からっぽのしょこ

読んだら書く!書いたら読む!同じ事は二度調べ(たく)ない

【Python】3.4.1:多次元ガウス分布の学習と予測:平均が未知の場合【緑ベイズ入門のノート】

はじめに

 『ベイズ推論による機械学習入門』の学習時のノートです。基本的な内容は「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」になります。「数式」と「プログラム」から理解するのが目標です。

 この記事は、3.4.1項の内容です。「尤度関数を平均が未知の多次元ガウス分布(多変量正規分布)」、「事前分布を多次元ガウス分布」とした場合の「パラメータの事後分布」と「未観測値の予測分布」の計算をPythonで実装します。

 省略してある内容等ありますので、本とあわせて読んでください。初学者な自分が理解できるレベルまで落として書き下していますので、分かる人にはかなりくどくなっています。同じような立場の人のお役に立てれば幸いです。

【数式読解編】

www.anarchive-beta.com

【他の節の内容】

www.anarchive-beta.com

【この節の内容】

・Pythonでやってみよう

 人工的に生成したデータを用いて、ベイズ推論を行ってみましょう。

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

# 3.4.1項で利用するライブラリ
import numpy as np
from scipy.stats import multivariate_normal # 多次元ガウス分布
import matplotlib.pyplot as plt

 SciPyライブラリのstatsから多次元ガウス分布のクラスmultivariate_normalを利用します。multivariate_normalの確率密度関数pdf()を使います。

・観測モデルの構築

 まずは、観測モデルを設定します。この例では、尤度$p(\mathbf{X} | \boldsymbol{\mu})$をデータごとに独立した多次元ガウス分布$\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}, \boldsymbol{\Lambda}^{-1})$とします。

 尤度のパラメータを設定します。この実装例では2次元のグラフで表現するため、$D = 2$のときのみ動作します。

# 真の平均パラメータを指定
mu_truth_d = np.array([25.0, 50.0])

# (既知の)分散共分散行列を指定
sigma2_dd = np.array([[600.0, 100.0], [100.0, 400.0]])

# (既知の精度)行列を計算
lambda_dd = np.linalg.inv(sigma2_dd)

 平均パラメータ$\boldsymbol{\mu} = (\mu_1, \cdots, \mu_D)$をmu_truth_dとして値を指定します。この例では未知の値であり、この値を求めるのが目的です。
 分散共分散行列$\boldsymbol{\Sigma} = (\sigma_{1,1}^2, \cdots, \sigma_{D,D}^2)$をsigma2_ddとして値を指定します。作図時に標準偏差に変換して利用します。$\sigma_{i,i}^2$は$i$次元方向の分散、$\sigma_{i,j}^2$は$i$次元と$j$次元方向の共分散です。
 分散共分散行列の逆行列を精度行列と呼びます。sigma2_ddから精度パラメータ(精度行列)$\boldsymbol{\Lambda} = \boldsymbol{\Sigma}^{-1}$を計算してlambda_ddとします。逆行列はnp.linalg.inv()で計算できます。

 グラフ用の点を作成します。

# 作図用のxのx軸の値を作成
x_0_line = np.linspace(
    mu_truth_d[0] - 3 * np.sqrt(sigma2_dd[0, 0]), 
    mu_truth_d[0] + 3 * np.sqrt(sigma2_dd[0, 0]), 
    num=500
)

# 作図用のxのx軸の値を作成
x_1_line = np.linspace(
    mu_truth_d[1] - 3 * np.sqrt(sigma2_dd[1, 1]), 
    mu_truth_d[1] + 3 * np.sqrt(sigma2_dd[1, 1]), 
    num=500
)

# 格子状のxの値を作成
x_0_grid, x_1_grid = np.meshgrid(x_0_line, x_1_line)

