概要

各分布の関係性をざっくり紹介
分布を計算・可視化する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の場合と比べて、誤差が目立つ

おわり