はじめに

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

　この記事では、平均が未知の1次元ガウス分布に対するベイズ推論をPythonでスクラッチ実装します。

【前節の内容】

【他の節の内容】

【この節の内容】

3.3.1 1次元ガウス分布のベイズ推論の実装：平均が未知の場合

　1次元ガウスモデル(Gaussian model)に対するベイズ推論(Bayesian inference)を実装して、人工的に生成したデータを用いて、パラメータの学習と未知変数の予測を行う。この記事では、生成分布の平均パラメータ(mean parameter)が未知の場合を扱う。平均が未知の1次元ガウスモデルでは、尤度関数を1次元ガウス分布(Gaussian distribution・一変量正規分布・Normal distribution)、事前分布を1次元ガウス分布とする。この記事では、Pythonを利用して実装する。
　1次元ガウスモデルについては「3.3.0：1次元ガウスモデルの生成モデルの導出【緑ベイズ入門のノート】 - からっぽのしょこ」、ベイズ推論については「3.3.1：1次元ガウス分布のベイズ推論の導出：平均が未知の場合【緑ベイズ入門のノート】 - からっぽのしょこ」、Rを利用する場合は「【R】3.3.1：1次元ガウス分布のベイズ推論の実装：平均が未知の場合【緑ベイズ入門のノート】 - からっぽのしょこ」を参照のこと。

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

# ライブラリを読込
import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt

ベイズ推論の実装

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

生成分布の設定

　データの生成分布(真の分布・ガウス分布) $p(x_n \mid \mu_{\mathrm{truth}}, \lambda) = \mathcal{N}(x_n \mid \mu_{\mathrm{truth}}, \lambda^{-1})$ のパラメータ(真のパラメータ) $\mu_{\mathrm{truth}}, \lambda$ を設定する。 $\lambda$ は既知とする。
　ガウス分布については「【Python】1次元ガウス分布の作図 - からっぽのしょこ」を参照のこと。

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

# 真のパラメータを指定
mu_truth = 25.0

# 既知のパラメータを指定
lmd = 0.01
print(mu_truth)
print(lmd)
print(1.0/np.sqrt(lmd))

25.0
0.01
10.0

　ガウス分布の平均パラメータ(実数) $\mu_{\mathrm{truth}}$ 、精度パラメータ(正の値) $\lambda \gt 0$ を指定する。
　 $\mu_{\mathrm{truth}}$ がガウスモデルにおける真のパラメータであり、この値を求めるのがここでの目的(学習)である。

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

# x軸の範囲を設定
#x_min = -15.0
#x_max = 65.0
k = 4.0
u = 5.0
x_size  = 1.0/np.sqrt(lmd) # 基準値を指定
x_size *= k # 定数倍
#x_size  = max(x_size, *np.abs(x_n-mu_truth)) # サンプルと比較
x_size  = np.ceil(x_size /u)*u # u単位で切り上げ
x_min   = mu_truth - x_size
x_max   = mu_truth + x_size
print(x_min, x_max)

# x軸の値を作成
x_vec = np.linspace(start=x_min, stop=x_max, num=1001)
print(x_vec[:5])

-15.0 65.0
[-15.   -14.92 -14.84 -14.76 -14.68]

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

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

# 生成分布の確率密度を計算:式(2.64)
model_dens_vec = norm.pdf(x=x_vec, loc=mu_truth, scale=1.0/np.sqrt(lmd))
print(model_dens_vec[:5])

[1.33830226e-05 1.38177629e-05 1.42657125e-05 1.47272413e-05
 1.52027287e-05]

　 $x$ の値ごとに、ガウス分布に従う確率密度 $\mathcal{N}(x \mid \mu_{\mathrm{truth}}, \lambda^{-1})$ を計算する。
　ガウス分布の確率密度関数 sp.stats.norm.pdf() の確率変数の引数 x に $x$ 、平均の引数 loc に $\mu_{\mathrm{truth}}$ 、標準偏差の引数 scale に $\sigma = \frac{1}{\sqrt{\lambda}}$ を指定する。

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

# 生成分布のラベルを作成
model_param_lbl = f'$\\mu_{{truth}} = {mu_truth:.3g}, \\lambda = {lmd:.3g}$'

