からっぽのしょこ

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

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

はじめに

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

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

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

【数式読解編】

www.anarchive-beta.com

【他の節の内容】

www.anarchive-beta.com

【この節の内容】

・Pythonでやってみよう

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

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

# 3.4.3項で利用するライブラリ
import numpy as np
from scipy.stats import multivariate_normal, multivariate_t # 多次元ガウス分布, 多次元スチューデントのt分布
import matplotlib.pyplot as plt

 この例では、SciPyライブラリのstatsから多次元ガウス分布の確率密度関数multivariate_normal.pdf()と多次元スチューデントのt分布の確率密度関数multivariate_t.pdf()を使います。ただし、multivariate_tはバージョン1.6.0で追加された?もののようなのでSciPyをアップデートする必要があるかもしれません。

・観測モデルの構築

 まずは、観測モデルを設定します。この例では、尤度$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])
sigma_truth_dd = np.array([[20.0, 15.0], [15.0, 30.0]])
lambda_truth_dd = np.linalg.inv(sigma_truth_dd**2)

 平均パラメータ$\boldsymbol{\mu} = (\mu_1, \cdots, \mu_D)$をmu_truth_dとして値を指定します。この項では未知の値であり、この値を求めるのが目的です。
 分散共分散行列$\boldsymbol{\Sigma} = (\sigma_{1,1}^2, \cdots, \sigma_{D,D}^2)$ではなく、標準偏差$\sigma_{i,i}$と相関係数$\sigma_{i,j}$の行列$(\sigma_{1,1}, \cdots, \sigma_{D,D})$をsigma_truth_ddとして値を指定します。これは作図時に標準偏差を利用するためです。
 sigma_truth_ddから精度パラメータ(精度行列)$\boldsymbol{\Lambda} = \boldsymbol{\Sigma}^{-1}$を計算してlambda_truth_ddとします。逆行列はnp.linalg.inv()で計算します。

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

# 作図用のxの点を作成
x_1_point = np.linspace(
    mu_truth_d[0] - 4 * sigma_truth_dd[0, 0], 
    mu_truth_d[0] + 4 * sigma_truth_dd[0, 0], 
    num=1000
)
x_2_point = np.linspace(
    mu_truth_d[1] - 4 * sigma_truth_dd[1, 1], 
    mu_truth_d[1] + 4 * sigma_truth_dd[1, 1], 
    num=1000
)
x_1_grid, x_2_grid = np.meshgrid(x_1_point, x_2_point)
x_point_arr = np.stack([x_1_grid.flatten(), x_2_grid.flatten()], axis=1)
x_dims = x_1_grid.shape
print(x_dims)
(1000, 1000)

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

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

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

 x_point_arrの値(の組み合わせ)ごとに確率密度を計算します。多次元ガウス分布の確率密度は、multivariate_normal.pdf()で計算できます。データの引数xx_point_arr、平均の引数meanmu_truth_d、分散共分散行列の引数covnp.linalg.inv(lambda_dd)を指定します。$\boldsymbol{\Sigma} = \boldsymbol{\Lambda}^{-1}$の計算をして精度行列から分散共分散行列に戻しています。sigma_truth_dd^2を指定することもできます。

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

# 確認
print(x_point_arr)
print(true_model)
[[-55.         -70.        ]
 [-54.83983984 -70.        ]
 [-54.67967968 -70.        ]
 ...
 [104.67967968 170.        ]
 [104.83983984 170.        ]
 [105.         170.        ]]
[2.52911750e-09 2.58863094e-09 2.64934710e-09 ... 2.64934710e-09
 2.58863094e-09 2.52911750e-09]

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

 尤度を作図します。

# 尤度を作図
plt.figure(figsize=(12, 9))
plt.contour(x_1_grid, x_2_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_truth_dd, 5)]) + '$', loc='left')
plt.colorbar()
plt.show()

