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

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

コインの裏表（乱数）
表が出た回数
表が出る確率
- 表が出たコインの枚数
正規分布

コインの裏表（乱数）

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

import random as rd
rd.randint(0, 1)  # randint(a,b) で a 以上 b 以下の整数をランダムに出力

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

import random as rd
def r(k): return rd.randint(0, 1)
R = []
for k in range(10):
      R.append(r(k))
print(R)

とします。たとえば

$\ \displaystyle{[1, 1, 1, 1, 0, 1, 0, 1, 0, 1]}$

のような結果が出力されます。もちろんいま乱数を使っているので、出力される結果は入力のたびに変わります。実際に何度か同じコマンドを実行してみると、毎回出力が変わることが確認できると思います。

表が出た回数

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

sum(R)

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

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

表が出る確率

表が出たコインの枚数

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

import random as rd
def r(k): return rd.randint(0, 1)
def l(t): return sum(R)
L = []
for t in range(10000): 
      R = []
      for k in range(20): 
            R.append(r(k))
      L.append(l(t))

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

ここで、上で作成したリストからヒストグラムを描画してみます。

import matplotlib.pyplot as plt
plt.hist(L,bins=20, range=(0, 20))
plt.show()

によって出力されるのが

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

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

正規分布

平均（mean value）標準偏差（standard deviation）を求めます：

import statistics as st
m = st.mean(L)
s = st.stdev(L)
print('平均：　{}'.format(m))
print('偏差：　{}'.format(s))

より

$\displaystyle{ \ {\rm Mean}=10.02=:\mu\\ \ {\rm Deviation}=2.24=:\sigma }$

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

$\displaystyle{ \ \mu\pm3\sigma=10.02\pm6.72}$

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

上記の値を正規分布

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

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

import sympy as sp
x = [t*20/10000 for t in range(10000)]
plt.plot(x, [1/(sp.sqrt(2*sp.pi*s**2))*sp.exp(-1/(2*s**2)*(t-m)**2)*10000 for t in x])
plt.hist(L,bins=20, range=(0, 20))
plt.show()

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

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

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

次の記事では、乱数を使った数値計算の応用として、円周率を求める問題を取り扱っています。あわせて読んでいただけると、よりよくPythonの使い方が学べると思います。

pianofisica.hatenablog.com

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

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

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