【Python】4.4.3：ガウス混合モデルの変分推論【緑ベイズ入門のノート】

はじめに

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

　この記事は多次元ガウス混合モデルにおける変分推論による近似推論をPythonで実装します。

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

【数式読解編】

www.anarchive-beta.com

【前節の内容】

www.anarchive-beta.com

【他の節一覧】

www.anarchive-beta.com

【この節の内容】

・Pythonでやってみる

　ガウス混合モデルに従い生成したデータを用いて、パラメータを推定してみましょう。

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

# 利用ライブラリ
import numpy as np
from scipy.stats import multivariate_normal
from scipy.special import psi
import matplotlib.pyplot as plt

　作図時にscipyライブラリのstatsモジュールから、多次元ガウス分布multivariate_normalクラスの確率密度の計算メソッドpdf()を使います。またパラメータの計算時にscipyライブラリのspecialモジュールから、ディガンマ関数psi()を使います。

・真の観測モデルの設定

　観測モデル$p(\mathbf{X} | \mathbf{S}, \boldsymbol{\mu}, \boldsymbol{\Lambda})$のパラメータを設定します。作図に関しては、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_{k,1,1}^2, \cdots, \sigma_{k,D,D}^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$の値をとります。このパラメータの値も求めます。

　以上が観測モデルの設定です。続いて、設定したガウス混合モデルに従ってデータを生成します。先に各データが持つクラスタ(潜在変数$\mathbf{S} = \{s_{1,1}, \cdots, s_{N,K}\}$)を生成します。

# 各データの真のクラスタを生成
N = 250
s_true_nk = np.random.multinomial(n=1, pvals=pi_true_k, size=N)

# 確認
print(s_true_nk[0:5])

[[0 0 1]
 [0 1 0]
 [0 0 1]
 [1 0 0]
 [0 0 1]]

　$\boldsymbol{\pi}$をパラメータとするカテゴリ分布に従い、潜在変数$\mathbf{S}$に1から$K$の値を割り当てます。カテゴリ分布の乱数は、多項分布の乱数生成関数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])

[2 1 2 0 2]

　s_true_nkは、行がデータ、列がクラスタに対応しているので、np.where()で行(データ)ごとに要素が1のインデックスを取得します。s_true_nの各要素が各データに対応し、0からK-1の値が1から$K$のクラスタに対応します。

　生成したクラスタに従い、観測データ$\mathbf{X} = \{x_{1,1}, \cdots, x_{N,D}\}$を生成します。

# 観測データを生成
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])

[[11.71152328  0.55523311]
 [-6.01880666 -2.75243653]
 [ 6.42341316  3.3940999 ]
 [ 3.89953429  2.16385938]
 [ 6.86726983  0.73660606]]

　各データに与えられたクラスタのパラメータ$\boldsymbol{\mu}_k,\ \boldsymbol{\Sigma}_k$を持つ多次元ガウス分布に従い、データを生成します。多次元ガウス分布の乱数は、numpy.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_{n,1}$がとり得る値)、第2引数にy軸の点($x_{n,2}$がとり得る値)を渡すことで直交する点(格子状の点)を作成します。

