モンテカルロ法による円周率の計算をRで行う

モンテカルロ法を使用して、息抜きにRで円周率の計算を行う。

有名問題なので背景はこの辺を参照。

L=1
count = 1000000
x = runif(count,0,L) #0からLの間の乱数を発生させる
y = runif(count,0,L) #0からLの間の乱数を発生させる
pi = 4*sum(x**2+y**2<L**2)/count

countの数を増やして繰り返しplotとしてみるようにした。
※）Rに慣れている人はこんなコードは書かないはず。非常に時間がかかるので注意。

L=1
totalCount = 1000

result = c()
cvec = 1:totalCount
for (count in cvec){
  x = runif(totalCount,0,L) #0からLの間の乱数を発生させる
  y = runif(totalCount,0,L) #0からLの間の乱数を発生させる
  pi = 4*sum(x**2+y**2<L**2)/totalCount
  result = c(result,pi)
}
plot(result)
result

totalCount = 10000として実行してみた。
※）実行には２０分ほどかかった気がするｗ

・・・・・・・・・・
[9946] 3.13496 3.14224 3.14996 3.14064 3.13712 3.14476 3.13808 3.13960 3.15128
[9955] 3.14856 3.13508 3.13692 3.14312 3.13912 3.13412 3.14264 3.14492 3.13232
[9964] 3.13880 3.13900 3.14716 3.13640 3.14712 3.14848 3.14032 3.15060 3.14208
[9973] 3.13532 3.14204 3.14852 3.14672 3.14236 3.14304 3.15292 3.13712 3.13504
[9982] 3.14880 3.14492 3.12244 3.14516 3.14272 3.13480 3.14568 3.14504 3.13816
[9991] 3.14020 3.13856 3.15320 3.14652 3.14520 3.13796 3.14524 3.14280 3.14516
[10000] 3.13400

こんなものだったか・・・。