# 生成分布を作図
fig, ax = plt.subplots(
    figsize=(9, 6), dpi=100, facecolor='white', 
    constrained_layout=True
)
fig.suptitle('Gaussian distribution', fontsize=20)
ax.plot(
    x_vec, model_dens_vec, 
    color='red', linewidth=1.0, 
    label='true model', zorder=10
) # 真の分布
ax.set_xlabel('$x$')
ax.set_ylabel('density')
ax.set_title(model_param_lbl, loc='left') # パラメータラベル
ax.legend()
ax.grid(zorder=0)
plt.show()

　真の分布(ガウス分布)を赤色の曲線で示す。

　真のパラメータ $\mu_{\mathrm{truth}}$ を求めることは、真の分布 $\mathcal{N}(x_n \mid \mu_{\mathrm{truth}}, \lambda^{-1})$ を求めることを意味する。

データの生成

　設定した生成分布(ガウス分布) $p(x_n \mid \mu_{\mathrm{truth}}, \lambda) = \mathcal{N}(x_n \mid \mu_{\mathrm{truth}}, \lambda^{-1})$ に従うデータ(観測データ) $\mathbf{X}$ を作成する。
　ガウスモデルのデータ生成については「【Python】1次元ガウス分布の乱数生成 - からっぽのしょこ」を参照のこと。

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

# データ数を指定
N = 100

# 観測データを生成
x_n = np.random.normal(loc=mu_truth, scale=1.0/np.sqrt(lmd), size=N)
print(x_n[:10].round(1))

[19.  18.9 29.6 18.8 31.9 13.5  9.2 43.  24.6 36.6]

　データ数 $N$ を指定して、ガウス分布に従う乱数 $\mathbf{X} = \{x_1, \cdots, x_N\}$ を生成する。
　ガウス分布の乱数生成関数 np.random.normal() の平均の引数 loc に $\mu_{\mathrm{truth}}$ 、標準偏差の引数 scale に $\sigma = \frac{1}{\sqrt{\lambda}}$ 、サンプルサイズの引数 size に $N$ を指定する。

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

# 階級数を指定
bin_num = 40

# 階級幅を計算
bin_size = (x_max - x_min) / bin_num

# 境界値の範囲を設定
bin_min = x_min - 0.5*bin_size
bin_max = x_max + 0.5*bin_size

# 度数を集計
obs_freq_vec, bin_vec = np.histogram(a=x_n, bins=bin_num+1, range=(bin_min, bin_max))
print(obs_freq_vec[:5])

# 密度を計算
obs_dens_vec = obs_freq_vec / (bin_size * N)
print(obs_dens_vec[:5])

# 階級値を作成
center_vec = bin_vec[:-1] + 0.5*bin_size
print(center_vec[:5])

[0 0 0 0 0]
[0. 0. 0. 0. 0.]
[-15. -13. -11.  -9.  -7.]

　階級数を指定して、階級幅 $w$ を設定する。
　 $x$ の範囲で階級値を作成して、観測データ x_n に含まれる要素数をカウントして、度数 $N_x$ と密度 $\frac{N_x}{w N}$ を求める。

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

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

# 密度軸の範囲を設定
u = 0.025
dens_max = np.max(obs_dens_vec)
dens_max = max(dens_max, model_dens_vec.max()) # 生成確率と比較
dens_max = np.ceil(dens_max /u)*u # u単位で切り上げ

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

ax.plot(
    x_vec, model_dens_vec, 
    color='red', linewidth=1.0, linestyle='--', 
    label='true model\n(generation distribution)', zorder=10
) # 生成分布
ax.bar(
    x=center_vec, height=obs_dens_vec, 
    width=bin_size, align='center', 
    color='hotpink', alpha=0.5, 
    label='observation data', zorder=11
) # 観測データ

dens_vals = ax.get_yticks()        # 密度
freq_vals = dens_vals * bin_size*N # 度数
ax2.set_yticks(ticks=freq_vals, labels=[f'{y:.1f}' for y in freq_vals]) # 度数軸目盛
ax.set_xlabel('$x$')
ax.set_ylabel('density')
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=dens_max) # (目盛の共通化用)
ax2.set_ylim(ymin=0.0, ymax=dens_max*bin_size*N) # (目盛の共通化用)

plt.show()

