はじめに

　『ベイズ推論による機械学習入門』(MLSシリーズ)の独学時のノートです。各種のモデルやアルゴリズムについて「数式・プログラム・図」を用いて解説します。
　本の補助として読んでください。

　この記事では、ポアソン分布に対するベイズ推論をPythonでスクラッチ実装します。

【前節の内容】

【他の節の内容】

【この節の内容】

3.2.3 ポアソン分布のベイズ推論の実装

　ポアソンモデル(Poisson model)に対するベイズ推論(Bayesian inference)を実装して、人工的に生成したデータを用いて、パラメータの学習と未知変数の予測を行う。ポアソンモデルでは、尤度関数をポアソン分布(Poisson distribution)、事前分布をガンマ分布(Gamma distribution)とする。この記事では、Pythonを利用して実装する。
　ポアソンモデルについては「3.2.3：ポアソンモデルの生成モデルの導出【緑ベイズ入門のノート】 - からっぽのしょこ」、ベイズ推論については「3.2.3：ポアソン分布のベイズ推論の導出【緑ベイズ入門のノート】 - からっぽのしょこ」、Rを利用する場合は「【R】3.2.3：ポアソン分布のベイズ推論の実装【緑ベイズ入門のノート】 - からっぽのしょこ」を参照のこと。

　利用するライブラリを読み込む。

# ライブラリを読込
import numpy as np
from scipy.stats import poisson, gamma, nbinom
import matplotlib.pyplot as plt

ベイズ推論の実装

　まずは、ポアソン分布に対するベイズ推論における一連の処理を確認する。
　観測モデル(ポアソンモデル)を設定して、モデルに従うデータ(トイデータ)を生成する。続いて、生成した観測データを用いて、事後分布の計算(パラメータの推定)を行う。さらに、事後分布のパラメータ(または観測データ)を用いて、予測分布の計算(未観測データの予測)を行う。

生成分布の設定

　データの生成分布(真の分布・ポアソン分布) $p(x_n \mid \lambda_{\mathrm{truth}}) = \mathrm{Poi}(x_n \mid \lambda_{\mathrm{truth}})$ のパラメータ(真のパラメータ) $\lambda_{\mathrm{truth}}$ を設定する。
　ポアソン分布については「【Python】ポアソン分布の作図 - からっぽのしょこ」を参照のこと。

　生成分布のパラメータ $\lambda_{\mathrm{truth}}$ を設定する。

# 真のパラメータを指定
lambda_truth = 4.0
print(lambda_truth)

4.0

　ポアソン分布の期待値パラメータ(正の値) $\lambda_{\mathrm{truth}} \gt 0$ を指定する。
　 $\lambda_{\mathrm{truth}}$ がポアソンモデルにおける真のパラメータであり、この値を求めるのがここでの目的(学習)である。

　生成分布の確率変数 $x$ の作図範囲を設定する。

# x軸の範囲を設定
x_min = 0
#x_max = 15
u = 5.0
k = 3.0
x_max = lambda_truth # 基準値を指定
x_max *= k # 定数倍
#x_max = max(x_max, x_n.max()) # サンプルと比較
x_max = np.ceil(x_max /u)*u # u単位で切り上げ
x_max = x_max.astype(dtype=np.int32) # 整数型に変換
print(x_min, x_max)

# x軸の値を作成
x_vec = np.arange(start=x_min, stop=x_max+1, step=1)
print(x_vec[:5])

0 15
[0 1 2 3 4]

　この例では、指定したパラメータ(または生成したデータ)を使って、範囲を設定している。

　生成分布の確率を計算する。

# 生成分布の確率を計算:式(2.37)
model_prob_vec = poisson.pmf(k=x_vec, mu=lambda_truth)
print(model_prob_vec[:5])

[0.01831564 0.07326256 0.14652511 0.19536681 0.19536681]

