はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。基本的な内容は「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」になります。「数式」と「プログラム」から理解するのが目標です。
この記事は、4.4.2項の内容です。「観測モデルを多次元ガウス混合分布(多変量正規混合分布)」、「事前分布をガウス・ウィシャート分布とディリクレ分布」とするガウス混合モデルに対するギブスサンプリングをPythonで実装します。
省略してある内容等ありますので、本とあわせて読んでください。初学者な自分が理解できるレベルまで落として書き下していますので、分かる人にはかなりくどくなっています。同じような立場の人のお役に立てれば幸いです。
【数式読解編】
【他の節一覧】
【この節の内容】
・Pythonでやってみる
人工的に生成したデータを用いて、ベイズ推論を行ってみましょう。アルゴリズム4.5を参考に実装します。
利用するライブラリを読み込みます。
# 4.4.2項で利用するライブラリ import numpy as np from scipy.stats import multivariate_normal, wishart, dirichlet # 多次元ガウス分布, ウィシャート分布, ディリクレ分布 import matplotlib.pyplot as plt
scipyライブラリのstatsモジュールから、多次元ガウス分布multivariate_normal、ウィシャート分布wishart、ディリクレ分布dirichletのクラスを利用します。確率密度の計算にはpdf()、乱数の生成にはrvs()を使います。
・観測モデルの設定
まずは、観測モデルを設定します。この節では、$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 0 1]
[1 0 0]
[0 0 1]]
これが本来は観測できない潜在変数であり、真のクラスタです。各行がデータを表し、値が1の列番号が割り当てられたクラスタ番号を表します。
s_truth_nkから各データのクラスタ番号を抽出します。
# クラスタ番号を抽出 _, s_truth_n = np.where(s_truth_nk == 1) # 確認 print(s_truth_n[:5])
[0 0 2 0 2]
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])
[[ 1.77629859e+01 4.72882578e+01]
[ 1.94567820e-02 1.76203249e+01]
[ 1.67164723e+01 -2.15832685e+01]
[ 1.65851789e+01 4.84626992e+01]
[ 2.78818169e+01 -1.21215116e+01]]
観測データの散布図を観測モデルに重ねて作図します。
# 観測データの散布図を作成 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項「多次元ガウス分布の学習と予測：精度が未知の場合」を参照してください。
・初期値の設定
設定した各パラメータの事前分布$p(\boldsymbol{\mu} | \boldsymbol{\Lambda}),\ p(\boldsymbol{\Lambda}),\ p(\boldsymbol{\pi})$に従いパラメータ$\boldsymbol{\mu},\ \boldsymbol{\Lambda},\ \boldsymbol{\pi}$を生成して初期値とします。
$\boldsymbol{\mu},\ \boldsymbol{\Lambda}$の事前分布から観測モデルのパラメータ$\boldsymbol{\mu},\ \boldsymbol{\Lambda}$を生成します。
# 観測モデルのパラメータをサンプル mu_kd = np.empty((K, D)) lambda_kdd = np.empty((K, D, D)) for k in range(K): # クラスタkの精度行列をサンプル lambda_kdd[k] = wishart.rvs(df=nu, scale=w_dd, size=1) # クラスタkの平均をサンプル mu_kd[k] = np.random.multivariate_normal( mean=m_d, cov=np.linalg.inv(beta * lambda_kdd[k]) ).flatten()
精度行列パラメータ$\boldsymbol{\Lambda}$の事前分布(ウィシャート分布)と平均パラメータ$\boldsymbol{\mu}$の事前分布(多次元ガウス分布)からそれぞれ$K$個のパラメータをサンプルして初期値とします。
ウィシャート分布に従う乱数は、SicPyライブラリのwishart.rvs()で生成できます。クラスタごとに生成するので、データ数の引数sizeに1を指定します。また、自由度の引数dfにnu、スケール行列の引数scaleにw_ddを指定します。
多次元ガウス分布の乱数生成関数np.random.multivariate_normal()の引数meanにm_d、covには今サンプルしたlambda_kdd[k]と係数betaの積をnp.linalg.inv()で逆行列に変換して指定します。2次元配列で出力されるため、flatten()で1次元配列に変換しておきます。
生成したパラメータを確認します。
# 確認 print(mu_kd) print(lambda_kdd)
[[ 41.29502182 23.02093784]
[-74.11919232 63.75529812]
[-35.84343453 25.75284128]]
[[[ 1.49649481e-04 -2.24077259e-04]
[-2.24077259e-04 1.41126515e-03]]
[[ 5.19813559e-04 -1.13218227e-04]
[-1.13218227e-04 6.80104601e-05]]
[[ 2.47702308e-04 -4.98632073e-04]
[-4.98632073e-04 2.22053690e-03]]]
$\boldsymbol{\pi}$の事前分布から$\boldsymbol{\pi}$を生成します。
# 混合比率をサンプル pi_k = dirichlet.rvs(alpha=alpha_k, size=1).flatten()
$\boldsymbol{\pi}$の事前分布(ディリクレ分布)に従い$\boldsymbol{\pi}$をサンプルして初期値とします。
ディリクレ分布に従う乱数は、SciPyライブラリの
のdirichlet.rvs()で生成できます。パラメータの引数alphaにalpha_kを指定します。2次元配列で出力されるため、flatten()で1次元配列に変換しておきます。
生成したパラメータを確認します。
# 確認 print(pi_k) print(np.sum(pi_k))
[0.12841012 0.04916972 0.82242016]
1.0
$K$個の要素の和が1になるのを確認できます。
パラメータの初期値が得られたので、グラフで確認しましょう。初期値を使った分布を計算します。
# 初期値による混合分布を計算 init_density = 0 for k in range(K): # クラスタkの分布の確率密度を計算 tmp_density = multivariate_normal.pdf( x=x_point, mean=mu_kd[k], cov=np.linalg.inv(beta * lambda_kdd[k]) ) # K個の分布の加重平均を計算 init_density += pi_k[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.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()
観測モデルのときと同様に作図します。サンプルを繰り返すことで真の分布に近付けていきます。
以上でモデルが整いました。次からは、推論を行っていきます。
・ギブスサンプリング
ギブスサンプリングでは、潜在変数$\mathbf{S}$のサンプリングと、パラメータ$\boldsymbol{\mu},\ \boldsymbol{\Lambda},\ \boldsymbol{\pi}$のサンプリングを交互に行います。
まずは、パラメータの初期値を用いて$\mathbf{S}$の条件付き分布(事後分布)$p(\mathbf{S} | \mathbf{X}, \boldsymbol{\mu}, \boldsymbol{\Lambda}, \boldsymbol{\pi})$を計算します。得られた条件付き分布から$\mathbf{S}$をサンプルします。次に、サンプルした$\mathbf{S}$を用いて各パラメータ$\boldsymbol{\mu},\ \boldsymbol{\Lambda},\ \boldsymbol{\pi}$の条件付き分布(事後分布)$p(\boldsymbol{\mu} | \boldsymbol{\Lambda}, \mathbf{X}, \mathbf{S}),\ p(\boldsymbol{\Lambda} | \mathbf{X}, \mathbf{S}),\ p(\boldsymbol{\pi} | \mathbf{X}, \mathbf{S})$を計算します。得られた条件付き分布からパラメータをサンプルします。ここまでが1回の試行です。
さらに次の試行では、前回サンプルしたパラメータを用いて$\mathbf{S}$の条件付き分布を更新となります。この処理を指定した回数繰り返します。
・コード全体
各パラメータの変数に関して、添字を使って代入するため予め変数を用意しておく必要があります。
また、後で各パラメータの更新値や分布の推移を確認するために、それぞれ試行ごとにリスト型のtrace_***に値(NumPy配列)を保存していきます。関心のあるパラメータを記録してください。
# 試行回数を指定 MaxIter = 150 # パラメータを初期化 eta_nk = np.zeros((N, K)) s_nk = np.zeros((N, K)) beta_hat_k = np.zeros(K) m_hat_kd = np.zeros((K, D)) w_hat_kdd = np.zeros((K, D, D)) nu_hat_k = np.zeros(K) alpha_hat_k = np.zeros(K) # 推移の確認用の受け皿を作成 trace_s_in = [np.repeat(np.nan, N)] trace_mu_ikd = [mu_kd.copy()] trace_lambda_ikdd = [lambda_kdd.copy()] trace_pi_ik = [pi_k.copy()] trace_beta_ik = [np.repeat(beta, K)] trace_m_ikd = [np.repeat(m_d.reshape((1, D)), K, axis=0)] trace_w_ikdd = [np.repeat(w_dd.reshape((1, D, D)), K, axis=0)] trace_nu_ik = [np.repeat(nu, K)] trace_alpha_ik = [alpha_k.copy()] # ギブスサンプリング for i in range(MaxIter): # 潜在変数の事後分布のパラメータを計算:式(4.94) for k in range(K): tmp_eta_n = np.diag( -0.5 * (x_nd - mu_kd[k]).dot(lambda_kdd[k]).dot((x_nd - mu_kd[k]).T) ).copy() # (何故か書き替え禁止になるのを防ぐためのcopy()) tmp_eta_n += 0.5 * np.log(np.linalg.det(lambda_kdd[k]) + 1e-7) tmp_eta_n += np.log(pi_k[k] + 1e-7) eta_nk[:, k] = np.exp(tmp_eta_n) eta_nk /= np.sum(eta_nk, axis=1, keepdims=True) # 正規化 # 潜在変数をサンプル：式(4.93) for n in range(N): s_nk[n] = np.random.multinomial(n=1, pvals=eta_nk[n], size=1).flatten() # 観測モデルのパラメータをサンプル for k in range(K): # muの事後分布のパラメータを計算：式(4.99) beta_hat_k[k] = np.sum(s_nk[:, k]) + beta m_hat_kd[k] = np.sum(s_nk[:, k] * x_nd.T, axis=1) m_hat_kd[k] += beta * m_d m_hat_kd[k] /= beta_hat_k[k] # lambdaの事後分布のパラメータを計算：式(4.103) tmp_w_dd = np.dot((s_nk[:, k] * x_nd.T), x_nd) tmp_w_dd += beta * np.dot(m_d.reshape(D, 1), m_d.reshape(1, D)) tmp_w_dd -= beta_hat_k[k] * np.dot(m_hat_kd[k].reshape(D, 1), m_hat_kd[k].reshape(1, D)) tmp_w_dd += np.linalg.inv(w_dd) w_hat_kdd[k] = np.linalg.inv(tmp_w_dd) nu_hat_k[k] = np.sum(s_nk[:, k]) + nu # lambdaをサンプル：式(4.102) lambda_kdd[k] = wishart.rvs(size=1, df=nu_hat_k[k], scale=w_hat_kdd[k]) # muをサンプル：式(4.98) mu_kd[k] = np.random.multivariate_normal( mean=m_hat_kd[k], cov=np.linalg.inv(beta_hat_k[k] * lambda_kdd[k]), size=1 ).flatten() # 混合比率のパラメータを計算：式(4.45) alpha_hat_k = np.sum(s_nk, axis=0) + alpha_k # piをサンプル：式(4.44) pi_k = dirichlet.rvs(size=1, alpha=alpha_hat_k).flatten() # 値を記録 _, s_n = np.where(s_nk == 1) trace_s_in.append(s_n.copy()) trace_mu_ikd.append(mu_kd.copy()) trace_lambda_ikdd.append(lambda_kdd.copy()) trace_pi_ik.append(pi_k.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)) + '%)')
・処理の解説
ギブスサンプリングで行う処理を確認していきます。
・潜在変数のサンプリング
潜在変数$\mathbf{S}$の条件付き分布(のパラメータ)を計算して、新たな$\mathbf{S}$をサンプルします。
$\mathbf{S}$の条件付き分布のパラメータを計算します。
# 潜在変数の事後分布のパラメータを計算:式(4.94) for k in range(K): tmp_eta_n = np.diag( -0.5 * (x_nd - mu_kd[k]).dot(lambda_kdd[k]).dot((x_nd - mu_kd[k]).T) ).copy() # (何故か書き替え禁止になるのを防ぐためのcopy()) tmp_eta_n += 0.5 * np.log(np.linalg.det(lambda_kdd[k]) + 1e-7) tmp_eta_n += np.log(pi_k[k] + 1e-7) eta_nk[:, k] = np.exp(tmp_eta_n) eta_nk /= np.sum(eta_nk, axis=1, keepdims=True) # 正規化
$n$番目のデータの潜在変数$\mathbf{s}_n$の条件付き分布(カテゴリ分布)のパラメータ$\boldsymbol{\eta}_n = (\eta_{n,1}, \eta_{n,2}, \cdots, \eta_{n,K})$の各項は、次の式で計算します。
$\eta_{n,k}$は、$n$番目の観測データ$\mathbf{x}_n$にクラスタ$k$が割り当てられる($s_{n,k} = 1$となる)確率を表します。
ただし、$N$個全てのデータを一度で処理するために、$(\mathbf{x}_n - \boldsymbol{\mu}_k)^{\top} \boldsymbol{\Lambda}_k (\mathbf{x}_n - \boldsymbol{\mu}_k)$の計算を1から$N$の全ての組み合わせで行います(tmp_eta_nの計算)。この計算結果は、$N \times N$のマトリクスになります。これは例えば、1行$N$列目の要素は$(\mathbf{x}_1 - \boldsymbol{\mu}_k)^{\top} \boldsymbol{\Lambda}_k (\mathbf{x}_N - \boldsymbol{\mu}_k)$の計算結果です。これは不要なので、np.diag()で対角成分(同じデータによる計算結果)のみを取り出します。
また、対数の計算において$\ln 0$とならないように、微小な値1e-7($10^{-7} = 0.0000001$)を加えています。
最後に、データごとの和が1となるように
で正規化します。
# 確認 print(np.round(eta_nk[:5], 3)) print(np.sum(eta_nk[:5], axis=1))
[[0. 0. 1. ]
[0.023 0.005 0.972]
[0.994 0.001 0.005]
[0. 0. 1. ]
[0.993 0. 0.007]]
[1. 1. 1. 1. 1.]
行(データ)ごとの和が1になるのを確認できます。
eta_nkを用いて、新たな$\mathbf{S}$を生成します。
# 潜在変数をサンプル：式(4.93) for n in range(N): s_nk[n] = np.random.multinomial(n=1, pvals=eta_nk[n], size=1).flatten()
真のクラスタのときと同様に、潜在変数をサンプルします。
ただし、ここではデータによってパラメータ(各クラスタの割り当て確率)が異なるため、for文で1データずつ生成します。
# 確認 print(s_nk[:5])
[[0. 0. 1.]
[0. 0. 1.]
[1. 0. 0.]
[0. 0. 1.]
[1. 0. 0.]]
これで$\mathbf{S}$が得られました。次は、更新した$\mathbf{S}$を用いて、パラメータ$\boldsymbol{\mu},\ \boldsymbol{\Lambda},\ \boldsymbol{\pi}$を更新します。
・観測モデルのパラメータのサンプリング
観測モデルのパラメータ$\boldsymbol{\mu},\ \boldsymbol{\Lambda}$の条件付き分布(のパラメータ)を計算し、新たな$\boldsymbol{\mu},\ \boldsymbol{\Lambda}$をサンプルします。
$\boldsymbol{\mu}_k$の条件付き分布のパラメータを計算します。
# muの事後分布のパラメータを計算：式(4.99) beta_hat_k[k] = np.sum(s_nk[:, k]) + beta m_hat_kd[k] = np.sum(s_nk[:, k] * x_nd.T, axis=1) m_hat_kd[k] += beta * m_d m_hat_kd[k] /= beta_hat_k[k]
$\boldsymbol{\mu}_k$の条件付き分布(多次元ガウス分布)の平均パラメータ$\hat{\mathbf{m}}_k = (\hat{m}_{k,1}, \hat{m}_{k,2}, \cdots, \hat{m}_{k,D})$、精度行列パラメータの係数$\hat{\beta}_k$を
で計算して、それぞれbeta_hat_k、m_hat_kdに代入します。
# 確認 print(beta_hat_k) print(m_hat_kd)
[ 79. 55. 119.]
[[ 28.02932663 -20.16180223]
[-19.57323076 -9.457385 ]
[ 6.28368833 33.89046545]]
ちなみに、$\sum_{n=1}^N s_{n,k} \mathbf{x}_n$の計算に注目すると
k = 0 # クラスタ番号 print(np.round(x_nd[:5], 1)) print(s_nk[:5, k])
[[ 17.8 47.3]
[ 0. 17.6]
[ 16.7 -21.6]
[ 16.6 48.5]
[ 27.9 -12.1]]
[0. 0. 1. 0. 1.]
print(np.round(s_nk[:5, k] * x_nd[:5].T, 1).T) print(np.sum(s_nk[:5, k] * x_nd[:5].T, axis=1))
[[ 0. 0. ]
[ 0. 0. ]
[ 16.7 -21.6]
[ 0. 0. ]
[ 27.9 -12.1]]
[ 44.59828916 -33.70478007]
対象としているクラスタのデータには1を掛けてそのままの値で、それ以外のクラスタのデータには0を掛けて値が0になっています。つまり、対象となるクラスタのデータのみの和をとっています。
同様に、$\sum_{n=1}^N s_{n,k}$は各クラスタが割り当てられたデータ数を表します。
$\boldsymbol{\lambda}_k$の条件付き分布のパラメータを計算します。
# lambdaの事後分布のパラメータを計算：式(4.103) tmp_w_dd = np.dot((s_nk[:, k] * x_nd.T), x_nd) tmp_w_dd += beta * np.dot(m_d.reshape(D, 1), m_d.reshape(1, D)) tmp_w_dd -= beta_hat_k[k] * np.dot(m_hat_kd[k].reshape(D, 1), m_hat_kd[k].reshape(1, D)) tmp_w_dd += np.linalg.inv(w_dd) w_hat_kdd[k] = np.linalg.inv(tmp_w_dd) nu_hat_k[k] = np.sum(s_nk[:, k]) + nu
$\boldsymbol{\lambda}_k$の条件付き分布(ウィシャート分布)のパラメータ$\hat{\mathbf{W}}_k = (\hat{w}_{k,1,1}, \cdots, \hat{w}_{k,D,D})$、自由度$\hat{\nu}_k$を
で計算して、結果をそれぞれnu_hat_k、w_hat_ddkに代入します。
# 確認 print(w_hat_kdd) print(nu_hat_k)
[[[ 4.55924169e-05 6.45009451e-06]
[ 6.45009451e-06 5.25854742e-05]]
[[ 1.59406207e-04 1.97749953e-05]
[ 1.97749953e-05 8.92096060e-05]]
[[ 2.85394262e-05 -5.12987474e-06]
[-5.12987474e-06 3.04729579e-05]]]
[ 80. 56. 120.]
nu_hat_k, w_hat_ddkを用いて、新たな$\boldsymbol{\lambda}_k$を生成します。
# lambdaをサンプル：式(4.102) lambda_kdd[k] = wishart.rvs(size=1, df=nu_hat_k[k], scale=w_hat_kdd[k])
パラメータ$\hat{\mathbf{W}}_k$を持つ自由度$\hat{\nu}_k$のウィシャート分布から$\boldsymbol{\lambda}_k$をサンプルします。
初期値のときと同様に、wishart.rvs()の自由度の引数dfにnu_hat_k[k]、スケール行列の引数scaleにw_hat_kdd[k]を指定します。
# 確認 print(lambda_kdd[k])
[[0.00354549 0.00078969]
[0.00078969 0.00452033]]
m_hat_kd, beta_hat_kとlambda_ddkを用いて、新たな$\boldsymbol{\mu}_k$を生成します。
# muをサンプル：式(4.98) mu_kd[k] = np.random.multivariate_normal( mean=m_hat_kd[k], cov=np.linalg.inv(beta_hat_k[k] * lambda_kdd[k]), size=1 ).flatten()
平均$\hat{\mathbf{m}}_k$、精度行列$\hat{\beta}_k \boldsymbol{\lambda}_k$の多次元ガウス分布から$\boldsymbol{\mu}_k$をサンプルします。
np.random.multivariate_normal()の平均の引数meanにm_hat_kd[k]、分散共分散行列の引数covに今サンプルしたlambda_kdd[k]と係数beta_hat_k[k]との積をnp.linalg.inv()で逆行列に変換して指定します。
# 確認 print(mu_kd[k])
[ 28.57478108 -22.40480772]
ここまでの処理をfor文ループで$K$回行うことで、$K$個のクラスタの超パラメータ$\hat{\mathbf{m}} = \{\hat{\mathbf{m}}_1, \hat{\mathbf{m}}_2, \cdots, \hat{\mathbf{m}}_K\}$、$\hat{\boldsymbol{\beta}} = \{\hat{\beta}_1, \hat{\beta}_2, \cdots, \hat{\beta}_k\}$、$\hat{\boldsymbol{\nu}} = \{\hat{\nu}_1, \hat{\nu}_2, \cdots, \hat{\nu}_K\}$、$\hat{\mathbf{W}} = \{\hat{\mathbf{W}}_1, \hat{\mathbf{W}}_2, \cdots, \hat{\mathbf{W}}_K\}$が得られます。
・混合比率のサンプリング
混合比率パラメータ$\boldsymbol{\pi}$の条件付き分布(のパラメータ)を計算し、新たな$\boldsymbol{\pi}$をサンプルします。
$\boldsymbol{\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. 56. 120.]
alpha_hat_kを用いて、新たな$\boldsymbol{\pi}$を生成します。
# piをサンプル：式(4.44) pi_k = dirichlet.rvs(size=1, alpha=alpha_hat_k).flatten()
パラメータ$\hat{\boldsymbol{\alpha}}$を持つディリクレ分布から$\boldsymbol{\pi}$をサンプルします。
事前分布からのサンプルのときと同様に、dirichlet.rvs()の引数alphaにalpha_hat_kを指定します。
# 確認 print(pi_k) print(np.sum(pi_k))
[0.28873963 0.23878639 0.47247398]
1.0
以上がギブスサンプリングによる近似推論で行うの個々の処理です。これを指定した回数繰り返し行うことで、各パラメータの値を徐々に真の値に近付けていきます。
・推論結果の確認
推論結果を確認していきます。
・最後のパラメータの確認
超パラメータm_hat_k, beta_hat_kとlambda_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] * lambda_ddk[k]) )
事前分布のときと同様に計算します。ただし、超パラメータが$K$個になったので、多次元ガウス分布も$K$個になります。$K$回分の計算結果を2次元配列に格納していきます。
クラスタ$k$の事後分布の平均は$\hat{\mathbf{m}}_k$、精度行列は$\hat{\beta}_k \boldsymbol{\Lambda}_k$です。
作成した配列は次のようになります。
# 確認 print(posterior_density_kg[:, :5])
[[0.00000000e+000 0.00000000e+000 0.00000000e+000 0.00000000e+000
0.00000000e+000]
[2.26628395e-207 6.47414410e-205 1.73031714e-202 4.32657629e-200
1.01213474e-197]
[0.00000000e+000 0.00000000e+000 0.00000000e+000 0.00000000e+000
0.00000000e+000]]
$\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(MaxIter) + ', N=' + str(N) + ', K=' + str(K), loc='left') plt.xlabel('$\mu_1$') plt.ylabel('$\mu_2$') plt.colorbar() plt.legend() plt.show()
事前分布のときと同様に作図できます。
$K$個のガウス分布がそれぞれ真の平均付近をピークとする分布を推定できています。
最後の試行でサンプルしたパラメータmu_kd, lambda_ddkとpi_kを用いて、混合分布を計算します。
# 最後にサンプルしたパラメータによる混合分布を計算 res_density = 0 for k in range(K): # クラスタkの分布の確率密度を計算 tmp_density = multivariate_normal.pdf( x=x_point, mean=mu_kd[k], cov=np.linalg.inv(lambda_kdd[k]) ) # K個の分布の加重平均を計算 res_density += tmp_density * pi_k[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.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) 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:Gibbs Sampling', fontsize=20) plt.title('$iter:' + str(MaxIter) + ', N=' + str(N) + ', \pi=[' + ', '.join([str(pi) for pi in np.round(pi_k, 3)]) + ']$', loc='left') plt.xlabel('$x_1$') plt.ylabel('$x_2$') 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('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('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('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('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('Gibbs Sampling', fontsize=20) plt.title('$\hat{\\alpha}$', loc='left') plt.legend() # 凡例 plt.grid() # グリッド線 plt.show()
・パラメータのサンプルの推移の確認
サンプルされた$\boldsymbol{\mu},\ \boldsymbol{\Lambda},\ \boldsymbol{\pi}$の値の推移を確認しましょう。
・コード(クリックで展開)
同様に作図します。それぞれ真の値***_truthを水平線plt.hlines()で描画しておきます。
# muの推移を作図 plt.figure(figsize=(12, 9)) plt.hlines(y=mu_truth_kd, xmin=0, xmax=MaxIter + 1, label='true val', color='red', linestyles='--') # 真の値 for k in range(K): for d in range(D): plt.plot(np.arange(MaxIter+1), np.array(trace_mu_ikd)[:, k, d], label='k=' + str(k + 1) + ', d=' + str(d + 1)) plt.xlabel('iteration') plt.ylabel('value') plt.suptitle('Gibbs Sampling', fontsize=20) plt.title('$\mu$', loc='left') plt.legend() # 凡例 plt.grid() # グリッド線 plt.show()
# lambdaの推移を作図 plt.figure(figsize=(12, 9)) plt.hlines(y=np.linalg.inv(sigma2_truth_kdd), xmin=0, xmax=MaxIter + 1, label='true val', color='red', linestyles='--') # 真の値 for k in range(K): for d1 in range(D): for d2 in range(D): plt.plot(np.arange(MaxIter + 1), np.array(trace_lambda_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('Gibbs Sampling', fontsize=20) plt.title('$\Lambda$', loc='left') plt.legend() # 凡例 plt.grid() # グリッド線 plt.show()
# piの推移を作図 plt.figure(figsize=(12, 9)) plt.hlines(y=pi_truth_k, xmin=0, xmax=MaxIter + 1, label='true val', color='red', linestyles='--') # 真の値 for k in range(K): plt.plot(np.arange(MaxIter + 1), np.array(trace_pi_ik)[:, k], label='k=' + str(k + 1)) plt.xlabel('iteration') plt.ylabel('value') plt.suptitle('Gibbs Sampling', fontsize=20) plt.title('$\pi$', loc='left') plt.legend() # 凡例 plt.grid() # グリッド線 plt.show()
試行回数が増えるに従って真の値に近づいているのを確認できます。
・おまけ：アニメーションによる推移の確認
animationモジュールを利用して、各分布の推移のアニメーション(gif画像)を作成するためのコードです。
・コード(クリックで展開)
# 追加ライブラリ import matplotlib.animation as animation
・パラメータの事後分布の推移
# 画像サイズを指定 fig = plt.figure(figsize=(12, 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_lambda_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 Mixture Model: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_2_Posterior.gif")
・サンプルしたパラメータによる分布とクラスタの推移
・コード(クリックで展開)
# 画像サイズを指定 fig = plt.figure(figsize=(12, 9)) # 作図処理を関数として定義 def update_model(i): # i回目のサンプルによる混合分布を計算 res_density = 0 for k in range(K): # クラスタkの分布の確率密度を計算 tmp_density = multivariate_normal.pdf( x=x_point, mean=trace_mu_ikd[i][k], cov=np.linalg.inv(trace_lambda_ikdd[i][k]) ) # K個の分布の加重平均を計算 res_density += trace_pi_ik[i][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) # クラスタ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:Gibbs Sampling', fontsize=20) plt.title('$iter:' + str(i) + ', N=' + str(N) + ', \pi=[' + ', '.join([str(pi) for pi in np.round(trace_pi_ik[i], 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_2_Model.gif")
(よく理解していないので、animationの解説は省略…とりあえずこれで動きます……)
上はplt.contour()、下はplt.contourf()で作成したものです。(アップロードサイズの上限を超えたためアニメーションを分割しました。左が0から100回目の結果まで、右が101から200回目の結果です。)
本当は、廃棄期間を設定した上でサンプル平均をとる操作が必要なはずですがここでは扱いません。
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
今回は先にPythonで組んでみました。そしてRではどこかミスっててまだうまく推定できてません。
この記事の投稿日の朝に公開されたモーニング娘。'20の新曲「ギューされたいだけなのに」です！
ウサギちゃんシンドローム♪
さらに前日に公開されたアンジュルムの新曲「SHAKA SHAKA TO LOVE」です！
歯磨けよ～
【次節の内容】
- 2021/05/14：加筆修正しました。