# xの点を作成
x_point_arr = np.stack([x_0_grid.flatten(), x_1_grid.flatten()], axis=1)
x_dims = x_0_grid.shape
print(x_dims)
(500, 500)

 作図用に、2次元ガウス分布に従う変数$\mathbf{x}_n = (x_{n,1}, x_{n,2})$がとり得る点(2値の組み合わせ)を作成します。
 $x_{n,1}$がとり得る値(x軸の値)をnp.linspace()で作成してx_0_lineとします。この例では、平均値を中心に標準偏差の3倍を範囲とします。np.linspace()を使うと指定した要素数で等間隔に切り分けます。np.arange()を使うと切り分ける間隔を指定できます。処理が重い場合は、この値を調整してください。$i$軸の標準偏差$\sigma_{i,i}$は、分散の平方根$\sigma_{i,i} = \sqrt{\sigma_{i,i}^2}$です。平方根はnp.sqrt()で計算できます。同様に、$x_{n,2}$がとり得る値(y軸の値)も作成してx_1_lineとします。
 x_0_linex_1_lineの要素の全ての組み合わせを持つ配列をnp.meshgrid()で作成してx_0_girdx_1_gridとします。これは、確率密度を等高線図にする際に格子状の点(2軸の全ての値が直交する点)を渡す必要があるためです。
 また、x_0_girdx_1_gridを1列に並べたものをnp.stack()で列方向に結合してx_point_arrとします。こちらは確率密度の計算に使います。

 作成した$\mathbf{x}$の点は次のようになります。

# 確認
print(x_point_arr[:5])
[[-48.48469228 -10.        ]
 [-48.19016446 -10.        ]
 [-47.89563663 -10.        ]
 [-47.60110881 -10.        ]
 [-47.30658098 -10.        ]]

 1列目がx軸の値、2列目がy軸の値に対応しています。

 尤度の確率密度を計算します。

# 尤度を計算:式(2.72)
true_model = multivariate_normal.pdf(
    x=x_point_arr, mean=mu_truth_d, cov=np.linalg.inv(lambda_dd)
)

 x_point_arrの値(の組み合わせ)ごとに確率密度を計算します。多次元ガウス分布の確率密度は、multivariate_normal.pdf()で計算できます。データの引数Xx_point_arr、平均の引数mumu_truth_d、分散共分散行列の引数sigmalambda_ddを逆行列$\boldsymbol{\Sigma} = \boldsymbol{\Lambda}^{-1}$に変換して指定します。勿論sigma2_ddを指定すればいいのですが、ここでは式に合わせてlambda_ddを使っています。

 計算結果を確認しましょう。

# 確認
print(x_point_arr)
print(true_model)
[[-48.48469228 -10.        ]
 [-48.19016446 -10.        ]
 [-47.89563663 -10.        ]
 ...
 [ 97.89563663 110.        ]
 [ 98.19016446 110.        ]
 [ 98.48469228 110.        ]]
[1.88323099e-07 1.94035442e-07 1.99890897e-07 ... 1.99890897e-07
 1.94035442e-07 1.88323099e-07]

 true_modelは、x_dimsを使ってx_1_grid, x_2_gridと同じ形状に変換してから作図に使います。

 尤度を作図します。

# 尤度を作図
plt.figure(figsize=(12, 9))
plt.contour(x_0_grid, x_1_grid, true_model.reshape(x_dims)) # 尤度
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
plt.suptitle('Multivariate Gaussian Distribution', fontsize=20)
plt.title('$\mu=[' + ', '.join([str(mu) for mu in mu_truth_d]) + ']' + 
          ', \Lambda=' + str([list(lmd_d) for lmd_d in np.round(lambda_dd, 5)]) + '$', 
          loc='left')
plt.colorbar()
plt.show()

f:id:anemptyarchive:20210518093651p:plain
尤度:多次元ガウス分布

 plt.contour()で等高線グラフを描画します。格子状の点を渡す必要があります。

 真のパラメータを求めることは、この真の分布を求めることを意味します。

・データの生成

 続いて、構築したモデルに従って観測データ$\mathbf{X} = \{\mathbf{x}_1, \mathbf{x}_2, \cdots, \mathbf{x}_N\}$を生成します。

 多次元ガウス分布に従う$N$個のデータをランダムに生成します。