　 $x$ の値ごとに、ポアソン分布に従う確率 $\mathrm{Poi}(x \mid \lambda_{\mathrm{truth}})$ を計算する。
　ポアソン分布の確率質量関数 sp.stats.poisson.pmf() の確率変数の引数 k に $x$ 、パラメータの引数 mu に $\lambda$ を指定する。

　生成分布のグラフを作成する。

# 生成分布のラベルを作成
model_param_lbl = f'$\\lambda_{{truth}} = {lambda_truth:.3g}$'

# 生成分布を作図
fig, ax = plt.subplots(
    figsize=(9, 6), dpi=100, facecolor='white', 
    constrained_layout=True
)
fig.suptitle('Poisson distribution', fontsize=20)
ax.bar(
    x=x_vec, height=model_prob_vec, 
    color='red', 
    label='true model', zorder=10
) # 真の分布
ax.set_xticks(ticks=x_vec) # x軸目盛
ax.set_xlabel('$x$')
ax.set_ylabel('probability')
ax.set_title(model_param_lbl, loc='left') # パラメータラベル
ax.legend()
ax.grid(zorder=0)
plt.show()

　真の分布を赤色の棒グラフで示す。

　真のパラメータ $\lambda_{\mathrm{truth}}$ を求めることは、真の分布 $\mathrm{Poi}(x_n \mid \lambda_{\mathrm{truth}})$ を求めることを意味する。

データの生成

　設定した生成分布(ポアソン分布) $p(x_n \mid \lambda_{\mathrm{truth}}) = \mathrm{Poi}(x_n \mid \lambda_{\mathrm{truth}})$ に従うデータ(観測データ) $\mathbf{X}$ を作成する。
　ポアソンモデルのデータ生成については「【Python】ポアソン分布の乱数生成 - からっぽのしょこ」を参照のこと。

　生成分布からデータ $\mathbf{X}$ を生成する。

# データ数を指定
N = 100

# 観測データを生成
x_n = np.random.poisson(lam=lambda_truth, size=N)
print(x_n[:10])

[2 6 5 4 5 3 3 6 5 5]

　データ数 $N$ を指定して、ポアソン分布に従う乱数 $\mathbf{X} = \{x_1, \cdots, x_N\}$ を生成する。
　ポアソン分布の乱数生成関数 np.random.poisson() のパラメータの引数 lamに $\lambda_{\mathrm{truth}}$ 、サンプルサイズの引数 size に $N$ を指定する。

　観測データ $\mathbf{X}$ を集計する。

# 度数を集計
obs_freq_vec = np.array([np.sum(x_n == x) for x in x_vec])
print(obs_freq_vec[:5])

# 相対度数を計算
obs_relfreq_vec = obs_freq_vec / N
print(obs_relfreq_vec[:5])

[ 1  8 14 19 20]
[0.01 0.08 0.14 0.19 0.2 ]

　 $x$ の値ごとに、観測データ x_n に含まれる要素数をカウントして、度数 $N_x$ と相対度数 $\frac{N_x}{N}$ を求める。

　観測データ $\mathbf{X}$ のグラフを作成する。

# 観測データのラベルを作成
obs_param_lbl  = f'$N = {N}, '
obs_param_lbl += f'\\lambda_{{truth}} = {lambda_truth:.3g}$'

# 相対度数軸の範囲を設定
u = 0.05
relfreq_max = np.max(obs_relfreq_vec)
relfreq_max = max(relfreq_max, model_prob_vec.max()) # 生成確率と比較
relfreq_max = np.ceil(relfreq_max /u)*u # u単位で切り上げ

# 観測データを作図
fig, ax = plt.subplots(
    figsize=(9, 6), dpi=100, facecolor='white', 
    constrained_layout=True
)
fig.suptitle('Poisson distribution', fontsize=20)
ax2 = ax.twinx() # 第2軸の設定用

