はじめに

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

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

【前節の内容】

www.anarchive-beta.com

【他の節の内容】

www.anarchive-beta.com

【この節の内容】

はじめに
3.3.3 1次元ガウス分布のベイズ推論の実装：平均と精度が未知の場合
- ベイズ推論の実装
- 学習の推移
  - データ数と分布の関係
  - 観測データと分布の関係
参考文献
おわりに

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

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

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

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

ベイズ推論の実装

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

生成分布の設定

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

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

# 真のパラメータを指定
mu_truth     = 5.0
lambda_truth = 0.25
print(mu_truth)
print(lambda_truth)
print(1.0/np.sqrt(lambda_truth))

5.0
0.25
2.0

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

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

# x軸の範囲を設定
#x_min = -5.0
#x_max = 15.0
k = 4.0
u = 5.0
x_size  = 1.0/np.sqrt(lambda_truth) # 基準値を指定
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])

-5.0 15.0
[-5.   -4.98 -4.96 -4.94 -4.92]

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

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

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

[7.43359757e-07 7.81433554e-07 8.21375294e-07 8.63272261e-07
 9.07215595e-07]

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

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

# 生成分布のラベルを作成
model_param_lbl = f'$\\mu_{{truth}} = {mu_truth:.3g}, \\lambda_{{truth}} = {lambda_truth:.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}}, \lambda_{\mathrm{truth}}$ を求めることは、真の分布 $\mathcal{N}(x_n \mid \mu_{\mathrm{truth}}, \lambda_{\mathrm{truth}}^{-1})$ を求めることを意味する。

データの生成

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

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

# データ数を指定
N = 100

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

[3.8 3.8 5.9 3.8 6.4 2.7 1.8 8.6 4.9 7.3]

　データ数 $N$ を指定して、ガウス分布に従う乱数 $\mathbf{X} = \{x_1, \cdots, x_N\}$ を生成する。
　ガウス分布の乱数生成関数 np.random.normal() の平均の引数 loc に $\mu_{\mathrm{truth}}$ 、標準偏差の引数 scale に $\sigma_{\mathrm{truth}} = \frac{1}{\sqrt{\lambda_{\mathrm{truth}}}}$ 、サンプルサイズの引数 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.]
[-5.  -4.5 -4.  -3.5 -3. ]

　階級数を指定して、階級幅 $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_{{truth}} = {lambda_truth:.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_{\mathrm{truth}}^{-1})$ に近付く。

事前分布の設定

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

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

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

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

0.0
1.0
1.0
1.0

　ガウス分布の平均パラメータ(実数) $m$ 、精度パラメータの係数(正の値) $\beta \gt 0$ と、ガンマ分布の形状パラメータ(正の値) $a \gt 0$ 、尺度パラメータ(正の値) $b \gt 0$ を指定する。

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

# μ軸の範囲を設定
mu_min = x_min
mu_max = x_max
print(mu_min, mu_max)

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

-5.0 15.0
[-5.  -4.8 -4.6 -4.4 -4.2]

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

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

# λ軸の範囲を設定
lambda_min = 0.0
#lambda_max = 1.0
k = 3.0
u = 0.5
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=101)
print(lambda_vec[:5])

0.0 1.0
[0.   0.01 0.02 0.03 0.04]

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

　格子点 $(\mu, \lambda)$ を作成する。

# 格子点を作成
mu_mat, lambda_mat = np.meshgrid(mu_vec, lambda_vec)
print(mu_mat[:5, :5])
print(mu_mat.shape)
print(lambda_mat[:5, :5])
print(lambda_mat.shape)

[[-5.  -4.8 -4.6 -4.4 -4.2]
 [-5.  -4.8 -4.6 -4.4 -4.2]
 [-5.  -4.8 -4.6 -4.4 -4.2]
 [-5.  -4.8 -4.6 -4.4 -4.2]
 [-5.  -4.8 -4.6 -4.4 -4.2]]
(101, 101)
[[0.   0.   0.   0.   0.  ]
 [0.01 0.01 0.01 0.01 0.01]
 [0.02 0.02 0.02 0.02 0.02]
 [0.03 0.03 0.03 0.03 0.03]
 [0.04 0.04 0.04 0.04 0.04]]
(101, 101)

　 $\mu, \lambda$ の格子状の点( mu_vec, lambda_vec の全ての組み合わせ)を np.meshgrid() で作成する。

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

