ガンマポアソンと指数

🎓 レベル：標準　|　重要度：A（必須）

📎 前提：ベータ二項モデル　|　数理：共役事前分布（統計）・ポアソン分布（統計）・指数分布・ガンマ分布・ベータ分布（統計）

要点（BLUF）

発生率 $\lambda$ の共役事前はガンマ分布。ポアソン尤度でも指数尤度でも事後は再びガンマになります。「比率の事前＝ベータ、カウント・待ち時間の事前＝ガンマ」と覚えます。
更新則は一見ちがって見えますが、 $\alpha$ には総発生数、 $\beta$ には総観測量（期間・時間）を足すと統一的に読めます。事後平均 $\alpha/\beta=$ 発生数 / 観測量。
ポアソンの事後予測は負の二項分布になり、 $\lambda$ の不確実性ぶん過分散（分散＞平均）が生じます。点推定のポアソン予測は分散を過小評価します。

1. ガンマ–ポアソン：カウントの共役

単位期間あたりの平均発生回数 $\lambda>0$ を推定します。観測 $x_1,\dots,x_n$ が独立にポアソン分布に従うとします（ポアソン分布）。

尤度のカーネル（ $\prod x_i!$ は $\lambda$ を含まない定数）：

p(x\mid\lambda)=\prod_{i=1}^{n}\frac{\lambda^{x_i}e^{-\lambda}}{x_i!}\;\propto\;\lambda^{\sum x_i}\,e^{-n\lambda}

事前のカーネル（ガンマ、レート表記 $\beta$ ）：

p(\lambda)=\frac{\beta^\alpha}{\Gamma(\alpha)}\lambda^{\alpha-1}e^{-\beta\lambda}\;\propto\;\lambda^{\alpha-1}e^{-\beta\lambda}

事後＝積。肩の指数を足すだけ：

p(\lambda\mid x)\propto\lambda^{\sum x_i}e^{-n\lambda}\cdot\lambda^{\alpha-1}e^{-\beta\lambda}=\lambda^{(\alpha+\sum x_i)-1}e^{-(\beta+n)\lambda}

これは $\mathrm{Gamma}(\alpha+\sum x_i,\ \beta+n)$ のカーネル。よって

\lambda\sim\mathrm{Gamma}(\alpha,\beta),\ x_i\sim\mathrm{Poisson}(\lambda)\ \Longrightarrow\ \lambda\mid x\sim\mathrm{Gamma}\!\Big(\alpha+\textstyle\sum x_i,\ \beta+n\Big)

import numpy as np
from scipy import stats
from scipy.integrate import trapezoid          # numpy 2.0+ で np.trapz は廃止

counts = np.array([3, 1, 4, 2, 3, 0, 2, 5, 1, 2])   # n=10 日ぶんの発生数
n = len(counts); S = counts.sum()
alpha, beta = 2.0, 1.0                                # 事前 Gamma(2,1)（レート表記）
a_post, b_post = alpha + S, beta + n                  # 事後 Gamma(α+Σx, β+n)

post = stats.gamma(a=a_post, scale=1/b_post)          # scipy は scale=1/rate
lo, hi = post.interval(0.95)
print(f"Σx={S}, n={n}  →  事後 Gamma({a_post:.0f}, {b_post:.0f})")
print(f"事前平均(発生率)={alpha/beta:.3f}  事後平均={post.mean():.4f} (=25/11)")
print(f"95%信用区間 = [{lo:.3f}, {hi:.3f}]")

# 閉じた式 vs 数値積分で事後一致を確認
lam = np.linspace(1e-4, 8, 4000)
like = np.exp(-n*lam) * lam**S                        # ポアソン尤度カーネル
pri  = stats.gamma(a=alpha, scale=1/beta).pdf(lam)
post_num = (like*pri) / trapezoid(like*pri, lam)
print(f"閉じた式 vs 数値積分 の最大差 = {np.max(np.abs(post_num - post.pdf(lam))):.2e}")

出力：