ax.bar(
    x=x_vec, height=model_prob_vec, 
    facecolor='none', edgecolor='red', linewidth=1.0, linestyle='--', 
    label='true model\n(generation distribution)', zorder=10
) # 生成分布
ax.bar(
    x=x_vec, height=obs_relfreq_vec, 
    color='hotpink', alpha=0.5, 
    label='observation data', zorder=11
) # 観測データ

ax.set_xticks(ticks=x_vec) # x軸目盛
relfreq_vals = ax.get_yticks()  # 相対度数
freq_vals    = relfreq_vals * N # 度数
ax2.set_yticks(ticks=freq_vals, labels=[f'{y:.1f}' for y in freq_vals]) # 度数軸目盛
ax.set_xlabel('$x$')
ax.set_ylabel('probability, relative frequency')
ax2.set_ylabel('frequency')
ax.set_title(obs_param_lbl, loc='left') # パラメータラベル
ax.legend()
ax.grid(zorder=0)
ax.set_ylim(ymin=0.0, ymax=relfreq_max) # (目盛の共通化用)
ax2.set_ylim(ymin=0.0, ymax=relfreq_max*N) # (目盛の共通化用)

plt.show()

　 $x$ の値ごとの生成確率(生成分布)を赤色(白抜き)の棒グラフ(破線)、観測データ $\mathbf{X}$ の度数 $N_x$ (相対度数 $\frac{N_x}{N}$ )を桃色(塗りつぶし)の棒グラフ(実線)で示す。

　データ数 $N$ が十分に大きいと、観測データ $\mathbf{X}$ のヒストグラムの形状が生成分布 $\mathrm{Poi}(x_n \mid \lambda_{\mathrm{truth}})$ に近付く。

事前分布の設定

　パラメータ $\lambda$ の事前分布(ガンマ分布) $p(\lambda \mid a, b) = \mathrm{Gam}(\lambda \mid a, b)$ のパラメータ(超パラメータ) $a, b$ を設定する。
　ガンマ分布については「【Python】ガンマ分布の作図 - からっぽのしょこ」を参照のこと。

　事前分布のパラメータ $a, b$ を設定する。

# 事前分布のパラメータを指定
a = 1.0
b = 1.0
print(a)
print(b)

1.0
1.0

　ガンマ分布の形状パラメータ(正の値) $a \gt 0$ 、尺度パラメータ(正の値) $b \gt 0$ を指定する。

　事前分布の確率変数 $\lambda$ の作図範囲を設定する。

# λ軸の範囲を設定
lambda_min = 0.0
#lambda_max = 15.0
u = 5.0
k = 3.0
lambda_max = lambda_truth # 基準値を指定
lambda_max *= k # 定数倍
lambda_max = np.ceil(lambda_max /u)*u # u単位で切り上げ
print(lambda_min, lambda_max)

# λ軸の値を作成
lambda_vec = np.linspace(start=lambda_min, stop=lambda_max, num=1001)
print(lambda_vec[:5])

0.0 15.0
[0.    0.015 0.03  0.045 0.06 ]

　この例では、指定したパラメータを使って、範囲を設定している。

　事前分布の確率密度を計算する。

# 事前分布の確率密度を計算:式(2.56)
prior_dens_vec = gamma.pdf(x=lambda_vec, a=a, scale=1.0/b)
print(prior_dens_vec[:5])

[1.         0.98511194 0.97044553 0.95599748 0.94176453]

　 $\lambda$ の値ごとに、ガンマ分布に従う確率密度 $\mathrm{Gam}(\lambda \mid a, b)$ を計算する。
　ガンマ分布の確率密度関数 sp.stats.gamma.pdf() の確率変数の引数 x に $\lambda$ 、形状パラメータの引数 a に $a$ 、尺度パラメータの引数 scale に $\frac{1}{b}$ を指定する。

　事前分布のグラフを作成する。