# μの事前分布の確率密度を計算:式(2.64)
N_dens_mat = norm.pdf(x=mu_mat, loc=m, scale=1.0/np.sqrt(beta*lambda_mat))
print(N_dens_mat[:5, :5])

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

# μ,λの結合事前分布の確率密度を計算:式(3.81)
prior_dens_mat = N_dens_mat * Gam_dens_mat
print(prior_dens_mat[:5, :5])

[[0.         0.         0.         0.         0.        ]
 [0.03520653 0.03555325 0.03588903 0.03621349 0.03652627]
 [0.04393913 0.04480883 0.04565921 0.04648851 0.04729503]
 [0.04749088 0.04890784 0.05030668 0.05168346 0.05303424]
 [0.04839414 0.05032887 0.05225726 0.05417279 0.05606876]]
[[1.         1.         1.         1.         1.        ]
 [0.99004983 0.99004983 0.99004983 0.99004983 0.99004983]
 [0.98019867 0.98019867 0.98019867 0.98019867 0.98019867]
 [0.97044553 0.97044553 0.97044553 0.97044553 0.97044553]
 [0.96078944 0.96078944 0.96078944 0.96078944 0.96078944]]
[[0.         0.         0.         0.         0.        ]
 [0.03485622 0.03519949 0.03553193 0.03585316 0.03616282]
 [0.04306908 0.04392156 0.04475509 0.04556798 0.04635853]
 [0.04608732 0.0474624  0.04881989 0.05015598 0.05146684]
 [0.04649658 0.04835545 0.05020822 0.05204865 0.05387027]]

　 $\mu, \lambda$ の値の組み合わせごとに、ガウス-ガンマ分布に従う確率密度 $\mathrm{NG}(\mu, \lambda \mid m, \beta, a, b) = \mathcal{N}(\mu \mid m, (\beta \lambda)^{-1}) \mathrm{Gam}(\lambda \mid a, b)$ を計算する。
　ガウス分布の確率密度関数 sp.stats.norm.pdf() の確率変数の引数 x に $\mu$ 、平均の引数 loc に $m$ 、標準偏差の引数 scale に $\sigma_{\mu} = \frac{1}{\sqrt{\beta \lambda}}$ を指定する。
　ガンマ分布の確率密度関数 sp.stats.gamma.pdf() の確率変数の引数 x に $\lambda$ 、形状の引数 a に $a$ 、尺度の引数 scale に $\frac{1}{b}$ を指定する。
　 $\mu$ の値ごとのガウス分布の確率密度 $\mathcal{N}(\mu \mid m, (\beta \lambda)^{-1})$ と、 $\lambda$ の値ごとのガンマ分布の確率密度 $\mathrm{Gam}(\lambda \mid a, b)$ の積を計算する。

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

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

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

ax.plot(
    0.0, 0.0, 
    color=cm.viridis(X=0.0), linewidth=1.0, linestyle='-', 
    label='prior distribution', zorder=10
) # (凡例表示用のダミー)
cs = ax.contourf(
    mu_mat, lambda_mat, prior_dens_mat, 
    cmap='viridis', alpha=0.5, 
    zorder=11
) # 事前分布
ax.axvline(
    x=mu_truth, 
    color='red', linewidth=1.0, linestyle='--', 
    label='true parameter', zorder=12
) # 真のパラメータ
ax.axhline(
    y=lambda_truth, 
    color='red', linewidth=1.0, linestyle='--', 
    zorder=12
) # 真のパラメータ

ax2x.set_xticks(ticks=[mu_truth], labels=['$\\mu_{truth}$']) # パラメータラベル
ax2y.set_yticks(ticks=[lambda_truth], labels=['$\\lambda_{truth}$']) # パラメータラベル
ax.set_xlabel('$\\mu$')
ax.set_ylabel('$\\lambda$')
ax.set_title(prior_param_lbl, loc='left') # パラメータラベル
fig.colorbar(cs, ax=ax, shrink=1.0, label='density')
ax.legend(prop={'size': 8})
ax.grid(zorder=0)
ax.set_xlim(xmin=mu_min, xmax=mu_max)   # (目盛の共通化用)
ax2x.set_xlim(xmin=mu_min, xmax=mu_max) # (目盛の共通化用)
ax.set_ylim(ymin=lambda_min, ymax=lambda_max)   # (目盛の共通化用)
ax2y.set_ylim(ymin=lambda_min, ymax=lambda_max) # (目盛の共通化用)

