Maximaで学ぶ確率・統計（モンテカルロ法）

数学の具体的な計算にMaximaを使って、数学もMaximaも同時に学んでしまいましょう。今回はMaximaを使った確率・統計の問題として、モンテカルロ法で円周率を求める問題を考えてみたいと思います。本記事のMaximaの学習事項としては、離散データの作成・処理の方法、散布図の作成を学びます。ただし、図の描画の確認はwxMaxima上でgnuplotを用いた場合で行っています。インターフェイスなどの動作環境の違いによって適宜変更点があるかもしれません。

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

乱数

Maximaには（擬似）乱数が用意されています：

random(1.0)$  /* 1.0 未満の負でない小数をランダムに出力します */

そこで２つの乱数 $(x,y)$ を生成します。1.0 未満の２つの小数をランダムに出力したとき、そのデータ点を多数（その数を $P$ としましょう）プロットしていくと、 $0\leq x\leq1, 0\leq y\leq1$ という領域をまんべんなく埋め尽くしていくはずです。

実際、1000個のデータ点をプロットしてみると

kill(all)$ P:1000$
p:0$ while p<P do(
x_point:random(1.0),
y_point:random(1.0),
xy[p]:[x_point,y_point],
p: p+1)$
XY:makelist(xy[p],p,0,P-1)$
plot2d ([discrete, XY],[style, points],
[color, blue],[xlabel, "x"],[ylabel, "y"]);

より

f:id:pianofisica:20200901115316p:plain

のようになって $0\leq x\leq1, 0\leq y\leq1$ の範囲にほぼ一様にデータ点が分布します。

モンテカルロ法

そこで、これらのデータ点のうち、 $x^2+y^2\leq1$ という条件を満たすデータ点の数（この数を $S$ としましょう）を求めてみます。これは第１象限内の半径１の円のなかにあるデータ点の数に等しいので、データ点の数の比 $S/P$ をとれば、半径１の円の４分の１におおよそ等しいはずです：

p:0$ while p<P do(
n[p]:0,
if (xy[p][1]^2+xy[p][2]^2)<=1 then n[p]:1 else n[p]:0,
p: p+1)$
S: sum(n[p],p,0,P-1);
bfloat(4*S/P); /* 円周率と見比べるために４倍しておきました */

もちろんデータ点が少なければ誤差が大きいわけですが、データ点を多くとり、試行回数を増やしていくほど良い近似値が得られるだろうと期待されます。これをMaximaを使ってみていくことにしましょう。そこで以下では、1,000個のデータ点をとり、試行回数を3,000回に設定して、上記の範囲内に入ったデータ点の個数を考えてみます。

円周率の計算

kill(all)$ T:3000$ P:1000$
t:1$ while t<T+1 do(
s[t]:0,
p:1, while p<P+1 do(
x_point:random(1.0),
y_point:random(1.0),
if (x_point^2+y_point^2)<=1 then s[t]:s[t]+1 else s[t]:s[t],
p: p+1),
t: t+1)$

s[t] は t 回目の試行で、1,000個のデータ点のうち、データ点が半径１の円の範囲内に入った場合の数です。

範囲内に入ったデータ点の個数ごとに、どれくらいの試行数でそれが実現したかを数えます：

h:0$ while h<P+1 do(
M[h]:0,
t:1, while t<T+1 do(
if s[t]=h then M[h]:M[h]+1 else M[h]:M[h],
t:t+1),
h:h+1)$

上の情報をもとにヒストグラムを作成します：

H:makelist([h,M[h]],h,0,P)$
plot2d ([discrete, H],[x, 7/10*P, 9/10*P],[style,points],
[xlabel, "number of score"],[ylabel, "number of events"]);

すると

f:id:pianofisica:20200901122248p:plain

のようなヒストグラムが得られます。

平均・分散・標準偏差を求めて

av: bfloat(sum(s[t],t,1,T)/T);
vr: bfloat(sum((s[t]-av)^2,t,1,T)/T);
dv: bfloat(sqrt(vr));

正規分布にあてはめてプロットすると

f(x):=1/sqrt(2*%pi*dv^2)*exp(-(x-av)^2/(2*dv^2))$
plot2d ([ [discrete, H], T*f(x)],[x, 7/10*P, 9/10*P],[style, points, lines],
[color, blue, red],[xlabel, "number of score"],[ylabel, "number of events"]);