# 事前分布のラベルを作成
prior_param_lbl  = f'$\\lambda_{{truth}} = {lambda_truth:.2f}, '
prior_param_lbl += f'a = {a:.3g}, b = {b:.3g}$'

# 事前分布を作図
fig, ax = plt.subplots(
    figsize=(9, 6), dpi=100, facecolor='white', 
    constrained_layout=True
)
fig.suptitle('Gamma distribution', fontsize=20)
ax2 = ax.twiny() # 第2軸の設定用

ax.axvline(
    x=lambda_truth, 
    color='red', linewidth=1.0, linestyle='--', 
    label='true parameter', zorder=10
) # 真のパラメータ
ax.plot(
    lambda_vec, prior_dens_vec, 
    color='purple', linewidth=1.0, 
    label='posterior distribution', zorder=11
) # 事前分布

ax2.set_xticks(ticks=[lambda_truth], labels=['$\\lambda_{truth}$']) # パラメータラベル
ax.set_xlabel('$\\lambda$')
ax.set_ylabel('density')
ax.set_title(prior_param_lbl, loc='left') # パラメータラベル
ax.legend()
ax.grid(zorder=0)
ax.set_xlim(xmin=lambda_min, xmax=lambda_max) # (目盛の共通化用)
ax2.set_xlim(xmin=lambda_min, xmax=lambda_max) # (目盛の共通化用)

plt.show()

　真のパラメータを赤色の垂直線(破線)、事前分布(ガンマ分布)を紫色の曲線(実線)で示す。

　真のパラメータ(真の分布のパラメータ) $\lambda_{\mathrm{truth}}$ と、パラメータ $\lambda$ の事前分布 $\mathrm{Gam}(\lambda \mid a, b)$ の位置関係を図で確認する。

事後分布の計算

　観測データ $\mathbf{X}$ から、パラメータ $\lambda$ の事後分布(ガンマ分布) $p(\lambda \mid \mathbf{X}, a, b) = \mathrm{Gam}(\lambda \mid \hat{a}, \hat{b})$ のパラメータ(超パラメータ) $\hat{a}, \hat{b}$ を求める(真のパラメータ $\lambda_{\mathrm{truth}}$ を分布推定する)。

　事後分布のパラメータ $\hat{a}, \hat{b}$ を計算する。

# 事後分布のパラメータを計算:式(3.38)
a_hat = np.sum(x_n) + a
b_hat = N + b
print(a_hat)
print(b_hat)

401.0
101.0

　事後分布のパラメータは、観測データ $\mathbf{X}$ を用いて、次の式で計算できる。

$\displaystyle \begin{aligned} \hat{a} &= \sum_{n=1}^N x_n + a \\ \hat{b} &= N + b \end{aligned} \tag{3.38}$

　事後分布の確率密度を計算する。

# 事後分布の確率密度を計算:式(2.56)
posterior_dens_vec = gamma.pdf(x=lambda_vec, a=a_hat, scale=1.0/b_hat)
print(posterior_dens_vec[:5])

[0. 0. 0. 0. 0.]

　事前分布のときと同様にして、ガンマ分布に従う確率密度 $\mathrm{Gam}(\lambda \mid \hat{a}, \hat{b})$ を計算する。

　事後分布のグラフを作成する。

# 事後分布のラベルを作成
posterior_param_lbl  = f'$N = {N}, '
posterior_param_lbl += f'\\lambda_{{truth}} = {lambda_truth:.3g}, '
posterior_param_lbl += f'\\hat{{a}} = {a_hat:.3g}, \\hat{{b}} = {b_hat:.3g}$'

# 事後分布を作図
fig, ax = plt.subplots(
    figsize=(9, 6), dpi=100, facecolor='white', 
    constrained_layout=True
)
fig.suptitle('Gamma distribution', fontsize=20)
ax2 = ax.twiny() # 第2軸の設定用