plt.show()

# 事前分布を作図
fig, ax = plt.subplots(
    figsize=(9, 6), dpi=100, facecolor='white', constrained_layout=True, 
    subplot_kw={'projection': '3d'}
)
fig.suptitle('Gaussian-Gamma distribution', fontsize=20)

ax.plot(
    [mu_truth, mu_truth], [lambda_min, lambda_max], [0.0, 0.0], 
    color='red', linewidth=1.0, linestyle='--', 
    label='true parameter', zorder=10
) # 真のパラメータ
ax.plot(
    [mu_min, mu_max], [lambda_truth, lambda_truth], [0.0, 0.0], 
    color='red', linewidth=1.0, linestyle='--', 
    zorder=10
) # 真のパラメータ
ax.contour(
    X=mu_mat, Y=lambda_mat, Z=prior_dens_mat, offset=0.0, 
    cmap='viridis', linewidths=1.0, linestyles=':', 
    zorder=11
) # 事前分布
ax.plot_surface(
    X=mu_mat, Y=lambda_mat, Z=prior_dens_mat, 
    cmap='viridis', alpha=0.6, 
    label='prior distribution', zorder=12
) # 事前分布

ax.set_xlabel('$\\mu$')
ax.set_ylabel('$\\lambda$')
ax.set_zlabel('density')
ax.set_title(prior_param_lbl, loc='left') # パラメータラベル
ax.legend(prop={'size': 8})

plt.show()

　真のパラメータを赤色の直線(破線)、事前分布(ガウス-ガンマ分布)を等高線や曲面(塗りつぶし)で示す。

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

事後分布の計算

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

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

# μの事後分布のパラメータを計算:式(3.83)
beta_hat = N + beta
m_hat    = (np.sum(x_n) + beta * m) / beta_hat
print(m_hat)
print(beta_hat)

5.130661252526928
101.0

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

$\displaystyle \begin{aligned} \hat{m} &= \frac{1}{\hat{\beta}} \left( \sum_{n=1}^N x_n + \beta m \right) \\ \hat{\beta} &= N + \beta \end{aligned} \tag{3.83}$

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

# λの事後分布のパラメータを計算:式(3.88)
a_hat = 0.5 * N + a
b_hat = 0.5 * (np.sum(x_n**2) + beta * m**2 - beta_hat * m_hat**2) + b
print(a_hat)
print(b_hat)

51.0
205.53419983252502

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

$\displaystyle \begin{aligned} \hat{a} &= \frac{N}{2} + a \\ \hat{b} &= \frac{1}{2} \left( \sum_{n=1}^N x_n^2 + \beta m^2 - \hat{\beta} \hat{m}^2 \right) + b \end{aligned} \tag{3.88}$

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

# μの事後分布の確率密度を計算:式(2.64)
N_dens_mat = norm.pdf(x=mu_mat, loc=m_hat, scale=1.0/np.sqrt(beta_hat*lambda_mat))
print(N_dens_mat[:5, :5])

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

# μ,λの結合事後分布の確率密度を計算:式(3.81)
posterior_dens_mat = N_dens_mat * Gam_dens_mat
print(posterior_dens_mat[:5, :5])

[[0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
  0.00000000e+00]
 [1.24258462e-23 9.42520186e-23 6.86609568e-22 4.80378482e-21
  3.22783809e-20]
 [5.44623265e-46 3.13347018e-44 1.66289163e-42 8.13975504e-41
  3.67508788e-39]
 [2.06726922e-68 9.02175653e-66 3.48777673e-63 1.19445490e-60
  3.62371762e-58]
 [7.39812434e-91 2.44895259e-87 6.89695102e-84 1.65253976e-80
  3.36871837e-77]]
[[0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
  0.00000000e+00]
 [3.81408498e-48 3.81408498e-48 3.81408498e-48 3.81408498e-48
  3.81408498e-48]
 [5.49878147e-34 5.49878147e-34 5.49878147e-34 5.49878147e-34
  5.49878147e-34]
 [4.48958010e-26 4.48958010e-26 4.48958010e-26 4.48958010e-26
  4.48958010e-26]
 [1.01512347e-20 1.01512347e-20 1.01512347e-20 1.01512347e-20
  1.01512347e-20]]