　 $x$ の値ごとの確率密度(生成分布)を赤色の曲線(破線)、観測データ $\mathbf{X}$ の度数 $N_x$ (密度 $\frac{N_x}{w N}$ )を桃色の棒グラフ(実線)で示す。

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

事前分布の設定

　パラメータ $\mu$ の事前分布(ガウス分布) $p(\mu \mid m, \lambda_{\mu}) = \mathcal{N}(\mu \mid m, \lambda_{\mu}^{-1})$ のパラメータ(超パラメータ) $m, \lambda_{\mu}$ を設定する。

　事前分布のパラメータ $m, \lambda_{\mu}$ を設定する。

# 事前分布のパラメータを指定
m         = 0.0
lambda_mu = 0.001
print(m)
print(lambda_mu)

0.0
0.001

　ガウス分布の平均パラメータ(実数) $m$ 、精度パラメータ(正の値) $\lambda_{\mu} \gt 0$ を指定する。

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

# μ軸の範囲を設定
mu_min = x_min
mu_max = x_max

# μ軸の値を作成
mu_vec = np.linspace(start=mu_min, stop=mu_max, num=1001)
print(mu_vec[:5])

[-15.   -14.92 -14.84 -14.76 -14.68]

　この例では、生成分布の確率変数 $x$ の範囲を使っている。

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

# 事前分布の確率密度を計算:式(2.64)
prior_dens_vec = norm.pdf(x=mu_vec, loc=m, scale=1.0/np.sqrt(lambda_mu))
print(prior_dens_vec[:5])

[0.01127332 0.01128682 0.01130027 0.01131365 0.01132698]

　 $\mu$ の値ごとに、ガウス分布に従う確率密度 $\mathcal{N}(\mu \mid m, \lambda_{\mu}^{-1})$ を計算する。
　ガウス分布の確率密度関数 sp.stats.norm.pdf() の確率変数の引数 x に $\mu$ 、平均の引数 loc に $m$ 、標準偏差の引数 scale に $\sigma_{\mu} = \frac{1}{\sqrt{\lambda_{\mu}}}$ を指定する。

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

# 事前分布のラベルを作成
prior_param_lbl  = f'$\\mu_{{truth}} = {mu_truth:.3g}, '
prior_param_lbl += f'm = {m:.3g}, \\lambda_{{\\mu}} = {lambda_mu:.3g}$'

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

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

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

plt.show()

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

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

事後分布の計算

　観測データ $\mathbf{X}$ から、パラメータ $\mu$ の事後分布(ガウス分布) $p(\mu \mid \mathbf{X}, m, \lambda_{\mu}) = \mathcal{N}(\mu \mid \hat{m}, \hat{\lambda}_{\mu}^{-1})$ のパラメータ(超パラメータ) $\hat{m}, \hat{\lambda}_{\mu}$ を求める(真のパラメータ $\mu_{\mathrm{truth}}$ を分布推定する)。

　事後分布のパラメータ $\hat{m}, \hat{\lambda}_{\mu}$ を計算する。

# 事後分布のパラメータを計算:式(3.53, 3.54)
lambda_mu_hat = N * lmd + lambda_mu
m_hat         = np.sum(x_n) * lmd + m * lambda_mu
m_hat        /= lambda_mu_hat
print(m_hat)
print(lambda_mu_hat)

25.883955369891094
1.001

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

$\displaystyle \begin{align} \hat{m} &= \frac{ \lambda \sum_{n=1}^N x_n + m \lambda_{\mu} }{ \hat{\lambda}_{\mu} } \tag{3.54}\\ \hat{\lambda}_{\mu} &= N \lambda + \lambda_{\mu} \tag{3.53} \end{align}$

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

# 事後分布の確率密度を計算:式(2.64)
posterior_dens_vec = norm.pdf(x=mu_vec, loc=m_hat, scale=1.0/np.sqrt(lambda_mu_hat))
print(posterior_dens_vec[:5])

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

　事前分布のときと同様にして、ガウス分布に従う確率密度 $\mathcal{N}(\mu \mid \hat{m}, \hat{\lambda}_{\mu}^{-1})$ を計算する。

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

# 事後分布のラベルを作成
posterior_param_lbl  = f'$N = {N}, '
posterior_param_lbl += f'\\mu_{{truth}} = {mu_truth:.3g}, '
posterior_param_lbl += f'\\hat{{m}} = {m_hat:.3g}, \\hat{{\\lambda}}_{{\\mu}} = {lambda_mu_hat:.3g}$'

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

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

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

