Pythonでポアソン分布と戯れる

はじめに

統計分布と戯れたシリーズ第三段ということでポアソン分布をお届けする。

ポアソン分布

• 書き方 Po(λ)
• 確率関数 $$ P(X=x) = \frac{λ^x}{x!}exp(-λ) $$
• 期待値 λ
• 分散 λ

確率関数を書いてみよう

def my_poisson(x, l):
    ret = (np.float_power(l, x) * np.exp(-l) / math.factorial(x)) 
    return ret

自前で定義した関数をそのまま使うケースは少ないだろうが、powerではなくfloat_powerを使わないと、xが大きくなった時に計算がおかしくなる。

確率関数を描画してみよう

定義した確率関数を使って、$Po(5)$ に従うポアソン分布を描画してみる。

# X～Po(5)
xs = []
ys = []
for x in range(20+1):
    xs.append(x)
    ys.append(my_poisson(x, 5))

fig = plt.figure(figsize=(8,5))    
ax = fig.add_subplot(1,1,1)
ax.bar(xs, ys)

こんな感じ。

なお、scipyでやるとこんな感じだ。検算するとよいだろう。

# Scipyの確率関数で検算する
from scipy.stats import poisson
rv = poisson(5)
ys_scipy = rv.pmf(xs)

fig = plt.figure(figsize=(8,5))    
ax = fig.add_subplot(1,1,1)
ax.bar(xs, ys_scipy)

再生性を確認しよう

$ X1～Po(λ1), X2～Po(λ2) $ である時、$X1 + X2～Po(λ1+λ2)$ になることが知られている。 これをポアソン分布の再生性という。

そこで、$X1～Po(3), X2～Po(5) $ に従う $X1, X2 $ の和が、$Po(8)$ に従うかやってみよう。

まずは、$X1＋X2$ の分布を求めてみる。

# X1～Po(3), X2～Po(5)
x1s = []
y1s = []
x2s = []
y2s = []
l1 = 3
l2 = 5

for x in range(15):
    x1s.append(x)
    y1s.append(my_poisson(x, l1))

for x in range(15):
    x2s.append(x)
    y2s.append(my_poisson(x, l2))
    
y1s /= np.array(y1s).sum()
y2s /= np.array(y2s).sum()

# 10000回X1, X2のサンプルを取得し、X1+X2の分布(出現回数)を取得する
sample_size = 10000
l1_l2 = [0] * 20
for i in range(sample_size):
    # X1から1つサンプルを取得
    x1_sample = np.random.choice(x1s, p=y1s)
    # X2から1つサンプルを取得
    x2_sample = np.random.choice(x2s, p=y2s)
    x1_x2[x1_sample + x2_sample] += 1

# 分布(出現回数）を標本数で割って確率を求める
x1_x2 = np.array(x1_x2)
x1_x2 = x1_x2/sample_size 

次に $Po(8)$ を求めよう。

xs = []
ys = []
for x in range(20):
    xs.append(x)
    ys.append(my_poisson(x, 8))

$X1＋X2$ の分布と $Po(8)$ の分布を重ねてみよう

fig = plt.figure(figsize=(8,5))    
ax = fig.add_subplot(1,1,1)
ax.bar(range(20), x1_x2[:20], color="r", alpha=0.3)
ax.bar(xs, ys, color="b", alpha=0.3)

紫の色のところが重なった部分になるが、ほぼ重なっていることが分かる。

中心極限定理を確認しよう

二項分布の時と同様に中心極限定理を確認しよう。
中心極限定理によれば、nが大きい時、ポアソン分布からの標本平均を用いた統計量 $\sqrt{n}(\bar{X_{n}}-λ)$は、正規分布$N(0,λ)$ に分布収束する。

$λ=5, n=100$ の条件でやってみよう。

# X～Po(5)の確率関数によるサンプリングの準備
n = 100
l = 5
x1s = [] 
y1s = []
for x in range(30):
    x1s.append(x)
    y1s.append(my_poisson(x, l))

y1s /= np.array(y1s).sum()    

# n=100で標本平均を求め、統計量を計算することを1000回繰り返す
traial_time = 1000
targets = []

for i in range(traial_time):
    x1_samples = []
    # n=100サンプリングする
    for j in range(n):
        x1_sample = np.random.choice(x1s, p=y1s)
        x1_samples.append(x1_sample)

    # n=100の標本平均を求める
    x1_samples_mean = np.array(x1_samples).mean()
    # print(x1_samples_mean)
    
    # 統計量を計算する
    target = math.sqrt(n) * (x1_samples_mean - l) 
    targets.append(target)

# 描画する
fig = plt.figure(figsize=(8,5))    
ax = fig.add_subplot(1,1,1)
ax.hist(targets, bins=200, color="b")

そうすると以下のような分布が得られる。

これが正規分布 $N(0,λ)$ に収束するということだから、これもグラフを描いて比較してみよう。
normal というのは正規分布の確率密度関数である。

# 正規分布を表示
fig = plt.figure(figsize=(8,5))
ax = fig.add_subplot(1,1,1)
x = np.linspace(-7, 7, 100)
y = normal(x, 0, l)
ax.plot(x, y, color="r")

上を実行すると下のようなグラフが得られる。
$\sqrt{n}(\bar{X_{n}}-λ)$ の分布が、下の正規分布に非常に近いということが分かる。

少数の法則を確認しよう

二項分布 $B(N, p)$ において $p=λ/N$ とし、$N$ を大きくするとポアソン分布 $Po(λ)$ に分布収束する。これを少数の法則という。

$N=100$, $λ=3$ で実験してみよう。

まずは、$X～Bin(N, λ/N)$ のからのサンプリング

# X～Bin(N, l/N)の確率関数によるサンプリングの準備
N = 100
l = 3
x1s_bin = [] 
y1s_bin = []

# X～Bin(N, l/N) の確率関数によるサンプリング
for x in range(N+1):
    x1s_bin.append(x)
    y1s_bin.append(binomial(x, N, l/N))

次に $X～λ(3)$ のサンプリング

x1s_po = []
y1s_po = []
for x in range(20):
    x1s_po.append(x)
    y1s_po.append(my_poisson(x, 3))

fig = plt.figure(figsize=(8,5))    
ax = fig.add_subplot(1,1,1)
ax.bar(x1s_po[:10], y1s_po[:10])

$X～Bin(N, λ/N)$ と $X～λ(3)$ を重ねてみよう。

fig = plt.figure(figsize=(8,5))    
ax = fig.add_subplot(1,1,1)
ax.bar(x1s_bin[:10], y1s_bin[:10], color="red", alpha=0.3)
ax.bar(x1s_po[:10], y1s_po[:10], color="blue", alpha=0.3)

以下の紫のところが重なった部分だが、きれいに重なっているのが分かる。