ポアソン分布

【Python統計学】ポアソン分布

この記事の内容（ポアソン分布）は二項分布の知識を前提とするので、簡単に復習しておきます (二項分布の詳細についてはこちらの記事を参照してください)。
 
$2$ 種類の事象 $A$ と $B$ が、それぞれ確率 $p,\ 1-p$ で起こる試行を $n$ 回繰り返すとき、$A$ が $k$ 回起こる確率は
 \[f(k)=P(X=k)={}_{n}\mathrm{C}_{k}\,p^{k}(1-p)^{n-k}\tag{1}\]
で与えられます。この式を使うと、コインを $10$ 枚投げて表が $3$ 枚出る確率などを計算できます。
 
しかし、$n$ に比べて $p$ が極端に小さな値をとるような（すなわち滅多に起こらないような）事象を扱う場合には少し扱いにくい公式です。

航空機事故の年間発生率

滅多に起こらない事象の代表例として航空機事故があります。国土交通省の運輸安全委員会によると、$2010$ 年から $2019$ 年の $10$ 年間に起きた大型機の事故は $29$ 件でした。$1$ 年間に約 $3$ 件の事故が起きていることになります。

この統計データをもとに、$1$ 年に $k$ 件の航空事故が起きる確率を考えてみましょう。$1$ 日当たりの航空事故発生率は
 \[p=\frac{3}{365}\tag{2}\]
です。ただし、同じ日に複数の事故は起こらないと仮定しています。$1$ 年 ($365$ 日) に $k$ 件の交通事故が発生する確率は
 \[f(k)={}_{365}\mathrm{C}_{k}\,p^k\,(1-p)^{365-k}\tag{3}\]
で与えられます。scipy.stats.binom() を使って、$1$ 年に $1$ 件の事故が発生する確率を計算してみます。

# POISSON DISTRIBUTION

# In[1]

import numpy as np
from scipy.stats import binom
from scipy.stats import poisson
import matplotlib.pyplot as plt

# 1日あたりの航空機事故発生率
p_day = 3 / 365

# 1年に1件の航空機事故が発生する確率を二項分布で計算する
p1 = binom.pmf(1, 365, p_day)

print("{:.5f}".format(p1))
# 0.14874

このように、現代ではコンピュータを使って二項分布の表式のままで計算できますが、数十年前（1980 年頃）は、コンピュータを使ったとしても $(1-p)^{364}$ のような計算を処理するのは大変でした。まして、ポアソン (Simeon Denis Poisson) が生きていた時代にはコンピュータなど存在しなかったので、$(1-p)^{364}$ などを計算するのは不可能であり、もう少し扱いやすい式が必要でした。

証明は省略しますが、二項分布の確率変数 $X$ の期待値 $\mu=np$ を一定に保ちながら $n\rightarrow \infty$ の極限をとると ($p\rightarrow 0$)、
 \[f(k)=\frac{\mu^k}{k!}e^{-\mu}\tag{4}\]
という式になることが知られており、これを ポアソン分布 (Poisson distribution) とよびます。$n$ と $p$ が期待値 $\mu$ に吸収されて表式から消えていることに注目してください。ポアソン分布の形は確率変数の期待値 $\mu$ だけで決まります。

# In[2]

# 確率変数
x = np.arange(0, 11)

# ポアソン分布
y = poisson.pmf(x, 3)

# FigureとAxes
fig = plt.figure()
ax = fig.add_subplot(111)
ax.grid()
ax.set_title("Poisson distribution", fontsize=18, pad=8)
ax.set_xlim(0, 10)
ax.set_ylim(0, 0.25)
ax.set_xlabel("k", fontsize = 15, labelpad=8)
ax.set_ylabel("P(X=k)", fontsize = 15, labelpad=8)

# データをプロット
ax.plot(x, y, marker="o", color="blue")

plt.show()

$1$ 年に $2$ 件あるいは $3$ 件の事故が発生する確率が最も高いことがわかります。事故が $0$ 件である確率は $5$ % なので、航空機事故が稀であるとはいえ、全くない年も珍しいことになります。実際、運輸安全委員会のデータによると、1974 年から 2019 年の間に大型機の事故件数が $0$ だった年は 2010 年だけでした。
 
scipy.stats.poisson.cdf() を使うと累積確率分布を計算できます。たとえば、ある年の大型機の事故が $3$ 件以下である確率 $P(X \leq 3)$ は次のように計算できます。

# In[3]

# ある年の大型機の事故が3件以下となる確率
p_leq3 = poisson.cdf(3, 3)

print("{:.5f}".format(p_leq3))
# 0.64723

$1$ から $P(X \leq 3)$ を引けば、ある年の大型機の事故が $4$ 件以上発生する確率 $P(X \geq 4)$ を計算できます。

# In[4]

# ある年の大型機の事故が4件以上となる確率
p_geq4 = 1 - p_leq3

print("{:.5f}".format(p_geq4))
# 0.35277

二項分布で $np$ 一定のまま極限をとったので、ポアソン分布の期待値も $np$ です。また、二項分布の分散
 \[\sigma^2=np(1-p)=\mu\left(1-\frac{\mu}{n}\right)\tag{5}\]
において、$n \rightarrow \infty$ の極限をとると、
 \[\sigma^2=\mu\tag{6}\]
となって、分散も期待値に等しいことがわかります。scipy.stats.poisson.var() で期待値 $3$ のポアソン分布にしたがう分散が $3$ になることを確認しておきましょう。

# In[5]