# (観測)データ数を指定
N = 50

# 多次元ガウス分布に従うデータを生成
x_nd = np.random.multivariate_normal(
    mean=mu_truth_d, cov=np.linalg.inv(lambda_dd), size=N
)

 生成するデータ数$N$をNとして値を指定します。

 多次元ガウス分布に従う乱数は、np.random.multivariate_normal()で生成できます。試行回数の引数sizeNを指定します。他の引数についてはmultivariate_normal.pdf()と同じです。生成したN個のデータをx_ndとします。

 観測したデータを確認しましょう。

# 確認
print(x_nd[:5])
[[ -2.71217414  18.05276797]
 [ 11.71327129  16.92522305]
 [ 57.41772843  53.95034246]
 [-13.68735434  17.50134914]
 [ 50.39252388  74.86825902]]


 観測データの散布図を尤度と重ねて確認します。

# 観測データの散布図を作成
plt.figure(figsize=(12, 9))
plt.scatter(x=x_nd[:, 0], y=x_nd[:, 1]) # 観測データ
plt.contour(x_0_grid, x_1_grid, true_model.reshape(x_dims)) # 真の分布
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
plt.suptitle('Multivariate Gaussian Distribution', fontsize=20)
plt.title('$N=' + str(N) + ', \mu=[' + ', '.join([str(mu) for mu in mu_truth_d]) + ']' + 
          ', \Lambda=' + str([list(lmd_d) for lmd_d in np.round(lambda_dd, 5)]) + '$', 
          loc='left')
plt.colorbar()
plt.show()

f:id:anemptyarchive:20210518093721p:plain
観測データ:多次元ガウス分布

 plt.scatter()で散布図を描画します。

・事前分布の設定

 次に、尤度に対する共役事前分布を設定します。多次元ガウス分布の平均パラメータ$\boldsymbol{\mu}$に対する事前分布$p(\boldsymbol{\mu})$として、多次元ガウス分布$\mathcal{N}(\boldsymbol{\mu} | \mathbf{m}, \boldsymbol{\Lambda}_{\boldsymbol{\mu}}^{-1})$を設定します。

 $\boldsymbol{\mu}$の事前分布のパラメータ(超パラメータ)を設定します。

# muの事前分布の平均パラメータを指定
m_d = np.array([0.0, 0.0])

# muの事前分布の分散共分散行列を指定
sigma2_mu_dd = np.array([[1000.0, 0.0], [0.0, 1000.0]])

# muの事前分布の精度行列を計算
lambda_mu_dd = np.linalg.inv(sigma2_mu_dd)

 尤度のときと同様に、多次元ガウス分布の平均パラメータ$\mathbf{m} = (m_1, \cdots, m_D)$をm_d、分散共分散行列$\boldsymbol{\Sigma}_{\boldsymbol{\mu}} = (\sigma_{(\mu),1,1}^2, \cdots, \sigma_{(\mu),D,D}^2)$をsigma2_mu_ddとして値を指定します。また、sigma2_mu_ddの逆行列を精度行列lambda_mu_ddとします。

 $\boldsymbol{\mu}$の事前分布もグラフで確認しましょう。作図用の$\boldsymbol{\mu}$の点を作成します。

# 作図用のmuのx軸の値を作成
mu_0_line = np.linspace(
    mu_truth_d[0] - 2 *  np.sqrt(sigma2_mu_dd[0, 0]), 
    mu_truth_d[0] + 2 *  np.sqrt(sigma2_mu_dd[0, 0]), 
    num=500
)

# 作図用のmuのy軸の値を作成
mu_1_line = np.linspace(
    mu_truth_d[1] - 2 *  np.sqrt(sigma2_mu_dd[1, 1]), 
    mu_truth_d[1] + 2 *  np.sqrt(sigma2_mu_dd[1, 1]), 
    num=500
)

# 格子状の点を作成
mu_0_grid, mu_1_grid = np.meshgrid(mu_0_line, mu_1_line)