f:id:anemptyarchive:20210410212337p: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_truth_dd), size=N
)

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

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

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

# 確認
print(x_nd[:5])
[[28.89434754 37.19128405]
 [12.02477961 52.19017493]
 [-0.58315191 78.15665593]
 [15.32763575 35.48928928]
 [19.75164614 39.96216006]]


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

# 観測データの散布図を作成
plt.figure(figsize=(12, 9))
plt.scatter(x=x_nd[:, 0], y=x_nd[:, 1]) # 観測データ
plt.contour(x_1_grid, x_2_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]) + ']' + 
          ', \Sigma=' + str([list(lmd_d) for lmd_d in np.round(np.sqrt(np.linalg.inv(lambda_truth_dd)), 1)]) + 
          '$', loc='left')
plt.colorbar()
plt.show()

f:id:anemptyarchive:20210410212352p:plain
観測データの散布図:多次元ガウス分布

 plt.scatter()で散布図を描画します。(なぜか$\Sigma_{\mu}$と表示されていますが$\Sigma$です。コードは直しましたが、画像の差し替えは非常に面倒なので見逃してください。)

・事前分布の設定

 次に、尤度に対する共役事前分布を設定します。多次元ガウス分布の平均パラメータ$\boldsymbol{\mu}$と精度パラメータ$\boldsymbol{\Lambda}$に対する事前分布$p(\boldsymbol{\mu}, \boldsymbol{\Lambda})$として、ガウス・ウィシャート分布$\mathrm{NW}(\boldsymbol{\mu}, \boldsymbol{\Lambda} | \mathbf{m}, \beta, \nu, \mathbf{W})$を設定します。

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

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

# lambdaの事前分布のパラメータを指定
w_dd = np.array([[0.0005, 0], [0, 0.0005]])
nu = 2

 多次元ガウス分布の平均パラメータ$\mathbf{m} = (m_1, \cdots, m_D)$をm_d、精度パラメータの係数$\beta$をbetaとして値を指定します。

 ウィシャート分布のパラメータ$\mathbf{W} = (w_{1,1}, \cdots, w_{D,D})$をW_dd、自由度$\nu$をnuとして値を指定します。$\mathbf{W}$は正定値行列、$\nu$は$D - 1$より大きい必要があります。

 事前分布のパラメータを設定できたので、これまでのように事前分布をグラフで確認しましょう。しかし、この例では精度パラメータ$\boldsymbol{\Lambda}$が未知の値なので、$\boldsymbol{\mu}$の事前分布$\mathcal{N}(\boldsymbol{\mu} | \mathbf{m}, (\beta \boldsymbol{\Lambda})^{-1})$の精度パラメータ$\beta \boldsymbol{\Lambda}$が求まりません。そこで、$\boldsymbol{\Lambda}$の事前分布$\mathcal{W}(\boldsymbol{\Lambda} | \nu, \mathbf{W})$から精度パラメータの期待値$\mathbb{E}[\boldsymbol{\Lambda}]$を求めて、それを利用することにします。

# lambdaの期待値を計算:式(2.89)
E_lambda_dd = nu * w_dd

 精度パラメータの期待値は、ウィシャート分布の期待値(2.89)より

$$ \mathbb{E}[\boldsymbol{\Lambda}] = \nu \mathbf{W} $$

で計算できます。

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

# 作図用のmuの点を作成
mu_1_point = np.linspace(mu_truth_d[0] - 100.0, mu_truth_d[0] + 100.0, num=1000)
mu_2_point = np.linspace(mu_truth_d[1] - 100.0, mu_truth_d[1] + 100.0, num=1000)
mu_1_grid, mu_2_grid = np.meshgrid(mu_1_point, mu_2_point)
mu_point_arr = np.stack([mu_1_grid.flatten(), mu_2_grid.flatten()], axis=1)
mu_dims = mu_1_grid.shape
print(mu_dims)

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

 尤度のときと同様に、多次元ガウス分布に従う変数$\boldsymbol{\mu}$がとり得る値を作成してmu_point_arrとします。この例では、真の値を中心に指定した範囲とします(本当は自動で良い感じの範囲になるようにしたかった)。

 m_dbeta, E_lambda_ddを用いて、尤度のときと同様にして確率密度を計算します。

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