Σx=23, n=10  →  事後 Gamma(25, 11)
事前平均(発生率)=2.000  事後平均=2.2727 (=25/11)
95%信用区間 = [1.471, 3.246]
閉じた式 vs 数値積分 の最大差 = 1.31e-14

出力の意味：事前 $\mathrm{Gamma}(2,1)$ （事前平均の発生率 $2.0$ ）に、10日で計23件のデータが入って、事後 $\mathrm{Gamma}(25,11)$ 、事後平均 $25/11=2.27$ 件/日。事前平均 $2.0$ とデータの平均 $2.3$ の間に落ち着きます。閉じた式とグリッド数値積分の差は $10^{-14}$ で、ベータ–二項と同じく完全一致です。

2. ガンマ–指数：待ち時間の共役

今度は事象の待ち時間（指数分布、レート $\lambda$ ）からレートを推定します。指数分布の尤度は

p(x\mid\lambda)=\prod_{i=1}^{n}\lambda e^{-\lambda x_i}=\lambda^{n}e^{-\lambda\sum x_i}\;\propto\;\lambda^{n}e^{-\lambda\sum x_i}

同じガンマ事前 $\propto\lambda^{\alpha-1}e^{-\beta\lambda}$ を掛けると、

p(\lambda\mid x)\propto\lambda^{n}e^{-\lambda\sum x_i}\cdot\lambda^{\alpha-1}e^{-\beta\lambda}=\lambda^{(\alpha+n)-1}e^{-(\beta+\sum x_i)\lambda}

すなわち $\mathrm{Gamma}(\alpha+n,\ \beta+\sum x_i)$ 。ポアソンとは足し先が入れ替わっています——指数では $n$ が形状 $\alpha$ に、 $\sum x_i$ がレート $\beta$ に足されます。

import numpy as np
from scipy import stats

waits = np.array([0.8, 2.1, 0.5, 1.7, 0.9, 1.3, 2.4, 0.3])  # n=8 の待ち時間
n = len(waits); T = waits.sum()
alpha, beta = 2.0, 2.0                                  # 事前 Gamma(2,2)
a_post, b_post = alpha + n, beta + T                    # 事後 Gamma(α+n, β+Σx)

post = stats.gamma(a=a_post, scale=1/b_post)
print(f"Σx={T:.1f}, n={n}  →  事後 Gamma({a_post:.0f}, {b_post:.1f})")
print(f"事後平均(レート λ) = {post.mean():.4f}")
print(f"平均待ち時間 1/λ の事後平均 = {b_post/(a_post-1):.4f}")  # E[1/λ]=β/(α-1)

出力：

Σx=10.0, n=8  →  事後 Gamma(10, 12.0)
事後平均(レート λ) = 0.8333
平均待ち時間 1/λ の事後平均 = 1.3333

出力の意味：8件の事象を観測するのに計 $10.0$ かかったので、事後 $\mathrm{Gamma}(10,12)$ 、発生レートの事後平均 $\lambda=10/12=0.833$ 件/単位時間。逆数の平均待ち時間は $E[1/\lambda]=\beta/(\alpha-1)=12/9=1.33$ 単位時間（ $1/\lambda$ は逆ガンマ分布に従い、その平均が $\beta/(\alpha-1)$ ）。