[[0.00000000e+000 0.00000000e+000 0.00000000e+000 0.00000000e+000
  0.00000000e+000]
 [4.73932332e-071 3.59485209e-070 2.61878724e-069 1.83220435e-068
  1.23112488e-067]
 [2.99476432e-079 1.72302677e-077 9.14387768e-076 4.47587341e-074
  2.02085051e-072]
 [9.28117076e-094 4.05038986e-091 1.56586530e-088 5.36260096e-086
  1.62689705e-083]
 [7.51000961e-111 2.48598924e-107 7.00125682e-104 1.67753189e-100
  3.41966506e-097]]

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

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

# 事後分布のラベルを作成
posterior_param_lbl  = f'$N = {N}, '
posterior_param_lbl += f'\\mu_{{truth}} = {mu_truth:.3g}, \\lambda_{{truth}} = {lambda_truth:.3g}, '
posterior_param_lbl += f'\\hat{{m}} = {m_hat:.3g}, \\hat{{\\beta}} = {beta_hat:.3g}, \\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('Gaussian-Gamma distribution', fontsize=20)
ax2x = ax.twiny() # 第2軸の設定用
ax2y = ax.twinx() # 第2軸の設定用

ax.plot(
    0.0, 0.0, 
    color=cm.viridis(X=0.0), linewidth=1.0, linestyle='-', 
    label='posterior distribution', zorder=10
) # (凡例表示用のダミー)
cs = ax.contourf(
    mu_mat, lambda_mat, posterior_dens_mat, 
    cmap='viridis', alpha=0.5, 
    zorder=11
) # 事後分布
ax.axvline(
    x=mu_truth, 
    color='red', linewidth=1.0, linestyle='--', 
    label='true parameter', zorder=12
) # 真のパラメータ
ax.axhline(
    y=lambda_truth, 
    color='red', linewidth=1.0, linestyle='--', 
    zorder=12
) # 真のパラメータ

ax2x.set_xticks(ticks=[mu_truth], labels=['$\\mu_{truth}$']) # パラメータラベル
ax2y.set_yticks(ticks=[lambda_truth], labels=['$\\lambda_{truth}$']) # パラメータラベル
ax.set_xlabel('$\\mu$')
ax.set_ylabel('$\\lambda$')
ax.set_title(posterior_param_lbl, loc='left') # パラメータラベル
fig.colorbar(cs, ax=ax, shrink=1.0, label='density')
ax.legend(prop={'size': 8})
ax.grid(zorder=0)
ax.set_xlim(xmin=mu_min, xmax=mu_max)   # (目盛の共通化用)
ax2x.set_xlim(xmin=mu_min, xmax=mu_max) # (目盛の共通化用)
ax.set_ylim(ymin=lambda_min, ymax=lambda_max)   # (目盛の共通化用)
ax2y.set_ylim(ymin=lambda_min, ymax=lambda_max) # (目盛の共通化用)

plt.show()

　2枚目の図は、パラメータの真の値から標準偏差1つ分の範囲を拡大したものである。

# 事後分布を作図
fig, ax = plt.subplots(
    figsize=(9, 6), dpi=100, facecolor='white', constrained_layout=True, 
    subplot_kw={'projection': '3d'}
)
fig.suptitle('Gaussian-Gamma distribution', fontsize=20)

ax.plot(
    [mu_truth, mu_truth], [lambda_min, lambda_max], [0.0, 0.0], 
    color='red', linewidth=1.0, linestyle='--', 
    label='true parameter', zorder=10
) # 真のパラメータ
ax.plot(
    [mu_min, mu_max], [lambda_truth, lambda_truth], [0.0, 0.0], 
    color='red', linewidth=1.0, linestyle='--', 
    zorder=10
) # 真のパラメータ
ax.contour(
    X=mu_mat, Y=lambda_mat, Z=posterior_dens_mat, offset=0.0, 
    cmap='viridis', linewidths=1.0, linestyles=':', 
    zorder=11
) # 事後分布
ax.plot_surface(
    X=mu_mat, Y=lambda_mat, Z=posterior_dens_mat, 
    cmap='viridis', alpha=0.6, 
    label='posterior distribution', zorder=12
) # 事後分布

