読者です 読者をやめる 読者になる 読者になる

モンテカルロ法で円周率の計算

Matlab/Octave


今年度、Octaveで簡単な統計計算をする授業をやっているのですが、その中で大数の法則中心極限定理の説明をするための一例として、モンテカルロ法で円周率を計算する練習問題を出題しました。

座標軸上(0,0)-(1,1)に1辺が長さ1の正方形があり、それに接するような原点を中心とした半径1の円(の4分の1)を考えます。ランダムに点を落としたとき、原点からの距離が1より小さければ円内、距離が1より大きければ円外にあるとできるので、円内・円外の確率から面積を概算し、それをもとに円周率を求めるというもの。

ランダムに点を落とす数(N)と、その実験を何回繰り返すか(K)によって、大数の法則中心極限定理がなんとなくわかるといいなぁ、という希望です。

以下がそのコード。

tic;

K = 1000;
N = 100000;
est_pi = zeros(K,1);

for k=1:K
  x = rand(N,1);
  y = rand(N,1);
  count = sum(x.*x + y.*y < 1);
  est_pi(k) = 4*(count / N);  
end
fprintf(1, 'K=%d N=%d ==> pi=%f\n', K, N, mean(est_pi));

toc;

試しにK=1000、N=1000000で計算したところ3.141577となりました。Matlab 2012aでの実行時間は42.26秒 (MacBook Pro 13, Mid 2009, 2.53GHz Core 2 Duo, 8GB RAM) でした。

こんどはJuliaで並列計算させてみようかな。