前回 : Pythonで学ぶ！ベイズ統計第4回 - ma-hi-ro Hobby

Pythonで学ぶ！ベイズ統計第5回

今回は、確率分布でのベイズ更新の例を紹介する
前回紹介した、ポワソン分布、ガンマ分布を使用する

Pythonで学ぶ！ベイズ統計第5回

確率分布のベイズ更新

推定対象のパラメータが $\theta$
$\theta$ についての事前分布が $P_1(\theta)$
$data$ が得られた場合の、尤度分布が $P_2(data | \theta)$
事後分布（ $data$ が得られた場合の、 $\theta$ の更新分布）は、 $P_3(\theta | data)$

これらから、ベイズ更新は下記になる

$\displaystyle P_3(\theta | data) = \frac{P_2(data | \theta) P_1(\theta)}{\int P_2(data | \theta) P_1(\theta) d \theta}$

事後分布の計算

交差点の年間の事故発生回数の予測（確率分布）を取得する

推定対象は、下記のポワソン分布の平均値 $\lambda$
事故発生回数の分布は、ポワソン分布 $f(x | \lambda)$ を仮定（尤度分布として使用）
事前分布は、ガンマ関数 $g(\lambda | \alpha, \beta)$ を仮定
事前分布は、例えば、「全国の交差点の $\lambda$ の分布」とかが分かれば定められる
得られたデータを $x$ とする（とある年の事故発生回数）

事後分布を計算すると

$\displaystyle \begin{aligned} P_3(\lambda | x) &= \frac{P_2(x | \lambda) P_1(\lambda)}{\int P_2(x | \lambda) P_1(\lambda) d \lambda} \\ &\propto f(x | \lambda) g(\lambda | \alpha, \beta) \\ &\propto \frac{\lambda ^x e^{-\lambda}}{x!} \frac{\beta ^{\alpha} \lambda^{\alpha - 1} e^{- \beta \lambda}}{\Gamma(\alpha)} \\ &\propto \lambda^{\alpha + x - 1} e^{- (\beta + 1) \lambda} \end{aligned}$

となる（途中から $\lambda$ に関係する部分のみ抽出）

上記から、事後分布もガンマ分布となり、各パラメータは、下記になる

$\displaystyle \begin{aligned} \alpha_1 &= \alpha_0 + x \\ \beta_1 &= \beta_0 + 1 \end{aligned}$

同様に、データが複数個得られた場合は、下記になる

$\displaystyle \begin{aligned} \alpha_1 &= \alpha_0 + \sum_i^n x_i \\ \beta_1 &= \beta_0 + n \end{aligned}$

具体例

事前分布のパラメータは、 $\alpha = 1.6, \beta = 1.0$ とし
データとして $x=4$ が得られた場合を考える

この場合、事後分布のパラメータは、

$\displaystyle \begin{aligned} \alpha_1 &= \alpha_0 + x &= 1.6 + 4 = 5.6 \\ \beta_1 &= \beta_0 + 1 &= 1.0 + 1 = 2.0 \end{aligned}$

平均値と分散は、下記になる

$\displaystyle \begin{aligned} \frac{\alpha_1}{\beta_1} &= 2.8 \\ \frac{\alpha_1}{\beta_1^2} &= 1.4 \\ \end{aligned}$

事前分布では、平均1.6、分散も1.6だった
追加されたデータが4のため平均値は大きくなり、データが追加されたことで分散は小さくなっている

さらに、過去のデータとして、5年分 $x = 3, 1, 2, 1, 5$ が得られた場合は

$\displaystyle \begin{aligned} \alpha_2 &= \alpha_1 + \sum_i^n x_i &= 5.6 + 12 &= 17.6 \\ \beta_2 &= \beta_1 + n &= 2.0 + 5 &= 7.0 \end{aligned}$

平均値と分散は、下記になる

$\displaystyle \begin{aligned} \frac{\alpha_2}{\beta_2} &\sim 2.51 \\ \frac{\alpha_2}{\beta_2^2} &\sim 0.36 \\ \end{aligned}$

分散がさらに小さくなり、推測が進んでいることがわかる

Pythonによるグラフ化

グラフで表現すると下図になる

import matplotlib.pyplot as plt
import numpy as np

from scipy import stats

x = np.arange(0, 10, step=0.01)

fig = plt.figure(figsize=(6, 5))
ax = fig.add_subplot(111)
data = stats.gamma.pdf(x, 1.6, scale=1.0/1.0)
ax.plot(x, data, '-', label=r'$\alpha =1.6, \beta =1.0$')
data = stats.gamma.pdf(x, 5.6, scale=1.0/2.0)
ax.plot(x, data, '-', label=r'$\alpha =5.6, \beta =2.0$')
data = stats.gamma.pdf(x, 17.6, scale=1.0/7.0)
ax.plot(x, data, '-', label=r'$\alpha =17.6, \beta =7.0$')
ax.set_xlim(-0.5, 9.5)
ax.grid()
ax.legend()