# muの点を作成
mu_point_arr = np.stack([mu_0_grid.flatten(), mu_1_grid.flatten()], axis=1)
mu_dims = mu_0_grid.shape
print(mu_dims)
(500, 500)

 $\mathbf{x}$の点のときと同様に、2次元ガウス分布に従う変数$\boldsymbol{\mu} = (\mu_1, \mu_2)$がとり得る値を作成してmu_point_matとします。この例では、真の値を中心に$\boldsymbol{\mu}$の事前分布の標準偏差の2倍を範囲とします。

 $\boldsymbol{\mu}$の事前分布の確率密度を計算します。

# muの事前分布を計算:式(2.72)
prior = multivariate_normal.pdf(
    x=mu_point_arr, mean=m_d, cov=np.linalg.inv(lambda_mu_dd)
)

 尤度のときと同様にして、mu_point_arrの値ごとに確率密度を計算します。

 計算結果は次のようになります。

# 確認
print(mu_point_arr)
print(prior)
[[-38.2455532  -13.2455532 ]
 [-37.99206401 -13.2455532 ]
 [-37.73857482 -13.2455532 ]
 ...
 [ 87.73857482 113.2455532 ]
 [ 87.99206401 113.2455532 ]
 [ 88.2455532  113.2455532 ]]
[7.01611475e-05 7.08423800e-05 7.15256308e-05 ... 5.56405454e-09
 5.44149678e-09 5.32129663e-09]


 $\boldsymbol{\mu}$の事前分布を作図します。

# muの事前分布を作図
plt.figure(figsize=(12, 9))
plt.scatter(x=mu_truth_d[0], y=mu_truth_d[1], color='red', s=100, marker='x') # 真の値
plt.contour(mu_0_grid, mu_1_grid, prior.reshape(mu_dims)) # muの事前分布
plt.xlabel('$\mu_1$')
plt.ylabel('$\mu_2$')
plt.suptitle('Multivariate Gaussian Distribution', fontsize=20)
plt.title('$m=[' + ', '.join([str(m) for m in m_d]) + ']' + 
          ', \Lambda_{\mu}=' + str([list(lmd_d) for lmd_d in np.round(lambda_mu_dd, 5)]) + '$', 
          loc='left')
plt.colorbar()
plt.show()

f:id:anemptyarchive:20210518093741p:plain
平均パラメータの事前分布:多次元ガウス分布

 $\boldsymbol{\mu}$の真の値を確率密度のグラフと重ねて表示します。

・事後分布の計算

 観測データ$\mathbf{X}$から平均パラメータ$\boldsymbol{\mu}$の事後分布$p(\boldsymbol{\mu} | \mathbf{X})$を求めます(平均パラメータ$\boldsymbol{\mu}$を分布推定します)。事後分布は多次元ガウス分布$\mathcal{N}(\boldsymbol{\mu} | \hat{\mathbf{m}}, \hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1})$になります。

 観測データx_ndを用いて、$\boldsymbol{\mu}$の事後分布のパラメータを計算します。

# muの事後分布のパラメータを計算:式(3.102),(3.103)
lambda_mu_hat_dd = N * lambda_dd + lambda_mu_dd
term_x_d = np.dot(lambda_dd, np.sum(x_nd, axis=0))
term_m_d = np.dot(lambda_mu_dd, m_d)
m_hat_d = np.dot(np.linalg.inv(lambda_mu_hat_dd), (term_x_d + term_m_d))

 事後分布のパラメータを

$$ \begin{align} \hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}} &= N \boldsymbol{\Lambda} + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \tag{3.102}\\ \hat{\mathbf{m}} &= \hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1} \left( \boldsymbol{\Lambda} \sum_{n=1}^N \mathbf{x}_n + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} \right) \tag{3.103} \end{align} $$

で計算して、結果をlambda_mu_hat_ddm_hat_dとします。

# 確認
print(m_hat_d)
print(lambda_mu_dd)
[23.79214395 49.72868583]
[[0.001 0.   ]
 [0.    0.001]]


 求めたパラメータを用いて、$\boldsymbol{\mu}$の事後分布の確率密度を計算します。