ax.set_xlabel('$\\mu$')
ax.set_ylabel('$\\lambda$')
ax.set_zlabel('density')
ax.set_title(posterior_param_lbl, loc='left') # パラメータラベル
ax.legend(prop={'size': 8})

plt.show()

　真のパラメータを赤色の直線(破線)、事後分布(ガウス-ガンマ分布)を等高線や曲面(塗りつぶし)で示す。

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

予測分布の計算

　観測データ $\mathbf{X}$ から、未観測データ $x_{*}$ の予測分布(スチューデントのt分布) $p(x_{*} \mid \mathbf{X}, m, \beta, a, b) = \mathrm{St}(x_{*} \mid \hat{\mu}_s, \hat{\lambda}_{s}, \hat{\nu}_s)$ のパラメータ $\hat{\mu}_s, \hat{\lambda}_{s}, \hat{\nu}_s$ を求める。
　スチューデントのt分布については「【Python】1次元スチューデントのt分布の作図 - からっぽのしょこ」を参照のこと。

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

# 予測分布のパラメータを計算:式(3.95')
mu_s_hat     = m_hat
lambda_s_hat = beta_hat / (1.0+beta_hat) * a_hat / b_hat
nu_s_hat     = 2.0 * a_hat
print(mu_s_hat)
print(lambda_s_hat)
print(nu_s_hat)

5.130661252526928
0.24570120223859973
102.0

# 予測分布のパラメータを計算:式(3.95')
mu_s_hat      = (np.sum(x_n) + beta_hat * m) / (N + beta)
tmp_term      = np.sum(x_n**2) + beta * m**2 - (np.sum(x_n) + beta * m)**2 / (N + beta)
lambda_s_hat  = (N + beta) / (N + 1 + beta)
lambda_s_hat *= (0.5 * N + a)
lambda_s_hat /= 0.5 * tmp_term + beta
nu_s_hat      = N + 2.0 * a
print(mu_s_hat)
print(lambda_s_hat)
print(nu_s_hat)

5.130661252526928
0.24570120223859945
102.0

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