# 確認
print(mu_point_arr)
print(prior_mu)
[[-75.        -50.       ]
 [-74.7997998 -50.       ]
 [-74.5995996 -50.       ]
 ...
 [124.5995996 150.       ]
 [124.7997998 150.       ]
 [125.        150.       ]]
[2.73841206e-06 2.77978389e-06 2.82166767e-06 ... 8.80609199e-13
 8.58897023e-13 8.37686604e-13]


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

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

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

 $\boldsymbol{\mu}$の真の値を確率密度のグラフと重ねて表示します。(こちらも$\Lambda_{\mu}$ではなく$E[\Lambda]$ですね。)

 $\boldsymbol{\Lambda}$についても事前分布(の確率密度)を計算してグラフで確認したいところですが、精度行列の分布はイメージしにくいため、平均パラメータの期待値$\mathbb{E}[\boldsymbol{\mu}]$と精度パラメータの期待値$\mathbb{E}[\boldsymbol{\Lambda}]$を用いた多次元ガウス分布を尤度と比較することにしましょう。ちなみに、ウィシャート分布の確率密度はwishart.pdf()で計算できます(たぶん)。

# muの期待値を計算:式(2.76)
E_mu_d = m_d

 平均パラメータの期待値は、ガウス分布の期待値(3.76)、より

$$ \mathbb{E}[\boldsymbol{\mu}] = \mathbf{m} $$

です。

 $\boldsymbol{\mu}$の事前分布の期待値(平均パラメータの期待値)と$\boldsymbol{\Lambda}$の事前分布の期待値(精度パラメータの期待値)を用いて、多次元ガウス分布の多次元ガウス分布の確率密度を計算します。

# 事前分布の期待値を用いた分布を計算:式(2.72)
prior_lambda = multivariate_normal.pdf(
    x=x_point_arr, mean=E_mu_d, cov=np.linalg.inv(beta * E_lambda_dd)
)

 E_mu_d, E_lambda_ddを用いて、尤度のときと同様にして計算します。

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

# 確認
print(prior_lambda)
[3.02641337e-06 3.05315107e-06 3.08004599e-06 ... 3.52209575e-13
 3.46349382e-13 3.40577958e-13]


 $\boldsymbol{\mu}$と$\boldsymbol{\Lambda}$の事前分布の期待値(平均と精度パラメータの期待値)を用いた分布を作図します。

# 事前分布の期待値を用いた分布を作図
plt.figure(figsize=(9, 9))
plt.contour(x_1_grid, x_2_grid, prior_lambda.reshape(x_dims)) # 事前分布の期待値を用いた分布
plt.contour(x_1_grid, x_2_grid, true_model.reshape(x_dims), 
            alpha=0.5, linestyles='--') # 尤度
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
plt.suptitle('Multivariate Gaussian Distribution', fontsize=20)
plt.title('$\\nu=' + str(nu) + 
          ', W=' + str([list(w_d) for w_d in np.round(w_dd, 5)]) + 
          '$', loc='left')
plt.colorbar()
plt.show()

f:id:anemptyarchive:20210410212514p:plain
事前分布の期待値を用いた分布:多次元ガウス分布

 尤度を破線で重ねて描画します。

・事後分布の計算

 観測データ$\mathbf{X}$からパラメータ$\boldsymbol{\mu},\ \boldsymbol{\Lambda}$の事後分布$p(\boldsymbol{\mu}, \boldsymbol{\Lambda} | \mathbf{X})$を求めます(パラメータ$\boldsymbol{\mu},\ \boldsymbol{\Lambda}$を分布推定します)。事後分布はガウス・ウィシャート分布$\mathrm{NW}(\boldsymbol{\mu} | \hat{\mathbf{m}}, \hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1}, \hat{\nu}, \hat{\mathbf{W}})$になります。

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