# muの事後分布を計算:式(2.72)
posterior = multivariate_normal.pdf(
    x=mu_point_arr, mean=m_hat_d, cov=np.linalg.inv(lambda_mu_hat_dd)
)

 更新した超パラメータm_hat_d, lambda_mu_hat_ddを用いて、事前分布のときと同様にして計算します。

 計算結果は次のようになります。

# 確認
print(posterior)
[2.60363564e-152 7.31735610e-152 2.04490731e-151 ... 2.11929381e-157
 7.20740505e-158 2.43731771e-158]


 $\boldsymbol{\mu}$の事後分布を作図します。

# muの事後分布を作図
plt.figure(figsize=(12, 9))
plt.scatter(x=mu_truth_d[0], y=mu_truth_d[1], color='red', s=100, marker='x') # 真の値
plt.contour(mu_0_grid, mu_1_grid, posterior.reshape(mu_dims)) # muの事後分布
plt.xlabel('$\mu_1$')
plt.ylabel('$\mu_2$')
plt.suptitle('Multivariate Gaussian Distribution', fontsize=20)
plt.title('$N=' + str(N) + 
          ', \hat{m}=[' + ', '.join([str(m) for m in np.round(m_hat_d, 1)]) + ']' + 
          ', \hat{\Lambda}_{\mu}=' + str([list(lmd_d) for lmd_d in np.round(lambda_mu_hat_dd, 5)]) + '$', 
          loc='left')
plt.colorbar()
plt.show()

f:id:anemptyarchive:20210518093807p:plain
平均パラメータの事後分布:多次元ガウス分布

 $\boldsymbol{\mu}$の真の値付近をピークとする分布を推定できています。分布が小さくなったわけではなく、かなり先の尖った形をしています。

・予測分布の計算

 最後に、観測データ$\mathbf{X}$から未観測のデータ$\mathbf{x}_{*}$の予測分布$p(\mathbf{x}_{*} | \mathbf{X})$を求めます。予測分布は多次元ガウス分布$\mathcal{N}(\mathbf{x}_{*} | \hat{\boldsymbol{\mu}}_{*}, \hat{\boldsymbol{\Lambda}}_{*}^{-1})$になります。

 $\boldsymbol{\mu}$の事後分布のパラメータを用いて、予測分布のパラメータを計算します。

# 予測分布のパラメータを計算:式(3.109'),(3.110')
lambda_star_hat_dd = np.linalg.inv(
    np.linalg.inv(lambda_dd) + np.linalg.inv(lambda_mu_hat_dd)
)
mu_star_hat_d = m_hat_d

 予測分布のパラメータを

$$ \begin{aligned} \hat{\boldsymbol{\Lambda}}_{*} &= ( \boldsymbol{\Lambda}^{-1} + \hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1} )^{-1} \\ \hat{\boldsymbol{\mu}}_{*} &= \hat{\mathbf{m}} \end{aligned} $$

で計算して、結果をlambda_star_hat_ddm_star_hat_dとします。

# 確認
print(lambda_star_hat_dd)
print(mu_star_hat_d)
[[ 0.00170541 -0.00042626]
 [-0.00042626  0.00255793]]
[23.79214395 49.72868583]

 $\mathbf{X}$から$\hat{\boldsymbol{\mu}}_{*},\ \hat{\boldsymbol{\Lambda}}_{*}$を学習しているのが式からも分かります。

 求めたパラメータを用いて、予測分布の確率密度を計算します。

# 予測分布を計算:式(2.72)
predict = multivariate_normal.pdf(
    x=x_point_arr, mean=mu_star_hat_d, cov=np.linalg.inv(lambda_star_hat_dd)
)

 尤度のときと同様に、x_point_arrの値ごとに確率密度を計算します。

 計算結果は次のようになります。

# 確認
print(predict)
[2.48593290e-07 2.55839313e-07 2.63257596e-07 ... 1.94062962e-07
 1.88378634e-07 1.82833756e-07]


 予測分布を尤度と重ねて作図します。

