Bakery Sales Dataset : (2) 日毎の販売回数をポアソン分布に従うとして、パラメタの事後分布を推定

 今回は、日毎の販売回数がポアソン分布に従うとして事後分布を推定していきます。

(1) 共役自然分布を使ったポアソン分布の事後分布推定の概要

ポアソン分布の尤度関数は次のように与えられる：

$$ p(X|\lambda) = \prod_{i=1}^{n} \frac{e^{-\lambda} \lambda^{x_i}}{x_i!} $$

ここで、$ X=(x_1,\dots,x_n) $ は観測データである。

ガンマ事前分布は次のように与えられる：

$$ p(\lambda|a,b) = \frac{b^a}{\Gamma(a)} \lambda^{a-1} e^{-b\lambda} $$

ベイズの定理により、

$$ p(\lambda|X) = \frac{p(\lambda,X)}{p(X)} = \frac{p(X|\lambda)p(\lambda|a,b)}{p(X)} \propto p(X|\lambda)p(\lambda|a,b) $$

$$ = \prod_{i=1}^{n} \frac{e^{-\lambda} \lambda^{x_i}}{x_i!} \cdot \frac{b^a}{\Gamma(a)} \lambda^{a-1} e^{-b\lambda} $$

$$ \propto \lambda^{\left(\sum_{i=1}^{n} x_i + a - 1\right)} e^{-(n+b)\lambda} $$

以上より、事後分布はガンマ分布となり、そのパラメータは次のようになる：

$$ a_{\text{post}} = \sum_{i=1}^{n} x_i + a $$

$$ b_{\text{post}} = n + b $$

したがって、$ \lambda $ の事後分布は次のようになる：

$$ p(\lambda|X) = \frac{b_{\text{post}}^{a_{\text{post}}}}{\Gamma(a_{\text{post}})} \lambda^{a_{\text{post}}-1} e^{-b_{\text{post}}\lambda} $$

つまり、ポアソン分布のパラメタ$\lambda$の事前分布をガンマ分布とおくと、$\lambda$の事後分布は同じくガンマ分布になる。

(2) 販売回数がポアソン分布に従うときのパラメタ推定

ここでは、1日ごとの販売回数を集計し、そのデータを用いてポアソン分布のパラメータのベイズ推定を行います。その後、事後分布を求めてみましょう。まずは、日付ごとの販売回数を集計します。

ここでは、事前分布としてガンマ分布を使用し、事後分布を計算してみましょう。ガンマ分布は、形状パラメータ（$\alpha$）と尺度パラメータ（$\beta$）を持ちます。これらのパラメータを適切に選ぶことが重要です。仮に、$\alpha=\beta=1$（一様事前分布に近い）として計算してみます。

ポアソン分布のパラメータ（$\lambda$）の事後分布をプロットしました。このグラフは、観測データとガンマ分布を事前分布とした場合の$\lambda$の事後確率密度を示しています。

この事後分布から、平均販売回数の最も確からしい値や、その不確実性を評価することができます。例えば、事後分布のピークは$\lambda$の最も可能性の高い値を示し、分布の幅は不確実性の程度を示します。

(3) 得られた$\lambda$の事後分布からのMAP（Maximum A Posteriori）推定
AP推定は事後分布を最大化する$\lambda$の値を見つけることです。ガンマ分布の場合、$a > 1$ のとき、$\lambda$のMAP推定値は次のようになります： 

$$ \hat{\lambda}_{\text{MAP}} = \arg\max_\lambda (p(X|\lambda) \cdot p(\lambda|a,b)) $$ 

$$ \hat{\lambda}_{\text{MAP}} = \frac{a_{\text{post}} - 1}{b_{\text{post}}} $$ 

ガンマ分布の最大値はそのモード（最頻値）となる性質があります。

実際求めてみます。

ポアソン分布のパラメータ（λ）の最大事後確率（MAP）推定値は約59.16です。これは、観測データをもとに、平均販売回数が一日あたり約59回であると推定されることを意味します。この値は、データと事前分布から導き出された最も確からしい販売回数の推定値です。

今度は、$\alpha=2, \beta=1$でやって、先ほどの結果と比較してみます。新しい事前分布に基づいた$\lambda$の最大事後確率（MAP）推定値は約59.16で同じでした。これは、事前分布を変更しても、観測データの影響が大きいため、MAP推定値に大きな変化が見られないことを示しています。この結果は、データが豊富であればあるほど、事前分布の影響が小さくなるというベイズ推定の特性を反映しています。

(4)事後予測分布

事後予測分布（Posterior Predictive Distribution）は、ベイズ統計学の枠組みで使用される重要な概念の1つで、ベイズ推論を行った結果から未知のデータポイントの分布を予測するための確率分布です。事後予測分布は以下のように表されます: 

$$ P(X_{\text{new}}|X) = \int P(X_{\text{new}}|\lambda)P(\lambda|X)d\lambda $$ 

ここで、ポアソン分布のパラメータ $\lambda$ がガンマ分布に従う場合、事後予測分布は負の二項分布（またはポリヤ分布）になります。 

負の二項分布の確率質量関数（PMF）は以下のように表されます: 

$$ P(x_{\text{new}}|a_{\text{post}}, b_{\text{post}}, X) = \frac{\Gamma(x_{\text{new}} + a_{\text{post}})}{x_{\text{new}}!\Gamma(a_{\text{post}})} \left(\frac{b_{\text{post}}}{b_{\text{post}} + 1}\right)^{a_{\text{post}}} \left(\frac{1}{b_{\text{post}} + 1}\right)^{x_{\text{new}}} $$

これをやろうとすると、先ほど求めたパラメタ$\alpha_post=9466, \beta=160 (\alpha=2, \beta=1)$となってしまい、累乗の計算が困難です。以下のように実施します。

最後にMCMCで事後予測分布を可視化します。

ここまでのソースコード