はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。基本的な内容は「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」になります。「数式」と「プログラム」から理解するのが目標です。
この記事は、4.3.3項の内容です。「観測モデルをポアソン混合分布」、「事前分布をガンマ分布とディリクレ分布」とするポアソン混合モデルに対する変分推論をPythonで実装します。
省略してある内容等ありますので、本とあわせて読んでください。初学者な自分が理解できるレベルまで落として書き下していますので、分かる人にはかなりくどくなっています。同じような立場の人のお役に立てれば幸いです。
【数式読解編】
【他の節の内容】
【この節の内容】
・Pythonでやってみる
人工的に生成したデータを用いて、ベイズ推論を行ってみましょう。アルゴリズム4.2を参考に実装します。
利用するライブラリを読み込みます。
# 4.3.3項で利用するライブラリ import numpy as np from scipy.stats import poisson, gamma # ポアソン分布, ガンマ分布 from scipy.special import psi # ディガンマ関数 import matplotlib.pyplot as plt
・モデルの設定
まずは、観測モデルを設定します。この例では、$K$個の観測モデルをポアソン分布$\mathrm{Poi}(x_n | \lambda_k)$とします。また、各データに対する$K$個の観測モデルの割り当てをカテゴリ分布$\mathrm{Cat}(\mathbf{s}_n | \boldsymbol{\pi})$によって行います。
観測モデルのパラメータと混合比率を設定します。
# K個の真のパラメータを指定 lambda_truth_k = np.array([10.0, 25.0, 40.0]) # 真の混合比率を指定 pi_truth_k = np.array([0.35, 0.25, 0.4]) # クラスタ数を取得 K = len(lambda_truth_k)
観測モデルの$K$個のパラメータ$\boldsymbol{\lambda} = \{\lambda_1, \lambda_2, \cdots, \lambda_K\}$をlambda_truth_kとして、$\lambda_k \geq 0$の値を指定します。
混合比率パラメータ$\boldsymbol{\pi} = (\pi_1, \pi_2, \cdots, \pi_K)$をpi_truth_kとして、$\sum_{k=1}^K \pi_k = 1$、$\pi_k \geq 0$の値を指定します。
この2つのパラメータの値を求めるのがここでの目的です。
設定したモデルをグラフで確認しましょう。グラフ用の点を作成します。
# 作図用のxの点を作成 x_line = np.arange(0, 2 * np.max(lambda_truth_k)) # 確認 print(x_line[:10])
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
作図用に、ポアソン分布に従う変数$x_n$がとり得る非負の整数をnp.arange()で作成してx_lineとします。この例では、0から$\boldsymbol{\lambda}$の最大値の2倍を範囲とします。
ポアソン混合分布の確率を計算します。
# 観測モデルを計算 model_prob = 0.0 for k in range(K): # クラスタkの分布の確率を計算 tmp_prob = poisson.pmf(k=x_line, mu=lambda_truth_k[k]) # K個の分布の加重平均を計算 model_prob += pi_truth_k[k] * tmp_prob
$K$個のパラメータlambda_truth_kを使って、x_lineの値ごとに次の式で確率を計算します。
$K$個のポアソン分布に対して、$\boldsymbol{\pi}$を使って加重平均をとります。
ポアソン分布の確率は、SciPyライブラリのstatsモジュールのpoisson.pmf()で計算できます。データの引数kにx_line、パラメータの引数muにlambda_truth_k[k]を指定します。
この計算をfor文でK回繰り返し行います。各パラメータによって求めた$K$回分の確率をpi_truth_k[k]で割り引いてmodel_probに加算していきます。
計算結果は次のようになります。
# 確認 print(np.round(model_prob[:10], 3))
[0. 0. 0.001 0.003 0.007 0.013 0.022 0.032 0.039 0.044]
観測モデルを作図します。
# 観測モデルを作図 plt.figure(figsize=(12, 9)) plt.bar(x=x_line, height=model_prob) # 真の分布 plt.xlabel('x') plt.ylabel('prob') plt.suptitle('Poisson Mixture Model', size = 20) plt.title('$\lambda=[' + ', '.join([str(lmd) for lmd in lambda_truth_k]) + ']' + ', \pi=[' + ', '.join([str(pi) for pi in pi_truth_k])+ ']$', loc='left') plt.show()
$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} = \{x_1, x_2, \cdots, x_N\}$を生成します。
$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にすることで生成できます。また、確率(パラメータ)の引数pvalsにpi_truth_k、試行回数の引数sizeにNを指定します。生成した$N$個の$K$次元ベクトルをs_truth_nkとします。
結果は次のようになります。
# 確認 print(s_truth_nk[:5])
[[1 0 0]
[1 0 0]
[1 0 0]
[0 0 1]
[0 0 1]]
これが本来は観測できない潜在変数であり、真のクラスタです。各行がデータを表し、値が1の列番号が割り当てられたクラスタ番号を表します。
s_truth_nkから各データのクラスタ番号を抽出します。
# 真のクラスタ番号を抽出 _, s_truth_n = np.where(s_truth_nk == 1) # 確認 print(s_truth_n[:5])
[0 0 0 2 2]
np.where()を使って行(データ)ごとに要素が1の列番号を抽出します。
各データのクラスタに従い$N$個のデータを生成します。
# (観測)データを生成
x_n = np.random.poisson(lam=lambda_truth_k[s_truth_n], size=N)
ポアソン分布に従う乱数は、np.random.poisson()で生成できます。パラメータの引数lamにlambda_truth_k[s_truth_n]、試行回数の引数sizeにNを指定します。各データのクラスタ番号s_truth_nを添字に使うことで、各データのクラスタに対応したパラメータ$\lambda_k$をlambda_truth_kから取り出すことができます。生成した$N$個のデータをx_nとします。式(4.28)を再現して、lamにapply(lambda_truth_k^t(s_truth_nk), 2, prod)を渡しても生成できます。
観測したデータを確認します。
# 確認 print(x_n[:10])
[ 9 12 8 37 46 13 40 45 5 9]
観測データが非負の整数になっているのを確認できます。
観測データをヒストグラムで確認しましょう。
# 観測データのヒストグラムを作成 plt.figure(figsize=(12, 9)) plt.bar(x=x_line, height=model_prob, label='true model', color='white', alpha=1, edgecolor='red', linestyle='--') # 真の分布 plt.bar(x=x_line, height=[np.sum(x_n == x) / len(x_n) for x in x_line], label='observation data') # 観測データ plt.xlabel('x') plt.ylabel('dens') plt.suptitle('Poisson Mixture Model', fontsize=20) plt.title('$N=' + str(N) + ', \lambda=[' + ', '.join([str(lmd) for lmd in lambda_truth_k]) + ']' + ', \pi=[' + ', '.join([str(pi) for pi in pi_truth_k]) + ']$', loc='left') plt.legend() plt.show()
真の分布を重ねて描画しています。
また、クラスタもヒストグラムで確認しましょう。
# 真のクラスタのヒストグラムを作成 plt.figure(figsize=(12, 9)) for k in range(K): plt.bar(x=x_line, height=[np.sum(x_n[s_truth_n == k] == x) for x in x_line], alpha=0.5, label='cluster:' + str(k + 1)) # 真のクラスタ plt.xlabel('x') plt.ylabel('count') plt.suptitle('Poisson Mixture Model', fontsize=20) plt.title('$N=' + str(N) + ', \lambda=[' + ', '.join([str(lmd) for lmd in lambda_truth_k]) + ']' + ', \pi=[' + ', '.join([str(pi) for pi in pi_truth_k]) + ']$', loc='left') plt.legend() plt.show()
どちらの分布もデータが十分に大きいと真の分布に近づきます。
・事前分布の設定
次に、観測モデル(観測データの分布)$p(\mathbf{X} | \boldsymbol{\lambda})$と潜在変数の分布$p(\mathbf{S} | \boldsymbol{\pi})$に対する共役事前分布を設定します。ポアソン分布のパラメータ$\boldsymbol{\lambda}$に対する事前分布$p(\boldsymbol{\lambda})$としてガンマ分布$\mathrm{Gam}(\boldsymbol{\lambda} | a, b)$、カテゴリ分布のパラメータ$\boldsymbol{\pi}$に対する事前分布$p(\boldsymbol{\pi})$としてディリクレ分布$p(\boldsymbol{\pi} | \boldsymbol{\alpha})$を設定します。
$\boldsymbol{\lambda}$の事前分布と$\boldsymbol{\pi}$の事前分布のパラメータ(超パラメータ)を設定します。
# lambdaの事前分布のパラメータを指定 a = 1.0 b = 1.0 # piの事前分布のパラメータを指定 alpha_k = np.repeat(2.0, K)
ガンマ分布のパラメータ$a,\ b$をそれぞれa, bとして、$a > 0,\ b > 0$の値を指定します。
ディリクレ分布のパラメータ$\boldsymbol{\alpha} = (\alpha_1, \alpha_2, \cdots, \alpha_K)$をalpha_kとして、$\alpha_k > 0$の値を指定します。
超パラメータを設定できたので、$\boldsymbol{\lambda}$の事前分布を確認しましょう。グラフ用の点を作成して、確率密度を計算します。
# 作図用のlambdaの点を作成 lambda_line = np.linspace(0.0, 2.0 * np.max(lambda_truth_k), num=1000) # lambdaの事前分布を計算 prior_lambda = gamma.pdf(x=lambda_line, a=a, scale=1 / b)
ガンマ分布に従う変数$\lambda_k$がとり得る0以上の値をnp.linspace()で作成してlambda_lineとします。この例では、0から$\boldsymbol{\lambda}$の最大値の2倍を範囲とします。np.linspace()を使うと指定した要素数で等間隔に切り分けます。np.arange()を使うと切り分ける間隔を指定できます。処理が重い場合は、この値を調整してください。
作成したlambda_lineの値ごとに確率密度を計算します。ガンマ分布の確率密度は、SciPyライブラリのstatsモジュールのgamma.pdf()で計算できます。データの引数xにlambda_line、a引数にa、scale引数にbを指定します。
計算結果は次のようになります。
# 確認 print(np.round(lambda_line[:10], 2)) print(np.round(prior_lambda[:10], 2))
[0. 0.08 0.16 0.24 0.32 0.4 0.48 0.56 0.64 0.72]
[1. 0.92 0.85 0.79 0.73 0.67 0.62 0.57 0.53 0.49]
$\boldsymbol{\lambda}$の事前分布を作図します。
# lambdaの事前分布を作図 plt.figure(figsize=(12, 9)) plt.plot(lambda_line, prior_lambda, label='prior', color='purple') # 事前分布 plt.vlines(x=lambda_truth_k, ymin=0.0, ymax=np.max(prior_lambda), label='true val', color='red', linestyle='--') # 真の値 plt.xlabel('$\lambda$') plt.ylabel('density') plt.suptitle('Gamma Distribution', fontsize=20) plt.title('a=' + str(a) + ', b=' + str(b), loc='left') plt.legend() plt.show()
真の値を垂直線plt.vlines()で描画しています。
$\pi$の事前分布(と事後分布)については省略します。詳しくは、3.2.2項「カテゴリ分布の学習と予測」をご参照ください。
・初期値の設定
潜在変数の近似事後分布の期待値$\mathbb{E}_{q(\mathbf{S})}[s_{n,k}]$をランダムに生成して初期値とします。
潜在変数の近似事後分布の期待値を初期化します。
# 潜在変数の近似事後分布の期待値を初期化 E_s_nk = np.random.rand(N, K) E_s_nk /= np.sum(E_s_nk, axis=1, keepdims=True)
一様乱数の生成関数np.random.rand()を使って0から1の値を生成し、N行K列の配列を作成してE_s_nkとします。ランダムに生成した値を、行ごとの和が1になるように正規化します。(正規化するので乱数の時点で0から1の値である必要はありませんが、確率であることを明示しておきます。)
生成したパラメータを確認します。
# 確認 print(E_s_nk[:5])
[[0.0206441 0.43149541 0.54786049]
[0.3678249 0.07025081 0.56192429]
[0.2804825 0.67591827 0.04359922]
[0.18476096 0.42784736 0.38739168]
[0.38808409 0.42529007 0.18662584]]
行がデータ、列がクラスタであり、それぞれのデータが各クラスタとなる確率を表しています。
E_s_nkを用いて、$\boldsymbol{\lambda}$の近似事後分布のパラメータの初期値を計算します。
# 初期値によるlambdaの近似事後分布のパラメータを計算:式(4.55) a_hat_k = np.sum(E_s_nk.T * x_n, axis=1) + a b_hat_k = np.sum(E_s_nk, axis=0) + b
$\boldsymbol{\lambda}$の近似事後分布のパラメータ$\hat{\mathbf{a}} = \{\hat{a}_1, \hat{a}_2, \cdots, \hat{a}_K\}$、$\hat{\mathbf{b}} = \{\hat{b}_1, \hat{b}_2, \cdots, \hat{b}_K\}$は
で計算して、結果をa_hat_k, b_hat_kとします。
# 確認 print(a_hat_k) print(b_hat_k)
[1994.56195769 2182.79524051 2016.64280181]
[82.6274252 89.37230443 81.00027038]
ちなみに、$\sum_{n=1}^N \mathbb{E}_{q(\mathbf{s}_n)} [s_{n,k}] x_n$は各データ$x_n$を各クラスタとなる確率で割り引いた上でのクラスタごとの総和で、$\sum_{n=1}^N \mathbb{E}_{q(\mathbf{s}_n)} [s_{n,k}]$は各クラスタとなる確率の総和です。
同様に、$\boldsymbol{\pi}$の近似事後分布のパラメータの初期値を計算します。
# 初期値によるpiのパラメータを計算:式(4.58) alpha_hat_k = np.sum(E_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)
[83.6274252 90.37230443 82.00027038]
超パラメータの初期値が得られたので、$\boldsymbol{\lambda}$の近似事後分布を確認しましょう。$K$個のガンマ分布を計算します。
# 初期値によるlambdaの近似事後分布を作図 posterior_lambda_kl = np.empty((K, len(lambda_line))) for k in range(K): posterior_lambda_kl[k] = gamma.pdf(x=lambda_line, a=a_hat_k[k], scale=1 / b_hat_k[k])
$K$個のガンマ分布を計算して、2次元配列に格納していきます。
作成した配列は次のようになります。
# 確認 print(posterior_lambda_kl)
[[0. 0. 0. ... 0. 0. 0.]
[0. 0. 0. ... 0. 0. 0.]
[0. 0. 0. ... 0. 0. 0.]]
$\boldsymbol{\lambda}$の近似事後分布の初期値を作図します。
# 初期値によるlambdaの近似事後分布を作図 plt.figure(figsize=(12, 9)) for k in range(K): plt.plot(lambda_line, posterior_lambda_kl[k], label='cluster:' + str(k + 1)) # 事後分布 plt.vlines(x=lambda_truth_k, ymin=0.0, ymax=np.max(posterior_lambda_kl), color='red', linestyle='--', label='true val') # 真の値 plt.xlabel('$\lambda$') plt.ylabel('density') plt.suptitle('Gamma Distribution', size=20) plt.title('$iter:' + str(0) + ', N=' + str(N) + ', \hat{a}=[' + ', '.join([str(a) for a in np.round(a_hat_k, 1)]) + ']' + ', \hat{b}=[' + ', '.join([str(b) for b in np.round(b_hat_k, 1)]) + ']$', loc='left') plt.legend() plt.show()
超パラメータ$\hat{a}_k,\ \hat{b}_k,\ \hat{\alpha}_k$からパラメータの期待値$\mathbb{E}[\lambda_k],\ \mathbb{E}[\pi_k]$を求めて、混合分布を計算します。
# lambdaの平均値を計算:式(2.59) E_lambda_k = a_hat_k / b_hat_k # piの平均値を計算:式(2.51) E_pi_k = alpha_hat_k / np.sum(alpha_hat_k) # 初期値による混合分布を計算 init_prob = 0.0 for k in range(K): # クラスタkの分布の確率を計算 tmp_prob = poisson.pmf(k=x_line, mu=E_lambda_k[k]) # K個の分布の加重平均を計算 init_prob += E_pi_k[k] * tmp_prob
$\mathbb{E}[\lambda_k]$は、ガンマ分布の期待値(2.59)より
で計算します。また、$\mathbb{E}[\pi_k]$は、ディリクレ分布の期待値(2.51)より
で計算します。
それぞれE_lambda_k, E_pi_kとして、観測モデルと同様にして混合分布を計算します。
超パラメータの初期値から求めた分布を作図します。
# 初期値による分布を作図 plt.figure(figsize=(12, 9)) plt.bar(x=x_line, height=model_prob, label='true model', color='white', alpha=1, edgecolor='red', linestyle='--') # 真の分布 plt.bar(x=x_line, height=init_prob, color='purple') # 初期値による分布 plt.xlabel('x') plt.ylabel('prob') plt.suptitle('Poisson Mixture Model', size = 20) plt.title('$iter:' + str(0) + ', E[\lambda]=[' + ', '.join([str(lmd) for lmd in np.round(E_lambda_k, 2)]) + ']' + ', E[\pi]=[' + ', '.join([str(pi) for pi in np.round(E_pi_k, 2)])+ ']$', loc='left') plt.legend() plt.show()
観測モデルのときと同様に作図します。
・変分推論
変分推論によりパラメータと潜在変数を推定します。
・コード全体
eta_nkについては添字を使って代入するため、予め変数を用意しておく必要があります。
また、後で各パラメータの更新値や分布の推移を確認するために、それぞれ試行ごとにtrace_***に値を保存していきます。
# 試行回数を指定 MaxIter = 50 # 変数を作成 eta_nk = np.zeros((N, K)) # 推移の確認用の受け皿を作成 trace_s_ink = [E_s_nk] trace_a_ik = [list(a_hat_k)] trace_b_ik = [list(b_hat_k)] trace_alpha_ik = [list(alpha_k)] # 変分推論 for i in range(MaxIter): # 潜在変数の近似事後分布のパラメータの各項を計算:式(4.60-4.62) E_lmd_k = a_hat_k / b_hat_k E_ln_lmd_k = psi(a_hat_k) - np.log(b_hat_k) E_ln_pi_k = psi(alpha_hat_k) - psi(np.sum(alpha_hat_k)) for n in range(N): # 潜在変数の近似事後分布のパラメータを計算:式(4.51) tmp_eta_k = np.exp(x_n[n] * E_ln_lmd_k - E_lmd_k + E_ln_pi_k) eta_nk[n] = tmp_eta_k / np.sum(tmp_eta_k) # 正規化 # 潜在変数の近似事後分布の期待値を計算:式(4.59) E_s_nk[n] = eta_nk[n].copy() # lambdaの近似事後分布のパラメータを計算:式(4.55) a_hat_k = np.sum(E_s_nk.T * x_n, axis=1) + a b_hat_k = np.sum(E_s_nk, axis=0) + b # 初期値によるpiのパラメータを計算:式(4.58) alpha_hat_k = np.sum(E_s_nk, axis=0) + alpha_k # 値を記録 trace_s_ink.append(E_s_nk) trace_a_ik.append(list(a_hat_k)) trace_b_ik.append(list(b_hat_k)) trace_alpha_ik.append(list(alpha_hat_k)) # 動作確認 #print(str(i + 1) + ' (' + str(np.round((i + 1) / MaxIter * 100, 1)) + '%)')
・処理の確認
続いて、変分推論で行う処理を確認していきます。
・潜在変数の近似事後分布$q(\boldsymbol{s}_n)$のパラメータの更新
$\mathbf{S}$の近似事後分布のパラメータを計算します。
# 潜在変数の近似事後分布のパラメータの各項を計算:式(4.60-4.62) E_lmd_k = a_hat_k / b_hat_k E_ln_lmd_k = psi(a_hat_k) - np.log(b_hat_k) E_ln_pi_k = psi(alpha_hat_k) - psi(np.sum(alpha_hat_k)) for n in range(N): # 潜在変数の近似事後分布のパラメータを計算:式(4.51) tmp_eta_k = np.exp(x_n[n] * E_ln_lmd_k - E_lmd_k + E_ln_pi_k) eta_nk[n] = tmp_eta_k / np.sum(tmp_eta_k) # 正規化 # 潜在変数の近似事後分布の期待値を計算:式(4.59) E_s_nk[n] = eta_nk[n].copy()
$n$番目のデータの潜在変数$\mathbf{s}_n$の近似事後分布(カテゴリ分布)のパラメータ$\boldsymbol{\eta}_n = (\eta_{n,1}, \eta_{n,2}, \cdots, \eta_{n,K})$を次の式で計算します。
右辺の計算結果をtmp_eta_kとします。各項は
で計算します。これらの項は、データによって値が変わらないため、ループの外で計算できます。さらに、総和が1になるように次の式で正規化します。
計算結果をeta_nkとします。
ちなみに、式(4.51)の$\exp(\cdot)$の括弧の中を$\tilde{\eta}_{n,k}$とすると、正規化の式は$\frac{\exp (\tilde{\eta}_{n,k})}{\sum_{k'=1}^K \exp (\tilde{\eta}_{n,k'})}$となります。これはソフトマックス関数の計算になります。オーバーフローが起きる場合は、「ソフトマックス関数のオーバーフロー対策【ゼロつく1のノート(数学)】 - からっぽのしょこ」を参考にしてください。
更新したパラメータeta_nk[n, ]を用いて、$\boldsymbol{s}_n$を事後分布の期待値を計算します。カテゴリ分布の期待値(2.31)より
結果は次のようになります。
# 確認 print(np.round(E_s_nk[:5], 3))
[[0.999 0.001 0. ]
[0.985 0.015 0. ]
[1. 0. 0. ]
[0. 0.05 0.95 ]
[0. 0.001 0.999]]
この処理をfor()ループで、$N$回繰り返すことで、全ての潜在変数の近似事後分布の期待値が得られます。
更新した$\mathbb{E}_{q(\mathbf{S})} [\mathbf{S}]$を用いて、$\boldsymbol{\lambda}$の近似事後分布と$\boldsymbol{\pi}$の近似事後分布のパラメータ$\hat{\mathbf{a}},\ \hat{\mathbf{b}},\ \hat{\boldsymbol{\alpha}}$を更新します。計算方法は初期値のときと同じです。
さらに、更新した$\hat{\mathbf{a}},\ \hat{\mathbf{b}},\ \hat{\boldsymbol{\alpha}}$を用いて、$\mathbb{E}_{q(\mathbf{S})} [\mathbf{S}]$を更新します。この処理を指定した回数繰り返します。
・推論結果の確認
・最後のパラメータの確認
超パラメータa_hat_kとb_hat_kの最後の更新値を用いて、$\boldsymbol{\lambda}$の近似事後分布を計算します。
# lambdaの近似事後分布を計算 posterior_lambda_kl = np.empty((K, len(lambda_line))) for k in range(K): posterior_lambda_kl[k] = gamma.pdf(x=lambda_line, a=a_hat_k[k], scale=1 / b_hat_k[k]) # lambdaの近似事後分布を作図 plt.figure(figsize=(12, 9)) for k in range(K): plt.plot(lambda_line, posterior_lambda_kl[k], label='cluster:' + str(k + 1)) # 事後分布 plt.vlines(x=lambda_truth_k, ymin=0.0, ymax=np.max(posterior_lambda_kl), color='red', linestyle='--', label='true val') # 真の値 plt.xlabel('$\lambda$') plt.ylabel('density') plt.suptitle('Gamma Distribution', size=20) plt.title('$iter:' + str(MaxIter) + ', N=' + str(N) + ', \hat{a}=[' + ', '.join([str(a) for a in np.round(a_hat_k, 1)]) + ']' + ', \hat{b}=[' + ', '.join([str(b) for b in np.round(b_hat_k, 1)]) + ']$', loc='left') plt.legend() plt.show()
初期値の設定のときと同様にして作図します。
3つのガンマ分布がそれぞれ3つの真の値付近をピークとする分布を推定できています。
また、a_hat_k, b_hat_kとalpha_hat_kを用いて、混合分布を計算します。
# lambdaの平均値を計算:式(2.59) E_lambda_k = a_hat_k / b_hat_k # piの平均値を計算:式(2.51) E_pi_k = alpha_hat_k / np.sum(alpha_hat_k) # 最後の推定値による混合分布を計算 res_prob = 0.0 for k in range(K): # クラスタkの分布の確率を計算 tmp_prob = poisson.pmf(k=x_line, mu=E_lambda_k[k]) # K個の分布の加重平均を計算 res_prob += E_pi_k[k] * tmp_prob # 最後の推定値による分布を作図 plt.figure(figsize=(12, 9)) plt.bar(x=x_line, height=model_prob, label='truth model', color='white', alpha=1, edgecolor='red', linestyle='--') # 真の分布 plt.bar(x=x_line, height=res_prob, color='purple') # 最後の推定値による分布 plt.xlabel('x') plt.ylabel('prob') plt.suptitle('Poisson Mixture Model:Variational Inference', size = 20) plt.title('$iter:' + str(MaxIter) + ', N=' + str(N) + ', E[\lambda]=[' + ', '.join([str(lmd) for lmd in np.round(E_lambda_k, 2)]) + ']' + ', E[\pi]=[' + ', '.join([str(pi) for pi in np.round(E_pi_k, 2)])+ ']$', loc='left') plt.legend() plt.show()
初期値の設定のときと同様にして作図します。
潜在変数の期待値E_s_nkの最後の更新値を用いて、クラスタのヒストグラムを作成します。作図用のデータフレームを作成します。
# 確率が最大のクラスタ番号を抽出 s_n = np.argmax(E_s_nk, axis=1) # 作図用の潜在変数の近似事後分布の期待値を計算:式(4.59) E_s_lk = np.empty((len(x_line), K)) for n in range(len(x_line)): # 潜在変数の近似事後分布のパラメータを計算:式(4.51) tmp_eta_k = np.exp(x_line[n] * E_ln_lmd_k - E_lmd_k + E_ln_pi_k) E_s_lk[n] = tmp_eta_k / np.sum(tmp_eta_k) # 正規化 # 作図用のxのクラスタとなる確率を抽出 E_s_line = E_s_lk[np.arange(len(x_line)), np.argmax(E_s_lk, axis=1)]
np.argmax()を使って行ごとに確率が最大の列番号を抽出して、各データのクラスタ番号s_nとします。
式(4.51)を見ると、観測データ$x_n$の値が同じであれば各クラスタとなる確率も同じ値になるのが分かります。よって、x_lineの値ごとに各クラスタとなる確率を求めてE_s_lkとします。
E_s_lkからは、行ごとに確率が最大の要素(値)を抽出してE_s_lineとします。
作成したデータフレームは次のようになります。
# 確認 print(s_n[:10]) print(x_line) print(np.round(E_s_line, 2))
[0 0 0 2 2 0 2 2 0 0]
[ 0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17.
18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35.
36. 37. 38. 39. 40. 41. 42. 43. 44. 45. 46. 47. 48. 49. 50. 51. 52. 53.
54. 55. 56. 57. 58. 59. 60. 61. 62. 63. 64. 65. 66. 67. 68. 69. 70. 71.
72. 73. 74. 75. 76. 77. 78. 79.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 0.99 0.98 0.96
0.92 0.82 0.65 0.57 0.76 0.88 0.94 0.97 0.97 0.97 0.96 0.94 0.9 0.85
0.78 0.69 0.59 0.53 0.64 0.74 0.82 0.88 0.92 0.95 0.97 0.98 0.99 0.99
0.99 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
1. 1. 1. 1. 1. 1. 1. 1. 1. 1. ]
s_nは各データのクラスタ番号です。E_s_lineはデータ$x_n$がとり得る値x_lineにおいて、各値が一番なる可能性が高いクラスタの確率値です。
確率が最大となったクラスタのヒストグラムを作成します。
# K個の色を指定 colormap_list = ['Reds', 'Greens', 'Blues'] color_list = ['red', 'green', 'blue'] # 最後のクラスタのヒストグラムを作成 plt.figure(figsize=(12, 9)) for k in range(K): plt.bar(x=x_line, height=[np.sum(x_n[s_truth_n == k] == x) for x in x_line], color='white', alpha=1, edgecolor=color_list[k], linestyle='--', label='true cluster:' + str(k + 1)) # 真のクラスタ for k in range(K): cm = plt.get_cmap(colormap_list[k]) # クラスタkのカラーマップを設定 plt.bar(x=x_line, height=[np.sum(x_n[s_n == k] == x) for x in x_line], color=[cm(p) for p in E_s_line], alpha=1, label='cluster:' + str(k + 1)) # 最後のクラスタ plt.xlabel('x') plt.ylabel('count') plt.suptitle('Poisson Mixture Model:Variational Inference', size=20) plt.title('$iter:' + str(MaxIter) + ', N=' + str(N) + ', E[\lambda]=[' + ', '.join([str(lmd) for lmd in np.round(a_hat_k / b_hat_k, 2)]) + ']' + ', E[\pi]=[' + ', '.join([str(pi) for pi in np.round(alpha_hat_k / np.sum(alpha_hat_k), 2)]) + ']$', loc='left') plt.legend() plt.show()
各データに割り当てられたクラスタとなる確率をplt.bar()の引数colorに指定することで、確率によって棒の色の濃淡を調整します。その際に、plt.get_cmap()にカラーマップ名を指定した作成したcm()に、x_lineの各値に対応した確率E_s_lineを指定します。
率の最大値に従って全てのデータのクラスタを決めているので、棒ごとに1つのクラスタになります。また、クラスタの境界付近のデータは、どちらのクラスタになるのか曖昧になるため色が薄くなる(確率が低くなる)のが分かります。
クラスタはランダムに決まるため、クラスタ番号は一致しません。
・超パラメータの推移の確認
超パラメータの更新値の推移を確認します。
・コード(クリックで展開)
各パラメータの学習ごとの値を格納したリストtrace_***をnp.array()でNumPy配列に変換して作図します。
# aの推移を作図 plt.figure(figsize=(12, 9)) for k in range(K): plt.plot(np.arange(MaxIter + 1), np.array(trace_a_ik).T[k], label='cluster:' + str(k + 1)) plt.xlabel('iteration') plt.ylabel('value') plt.suptitle('Variational Inference', size=20) plt.title('$\hat{a}$', loc='left') plt.legend() plt.grid() plt.show()
# bの推移を作図 plt.figure(figsize=(12, 9)) for k in range(K): plt.plot(np.arange(MaxIter + 1), np.array(trace_b_ik).T[k], label='cluster:' + str(k + 1)) plt.xlabel('iteration') plt.ylabel('value') plt.suptitle('Variational Inference', size=20) plt.title('$\hat{b}$', 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).T[k], label='cluster:' + str(k + 1)) plt.xlabel('iteration') plt.ylabel('value') plt.suptitle('Variational Inference', size=20) plt.title('$\hat{\\alpha}$', 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): # 前フレームのグラフを初期化 plt.cla() # i回目のlambdaの近似事後分布を計算 posterior_lambda_k = np.empty((K, len(lambda_line))) for k in range(K): posterior_lambda_k[k] = gamma.pdf(x=lambda_line, a=trace_a_ik[i][k], scale=1 / trace_b_ik[i][k]) # i回目のlambdaの近似事後分布を作図 for k in range(K): plt.plot(lambda_line, posterior_lambda_k[k], label='cluster:' + str(k + 1)) # 事後分布 plt.vlines(x=lambda_truth_k, ymin=0.0, ymax=np.max(posterior_lambda_k), color='red', linestyle='--', label='true val') # 真の値 plt.xlabel('$\lambda$') plt.ylabel('density') plt.suptitle('Gamma Distribution:Variational Inference', size=20) plt.title('$iter:' + str(i) + ', N=' + str(N) + ', \hat{a}=[' + ', '.join([str(a) for a in np.round(trace_a_ik[i], 1)]) + ']' + ', \hat{b}=[' + ', '.join([str(b) for b in np.round(trace_b_ik[i], 1)]) + ']$', loc='left') plt.legend() # gif画像を作成 posterior_anime = animation.FuncAnimation(fig, update_posterior, frames=MaxIter + 1, interval=100) posterior_anime.save("ch4_3_3_Posterior.gif")
・パラメータの期待値による分布の推移
・コード(クリックで展開)
# 画像サイズを指定 fig = plt.figure(figsize=(12, 9)) # 作図処理を関数として定義 def update_model(i): # 前フレームのグラフを初期化 plt.cla() # i回目のパラメータの平均値を計算:式(2.59,2.51) E_lambda_k = np.array(trace_a_ik[i]) / np.array(trace_b_ik[i]) E_pi_k = np.array(trace_alpha_ik[i]) / np.sum(trace_alpha_ik [i]) # i回目の推定値による混合分布を計算 res_prob = 0.0 for k in range(K): # クラスタkの分布の確率を計算 tmp_prob = poisson.pmf(k=x_line, mu=E_lambda_k[k]) # K個の分布の加重平均を計算 res_prob += E_pi_k[k] * tmp_prob # i回目のサンプルによる分布を作図 plt.bar(x=x_line, height=model_prob, label='true model', color='white', alpha=1, edgecolor='red', linestyle='--') # 真の分布 plt.bar(x=x_line, height=res_prob, color='purple') # 推定値による分布 plt.xlabel('x') plt.ylabel('prob') plt.suptitle('Poisson Mixture Model:Variational Inference', size = 20) plt.title('$iter:' + str(i) + ', N=' + str(N) + ', E[\lambda]=[' + ', '.join([str(lmd) for lmd in np.round(E_lambda_k, 2)]) + ']' + ', E[\pi]=[' + ', '.join([str(pi) for pi in np.round(E_pi_k, 2)]) + ']$', loc='left') plt.ylim(0.0, 0.1) plt.legend() # gif画像を作成 model_anime = animation.FuncAnimation(fig, update_model, frames=MaxIter + 1, interval=100) model_anime.save("ch4_3_3_Model.gif")
・潜在変数の近似事後分布の期待値によるヒストグラムの推移
・コード(クリックで展開)
# K個の色を指定 colormap_list = ['Reds', 'Greens', 'Blues'] color_list = ['red', 'green', 'blue'] # 画像サイズを指定 fig = plt.figure(figsize=(12, 9)) # 作図処理を関数として定義 def update_cluster(i): # 超パラメータを取得 a_hat_k = np.array(trace_a_ik[i]) b_hat_k = np.array(trace_b_ik[i]) alpha_hat_k = np.array(trace_alpha_ik[i]) # 潜在変数の近似事後分布のパラメータの各項を計算:式(4.60-4.62) E_lmd_k = a_hat_k / b_hat_k E_ln_lmd_k = psi(a_hat_k) - np.log(b_hat_k) E_ln_pi_k = psi(alpha_hat_k) - psi(np.sum(alpha_hat_k)) # 作図用の潜在変数の近似事後分布の期待値を計算:式(4.59) E_s_lk = np.empty((len(x_line), K)) for n in range(len(x_line)): # 潜在変数の近似事後分布のパラメータを計算:式(4.51) tmp_eta_k = np.exp(x_line[n] * E_ln_lmd_k - E_lmd_k + E_ln_pi_k) E_s_lk[n] = tmp_eta_k / np.sum(tmp_eta_k) # 正規化 # 作図用のxのクラスタとなる確率を抽出 E_s_line = E_s_lk[np.arange(len(x_line)), np.argmax(E_s_lk, axis=1)] # 確率が最大のクラスタ番号を抽出 s_n = np.argmax(trace_s_ink[i], axis=1) # 前フレームのグラフを初期化 plt.cla() # i回目のクラスタの散布図を作成 for k in range(K): plt.bar(x=x_line, height=[np.sum(x_n[s_truth_n == k] == x) for x in x_line], color='white', alpha=1, edgecolor=color_list[k], linestyle='--', label='true:' + str(k + 1)) # 真のクラスタ for k in range(K): # k番目のグラフのカラーマップを設定 cm = plt.get_cmap(colormap_list[k]) plt.bar(x=x_line, height=[np.sum(x_n[s_n == k] == x) for x in x_line], color=[cm(p) for p in E_s_line], alpha=1, label='cluster:' + str(k + 1)) # 推定したクラスタ plt.xlabel('x') plt.ylabel('count') plt.suptitle('Poisson Mixture Model:Variational Inference', size=20) plt.title('$iter:' + str(i) + ', N=' + str(N) + ', E[\lambda]=[' + ', '.join([str(lmd) for lmd in np.round(a_hat_k / b_hat_k, 2)]) + ']' + ', E[\pi]=[' + ', '.join([str(pi) for pi in np.round(alpha_hat_k / np.sum(alpha_hat_k), 2)]) + ']$', loc='left') plt.legend() # gif画像を作成 cluster_anime = animation.FuncAnimation(fig, update_cluster, frames=MaxIter + 1, interval=100) cluster_anime.save("ch4_3_3_Cluster.gif")
(よく理解していないので、animationの解説は省略…とりあえずこれで動きます……)
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
もう少しPyPlot力を上げたい。plt.関数名()のやり方だけでどこまでいけるのか！？
【次節の内容】