# 予測分布を作図
plt.figure(figsize=(12, 9))
plt.contour(x_0_grid, x_1_grid, true_model.reshape(x_dims), 
            alpha=0.5, linestyles='--') # 真の分布
plt.scatter(x=mu_truth_d[0], y=mu_truth_d[1], color='red', s=100, marker='x') # 真のmu
plt.scatter(x=x_nd[:, 0], y=x_nd[:, 1]) # 観測データ
plt.contour(x_0_grid, x_1_grid, predict.reshape(x_dims)) # 予測分布
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
plt.suptitle('Multivariate Gaussian Distribution', fontsize=20)
plt.title('$N=' + str(N) + 
          ', \hat{\mu}_{*}=[' + ', '.join([str(mu) for mu in np.round(mu_star_hat_d, 1)]) + ']' + 
          ', \hat{\Lambda}_{*}=' + str([list(lmd_d) for lmd_d in np.round(lambda_star_hat_dd, 5)]) + '$', 
          loc='left')
plt.colorbar()
plt.show()

f:id:anemptyarchive:20210518093923p:plain
未観測データの予測分布:多次元ガウス分布

 観測データが増えると、予測分布が真の分布に近づきます。

・おまけ:アニメーションで推移の確認

 animationモジュールを利用して、事後分布と予測分布の推移をアニメーション(gif画像)で確認するためのコードです。

・コード(クリックで展開)

 異なる点のみを簡単に解説します。

・ライブラリの読み込み

# 3.4.1項で利用するライブラリ
import numpy as np
from scipy.stats import multivariate_normal # 多次元ガウス分布
import matplotlib.pyplot as plt
import matplotlib.animation as animation


・モデルの設定

# 真の平均パラメータを指定
mu_truth_d = np.array([25.0, 50.0])

# (既知の)分散共分散行列を指定
sigma2_dd = np.array([[600.0, 100.0], [100.0, 400.0]])

# (既知の精度)行列を計算
lambda_dd = np.linalg.inv(sigma2_dd)

# muの事前分布の平均パラメータを指定
m_d = np.array([0.0, 0.0])

# muの事前分布の分散共分散行列を指定
sigma2_mu_dd = np.array([[1000.0, 0.0], [0.0, 1000.0]])

# muの事前分布の精度行列を計算
lambda_mu_dd = np.linalg.inv(sigma2_mu_dd)

# 初期値による予測分布のパラメータを計算:式(3.109),(3.110)
lambda_star_dd = np.linalg.inv(
    np.linalg.inv(lambda_dd) + np.linalg.inv(lambda_mu_dd)
)
mu_star_d = m_d

 事前分布のパラメータを用いて、予測分布のパラメータを計算しておきます。

・作図用の点の作成

# 作図用のmuのx軸の値を作成
mu_0_line = np.linspace(
    mu_truth_d[0] - 2 *  np.sqrt(sigma2_mu_dd[0, 0]), 
    mu_truth_d[0] + 2 *  np.sqrt(sigma2_mu_dd[0, 0]), 
    num=500
)

# 作図用のmuのy軸の値を作成
mu_1_line = np.linspace(
    mu_truth_d[1] - 2 *  np.sqrt(sigma2_mu_dd[1, 1]), 
    mu_truth_d[1] + 2 *  np.sqrt(sigma2_mu_dd[1, 1]), 
    num=500
)

# 格子状の点を作成
mu_0_grid, mu_1_grid = np.meshgrid(mu_0_line, mu_1_line)

# muの点を作成
mu_point_arr = np.stack([mu_0_grid.flatten(), mu_1_grid.flatten()], axis=1)
mu_dims = mu_0_grid.shape
# 作図用のxのx軸の値を作成
x_0_line = np.linspace(
    mu_truth_d[0] - 3 * np.sqrt(sigma2_dd[0, 0]), 
    mu_truth_d[0] + 3 * np.sqrt(sigma2_dd[0, 0]), 
    num=500
)