なお、指数尤度の事後予測（次の待ち時間 $x'$ ）は $\lambda$ を事後ガンマで積分するとロマックス分布（パレート第II種）になります。 $p(x'\mid D)=\int\mathrm{Exp}(x'\mid\lambda)\mathrm{Gamma}(\lambda\mid\alpha',\beta')d\lambda$ の結果で、指数分布より裾が重い（まれに非常に長い待ち時間が出る）。ポアソンの予測が過分散の負の二項になったのと同じ理屈で、レートの不確実性が予測の裾を伸ばします。

3. 統一的な読み方：α＝総発生数、β＝総観測量

ポアソンとガンマで足し先が入れ替わって見えるのは、「1観測が何を運ぶか」が違うだけです。どちらも

\boxed{\ \alpha_{\text{post}}=\alpha+(\text{観測した総発生数}),\qquad \beta_{\text{post}}=\beta+(\text{観測した総観測量・時間})\ }

と読めます。

	ポアソン（カウント）	指数（待ち時間）
1観測 $x_i$ が運ぶもの	発生数 $x_i$ 件・期間 1	発生 1 件・時間 $x_i$
総発生数 → $\alpha$ に加算	$\sum x_i$	$n$
総観測量 → $\beta$ に加算	$n$	$\sum x_i$
事後平均 $\alpha/\beta$	$\dfrac{\alpha+\sum x_i}{\beta+n}$	$\dfrac{\alpha+n}{\beta+\sum x_i}$

要するに発生率＝総発生数 / 総観測量であり、ハイパーパラメータ $\alpha,\beta$ は「事前に見込んだ仮想の発生数・観測量」（擬似観測）。 $\alpha/\beta$ が事前平均レート、 $\beta$ が大きいほど事前が強い、とベータ二項モデルと同じ構図で読めます。

4. 事後予測は負の二項分布：過分散が生まれる

ポアソンの事後予測（次の期間の発生数 $x'$ ）は、 $\lambda$ を事後で積分すると閉じた形の負の二項分布になります。

p(x'\mid D)=\int_0^\infty \mathrm{Poisson}(x'\mid\lambda)\,\mathrm{Gamma}(\lambda\mid\alpha',\beta')\,d\lambda=\mathrm{NB}\!\left(x'\,\middle|\,r=\alpha',\ p=\tfrac{\beta'}{\beta'+1}\right)

ポアソンは「分散＝平均」（等分散）ですが、 $\lambda$ がまだ不確かなぶん、予測には余分なばらつきが乗って過分散（分散＞平均）になります。これは現実のカウントデータがしばしばポアソンより散らばる理由の一つです。

import numpy as np
from scipy import stats

a_post, b_post = 25, 11                                # §1 の事後 Gamma(25,11)
# 事後予測 = 負の二項 NB(r=α', p=β'/(β'+1))
pp   = stats.nbinom(n=a_post, p=b_post/(b_post + 1))
# プラグイン：λ を事後平均に固定した素朴なポアソン予測
plug = stats.poisson(a_post / b_post)
print(f"事後予測 NB:      平均={pp.mean():.3f}, 分散={pp.var():.3f}")
print(f"プラグイン Poisson: 平均={plug.mean():.3f}, 分散={plug.var():.3f}")

出力：

事後予測 NB:      平均=2.273, 分散=2.479
プラグイン Poisson: 平均=2.273, 分散=2.273

出力の意味：平均はどちらも $2.27$ ですが、事後予測の分散 $2.48$ は、プラグインのポアソン分散 $2.27$ （＝平均、等分散）より大きい。差の $0.21$ が「 $\lambda$ をまだ確定できていない不確実性」のぶんです。点推定のポアソンで予測区間を作ると、この過分散を取りこぼして過信します（信用区間と事後予測分布と同じ構図がカウントでも起こる）。

⚠️ よくある誤解

「ポアソンも指数も更新は同じ」ではない：足し先が逆です。ポアソンは $\sum x_i\to\alpha,\ n\to\beta$ 、指数は $n\to\alpha,\ \sum x_i\to\beta$ 。「 $\alpha$ ＝総発生数、 $\beta$ ＝総観測量」で覚えれば取り違えません。
レート $\beta$ とスケール $1/\beta$ の取り違え：ガンマには「レート $\beta$ 」と「スケール $1/\beta$ 」の2流儀。更新則 $\beta\to\beta+n$ はレート表記。SciPy の stats.gamma は scale（＝ $1/\beta$ ）を取るので、scale=1/rate と渡します。
「事後予測はポアソンのまま」ではない： $\lambda$ の不確実性を積分すると負の二項になり、過分散します。等分散のポアソンで予測区間を作るのは過信。
「 $\lambda$ が事後平均なら待ち時間も逆数でよい」ではない： $1/\lambda$ の事後平均は $1/E[\lambda]$ ではなく $E[1/\lambda]=\beta/(\alpha-1)$ 。非線形変換の平均は取り違えやすい点です。