この記事は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パッケージを library か require でロードしたのち、ユーザ定義関数を呼び出します。
> library('animation')
> buffon.needle()
buffon.needle 関数は引数を与えることもできます。
> buffon.needle = function(n = 100, sec = 0.1, l = 1.0, d = 1.0)
結果
アニメーション
針をプロットするのと、モンテカルロ法のアニメーションを見ることができます。
ビュフォンの針では下記の関数の積分値が解となるので、その積分＝面積をモンテカルロ法で近似解を得ている。
\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を使っていってほしいです。
僕も大学生のときにちょっと使っていた程度ですが、もっと使っていれば楽にできたことたくさんあったよなぁ。と今では思います。