# 作図用のxのx軸の値を作成
x_1_line = np.linspace(
    mu_truth_d[1] - 3 * np.sqrt(sigma2_dd[1, 1]), 
    mu_truth_d[1] + 3 * np.sqrt(sigma2_dd[1, 1]), 
    num=500
)

# 格子状のxの値を作成
x_0_grid, x_1_grid = np.meshgrid(x_0_line, x_1_line)

# xの点を作成
x_point_arr = np.stack([x_0_grid.flatten(), x_1_grid.flatten()], axis=1)
x_dims = x_0_grid.shape


・推論処理

 各試行の結果をリストに格納していく必要があります。$\boldsymbol{\mu}$の事後分布をtrace_posterior、予測分布をtrace_predict、各パラメータをtrace_***として、初期値の結果を持つように作成しておきます。

# データ数(試行回数)を指定
N = 100

# 観測データの受け皿を作成
x_nd = np.empty((N, 2))

# 推移の記録用の受け皿を初期化
trace_m = [m_d]
trace_lambda_mu = [lambda_mu_dd]
trace_posterior = [
    multivariate_normal.pdf(
        x=mu_point_arr, mean=m_d, cov=np.linalg.inv(lambda_mu_dd)
    )
]
trace_mu_star = [mu_star_d]
trace_lambda_star = [lambda_star_dd]
trace_predict = [
    multivariate_normal.pdf(
        x=x_point_arr, mean=mu_star_d, cov=np.linalg.inv(lambda_star_dd)
    )
]

# ベイズ推論
for n in range(N):
    # 多次元ガウス分布に従うデータを生成
    x_nd[n] = np.random.multivariate_normal(
        mean=mu_truth_d, cov=np.linalg.inv(lambda_dd), size=1
    ).flatten()
    
    # muの事後分布のパラメータを更新:式(3.102),(3.102)
    old_lambda_mu_dd = lambda_mu_dd.copy()
    lambda_mu_dd += lambda_dd
    term_m_d = np.dot(lambda_dd, x_nd[n]) + np.dot(old_lambda_mu_dd, m_d)
    m_d = np.dot(np.linalg.inv(lambda_mu_dd), term_m_d)
    
    # muの事後分布(多次元ガウス分布)を計算:式(2.72)
    trace_posterior.append(
        multivariate_normal.pdf(
            x=mu_point_arr, mean=m_d, cov=np.linalg.inv(lambda_mu_dd)
        )
    )
    
    # 予測分布のパラメータを計算:式(3.109),(3.110)
    lambda_star_dd = np.linalg.inv(
        np.linalg.inv(lambda_dd) + np.linalg.inv(lambda_mu_dd)
    )
    mu_star_d = m_d
    
    # 予測分布を計算:式(2.72)
    trace_predict.append(
        multivariate_normal.pdf(
            x=x_point_arr, mean=mu_star_d, cov=np.linalg.inv(lambda_star_dd)
        )
    )
    
    # n回目の結果を記録
    trace_m.append(m_d)
    trace_lambda_mu.append(lambda_mu_dd)
    trace_mu_star.append(mu_star_d)
    trace_lambda_star.append(lambda_star_dd)
    
    # 動作確認
    #print('n=' + str(n + 1) + ' (' + str(np.round((n + 1) / N * 100, 1)) + '%)')

 観測された各データによってどのように学習する(分布が変化する)のかを確認するため、for文で1データずつ処理します。よって、データ数Nがイタレーション数になります。

 一度の処理で事後分布のパラメータを計算するのではなく、事前分布(1ステップ前の事後分布)に対して繰り返し観測データの情報を与えることでパラメータを更新(上書き)していきます。
 それに伴い、事後分布のパラメータの計算式(3.102-103)の$N$と$\sum_{n=1}^N \mathbf{x}_n$の計算は、ループ処理によって$N$回繰り返し1x_nd[n]を加えることで行います。$n$回目のループ処理のときには、$n-1$回分の1x_nd[n]が既にm_dlambda_mu_ddに加えられているわけです。
 ただし、事後分布のパラメータの計算において、更新前(1ステップ前)のパラメータを使うため、old_lambda_mu_ddとして値を一時的に保存しておきます。