$\displaystyle \begin{aligned} \hat{\mu}_s &= \hat{m} \\ &= \frac{ \sum_{n=1}^N x_n + \beta m }{ N + \beta } \\ \hat{\lambda}_s &= \frac{ \hat{\beta} \hat{a} }{ (1 + \hat{\beta}) \hat{b} } \\ &= \frac{ (N + \beta) \left( \frac{N}{2} + a \right) }{ (N + 1 + \beta) \left[ \frac{1}{2} \left\{ \sum_{n=1}^N x_n^2 + \beta m^2 - \frac{1}{N + \beta} \Bigl( \sum_{n=1}^N x_n + \beta m \Bigr)^2 \right\} + b \right] } \\ \hat{\nu}_s &= 2 \hat{a} \\ &= N + 2 a \end{aligned} \tag{3.95'}$

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

# 予測分布の確率を計算:式(3.76)
predict_dens_vec = t.pdf(x=x_vec, df=nu_s_hat, loc=mu_s_hat, scale=1.0/np.sqrt(lambda_s_hat))
print(predict_dens_vec[:5])

[2.25938682e-06 2.35225721e-06 2.44882700e-06 2.54923825e-06
 2.65363820e-06]

　 $x$ の値ごとに、スチューデントのt分布に従う確率密度 $\mathrm{St}(x \mid \mu_s, \hat{\lambda}_s, \hat{\nu}_s)$ を計算する。
　スチューデントのt分布の確率密度関数 sp.stats.t.pdf() の確率変数の引数 x に $x$ 、自由度の引数 nu に $\hat{\nu}_s$ 、中心の引数 mu に $\hat{\mu}_s$ 、尺度の引数 sigma に $\hat{\sigma}_s = \frac{1}{\sqrt{\hat{\lambda}_s}}$ を指定する。

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

# 予測分布のラベルを作成
predict_param_lbl  = f'$N = {N}, '
predict_param_lbl += f'\\mu_{{truth}} = {mu_truth:.3g}, \\lambda_{{truth}} = {lambda_truth:.3g}, '
predict_param_lbl += f'\\hat{{\\mu}}_s = {mu_s_hat:.3g}, \\hat{{\\lambda}}_s = {lambda_s_hat:.3g}, \\hat{{\\nu}}_s = {nu_s_hat:.3g}$'

# 予測分布を作図
fig, ax = plt.subplots(
    figsize=(9, 6), dpi=100, facecolor='white', 
    constrained_layout=True
)
fig.suptitle("Student's t 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()

　真の分布(ガウス分布)を赤色の曲線(破線)、予測分布(スチューデントのt分布)を紫色の曲線(実線)で示す。

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

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

スポンサードリンク

学習の推移

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

データ数と分布の関係

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

　 $n$ 個のデータから求めた( $n$ 回更新した)事後分布(ガウス-ガンマ分布) $\mathcal{N}(\mu, \lambda \mid \hat{m}, \hat{\beta}, \hat{a}, \hat{b})$ を等高線や曲面、 $n$ 番目の観測データ $x_n$ に対応する位置 $(\mu, \lambda) = (x_n, (\frac{\sigma}{x_n - \mu})^2)$ を桃色の点で示す。

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

$\displaystyle \begin{aligned} m^{(n)} &\leftarrow \frac{1}{\hat{\beta}^{(n)}} \left( x_n + \beta^{(n-1)} m^{(n-1)} \right) \\ \beta^{(n)} &\leftarrow 1 + \beta^{(n-1)} \\ a^{(n)} &\leftarrow \frac{1}{2} + a^{(n-1)} \\ b^{(n)} &\leftarrow \frac{1}{2} \left( x_n^2 + \beta^{(n-1)} (m^{(n-1)})^2 - \beta^{(n)} (m^{(n)})^2 \right) + b^{(n-1)} \end{aligned}$

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

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

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

　 $n$ 個のデータから求めた( $n$ 回更新した)予測分布(スチューデントのt分布) $\mathrm{St}(x_{*} \mid \hat{\mu}_s, \hat{\lambda}_s^{-1}, \hat{\nu}_s)$ を紫色の曲線(実線)、 $n$ 番目の観測データ $x_n$ を桃色の点で示す。

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

$\displaystyle \begin{aligned} \mu_s^{(n)} &= m^{(n-1)} \\ \lambda_s^{(n)} &= \frac{\beta^{(n)} a^{(n)}}{(1 + \beta^{(n)}) b^{(n)}} \\ \nu_s^{(n)} &\leftarrow 1 + \nu_s^{(n-1)} \end{aligned}$

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

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

観測データと分布の関係

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

　真の分布の期待値 $\mathbb{E}[x_n]$ と真のパラメータ $\mu_{\mathrm{truth}}, \lambda_{\mathrm{truth}}$ の各図(分布)に対応する位置を赤色の垂直線など(破線)、観測データの標本平均 $\bar{x}$ を桃色の垂直線(破線)、事後分布の期待値 $\mathbb{E}[\mu], \mathbb{E}[\lambda]$ と予測分布の期待値 $\mathbb{E}[x_{*}]$ を紫色の垂直線など(破線)で示す。
　また、 $\sigma$ 軸と $\lambda$ 軸、 $\lambda$ 軸と $p(x)$ 軸の対応関係を恒等関数や変換曲線で示す。

　生成分布(ガウス分布)、事後分布(ガウス-ガンマ分布)、予測分布(スチューデントのt分布)の期待値は、それぞれ次の式で計算できる。

$\displaystyle \begin{aligned} \mathbb{E}_{\mathcal{N}(x_n \mid \mu_{\mathrm{truth}}, \lambda_{\mathrm{truth}}^{-1})}[x] &= \mu_{\mathrm{truth}} \\ \mathbb{E}_{\mathrm{NG}(\mu, \lambda \mid \hat{m}, \hat{\beta}, \hat{a}, \hat{b})}[\mu] &= \hat{m} \\ \mathbb{E}_{\mathrm{NG}(\mu, \lambda \mid \hat{m}, \hat{\beta}, \hat{a}, \hat{b})}[\lambda] &= \frac{\hat{a}}{\hat{b}} \\ &= \frac{1 + \hat{\beta}}{\hat{\beta}} \hat{\lambda}_s \\ \mathbb{E}_{\mathrm{St}(x_{*} \mid \hat{\mu}_s, \hat{\lambda}_s^{-1}, \hat{\nu}_s)}[x] &= \hat{\mu}_s \\ &= \hat{m} \end{aligned}$

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

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

　この記事では、1次元ガウス分布に対するベイズ推論を扱った。次の記事では、多次元ガウス分布を扱う。