ax.axvline(
    x=lambda_truth, 
    color='red', linewidth=1.0, linestyle='--', 
    label='true parameter', zorder=10
) # 真のパラメータ
ax.plot(
    lambda_vec, posterior_dens_vec, 
    color='purple', linewidth=1.0, 
    label='posterior distribution', zorder=11
) # 事後分布

ax2.set_xticks(ticks=[lambda_truth], labels=['$\\lambda_{truth}$']) # パラメータラベル
ax.set_xlabel('$\\lambda$')
ax.set_ylabel('density')
ax.set_title(posterior_param_lbl, loc='left') # パラメータラベル
ax.legend()
ax.grid(zorder=0)
ax.set_xlim(xmin=lambda_min, xmax=lambda_max) # (目盛の共通化用)
ax2.set_xlim(xmin=lambda_min, xmax=lambda_max) # (目盛の共通化用)

plt.show()

　真のパラメータを赤色の垂直線(破線)、事後分布(ガンマ分布)を紫色の曲線(実線)で示す。

　データ数 $N$ が十分に大きいと、パラメータ $\lambda$ の事後分布 $\mathrm{Gam}(\lambda \mid \hat{a}, \hat{b})$ のピークが真のパラメータ $\lambda_{\mathrm{truth}}$ に近付く。

予測分布の計算

　観測データ $\mathbf{X}$ から、未観測データ $x_{*}$ の予測分布(負の二項分布) $p(x_{*} \mid \mathbf{X}, a, b) = \mathrm{NB}(x_{*} \mid \hat{r}, \hat{p})$ のパラメータ $\hat{r}, \hat{p}$ を求める。
　負の二項分布については「準備中」を参照のこと。

　予測分布のパラメータ $\hat{r}, \hat{p}$ を計算する。

# 事後分布のパラメータにより予測分布のパラメータを計算:式(3.44')
r_hat = a_hat
q_hat = 1.0 / (1.0 + b_hat)
p_hat = b_hat / (1.0 + b_hat)
print(r_hat)
print(q_hat)
print(p_hat)

401.0
0.00980392156862745
0.9901960784313726

# 観測データにより予測分布のパラメータを計算:式(3.44')
r_hat = np.sum(x_n) + a
q_hat = 1.0 / (N + 1.0 + b)
p_hat = (N + b) / (N + 1.0 + b)
print(r_hat)
print(q_hat)
print(p_hat)

401.0
0.00980392156862745
0.9901960784313726

# 失敗確率パラメータから成功確率パラメータを計算
p_hat = 1.0 - q_hat
print(p_hat)

0.9901960784313726

　予測分布のパラメータは、事後分布のパラメータ $\hat{a}, \hat{b}$ または観測データ $\mathbf{X}$ を用いて、次の式で計算できる。