# 平均3のポアソン分布の分散
var = poisson.var(3)

print(var)
# 3.0

参考資料：運輸安全委員会ホームページ
https://jtsb.mlit.go.jp/jtsb/aircraft/air-accident-toukei.php

numpy.random.poisson()

numpy.random.poisson() は、ポアソン分布にしたがう乱数を生成します。

numpy.random.poisson(lam=1.0, size=None)

lam で分布の期待値、size で生成する乱数の個数を指定します。
 
以下のコードを実行すると、期待値 $\mu=1.5$ のポアソン分布にしたがう乱数を $1000$ 個生成し、各数字の発生比率をヒストグラムで表示します。

# NUMPY_POISSON

import numpy as np
import matplotlib.pyplot as plt

# 疑似乱数ジェネレータのシードを設定
np.random.seed(0)

# μ=1.5のポアソン分布に従う1000個の乱数
k = np.random.poisson(1.5, 1000)

# 階級値の設定
bins = np.arange(0, 10)-0.5

# FigureとAxes
fig = plt.figure()
ax = fig.add_subplot(111)
ax.grid()
ax.set_title("Poisson distribution", fontsize=18, pad=8)
ax.set_ylim(0, 0.4)
ax.set_xlabel("k", fontsize = 15, labelpad=6)
ax.set_ylabel("P(X=k)", fontsize = 15, labelpad=6)

# 乱数の生成比率をプロット
ax.hist(k, bins, rwidth=0.6, density=True, color="darkorange")
plt.show()

上のコードでは Axes.hist() に density=True で正規化して各乱数の割合をプロットしています。このキーワードを省くか False を指定すれば各乱数の個数をプロットします。

scipy.stats.poisson()

scipy.stats.poisson() はポアソン分布オブジェクトを生成します。

scipy.stats.poisson(*args, **kwds)

オブジェクトに備わるメソッドを使って乱数や確率密度関数、各種統計量を取得できます。たとえば、poisson.rvs() はポアソン分布に従う乱数を生成します。

# SCIPY_POISSON

# In[1]

import numpy as np
from scipy.stats import poisson
import matplotlib.pyplot as plt

# μ=2のポアソン分布に従う乱数を1000個生成
# 乱数ジェネレータのシードを0に設定
k = poisson.rvs(2, size=1000, random_state=0)

# 階級値の設定
bins = np.arange(0, 10)-0.5

# FigureとAxes
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_title("Poisson distribution", fontsize=18, pad=8)
ax.set_xlabel("k", fontsize = 15, labelpad=6)
ax.set_ylabel("N(k)", fontsize = 15, labelpad=6)

# 各乱数の個数をヒストグラムで表示
ax.hist(k, bins, rwidth=0.6, color="mediumblue")

plt.show()

poisson.pmf() はポアソン分布の確率密度関数を返します。

# In[2]

# 確率変数
x = np.arange(0, 11)

# ポアソン分布(μ=2.5)
y = poisson.pmf(x, 2.5)

# FigureとAxes
fig = plt.figure()
ax = fig.add_subplot(111)
ax.grid()
ax.set_title("Poisson distribution", fontsize=18, pad=8)
ax.set_xlim(0, 10)
ax.set_xlabel("k", fontsize = 15, labelpad=8)
ax.set_ylabel("P(X=k)", fontsize = 15, labelpad=8)

# データをプロット
ax.plot(x, y, marker="o", color="orangered")

plt.show()

poisson.cdf(k, mu) は期待値 mu のポアソン分布の累積分布関数、すなわち事象発生回数が k 以下であるような確率分布を返します。

# In[3]

# 確率変数
k = np.arange(0, 11)

# ポアソン累積分布関数(μ=3)
cdf = poisson.cdf(k, 3)

# FigureとAxes
fig, ax = plt.subplots(1, 1)
ax.grid()
ax.set_title("Poisson Cumulative Distribution", fontsize=15, pad=8)
ax.set_xlim(0, 10)
ax.set_ylim(0, 1.2)
ax.set_xlabel("k", fontsize = 15, labelpad=8)
ax.set_ylabel("P(X<k)", fontsize = 15, labelpad=8)

# データをプロット
ax.plot(k, cdf, marker="o", color="red")

plt.show()

mean(), median(), var(), std() メソッドで、ポアソン分布の期待値、中央値、分散、標準偏差を取得できます。

# In[4]

# ポアソン分布の期待値を設定
mu = 2.5

# 期待値
mean = poisson.mean(mu)

# 中央値
median = poisson.median(mu)

# 分散
var = poisson.var(mu)

# 標準偏差
std = poisson.std(mu)

print("期待値: {:.3f}".format(mean))
print("中央値: {:.3f}".format(median))
print("分散: {:.3f}".format(var))
print("標準偏差: {:.3f}".format(std))

# 期待値: 2.500
# 中央値: 2.000
# 分散: 2.500
# 標準偏差: 1.581

stats() メソッドは期待値、分散、歪度、尖度の $4$ 種の統計量をまとめて取得します。

# In[5]

# μ=2.5のポアソン分布の期待値、分散、尖度、歪度
mean, var, skew, kurt = poisson.stats(mu, moments="mvsk")

print("期待値: {:.3f}".format(mean))
print("分散: {:.3f}".format(var))
print("尖度: {:.3f}".format(skew))
print("歪度: {:.3f}".format(kurt))

# 期待値: 2.500
# 分散: 2.500
# 尖度: 0.632
# 歪度: 0.400