11
13

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

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

Posted at

この記事は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を使っていってほしいです。
僕も大学生のときにちょっと使っていた程度ですが、もっと使っていれば楽にできたことたくさんあったよなぁ。と今では思います。

11
13
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
11
13

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?