モンテカルロシミュレーションによる円周率の推定

簡単なモンテカルロシミュレーション

円周率をモンテカルロシミュレーションで求められるって話は何度か聞いてたけど実際に学習したことなかったからやってみた.
円周率は直径に対する比率のことらしい. 原点より0<=x,y<=1の範囲でランダムにプロットした際に半径1の円に内包した率を考える.
これは, 一辺が2の正方形(原点から+-方向に1だから)の面積に対する半径1の円の面積が占める割合と等価である.
したがって、半径1であることから円の面積と円周率が同じであるため上述の割合を計算することで確率論的に円周率を求められるということである.

実験

Githubのリポジトリhttps://github.com/totomo/math/monte_carloにあげたプログラムで数値実験を行った. 言語はCを採用した. 乱数発生器については面倒なのでmath.hに標準で定義されているrand関数を用いた.

結果

下記に実行結果を示す.

$ make clean
$ make
$ make run
./estpi 10000
result: 3.1712000

estpi(実行ファイル)の第一引数には試行回数を設定してある. 中心極限定理より試行回数にしたがってシミュレーション結果も円周率に漸近するはずである. 検証のために下記のようなシェルスクリプトでestpiに与える引数を変更してファイルに書き出した.

rm -f result.dat; 
for((i=100;i<100^5;i+=100)) 
do ./estpi $i >> result.dat; 
done

以下に横軸を試行回数、縦軸を推定した円周率としてプロットしたグラフを示す.
グラフから試行回数にしたがって円周率に漸近していることがわかる.