はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。基本的な内容は「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」になります。「数式」と「プログラム」から理解するのが目標です。
この記事は4.4.2項の内容です。ガウス混合モデルにおけるギブスサンプリングによる近似推論をPythonで実装します。
省略してある内容等ありますので、本と併せて読んでください。初学者な自分が理解できるレベルまで落として書き下していますので、分かる人にはかなりくどくなっています。同じような立場の人のお役に立てれば幸いです。
【数式読解編】
【前節の内容】
まだ
【他の節一覧】
【この節の内容】
・Pythonでやってみる
ガウス混合モデルに従い生成したデータを用いて、パラメータを推定してみましょう。
利用するライブラリを読み込みます。
# 利用ライブラリ 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()を使います。
・真の観測モデルの設定
観測モデルのパラメータを設定します。作図に関しては、2次元のグラフで表現するため$D = 2$のときのみ動作します。
# 真の観測モデルのパラメータを指定 D = 2 K = 3 mu_true_kd = np.array( [[0.0, 4.0], [-5.0, -5.0], [5.0, -2.5]] ) sigma2_true_kdd = np.array( [[[8.0, 0.0], [0.0, 8.0]], [[4.0, -2.5], [-2.5, 4.0]], [[6.5, 4.0], [4.0, 6.5]]] ) # 確認 print(mu_true_kd) print(sigma2_true_kdd)
[[ 0. 4. ]
[-5. -5. ]
[ 5. -2.5]]
[[[ 8. 0. ]
[ 0. 8. ]]
[[ 4. -2.5]
[-2.5 4. ]]
[[ 6.5 4. ]
[ 4. 6.5]]]
真の観測モデルにおけるクラスタごとの平均パラメータ$\boldsymbol{\mu} = \{\boldsymbol{\mu_1}, \cdots, \boldsymbol{\mu}_K\}$、$\boldsymbol{\mu}_k = \{\mu_{k1}, \cdots, \mu_{kD}\}$をmu_true_kdとし、分散共分散行列(精度行列の逆行列)パラメータ$\boldsymbol{\Lambda}^{-1} = \{\boldsymbol{\Lambda}_1^{-1}, \cdots, \boldsymbol{\Lambda}_K^{-1}\}$、$\boldsymbol{\Lambda}_k^{-1} = \boldsymbol{\Sigma}_k = \{\sigma_{k11}^2, \cdots, \sigma_{kDD}^2\}$をsigma2_true_kddとして、値を指定します。ただし$\boldsymbol{\Sigma}_k$は、正定値行列である必要があります。
この2つのパラメータの値を観測データから求めることがこの例の目的となります。
続いて混合比率を設定します。
# 真の混合比率を指定 pi_true_k = np.array([0.5, 0.2, 0.3])
混合比率$\boldsymbol{\pi} = \{\pi_1, \cdots, \pi_K\}$をpi_true_kとして、値を指定します。ここで$\pi_k$は各データ$\mathbf{x}_n$のクラスタ(各潜在変数$\mathbf{s}_n$)が$k$となる確率であり、$0 \leq \pi_k \leq 1,\ \sum_{k=1}^K \pi_k = 1$の値をとります。このパラメータの値も求めます。
以上が観測モデルの設定です。続いて、設定したガウス混合モデルに従ってデータを生成します。まずは各データが持つクラスタを割り当てます。
# 各データの真のクラスタを生成 N = 250 s_true_nk = np.random.multinomial(n=1, pvals=pi_true_k, size=N) # 確認 print(s_true_nk[0:5])
[[1 0 0]
[1 0 0]
[1 0 0]
[0 0 1]
[0 1 0]]
$\boldsymbol{\pi}$をパラメータとするカテゴリ分布に従い、潜在変数$\mathbf{S} = \{s_{11}, \cdots, s_{NK}\}$を生成します。カテゴリ分布の乱数は、多項分布の乱数生成関数np.random.multinomial()の引数nに1を指定することで生成できます。確率値引数pvalsにpi_true_k、サンプルサイズ引数sizeにNを指定します。
出力は、各列がクラスタ1から$K$に対応していて、各行(各データ)に割り当てられたクラスタが1でそれ以外は0となります。これが(本来は観測できない)真のクラスタとなります。
処理の都合上、割り当てられたクラスタを抽出しておきます。
# 各データのクラスタを抽出 _, s_true_n = np.where(s_true_nk == 1) # 確認 print(s_true_n[0:5])
[0 0 0 2 1]
s_true_nkは、行がデータ、列がクラスタに対応しているので、np.where()で行(データ)ごとに要素が1のインデックスを取得します。
生成したクラスタに従い観測データ$\mathbf{X} = \{x_{11}, \cdots, x_{ND}\}$を生成します。
# 観測データを生成 x_nd = np.array([ np.random.multivariate_normal( mean=mu_true_kd[k], cov=sigma2_true_kdd[k], size=1 ).flatten() for k in s_true_n ]) # 確認 print(x_nd[0:5])
[[-0.11453912 8.88713788]
[-1.63426849 9.43266429]
[-0.41770276 2.37368869]
[ 5.89613011 -0.87185337]
[-5.34799934 -5.55523625]]
各クラスタに対応するパラメータ$\boldsymbol{\mu}_k,\ \boldsymbol{\Sigma}_k$を持つ多次元ガウス分布に従い乱数(人工データ)を生成します。多次元ガウス分布の乱数は、np.random.multivariate_normal()で生成できます。
データごとにクラスタが異なるので、for文で1データずつs_true_nからクラスタの値を取り出し、そのクラスタに従ってデータ$\mathbf{x}_n$を生成します。平均引数meanにmu_true_kd[k]、分散共分散行列引数covにsigma2_true_kdd[k]を指定します。また1データずつ生成するので、サンプルサイズ引数sizeは1です。
グラフを作成して観測モデルと観測データを確認しましょう。
# 作図用の格子状の点を作成 X_line, Y_line = np.meshgrid(np.arange(-10.0, 10.0, 0.1), np.arange(-10.0, 10.0, 0.1)) x_line = X_line.flatten() y_line = Y_line.flatten() # 観測モデルを計算 z_kline = np.empty((K, len(x_line))) for k in range(K): tmp_z_line = [ multivariate_normal.pdf( x=(x, y), mean=mu_true_kd[k], cov=sigma2_true_kdd[k] ) for x, y in zip(x_line, y_line) ] z_kline[k] = tmp_z_line.copy() Z_true_kline = z_kline.reshape((K, *X_line.shape)) # 確認 print(Z_true_kline.shape)
(3, 200, 200)
等高線グラフ用に格子状の点を用意する必要があります。np.meshgrid()の第1引数にx軸の点($x_1$がとり得る値)、第2引数にy軸の点($x_2$がとり得る値)を渡すことで直交する点(格子状の点)を作成します。
各交点$(x_1, x_2)$の確率密度をmultivariate_normal.pdf()で計算します。
観測モデルと観測データのグラフを作成します。
# 観測データの散布図を作成 fig = plt.figure(figsize=(10, 10)) for k in range(K): plt.contour(X_line, Y_line, Z_true_kline[k]) # 真の観測モデル k_idx = np.where(s_true_n == k) plt.scatter(x_nd[k_idx, 0], x_nd[k_idx, 1], label='cluster'+str(k+1)) # 真の観測データ plt.scatter(mu_true_kd[:, 0], mu_true_kd[:, 1], marker='+', s=100) # 真の平均 plt.suptitle('Gaussian Mixture Model', fontsize=20) plt.title('K=' + str(K) + ', N=' + str(N), loc='left', fontsize=20) plt.xlabel('$x_1$') plt.ylabel('$x_2$') plt.show()
等高線グラフはplt.contour()で、散布図はplt.scatter()で描画できます。
作図時にnp.where(s_true_n == k)で、各クラスタが割り当てられたデータのインデックスを取得しておきます。そのインデックスを添字として用いてx_ndから各クラスタのデータを取り出し、クラスタごとに色を変えた散布図を描きます。
ここまでが観測モデル(ガウス混合モデル)に関する設定です。次に事前分布のパラメータの設定を行います。
・パラメータの設定
観測モデルのパラメータは本来知り得ないものです。そこで初期値を設定し、ギブスサンプリングにより真の値に近付けていきます。
# 観測モデルのパラメータの初期値を指定 mu_kd = np.zeros((K, D)) lambda_kdd = np.array([ np.linalg.inv(np.identity(D) * 100) for _ in range(K) ]) # 混合比率の初期値を指定 pi_k = np.random.choice(np.arange(0.0, 1.0, 0.01), size=K) pi_k /= np.sum(pi_k)
$\boldsymbol{\mu},\ \boldsymbol{\Lambda},\ \boldsymbol{\pi}$をそれぞれmu_kd、lambda_kdd、pi_kとします。
各事前分布のパラメータを設定します。
# muの事前分布のパラメータを指定 beta = 1 m_d = np.array([0.0, 0.0]) # lambdaの事前分布のパラメータを指定 nu = D w_dd = np.identity(D) * 10 # piの事前分布のパラメータを指定 alpha_k = np.repeat(1, K)
観測モデルの平均パラメータ$\boldsymbol{\mu}_k$の事前分布(多次元ガウス分布)の平均パラメータ$\mathbf{m} = \{m_1, \cdots, m_D\}$、精度行列パラメータの係数$\beta$を、それぞれm_d、betaとして値を指定します。
観測モデルの精度行列パラメータ$\boldsymbol{\Lambda}_k$の事前分布(ウィシャート分布)の自由度$\nu$、パラメータ$\mathbf{W} = \{w_{11}, \cdots, w_{DD}\}$を、それぞれnu、w_ddとして値を指定します。ただし、$\nu > D - 1$の値をとり、$\mathbf{W}$は正定値行列です。
混合比率$\boldsymbol{\pi}$の事前分布(ディリクレ分布)のパラメータ$\boldsymbol{\alpha} = \{\alpha_1, \cdots, \alpha_K\}$をalpha_kとして、$\alpha_k > 0$の値を指定します。
以上で推論に必要なモデルに関する設定は完了です。次はギブスサンプリングを実装します。
・ギブスサンプリング
潜在変数(各データのクラスタ)$\mathbf{S}$とパラメータ$\boldsymbol{\mu},\ \boldsymbol{\Lambda},\ \boldsymbol{\pi}$のサンプリングを交互に行います。またサンプルされたパラメータを用いて、繰り返し事後分布のパラメータを計算(事前分布のパラメータを更新)します。
・コード全体
本では縦ベクトルとしているところを、この例では1次元配列(横ベクトル)として扱うため、転置などの処理が異なっている点に注意してください。
・コード(クリックで展開)
# 試行回数を指定 MaxIter = 500 # 推移の確認用の受け皿 trace_s_in = np.zeros((MaxIter, N)) trace_mu_ikd = np.zeros((MaxIter+1, K, D)) trace_lambda_ikdd = np.zeros((MaxIter+1, K, D, D)) trace_mu_ikd[0] = mu_kd.copy() trace_lambda_ikdd[0] = lambda_kdd.copy() # ギブスサンプリング for i in range(MaxIter): # 初期化 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)) nu_hat_k = np.zeros(K) w_hat_kdd = np.zeros((K, D, D)) alpha_hat_k = np.zeros(K) # 潜在変数のパラメータを計算：式(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() 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) nu_hat_k[k] = np.sum(s_nk[:, k]) + nu 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) # 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() # 値を記録 _, trace_s_in[i] = np.where(s_nk == 1) trace_mu_ikd[i+1] = mu_kd.copy() trace_lambda_ikdd[i+1] = lambda_kdd.copy() # 動作確認 print(str(i+1) + ' (' + str(np.round((i + 1) / MaxIter * 100, 1)) + '%)')
1 (0.2%)
2 (0.4%)
3 (0.6%)
4 (0.8%)
5 (1.0%)
(省略)
496 (99.2%)
497 (99.4%)
498 (99.6%)
499 (99.8%)
500 (100.0%)
trace_***は、パラメータや分布の推移を確認するためのオブジェクトです。推論自体には不要です。
・処理の解説
始めに、パラメータ$\boldsymbol{\eta}_n = \{\eta_{n,1}, \cdots, \eta_{n,K}\}$を持つカテゴリ分布に従い各潜在変数$\mathbf{s}_n$をサンプルします(各データにクラスタを割り当てます)。
まずは各データ$\mathbf{x}_n$のクラスタが$k$である確率を表すパラメータ$\eta_{n,k}$を計算します。
# 初期化 eta_nk = np.zeros((N, K)) # 潜在変数のパラメータを計算：式(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() 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) # 確認 print(np.round(eta_nk[0:5], 4))
[[0.0185 0. 0. ]
[0.0001 0.0539 0. ]
[0.0001 0. 0.0161]
[0.0522 0. 0. ]
[0. 0. 0.0052]]
$\eta_{n,k}$は、次の式で計算します。
ただし全てのデータ($n = 1, \cdots, N$)を一度に処理するために、$(\mathbf{x}_n - \boldsymbol{\mu}_k)^{\top} \boldsymbol{\Lambda}_k (\mathbf{x}_n - \boldsymbol{\mu}_k)$の計算を1から$N$の全ての組み合わせで行います(for文の1行目)。この計算結果は、$N \times N$の配列になります。これは例えば、最初の行の最後の要素(Pythonに合わせて言うと0行目の$N-1$番目の要素)は、$(\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$)を加えています。
さらに$\sum_{k=1}^K \eta_{n,k} = 1$となるように、1から$K$の和をとったもので割って正規化する必要があります。
# 正規化 eta_nk /= np.sum(eta_nk, axis=1, keepdims=True) # 確認 print(np.round(eta_nk[0:5], 3)) print(np.sum(eta_nk[0:5], axis=1))
[[1. 0. 0. ]
[0.001 0.999 0. ]
[0.007 0. 0.993]
[0.999 0. 0.001]
[0. 0. 1. ]]
[1. 1. 1. 1. 1.]
正規化の計算を式にすると次のようになります。
パラメータが得られたので、カテゴリ分布に従い潜在変数$\mathbf{S}$を生成します。
# 初期化 s_nk = np.zeros((N, K)) # 潜在変数をサンプル：式(4.93) for n in range(N): s_nk[n] = np.random.multinomial(n=1, pvals=eta_nk[n], size=1).flatten() # 確認 print(s_nk[0:5])
[[1. 0. 0.]
[0. 1. 0.]
[0. 0. 1.]
[1. 0. 0.]
[0. 0. 1.]]
真のクラスタのときと同様に、潜在変数をサンプルします。ただしデータによってどのクラスタが割り当てられるのかの確率が異なるため、for文で1データずつ生成します。
続いて、それぞれ仮定した分布に従いパラメータ$\boldsymbol{\mu},\ \boldsymbol{\Lambda},\ \boldsymbol{\pi}$をサンプルします。ここからの内容は、for文でクラスタごとに行います。
$\boldsymbol{\mu}_k$の事後分布(多次元ガウス分布)の平均パラメータ$\hat{\mathbf{m}}_k$、精度行列パラメータの係数$\hat{\beta}_k$を計算します。
# クラスタを指定 k = 0 # 初期化 beta_hat_k = np.zeros(K) m_hat_kd = np.zeros((K, D)) # 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] # 確認 print(beta_hat_k) print(m_hat_kd)
[114. 0. 0.]
[[-0.37745058 4.43354391]
[ 0. 0. ]
[ 0. 0. ]]
$\hat{\mathbf{m}}_k,\ \hat{\beta}_k$は、それぞれ次の式で計算します。
計算結果をそれぞれbeta_hat_k、m_hat_kdとします。
ちなみに$\sum_{n=1}^N s_{n,k} \mathbf{x}_n$の計算を確認してみると
print(np.round(x_nd[0:5], 2)) print(s_nk[0:5, k]) print(np.round(s_nk[:, k] * x_nd.T, 2).T[0:5])
[[-2.36 0.93]
[-4.59 -4.35]
[ 8.55 -1.12]
[ 1.15 3.53]
[ 5.4 -7.26]]
[1. 0. 0. 1. 0.]
[[-2.36 0.93]
[-0. -0. ]
[ 0. -0. ]
[ 1.15 3.53]
[ 0. -0. ]]
対象としているクラスタのデータ(この例だと$\mathbf{x}_1,\ \mathbf{x}_4$)には1を掛けてそのままの値で、それ以外のクラスタのデータには0を掛けて値が0になっています。つまり対象となるクラスタのデータのみの和をとっていることが分かります。
$\boldsymbol{\lambda}_k$の事後分布(ウィシャート分布)のパラメータ$\hat{\mathbf{W}}_k$、自由度$\hat{\nu}_k$も計算します。
# 初期化 nu_hat_k = np.zeros(K) w_hat_kdd = np.zeros((K, D, D)) # lambdaの事後分布のパラメータを計算：式(4.103) nu_hat_k[k] = np.sum(s_nk[:, k]) + nu 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) # 確認 print(nu_hat_k) print(np.round(w_hat_kdd, 4))
[115. 0. 0.]
[[[0.0011 0. ]
[0. 0.0015]]
[[0. 0. ]
[0. 0. ]]
[[0. 0. ]
[0. 0. ]]]
$\hat{\mathbf{W}}_k,\ \hat{\nu}_k$は、それぞれ次の式で計算します。
計算結果をそれぞれnu_hat_k、w_hat_kddとします。
パラメータが得られたので、ウィシャート分布に従い$\boldsymbol{\lambda}_k$をサンプルします。
# lambdaをサンプル：式(4.102) smp_lambda = wishart.rvs(size = 1, df=nu_hat_k[k], scale=w_hat_kdd[k]) # 確認 print(smp_lambda)
[[0.12919771 0.00601156]
[0.00601156 0.14763307]]
wishart.rvs()で$\boldsymbol{\lambda}_k$を生成(サンプル)します。クラスタごとに生成するので、size引数に1を指定します。また自由度引数dfにnu_hat_k[k]、スケール行列引数scaleにw_hat_kdd[k]を指定します。
多次元ガウス分布に従い$\boldsymbol{\mu}_k$もサンプルします。
# muをサンプル：式(4.98) smp_mu = np.random.multivariate_normal( mean=m_hat_kd[k], cov=np.linalg.inv(beta_hat_k[k] * lambda_kdd[k]), size=1 ) # 確認 print(smp_mu)
[[-0.46719943 4.39865598]]
np.random.multivariate_normal()で$\boldsymbol{\mu}_k$を生成(サンプル)します。こちらもクラスタごとに生成するので、size引数に1を指定します。また平均引数meanにm_hat_kd[k]、分散共分散行列引数covに先ほどサンプルしたlambda_kdd[k]と係数beta_hat_k[k]との積をnp.linalg.inv()で逆行列に変換して指定します。
サイズが1でも2次元配列で出力されるため、flatten()で1次元に変換してから代入しています。
ここまでの処理を全てのクラスタ($k = 1, \cdots, K$)で行います。
$\boldsymbol{\pi}$の事後分布(ディリクレ分布)のパラメータ$\hat{\alpha}_k$を計算します。
# 初期化 alpha_hat_k = np.zeros(K) # 混合比率のパラメータを計算：式(4.45) alpha_hat_k = np.sum(s_nk, axis=0) + alpha_k # 確認 print(alpha_hat_k)
[114. 70. 69.]
$\hat{\alpha}_k$は、次の式で計算します。
計算結果をalpha_hat_kとします。
パラメータが得られたので、ディリクレ分布に従い$\boldsymbol{\pi}$をサンプルします。
# piをサンプル：式(4.44) smp_pi = dirichlet.rvs(size=1, alpha=alpha_hat_k) # 確認 print(smp_pi)
[[0.28038645 0.5373595 0.18225405]]
dirichlet.rvs()で$\boldsymbol{\pi}$を生成(サンプル)します。これは全データに対して共通するパラメータなので、size引数に1を指定します。またパラメータ引数alphaにalpha_hat_kを指定します。
こちらも2次元配列で出力されるため、flatten()で1次元に変換してから代入しています。
以上がギブスサンプリングによる近似推論で行うの個々の処理です。これを指定した回数繰り返し行うことで、各パラメータの値を徐々に真の値に近付けていきます。
・結果の確認
最後にサンプルしたパラメータを用いて、観測モデルの近似事後分布(の確率密度)を計算し、作図します。
# 近似事後分布を計算 z_kline = np.empty((K, len(x_line))) for k in range(K): tmp_z_line = [ multivariate_normal.pdf( x=(x, y), mean=mu_kd[k], cov=np.linalg.inv(lambda_kdd[k]) ) for x, y in zip(x_line, y_line) ] z_kline[k] = tmp_z_line.copy() Z_kline = z_kline.reshape((K, *X_line.shape))
# 各データのクラスタを抽出 _, s_n = np.where(s_nk == 1) # サンプルしたパラメータによる近似事後分布を作図 fig = plt.figure(figsize=(10, 10)) for k in range(K): plt.contour(X_line, Y_line, Z_kline[k]) # 近似事後分布 plt.contour(X_line, Y_line, Z_true_kline[k], linestyles='dotted', alpha=0.5) # 真の観測モデル k_idx = np.where(s_n == k) # クラスタkのインデックスを取得 plt.scatter(x_nd[k_idx, 0], x_nd[k_idx, 1], label='cluster'+str(k+1)) # 観測データ plt.scatter(mu_true_kd[:, 0], mu_true_kd[:, 1], marker='+', s=100, alpha=0.5) # 真の平均 plt.suptitle('Gibbs Sampling', fontsize=20) plt.title('K=' + str(K) + ', N=' + str(N) + ', iter:' + str(MaxIter), loc='left', fontsize=20) plt.xlabel('$x_{1}$') plt.ylabel('$x_{2}$') plt.show()
うまく推定できていることが確認できます。ただしクラスタの順番についてはランダムに決まります。
実際には、始めの頃のサンプルは初期値に依存しているため破棄する期間(burn-in period)を設ける必要がありますが、この資料ではアルゴリズムの理解を焦点にしているためこれ以上は扱いません。
サンプルされた$\boldsymbol{\mu},\ \boldsymbol{\Lambda}$の値の推移を確認しましょう。
# muの推移を確認 fig = plt.figure(figsize=(15, 10)) for k in range(K): for d in range(D): plt.plot(np.arange(MaxIter+1), trace_mu_ikd[:, k, d], label='k=' + str(k+1) + ', d=' + str(d+1)) plt.xlabel('iteration') plt.ylabel('value') plt.legend() plt.show()
# lambdaの推移を確認 fig = plt.figure(figsize=(15, 10)) for k in range(K): for d1 in range(D): for d2 in range(D): plt.plot(np.arange(MaxIter+1), trace_lambda_ikdd[:, k, d1, d2], label='k=' + str(k+1) + ', d=' + str(d1+1) + ', d''=' + str(d2+1)) plt.xlabel('iteration') plt.ylabel('value') plt.legend() plt.show()
以上がギブスサンプリングによる近似推論の処理です。
・おまけ
matplotlib.animationモジュールを利用して、gif画像を作成するコードです。
・コード(クリックで展開)
# 追加ライブラリ import matplotlib.animation as animation
イタレーションごとに、サンプルしたパラメータを用いて近似事後分布の確率密度を計算します。
# 近似事後分布を計算 Z_ikline = np.zeros((MaxIter+1, K, *X_line.shape)) for i in range(MaxIter + 1): z_kline = np.empty((K, len(x_line))) for k in range(K): tmp_z_line = [ multivariate_normal.pdf( x=(x, y), mean=trace_mu_ikd[i, k], cov=np.linalg.inv(trace_lambda_ikdd[i, k]) ) for x, y in zip(x_line, y_line) ] z_kline[k] = tmp_z_line.copy() Z_ikline[i] = z_kline.reshape((K, *X_line.shape)) # 動作確認 print(str(i) + ' (' + str(np.round((i) / (MaxIter) * 100, 1)) + '%)')
0 (0.0%)
1 (0.2%)
2 (0.4%)
3 (0.6%)
4 (0.8%)
5 (1.0%)
(省略)
496 (99.2%)
497 (99.4%)
498 (99.6%)
499 (99.8%)
500 (100.0%)
処理に時間がかかる場合は、描画範囲や点の間隔(X_line, Y_line作成時のnp.meshgrid()の引数の値)を調整してください。
gif画像を作成して、近似事後分布の推移を確認します。
# グラフを初期化 #plt.cla() # グラフを作成 fig = plt.figure(figsize=(12, 12)) fig.suptitle('Gibbs Sampling', fontsize=20) ax = fig.add_subplot(1, 1, 1) # 作図処理を関数として定義 def update(i): # 前フレームのグラフを初期化 ax.cla() # nフレーム目のグラフを描画 for k in range(K): ax.contour(X_line, Y_line, Z_ikline[i, k]) # 近似事後分布 ax.contour(X_line, Y_line, Z_true_kline[k], linestyles='dotted', alpha=0.5) # 真の観測モデル if i > 0: # 初期値以外のとき k_idx = np.where(trace_s_in[i-1] == k) # クラスタkのインデックスを取得 ax.scatter(x_nd[k_idx, 0], x_nd[k_idx, 1], label='cluster'+str(k+1)) # クラスタkの観測データ if i == 0: # 初期値のとき ax.scatter(x_nd[:, 0], x_nd[:, 1]) # 全ての観測データ ax.scatter(mu_true_kd[:, 0], mu_true_kd[:, 1], marker='+', s=100, alpha=0.5) # 真の平均 # グラフの設定 ax.set_title('K=' + str(K) + ', N=' + str(N) + ', iter:' + str(i), loc='left', fontsize=20) ax.set_xlabel('$x_{1}$') ax.set_ylabel('$x_{2}$') # gif画像を作成 ani = animation.FuncAnimation(fig, update, frames=MaxIter + 1, interval=100) ani.save("filename.gif")
詳細をよく理解していないので解説は省略します…。
(重くて貼れなかったので100回までの画像に差し替えました。)
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
今回は先にPythonで組んでみました。そしてRではどこかミスっててまだうまく推定できてません。
この記事の投稿日の朝に公開されたモーニング娘。'20の新曲「ギューされたいだけなのに」です！
ウサギちゃんシンドローム♪
さらに前日に公開されたアンジュルムの新曲「SHAKA SHAKA TO LOVE」です！
歯磨けよ～
【次節の内容】