# muの事後分布のパラメータを計算:式(3.129)
beta_hat = N + beta
m_hat_d = (np.sum(x_nd, axis=0) + beta * m_d) / beta_hat

 事後分布のパラメータを

$$ \begin{aligned} \hat{\beta} &= N + \beta \\ \hat{\mathbf{m}} &= \frac{1}{\hat{\beta}} \left( \sum_{n=1}^N \mathbf{x}_n + \beta \mathbf{m} \right) \end{aligned} \tag{3.129} $$

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

# 確認
print(beta_hat)
print(m_hat_d)
51
[22.12609475 52.19317486]


 $\boldsymbol{\Lambda}$の事後分布のパラメータも計算します。

# lambdaの事後分布のパラメータを計算:式(3.133)
term_x_dd = np.dot(x_nd.T, x_nd)
term_m_dd = beta * np.dot(m_d.reshape([2, 1]), m_d.reshape([1, 2]))
term_m_hat_dd = beta_hat * np.dot(m_hat_d.reshape([2, 1]), m_hat_d.reshape([1, 2]))
w_hat_dd = np.linalg.inv(
  term_x_dd + term_m_dd - term_m_hat_dd + np.linalg.inv(w_dd)
)
nu_hat = N + nu

 事後分布のパラメータを

$$ \begin{aligned} \hat{\mathbf{W}}^{-1} &= \sum_{n=1}^N \mathbf{x}_n \mathbf{x}_n^{\top} + \beta \mathbf{m} \mathbf{m}^{\top} - \hat{\beta} \hat{\mathbf{m}} \hat{\mathbf{m}}^{\top} + \mathbf{W}^{-1} \\ \hat{\nu} &= N + \nu \end{aligned} \tag{3.133} $$

で計算して、結果をw_hat_ddnu_hatとします。ただし、$\sum_{n=1}^N \mathbf{x}_n \mathbf{x}_n^{\top}$の計算を効率よく処理するために、$\mathbf{X}^{\top} \mathbf{X}$で計算しています。

# 確認
print(w_hat_dd)
print(nu_hat)
[[ 7.31498011e-05 -2.10869750e-05]
 [-2.10869750e-05  3.34901788e-05]]
52


 事前分布のときと同様に、求めたパラメータを用いて、学習後の精度パラメータの期待値$\mathbb{E}[\hat{\boldsymbol{\Lambda}}]$を計算します。

# lambdaの期待値を計算:式(2.89)
E_lambda_hat_dd = nu_hat * w_hat_dd

# 確認
print(E_lambda_hat_dd)
[[ 0.00380379 -0.00109652]
 [-0.00109652  0.00174149]]


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

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

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

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

# 確認
print(posterior_mu)
[0. 0. 0. ... 0. 0. 0.]


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

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

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

 $\boldsymbol{\mu}$の真の値付近をピークとする分布を推定できています。

 グラフの表示範囲が狭めて確認してみましょう。

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

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

 先のとがった分布になっているのが分かります。

 $\boldsymbol{\Lambda}$の事後分布についても事前分布のときと同様に、パラメータの期待値を用いた分布を確認しましょう。学習後の平均パラメータの期待値$\mathbb{E}[\hat{\boldsymbol{\mu}}]$を計算します。

# muの期待値を計算:式(2.76)
E_mu_hat_d = m_hat_d

# 確認
print(E_mu_hat_d)
[22.12609475 52.19317486]


 $\boldsymbol{\mu}$と$\boldsymbol{\Lambda}$の事後分布の期待値(平均と精度パラメータの期待値)を用いて、多次元ガウス分布の確率密度を計算します。

