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

今年度、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で並列計算させてみようかな。