plt.show()

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

　データ数 $N$ が十分に大きいと、パラメータ $\mu$ の事後分布 $\mathcal{N}(\mu \mid \hat{m}, \hat{\lambda}_{\mu}^{-1})$ のピークが真のパラメータ $\mu_{\mathrm{truth}}$ に近付く。

予測分布の計算

　観測データ $\mathbf{X}$ から、未観測データ $x_{*}$ の予測分布(ガウス分布) $p(x_{*} \mid \mathbf{X}, m, \lambda_{\mu}) = \mathcal{N}(x_{*} \mid \hat{\mu}_{*}, \hat{\lambda}_{*}^{-1})$ のパラメータ $\hat{\mu}_{*}, \hat{\lambda}_{*}$ を求める。

　予測分布のパラメータ $\hat{\mu}_{*}, \hat{\lambda}_{*}$ を計算する。

# 事後分布のパラメータにより予測分布のパラメータを計算:式(3.62')
mu_star_hat     = m_hat
lambda_star_hat = lmd * lambda_mu_hat / (lmd + lambda_mu_hat)
print(mu_star_hat)
print(lambda_star_hat)

25.883955369891094
0.009901088031651831

# 観測データにより予測分布のパラメータを計算:式(3.62')
mu_star_hat      = lmd * np.sum(x_n) + m * lambda_mu
mu_star_hat     /= N * lmd + lambda_mu
lambda_star_hat  = (N * lmd + lambda_mu) * lmd
lambda_star_hat /= (N + 1) * lmd + lambda_mu
print(mu_star_hat)
print(lambda_star_hat)

25.883955369891094
0.009901088031651831

# 分散パラメータから精度パラメータを計算:式(3.63')
sigma2_star_hat = 1.0/lmd + 1.0/lambda_mu_hat
lambda_star_hat = 1.0/sigma2_star_hat
print(sigma2_star_hat)
print(lambda_star_hat)

100.999000999001
0.00990108803165183

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

