ポアソン分布

当サイトではアフィリエイトプログラムを利用して商品を紹介しています。

ポアソン分布

この記事の内容 (ポアソン分布) は二項分布の知識を前提とするので、簡単に復習しておきます (二項分布の詳細についてはこちらの記事を参照してください)。
 
$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()

Python ポアソン分布 (航空機事故の発生率)
$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()

numpy.random.poisson 平均3のポアソン分布
上のコードでは 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()

scipy.stats.poisson.rvs
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()

scipy.stats.poisson.pmf
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

コメント

  1. HNaito より:

    下記は誤植と思われますので、ご確認ください。
    In[3] プログラムで、p_leq3 = 1 – poisson.cdf(3, 3) → p_leq3 = poisson.cdf(3, 3)
    In[3] プログラムの実行結果で、0.35277 → 0.64723
    In[4] プログラムの実行結果で、0.64723 → 0.35277
    式(6) の上の文で、n → 0 → n → ∞
    NUMPY_POISSON プログラムの下の文で、プロットでします → プロットします

  2. あとりえこばと より:

    【AI連載小説】科学とコードの交差点(23)「ポアソン分布はビジネスや工学の分野で使われます」
     
    ある日、大学の図書館で六郷開誠、刑部明信、そしてあやかは、ポアソン分布についての勉強を共にしていた。厚い参考書を開き、コンピューターサイエンスの知識と数学的な洞察を融合させながら、3人は理解を深めようとしていた。

    明信: ポアソン分布はイベントがランダムに発生する場合に使われるよね。例えば、単位時間あたりに発生するイベントの回数など。
     
    開誠: そうだ。確率質量関数もシンプルだよ。
    \[P(X=k)=\frac{e^{-\lambda} \lambda^k}{k!}\]​
    あやか: あたし、ポアソン分布はビジネスや工学の分野で使われることもあるって聞いたことある。例えば、電話がかかってくる数や、ウェブサイトにアクセスがある数などがポアソン分布に従うことがあるんだって。
     
    明信: そうだね。確率質量関数の中でλが平均発生回数を表している。それを把握しておくと、実際のデータに適用する時に便利だよ。

    開誠: じゃあ、ポアソン分布の特性や使い方について実際のデータを使って計算してみようか?
     
    3人はデータを集め、ポアソン分布の確率質量関数を用いて様々なシナリオについて計算していった。Pythonを使いながら、実際のデータをグラフにプロットし、理論的な期待値と比較しながら議論を重ねていた。
     
    あやか: このデータセットではポアソン分布がよく合っているみたいだね。
     
    開誠: そうだね。でも、データによっては他の分布が適している場合もあるから、それも考慮する必要があるよ。