Maximaで学ぶ確率・統計（ヒストグラムの作成）

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

コインの裏表（乱数）
表が出た回数
表が出る確率
- 表が出たコインの枚数
- ヒストグラム
正規分布

コインの裏表（乱数）

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

random(n)$  /* n 未満の非負整数をランダムに出力します */

そこで、 $n=2$ の場合を考え、１を表に、０を裏に対応させることで、ランダムな出力をコイントスだと思うことができます。たとえば、コイントスを10回だけ試行する場合には

makelist(random(2), i, 1, 10);

とします。出力される結果が、入力のたびに変わることが確認できると思います。

表が出た回数

いま考えている問題では、10回試行したリストの各要素の和がそのまま、10回のコイントスのうちで表が出た回数に等しくなります：

sum(makelist(random(2), i, 1, 10)[n], n, 1, 10);

これもまた、入力のたびに出力が変わります。この量はまた、10枚のコインを１度に投げたときに、表が出たコインの枚数だとも思えます。

以下の説明ではこのようにみなすことにします。

表が出る確率

表が出たコインの枚数

１度に投げるコインの枚数を10枚から20枚に増やし、20枚のコインを投げるという試行を10,000回くりかえしてみましょう：

kill(all)$
L:makelist(sum(makelist(random(2), i, 1, 20)[n], n, 1, 20), k, 1, 10000)$

ここで作成したリストの各要素は、各試行ごとの投げた20枚のコインのうち表が出たコインの枚数です。

length(L);

から、リストの個数（試行数）は10,000だけあることがわかります。以下のようにして、どれだけの試行数で何枚のコインが表になったか数えます。

k:0$ while k<21 do(
N[k]:0,
i:1, while i<10001 do(
if L[i]=k then N[k]:N[k]+1 else N[k]:N[k],
i:i+1),
k:k+1)$

ここで N[k] は20枚のコインのうち k 枚のコインだけ表が出た試行数です。

ヒストグラム

表が出たコインの枚数ごとに、その事象が何回だけ出現したかを２次元データとしてプロットしたものがヒストグラム（散布図）です：

NL:makelist([k,N[k]],k,0,20)$
plot2d ([discrete, NL],[style,points],
[xlabel, "number of coins showing its face"],[ylabel, "number of events"]);

以上のようにして得られたヒストグラムが以下です：

f:id:pianofisica:20200830171219p:plain

10,000回の試行のうち、20枚中3枚以下、17枚以上のコインが表になったことはほとんどなく、10枚表が出たことが1,800回近くにのぼることが観測できます。

このような分布は正規分布（ガウス分布）としてよく知られています。

正規分布

平均（average）分散（variance）偏差（deviation）を求めます：

av: bfloat(sum(L[k],k,1,10000)/10000);
vr: bfloat(sum((L[k]-av)^2,k,1,10000)/10000);
dv: bfloat(sqrt(vr));

より

$\displaystyle{ \ {\rm Average}=10.04=:\mu\\ \ {\rm Variance}=5.092\\ \ {\rm Deviation}=2.257=:\sigma }$

です。３シグマが $3\sigma=6.77$ なので

$\displaystyle{ \ \mu\pm3\sigma=10.04\pm6.77}$

から、20枚中3以下、17枚以上のコインが表になる事象が３シグマ以上のまれなイベントであることがわかります。

上記の値を正規分布

$\displaystyle{ \ \rho(x)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right) }$

に当てはめ、ヒストグラムと合わせて図示してみます：

f(x):=1/sqrt(2*%pi*dv^2)*exp(-(x-av)^2/(2*dv^2))$
plot2d ([ [discrete, NL], 10000*f(x)], [x,0,20], [style, points, lines],
[color, blue, red],
[xlabel, "number of coins showing its face"],[ylabel, "number of events"]);

と入力して得られたのが次の図です：

f:id:pianofisica:20200830175321p:plain

ガウス分布関数は事象の出現確率を与えるので、ヒストグラムと比較するためには全試行回数を掛ける必要があることに注意します。

ヒストグラムが正規分布とよく合っていることがわかると思います。

今回は、乱数を使ってコンピュータ上で仮想的にコイントスを実現しました。

実際に、20枚のコインを一度に投げて表が出た枚数を数えることはできますが、それを10,000回もくりかえすことは到底できません。コンピュータを使えばそれが（仮想的にではありますが）瞬く間にできてしまいます。計算機の威力を実感できる題材だったのではないでしょうか？

キーワード：Maxima、確率・統計、乱数、データ作成、プロット