　生成した各交点$(x_{n,1}, x_{n,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_{n,1}$')
plt.ylabel('$x_{n,2}$')
plt.show()

　等高線グラフはplt.contour()で、散布図はplt.scatter()で描画できます。
　作図時にnp.where(s_true_n == k)で、各クラスタが割り当てられたデータのインデックスを取得しておきます。そのインデックスを添字として用いてx_ndから各クラスタのデータを取り出し、クラスタごとに色を変えた散布図を描きます。

ガウス混合モデル

　ここまでが観測モデル(ガウス混合モデル)に関する設定です。次に事前分布のパラメータの設定を行います。

・事前分布のパラメータの設定

　観測モデルのパラメータは本来知り得ないものです。そこで事前分布を設定し、事後分布を求めます(分布推定します)。

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

# lambdaの事前分布のパラメータを指定
nu = D
w_dd = np.identity(D) * 0.05

# piの事前分布のパラメータを指定
alpha_k = np.repeat(1, K)

# 確認
print(w_dd)
print(np.linalg.inv(nu * w_dd)) # 分散

[[0.05 0.  ]
 [0.   0.05]]
[[10.  0.]
 [ 0. 10.]]

　観測モデルの平均パラメータ$\boldsymbol{\mu}_k$の事前分布(多次元ガウス分布)の平均パラメータ$\mathbf{m} = \{m_1, \cdots, m_D\}$、精度行列パラメータの係数$\beta$を、それぞれm_d、betaとして値を指定します。
　観測モデルの精度行列パラメータ$\boldsymbol{\Lambda}_k$の事前分布(ウィシャート分布)の自由度$\nu$、パラメータ$\mathbf{W} = \{w_{1,1}, \cdots, w_{D,D}\}$を、それぞれnu、w_ddとして値を指定します。ただし、$\nu > D - 1$の値をとり、$\mathbf{W}$は正定値行列です。
　混合比率$\boldsymbol{\pi}$の事前分布(ディリクレ分布)のパラメータ$\boldsymbol{\alpha} = \{\alpha_1, \cdots, \alpha_K\}$をalpha_kとして、$\alpha_k > 0$の値を指定します。

　近似事後分布のパラメータの初期値をランダムに設定します。

# mu近似事後分布のパラメータの初期値をランダムに設定
beta_hat_k = np.random.choice(np.arange(0.1, 10.0, 0.1), size=K, replace=True)
m_hat_kd = np.random.choice(
    np.arange(np.min(np.round(x_nd, 1)), np.max(np.round(x_nd, 1)), 0.1), size=(K, D), replace=True
)

# lambda近似事後分布のパラメータの初期値を設定
nu_hat_k = np.repeat(nu, K)
w_hat_kdd = np.array([w_dd for _ in range(K)])

# pi近似事後分布のパラメータの初期値をランダムに設定
alpha_hat_k = np.random.choice(np.arange(0.1, 10.0, 0.1), size=K, replace=True)

# 確認
print(m_hat_kd)

[[ 4.8  1.5]
 [-1.4  1.4]
 [-5.4 -0.8]]

　$\boldsymbol{\mu}_k$の近似事後分布(多次元ガウス分布)の平均パラメータ$\hat{\mathbf{m}} = \{\hat{m}_{1,1}, \cdots, \hat{m}_{K,D}\}$、精度行列パラメータの係数$\hat{\boldsymbol{\beta}} = \{\hat{\beta}_1, \cdots, \hat{\beta}_K\}$を、それぞれm_hat_kd、beta_hat_kとします。
　$\boldsymbol{\Lambda}_k$の近似事後分布(ウィシャート分布)の自由度$\hat{\boldsymbol{\nu}} = \{\hat{\nu}_1, \cdots, \hat{\nu}_K\}$、パラメータ$\hat{\mathbf{W}} = \{\hat{\mathbf{W}}_1, \cdots, \hat{\mathbf{W}}_K\}$、$\hat{\mathbf{W}}_k = \{w_{k,1,1}, \cdots, w_{k,D,D}\}$を、それぞれnu_hat_k、w_hat_kddとします。
　$\boldsymbol{\pi}$の近似事後分布(ディリクレ分布)のパラメータ$\hat{\boldsymbol{\alpha}} = \{\hat{\alpha}_1, \cdots, \hat{\alpha}_K\}$をalpha_hat_kとします。

　m_hat_kdは、観測データx_ndの最小値から最大値の範囲でランダムに値を決めています。beta_hat_kとalpha_hat_kは、範囲を指定してランダムに値を決めています。w_hat_ddkは、(値をランダムに生成するのが面倒なので)nu_hat_kと共に事前分布のパラメータの値を複製しています。

　事前分布と事後分布のパラメータのことを超パラメータ(ハイパーパラメータ)と呼びます。以上で推論に必要なモデルに関する設定は完了です。次は変分推論を実装します。

・変分推論の実装

　各データのクラスタ$\mathbf{S}$の近似事後分布$q(\mathbf{S})$のパラメータの計算(超パラメータの更新)と、パラメータ$\boldsymbol{\mu},\ \boldsymbol{\Lambda},\ \boldsymbol{\pi}$の近似事後分布のパラメータの計算(超パラメータの更新)を交互に行います。

　本では縦ベクトルとしているところを、この例では横ベクトルとして扱うため、転置などの処理が異なっている点に注意してください。(またNotebookの都合上、先にコード全体を載せて後で個々の処理について解説します。)

・コード全体

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

　$q(\mathbf{S})$の計算時に使用する中間変数(オブジェクト)に関して、クラスタごとの計算結果を添字を用いて代入していくため、最初に受け皿としてのオブジェクトを作成しておきます。

# 試行回数を指定
MaxIter = 100

# 途中計算に用いる項の受け皿を作成
ln_eta_nk = np.zeros((N, K))
E_lmd_kdd = np.zeros((K, D, D))
E_ln_det_lmd_k = np.zeros(K)
E_lmd_mu_kd = np.zeros((K, D))
E_mu_lmd_mu_k = np.zeros(K)
E_ln_pi_k = np.zeros(K)

# 推移の確認用の受け皿
trace_E_s_ink = np.zeros((MaxIter, N, K))
trace_E_mu_ikd = np.zeros((MaxIter+1, K, D))
trace_E_lambda_ikdd = np.zeros((MaxIter+1, K, D, D))
trace_E_mu_ikd[0] = m_hat_kd.copy()
trace_E_lambda_ikdd[0] = np.repeat(nu_hat_k, D * D).reshape(K, D, D) * w_hat_kdd

# 変分推論
for i in range(MaxIter):
    
    # Sの近似事後分布のパラメータを計算:式(4.109)
    for k in range(K):
        E_lmd_kdd[k] = nu_hat_k[k] * w_hat_kdd[k]
        E_ln_det_lmd_k[k] = np.sum(psi(0.5 * (nu_hat_k[k] - np.arange(D))))
        E_ln_det_lmd_k[k] += D * np.log(2) + np.log(np.linalg.det(w_hat_kdd[k]))
        E_lmd_mu_kd[k] = np.dot(E_lmd_kdd[k], m_hat_kd[k:k+1].T).flatten()
        E_mu_lmd_mu_k[k] = np.dot(m_hat_kd[k], E_lmd_mu_kd[k]) + D / beta_hat_k[k]
        E_ln_pi_k[k] = psi(alpha_hat_k[k]) - psi(np.sum(alpha_hat_k))
        ln_eta_nk[:, k] = -0.5 * np.diag(
            x_nd.dot(E_lmd_kdd[k]).dot(x_nd.T)
        )
        ln_eta_nk[:, k] += np.dot(x_nd, E_lmd_mu_kd[k:k+1].T).flatten()
        ln_eta_nk[:, k] -= 0.5 * E_mu_lmd_mu_k[k]
        ln_eta_nk[:, k] += 0.5 * E_ln_det_lmd_k[k] + E_ln_pi_k[k]
    tmp_eta_nk = np.exp(ln_eta_nk)
    eta_nk = ((tmp_eta_nk.T + 1e-7) / np.sum(tmp_eta_nk + 1e-7, axis=1)).T # 正規化
    
    # Sの近似事後分布の期待値を計算:式(4.59)
    E_s_nk = eta_nk.copy()
    
    for k in range(K):
        
        # muの近似事後分布のパラメータを計算:式(4.114)
        beta_hat_k[k] = np.sum(E_s_nk[:, k]) + beta
        m_hat_kd[k] = np.sum(E_s_nk[:, k] * x_nd.T, axis=1) + beta * m_d
        m_hat_kd[k] /= beta_hat_k[k]
        
        # lambdaの近似事後分布のパラメータを計算:式(4.118)
        nu_hat_k[k] = np.sum(E_s_nk[:, k]) + nu
        tmp_w_dd = np.dot(E_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)
        
    # piの近似事後分布のパラメータを計算:式(4.58)
    alpha_hat_k = np.sum(E_s_nk, axis=0) + alpha_k
    
    # 観測モデルのパラメータの期待値を記録
    trace_E_s_ink[i] = E_s_nk.copy()
    trace_E_mu_ikd[i+1] = m_hat_kd.copy()
    trace_E_lambda_ikdd[i+1] = np.repeat(nu_hat_k, D * D).reshape(K, D, D) * w_hat_kdd
    
    # 動作確認
    print(str(i+1) + ' (' + str(np.round((i + 1) / MaxIter * 100, 1)) + '%)')

1 (1.0%)
2 (2.0%)
3 (3.0%)
4 (4.0%)
5 (5.0%)
(省略)
96 (96.0%)
97 (97.0%)
98 (98.0%)
99 (99.0%)
100 (100.0%)

　trace_***は、パラメータや分布の推移を確認するためのオブジェクトです。推論自体には不要です。\ 　$q(\boldsymbol{\mu} | \boldsymbol{\Lambda})$の期待値$\mathbb{E}_{q(\boldsymbol{\mu} | \boldsymbol{\Lambda})}[\boldsymbol{\mu}] = \hat{\mathbf{m}}$、$q(\boldsymbol{\Lambda})$の期待値$\mathbb{E}_{q(\boldsymbol{\Lambda})}[\boldsymbol{\Lambda}] = \hat{\boldsymbol{\nu}} \hat{\mathbf{W}}$、$q(\mathbf{S})$の期待値$\mathbb{E}_{q(\mathbf{S})}[\mathbf{S}] = \boldsymbol{\eta}$を記録しておきます。

・処理の解説

　まずは、潜在変数$\mathbf{s}_n$の近似事後分布(カテゴリ分布)のパラメータ$\boldsymbol{\eta}_n = \{\eta_{n,1}, \cdots, \eta_{n,K}\}$を計算します。$\boldsymbol{\eta}_n$の計算式(4.109)の期待値に関する項の計算から行います。

# Sの近似事後分布のパラメータを計算:式(4.109)
for k in range(K):
    E_lmd_kdd[k] = nu_hat_k[k] * w_hat_kdd[k]
    E_ln_det_lmd_k[k] = np.sum(psi(0.5 * (nu_hat_k[k] - np.arange(D))))
    E_ln_det_lmd_k[k] += D * np.log(2) + np.log(np.linalg.det(w_hat_kdd[k]))
    E_lmd_mu_kd[k] = np.dot(E_lmd_kdd[k], m_hat_kd[k:k+1].T).flatten()
    E_mu_lmd_mu_k[k] = np.dot(m_hat_kd[k], E_lmd_mu_kd[k]) + D / beta_hat_k[k]
    E_ln_pi_k[k] = psi(alpha_hat_k[k]) - psi(np.sum(alpha_hat_k))

　それぞれ次の式で計算します。

$$ \begin{align} \mathbb{E}_{q(\boldsymbol{\Lambda}_k)} [ \boldsymbol{\Lambda}_k ] &= \hat{\nu} \hat{\mathbf{W}}_k \tag{4.119}\\ \mathbb{E}_{q(\boldsymbol{\Lambda}_k)} [ \ln |\boldsymbol{\Lambda}_k| ] &= \sum_{d=1}^D \psi \Bigl( \frac{\hat{\nu}_k + 1 - d}{2} \Bigr) + D \ln 2 + \ln |\hat{\mathbf{W}}_k| \tag{4.120}\\ \mathbb{E}_{q(\boldsymbol{\mu}_k, \boldsymbol{\Lambda}_k)} [ \boldsymbol{\Lambda}_k \boldsymbol{\mu}_k ] &= \hat{\nu} \hat{\mathbf{W}}_k \hat{\mathbf{m}}_k \tag{4.121}\\ &= \mathbb{E}_{q(\boldsymbol{\Lambda}_k)} [ \boldsymbol{\Lambda}_k ] \hat{\mathbf{m}}_k \\ \mathbb{E}_{q(\boldsymbol{\mu}_k, \boldsymbol{\Lambda}_k)} [ \boldsymbol{\mu}_k^{\top} \boldsymbol{\Lambda}_k \boldsymbol{\mu}_k ] &= \hat{\nu} \hat{\mathbf{m}}_k^{\top} \hat{\mathbf{W}}_k \hat{\mathbf{m}}_k + \frac{D}{\hat{\beta}_k} \tag{4.122}\\ &= \hat{\mathbf{m}}_k^{\top} \mathbb{E}_{q(\boldsymbol{\mu}_k, \boldsymbol{\Lambda}_k)} [ \boldsymbol{\Lambda}_k \boldsymbol{\mu}_k ] + \frac{D}{\hat{\beta}_k} \\ \mathbb{E}_{q(\boldsymbol{\pi})} [ \ln \pi_k ] &= \psi(\hat{\alpha}_k) - \psi \Bigl( \sum_{k=1}^K \hat{\alpha}_k \Bigr) \tag{4.62} \end{align} $$

　ここで$\psi(\cdot)$はディガンマ関数です。

　これらを用いて$\boldsymbol{\eta} = \{\boldsymbol{\eta}_1, \cdots, \boldsymbol{\eta}_N\}$を計算します。

# Sの近似事後分布のパラメータを計算:式(4.109)
for k in range(K):
    ln_eta_nk[:, k] = -0.5 * np.diag(
        x_nd.dot(E_lmd_kdd[k]).dot(x_nd.T)
    )
    ln_eta_nk[:, k] += np.dot(x_nd, E_lmd_mu_kd[k:k+1].T).flatten()
    ln_eta_nk[:, k] -= 0.5 * E_mu_lmd_mu_k[k]
    ln_eta_nk[:, k] += 0.5 * E_ln_det_lmd_k[k] + E_ln_pi_k[k]
tmp_eta_nk = np.exp(ln_eta_nk)

# 確認
print(np.round(tmp_eta_nk[0:5], 3))

[[0.003 0.    0.   ]
 [0.    0.    0.025]
 [0.004 0.002 0.   ]
 [0.003 0.015 0.   ]
 [0.026 0.001 0.   ]]

　$\eta_{n,k}$は、次の式で計算します。

$$ \begin{align} \eta_{n,k} &\propto \exp \Biggl\{ - \frac{1}{2} \mathbf{x}_n^{\top} \mathbb{E}_{q(\boldsymbol{\Lambda}_k)} [ \boldsymbol{\Lambda}_k ] \mathbf{x}_n - \mathbf{x}_n^{\top} \mathbb{E}_{q(\boldsymbol{\mu}_k, \boldsymbol{\Lambda}_k)} [ \boldsymbol{\Lambda}_k \boldsymbol{\mu}_k ] + \frac{1}{2} \mathbb{E}_{q(\boldsymbol{\mu}_k, \boldsymbol{\Lambda}_k)} [ \boldsymbol{\mu}_k^{\top} \boldsymbol{\Lambda}_k \boldsymbol{\mu}_k ] \Biggr.\\ &\qquad \Biggl. + \frac{1}{2} \mathbb{E}_{q(\boldsymbol{\Lambda}_k)} [ \ln |\boldsymbol{\Lambda}_k| ] + \mathbb{E}_{q(\boldsymbol{\pi})} [ \ln \pi_k ] \Biggr\} \tag{4.109} \end{align} $$

　ただし全てのデータ($n = 1, \cdots, N$)を一度に処理するために、$(\mathbf{x}_n - \boldsymbol{\mu}_k)^{\top} \boldsymbol{\Lambda}_k (\mathbf{x}_n - \boldsymbol{\mu}_k)$の計算を1から$N$の全ての組み合わせで行います(for文の2行目)。この計算結果は、$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$)を加えています。

　さらに$\sum_{k=1}^K \eta_{n,k} = 1$となるように、1から$K$の和をとったもので割って正規化する必要があります。

# 正規化
eta_nk = ((tmp_eta_nk.T + 1e-7) / np.sum(tmp_eta_nk + 1e-7, axis=1)).T

# 確認
print(np.round(eta_nk[0:5], 3))
print(np.sum(eta_nk[0:5], axis=1))

[[1.    0.    0.   ]
 [0.    0.008 0.992]
 [0.644 0.356 0.   ]
 [0.191 0.809 0.   ]
 [0.974 0.026 0.   ]]
[1. 1. 1. 1. 1.]

　正規化の計算を式にすると次のようになります。

$$ \eta_{n,k} \leftarrow \frac{ \eta_{n,k} }{ \sum_{k'=1}^K \eta_{n,k'} } $$

　$q(\mathbf{S})$のパラメータが得られたので、$q(\mathbf{S})$の期待値を計算します。

# Sの近似事後分布の期待値を計算:式(4.59)
E_s_nk = eta_nk.copy()

　カテゴリ分布の期待値(2.31)より、パラメータの値そのままです。

$$ \mathbb{E}_{q(\mathbf{s}_n)}[s_{n,k}] = \eta_{n,k} \tag{4.59} $$

　$\mathbf{S}$の近似事後分布の期待値が求まりました。次はこれを用いて、各パラメータの近似事後分布(のパラメータ)を求めます。ここからの内容は、for文でクラスタごとに行います。

　$\boldsymbol{\mu}_k$の近似事後分布(多次元ガウス分布)の平均パラメータ$\hat{\mathbf{m}}_k$、精度行列パラメータの係数$\hat{\beta}_k$を計算します。

# muの近似事後分布のパラメータを計算:式(4.114)
beta_hat_k[k] = np.sum(E_s_nk[:, k]) + beta
m_hat_kd[k] = np.sum(E_s_nk[:, k] * x_nd.T, axis=1) + beta * m_d
m_hat_kd[k] /= beta_hat_k[k]

# 確認
print(np.round(beta_hat_k, 3))
print(np.round(m_hat_kd, 3))

[ 79.019 118.949  55.063]
[[ 5.306 -2.288]
 [-0.053  4.141]
 [-4.806 -4.982]]

　$\hat{\beta}_k,\ \hat{\mathbf{m}}_k$は、それぞれ次の式で計算します。

$$ \begin{align} \hat{\beta}_k &= \sum_{n=1}^N \mathbb{E}_{q(\mathbf{s}_n)} [s_{n,k}] + \beta \\ \hat{\mathbf{m}}_k &= \frac{ \sum_{n=1}^N \mathbb{E}_{q(\mathbf{s}_n)} [s_{n,k}] \mathbf{x}_n + \beta \mathbf{m} }{ \hat{\beta}_k } \tag{4.114} \end{align} $$

　$\boldsymbol{\lambda}_k$の近似事後分布(ウィシャート分布)のパラメータ$\hat{\mathbf{W}}_k$、自由度$\hat{\nu}_k$も計算します。

# lambdaの近似事後分布のパラメータを計算:式(4.118)
nu_hat_k[k] = np.sum(E_s_nk[:, k]) + nu
tmp_w_dd = np.dot(E_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, 3))

