単位時間の間に発生する事象の平均回数を$\lambda$とした時、ポアソン分布を用いることで、一定期間内に事象が起こる回数を確率的に表現することができます。
今回はこのポアソン分布に対してベイズ推定のアプローチから、一定期間内に事象が発生する確率（事後確率）を導出していきたいと思います。
ポアソン分布に対しては、そのパラメータの事前分布に共役事前分布であるガンマ分布を用いることで、解析的にポアソン分布のベイズ学習とベイズ予測を行うことができます。
- ポアソン分布と共役事前分布であるガンマ分布を用いて解析的にベイズ推定を行う
今回扱いたい問題設定
まず、実際の実務を想定して、実際に予測したい問題を定式化します。
今回は、具体例を通してベイズ推定を学んでいきます。
北海道は秋になると鮭が遡上するため、8月から9月頃にかけて鮭釣りを楽しむことができます。
あなたは東京に住んでいるため、過去の北海道の鮭釣りの釣れた記録から、今年どれだけの鮭が釣れるかを予測するため、ポアソン分布を用いようと考えています。
過去に鮭釣りのために、北海道には20日間滞在しており、それぞれ、[0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 2, 1, 1, 1, 2, 3, 1 ,2, 0, 0, 0, 2, 1, 0, 2, 1, 0] 釣れたというデータがあります。
このデータを用いて、ポアソン分布を用いて、次回釣りにいった時、1日で釣れる鮭の匹数を予測したいと考えています。
過去に釣ることができた鮭の数から、今後鮭を釣れる数を予測する問題です。
今回の問題のように、一定期間内にイベントがある一定数起きるような事象は、ポアソン分布を用いて表現することができます。
今回はベイズ的な手法を用いて、データを用いてモデルを学習し、学習したモデルを利用して、ベイズ予測を行います。
ベイズ学習を行う
まず、ベイズ学習をするためにモデル化をします。今回、予測したいのは、北海道に出向いた時に、1日あたり鮭を何匹釣れることができるのか、です。
このように、一定期間の間に、非負(ゼロ以上の整数)回の事象を予測するには、ポアソン分布を用います。ポアソン分布は、パラメータ$\lambda$を用いて下記のように定義されています。
\begin{equation}
\begin{split}
Poi(x| \lambda) = \frac{\lambda^{x}}{x !} e^{- \lambda}
\end{split}
\end{equation}これが、ベイズ推定における尤度関数に相当する式です。ここで、ベイズ推定では、(1)に登場するパラメータも確率的に表すので、$\lambda$も確率的に表現してみます。
(1)式のポアソン分布のパラメータ$\lambda$は、0以上の実数の範囲を取ります。
よって、0以上の実数の範囲を表すことができ、ポアソン分布の共役事前分布であるガンマ分布を用いてこのパラメータ$lambda$を確率的に表現してみます。
\begin{equation}
\begin{split}
p(\lambda | a, b) &= Gam(\lambda | a, b) \\
&= \frac{b^a}{\Gamma(a)} \lambda^{a-1}e^{-b \lambda}
\end{split}
\end{equation}ガンマ分布を用いて、(2)式のようにポアソン分布のパラメータを表現することができました。
ここで、ベイズ推定の学習では、データ$\bm{X}$を元に、$\lambda$の事後分布 $p(\lambda | \bm{X})$を求めます。これは条件付き確率の式をそのまま用いると、下記のように表現できます。
\begin{equation}
\begin{split}
Poi(\lambda | \bm{X}) &= \frac{p(\bm{X} | \lambda) p(\lambda)}{p(\bm{X})} \\
& \propto
p(\bm{X} | \lambda) p(\lambda) \\
&=
\prod_{n=1}^{N} p(x_n | \lambda) p(\lambda) \\
&=
\prod_{n=1}^{N} Poi(x_n | \lambda) Gam(\lambda | a, b)
\end{split}
\end{equation}ここで、(2)式目の段階で$\propto$ をやっていますが、これは分母の$p(\bm{X})$は$\lambda$に関係ないので考えないようにしています。
今は$\lambda$の大まかな形状について考えているので、$\lambda$に関係ない項は考えないようにしています。現時点でこれが理解できなくても、今回の事例でベイズ学習とベイズ推論を考える上で、このように考えて問題ないことは、記事を最後まで見れば納得できると思います。今はそういうものか、という感じで考えてください。
(3)式は積の形になって変形しにくいので、(3)式の対数を考えます。
\begin{equation}
\begin{split}
ln Poi(\lambda | \bm{X}) &=
ln \prod_{n=1}^{N} Poi(x_n | \lambda) Gam(\lambda | a, b) \\
&=
\sum_{n=1}^{N} ln Poi(x_n | \lambda) + ln Gam(\lambda | a, b) + const. \\
&=
\sum_{n=1}^{N} \{ x_n ln \lambda - ln x_n! - \lambda \} + (a-1)ln\lambda -b\lambda + const. \\
&=
(\sum_{n=1}^{N} x_n +a -1)ln \lambda - (N-b)\lambda + const.
\end{split}
\end{equation}上の式変形では、$\lambda$に関係ない項は、全て、const.の定数の項に押し込めています。
ここで(4)式の事後分布$Poi(\lambda | \bm{X})$ですが、これはガンマ分布と同じ形をしているのがわかります。実際に、(2)式のガンマ分布も対数をとってみると、
\begin{equation}
\begin{split}
ln Gam(\lambda | a, b) &= \frac{b^a}{\Gamma(a)} \lambda^{a-1}e^{-b \lambda} \\
&= (a-1)ln\lambda -b\lambda
\end{split}
\end{equation}ここで、(4)式と(5)式を比較すると、新しい事後分布のガンマ分布のパラメタを$\hat{a},\hat{b}$とすると、
\begin{equation}
\begin{split}
\hat{a} &= a + \sum_{n=1}^{N} x_n \\
\hat{b} &= b + N
\end{split}
\end{equation}となって、ガンマ関数のパラメータを更新することができました。
よって、データによって学習した事後分布 $Poi(\lambda | \bm{X})$は、下記のようになります。
\begin{equation}
\begin{split}
Poi(\lambda | \bm{X}) &= Gam(\lambda | \hat{a}, \hat{b})
\end{split}
\end{equation}これが、ベイズ学習によって得られた事後分布になります。ベイズ推定の学習では、パラメータの事後分布が得ることができれば学習は完了です。続いて、データによって学習した事後分布を用いて、ベイズ予測をしてみましょう。
事後分布をもとにベイズ予測を行う
ベイズ予測は、ベイズ学習によって得られた事後分布をもとに、未知の値$\hat{x}$を予測することです。ベイズ予測を一般化した式を下記に記載します。
\begin{equation}
\begin{split}
p(\hat{x} | \bm{X}) = \int p(\hat{x}| \theta)p(\theta | \bm{X})d\theta
\end{split}
\end{equation}(8)式をもとに今回の場合で考えると、
\begin{equation}
\begin{split}
p(\hat{x} | \bm{X}) &= \int p(\hat{x}| \lambda)p(\lambda | \bm{X})d\lambda \\
&=
\int \frac{\lambda^{\hat{x}}}{\hat{x} !} e^{- \lambda} Gam(\lambda | \hat{a}, \hat{b}) d \lambda \\
&=
\int \frac{\lambda^{\hat{x}}}{\hat{x} !} e^{- \lambda} \frac{ \hat{b}^{\hat{a}}}{\Gamma(\hat{a})} \lambda^{\hat{a}-1}e^{-\hat{b} \lambda} d\lambda \\
&=
\frac{\hat{b}^{\hat{a}}}{\hat{x} ! \Gamma(\hat{a})}
\int
\frac{}{} \lambda^{\hat{x} + \hat{a} -1} e^{-(1+\hat{b})\lambda}d \lambda \\
&=
\frac{\hat{b}^{\hat{a}}}{\hat{x} ! \Gamma(\hat{a})}
\frac{\Gamma(\hat{x} + \hat{a})}{(1 + \hat{b})^{\hat{x} + \hat{a}}} \\
&=
\frac{\Gamma(\hat{x} + \hat{a})}{\hat{x} ! \Gamma(\hat{a})}
\frac{\hat{b}^{\hat{a}}}{(1 + \hat{b})^{\hat{x} + \hat{a}}}
\end{split}
\end{equation}よって、予測分布を求めることができた。
とはいえ、予測分布を解析的に求めることができただけでは、なんとなく理解も十分ではないと思うので、今回はこれを実装してみます。
ポアソン分布のベイズ推定を実装してみる
実装というのは大げさですが、Poisson分布のパラメータ $\lambda$の分布が、事前分布と事後分布でどのように変わったのかみてみましょう。
from scipy.stats import gamma
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0, 10,100)
# 事前分布のパラメータ
a = 1.0
b = 1.0
# データ
data = np.array([0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 2, 1, 1, 1, 2, 3, 1 ,2, 0, 0, 0, 2, 1, 0, 2, 1, 0])
# データから事後分布のパラメータを算出
def possison_posterior(data, a, b):
N = data.shape[0]
hat_a = np.sum(data) + a
hat_b = N + b
return hat_a, hat_b
hat_a, hat_b = possison_posterior(data, 1.0, 1.0)
plt.plot(x, gamma.pdf(x,a,scale=1./b), label="prior", color="red")
plt.plot(x, gamma.pdf(x,hat_a,scale=1./hat_b),label="posterior", color="blue")
plt.xlabel('$\lambda$')
plt.ylabel('$p(\lambda |a,b)$')
plt.legend()
plt.grid()
plt.show()コードを実行すると、事前分布と解析的に導出した事後分布を図示することができました。
