Help us understand the problem. What is going on with this article?

ビュフォンの針をRのanimationで可視化モンテカルロ法

More than 3 years have passed since last update.

この記事はR Advent Calendar 2015の3日目の記事です。

はじめに

数値解析の分野では、Rを使用せずに違う言語で描画プログラムを書いた経験のある方が多数いらっしゃると思います。
Rにはanimationパッケージというものがあり、プロットを結合してアニメーションとして出力することができます。

今回はそれを使った数値解析・モンテカルロ法の話です。

ビュフォンの針とは

ビュフォンの針(Buffon's needle problem)というのをご存知でしょうか。

長さ a の針を等間隔 d で引かれた平行線群のなかに落としたとき、針と平行線の交わる確率がどのようになるかという問題です。

このとき

  • 針の中心は平行線の存在する空間に一様に分布
  • 針と平行線のなす角度の分布は区間 (0, π) 上の一様分布

を満たすことができれば

\frac{2a}{\pi d}

となります。つまり、 d=2a とおけば

\frac{2a}{\pi d} = \frac{1}{\pi}

となり、円周率(の逆数)を得ることができます。

a=l と記述されるのが一般ですが、見難いため a にしています。

今回は、この問題をRで実装し、シミュレートして確認します。

実装

Rコード

buffon.needle = function(n = 100, sec = 0.1, l = 1.0, d = 1.0) {
  x <- y <- xd <- yd <- phi <- ctr <- numeric(0)
  sum <- 0; PI <- rep(NA, n)
  layout(matrix(c(1, 2), 2))
  for (j in 1:n) {
    # Draw parallel lines
    xlim <- c(-0.5 * l, 1.5 * l); ylim <- c(0, 2 * d)
    plot(1, xlim = xlim, ylim = ylim, type = 'n',
         xlab = '', ylab = '', axes = FALSE)
    abline(h = c(0.5 * d, 1.5 * d), lwd = 2, col = 'blue')

    # Simulate
    phi <- c(phi, runif(1, 0, pi))
    ctr <- c(ctr, runif(1, 0, 0.5 * d))
    x = c(x, runif(1, 0, l))
    y = c(y, sample(c(0.5 * d + ctr[j], 1.5 * d - ctr[j]), 1))
    xd = c(xd, 0.5 * l * cos(phi[j]))
    yd = c(yd, 0.5 * l * sin(phi[j]))

    # Draw Needle
    segments(x - xd, y - yd, x + xd, y + yd, col = 'gray')

    # Draw Mote Carlo
    xx = seq(0, pi, length = 200)
    plot(xx, 0.5 * l * sin(xx), ylim = c(0, 0.5 * d), type = 'l', bty = 'l',
         xlab = '', ylab = '', col = 'lightgray')
    idx = as.numeric(ctr > 0.5 * l * sin(phi)) + 1
    points(phi, ctr, col = c('red', 'blue')[idx])
    sum = sum + (ctr[j] <= 0.5 * l * sin(phi[j]))
    if (n > 0) {
      PI[j] = 2 * l * j/(d * sum)
    }
    ani.pause(sec)
  }

  return(PI)
}

animationの話だからといってやたらめったら変わったことはしていません。
今回は ani.pause() を使用して、段階的に表示することだけをしています。
実際の現場では ani.start(), ani.stop(), ani.record(), ani.replay() を使って録画することもあると思いますので、適したときに使用してください。

実行

animationパッケージを libraryrequire でロードしたのち、ユーザ定義関数を呼び出します。

> library('animation')
> buffon.needle()

buffon.needle 関数は引数を与えることもできます。

> buffon.needle = function(n = 100, sec = 0.1, l = 1.0, d = 1.0)

結果

アニメーション

針をプロットするのと、モンテカルロ法のアニメーションを見ることができます。

output2.gif

ビュフォンの針では下記の関数の積分値が解となるので、その積分=面積をモンテカルロ法で近似解を得ている。

\frac{L}{2}sin \theta

出力

a=l=1, d=1として関数を実行しているため、近似解は円周率に収束していくはずとなります。
その収束していく流れがわかるようにステップ毎の近似解を返却しています。

> buffon.needle(30)
 [1] 2.000000 2.000000 2.000000 2.666667 2.500000 2.400000 2.800000 3.200000 3.000000 2.857143 3.142857
[12] 3.000000 2.888889 2.800000 2.727273 2.666667 2.833333 2.769231 2.714286 2.857143 2.800000 2.750000
[23] 2.705882 2.823529 2.941176 3.058824 3.176471 3.111111 3.052632 3.157895

おわりに

Rというより数値解析の話が強かったかもしれませんが、数値解析方面をやっている方はどんどんRを使っていってほしいです。
僕も大学生のときにちょっと使っていた程度ですが、もっと使っていれば楽にできたことたくさんあったよなぁ。と今では思います。

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away