# 事後分布の期待値を用いた分布を計算:式(2.72)
posterior_lambda = multivariate_normal.pdf(
    x=x_point_arr, mean=E_mu_hat_d, cov=np.linalg.inv(E_lambda_hat_dd)
)

 E_mu_hat_d, E_lambda_hat_ddを用いて、尤度のときと同様にして計算します。

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

# 確認
print(posterior_lambda)
[3.14195783e-10 3.22303763e-10 3.30588714e-10 ... 2.10298884e-10
 2.04154530e-10 1.98170360e-10]


 $\boldsymbol{\mu}$と$\boldsymbol{\Lambda}$の事後分布の期待値(平均と精度パラメータの期待値)を用いた分布を作図します。

# 事後分布の期待値を用いた分布を作図
plt.figure(figsize=(12, 9))
plt.contour(x_1_grid, x_2_grid, posterior_lambda.reshape(x_dims)) # 事後分布の期待値を用いた分布
plt.contour(x_1_grid, x_2_grid, true_model.reshape(x_dims), 
            alpha=0.5, linestyles='--') # 尤度
#plt.scatter(x=x_nd[:, 0], y=x_nd[:, 1]) # 観測データ
plt.scatter(x=mu_truth_d[0], y=mu_truth_d[1], color='red', s=100, marker='x') # 真のmu
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
plt.suptitle('Multivariate Gaussian Distribution', fontsize=20)
plt.title('$N=' + str(N) + ', \hat{\\nu}=' + str(nu_hat) + 
          ', \hat{W}=' + str([list(w_d) for w_d in np.round(w_hat_dd, 5)]) + 
          '$', loc='left')
plt.colorbar()
plt.show()

f:id:anemptyarchive:20210410212642p:plain
事後分布の期待値を用いた分布:多次元ガウス分布

 分布の形状が近づいているのを確認できます。

・予測分布の計算

 最後に、観測データ$\mathbf{X}$から未観測のデータ$\mathbf{x}_{*}$の予測分布$p(\mathbf{x}_{*} | \mathbf{X})$を求めます。予測分布はスチューデントのt分布$p(\mathbf{x}_{*} | \hat{\boldsymbol{\mu}}_s, \hat{\boldsymbol{\Lambda}}_s, \hat{\nu}_s)$になります。

 事後分布のパラメータ、または観測データと事前分布のパラメータを用いて、予測分布のパラメータを計算します。

# 次元数:(固定)
D = 2.0

# 予測分布のパラメータを計算:式(3.140')
mu_s_d = m_hat_d
lambda_s_hat_dd = (1.0 - D + nu_hat) * beta_hat / (1 + beta_hat) * w_hat_dd
nu_s_hat = 1.0 - D + nu_hat

 予測分布のパラメータを

$$ \begin{aligned} \hat{\boldsymbol{\mu}}_s &= \hat{\mathbf{m}} \\ &= \frac{1}{N + \beta} \left( \sum_{n=1}^N \mathbf{x}_n + \beta \mathbf{m} \right) \\ \hat{\boldsymbol{\Lambda}}_s &= \frac{ (1 - D + \hat{\nu}) \hat{\beta} }{ 1 + \hat{\beta} } \hat{\mathbf{W}} \\ &= \frac{ (1 - D + \hat{\nu}) (N + \beta) }{ N + 1 + \beta } \hat{\mathbf{W}} \\ \hat{\nu}_s &= 1 - D + \hat{\nu} \\ &= N + 1 - D + \nu \end{aligned} $$

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

 それぞれ上の式だと、事後分布のパラメータm_hat_d, beta_hatw_hat_dd, nu_hatで計算できます。下の式だと、観測データx_ndと事前分布のパラメータm_d, betaw_dd, nuで計算できます。

# 確認
print(mu_s_d)
print(lambda_s_hat_dd)
print(nu_s_hat)
[22.12609475 52.19317486]
[[ 0.0036589  -0.00105475]
 [-0.00105475  0.00167515]]