[ 80 119  56]
[[[ 0.002 -0.001]
  [-0.001  0.002]]

 [[ 0.001 -0.   ]
  [-0.     0.001]]

 [[ 0.004  0.001]
  [ 0.001  0.004]]]

　$\hat{\nu}_k,\ \hat{\mathbf{W}}_k$は、それぞれ次の式で計算します。

$$ \begin{align} \hat{\nu}_k &= \sum_{n=1}^N \mathbb{E}_{q(\mathbf{s}_n)} [s_{n,k}] + \nu \\ \hat{\mathbf{W}}_k^{-1} &= \sum_{n=1}^N \mathbb{E}_{q(\mathbf{s}_n)} [s_{n,k}] \mathbf{x}_n \mathbf{x}_n^{\top} + \beta \mathbf{m} \mathbf{m}^{\top} - \hat{\beta}_k \hat{\mathbf{m}}_k \hat{\mathbf{m}}^{\top} + \mathbf{W}^{-1} \tag{4.103} \end{align} $$

　ここまでの処理を全てのクラスタ($k = 1, \cdots, K$)で行います。

　最後に、$\boldsymbol{\pi}$の近似事後分布(ディリクレ分布)のパラメータ$\hat{\alpha}_k$を計算します。

# piの近似事後分布のパラメータを計算:式(4.58)
alpha_hat_k = np.sum(E_s_nk, axis=0) + alpha_k

