モンテカルロ法で円周率を求める
モンテカルロ法で円周率を求める
投稿日2017.8.28
対象となるある現象がいくつかの確率変数の関数で表現できる場合（数値モデル）、各変数で仮定される確率分布に沿った標本値を大量に生成（乱数）して、その計算結果の分布を推定（結果の推定）する方法です。
数値モデル → 乱数挿入 → 結果の推定
モンテカルロ法を使用して円周率 (π) を求めてみます
上図の四角の中にランダムに点を打ってみます
四角の面積は 2×2=4、円の面積は π × 1^2 = π
打った点の個数をn、円の中に入った個数をmとします
非常に多くの点を打つことで四角の中は塗りつぶされます
そこで次のような比率を考えます
点を多く打つことで、全体の点の数と円内の点の数との比が
四角の面積と円の面積の比に等しくなるはずです．
n : m = 4 : π
π＝4 × m / n （数値モデル）
Rを使用してnの乱数を多く作成すれば円周率が求められることになります
ではやってみます
一様乱数から(x,y)を作成します （乱数）
-1<x<1, -1<y<1
試しに5組分だけ作成してみます
dat<-matrix(runif(5*2,-1,1),2)
イメージ
一行目をｘ、二行目をｙとして x^2 + y^2 < 1 となる点の個数が
円内の個数mということになります
m<-colSums(dat^2) < 1
全ての点をプロットしてｍは赤で示します
plot( dat [1,], dat [2,], col = (2-m) , pch = 20 )
m はcolSums(dat^2) < 1の場合にはTRUE＝1、それいがいはFALSE＝0 ここがポイント
col=(2-m)とした場合には円内に入ったら1（黒）、円外では2（赤）となります
これまでの式をまとめて100個打ってみます
dat<-matrix(runif(100*2,-1,1),2)
m<-colSums(dat^2) < 1
plot( dat [1,], dat [2,], col = (2-m) , pch = 20 )
分かりにくいので10000個打ってみます
dat<-matrix(runif(10000*2,-1,1),2)
m<-colSums(dat^2) < 1
plot( dat [1,], dat [2,], col = (2-m) , pch = 20 )
おっ！円が見えてきました！！！
ちょっと感動
dat<-matrix(runif(n*2,-1,1),2)
m<-colSums(dat^2) < 1
plot( dat [1,], dat [2,], col = (2-m) , pch = 20 )
結果の推定
π＝4 × m / n を推定します
上記の結果より
π＝4 × m / 10000 = 3.1388 (乱数なので結果は異なります)
まとめて関数にしてみます
pai<-function(n){
dat<-matrix ( runif ( n*2 , -1 , 1 ) , 2 )
m<-colSums( dat^2 ) < 1
pl<-plot( dat [1,] , dat [2,] , col = (2-m) , pch = 20 , xlab = n , ylab = "" )
return( ( 4*sum( m ) ) / n )
}
nを100,500,1000,5000と変更して図を描いてみます
par(mfrow = c(2,2))
pai(100) ;pai(500) ;pai(1000) ;pai(5000)
dev.off()
n=100 π＝ 3.4
n=500 π＝ 3.128
n=1000 π＝ 3.144
n=5000 π＝ 3.1352