・事後分布の推移

# 画像サイズを指定
fig = plt.figure(figsize=(9, 9))

# 作図処理を関数として定義
def update_posterior(n):
    # 前フレームのグラフを初期化
    plt.cla()
    
    # nフレーム目のmuの事後分布を作図
    plt.scatter(x=mu_truth_d[0], y=mu_truth_d[1], 
                color='red', s=100, marker='x') # 真の値
    plt.contour(mu_0_grid, mu_1_grid, trace_posterior[n].reshape(mu_dims)) # muの事後分布
    plt.xlabel('$\mu_1$')
    plt.ylabel('$\mu_2$')
    plt.suptitle('Multivariate Gaussian Distribution', fontsize=20)
    plt.title('$N=' + str(n) + 
              ', \hat{m}=[' + ', '.join([str(m) for m in np.round(trace_m[n], 1)]) + ']' + 
              ', \hat{\Lambda}_{\mu}=' + str([list(lmd_d) for lmd_d in np.round(trace_lambda_mu[n], 5)]) + '$', 
              loc='left')

# gif画像を作成
posterior_anime = animation.FuncAnimation(fig, update_posterior, frames=N + 1, interval=100)
posterior_anime.save("ch3_4_1_Posterior.gif")

 各フレーム(各試行)におけるパラメータの値をタイトルとして表示しています。ややこしければ、plt.title('n=' + str(n), loc='left')として試行回数だけ表示するだけでもそれっぽくなります。

 初期値(事前分布)を含むため、フレーム数の引数nframesN + 1です。

・予測分布の推移

# 尤度を計算:式(2.72)
true_model = multivariate_normal.pdf(
    x=x_point_arr, mean=mu_truth_d, cov=np.linalg.inv(lambda_dd)
)

# 画像サイズを指定
fig = plt.figure(figsize=(9, 9))

# 作図処理を関数として定義
def update_predict(n):
    # 前フレームのグラフを初期化
    plt.cla()
    
    # nフレーム目の予測分布を作図
    plt.contour(x_0_grid, x_1_grid, true_model.reshape(x_dims), 
                alpha=0.5, linestyles='--') # 真の分布
    plt.scatter(x=mu_truth_d[0], y=mu_truth_d[1], 
                color='red', s=100, marker='x') # 真のmu
    plt.scatter(x=x_nd[:n, 0], y=x_nd[:n, 1]) # 観測データ
    plt.contour(x_0_grid, x_1_grid, trace_predict[n].reshape(x_dims)) # 予測分布
    plt.xlabel('$x_1$')
    plt.ylabel('$x_2$')
    plt.suptitle('Multivariate Gaussian Distribution', fontsize=20)
    plt.title('$N=' + str(n) + 
              ', \hat{\mu}_{*}=[' + ', '.join([str(mu) for mu in np.round(trace_mu_star[n], 1)]) + ']' + 
              ', \hat{\Lambda}_{*}=' + str([list(lmd_d) for lmd_d in np.round(trace_lambda_star[n], 5)]) + '$', 
              loc='left')

# gif画像を作成
predict_anime = animation.FuncAnimation(fig, update_predict, frames=N + 1, interval=100)
predict_anime.save("ch3_4_1_Predict.gif")

 (よく理解していないので、animationの解説は省略…とりあえずこれで動きます……)


f:id:anemptyarchive:20210518094446g:plain
平均パラメータの事後分布の推移:多次元ガウス分布

 データが増えるに従って先が尖っていきます。

f:id:anemptyarchive:20210518094657g:plain
未観測データの予測分布の推移:多次元ガウス分布

 新たなデータによって平均(分布の中心)が推移しているのを確認できます。

参考文献

  • 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.

おわりに

 Pythonにもかなり慣れてきたので、そろそろNumPyとMatplotlibのpyplot以外のライブラリ?モジュール?にも触ってみたい。

【次節の内容】

www.anarchive-beta.com