$\displaystyle \begin{align} \hat{\mu}_{*} &= \hat{m} \\ &= \frac{ \lambda \sum_{n=1}^N x_n + m \lambda_{\mu} }{ N \lambda + \lambda_{\mu} } \tag{3.62'}\\ \hat{\lambda}_{*} &= \frac{ \lambda \hat{\lambda}_{\mu} }{ \lambda + \hat{\lambda}_{\mu} } \\ &= \frac{ (N \lambda + \lambda_{\mu}) \lambda }{ (N + 1) \lambda + \lambda_{\mu} } \tag{3.62'}\\ \hat{\lambda}_{*}^{-1} &= \lambda^{-1} + \hat{\lambda}_{\mu}^{-1} \tag{3.63'} \end{align}$

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

# 予測分布の確率を計算:式(2.64)
predict_dens_vec = norm.pdf(x=x_vec, loc=mu_star_hat, scale=1.0/np.sqrt(lambda_star_hat))
print(predict_dens_vec[:5])

[1.01167579e-05 1.04494068e-05 1.07923097e-05 1.11457588e-05
 1.15100540e-05]

　生成分布のときと同様にして、ガウス分布に従う確率密度 $\mathcal{N}(x \mid \hat{\mu}_{*}, \hat{\lambda}_{*}^{-1})$ を計算する。

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

# 予測分布のラベルを作成
predict_param_lbl  = f'$N = {N}, '
predict_param_lbl += f'\\mu_{{truth}} = {mu_truth:.3g}, \\lambda = {lmd:.3g}, '
predict_param_lbl += f'\\hat{{\\mu}}_{{*}} = {mu_star_hat:.3g}, \\hat{{\\lambda}}_{{*}} = {lambda_star_hat:.3g}$'

# 予測分布を作図
fig, ax = plt.subplots(
    figsize=(9, 6), dpi=100, facecolor='white', 
    constrained_layout=True
)
fig.suptitle('Gaussian distribution', fontsize=20)
ax.plot(
    x_vec, model_dens_vec, 
    color='red', linewidth=1.0, linestyle='--', 
    label='true model', zorder=10
) # 真の分布
ax.plot(
    x_vec, predict_dens_vec, 
    color='purple', linewidth=1.0, 
    label='predict distribution', zorder=11
) # 予測分布
ax.set_xlabel('$x$')
ax.set_ylabel('density')
ax.set_title(predict_param_lbl, loc='left') # パラメータラベル
ax.legend()
ax.grid(zorder=0)
plt.show()

　真の分布(ガウス分布)を赤色の曲線(破線)、予測分布(ガウス分布)を紫色の曲線(実線)で示す。

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

　以上で、平均が未知の1次元ガウスモデルのベイズ推論を実装した。

スポンサードリンク

学習の推移

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

データ数と分布の関係

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

　 $n$ 個のデータから求めた( $n$ 回更新した)事後分布(ガウス分布) $\mathcal{N}(\mu \mid m^{(n)}, \lambda_{\mu}^{(n)})$ を紫色の曲線(実線)、 $n$ 番目の観測データ $x_n$ に対応する位置 $\mu = x_n$ を桃色の点で示す。

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

$\displaystyle \begin{aligned} m^{(n)} &\leftarrow \frac{ x_n \lambda + m^{(n-1)} \lambda_{\mu}^{(n-1)} }{ \lambda_{\mu}^{(n)} } \\ \lambda_{\mu}^{(n)} &\leftarrow \lambda + \lambda_{\mu}^{(n-1)} \end{aligned}$

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

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

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

　 $n$ 個のデータから求めた( $n$ 回更新した)予測分布(ガウス分布) $\mathcal{N}(x_{*} \mid \hat{\mu}_{*}, \hat{\lambda}_{*}^{-1})$ を紫色の曲線(実線)、 $n$ 番目の観測データ $x_n$ を桃色の点で示す。

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

$\displaystyle \begin{aligned} \mu_{*}^{(n)} &= m^{(n)} \\ \lambda_{*}^{(n)} &= \frac{ \lambda \lambda_{\mu}^{(n)} }{ \lambda + \lambda_{\mu}^{(n)} } \end{aligned}$

　超パラメータの初期値 $m^{(0)}, \lambda^{(0)}$ を用いて求めたパラメータを $\mu_{*}^{(0)}, \lambda_{*}^{(0)}$ として、 $\mu_{*}^{(n-1)}, \lambda_{*}^{(n-1)}$ は $n-1$ 回目の更新値を表す。

　データ数 $N$ が大きくなるのに応じて、未観測データ $x_{*}$ の予測分布 $\mathcal{N}(x_{*} \mid \hat{\mu}_{*}, \hat{\lambda}_{*}^{-1})$ の形状が真の分布 $\mathcal{N}(x_n \mid \mu_{\mathrm{truth}}, \lambda^{-1})$ に近付いていくのを確認できる。

観測データと分布の関係

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

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

　生成分布(ガウス分布)、事後分布(ガウス分布)、予測分布(ガウス分布)の期待値は、それぞれ次の式で計算できる。

$\displaystyle \begin{aligned} \mathbb{E}_{\mathcal{N}(x_n \mid \mu_{\mathrm{truth}}, \lambda^{-1})}[x] &= \mu_{\mathrm{truth}} \\ \mathbb{E}_{\mathcal{N}(\mu \mid \hat{m}, \hat{\lambda}_{\mu}^{-1})}[\mu] &= \hat{m} \\ \mathbb{E}_{\mathcal{N}(x_{*} \mid \hat{\mu}_{*}, \hat{\lambda}_{*}^{-1})}[x] &= \hat{\mu}_{*} \\ &= \hat{m} \end{aligned}$

　真の分布 $\mathcal{N}(x_n \mid \mu_{\mathrm{truth}}, \lambda^{-1})$ や真のパラメータ $\mu_{\mathrm{truth}}$ と、観測データ $\mathbf{X}$ 、事後分布 $\mathcal{N}(\mu \mid \hat{m}, \hat{\lambda}_{\mu}^{-1})$ 、予測分布 $\mathcal{N}(x_{*} \mid \hat{\mu}_{*}, \hat{\lambda}_{*}^{-1})$ 、またそれぞれの統計量の対応関係を確認できる。

　以上で、平均が未知の1次元ガウスモデルのベイズ推論における学習推移を確認した。

　この記事では、平均が未知の場合の1次元ガウス分布に対するベイズ推論を扱った。次の記事では、精度が未知の場合を扱う。