統計のお勉強です。今日はポアソン分布 について。
ポアソン分布とは ポアソン分布 (Poisson distribution) は、「稀にしか発生しない事象が、一定時間内にどれくらい発生するか」 をモデル化するときに使える確率分布です。
こういった事例は世の中にいくらでもあって、例えば、
雨の降る日数 来客数 注文数 事故発生数 不良品の個数 サッカーの試合でのチーム別得点数 などがあります。
この確率分布の具体的な式は次のようになります。
$$
P(\textcolor{#ff0000}{x}) = \frac{e^{-\textcolor{#009900}{\lambda}} \cdot \textcolor{#009900}{\lambda}^\textcolor{#ff0000}{x}}{\textcolor{#ff0000}{x}!}
$$
\( \textcolor{#009900}{\lambda} \): 単位時間内の平均発生回数 (> 0)
この \( P(x) \) を計算すると、「単位時間内に事象が \( x \) 回発生する確率」 が分かってしまうんです。
不思議ですね！
ネイピア数 \( e \) が出てくるところがロマンチックです。
この式を使うには、単位時間あたりの平均発生回数を示す \( \lambda \gt 0 \) を与えてやる必要があります。
それが唯一の情報なので、それくらいは与えてやりましょう。
☝️ P(x) の読み方
上記の \( P(x) \) の式は、正確にはポアソン分布の 確率「質量」関数 と呼びます。
「質量」という用語は、二項分布やポアソン分布のように、離散型の確率変数を使う分布の場合に使います。
つまり、\( x \) の具体的な値が 0、1、2、3 という整数であるということです（平均回数を表す \( \lambda \) は整数でなくてもよいので混同しないよう注意）。
一方、正規分布の式のように、連続型の確率変数を使う分布の場合は、確率「密度」関数 と呼びます。
\( P(x) \) の代わりに \( P(X=k) \)（あるいは \( P(X=x) \)）と表現されていることもあります。
これは、回数を表す確率変数 \( X \) の具体的な値が \( k \) であることを示しているだけで意味は同じなので、\( k \) を \( x \) と読み替えれば OK です。
確率変数 \( X \) が上記のような確率分布を持つとき、確率変数 \( X \) はポアソン分布 \( Po(\lambda) \) に従う といい、\( X \sim Po(\lambda) \) と表現します。
ポアソン分布の例題 具体的な問題を見ると、ポアソン分布の使い方は簡単に分かります。
\( P(x) \) は「単位時間内に事象が \( x \) 回発生する確率」を求める関数ですが、この \( P(x) \) が扱おうとしている「単位時間」に対して、\( \lambda \)（単位時間あたりの平均発生回数）の方の「単位時間」を揃えてやるのがポイント です。
簡単な問題から順番に見ていきます。
問題 1 1 週間に平均 2 日雨が降るとき、1 週間に 3 日雨が降る確率は？
この問題では「単位時間」が 1 週間で揃っているので、\( \lambda \) の値としては 2（1 週間あたりの平均雨天日数）をそのまま使ってやれば OK です。
\( \lambda = 2 \) … （1 週間あたりの）平均雨天日数 \( x = 3 \) … （1 週間あたりの）雨天日数が 3 になる確率を求めたい $$
P(3) = \frac{e^{-2} \cdot 2^3}{3!} = \frac{1}{e^2} \cdot \frac{8}{6} \approx 0.1804
$$
よって、3 日雨が降る確率は、約 18.04 ％です。
同様に、\( P(0) \)、\( P(1) \)、\( P(2) \) を計算することで、雨がまったく降らない確率、1 日だけ降る確率、2 日降る確率を求めることができます。
雨の降る日数が 3 日以下である確率を求めたいときは、\( P(0) \) から \( P(3) \) までの確率を足し合わせる必要があります。
問題 2 1 時間あたり平均 2 人の客が来るとき、3 時間で 5 人の客が来る確率は？
このケースでは、求めるべき確率 \( P(5) \) の単位時間が 3 時間なので、\( \lambda \) の単位時間も 3 時間に合わせる 必要があります。
1 時間あたりの平均客数が 2 人なので、3 時間あたりの平均客数は \( \lambda = 6 \) 人となります。
\( \lambda = 6 \) … （3 時間あたりの）平均客数 \( x = 5 \) … （3 時間あたりの）客数が 5 になる確率を求めたい
$$
P(5) = \frac{e^{-6} \cdot 6^5}{5!} \approx 0.1606
$$ よって、5 人の客が来る確率は、約 16.06 ％です。
問題 3 不良品率が 0.1 ％のとき、1 日の生産が 1000 個で 1 日の不良品数が 2 個以下になる確率は？
この問題では、\( \lambda \) の値として、「1 日あたりの平均不良品数」を設定します（％ではなく数にする）。
1 日あたりの生産数が 1000 個で、不良品率が 0.1 ％ということは、平均不良品数は \( \lambda = 1000 \cdot 0.001 = 1 \) となります。
不良品数が 2 個以下になる確率を求めたいときは、\( P(0) \)、\( P(1) \)、\( P(2) \) を全て求めて足し合わせる必要があります。
ポアソン分布の式は便利ですが、このようにひとつずつ足し合わせないといけないのは面倒なところですね。
\( \lambda = 1 \) … （1 日あたりの）平均不良品数 \( x = 0, 1, 2 \) … （1 日あたりの）不良品数が 0 or 1 or 2 になる確率を求めたい $$
\begin{align}
P(0) &= \frac{e^{-1} \cdot 1^0}{0!} \approx 0.3679 \\
P(1) &= \frac{e^{-1} \cdot 1^1}{1!} \approx 0.3679 \\
P(2) &= \frac{e^{-1} \cdot 1^2}{2!} \approx 0.1839 \\
\end{align}
$$
$$
P(0) + P(1) + P(2) \approx 0.9197
$$
よって、不良品数が 2 個以下になる確率は、約 91.97 ％です。
☝️ ワンポイント
これまでの例から分かるように、ポアソン分布を使う場合は、\( \lambda \) の値は小さい値になります。
1 以下の小数点数（例: 0.23）になっても構いません。
問題 4 過去 5 年間に大きな地震が 7 回起こっています。これから 3 日の間に大きな地震が 1 回起こる確率は？
\( \lambda \) の値として、「3 日間の平均地震発生回数」を計算して代入するのがポイントです。
\( \lambda = 7 \times \frac{3}{365 \times 5} \approx 0.0115 \) … （3 日間の）平均発生回数 \( x = 1 \) … （3 日間の）発生回数が 1 回になる確率を求めたい $$
P(1) = \frac{e^{-0.0115} \cdot 0.0115^1}{1!} \approx 0.01137
$$
よって、地震が 1 回発生する確率は、約 1.137 ％です。
問題 5 日本の人口は 1 億 2737 万人とします。日本全体で、過去 5 年間に平均で 20204 人/年の食中毒が発生しました。千葉県松戸市（人口 484600 人）において食中毒にかかる人が 1 人も発生しない日の確率は？
数字がいっぱい出てきますが、落ち着いて考えればこれまでの問題と同じように計算できます。
「1 人も発生しない」というのは発生回数が 0 回ということなので、\( x = 0 \) です。
よって、確率 \( P(0) \) を求める問題です。
この確率は「1 日あたり」の確率なので、平均発生回数 \( \lambda \) の単位時間も 1 日で考えます (20204 / 365)。
さらに、この平均値は日本全体での値なので、松戸市内だけの値として人数の割合 (484600 / 127370000) を掛け合わせてやる必要があります。
\( \lambda = \frac{20204}{365} \times \frac{484600}{127370000} \approx 0.2106 \) … （松戸市での 1 日あたりの）平均発生回数 \( x = 0 \) … （松戸市での 1 日あたりの）発生回数が 0 になる確率を求めたい $$
P(0) = \frac{e^{-0.2106} \cdot 0.2106^0}{0!} = e^{-0.2106} \approx 0.8101
$$
よって、松戸市で 1 日に 1 人も食中毒が発生しない確率は、約 81.01 ％です。
問題 6 ある機械の単位時間あたりの修理回数を確率変数 \( X \) で表すと、\( X \) は平均回数 0.8 のポアソン分布に従います。運用コストとして \( 100 + 50 X^2 \) かかるとき、運用コストの期待値は？
この問題を解くには、ポアソン分布に従う確率変数 \( X \) の 期待値 \( \text{E}[X] \) と 分散 \( \text{V}[X] \) がどのような値になるかを知っている必要があります。
\( X \) の期待値は平均回数のことなので \( \lambda = 0.8 \) ですね。
結論から言うと、\( X \) の分散も \( \lambda \) になります（計算は省略）。
式で表現すると次のように書けます。
$$
\begin{align}
\text{E}[X] &= \textcolor{#ff0000}{\lambda} \\
\text{V}[X] &= \textcolor{#ff0000}{\lambda}
\end{align}
$$
平均回数と分散が等しくなるというのは、ポアソン分布の面白いところです。
この問題では、運用コストが \( 100 + 50 X^2 \) になると言っているので、求めるべき運用コストの期待値は \( \text{E}[100 + 50 X^2] \) と書けます。
あとはこれを計算するだけです。
$$
\begin{align}
\text{E}[100 + 50 X^2] &= 100 + 50 \text{E}[X^2] \\
&= 100 + 50 (\text{V}[X] + \text{E}[X]^2) \\
&= 100 + 50 (\lambda + \lambda^2) \\
&= 100 + 50 (0.8 + 0.8^2) \\
&= 172
\end{align}
$$
よって、運用コストの期待値は 172 になります。
☝️ E[X^2] の計算方法
上記の計算では、\( \text{E}[X^2] = \text{V}[X] + \text{E}[X]^2 \) という展開がポイントになります。
この式は、確率変数 \( X \) の分散の定義式から導かれます。
$$
\begin{align}
\textcolor{green}{\text{V}[X]}
&= \text{E} \left[ (X - \text{E}[X])^2 \right] \\
&= \text{E} \left[ X^2 - 2 \text{E}[X] \cdot X + \text{E}[X]^2 \right] \\
&= \text{E}[X^2] - 2\text{E}[X] \cdot \text{E}[X] + \text{E}[X]^2 \\
&= \textcolor{green}{\text{E}[X^2] - \text{E}[X]^2}
\end{align}
$$
\( \text{E}[\text{E}[X]] \) のように入れ子になっている部分が若干ややこしいですが、内側の \( \text{E}[X] \) は定数扱いの期待値（平均値）なので、そのまま外側の \( \text{E}[~] \) を外すことができます。
2 という定数の期待値 \( \text{E}[2] \) が 2 であるのと同じことです。
上記のように、確率変数 \( X \) や期待値 \( \text{E} \) の記号を使うと混乱するという人は、次のように具体的な値 \( x_i \) と平均値 \( \bar{x} \) を使って記述すると理解しやすいかもしれません。
$$
\begin{align}
\textcolor{green}{\text{V}[X]}
&= \frac{1}{n} \sum (x_i - \bar{x})^2 \\
&= \frac{1}{n} \sum (x_i^2 - 2 \bar{x} \cdot x_i + \bar{x}^2) \\
&= \frac{1}{n} \sum x_i^2 - 2 \bar{x} \cdot \frac{1}{n} \sum x_i + \bar{x}^2 \\
&= \frac{1}{n} \sum x_i^2 - 2 \bar{x} \cdot \bar{x} + \bar{x}^2 \\
&= \frac{1}{n} \sum x_i^2 - \bar{x}^2 \\
&= \textcolor{green}{\text{E}[X^2] - \text{E}[X]^2}
\end{align}
$$
何度も出てくる \( \frac{1}{n} \sum x_i \) という部分は平均値を示しているので、期待値 \( \text{E}[X] \) と同じ意味です。
そう考えると、\( \text{E}[X] \) という記号を使った方がスッキリ書けますね。
（おまけ）Python や Excel でポアソン分布の計算 Python の場合 Python でポアソン分布の確率質量関数 poisson() を定義する例です。
def poisson ( * , lam , k ):
"""
Poisson distribution function.
Args:
lam (float): the average number of events in a unit of time.
k (int): the number of events.
Returns:
float: the probability of k events occurring in a unit of time.
"""
from math import exp , factorial
return exp ( - lam ) * lam ** k / factorial ( k )
# 使用例
print ( " {:.4f} " . format ( poisson ( lam = 2 , k = 3 ))) # 0.1804
print ( " {:.4f} " . format ( poisson ( lam = 6 , k = 5 ))) # 0.1606
あるいは、scipy.stats モジュールが提供する poisson() 関数を使う方法もあります。
サンプルコードでは、poisson() 関数の最初の引数を * とすることで、名前付き引数の形で呼び出すことを強制しています。
こうすることで、呼び出し時に引数の順番を間違えてしまうミスをなくせます。
Excel の場合 Excel では、セルに =POISSON.DIST(イベント数, 平均, 関数形式)
イベント数 … \( x \) のこと。単位時間にこの数だけ発生する確率を求めたい。平均 … \( \lambda \) のこと。単位時間の平均発生回数をセットする。関数形式 … TRUE なら \( x = 0～イベント数 \) の累積確率、FALSE なら \( x = イベント数 \) の確率。今回は FALSE をセットすれば OK。例えば、
=POISSON.DIST(3, 2, FALSE)
のように入力すると、平均回数 \( \lambda = 2 \) のときに \( x = 3 \) 回発生する確率 \( P(3) \) が計算され、0.180447044 と表示されます。
関連記事 
