ポワソン分布は単位時間あたりλ回起こる出来事の確率分布と説明されることが多い。

これをPythonのシミュレーションで確かめてみる。

import random
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from pandas import DataFrame, Series
from scipy.stats import poisson

%matplotlib inline

シミュレーションを行うには単位時間あたりの長さのほかにシミュレートする時間の長さが必要だ。

単位時間あたり4回起こる出来事を長さが2の時間の中にランダムに配置するとすると、8個のイベントを配置するということになる。

 [random.random() * 2 for _ in range(0, 8)]

[1.7265978877690655,
 0.05240790096914916,
 1.0092316367962268,
 1.7811176759439296,
 1.5008611575306552,
 1.4178547662306291,
 0.3848485501520824,
 1.9601703946981495]

単位時間1の中にいくつ入ったのかを数える。

[random.random() * 2 < 1 for _ in range(0, 8)].count(True)

3

これを10000回繰り返してヒストグラムを描く。

simulations = [
    [random.random() * 2 < 1 for _ in range(0, 8)].count(True)
    for _ in range (0, 10000)]
plt.hist(simulations, bins=range(0, 8+2), normed=True, color='gray')
plt.show()

これは実質的には試行回数が8で確率が1/2の二項分布のシミュレーションを行っていることになる。 8個のイベントしか配置していないから8回よりも多くなることはない。

一方ポワソン分布のパラメタは単位時間あたりの期待値だけで、8のような上限はない。そこで期待値が4に固定されるようにして、長さとイベント数を増やしてみる。

さらにシミュレーション結果の隣にポワソン分布の理論値のヒストグラムを並べるようにする。

lam = 4

lengths = [2, 4, 8, 16]

fig_columns = 2
fig_rows = len(lengths) / fig_columns + (1 if len(lengths) % fig_columns > 0 else 0)
fig, axes = plt.subplots(fig_rows, fig_columns)

for i, length in enumerate(lengths):
    num_events = int(length * lam) # 単位長さ当たりのイベント数num_events/length=lamになるようにする

    simulations = [
        # 長さlengthの中にnum_events個のイベントを配置し、1未満に入ったイベントの数を数える
        [random.random() * length < 1 for _ in range(0, num_events)].count(True)
        for _ in range (0, 10000)] # 10000回シミュレートする

    ax = axes[i / fig_columns, i % fig_columns]
 
    x_range = Series(range(0, int(poisson.ppf(0.99999, lam))))

    (_, _, bar_sim) = ax.hist(simulations, bins=x_range, normed=True, color='gray', width=0.5)
    bar_poi = ax.bar(x_range + 0.5, poisson.pmf(x_range, lam), color='white', width=0.5)

    ax.set_xlabel('count')
    ax.set_ylabel('probability')
    ax.set_title('length=%d, # of events=%d' % (length, num_events))

fig.legend((bar_sim[0], bar_poi[0]), ('simulation', 'poisson'), 'upper right', bbox_to_anchor=(1.2, 1))
fig.subplots_adjust(wspace=0.5, hspace=1)

plt.show()

長さとイベント数が多くなるにつれてシミュレーションの分布がポワソン分布に近づいていることがわかる。

参考にした文献：「ポワソン分布・ポワソン回帰・ポワソン過程」島谷健一郎（近代科学社）

wonderful cool something

ポワソン分布のシミュレーション