# 確認
print(np.round(alpha_hat_k, 3))

[ 78.96  118.977  55.063]

　$\hat{\alpha}_k$は、次の式で計算します。

$$ \hat{\alpha}_k = \sum_{n=1}^N \mathbb{E}_{q(\boldsymbol{s}_n)} [s_{n,k}] + \alpha_k \tag{4.58} $$

　以上で各パラメータの近似事後分布(のパラメータ)が求まりました。次の試行ではこれを用いて潜在変数の近似事後分布(のパラメータ)を更新します。

　以上が変分推論による近似推論で行う個々の処理です。これを指定した回数くり返し行うことで、徐々に各パラメータの値を真の値に近付けていきます。

・結果の確認

　最後の試行で求めたパラメータの期待値を用いて、観測モデルの近似事後分布(の確率密度)の期待値を計算します。

# 観測モデルの確率密度を計算
z_kline = np.empty((K, len(x_line)))
for k in range(K):
    tmp_z_line = [
        multivariate_normal.pdf(
            x=(x, y), mean=m_hat_kd[k], cov=np.linalg.inv(nu_hat_k[k] * w_hat_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))

　真の観測モデルと重ねて近似事後分布(の期待値)を描画します。

# カラーマップを指定
cmap_list = ['Blues', 'Oranges', 'Greens']

# 各データのクラスタを抽出
max_p_idx = np.argmax(E_s_nk, axis=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(max_p_idx == k) # クラスタkのインデックスを取得
    plt.scatter(x_nd[k_idx, 0], x_nd[k_idx, 1], c=E_s_nk[k_idx, k], cmap=cmap_list[k], label='cluster'+str(k+1)) # 観測データ
plt.scatter(mu_true_kd[:, 0], mu_true_kd[:, 1], marker='+', s=100, alpha=0.5) # 真の平均
plt.suptitle('Variational Inference', fontsize=20)
plt.title('K=' + str(K) + ', N=' + str(N) + ', iter:' + str(MaxIter), loc='left', fontsize=20)
plt.xlabel('$x_{n,1}$')
plt.ylabel('$x_{n,2}$')
plt.show()

　各データのクラスタ割り当て確率の期待値E_s_nk[n]が最も高いクラスタを、そのデータのクラスタとみなすことにします。np.argmax(E_s_nk, axis=1)で、行(データ)ごとに最大値のインデックスを返します。この出力は0からK-1の値となり、データごとのクラスタ番号に対応します。これを用いて、真の散布図のときと同様にクラスタごとに色分けします。
　ここでは更に、各データがそのクラスタとなる確率をグラデーションで表現します。散布図メソッドplt.scatter()のカラー引数cにE_s_nkを、カラーマップ引数cmapにグラデーションの色を指定します。ただしクラスタごとに別の色にするために、クラスタ数と同じ数の色を設定する必要があります。cmapに指定できる色は「Choosing Colormaps in Matplotlib」を参照してください。

変分推論によるガウス混合モデルの近似事後分布

　出力されたグラフを見ると、クラスタの境界に近づくほど観測データの色が薄く(そのクラスタとなる確率が低く)なっていることが確認できます。各クラスタが持つ平均$\mathbb{E}[{\boldsymbol{\mu}_k}] = \hat{\mathbf{m}}_k$(等高線の中心)から離れたデータであっても、他のクラスタの分布からも離れていれば色が濃く(そのクラスタとなる確率が高く)なっていることも確認できます。
　ただしクラスタの順番はランダムに決まるため、真のクラスタの番号(色)とは異なります。

　$\boldsymbol{\mu},\ \boldsymbol{\Lambda}$の近似事後分布の期待値の推移を確認しましょう。

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

　trace_E_mu_ikd、trace_E_lambda_ikddから各クラスタと次元の値をそれぞれ取り出して、折れ線グラフを描画します。

# muの推移を確認
fig = plt.figure(figsize=(15, 10))
for k in range(K):
    for d in range(D):
        plt.plot(np.arange(MaxIter+1), trace_E_mu_ikd[:, k, d], 
                 label='k=' + str(k+1) + ', d=' + str(d+1))
plt.suptitle('Variational Inference', fontsize=20)
plt.title('$\mathbf{\mu}$' + ': N=' + str(N), loc='left', fontsize=20)
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_E_lambda_ikdd[:, k, d1, d2], 
                 label='k=' + str(k+1) + ', d=' + str(d1+1) + ', d''=' + str(d2+1))
plt.suptitle('Variational Inference', fontsize=20)
plt.title('$\mathbf{\Lambda}$' + ': N=' + str(N), loc='left', fontsize=20)
plt.xlabel('iteration')
plt.ylabel('value')
plt.legend()
plt.show()

$\boldsymbol{\mu}$と$\boldsymbol{\Lambda}$の推移

　収束していることが確認できます。

　以上が変分推論による近似推論の処理です。

・おまけ

　matplotlib.animationモジュールを利用して、gif画像を作成するコードです。

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

# 追加ライブラリ
import matplotlib.animation as animation

　更新回数(試行回数)が増えるに従って、近似事後分布(の平均)が真の観測モデルに近づいていく様子を確認します。

# 描画する回数を指定(変更)
#MaxIter = 100

# 観測モデルを計算
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_E_mu_ikd[i, k], cov=np.linalg.inv(trace_E_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 (1.0%)
2 (2.0%)
3 (3.0%)
4 (4.0%)
5 (5.0%)
(省略)
96 (96.0%)
97 (97.0%)
98 (98.0%)
99 (99.0%)
100 (100.0%)

　処理に時間がかかる場合は、描画する試行回数(MaxIterの値)、または描画範囲や点の間隔(X_line, Y_line作成時のnp.meshgrid()の引数の値)を調整してください。

　gif画像を作成して、分布の推移を確認します。

# カラーマップを指定
cmap_list = ['Blues', 'Oranges', 'Greens']

# グラフを初期化
plt.cla()

# グラフを作成
fig = plt.figure(figsize=(9, 9))
fig.suptitle('Variational Inference', fontsize=20)
ax = fig.add_subplot(1, 1, 1)

# 作図処理を関数として定義
def update(i):
    
    # 前フレームのグラフを初期化
    ax.cla()
    
    # 各データのクラスタを抽出
    max_p_idx = np.argmax(trace_E_s_ink[i], axis=1)
    
    # 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(max_p_idx == k) # クラスタkのインデックスを取得
            ax.scatter(x_nd[k_idx, 0], x_nd[k_idx, 1], c=trace_E_s_ink[i, k_idx, k], 
                       cmap=cmap_list[k], 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_{n,1}$')
    ax.set_ylabel('$x_{n,2}$')

# gif画像を作成
ani = animation.FuncAnimation(fig, update, frames=MaxIter + 1, interval=100)
ani.save("filename")

　詳細をよく理解していないので解説は省略します…。(Notebookだとエラーが表示されますが、たぶんフォルダには出力されています。)

変分推論による近似事後分布の推移

参考文献

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

おわりに

　色々と言葉がごちゃごちゃしてくる。観測モデルのパラメータ$\boldsymbol{\mu},\ \boldsymbol{\Lambda},\ \boldsymbol{\pi}$を求めたい。そこでそれぞれ(近似)事後分布を求める(分布推定する)。パラメータの分布自体を描くのは省略して各分布の期待値(平均？)を求めて、それを用いて推定した観測モデル(？)を求めた(描画した)。この最後のグラフも近似事後分布と言っちゃったけど、本当は何と呼べばいいんだ？

　Rの方ではやり方が分からなかった散布図を確率値でグラデーションにするのができて満足。でもクラスタ数に応じて色の設定を手打ちしなきゃいけないのがまだ不満。

　次で4章ラストです！

【次節の内容】続く