Pyhton を使って確率の計算をやってみる：中心極限定理と正規分布

簡単な確率の計算をPythonを使ってやってみましょう。

確率は何とも概念がつかみにくいものがあるのですが、実際に計算してみることで概念をつかんでみたいと思います。

この記事はこんな人におすすめです。

Pythonで確率の計算をしてみたい。
中心極限定理について知りたい

とりあえず簡単な問題からやっていきましょう。

プログラミング無料体験はこちら↓↓↓

確率の計算を行ってみる

サイコロを使った確率の計算

確率の定番、サイコロを使って確率を求めてみたいと思います。

今回のお題は以下の「演習で身につける統計学入門」からの出題です。

公平なサイコロAとBを振って出た目の和が10以上になる事象Dの起きる確率を求めなさい。

Pythonなどを使わずに普通に答えを求めた場合、まずサイコロを２回振って出た数の和が10以上になる組み合わせを見つける必要があります。

その組み合わせは、|4, 6|、|5, 5|、|5, 6|、|6, 4|、|6, 5|、|6, 6| の6通り。

つまり、6つの事象があるということ。

根本事象は 6 X 6 = 36 なので　 $P(A) = \frac{n(A)}{n(S)}$ 　から、 $P(D) = \frac{6}{36} = \frac{1}{6} = 0.16666$ 　、つまり16％ほどになります。

しかし毎回事象を数えるのは面倒なので、プログラムで簡単に計算できるようにしてみます。

count = 0 

# for 文を2回繰り返し、二つのサイコロを振って出た目の和が10以上になる組み合わせを求める。

for i in range (1, 7): 　# サイコロ Aを振る。1 以上 7 未満、つまり6
    for j in range (1, 7): 　# サイコロBを振る。1 以上 7 未満、つまり6
        if i + j >= 10: 　# もしA と Bの合計が10以上になる場合
           count += 1 　# countに1を足して代入

print(count/36) 　# countを根本事象36で割ると
0.16666666666666666

つまり確率は16％ほどとなります。
この計算の仕方だと、根本事象である36がわかっている必要がありました。

しかし根本事象がわからない場合もありますので、わからない場合を前提にしたプログラムを書くと以下。

A = 0
B = 0
# for 文を2回繰り返し、二つのサイコロを振って出た目の和が10以上になる組み合わせを求める。

for i in range (1, 7): 　# サイコロ Aを振る。1 以上 7 未満、つまり6
    for j in range (1, 7): 　# サイコロBを振る。1 以上 7 未満、つまり6
        A += 1 　# 根本事象を求める
        if i + j >= 10 or j + i >= 10:  #もしA と Bの合計が10以上になる場合
           B += 1  # B に1を足して代入

print(B/A)  # Bを根本事象A で割ると
0.16666666666666666

確率は同じように16％ほどとなります。

コインを使った確率の計算

では次に、これも確率の勉強では定番のコインを使って確率を求めてみます。

例題は同じように「演習で身につける統計学入門」からの出題です。

公平な効果を3回トスしたとき、

裏が1回だけ出る確率を求めなさい。

裏が1回以上出る確率を求めなさい。

という問題を解いてみましょう。

この問題を解く場合、表か裏かの離散確率分布になるので二項分布を使って確率を計算するのが妥当です。

考え方として表が0裏が1とし、どちらかが出る確率は1/2、つまり0.5で50％になりますね。

表が１回だけ出る確率

この考え方決めたうえで、scipy.statsというモジュールを使って確率を求めてみます。

from scipy.stats import binom  # scipy.stats モジュールから binom(二項分布)をインポート
win = binom.pmf(1, 3, 0.5)  # 3回投げて1回裏が出る確率を求める。

print(win)
0.375

実行すると、裏が一回だけで出る確率は37.5％になります。

裏が一回以上現れる確率

では、裏が1回以上現れる確率を求めます。

考え方としては、裏が一度も出ない確率を求めて100％から引くと、残りは裏が1回以上出る確率になります。

from scipy.stats import binom  # scipy.stats モジュールから binom(二項分布)をインポート
win2 = binom.pmf(0, 3, 0.5) # 3回投げて裏が一度も出ない確率を求める。 

print(1-win2)  # 1から裏が一度も出ない確率を引く
0.875

というわけで、コインを３回投げて裏が１回以上出る確率は87.5％になります。

中心極限定理とは

正規分布から確率を求める方法はいくつか書きましたが、今回は確率を勉強するうえで覚えておきたい言葉があります。