$\displaystyle \begin{aligned} \hat{r} &= \hat{a} \\ &= \sum_{n=1}^N x_n + a \\ \hat{q} &= \frac{1}{1 + \hat{b}} \\ &= \frac{1}{N + 1 + b} \\ \hat{p} &= 1 - \hat{q} \\ &= \frac{\hat{b}}{1 + \hat{b}} \\ &= \frac{N + b}{N + 1 + b} \end{aligned} \tag{3.44'}$

　予測分布の確率を計算する。

# 予測分布の確率密度を計算:式(3.43)
predict_prob_vec = nbinom.pmf(k=x_vec, n=r_hat, p=p_hat)
print(predict_prob_vec[:5])

[0.01923986 0.07563908 0.14905347 0.19630245 0.19437791]

　 $x$ の値ごとに、負の二項分布に従う確率 $\mathrm{NB}(x \mid \hat{r}, \hat{p})$ を計算する。
　負の二項分布の確率質量関数 sp.stats.nbinom.pmf() の確率変数の引数 k に $x$ 、成功回数の引数 n に $\hat{r}$ 、成功確率の引数 p に $\hat{p}$ を指定する。

　予測分布のグラフを作成する。

# 予測分布のラベルを作成
predict_param_lbl  = f'$N = {N}, '
predict_param_lbl += f'\\lambda_{{truth}} = {lambda_truth:.3g}, '
predict_param_lbl += f'\\hat{{r}} = {r_hat:.3g}, \\hat{{p}} = {p_hat:.3g}$'

# 予測分布を作図
fig, ax = plt.subplots(
    figsize=(9, 6), dpi=100, facecolor='white', 
    constrained_layout=True
)
fig.suptitle('Negative Binomial distribution', fontsize=20)
ax.bar(
    x=x_vec, height=model_prob_vec, 
    facecolor='none', edgecolor='red', linewidth=1.0, linestyle='--', 
    label='true model', zorder=10
) # 真の分布
ax.bar(
    x=x_vec, height=predict_prob_vec, 
    color='purple', alpha=0.5, 
    label='predict distribution', zorder=11
) # 予測分布
ax.set_xticks(ticks=x_vec) # x軸目盛
ax.set_xlabel('$x$')
ax.set_ylabel('probability')
ax.set_title(predict_param_lbl, loc='left') # パラメータラベル
ax.legend()
ax.grid(zorder=0)
plt.show()

　真の分布を赤色(白抜き)の棒グラフ(破線)、予測分布(負の二項分布)を紫色(塗りつぶし)の棒グラフ(実線)で示す。

　データ数 $N$ が十分に大きいと、未観測データ $x_{*}$ の予測分布 $\mathrm{NB}(x_{*} \mid \hat{r}, \hat{p})$ の形状が真の分布 $\mathrm{Poi}(x_{n} \mid \lambda_{\mathrm{truth}})$ に近付く。

　以上で、ポアソンモデルのベイズ推論を実装した。

スポンサードリンク

学習の推移

　次は、ポアソン分布に対するベイズ推論を図で確認する。
　データ数を増やして分布の変化をアニメーションで確認する。
　作図コードについては「Suyama-Bayes/code/poisson_model/bayesian_inference/plot_parameter_updates.py at master · anemptyarchive/Suyama-Bayes · GitHub」を参照のこと。

データ数と分布の関係

　データ数 $N$ を変化させたときの事後分布 $p(\lambda \mid \mathbf{X}, a, b)$ の変化をアニメーションにする。

　 $n$ 個のデータから求めた( $n$ 回更新した)事後分布(ガンマ分布) $\mathrm{Gam}(\lambda \mid a^{(n)}, b^{(n)})$ を紫色の曲線(実線)、 $n$ 番目の観測データ $x_n$ に対応する位置 $\lambda = x_n$ を桃色の点で示す。

　「ベイズ推論の実装」では、 $N$ 個(複数)のデータ $\mathbf{X}$ を用いて、事後分布のパラメータ $\hat{a}, \hat{b}$ を一括更新した。
　 $n$ 番目(1つ)のデータ $x_n$ を用いて逐次更新する場合、 $n$ 回目の事後分布のパラメータ $a^{(n)}, b^{(n)}$ は、次の式で計算できる。

$\displaystyle \begin{aligned} a^{(n)} &\leftarrow x_n + a^{(n-1)} \\ b^{(n)} &\leftarrow 1 + b^{(n-1)} \end{aligned}$

　超パラメータの初期値(1回目の更新における事前分布のパラメータ)を $a^{(0)}, b^{(0)}$ として、 $a^{(n-1)}, b^{(n-1)}$ は、 $n-1$ 回目の更新における事後分布のパラメータであり、また $n$ 回目の更新における事前分布のパラメータを表す。

　データ数 $N$ が大きくなるのに応じて、パラメータ $\lambda$ の事後分布 $\mathrm{Gam}(\lambda \mid \hat{a}, \hat{b})$ のピークが真のパラメータ $\lambda_{\mathrm{truth}}$ に近付いていくのを確認できる。

　データ数 $N$ を変化させたときの予測分布 $p(x_{*} \mid \mathbf{X}, a, b)$ の変化をアニメーションにする。

　 $n$ 個のデータから求めた( $n$ 回更新した)予測分布(負の二項分布) $\mathrm{NB}(x_{*} \mid r^{(n)}, p^{(n)})$ を紫色(塗りつぶし)の棒グラフ(実線)、 $n$ 番目の観測データ $x_n$ を桃色の点で示す。

　「ベイズ推論の実装」では、 $N$ 個(複数)のデータ $\mathbf{X}$ を用いて、予測分布のパラメータ $\hat{r}, \hat{p}$ を一括更新した。
　 $n$ 番目(1つ)のデータ $x_n$ を用いて逐次更新する場合、 $n$ 回目の予測分布のパラメータ $r^{(n)}, p^{(n)}$ は、次の式で計算できる。

$\displaystyle \begin{aligned} r^{(n)} &\leftarrow x_n + r^{(n-1)} \\ q^{(n)} &\leftarrow \frac{1}{1 + \frac{1}{q^{(n-1)}}} \\ p^{(n)} &= 1 - q^{(n)} \end{aligned}$

　超パラメータの初期値 $a^{(0)}, b^{(0)}$ を用いて求めたパラメータを $r^{(0)}, p^{(0)}$ として、 $r^{(n-1)}, p^{(n-1)}$ は $n-1$ 回目の更新値を表す。

　データ数 $N$ が大きくなるのに応じて、未観測データ $x_{*}$ の予測分布 $\mathrm{NB}(x_{*} \mid \hat{r}, \hat{p})$ の形状が真の分布 $\mathrm{Poi}(x_{n} \mid \lambda_{\mathrm{truth}})$ に近付いていくのを確認できる。

観測データと分布の関係

　データ数 $N$ の変化による観測データと事後分布、予測分布の関係をアニメーションにする。

　真の分布の期待値 $\mathbb{E}[x_n]$ と真のパラメータ $\lambda_{\mathrm{truth}}$ の各図(分布)に対応する位置を赤色の垂直線(破線)、観測データの標本平均 $\bar{x}$ を桃色の垂直線(破線)、事後分布の期待値 $\mathbb{E}[\lambda]$ と予測分布の期待値 $\mathbb{E}[x_{*}]$ を紫色の垂直線(破線)で示す。

　生成分布(ポアソン分布)、事後分布(ガンマ分布)、予測分布(負の二項分布)の期待値は、それぞれ次の式で計算できる。

$\displaystyle \begin{align} \mathbb{E}_{\mathrm{Poi}(x_n \mid \lambda_{\mathrm{truth}})}[x] &= \lambda_{\mathrm{truth}} \\ \mathbb{E}_{\mathrm{Gam}(\lambda \mid \hat{a}, \hat{b})}[\lambda] &= \frac{\hat{a}}{\hat{b}} \\ \mathbb{E}_{\mathrm{NB}(x_{*} \mid \hat{r}, \hat{p})}[x] &= \frac{\hat{r} (1 - \hat{p})}{\hat{p}} \\ &= \frac{\hat{a}}{\hat{b}} \end{align}$

　真の分布 $\mathrm{Poi}(x_n \mid \lambda_{\mathrm{truth}})$ や真のパラメータ $\lambda_{\mathrm{truth}}$ と、観測データ $\mathbf{X}$ 、事後分布 $\mathrm{Gam}(\lambda \mid \hat{a}, \hat{b})$ 、予測分布 $\mathrm{NB}(x_{*} \mid \hat{r}, \hat{p})$ 、またそれぞれの統計量の対応関係を確認できる。

　以上で、ポアソンモデルのベイズ推論における学習推移を確認した。

　この記事では、ポアソン分布に対するベイズ推論を実装した。次の記事では、1次元ガウス分布を扱う。