11 ポワソン分布

Author

Ryo Nakagami

Published

2024-09-25

Modified

2024-10-19

ポワソン分布の性質

稀な現象の大量観測（二項分布に置き換えるならば \(n\) がおおきく，\(p\) が十分小さい）において発生する個数の分布を表すのにポワソン分布が用いられます．例として，

ある都市の１日に起こる交通事故の件数の分布
ある都市で１年間に肺がんによりなくなる人数の分布

Def: ポワソン分布(Poisson distribution)

確率変数 \(X\) がパラメーター \(\lambda > 0\) のポワソン分布に従うとき，標本空間は \(\mathcal{X} = \{0, 1, \cdots, n\}\)，確率関数 \(f_X(x)\) は

\[ f_X(x) = \bigg\{\begin{array}{c} \frac{\lambda^x}{x!}\exp(-\lambda) & x \in \mathcal{X}\\0 & \text{otherwise}\end{array} \]

このとき，\(X\sim \operatorname{Po}(\lambda)\) と表す．

ポワソン分布の \(\lambda > 0\) は強度もしくは生起率と呼ばれるパラメーターです．

▶ 確率関数の和

\[ \begin{align*} \sum_{x=0}^\infty f_X(x) &= \sum_{x=0}^\infty\frac{\lambda^x}{x!}\exp(-\lambda)\\ &= \exp(-\lambda)\sum_{x=0}^\infty\frac{\lambda^x}{x!}\\ &= \exp(-\lambda)\left[1 + \frac{\lambda}{1!} + \frac{\lambda^2}{2!} + \cdots\right]\\ &= \exp(-\lambda)\exp(\lambda) \quad\because\text{マクローリン展開より}\\ &= 1 \end{align*} \]

▶ 期待値の導出

確率変数 \(X\sim\operatorname{Po}(\lambda)\) について，

\[ \begin{align*} \mathbb E[X] &= \sum_{x=0}^\infty x \frac{\lambda^x}{x!}\exp(-\lambda)\\ &= \lambda\exp(-\lambda)\sum_{x=1}^\infty\frac{\lambda^{x-1}}{(x-1)!}\\ &= \lambda\exp(-\lambda)\sum_{k=0}^\infty\frac{\lambda^k}{k!}\\ &= \lambda\exp(-\lambda)\exp(\lambda)\\ &= \lambda \end{align*} \]

▶ 分散の導出

\[ \begin{align*} \mathbb E[X(X-1)] &= \sum_{x=0}^\infty x(x-1) \frac{\lambda^x}{x!}\exp(-\lambda)\\ &= \lambda^2\exp(-\lambda)\sum_{x=2}^\infty\frac{\lambda^{x-2}}{(x-2)!}\\ &= \lambda^2 \end{align*} \]

従って，

\[ \begin{align*} \operatorname{Var}(X) &= \mathbb E[X(X- 1)] + \mathbb E[X](1 - \mathbb E[X])\\ &= \lambda^2 + \lambda(1 - \lambda)\\ &= \lambda \end{align*} \]

▶ 確率母関数の計算

確率変数 \(X\sim\operatorname{Po}(\lambda)\) としたとき，確率母関数 \(G_X(s)\) は

\[ \begin{align*} G_X(s) &= \sum_{x=0}^\infty s^x\frac{\lambda^x}{x!}\exp(-\lambda)\\ &= \exp(-\lambda)\sum_{x=0}^\infty \frac{(s\lambda)^x}{x!}\\ &= \exp(-\lambda)\exp(s\lambda)\\ &= \exp((s-1)\lambda) \end{align*} \]

期待値は

\[ \begin{align*} \mathbb E[X] &= G_X^\prime(1)\\ &= \lambda \exp((1-1)\lambda)\\ &= \lambda \end{align*} \]

分散は \(\operatorname{Var}(X) = \mathbb E[X(X-1)] + \mathbb E[X](1 - \mathbb E[X])\) なので

\[ \begin{align*} \mathbb E[X(X-1)] &= G_X^{\prime\prime}(1)\\ &= \lambda \lambda\exp((1-1)\lambda)\\ &= \lambda^2 \end{align*} \]

従って，

\[ \begin{align*} \operatorname{Var}(X) &= \lambda^2 + \lambda(1-\lambda)\\ &= \lambda \end{align*} \]

▶ 積率母関数の計算

確率変数 \(X\sim\operatorname{Po}(\lambda)\) としたとき，積率母関数 \(M_X(t)\) は

\[ \begin{align*} M_X(t) &= \mathbb E[\exp(tX)]\\ &= \sum_{x=0}^\infty\exp(-\lambda)\exp(tX)\frac{\lambda^x}{x!}\\ &= \exp(-\lambda)\sum_{x=0}^\infty\frac{(\exp(t)\lambda)^x}{x!}\\ &= \exp(-\lambda)\exp(e^t\lambda)\\[5pt] &= \exp[\lambda(e^t - 1)] \end{align*} \]

Example 11.1 : 19-20シーズンのMan Utdのゴール数とポワソン分布

Man Utd English Premier League 2019/20 fixture and resultsというサイトではプレミアリーグの各シーズン及び各チームの試合結果が保存されています. 上記サイトのURLの構造は

https://fixturedownload.com/results/epl-<シーズン>/<team名>

19-20シーズンのMan Utdのデータを抽出したい場合は次のようにパラメーターを設定します．

parameter	value	意味
team名	`man-utd`	Man Utd, すべて小文字、スペースはハイフンとなる
シーズン	`2019`	2019-2020シーズン

抽出結果は以下のようになります

Code

import pandas as pd
target_team = 'man-utd'
TARGET_TEAMNAME = target_team.replace('-', ' ').title()
URL_PATH = 'https://fixturedownload.com/results/epl-2019/'+target_team

df = pd.read_html(URL_PATH, flavor="bs4")[0]
df.head()

	Round Number	Date	Location	Home Team	Away Team	Result
0	1	11/08/2019 16:30	Old Trafford	Man Utd	Chelsea	4 - 0
1	2	19/08/2019 20:00	Molineux Stadium	Wolves	Man Utd	1 - 1
2	3	24/08/2019 15:00	Old Trafford	Man Utd	Crystal Palace	1 - 2
3	4	31/08/2019 12:30	St. Mary's Stadium	Southampton	Man Utd	1 - 1
4	5	14/09/2019 15:00	Old Trafford	Man Utd	Leicester	1 - 0

ここからMan Utdのゲームごとのゴール数配列を抽出し，ポワソン分布比較してみます．

Code

import numpy as np
from plotly import express as px
import plotly.graph_objects as go
from scipy.stats import poisson

from regmonkey_style import stylewizard as sw 

sw.set_templates("regmonkey_scatter")


def extract_goal(data, target_team):
    score_list = data["Result"].str.split(r"\s-\s", expand=True)
    score_list.columns = ["home_score", "away_score"]
    score_list = score_list.astype({"home_score": "int64", "away_score": "int64"})

    data = pd.concat([data, score_list], axis=1)
    return np.where(
        df["Home Team"] == target_team, data["home_score"], data["away_score"]
    )


goal_array = extract_goal(df.copy(), target_team=TARGET_TEAMNAME)

## 可視化
x, y = np.unique(goal_array, return_counts=True)

fig = go.Figure()
fig.add_trace(
    go.Bar(
        x=x,
        y=y / sum(y),  # Normalizing the frequency
        name="Goals per games", 
        width=0.3
    )
)

fig.add_trace(
    go.Scatter(
        x=x,
        y=poisson(np.mean(goal_array)).pmf(x),
        mode="lines+markers",
        name="Poisson",
        line=dict(color="gray", dash="dot"),
    )
)

fig.update_layout(
    title="19-20シーズン Man Utd Goal Distribution",
    xaxis_title="Goals",  # X-axis label
    yaxis_title="Frequency",  # Y-axis label
    legend_title="Legend",  # Legend title
    showlegend=True,  # Ensure legend is displayed
)

fig.show()

Fittingはそんなに悪くないが，scoreがゼロのところが理論頻度より実現値のほうが頻度が多い傾向がわかります
0が過剰に含まれたデータに対してはゼロ過剰ポアソン分布(Zero-inflated Poisson Model, ZIP)を用いて対処することが考えられます

ポワソン分布の再生性

Theorem 11.1 : ポワソン分布の再生性

互いに独立な確率変数 \(X\sim\operatorname{Po}(\lambda_x), Y\sim\operatorname{Po}(\lambda_y)\) について, \(X + Y\) もパラメーター \(\lambda_x + \lambda_y\) のポワソン分布に従う．つまり，

\[ X + Y \sim \operatorname{Po}(\lambda_x + \lambda_y) \]

Proof: PGFを用いた証明

\[ \begin{align*} G_{X+Y}(s) &= G_{X}(s)G_{Y}(s)\\ &= \exp(-\lambda_x)\exp(s\lambda_x)\exp(-\lambda_y)\exp(s\lambda_y)\\ &= \exp(-(\lambda_x+\lambda_y))\exp(s(\lambda_x+\lambda_y)) \end{align*} \]

これは，パラメーター \(\lambda_x + \lambda_y\) のポワソン分布のPGFに一致しています．

Proof: MGFを用いた証明

\[ \begin{align*} M_{X+Y}(t) &= M_{X}(t)M_{X}(t)\\ &= \exp[\lambda_x(e^t - 1)]\exp[\lambda_y(e^t - 1)]\\ &= \exp[\lambda_x(e^t - 1) + \lambda_y(e^t - 1)]\\ &= \exp[(\lambda_x + \lambda_y)(e^t - 1)] \end{align*} \]

これは，パラメーター \(\lambda_x + \lambda_y\) のポワソン分布のMGFに一致しています．

Proof

\(Z = X +Y\) としたとき，

\[ \begin{align*} P_{Z}(z) &= \sum_{j=0}^z\frac{\lambda_x^j}{j!}\exp(-\lambda_x)\frac{\lambda_y^{z-j}}{(z-j)!}\exp(-\lambda_y)\\ &= \exp(-\lambda_x-\lambda_y)\sum_{j=0}^z\frac{\lambda_x^j\lambda_y^{z-j}}{j!(z-j)!}\frac{z!}{z!}\\ &= \frac{1}{{z!}}\exp(-\lambda_x-\lambda_y)\underbrace{\sum_{j=0}^z\frac{z!\lambda_x^j\lambda_y^{z-j}}{j!(z-j)!}}_{=(\lambda_x + \lambda_y)^z}\\ &= \frac{(\lambda_x + \lambda_y)^z}{z!}\exp(-\lambda_x-\lambda_y) \end{align*} \]

📘 REMARKS

「互いに独立」という仮定が成立しないと，再生性の性質は成立しません．実際に

\[ \begin{align*} X_1 &\sim \operatorname{Po}(\lambda)\\ X_1 &= X_2 \end{align*} \]

とした場合，\(Y = X_1 + X_2\) という確率変数の確率関数は

\[ \begin{align*} f_Y(y) = \left\{\begin{array}{cc} \frac{\lambda^{y/2}}{(y/2)!}\exp(-\lambda) & y \text{: even}\\ 0 & y \text{: odd} \end{array}\right. \end{align*} \]

となり，\(Y\) はポワソン分布に従う訳ではないことが確認できます．

他の確率分布との関係

二項分布との関係

確率変数 \(X\) について

\[ X \sim \operatorname{Binom}(n, p) \]

を考えます．このとき，\(p\) は十分小さく

\[ \lim_{n\to\infty} np = \lambda \]

が成立するとします．\(X\) についての確率関数を求めると

\[ \begin{align*} P(X = x) &= {}_nC_x p^x(1-p)^{n-x}\\ &= \frac{n(n-1)\cdots(n-x+1)}{x!}\left(\frac{\lambda}{n}\right)^x\left(1-\frac{\lambda}{n}\right)^{n-x}\\ &\to \frac{\lambda^x}{x!}\exp(-\lambda) \quad (\text{as } n\to \infty) \end{align*} \]

Estimation

Ommited variable and overdispersion

Def: over-dispersion

観察データについて，期待される分散よりも大きい分散(variation)が確認される状があるとき，overdispersionであるという．

ポワソン分布において，上で確認したように分散と平均は一致するはずですが，ommitted variable biasを原因として overdispersionが発生するケースが多くあります．

\[ \begin{align*} y &\sim \operatorname{Po}(\lambda)\\ \log(\lambda) &= \beta_0 + \beta_1 x_1 + \beta_2 x_2 \end{align*} \]

という確率変数を考えます．このとき次のような条件を想定します：

観察データにおいて, \(x_2\) は観察不能(=missing variable)
\(x_1 \perp x_2\)

このとき，\(y\) について \(x_1\) を条件づけたときの期待値を計算すると

\[ \begin{align*} \mathbb E[y\vert x_1] &= \mathbb E[\mathbb E[y\vert x_1, x_2]\vert x_1]\\ &= \mathbb E[\exp(\beta_0 + \beta_1 x_1 + \beta_2 x_2)\vert x_1]\\ &= \exp(\beta_0 + \beta_1 x_1)\mathbb E[\exp(\beta_2 x_2)\vert x_1]\\ &= \exp(\beta_0 + \beta_1 x_1)\mathbb E[\exp(\beta_2 x_2)] \because\text{独立性}\\ &= \exp(\tilde\beta_0+ \beta_1 x_1) \end{align*} \]

と，切片には影響を与えますが，\(\beta_1\) についてはunbiasedに推定できることがわかります．一方，条件付き分散を見てみると

\[ \begin{align*} \operatorname{Var}(y\vert x_1) =& \mathbb E[\operatorname{Var}(y\vert x_1, x_2)\vert x_1] + \operatorname{Var}(\mathbb E[y\vert x_1, x_2]\vert x_1)\\[3pt] =& \mathbb E[\exp(\beta_0 + \beta_1 x_1 + \beta_2 x_2)\vert x_1] \\ &+ \operatorname{Var}(\exp(\beta_0 + \beta_1 x_1 + \beta_2 x_2)\vert x_1)\\ =& \exp(\tilde\beta_0+ \beta_1 x_1) \\ & + [\exp(\beta_0 + \beta_1x_1)]^2\operatorname{Var}(\exp(\beta_2 x_2))\\ >& \exp(\tilde\beta_0+ \beta_1 x_1) = \mathbb E[y\vert x_1] \end{align*} \]

以上より，OVBのとき，ovberdispersionが発生してしまうことがわかります．

📘 REMARKS

ポワソン回帰においてover-dispersionが発生した場合は以下のケースが考えられます：

重要な特徴量が欠落変数(ommited variable)になってしまっている
説明変数，被説明変数(response variable)についてmeasurement errorが発生してしまっている
\(\log(\lambda)\) と特徴量ベクトル \(\mathbf x\) のモデル特定に誤りがある
Outlierの存在
response variableが複数の確率分布の混合(mixture)に基づいている場合