それは「中心極限定理」ですね。

確率を正規分布を使って求めるのですが、その時に知っておきたいのが「中心極限定理」です。

面倒くさそうな文字ですが、知っておくと確率の考え方を理解するのに役に立ちますので、ざっくりと覚えましょう。

簡単に言うと「母集団から取り出した標本の平均は、標本数を十分大きくしたとき正規分布に近づく（演習で身につける統計学より）」という定理です。

ではなぜ正規分布が都合がいいのかというと、以下の３点があるから（プログラミングのための確率統計）。

ある値 $\mu$ を中心として
その付近の値は出やすく
そこから外れた値ほど出にくくなる

世の中はそんな正規分布にあふれているので、使い勝手がいいんですね。

とりあえず、中心極限定理の数式を載せておきます。

１．標本平均の期待値は母平均に等しい

$E\left [ \bar{X} \right ] = \mu$

2. 標本平均の分散は母分散を標本数 n で割ったものに等しい。

$E\left [ (\bar{X}-\mu )^{2} \right ] = \frac{\sigma ^{2}}{n}$

というわけで、標本採集の回数を増やすとどうなるか見てみましょう。

標本採取の数を増やすシミュレーションを行う

下のコードは、6,7,8,9,10という数字の中から4個をランダムに取り出し、その平均を求める操作を40000回繰り返すものです。

import numpy as np
from matplotlib import pyplot as plt

t = list()
for i in range(40000):
　　s = np.random.randint(6, 11, 4)
　　t+= [sum(s)/len(s)]

plt.hist(t, bins=16)

この操作を行って得られた平均をグラフ化すると、以下のようなきれいな正規分布になりますね。

中心極限定理と標準化を使って確率を求める

では、実際に中心極限定理の概念と標準化を使って、確率を求めてみたいと思います。

例えば以下のような場合。

果物店の果物一個当たりの重量が平均257g で分散が63の場合、35個取り出した時の標本平均が260g以上になる確率を求めてみます。

Pythonにはこのような計算を簡単に行うためのパッケージが用意されていて、scipyもその一つになります。

以下コードになります。

from scipy.stats import norm 　# scipyのstatsからnorm(正規分布）読み込む。
p = norm.cdf(x = 260, loc = 257, scale = ((63/35)**0.5))

print(1-p)
0.0126736593387341

以上より、果物が260g以上になる確率は0.0126、つまり1.26％ですね。

もう一つやってみましょう。

英語の成績が平均が67点で分散が28の場合、40人を無作為に取り出した時の平均が65点から69点の範囲にある確率を求める。

この場合、65点のまでの確率と69点までの確率を求め、69点までの確率から65点までの確率を引くと、その範囲にある確率を求めることができます。

from scipy.stats import norm
p1 = norm.cdf(x = 65, loc = 67, scale = ((28/40)**0.5))
p2 = norm.cdf(x = 69, loc = 67, scale = ((28/40)**0.5))

print(p2-p1)
0.9831725905172433

以上の計算により、確率は98％になります。

しかし、以下のような場合は別のやり方が必要になります。

例えば、

英語の試験を受けた生徒の中から無作為に選んだ37人の平均と分散は63と36。

この時、母集団の平均が60未満になる確率を求める。

この場合、母分散 $\sigma ^{2}$ には以下の関係があります。

$\sigma ^{2} = S^{2}n/(n-1)$

これを使って、中心極限定理に当てはめ標準化すると以下のようになります。

$Z = \frac{\bar{X}-\mu }{S/\sqrt{n-1}}$

といっても、scipyを使っても計算できるので、サクッと計算してしまいましょう。

と、その前に、標本データの平均は母集団の平均と違う可能性が大きいです。

そのため、母集団の分散は標本データの分散よりも大きくなるので、補正する必要があるんですね。

それが（n – 1)の意味になり、今回の問いの場合は（37－1）というわけです。

scipyでも同じ考え方で、確率を求めます。

from scipy.stats import norm 
p3 = norm.cdf(x = 60, loc = 63, scale = ((36/(37-1))**0.5)) 
print(p3) 0.0013498980316300933

計算してみると確率は0.0013、つまり0.13%ということになります。

まとめ

この記事では、Pythonを使って簡単な確率の計算を行います。確率の計算に慣れるため、サイコロやコインを使った問題から始め、中心極限定理についても説明しました。

Pythonを使うことで、確率の計算が簡単にできることが分かります。

また中心極限定理を使うことで、標本の平均が母集団の平均にどのように近づくかを理解することもできます。