51.0

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

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

# 予測分布を計算:式(3.121)
predict = multivariate_t.pdf(
    x=x_point_arr, loc=mu_s_d, shape=np.linalg.inv(lambda_s_hat_dd), df=nu_s_hat
)

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

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

# 確認
print(predict)
[4.75852679e-09 4.83856513e-09 4.91968684e-09 ... 3.66368309e-09
 3.59397760e-09 3.52542480e-09]


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

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

f:id:anemptyarchive:20210410212741p:plain
予測分布:多次元スチューデントのt分布

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

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

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

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

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

# 利用するライブラリ
import numpy as np
from scipy.stats import multivariate_normal, multivariate_t # 多次元ガウス分布, 多次元スチューデントのt分布
import matplotlib.pyplot as plt
import matplotlib.animation as animation


・モデルの設定

# 真のパラメータを指定
mu_truth_d = np.array([25.0, 50.0])
sigma_truth_dd = np.array([[20.0, 15.0], [15.0, 30.0]])
lambda_truth_dd = np.linalg.inv(sigma_truth_dd**2)

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

# lambdaの事前分布のパラメータを指定
w_dd = np.array([[0.0005, 0], [0, 0.0005]])
inv_w_dd = np.linalg.inv(w_dd)
nu = 2

# lambdaの期待値を計算:式(2.89)
E_lambda_dd = nu * w_dd

# 初期値による予測分布のパラメータを計算:式(3.140)
mu_s_d = m_d
lambda_s_dd = (nu - 1.0) * beta / (1 + beta) * w_dd
nu_s = nu - 1.0


・推論処理

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

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

# 作図用のmuの点を作成
mu_1_point = np.linspace(mu_truth_d[0] - 100.0, mu_truth_d[0] + 100.0, num=1000)
mu_2_point = np.linspace(mu_truth_d[1] - 100.0, mu_truth_d[1] + 100.0, num=1000)
mu_1_grid, mu_2_grid = np.meshgrid(mu_1_point, mu_2_point)
mu_point_arr = np.stack([mu_1_grid.flatten(), mu_2_grid.flatten()], axis=1)
mu_dims = mu_1_grid.shape

# 作図用のxの点を作成
x_1_point = np.linspace(mu_truth_d[0] - 4 * sigma_truth_dd[0, 0], mu_truth_d[0] + 4 * sigma_truth_dd[0, 0], num=300)
x_2_point = np.linspace(mu_truth_d[1] - 4 * sigma_truth_dd[1, 1], mu_truth_d[1] + 4 * sigma_truth_dd[1, 1], num=300)
x_1_grid, x_2_grid = np.meshgrid(x_1_point, x_2_point)
x_point_arr = np.stack([x_1_grid.flatten(), x_2_grid.flatten()], axis=1)
x_dims = x_1_grid.shape

# 推移の記録用の受け皿を初期化
x_nd = np.empty((N, 2))
trace_beta = [beta]
trace_m = [list(m_d)]
trace_posterior = [
    multivariate_normal.pdf(
        x=mu_point_arr, mean=m_d, cov=np.linalg.inv(beta * E_lambda_dd)
    )
]
trace_mu_s = [list(mu_s_d)]
trace_lambda_s = [[list(lmd_d) for lmd_d in lambda_s_dd]]
trace_nu_s = [nu_s]
trace_predict = [
    multivariate_t.pdf(
        x=x_point_arr, loc=mu_s_d, shape=np.linalg.inv(lambda_s_dd), df=nu_s
    )
]

