はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。基本的な内容は「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」になります。「数式」と「プログラム」から理解するのが目標です。
この記事は、4.4.4項の内容です。「観測モデルを多次元ガウス混合分布(多変量正規混合分布)」、「事前分布をガウス・ウィシャート分布とディリクレ分布」とするガウス混合モデルに対する崩壊型ギブスサンプリング(周辺化ギブスサンプリング)をPythonで実装します。
省略してある内容等ありますので、本とあわせて読んでください。初学者な自分が理解できるレベルまで落として書き下していますので、分かる人にはかなりくどくなっています。同じような立場の人のお役に立てれば幸いです。
【数式読解編】
【他の節の内容】
【この節の内容】
・Pythonでやってみる
人工的に生成したデータを用いて、ベイズ推論を行ってみましょう。アルゴリズム4.7を参考に実装します。
利用するライブラリを読み込みます。
# 4.4.4項で利用するライブラリ import numpy as np from scipy.stats import multivariate_normal, multivariate_t # 多次元ガウス分布, 多次元スチューデントのt分布 import matplotlib.pyplot as plt
scipyライブラリのstatsモジュールから、多次元ガウス分布multivariate_normal、多次元スチューデントのt分布multivariate_tのクラスを利用します。確率密度の計算にpdf()を使います。
・観測モデルの設定
まずは、観測モデルを設定します。この節では、$K$個の観測モデルを多次元ガウス分布$\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Lambda}_k^{-1})$とします。また、各データに対する観測モデルの割り当てをカテゴリ分布$\mathrm{Cat}(\mathbf{s}_n | \boldsymbol{\pi})$によって行います。
観測モデルのパラメータを設定します。この実装例では2次元のグラフで表現するため、$D = 2$のときのみ動作します。
# 次元数を設定:(固定) D = 2 # クラスタ数を指定 K = 3 # K個の真の平均を指定 mu_truth_kd = np.array( [[5.0, 35.0], [-20.0, -10.0], [30.0, -20.0]] ) # K個の真の分散共分散行列を指定 sigma2_truth_kdd = np.array( [[[250.0, 65.0], [65.0, 270.0]], [[125.0, -45.0], [-45.0, 175.0]], [[210.0, -15.0], [-15.0, 250.0]]] )
$K$の平均パラメータ$\boldsymbol{\mu} = \{\boldsymbol{\mu}_1, \boldsymbol{\mu}_2, \cdots, \boldsymbol{\mu}_K\}$、$\boldsymbol{\mu}_k = (\mu_{k,1}, \mu_{k,2}, \cdots, \mu_{k,D})$をmu_truth_kd、分散共分散行列パラメータ$\boldsymbol{\Sigma} = \{\boldsymbol{\Sigma}_1, \boldsymbol{\Sigma}_2, \cdots, \boldsymbol{\Sigma}_K\}$、$\boldsymbol{\Sigma}_k = (\sigma_{k,1,1}^2, \cdots, \sigma_{k,D,D}^2)$をsigma2_truth_ddkとして、値を指定します。$\boldsymbol{\Sigma}_k$は正定値行列です。また、分散共分散行列の逆行列$\boldsymbol{\Lambda}_k = \boldsymbol{\Sigma}_k^{-1}$を精度行列と呼びます。
設定したパラメータは次のようになります。
# 確認 print(mu_truth_kd) print(sigma2_truth_kdd)
[[ 5. 35.]
[-20. -10.]
[ 30. -20.]]
[[[250. 65.]
[ 65. 270.]]
[[125. -45.]
[-45. 175.]]
[[210. -15.]
[-15. 250.]]]
混合比率を設定します。
# 真の混合比率を指定 pi_truth_k = np.array([0.45, 0.25, 0.3])
混合比率パラメータ$\boldsymbol{\pi} = (\pi_1, \pi_2, \cdots, \pi_K)$をpi_true_kとして、$\sum_{k=1}^K \pi_k = 1$、$\pi_k \geq 0$の値を指定します。ここで、$\pi_k$は各データがクラスタ$k$となる($s_{n,k} = 1$となる)確率を表します。
この3つのパラメータの値を観測データから求めるのがここでの目的です。
設定した観測モデルをグラフで確認しましょう。作図用の点を作成します。
# 作図用のx軸のxの値を作成 x_1_line = np.linspace( np.min(mu_truth_kd[:, 0] - 3 * np.sqrt(sigma2_truth_kdd[:, 0, 0])), np.max(mu_truth_kd[:, 0] + 3 * np.sqrt(sigma2_truth_kdd[:, 0, 0])), num=300 ) # 作図用のy軸のxの値を作成 x_2_line = np.linspace( np.min(mu_truth_kd[:, 1] - 3 * np.sqrt(sigma2_truth_kdd[:, 1, 1])), np.max(mu_truth_kd[:, 1] + 3 * np.sqrt(sigma2_truth_kdd[:, 1, 1])), num=300 ) # 作図用の格子状の点を作成 x_1_grid, x_2_grid = np.meshgrid(x_1_line, x_2_line) # 作図用のxの点を作成 x_point = np.stack([x_1_grid.flatten(), x_2_grid.flatten()], axis=1) # 作図用に各次元の要素数を保存 x_dim = x_1_grid.shape print(x_dim)
(300, 300)
作図用に、2次元ガウス分布に従う変数$\mathbf{x}_n = (x_{n,1}, x_{n,2})$がとり得る点(値の組み合わせ)を作成します。$x_{n,1}$がとり得る値(x軸の値)をnp.linspace()で作成してx_1_lineとします。この例では、$K$個のクラスタそれぞれで平均値から標準偏差の3倍を引いた値と足した値を計算して、その最小値から最大値までを範囲とします。
np.linspace()を使うと指定した要素数で等間隔に切り分けます。np.arange()を使うと切り分ける間隔を指定できます。処理が重い場合は、この値を調整してください。
同様に、$x_{n,2}$がとり得る値(y軸の値)も作成してx_2_lineとします。
np.meshgrid()を使って、x_1_lineとx_2_lineの要素の全ての組み合わせを持つような2次元配列に変換してそれぞれx_1_grid, x_2_gridとします。
これは、確率密度を等高線図にする際に格子状の点(2軸の全ての値が直交する点)を渡す必要があるためです。
また、x_1_grid, x_2_gridの要素を列として持つ2次元配列を作成してx_pointとします。これは、確率密度の計算に使います。
計算結果の確率密度は1次元配列になります。作図時にx_1_grid, x_2_gridと同じ形状にする必要があるため、形状をx_dimとして保存しておきます。
作成した$\mathbf{x}$の点は次のようになります。
# 確認 print(x_point[:5])
[[-53.54101966 -67.4341649 ]
[-53.11621983 -67.4341649 ]
[-52.69142 -67.4341649 ]
[-52.26662017 -67.4341649 ]
[-51.84182033 -67.4341649 ]]
1列目がx軸の値、2列目がy軸の値に対応しています。
観測モデルの確率密度を計算します。
# 観測モデルを計算 true_model = 0 for k in range(K): # クラスタkの分布の確率密度を計算 tmp_density = multivariate_normal.pdf( x=x_point, mean=mu_truth_kd[k], cov=sigma2_truth_kdd[k] ) # K個の分布の加重平均を計算 true_model += pi_truth_k[k] * tmp_density
$K$個のパラメータmu_truth_kd, sigmga2_truth_ddkを使って、x_point_matの点ごとに、次の混合ガウス分布の式で確率密度を計算します。
$K$個の多次元ガウス分布に対して、$\boldsymbol{\pi}$を使って加重平均をとります。
多次元ガウス分布の確率密度は、SciPyライブラリのmultivariate_normal.pdf()で計算できます。$k$番目の分布の場合は、データの引数Xにx_point、平均の引数muにmu_truth_kd[k]、分散共分散行列の引数sigmaにsigma2_truth_kdd[k]を指定します。
この計算をfor文でK回繰り返し行います。各パラメータによって求めた$K$回分の確率密度をそれぞれpi_truth_k[k]で割り引いてmodel_probに加算していきます。
計算結果は次のようになります。
# 確認 print(true_model[:5])
[9.06768334e-13 1.07345258e-12 1.27033352e-12 1.50260830e-12
1.77630256e-12]
観測モデルを作図します。
# 観測モデルを作図 plt.figure(figsize=(12, 9)) plt.contour(x_1_grid, x_2_grid, true_model.reshape(x_dim)) # 真の分布 plt.suptitle('Gaussian Mixture Model', fontsize=20) plt.title('K=' + str(K), loc='left') plt.xlabel('$x_1$') plt.ylabel('$x_2$') plt.colorbar() plt.show()
等高線グラフはplt.contour()で描画します。
$K$個の分布を重ねたような分布になります。ただし、積分して1となるようにしているため、$K$個の分布そのものよりも1つ１つの山が小さく(確率密度が小さく)なっています。
・データの生成
続いて、構築した観測モデルに従ってデータを生成します。先に、潜在変数$\mathbf{S} = \{\mathbf{s}_1, \mathbf{s}_2, \cdots, \mathbf{s}_N\}$、$\mathbf{s}_n = (s_{n,1}, s_{n,2}, \cdots, s_{n,K})$を生成し、各データにクラスタを割り当てます。次に、割り当てられたクラスタに従い観測データ$\mathbf{X} = \{\mathbf{x}_1, \mathbf{x}_2, \cdots, \mathbf{x}_N\}$、$\mathbf{x}_n = (x_{n,1}, x_{n,2}, \cdots, x_{n,D})$を生成します。
$N$個の潜在変数$\mathbf{S}$を生成します。
# (観測)データ数を指定 N = 250 # 潜在変数を生成 s_truth_nk = np.random.multinomial(n=1, pvals=pi_truth_k, size=N)
生成するデータ数$N$をNとして値を指定します。
カテゴリ分布に従う乱数は、多項分布に従う乱数生成関数np.random.multinomial()のn引数を1にすることで生成できます。また、確率(パラメータ)の引数probにpi_truth_k、試行回数の引数sizeにNを指定します。生成した$N$個の$K$次元ベクトルをs_truth_nkとします。
結果は次のようになります。
# 確認 print(s_truth_nk[:5])
[[1 0 0]
[1 0 0]
[0 1 0]
[0 1 0]
[1 0 0]]
これが本来は観測できない潜在変数であり、真のクラスタです。各行がデータを表し、値が1の列番号が割り当てられたクラスタ番号を表します。
s_truth_nkから各データのクラスタ番号を抽出します。
# クラスタ番号を抽出 _, s_truth_n = np.where(s_truth_nk == 1) # 確認 print(s_truth_n[:5])
[0 0 1 1 0]
np.where()を使って行(データ)ごとに要素が1の列番号を抽出します。
各データに割り当てられたクラスタに従い$N$個のデータを生成します。
# (観測)データを生成 x_nd = np.array([ np.random.multivariate_normal( mean=mu_truth_kd[k], cov=sigma2_truth_kdd[k], size=1 ).flatten() for k in s_truth_n ])
各データ$\mathbf{x}_n$に割り当てられたクラスタ($s_{n,k} = 1$である$k$)のパラメータ$\boldsymbol{\mu}_k,\ \boldsymbol{\Sigma}_k$を持つ多次元ガウス分布に従いデータ(値)を生成します。多次元ガウス分布の乱数は、np.random.multivariate_normaln()で生成できます。
データごとにクラスタが異なるので、for文(リスト内包表記)で1データずつ処理します。$n$番目のデータのクラスタ番号s_true_n[n]を取り出してkとします。そのクラスタに従ってデータ$\mathbf{x}_n$を生成します。平均の引数meanにmu_true_kd[k]、分散共分散行列の引数covにsigma2_true_kdd[k]を指定します。また1データずつ生成するので、データ数の引数sizeは1です。2次元配列で出力されるため、flatten()で1次元配列に変換しておきます。
生成した観測データは次のようになります。
# 確認 print(x_nd[:5])
[[ -9.02233866 5.09160069]
[-14.81289351 26.69130663]
[ 4.35088712 -36.27371768]
[-22.60485994 -41.78585905]
[-25.94071671 10.99908115]]
観測データの散布図を観測モデルに重ねて作図します。
# 観測データの散布図を作成 plt.figure(figsize=(12, 9)) for k in range(K): k_idx, = np.where(s_truth_n == k) # クラスタkのデータのインデック plt.scatter(x=x_nd[k_idx, 0], y=x_nd[k_idx, 1], label='cluster:' + str(k + 1)) # 観測データ plt.contour(x_1_grid, x_2_grid, true_model.reshape(x_dim), linestyles='--') # 真の分布 plt.suptitle('Gaussian Mixture Model', fontsize=20) plt.title('$N=' + str(N) + ', K=' + str(K) + ', \pi=[' + ', '.join([str(pi) for pi in pi_truth_k]) + ']$', loc='left') plt.xlabel('$x_1$') plt.ylabel('$x_2$') plt.colorbar() plt.show()
散布図はplt.scatter()で描画します。各クラスタが割り当てられたデータのインデックスをnp.where(s_true_n == k)で抽出してk_idxとします。k_idxを使ってクラスタ$k$が割り当てられたデータを取り出して、クラスタごとに色を変えた散布図を描きます。
観測データから各データのクラスタも推定します。
・事前分布の設定
次に、観測モデル(観測データの分布)$p(\mathbf{X} | \mathbf{S}, \boldsymbol{\mu}, \boldsymbol{\Lambda})$と潜在変数の分布$p(\mathbf{S} | \boldsymbol{\pi})$に対する共役事前分布を設定します。多次元ガウス分布のパラメータ$\boldsymbol{\mu},\ \boldsymbol{\Lambda}$に対する事前分布$p(\boldsymbol{\mu}, \boldsymbol{\Lambda})$としてガウス・ウィシャート分布$\mathrm{NW}(\boldsymbol{\mu}, \boldsymbol{\Lambda} | \mathbf{m}, \beta, \nu, \mathbf{W})$、カテゴリ分布のパラメータ$\boldsymbol{\pi}$に対する事前分布$p(\boldsymbol{\pi})$としてディリクレ分布$p(\boldsymbol{\pi} | \boldsymbol{\alpha})$を設定します。
各事前分布のパラメータ(超パラメータ)を設定します。
# muの事前分布のパラメータを指定 beta = 1.0 m_d = np.repeat(0.0, D) # lambdaの事前分布のパラメータを指定 w_dd = np.identity(D) * 0.0005 nu = D # piの事前分布のパラメータを指定 alpha_k = np.repeat(2.0, K)
$K$個の観測モデルの平均パラメータ$\boldsymbol{\mu}$の共通の事前分布(多次元ガウス分布)の平均パラメータ$\mathbf{m} = (m_1, m_2, \cdots, m_D)$をm_d、精度行列パラメータの係数$\beta$をbetaとして値を指定します。
$K$個の観測モデルの精度行列パラメータ$\boldsymbol{\Lambda}$の共通の事前分布(ウィシャート分布)のパラメータ$\mathbf{W} = (w_{1,1}, \cdots, w_{D,D})$をw_dd、自由度$\nu$をnuとして値を指定します。$\mathbf{W}$は正定値行列であり、自由度は$\nu > D - 1$の値をとります。
混合比率パラメータ$\boldsymbol{\pi}$の事前分布(ディリクレ分布)のパラメータ$\boldsymbol{\alpha} = (\alpha_1, \alpha_2, \cdots, \alpha_K)$をalpha_kとして、$\alpha_k > 0$の値を指定します。
パラメータを設定できたので、$\boldsymbol{\mu}$の事前分布をグラフで確認しましょう。作図用の$\boldsymbol{\mu}$の点を作成します。
# muの事前分布の標準偏差を計算 sigma_mu_d = np.sqrt( np.linalg.inv(beta * nu * w_dd).diagonal() ) # 作図用のx軸のmuの値を作成 mu_0_line = np.linspace( np.min(mu_truth_kd[:, 0]) - sigma_mu_d[0], np.max(mu_truth_kd[:, 0]) + sigma_mu_d[0], num=300 ) # 作図用のy軸のmuの値を作成 mu_1_line = np.linspace( np.min(mu_truth_kd[:, 1]) - sigma_mu_d[1], np.max(mu_truth_kd[:, 1]) + sigma_mu_d[1], num=300 ) # 作図用の格子状の点を作成 mu_0_grid, mu_1_grid = np.meshgrid(mu_0_line, mu_1_line) # 作図用のmuの点を作成 mu_point = np.stack([mu_0_grid.flatten(), mu_1_grid.flatten()], axis=1) # 作図用に各次元の要素数を保存 mu_dim = mu_0_grid.shape print(mu_dim)
(300, 300)
$\mathbf{x}$の点のときと同様に作成します。各軸の描画範囲は、真の平均の最小値から事前分布の標準偏差を引いた値から、真の平均の最大値に標準偏差を足した値までとします。事前分布の標準偏差は$(\beta \boldsymbol{\Lambda}_k)^{-1}$の対角成分の平方根です。ただしこの時点では$\boldsymbol{\Lambda}$が未知の値なので、代わりに$\boldsymbol{\Lambda}$の期待値$\mathbb{E}[\boldsymbol{\Lambda}]$を使うことにします。ウィシャート分布の期待値(2.89)より
です。
また、逆行列の計算はnp.linalg.inv()、対角成分の抽出はdiagonal()、平方根の計算はnp.sqrt()で行えます。
面倒でしたら適当な値を直接指定してください。
$\boldsymbol{\mu}$の事前分布の確率密度を計算します。
# muの事前分布を計算
prior_mu_density = multivariate_normal.pdf(
x=mu_point, mean=m_d, cov=np.linalg.inv(beta * nu * w_dd)
)
観測モデルのときと同様に、多次元ガウス分布の確率密度を計算します。ここでも、$\boldsymbol{\mu}$の事前分布の分散共分散行列を$(\beta \mathbb{E}[\boldsymbol{\Lambda}])^{-1}$とします。
計算結果は次のようになります。
# 確認 print(prior_mu_density[:5])
[1.10780663e-05 1.12959859e-05 1.15165400e-05 1.17397162e-05
1.19655008e-05]
$\boldsymbol{\mu}$の事前分布を作図します。
# muの事前分布を作図 plt.figure(figsize=(12, 9)) plt.scatter(x=mu_truth_kd[:, 0], y=mu_truth_kd[:, 1], label='true val', color='red', s=100, marker='x') # 真の平均 plt.contour(mu_0_grid, mu_1_grid, prior_mu_density.reshape(mu_dim)) # 事前分布 plt.suptitle('Gaussian Mixture Model', fontsize=20) plt.title('iter:' + str(0) + ', K=' + str(K), loc='left') plt.xlabel('$\mu_1$') plt.ylabel('$\mu_2$') plt.colorbar() plt.legend() plt.show()
真の平均値をプロットしておきます。この真の平均を分布推定します。
$\boldsymbol{\Lambda}$と$\boldsymbol{\pi}$の事前分布(と事後分布)については省略します。詳しくは、3.2.2 項「カテゴリ分布の学習と予測」と3.4.2項「多次元ガウス分布の学習と予測：精度が未知の場合」を参照してください。
・初期値の設定
設定した$\boldsymbol{\pi}$の事前分布$p(\boldsymbol{\pi})$に従い潜在変数$\mathbf{S}$を生成して初期値とします。また、生成した$\mathbf{S}$を用いて各パラメータ$\boldsymbol{\mu},\ \boldsymbol{\Lambda},\ \boldsymbol{\pi}$の($N$個全てのデータを用いた)事後分布$p(\boldsymbol{\mu} | \boldsymbol{\Lambda}, \mathbf{X}, \mathbf{S}),\ p(\boldsymbol{\Lambda} | \mathbf{X}, \mathbf{S}),\ p(\boldsymbol{\pi} | \mathbf{X}, \mathbf{S})$のパラメータ(超パラメータ)$\hat{\mathbf{m}},\ \hat{\boldsymbol{\beta}},\ \hat{\mathbf{W}},\ \hat{\boldsymbol{\nu}},\ \hat{\boldsymbol{\alpha}}$を計算して初期値とします。
alpha_kを用いて、$\mathbf{S}$の初期値を生成します。
# 潜在変数を初期化 s_nk = np.random.multinomial(n=1, pvals=alpha_k / np.sum(alpha_k), size=N)
真のクラスタのときと同様に生成してs_nkとします。ただし、ここでは各クラスタとなる確率として、$\boldsymbol{\alpha}$から求めた混合比率の期待値($\boldsymbol{\pi}$の事前分布の期待値)$\mathbb{E}[\boldsymbol{\pi}]$を使うことにします。クラスタ$k$となる確率の期待値$\mathbb{E}[\pi_k]$は、ディリクレ分布の期待値(2.51)より
です。
生成した潜在変数は次のようになります。
# 確認 print(s_nk[:5])
[[1 0 0]
[0 0 1]
[0 1 0]
[0 0 1]
[1 0 0]]
s_nkを用いて、$\boldsymbol{\mu}$の事後分布のパラメータを計算します。
# 初期値によるmuの事後分布のパラメータを計算:式(4.99) beta_hat_k = np.sum(s_nk, axis=0) + beta m_hat_kd = (np.dot(s_nk.T, x_nd) + beta * m_d) / beta_hat_k.reshape(K, 1)
$\boldsymbol{\mu}$の条件付き分布($K$個の多次元ガウス分布)の平均パラメータ$\hat{\mathbf{m}} = \{\hat{\mathbf{m}}_1, \hat{\mathbf{m}}_2, \cdots, \hat{\mathbf{m}}_K\}$、$\hat{\mathbf{m}}_k = (\hat{m}_{k,1}, \hat{m}_{k,2}, \cdots, \hat{m}_{k,D})$、精度行列パラメータの係数$\boldsymbol{\hat{\beta}} = \{\hat{\beta}_1, \hat{\beta}_2, \cdots, \hat{\beta}_K\}$の各項を
で計算して、それぞれbeta_hat_k、m_hat_kdとします。この時点では、$N$個全てのデータを使って計算しているため、(崩壊型でない)ギブスサンプリング(4.4.2項)の式(4.99)です。
# 確認 print(beta_hat_k) print(m_hat_kd)
[79. 91. 83.]
[[ 4.53341639 10.79193878]
[ 8.21529036 6.68526828]
[ 5.39221777 7.9168154 ]]
同様に、$\boldsymbol{\lambda}$の事後分布のパラメータを計算します。
# 初期値によるlambdaの事後分布のパラメータを計算:式(4.103) w_hat_kdd = np.zeros((K, D, D)) for k in range(K): inv_w_dd = np.dot(s_nk[:, k] * x_nd.T, x_nd) inv_w_dd += beta * np.dot(m_d.reshape((D, 1)), m_d.reshape((1, D))) inv_w_dd -= beta_hat_k[k] * np.dot(m_hat_kd[k].reshape((D, 1)), m_hat_kd[k].reshape((1, D))) inv_w_dd += np.linalg.inv(w_dd) w_hat_kdd[k] = np.linalg.inv(inv_w_dd) nu_hat_k = np.sum(s_nk, axis=0) + nu
$\boldsymbol{\lambda}$の条件付き分布($K$個のウィシャート分布)のパラメータ$\hat{\mathbf{W}} = \{\hat{\mathbf{W}}_1, \hat{\mathbf{W}}_2, \cdots, \hat{\mathbf{W}}_K\}$、$\hat{\mathbf{W}}_k = (\hat{w}_{k,1,1}, \cdots, \hat{w}_{k,D,D})$、自由度$\boldsymbol{\hat{\nu}} = \{\hat{\nu}_1, \hat{\nu}_2, \cdots, \hat{\nu}_K\}$を
で計算して、それぞれnu_hat_k、w_hat_ddkとします。
# 確認 print(w_hat_kdd) print(nu_hat_k)
[[[ 2.71849829e-05 -1.03335458e-06]
[-1.03335458e-06 1.19714702e-05]]
[[ 1.78771769e-05 1.53772085e-06]
[ 1.53772085e-06 1.28331770e-05]]
[[ 2.25146259e-05 4.28859168e-06]
[ 4.28859168e-06 1.20031226e-05]]]
[80 92 84]
$\boldsymbol{\pi}$の事後分布のパラメータを計算します。
# 初期値によるpiの事後分布のパラメータを計算:式(4.45) alpha_hat_k = np.sum(s_nk, axis=0) + alpha_k
$\boldsymbol{\pi}$の条件付き分布(ディリクレ分布)のパラメータ$\hat{\boldsymbol{\alpha}} = (\hat{\alpha}_1, \hat{\alpha}_2, \cdots, \hat{\alpha}_K)$の各項を
で計算して、alpha_hat_kとします。
# 確認 print(alpha_hat_k)
[80. 92. 84.]
超パラメータの初期値が得られたので、$\boldsymbol{\mu}$の事後分布を確認しましょう。$K$個の事後分布を計算します。
# 初期値によるmuの事後分布を計算 posterior_density_kg = np.empty((K, mu_point.shape[0])) for k in range(K): # クラスタkのmuの事後分布を計算 posterior_density_kg[k] = multivariate_normal.pdf( x=mu_point, mean=m_hat_kd[k], cov=np.linalg.inv(beta_hat_k[k] * nu_hat_k[k] * w_hat_kdd[k]) )
事前分布のときと同様に計算します。ただし、超パラメータが$K$個になったので、多次元ガウス分布も$K$個になります。$K$回分の計算結果を2次元配列に格納していきます。
クラスタ$k$の事後分布の平均は$\hat{\mathbf{m}}$で、精度行列として$\hat{\beta}_k \mathbb{E}[\boldsymbol{\Lambda}_k]$を使います。事前分布のときと同様に、$\mathbb{E}[\boldsymbol{\Lambda}_k] = \hat{\nu}_k \hat{\mathbf{W}}_k$です。
作成した配列は次のようになります。
# 確認 print(posterior_density_kg[:, :5])
[[3.51756671e-174 1.15039362e-172 3.67068377e-171 1.14273063e-169
3.47086043e-168]
[1.28177183e-217 5.00895620e-216 1.91584129e-214 7.17212097e-213
2.62791592e-211]
[8.73842381e-222 5.02919979e-220 2.82999329e-218 1.55701449e-216
8.37569075e-215]]
$\boldsymbol{\mu}$の事後分布を作図します。
# 初期値によるmuの事後分布を作図 plt.figure(figsize=(12, 9)) plt.scatter(x=mu_truth_kd[:, 0], y=mu_truth_kd[:, 1], label='true val', color='red', s=100, marker='x') # 真の平均 for k in range(K): plt.contour(mu_0_grid, mu_1_grid, posterior_density_kg[k].reshape(mu_dim)) # 事後分布 plt.suptitle('Gaussian Distribution', fontsize=20) plt.title('iter:' + str(0) + ', K=' + str(K), loc='left') plt.xlabel('$\mu_1$') plt.ylabel('$\mu_2$') plt.colorbar() plt.legend() plt.show()
事前分布のときと同様に作図できます。
m_hat_k, beta_hat_kとnu_hat_k, w_hat_ddkから求めた各パラメータの期待値を用いて、混合分布を計算します。
# 初期値による混合分布を計算 init_density = 0 for k in range(K): # クラスタkの分布の確率密度を計算 tmp_density = multivariate_normal.pdf( x=x_point, mean=m_hat_kd[k], cov=np.linalg.inv(nu_hat_k[k] * w_hat_kdd[k]) ) # K個の分布の加重平均を計算 init_density += alpha_hat_k[k] / np.sum(alpha_hat_k) * tmp_density
観測モデルのときと同様に、多次元ガウス混合分布の確率密度を計算します。ただし、パラメータ$\boldsymbol{\mu},\ \boldsymbol{\Lambda},\ \boldsymbol{\pi}$は未知なので、代わりに各パラメータの期待値$\mathbb{E}[\boldsymbol{\mu}],\ \mathbb{E}[\boldsymbol{\Lambda}],\ \mathbb{E}[\boldsymbol{\pi}]$を使うことにします。
クラスタ$k$における平均パラメータの期待値$\mathbb{E}[\boldsymbol{\mu}_k]$は、多次元ガウス分布の期待値(2.76)より
です。
また、混合比率パラメータの期待値は、先ほどと同様に$\mathbb{E}[\pi_k] = \frac{\hat{\alpha}_k}{\sum_{k'=1}^K \hat{\alpha}_{k'}}$です。
超パラメータの初期値から求めた分布を作図します。
# 初期値による分布を作図 plt.figure(figsize=(12, 9)) plt.contour(x_1_grid, x_2_grid, true_model.reshape(x_dim), alpha=0.5, linestyles='dashed') # 真の分布 plt.contour(x_1_grid, x_2_grid, init_density.reshape(x_dim)) # 期待値による分布:(等高線) plt.suptitle('Gaussian Mixture Model', fontsize=20) plt.title('iter:' + str(0) + ', K=' + str(K), loc='left') plt.xlabel('$x_1$') plt.ylabel('$x_2$') plt.colorbar() plt.show()
観測モデルのときと同様に作図します。
ここに、クラスタの初期値も重ねて確認しましょう。
# クラスタ番号を抽出 _, s_n = np.where(s_nk == 1) # クラスタの散布図を作成 plt.figure(figsize=(12, 9)) plt.contour(x_1_grid, x_2_grid, true_model.reshape(x_dim), alpha=0.5, linestyles='dashed') # 真の分布 #plt.contour(x_1_grid, x_2_grid, init_density.reshape(x_dim)) # 期待値による分布:(等高線) plt.contourf(x_1_grid, x_2_grid, init_density.reshape(x_dim), alpha=0.5) # 期待値による分布:(塗りつぶし) for k in range(K): k_idx, = np.where(s_n == k) # クラスタkのデータのインデックス plt.scatter(x=x_nd[k_idx, 0], y=x_nd[k_idx, 1], label='cluster:' + str(k + 1)) # サンプルしたクラスタ plt.suptitle('Gaussian Mixture Model', fontsize=20) plt.title('iter:' + str(0) + ', K=' + str(K), loc='left') plt.xlabel('$x_1$') plt.ylabel('$x_2$') plt.legend() plt.show()
データの生成のときと同様にして作図します。
・崩壊型ギブスサンプリング
崩壊型ギブスサンプリングでは、潜在変数$\mathbf{s}_n$のサンプリングと、各パラメータ$\boldsymbol{\mu},\ \boldsymbol{\Lambda},\ \boldsymbol{\pi}$の事後分布$p(\boldsymbol{\mu} | \boldsymbol{\Lambda}, \mathbf{X}_{\backslash n}, \mathbf{S}_{\backslash n}),\ p(\boldsymbol{\Lambda} | \mathbf{X}_{\backslash n}, \mathbf{S}_{\backslash n}),\ p(\boldsymbol{\pi} | \mathbf{X}_{\backslash n}, \mathbf{S}_{\backslash n})$のパラメータ(超パラメータ)$\hat{\mathbf{m}}_{\backslash n},\ \hat{\boldsymbol{\beta}}_{\backslash n},\ \hat{\mathbf{W}}_{\backslash n},\ \hat{\boldsymbol{\nu}}_{\backslash n},\ \hat{\boldsymbol{\alpha}}_{\backslash n}$の更新を交互に行います。
まずは、超パラメータの初期値を用いて1番目の潜在変数$\mathbf{s}_1$の事後分布$p(\mathbf{s}_1 | \mathbf{X}, \mathbf{S}_{\backslash 1})$を計算します。得られた事後分布から$\mathbf{s}_1$サンプルします。次に、1番目を更新した$\mathbf{S}$を用いて超パラメータを更新します。さらに、更新した超パラメータを用いて$\mathbf{s}_2$の事後分布$p(\mathbf{s}_2 | \mathbf{X}, \mathbf{S}_{\backslash 2})$を計算して$\mathbf{s}_2$をサンプルし、超パラメータを更新します。これを$N$個全てのデータで行います。ここまでが1回の試行です。
さらに次の試行では、前回の最後に更新した超パラメータを用いて$\mathbf{s}_1$の事後分布の計算となります。この処理を指定した回数繰り返します。
・コード全体
スチューデントのt分布の確率密度の変数density_st_kに関して、添字を使って代入するため予め変数を用意しておく必要があります。
また、後で各パラメータの更新値や分布の推移を確認するために、それぞれ試行ごとにtrace_***に値を保存していきます。関心のあるパラメータを記録してください。
# 試行回数を指定 MaxIter = 150 # 途中計算に用いる項の受け皿を作成 density_st_k = np.zeros(K) # 推移の確認用の受け皿を作成 _, s_n = np.where(s_nk == 1) # クラスタ番号を抽出 trace_s_in = [s_n.copy()] trace_beta_ik = [beta_hat_k.copy()] trace_m_ikd = [m_hat_kd.copy()] trace_w_ikdd = [w_hat_kdd.copy()] trace_nu_ik = [nu_hat_k.copy()] trace_alpha_ik = [alpha_hat_k.copy()] # 崩壊型ギブスサンプリング for i in range(MaxIter): for n in range(N): # n番目のデータの潜在変数を初期化 s_nk[n] = np.repeat(0, K) # muの事後分布のパラメータを計算:式(4.128) beta_hat_k = np.sum(s_nk, axis=0) + beta m_hat_kd = (np.dot(s_nk.T, x_nd) + beta * m_d) / beta_hat_k.reshape(K, 1) # lambdaの事後分布のパラメータを計算:式(4.128) term_m_dd = beta * np.dot(m_d.reshape((D, 1)), m_d.reshape((1, D))) for k in range(K): inv_w_dd = np.dot(s_nk[:, k] * x_nd.T, x_nd) inv_w_dd += term_m_dd inv_w_dd -= beta_hat_k[k] * np.dot(m_hat_kd[k].reshape((D, 1)), m_hat_kd[k].reshape((1, D))) inv_w_dd += np.linalg.inv(w_dd) w_hat_kdd[k] = np.linalg.inv(inv_w_dd) nu_hat_k = np.sum(s_nk, axis=0) + nu # piの事後分布のパラメータを計算:式(4.73) alpha_hat_k = np.sum(s_nk, axis=0) + alpha_k # スチューデントのt分布のパラメータを計算 nu_st_hat_k = 1 - D + nu_hat_k mu_st_hat_kd = m_hat_kd.copy() term_lmd_k = nu_st_hat_k * beta_hat_k / (1 + beta_hat_k) lambda_st_hat_kdd = term_lmd_k.reshape((K, 1, 1)) * w_hat_kdd # スチューデントのt分布の確率密度を計算:式(4.129) for k in range(K): density_st_k[k] = multivariate_t.pdf( x=x_nd[n], loc=mu_st_hat_kd[k], shape=np.linalg.inv(lambda_st_hat_kdd[k]), df=nu_st_hat_k[k] ) # カテゴリ分布のパラメータを計算:式(4.75) eta_k = alpha_hat_k / np.sum(alpha_hat_k) # 潜在変数のサンプリング確率を計算:式(4.124) tmp_p_k = density_st_k * eta_k prob_s_k = tmp_p_k / np.sum(tmp_p_k) # 正規化 # n番目のデータの潜在変数をサンプル s_nk[n] = np.random.multinomial(n=1, pvals=prob_s_k, size=1).flatten() # muの事後分布のパラメータを計算:式(4.99) beta_hat_k = np.sum(s_nk, axis=0) + beta m_hat_kd = (np.dot(s_nk.T, x_nd) + beta * m_d) / beta_hat_k.reshape(K, 1) # lambdaの事後分布のパラメータを計算:式(4.103) term_m_dd = beta * np.dot(m_d.reshape((D, 1)), m_d.reshape((1, D))) for k in range(K): inv_w_dd = np.dot(s_nk[:, k] * x_nd.T, x_nd) inv_w_dd += term_m_dd inv_w_dd -= beta_hat_k[k] * np.dot(m_hat_kd[k].reshape((D, 1)), m_hat_kd[k].reshape((1, D))) inv_w_dd += np.linalg.inv(w_dd) w_hat_kdd[k] = np.linalg.inv(inv_w_dd) nu_hat_k = np.sum(s_nk, axis=0) + nu # piの事後分布のパラメータを計算:式(4.45) alpha_hat_k = np.sum(s_nk, axis=0) + alpha_k # i回目のパラメータを記録 _, s_n = np.where(s_nk == 1) # クラスタ番号を抽出 trace_s_in.append(s_n.copy()) trace_beta_ik.append(beta_hat_k.copy()) trace_m_ikd.append(m_hat_kd.copy()) trace_w_ikdd.append(w_hat_kdd.copy()) trace_nu_ik.append(nu_hat_k.copy()) trace_alpha_ik.append(alpha_hat_k.copy()) # 動作確認 #print(str(i + 1) + ' (' + str(np.round((i + 1) / MaxIter * 100, 1)) + '%)')
・処理の解説
崩壊型ギブスサンプリングで行う処理を確認していきます。
・$n$番目のデータの除去
$n$番目の潜在変数$\mathbf{s}_n$を含めない超パラメータ$\mathbf{S}_{\backslash n}$を用いて、パラメータ$\hat{\mathbf{m}}_{\backslash n},\ \hat{\boldsymbol{\beta}}_{\backslash n},\ \hat{\mathbf{W}}_{\backslash n},\ \hat{\boldsymbol{\nu}}_{\backslash n},\ \hat{\boldsymbol{\alpha}}_{\backslash n}$を計算します。
n番目の潜在変数を0に置き換えます。
# n番目のデータの潜在変数を初期化 s_nk[n] = np.repeat(0, K)
s_nkのn行目を全て0にした上で、初期値のときと同じ計算すると、式(4.128)と式(4.73)
の計算となります。(本当は本にあるように$n$番目のデータに関して足し引きする処理をしたかったのですが、実装できませんでした。)
・潜在変数のサンプリング
$\mathbf{s}_n$の事後分布(のパラメータ)を計算して、新たな$\mathbf{s}_n$をサンプルします。
$\mathbf{s}_n$の事後分布(サンプリング確率)は
を正規化することで求まります。サンプリングに用いるこの2つの分布を計算しましょう。
m_hat_kd, beta_hat_kとw_hat_kkd, nu_hat_kを用いて、式(4.124)の前の項のパラメータを計算します。
# スチューデントのt分布のパラメータを計算 nu_st_hat_k = 1 - D + nu_hat_k mu_st_hat_kd = m_hat_kd.copy() term_lmd_k = nu_st_hat_k * beta_hat_k / (1 + beta_hat_k) lambda_st_hat_kdd = term_lmd_k.reshape((K, 1, 1)) * w_hat_kdd
$K$個のスチューデントのt分布のパラメータ$\hat{\boldsymbol{\nu}}^{(s)}_{\backslash n} = \{\hat{\nu}^{(s)}_{1 \backslash n}, \hat{\nu}^{(s)}_{2 \backslash n}, \cdots, \hat{\nu}^{(s)}_{K \backslash n}\}$、$\hat{\boldsymbol{\mu}}^{(s)}_{\backslash n} = \{\hat{\boldsymbol{\mu}}^{(s)}_{1 \backslash n}, \hat{\boldsymbol{\mu}}^{(s)}_{2 \backslash n}, \cdots, \hat{\boldsymbol{\mu}}^{(s)}_{K \backslash n}\}$、$\hat{\boldsymbol{\mu}}^{(s)}_{k \backslash n} = (\hat{\mu}^{(s)}_{k,1 \backslash n}, \hat{\mu}^{(s)}_{k,2 \backslash n}, \cdots, \hat{\mu}^{(s)}_{k,D \backslash n})$、$\hat{\boldsymbol{\Lambda}}^{(s)}_{\backslash n} = \{\hat{\boldsymbol{\Lambda}}^{(s)}_{1 \backslash n}, \hat{\boldsymbol{\Lambda}}^{(s)}_{2 \backslash n}, \cdots, \hat{\boldsymbol{\Lambda}}^{(s)}_{K \backslash n}\}$、$\hat{\boldsymbol{\Lambda}}^{(s)}_{k \backslash n} = (\hat{\lambda}^{(s)}_{k,1,1 \backslash n}, \cdots, \hat{\lambda}^{(s)}_{k,D,D \backslash n})$の各項を
で計算して、それぞれnu_st_hat_k, mu_st_hat_kd, lambda_st_hat_kddとします。
# 確認 print(nu_st_hat_k) print(mu_st_hat_kd) print(lambda_st_hat_kdd)
[110 69 73]
[[ 5.61440921 38.30120363]
[-18.36582238 -12.82442315]
[ 30.29981635 -17.36361393]]
[[[ 0.00402325 -0.00035665]
[-0.00035665 0.00311277]]
[[ 0.00559949 0.00256636]
[ 0.00256636 0.00539429]]
[[ 0.00539537 0.00127615]
[ 0.00127615 0.00374864]]]
nu_st_hat_k, mu_st_hat_kd, lambda_st_hat_kddを用いて、多次元スチューデントのt分布の確率密度を計算します。
# スチューデントのt分布の確率密度を計算:式(4.129) density_st_k = np.zeros(K) for k in range(K): density_st_k[k] = multivariate_t.pdf( x=x_nd[n], loc=mu_st_hat_kd[k], shape=np.linalg.inv(lambda_st_hat_kdd[k]), df=nu_st_hat_k[k] )
クラスタごとに確率密度を計算して、density_st_kに代入していきます。多次元スチューデントのt分布の確率密度は、SciPyライブラリのmultivariate_t.pdf()で計算できます。データの引数xにn番目のデータx_nd[n]、平均の引数locにmu_st_hat_kd[k]、スケール行列のshapeにlambda_st_hat_kdd[k]の逆行列、自由度の引数dfにnu_st_hat_k[k]を指定します。
# 確認 print(density_st_k)
[1.87135293e-04 1.80438926e-09 3.90827085e-08]
alpha_hat_kを用いて、式(4.124)の後の項のパラメータを計算します。
# カテゴリ分布のパラメータを計算:式(4.75)
eta_k = alpha_hat_k / np.sum(alpha_hat_k)
カテゴリ分布のパラメータ$\boldsymbol{\eta}_{\backslash n} = (\eta_{1 \backslash n}, \eta_{2 \backslash n}, \cdots, \eta_{K \backslash n})$の各項を
で計算して、eta_kとします。
# 確認 print(eta_k) print(np.sum(eta_k))
[0.43529412 0.2745098 0.29019608]
1.0
和が1になります。
式(4.124)の計算をして、結果を正規化し$\mathbf{s}_n$のサンプリング確率を求めます。
# 潜在変数のサンプリング確率を計算:式(4.124) tmp_p_k = density_st_k * eta_k prob_s_k = tmp_p_k / np.sum(tmp_p_k) # 正規化
多次元スチューデントのt分布の確率密度とカテゴリ分布のパラメータの積を計算します。(式(4.75)の計算にて正規化を行っているため、tmp_prob_kの時点で和が1になります。なので、ここでの正規化の必要性はありません。)
# 確認 print(prob_s_k) print(np.sum(prob_s_k))
[9.99854709e-01 6.07976074e-06 1.39211332e-04]
1.0
prob_s_kを用いて、新たな$\mathbf{s}_n$を生成します。
# n番目のデータの潜在変数をサンプル s_nk[n] = np.random.multinomial(n=1, pvals=prob_s_k, size=1).flatten()
prob_s_kをパラメータとして持つカテゴリ分布(4.124)から$\mathbf{s}_n$をサンプルします。
出力は、全ての要素を0にしたs_nk[n]に出力を代入します。
これで$n$番目の潜在変数$\mathbf{s}_n$を更新できました。for文で$N$個全てのデータに関して同じ処理を行い全ての潜在変数$\mathbf{S}$を更新します。
次の$n + 1$番目のデータに関する処理では、始めに$\mathbf{S}_{\backslash n+1}$を用いて超パラメータ$\hat{\mathbf{m}}_{\backslash n+1},\ \hat{\boldsymbol{\beta}}_{\backslash n+1},\ \hat{\mathbf{W}}_{\backslash n+1},\ \hat{\boldsymbol{\nu}}_{\backslash n+1},\ \hat{\boldsymbol{\alpha}}_{\backslash n+1}$を計算します。このとき、更新した$\mathbf{s}_n$を含め更新前の$\mathbf{s}_{n+1}$を除くことで超パラメータを更新しています。
最後に、初期値のときと同様にして、更新した$\mathbf{S}$を用いた超パラメータ$\hat{\mathbf{m}},\ \hat{\boldsymbol{\beta}},\ \hat{\mathbf{W}},\ \hat{\boldsymbol{\nu}},\ \hat{\boldsymbol{\alpha}}$を更新します。これは$N$個全てのデータを用いた超パラメータです。
以上が崩壊型ギブスサンプリングによる近似推論で行うの個々の処理です。これを指定した回数繰り返し行うことで、各パラメータの値を徐々に真の値に近付けていきます。
・推論結果の確認
推論結果を確認していきます。
・最後のパラメータの確認
超パラメータm_hat_k, beta_hat_kとnu_hat_k, w_hat_ddkの最後の更新値を用いて、$\boldsymbol{\mu}$の事後分布を計算します。
# muの事後分布を計算 posterior_density_kg = np.empty((K, mu_point.shape[0])) for k in range(K): # クラスタkのmuの事後分布を計算 posterior_density_kg[k] = multivariate_normal.pdf( x=mu_point, mean=m_hat_kd[k], cov=np.linalg.inv(beta_hat_k[k] * nu_hat_k[k] * w_hat_kdd[k]) )
# muの事後分布を作図 plt.figure(figsize=(12, 9)) plt.scatter(x=mu_truth_kd[:, 0], y=mu_truth_kd[:, 1], label='true val', color='red', s=100, marker='x') # 真の平均 for k in range(K): plt.contour(mu_0_grid, mu_1_grid, posterior_density_kg[k].reshape(mu_dim)) # 事後分布 plt.suptitle('Gaussian Distribution', fontsize=20) plt.title('iter:' + str(MaxIter) + ', N=' + str(N) + ', K=' + str(K), loc='left') plt.xlabel('$\mu_1$') plt.ylabel('$\mu_2$') plt.colorbar() plt.legend() plt.show()
変数名やタイトル以外は「初期値の設定」のときと同じコードです。
$K$個のガウス分布がそれぞれ真の平均付近をピークとする分布を推定できています。
m_hat_k, beta_hat_kとnu_hat_k, w_hat_ddkから求めた各パラメータの期待値$\mathbb{E}[\boldsymbol{\mu}],\ \mathbb{E}[\boldsymbol{\Lambda}],\ \mathbb{E}[\boldsymbol{\pi}]$を用いて、混合分布を計算します。
# 最後に更新したパラメータの期待値による混合分布を計算 res_density = 0 for k in range(K): # クラスタkの分布の確率密度を計算 tmp_density = multivariate_normal.pdf( x=x_point, mean=m_hat_kd[k], cov=np.linalg.inv(nu_hat_k[k] * w_hat_kdd[k]) ) # K個の分布の加重平均を計算 res_density += alpha_hat_k[k] / np.sum(alpha_hat_k) * tmp_density
# 最後に更新したパラメータの期待値による分布を作図 plt.figure(figsize=(12, 9)) plt.contour(x_1_grid, x_2_grid, true_model.reshape(x_dim), alpha=0.5, linestyles='dashed') # 真の分布 plt.scatter(x=mu_truth_kd[:, 0], y=mu_truth_kd[:, 1], color='red', s=100, marker='x') # 真の平均 plt.contourf(x_1_grid, x_2_grid, res_density.reshape(x_dim), alpha=0.5) # 期待値による分布:(塗りつぶし) for k in range(K): k_idx, = np.where(s_n == k) # クラスタkのデータのインデックス plt.scatter(x=x_nd[k_idx, 0], y=x_nd[k_idx, 1], label='cluster:' + str(k + 1)) # サンプルしたクラスタ #plt.contour(x_1_grid, x_2_grid, res_density.reshape(x_dim)) # 期待値による分布:(等高線) #plt.colorbar() # 等高線の値:(等高線用) plt.suptitle('Gaussian Mixture Model:Collapsed Gibbs Sampling', fontsize=20) plt.title('$iter:' + str(MaxIter) + ', N=' + str(N) + ', \pi=[' + ', '.join([str(pi) for pi in np.round(alpha_hat_k / np.sum(alpha_hat_k), 3)]) + ']$', loc='left') plt.xlabel('$x_1$') plt.ylabel('$x_2$') plt.legend() plt.show()
タイトル以外は「初期値の設定」のときと同じコードです。上の図はplt.contour()で、下の図はplt.contourf()で作成しました。
真の分布に近い分布になっています。
また、クラスタリングできていることが確認できます。ただし、クラスタの番号についてはランダムに決まるため真のクラスタの番号と一致しません。
・超パラメータの推移の確認
超パラメータ$\hat{\mathbf{m}},\ \hat{\boldsymbol{\beta}},\ \hat{\boldsymbol{\nu}}, \hat{\mathbf{W}},\ \hat{\boldsymbol{\alpha}}$の更新値の推移を確認します。
・コード(クリックで展開)
学習ごとの値を格納した各パラメータのリストtrace_***をnp.arange()でNumPy配列に変換して作図します。
# mの推移を作図 plt.figure(figsize=(12, 9)) for k in range(K): for d in range(D): plt.plot(np.arange(MaxIter+1), np.array(trace_m_ikd)[:, k, d], label='k=' + str(k + 1) + ', d=' + str(d + 1)) plt.xlabel('iteration') plt.ylabel('value') plt.suptitle('Collapsed Gibbs Sampling', fontsize=20) plt.title('$\hat{m}$', loc='left') plt.legend() # 凡例 plt.grid() # グリッド線 plt.show()
# betaの推移を作図 plt.figure(figsize=(12, 9)) for k in range(K): plt.plot(np.arange(MaxIter + 1), np.array(trace_beta_ik)[:, k], label='k=' + str(k + 1)) plt.xlabel('iteration') plt.ylabel('value') plt.suptitle('Collapsed Gibbs Sampling', fontsize=20) plt.title('$\hat{\\beta}$', loc='left') plt.legend() # 凡例 plt.grid() # グリッド線 plt.show()
# nuの推移を作図 plt.figure(figsize=(12, 9)) for k in range(K): plt.plot(np.arange(MaxIter + 1), np.array(trace_nu_ik)[:, k], label='k=' + str(k + 1)) plt.xlabel('iteration') plt.ylabel('value') plt.suptitle('Collapsed Gibbs Sampling', fontsize=20) plt.title('$\hat{\\nu}$', loc='left') plt.legend() # 凡例 plt.grid() # グリッド線 plt.show()
# wの推移を作図 plt.figure(figsize=(12, 9)) for k in range(K): for d1 in range(D): for d2 in range(D): plt.plot(np.arange(MaxIter + 1), np.array(trace_w_ikdd)[:, k, d1, d2], alpha=0.5, label='k=' + str(k + 1) + ', d=' + str(d1 + 1) + ', d''=' + str(d2 + 1)) plt.xlabel('iteration') plt.ylabel('value') plt.suptitle('Collapsed Gibbs Sampling', fontsize=20) plt.title('$\hat{w}$', loc='left') plt.legend() # 凡例 plt.grid() # グリッド線 plt.show()
# alphaの推移を作図 plt.figure(figsize=(12, 9)) for k in range(K): plt.plot(np.arange(MaxIter + 1), np.array(trace_alpha_ik)[:, k], label='k=' + str(k + 1)) plt.xlabel('iteration') plt.ylabel('value') plt.suptitle('Collapsed Gibbs Sampling', fontsize=20) plt.title('$\hat{\\alpha}$', loc='left') plt.legend() # 凡例 plt.grid() # グリッド線 plt.show()
・おまけ：アニメーションによる推移の確認
animationモジュールを利用して、各分布の推移のアニメーション(gif画像)を作成するためのコードです。
# 追加ライブラリ import matplotlib.animation as animation
・平均パラメータの事後分布の推移
・コード(クリックで展開)
# 画像サイズを指定 fig = plt.figure(figsize=(9, 9)) # 作図処理を関数として定義 def update_posterior(i): # i回目のmuの事後分布を計算 posterior_density_kg = np.empty((K, mu_point.shape[0])) for k in range(K): # クラスタkのmuの事後分布を計算 posterior_density_kg[k] = multivariate_normal.pdf( x=mu_point, mean=trace_m_ikd[i][k], cov=np.linalg.inv(trace_beta_ik[i][k] * trace_nu_ik[i][k] * trace_w_ikdd[i][k]) ) # 前フレームのグラフを初期化 plt.cla() # i回目のmuの事後分布を作図 plt.scatter(x=mu_truth_kd[:, 0], y=mu_truth_kd[:, 1], label='true val', color='red', s=100, marker='x') # 真の平均 for k in range(K): plt.contour(mu_0_grid, mu_1_grid, posterior_density_kg[k].reshape(mu_dim)) # 事後分布 plt.suptitle('Gaussian Distribution:Collapsed Gibbs Sampling', fontsize=20) plt.title('iter:' + str(i) + ', N=' + str(N) + ', K=' + str(K), loc='left') plt.xlabel('$\mu_1$') plt.ylabel('$\mu_2$') plt.legend() # gif画像を作成 posterior_anime = animation.FuncAnimation(fig, update_posterior, frames=MaxIter + 1, interval=100) posterior_anime.save("ch4_4_4_Posterior.gif")
・パラメータの期待値による分布とクラスタの推移
・コード(クリックで展開)
# 画像サイズを指定 fig = plt.figure(figsize=(9, 9)) # 作図処理を関数として定義 def update_model(i): # i回目の混合比率を計算 pi_hat_k = trace_alpha_ik[i] / np.sum(trace_alpha_ik[i]) # i回目のサンプルによる混合分布を計算 res_density = 0 for k in range(K): # クラスタkの分布の確率密度を計算 tmp_density = multivariate_normal.pdf( x=x_point, mean=trace_m_ikd[i][k], cov=np.linalg.inv(trace_nu_ik[i][k] * trace_w_ikdd[i][k]) ) # K個の分布の加重平均を計算 res_density += pi_hat_k[k] * tmp_density # 前フレームのグラフを初期化 plt.cla() # i回目のサンプルによる分布を作図 plt.contour(x_1_grid, x_2_grid, true_model.reshape(x_dim), alpha=0.5, linestyles='dashed') # 真の分布 plt.scatter(x=mu_truth_kd[:, 0], y=mu_truth_kd[:, 1], color='red', s=100, marker='x') # 真の平均 plt.contour(x_1_grid, x_2_grid, res_density.reshape(x_dim)) # 期待値による分布:(等高線) #plt.contourf(x_1_grid, x_2_grid, res_density.reshape(x_dim), alpha=0.5) # 期待値による分布:(塗りつぶし) for k in range(K): k_idx, = np.where(trace_s_in[i] == k) plt.scatter(x=x_nd[k_idx, 0], y=x_nd[k_idx, 1], label='cluster:' + str(k + 1)) # サンプルしたクラスタ plt.suptitle('Gaussian Mixture Model:Collapsed Gibbs Sampling', fontsize=20) plt.title('$iter:' + str(i) + ', N=' + str(N) + ', \pi=[' + ', '.join([str(pi) for pi in np.round(pi_hat_k, 3)]) + ']$', loc='left') plt.xlabel('$x_1$') plt.ylabel('$x_2$') plt.legend() # gif画像を作成 model_anime = animation.FuncAnimation(fig, update_model, frames=MaxIter + 1, interval=100) model_anime.save("ch4_4_4_Model.gif")
(よく理解していないので、animationの解説は省略…とりあえずこれで動きます……)
上の図はplt.contour()で、下の図はplt.contourf()で作成しました。
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
4章完了！そして5章は読むだけにしたので、この本完了です。ずっとこの本だけをやってたわけではないですが、1年以上かかりました。
RとPythonの両方で実装してみて、ggplot2とPyPlotで一方では簡単にできることが他方でははてこずったり、データフレーム形式にまとめるのが楽だったり手間だったり、色々似てるようで似てない点の雰囲気を味わえました。2つの言語で同じものを組むのもいい経験でいした。
しかし、2言語分の記事を書くのは初めてだったのですが、予想外に大変でした。書く際はコピペするので純粋に2つ記事を書くよりは楽ではあるのですが、誤字脱字修正しようとなると労力が2倍になるのが辛かった。さらにブログ記事と手元資料でもコピペしているので、ほぼ同じテキストが4つあって、、、さらにさらに、この本のいいところがどの節でも同じ手順で進めることなので、このブログでも構成や表現を統一することを意識していて、複数の記事を連動して修正しないといけなかったり。今後はヘンに拘りを持つのは止めようと思います。でも書き切れて良かった！
最後まで読んでいただきありがとうございました。それなりに苦労したので、誰かの役に立てば幸いです。
【次節の内容】
おわり。
