概要
- 各分布の関係性をざっくり紹介
- 分布を計算・可視化するPythonコード作成
- 条件を変えて実行し、近似誤差の大きい場合と小さい場合があることを確認
参考書籍
各分布の関係をざっくりと
超幾何分布と二項分布
超幾何分布では
- 赤青のボールが20個ずつ入った箱から、ボールを10個取り出した時、赤が3個である確率はいくらか？
といった問題を考えます。1回目に赤が出る確率は50%ですが、1回目の影響で残りの赤青の個数が変わるので、2回目は確率が変わります。これは非復元抽出と呼ばれ、各試行が独立ではありません
（参考：biostatistics - 超幾何分布）
一方、二項分布では
- 赤青のボールが1個ずつ入った箱から、ボールを1個取り出して戻す試行を10回繰り返した時、赤が3回出る確率はいくらか？
といった問題を考えます。各試行で元の状態に戻すので、各試行は独立であり、復元抽出の問題です。
よくある例だと「コインを10回投げて表が3回出る確率はいくらか？」というのと同じです
（参考：統計WEB - 二項分布）
超幾何分布は、全体の個数（箱の中のボールの数）が抽出数よりも十分に大きい場合、二項分布で良く近似できます
二項分布とポアソン分布と正規分布
二項分布は、
- 試行回数（コインを投げる回数）が十分多い場合、正規分布で良く近似できます
- 1回の試行の確率が極めて小さい場合（例えば、歪んでいて1000回に1回しか表が出ないコイン）、ポアソン分布で良く近似できます
（近似するのは計算が楽になるからです）
分布の計算と近似誤差の確認
Pythonコード（関数部）
import scipy as sc
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib
# functions ------------------------
def binomial(n, x, p):
'''
二項分布
確率pの事象が、n回の試行中でx回起こる
各試行は独立
例：サイコロで3が出る事象（確率p=1/6）が、サイコロをn=10回振った時にx=5回起こる）
nCx * p^x * (1-p)^(n-x)
'''
return (
sc.special.comb(n, x)
* np.power(p, x)
* np.power(1 - p, n - x)
)
def poisson(n, x, p):
'''
ポアソン分布
p->0のとき、二項分布はポアソン分布となる
'''
mean = n * p
return (
np.power(mean, x)
* np.exp(-mean)
/ sc.special.factorial(x)
)
def gaussian(n, x, p):
'''
正規分布
n->ooのとき、二項分布は正規分布となる
'''
mean = n * p
variance = np.power(n * p * (1 - p), 1/2)
return (
1.0
/ (np.power(2 * np.pi, 1/2) * variance)
/ np.exp(
np.power(x - mean, 2)
/ (2 * np.power(variance, 2))
)
)
def hypergeometric(n, x, p, N):
'''
超幾何分布
各試行は独立ではない
非復元抽出のため全数Nを考える必要あり
'''
target_size = N * p # 全体中の当たりの数
return (
sc.special.comb(target_size, x) # 当たりx個の組み合わせ
* sc.special.comb(N - target_size, n - x) # ハズレn-x個の組み合わせ
/ sc.special.comb(N, n) # n個の全組み合わせ
)
def calc_distributions(n, p, cols, N=0):
'''
ex.:
n = 30 # 試行回数
p = 0.16 # 当たりの割合
N = 1000 # 全数（超幾何分布でのみ使用）
'''
# 結果格納の配列を用意
bi_arr = np.zeros(n + 1, dtype=np.float)
po_arr = np.zeros(n + 1, dtype=np.float)
ga_arr = np.zeros(n + 1, dtype=np.float)
for x in np.arange(n + 1):
# 結果を配列に格納
bi_arr[x] = binomial(n, x, p)
po_arr[x] = poisson(n, x, p)
ga_arr[x] = gaussian(n, x, p)
# 結果をデータフレームに格納
df = pd.DataFrame()
df['二項分布'] = bi_arr
df['ポアソン分布'] = po_arr
df['正規分布'] = ga_arr
if N > 0:
hy_arr = np.zeros(n + 1, dtype=np.float)
for x in np.arange(n + 1):
hy_arr[x] = hypergeometric(n, x, p, N)
df['超幾何分布'] = hy_arr
return df[cols]
def visualize(df,n, p, N=0):
plot_params = {
'linestyle':'-',
'markersize':5,
'marker':'o'
}
colors = [
'black',
'blue',
'green',
'purple',
]
plt.close(1)
figure_ = plt.figure(1, figsize=(8,4))
axes_ = figure_.add_subplot(111) # Axes作成
for i, col in enumerate(df.columns):
if i == 0:
plot_params['linewidth'] = 2
plot_params['alpha'] = 1.0
else:
plot_params['linewidth'] = 1
plot_params['alpha'] = 0.4
plot_params['color'] = colors[i]
axes_.plot(
df.index.values, df[col].values,
label = col,
**plot_params,
)
plt.legend()
plt.xlabel('当たりが出る回数x')
plt.ylabel('n回中x回で当たりが出る確率')
title = '当たりの割合p:{p:.2f}, 試行回数n:{n}'.format(p=p, n=n)
if '超幾何分布' in df.columns:
title += ', 全数N:{N}'.format(N=N)
plt.title(title)
xmax = n * p * 2
if xmax<10: xmax = 10;
plt.xlim([0, xmax])
# （非復元抽出を考えると「回数」という表現はイマイチ）
条件を変えて実行、誤差確認
サイコロを100回振って、1が何回出るか
n = 100; p = 0.167
df = calc_distributions(n, p, cols=['二項分布', 'ポアソン分布', '正規分布'])
visualize(df, n, p)
二項分布と正規分布が良く一致
100面体を100回振って、1が何回出るか
n = 100; p = 0.01
df = calc_distributions(n, p, cols=['二項分布', 'ポアソン分布', '正規分布'])
visualize(df, n, p)
二項分布とポアソン分布が良く一致
5000人の住民のうち10%が子ども。100人選んだ中に子どもは何人いるか？
n = 100; p = 0.1; N = 5000
df = calc_distributions(n, p, cols=['超幾何分布','二項分布'], N=N)
visualize(df, n, p, N)
超幾何分布と二項分布は良く一致
500人の住民のうち10%が子ども。100人選んだ中に子どもは何人いるか？
n = 100; p = 0.1; N = 500
df = calc_distributions(n, p, cols=['超幾何分布','二項分布'], N=N)
visualize(df, n, p, N)
N=5000の場合と比べて、誤差が目立つ
おわり