# ベイズ推論
for n in range(N):
    # 多次元ガウス分布に従うデータを生成
    x_nd[n] = np.random.multivariate_normal(
        mean=mu_truth_d, cov= np.linalg.inv(lambda_truth_dd)
    ).flatten()
    
    # muの事後分布のパラメータを:式(3.129)
    old_beta = beta
    old_m_d = m_d.copy()
    beta += 1
    m_d = (x_nd[n] + old_beta * m_d) / beta
    
    # lambdaの事後分布のパラメータを更新:式(1.33)
    term_x_dd = np.dot(x_nd[n].reshape([2, 1]), x_nd[n].reshape([1, 2]))
    old_term_m_dd = old_beta * np.dot(old_m_d.reshape([2, 1]), old_m_d.reshape([1, 2]))
    term_m_dd = beta * np.dot(m_d.reshape([2, 1]), m_d.reshape([1, 2]))
    inv_w_dd += term_x_dd + old_term_m_dd - term_m_dd
    nu += 1
    
    # lambdaの期待値を計算:式(2.89)
    E_lambda_dd = nu * np.linalg.inv(inv_w_dd)
    
    # muの事後分布を計算:式(2.72)
    trace_posterior.append(
        multivariate_normal.pdf(
            x=mu_point_arr, mean=m_d, cov=np.linalg.inv(beta * E_lambda_dd)
        )
    )
    
    # 予測分布のパラメータを更新:式(3.140)
    mu_s_d = m_d
    lambda_s_dd = (nu - 1.0) * beta / (1 + beta) * np.linalg.inv(inv_w_dd)
    nu_s = nu - 1.0
    
    # 予測分布を計算:式(3.121)
    trace_predict.append(
        multivariate_t.pdf(
            x=x_point_arr, loc=mu_s_d, shape=np.linalg.inv(lambda_s_dd), df=nu_s
        )
    )
    
    # 超パラメータを記録
    trace_beta.append(beta)
    trace_m.append(list(m_d))
    trace_mu_s.append(list(mu_s_d))
    trace_lambda_s.append([list(lmd_d) for lmd_d in lambda_s_dd])
    trace_nu_s.append(nu_s)
    
    # 動作確認
    #print('n=' + str(n + 1) + ' (' + str(np.round((n + 1) / N * 100, 1)) + '%)')

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

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

・事後分布の推移

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

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

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

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

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

・予測分布の推移

## 予測分布の推移をgif画像化

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

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

# 作図処理を関数として定義
def update_predict(n):
    # 前フレームのグラフを初期化
    plt.cla()
    
    # nフレーム目の予測分布を作図
    plt.contour(x_1_grid, x_2_grid, np.array(trace_predict[n]).reshape(x_dims)) # 予測分布
    plt.contour(x_1_grid, x_2_grid, true_model.reshape(x_dims), 
                alpha=0.5, linestyles='--') # 尤度
    plt.scatter(x=x_nd[:n, 0], y=x_nd[:n, 1]) # 観測データ
    plt.xlabel('$x_1$')
    plt.ylabel('$x_2$')
    plt.suptitle("Multivariate Student's t Distribution", fontsize=20)
    plt.title('$N=' + str(n) + ', \hat{\\nu}_s=' + str(trace_nu_s[n]) + 
              ', \hat{\mu}_s=[' + ', '.join([str(mu) for mu in np.round(trace_mu_s[n], 1)]) + ']' + 
              ', \hat{\Lambda}_s=' + str([list(lmd_d) for lmd_d in np.round(trace_lambda_s[n], 5)]) + 
              '$', loc='left')

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

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


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

f:id:anemptyarchive:20210410213958g:plain
予測分布の推移:多次元スチューデントのt分布

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

参考文献

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

おわりに

 この節の修正とPythonでの実装をしてる間に、SciPyに多次元版のスチューデントのt分布のモジュールが実装されて歓喜した。

 2021年4月12日は、元℃-ute・元Buono!の鈴木愛理さんの27歳のお誕生日です!おめでとうございます!

 ほんと素敵な方です。ぜひ見て!聴いて!ください。

【次節の内容】

www.anarchive-beta.com