はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。基本的な内容は「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」になります。「数式」と「プログラム」から理解するのが目標です。
この記事は、4.4.4項の内容です。「観測モデルを多次元ガウス混合分布(多変量正規混合分布)」、「事前分布をガウス・ウィシャート分布とディリクレ分布」とするガウス混合モデルに対する崩壊型ギブスサンプリング(周辺化ギブスサンプリング)をR言語で実装します。
省略してある内容等ありますので、本とあわせて読んでください。初学者な自分が理解できるレベルまで落として書き下していますので、分かる人にはかなりくどくなっています。同じような立場の人のお役に立てれば幸いです。
【数式読解編】
【他の節の内容】
【この節の内容】
・Rでやってみる
人工的に生成したデータを用いて、ベイズ推論を行ってみましょう。アルゴリズム4.7を参考に実装します。
利用するパッケージを読み込みます。
# 利用パッケージ library(tidyverse) library(mvnfast)
mvnfastは、多次元ガウス分布に関するパッケージです。多次元ガウス分布に従う乱数生成関数rmvn()と確率密度関数dmvn()、多次元スチューデントのt分布の確率密度関数dmvt()を使います。
・観測モデルの設定
まずは、観測モデルを設定します。この節では、$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 <- matrix( c(5, 35, -20, -10, 30, -20), nrow = K, ncol = D, byrow = TRUE ) # K個の真の分散共分散行列を指定 sigma2_truth_ddk <- array( c(250, 65, 65, 270, 125, -45, -45, 175, 210, -15, -15, 250), dim = c(D, D, K) )
$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}$を精度行列と呼びます。
mu_true_kdの作成の際に、matrix()の引数にbyrow = TRUEを指定すると行方向に要素を並べるため、値を指定しやすくなります。
設定したパラメータは次のようになります。
# 確認
mu_truth_kd
sigma2_truth_ddk
## [,1] [,2] ## [1,] 5 35 ## [2,] -20 -10 ## [3,] 30 -20 ## , , 1 ## ## [,1] [,2] ## [1,] 250 65 ## [2,] 65 270 ## ## , , 2 ## ## [,1] [,2] ## [1,] 125 -45 ## [2,] -45 175 ## ## , , 3 ## ## [,1] [,2] ## [1,] 210 -15 ## [2,] -15 250
混合比率を設定します。
# 真の混合比率を指定 pi_truth_k <- c(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_vec <- seq( min(mu_truth_kd[, 1] - 3 * sqrt(sigma2_truth_ddk[1, 1, ])), max(mu_truth_kd[, 1] + 3 * sqrt(sigma2_truth_ddk[1, 1, ])), length.out = 300) # 作図用のy軸のxの値を作成 x_2_vec <- seq( min(mu_truth_kd[, 2] - 3 * sqrt(sigma2_truth_ddk[2, 2, ])), max(mu_truth_kd[, 2] + 3 * sqrt(sigma2_truth_ddk[2, 2, ])), length.out = 300 ) # 作図用のxの点を作成 x_point_mat <- cbind( rep(x_1_vec, times = length(x_2_vec)), rep(x_2_vec, each = length(x_1_vec)) )
作図用に、2次元ガウス分布に従う変数$\mathbf{x}_n = (x_{n,1}, x_{n,2})$がとり得る点(値の組み合わせ)を作成します。$x_{n,1}$がとり得る値(x軸の値)をseq()で作成してx_1_vecとします。この例では、$K$個のクラスタそれぞれで平均値から標準偏差の3倍を引いた値と足した値を計算して、その最小値から最大値までを範囲とします。length.out引数を使うと指定した要素数で等間隔に切り分けます。by引数を使うと切り分ける間隔を指定できます。処理が重い場合は、この値を調整してください。同様に、$x_{n,2}$がとり得る値(y軸の値)も作成してx_2_vecとします。
x_1_vecとx_2_vecの要素の全ての組み合わせを持つようにマトリクスを作成してx_point_matとします。これは、確率密度を等高線図にする際に格子状の点(2軸の全ての値が直交する点)を渡す必要があるためです。
作成した$\mathbf{x}$の点は次のようになります。
# 確認 x_point_mat[1:5, ]
## [,1] [,2] ## [1,] -53.54102 -67.43416 ## [2,] -53.11622 -67.43416 ## [3,] -52.69142 -67.43416 ## [4,] -52.26662 -67.43416 ## [5,] -51.84182 -67.43416
1列目がx軸の値、2列目がy軸の値に対応しています。
観測モデルの確率密度を計算します。
# 観測モデルを計算 model_density <- 0 for(k in 1:K) { # クラスタkの分布の確率密度を計算 tmp_density <- mvnfast::dmvn( X = x_point_mat, mu = mu_truth_kd[k, ], sigma = sigma2_truth_ddk[, , k] ) # K個の分布の加重平均を計算 model_density <- model_density + pi_truth_k[k] * tmp_density }
$K$個のパラメータmu_truth_kd, sigmga2_truth_ddkを使って、x_point_matの点ごとに、次の混合ガウス分布の式で確率密度を計算します。
$K$個の多次元ガウス分布に対して、$\boldsymbol{\pi}$を使って加重平均をとります。
多次元ガウス分布の確率密度は、mvnfastパッケージのdmvn()で計算できます。$k$番目の分布の場合は、データの引数Xにx_point_mat、平均の引数muにmu_truth_kd[k, ]、分散共分散行列の引数sigmaにsigma2_truth_ddk[, , k]を指定します。
この計算をfor()ループでK回繰り返し行います。各パラメータによって求めた$K$回分の確率密度をそれぞれpi_truth_k[k]で割り引いてmodel_probに加算していきます。
計算結果をデータフレームに格納します。
# 観測モデルをデータフレームに格納 model_df <- tibble( x_1 = x_point_mat[, 1], x_2 = x_point_mat[, 2], density = model_density ) # 確認 head(model_df)
## # A tibble: 6 x 3 ## x_1 x_2 density ## <dbl> <dbl> <dbl> ## 1 -53.5 -67.4 9.07e-13 ## 2 -53.1 -67.4 1.07e-12 ## 3 -52.7 -67.4 1.27e-12 ## 4 -52.3 -67.4 1.50e-12 ## 5 -51.8 -67.4 1.78e-12 ## 6 -51.4 -67.4 2.10e-12
ggplot2パッケージを利用して作図するには、データフレームを渡す必要があります。
観測モデルを作図します。
# 観測モデルを作図 ggplot(model_df, aes(x = x_1, y = x_2, z = density, color = ..level..)) + geom_contour() + # 真の分布 labs(title = "Gaussian Mixture Model", subtitle = paste0("K=", K), x = expression(x[1]), y = expression(x[2]))
等高線グラフはgeom_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 <- rmultinom(n = N, size = 1, prob = pi_truth_k) %>% t()
生成するデータ数$N$をNとして値を指定します。
カテゴリ分布に従う乱数は、多項分布に従う乱数生成関数rmultinom()のsize引数を1にすることで生成できます。また、試行回数の引数nにN、確率(パラメータ)の引数probにpi_truth_kを指定します。生成した$N$個の$K$次元ベクトルをs_truth_nkとします。ただし、クラスタが行でデータが列方向に並んで出力されるので、t()で転置しておきます。
結果は次のようになります。
# 確認 s_truth_nk[1:5, ]
## [,1] [,2] [,3] ## [1,] 0 0 1 ## [2,] 1 0 0 ## [3,] 1 0 0 ## [4,] 1 0 0 ## [5,] 1 0 0
これが本来は観測できない潜在変数であり、真のクラスタです。各行がデータを表し、値が1の列番号が割り当てられたクラスタ番号を表します。
s_truth_nkから各データのクラスタ番号を抽出します。
# クラスタ番号を取得 s_truth_n <- which(t(s_truth_nk) == 1, arr.ind = TRUE) %>% .[, "row"] # 確認 s_truth_n[1:5]
## [1] 3 1 1 1 1
which()を使って行(データ)ごとに要素が1の列番号を抽出します。
各データに割り当てられたクラスタに従い$N$個のデータを生成します。
# (観測)データを生成 x_nd <- matrix(0, nrow = N, ncol = D) for(n in 1:N) { # n番目のデータのクラスタ番号を取得 k <- s_truth_n[n] # n番目のデータを生成 x_nd[n, ] = mvnfast::rmvn(n = 1, mu = mu_truth_kd[k, ], sigma = sigma2_truth_ddk[, , k]) }
各データ$\mathbf{x}_n$に割り当てられたクラスタ($s_{n,k} = 1$である$k$)のパラメータ$\boldsymbol{\mu}_k,\ \boldsymbol{\Sigma}_k$を持つ多次元ガウス分布に従いデータ(値)を生成します。多次元ガウス分布の乱数は、mvnfastパッケージのrmvn()で生成できます。
データごとにクラスタが異なるので、for()で1データずつ処理します。$n$番目のデータのクラスタ番号s_true_n[n]を取り出してkとします。そのクラスタに従ってデータ$\mathbf{x}_n$を生成します。平均の引数muにmu_true_kd[k, ]、分散共分散行列の引数sigmaにsigma2_true_ddk[, , k]を指定します。また1データずつ生成するので、データ数の引数nは1です。
観測データもグラフで確認しましょう。作図用のデータフレームを作成します。
# 観測データと真のクラスタ番号をデータフレームに格納 x_df <- tibble( x_n1 = x_nd[, 1], x_n2 = x_nd[, 2], cluster = as.factor(s_truth_n) ) # 確認 head(x_df)
## # A tibble: 6 x 3 ## x_n1 x_n2 cluster ## <dbl> <dbl> <fct> ## 1 18.5 -9.83 3 ## 2 -10.3 47.1 1 ## 3 17.4 64.6 1 ## 4 18.0 31.5 1 ## 5 -6.63 35.6 1 ## 6 32.8 48.9 1
2次元の観測データに1つのクラスタが割り当てられています。
観測データの散布図を観測モデルに重ねて作図します。
# 観測データの散布図を作成 ggplot() + geom_contour(data = model_df, aes(x = x_1, y = x_2, z = density), linetype = "dashed") + # 真の分布 #geom_contour_filled(data = model_df, aes(x = x_1, y = x_2, z = density, fill = ..level..), # alpha = 0.6, linetype = "dashed") + # 真の分布:(塗りつぶし) geom_point(data = x_df, aes(x = x_n1, y = x_n2, color = cluster)) + # 真のクラスタ labs(title = "Gaussian Mixture Model", subtitle = paste0('N=', N, ", pi=(", paste0(pi_truth_k, collapse = ", "), ")"), x = expression(x[1]), y = expression(x[2]))
散布図はgeom_point()で描画します。クラスタごとに色分けしています。観測データから各データのクラスタも推定します。
・事前分布の設定
次に、観測モデル(観測データの分布)$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 m_d <- rep(0, D) # lambdaの事前分布のパラメータを指定 w_dd <- diag(D) * 0.0005 nu <- D # piの事前分布のパラメータを指定 alpha_k <- rep(2, 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 <- solve(beta * nu * w_dd) %>% diag() %>% sqrt() # 作図用のx軸のmuの値を作成 mu_1_vec <- seq( min(mu_truth_kd[, 1]) - sigma_mu_d[1], max(mu_truth_kd[, 1]) + sigma_mu_d[1], length.out = 250) # 作図用のy軸のmuの値を作成 mu_2_vec <- seq( min(mu_truth_kd[, 2]) - sigma_mu_d[2], max(mu_truth_kd[, 2]) + sigma_mu_d[2], length.out = 250 ) # 作図用のmuの点を作成 mu_point_mat <- cbind( rep(mu_1_vec, times = length(mu_2_vec)), rep(mu_2_vec, each = length(mu_1_vec)) )
$\mathbf{x}$の点のときと同様に作成します。各軸の描画範囲は、真の平均の最小値から事前分布の標準偏差を引いた値から、真の平均の最大値に標準偏差を足した値までとします。事前分布の標準偏差は$(\beta \boldsymbol{\Lambda}_k)^{-1}$の対角成分の平方根です。ただしこの時点では$\boldsymbol{\Lambda}$が未知の値なので、代わりに$\boldsymbol{\Lambda}$の期待値$\mathbb{E}[\boldsymbol{\Lambda}]$を使うことにします。ウィシャート分布の期待値(2.89)より
です。
また、逆行列の計算はsolve()、対角成分の抽出はdiag()、平方根の計算はsqrt()で行えます。
面倒でしたら適当な値を直接指定してください。
$\boldsymbol{\mu}$の事前分布の確率密度を計算して、作図用のデータフレームを作成します。
# muの事前分布をデータフレームに格納 prior_mu_df <- tibble( mu_1 = mu_point_mat[, 1], mu_2 = mu_point_mat[, 2], density = mvnfast::dmvn(X = mu_point_mat, mu = m_d, sigma = solve(beta * nu * w_dd)) )
観測モデルのときと同様に、多次元ガウス分布の確率密度を計算します。ここでも、$\boldsymbol{\mu}$の事前分布の分散共分散行列を$(\beta \mathbb{E}[\boldsymbol{\Lambda}])^{-1}$とします。
作成したデータフレームは次のようになります。
# 確認 head(prior_mu_df)
## # A tibble: 6 x 3 ## mu_1 mu_2 density ## <dbl> <dbl> <dbl> ## 1 -51.6 -51.6 0.0000111 ## 2 -51.2 -51.6 0.0000113 ## 3 -50.7 -51.6 0.0000116 ## 4 -50.3 -51.6 0.0000119 ## 5 -49.8 -51.6 0.0000121 ## 6 -49.3 -51.6 0.0000124
$\boldsymbol{\mu}$の事前分布を作図します。
# 真の平均をデータフレームに格納 mu_df <- tibble( x_1 = mu_truth_kd[, 1], x_2 = mu_truth_kd[, 2] ) # muの事前分布を作図 ggplot() + geom_contour(data = prior_mu_df, aes(x = mu_1, y = mu_2, z = density, color = ..level..)) + # 事前分布 geom_point(data = mu_df, aes(x = x_1, y = x_2), color = "red", shape = 4, size = 5) + # 真の値 labs(title = "Gaussian Distribution", subtitle = paste0("m=(", paste0(m_d, collapse = ", "), ")", ", sigma2_mu=(", paste0(solve(beta * nu * w_dd), collapse = ", "), ")"), x = expression(mu[1]), y = expression(mu[2]))
真の平均値をプロットしておきます。この真の平均を分布推定します。
$\boldsymbol{\Lambda}$と$\boldsymbol{\pi}$の事前分布(と事後分布)については省略します。詳しくは、3.2.2 項「カテゴリ分布の学習と予測」と3.4.2項「多次元ガウス分布の学習と予測：精度が未知の場合」を参照してください。
・初期値の設定
設定した$\boldsymbol{\pi}$の事前分布$p(\boldsymbol{\pi})$に従い潜在変数$\mathbf{S}$を生成して初期値とします。また、生成した$\mathbf{S}$を用いて各パラメータ$\boldsymbol{\mu},\ \boldsymbol{\Lambda},\ \boldsymbol{\pi}$の($N$個全てのデータを用いた)事後分布$p(\boldsymbol{\mu} | \boldsymbol{\Lambda}, \mathbf{X}, \mathbf{S}),\ p(\boldsymbol{\Lambda} | \mathbf{X}, \mathbf{S}),\ p(\boldsymbol{\pi} | \mathbf{X}, \mathbf{S})$のパラメータ(超パラメータ)$\hat{\mathbf{m}},\ \hat{\boldsymbol{\beta}},\ \hat{\mathbf{W}},\ \hat{\boldsymbol{\nu}},\ \hat{\boldsymbol{\alpha}}$を計算して初期値とします。
alpha_kを用いて、$\mathbf{S}$の初期値を生成します。
# 潜在変数を初期化 s_nk <- rmultinom(n = N, size = 1, prob = alpha_k / sum(alpha_k)) %>% t()
真のクラスタのときと同様に生成してs_nkとします。ただし、ここでは各クラスタとなる確率として、$\boldsymbol{\alpha}$から求めた混合比率の期待値($\boldsymbol{\pi}$の事前分布の期待値)$\mathbb{E}[\boldsymbol{\pi}]$を使うことにします。クラスタ$k$となる確率の期待値$\mathbb{E}[\pi_k]$は、ディリクレ分布の期待値(2.51)より
です。
生成した潜在変数は次のようになります。
# 確認 s_nk[1:5, ]
## [,1] [,2] [,3] ## [1,] 0 0 1 ## [2,] 0 0 1 ## [3,] 0 0 1 ## [4,] 1 0 0 ## [5,] 0 1 0
s_nkを用いて、$\boldsymbol{\mu}$の事後分布のパラメータを計算します。
# 初期値によるmuの事後分布のパラメータを計算:式(4.99) beta_hat_k <- colSums(s_nk) + beta m_hat_kd <- t(t(x_nd) %*% s_nk + beta * m_d) / beta_hat_k
$\boldsymbol{\mu}$の条件付き分布($K$個の多次元ガウス分布)の平均パラメータ$\hat{\mathbf{m}} = \{\hat{\mathbf{m}}_1, \hat{\mathbf{m}}_2, \cdots, \hat{\mathbf{m}}_K\}$、$\hat{\mathbf{m}}_k = (\hat{m}_{k,1}, \hat{m}_{k,2}, \cdots, \hat{m}_{k,D})$、精度行列パラメータの係数$\boldsymbol{\hat{\beta}} = \{\hat{\beta}_1, \hat{\beta}_2, \cdots, \hat{\beta}_K\}$の各項を
で計算して、それぞれbeta_hat_k、m_hat_kdとします。この時点では、$N$個全てのデータを使って計算しているため、(崩壊型でない)ギブスサンプリング(4.4.2項)の式(4.99)です。
# 確認
beta_hat_k
m_hat_kd
## [1] 81 92 80 ## [,1] [,2] ## [1,] 6.547434 9.463665 ## [2,] 4.291491 5.433040 ## [3,] 10.695865 5.902112
同様に、$\boldsymbol{\lambda}$の事後分布のパラメータを計算します。
# 初期値によるlambdaの事後分布のパラメータを計算:式(4.103) w_hat_ddk <- array(0, dim = c(D, D, K)) term_m_dd <- beta * matrix(m_d) %*% t(m_d) for(k in 1:K) { term_x_dd <- t(s_nk[, k] * x_nd) %*% x_nd term_m_hat_dd <- beta_hat_k[k] * matrix(m_hat_kd[k, ]) %*% t(m_hat_kd[k, ]) w_hat_ddk[, , k] <- solve(term_x_dd + term_m_dd - term_m_hat_dd + solve(w_dd)) } nu_hat_k <- colSums(s_nk) + nu
$\boldsymbol{\lambda}$の条件付き分布($K$個のウィシャート分布)のパラメータ$\hat{\mathbf{W}} = \{\hat{\mathbf{W}}_1, \hat{\mathbf{W}}_2, \cdots, \hat{\mathbf{W}}_K\}$、$\hat{\mathbf{W}}_k = (\hat{w}_{k,1,1}, \cdots, \hat{w}_{k,D,D})$、自由度$\boldsymbol{\hat{\nu}} = \{\hat{\nu}_1, \hat{\nu}_2, \cdots, \hat{\nu}_K\}$を
で計算して、それぞれnu_hat_k、w_hat_ddkとします。
# 確認
w_hat_ddk
nu_hat_k
## , , 1 ## ## [,1] [,2] ## [1,] 1.908850e-05 1.517377e-06 ## [2,] 1.517377e-06 1.662394e-05 ## ## , , 2 ## ## [,1] [,2] ## [1,] 2.085105e-05 1.780083e-06 ## [2,] 1.780083e-06 1.173980e-05 ## ## , , 3 ## ## [,1] [,2] ## [1,] 2.381185e-05 2.880206e-06 ## [2,] 2.880206e-06 1.184004e-05 ## [1] 82 93 81
$\boldsymbol{\pi}$の事後分布のパラメータを計算します。
# 初期値によるpiの事後分布のパラメータを計算:式(4.45) alpha_hat_k <- colSums(s_nk) + alpha_k
$\boldsymbol{\pi}$の条件付き分布(ディリクレ分布)のパラメータ$\hat{\boldsymbol{\alpha}} = (\hat{\alpha}_1, \hat{\alpha}_2, \cdots, \hat{\alpha}_K)$の各項を
で計算して、alpha_hat_kとします。
# 確認
alpha_hat_k
## [1] 82 93 81
超パラメータの初期値が得られたので、$\boldsymbol{\mu}$の事後分布を確認しましょう。$K$個の事後分布を計算します。
# 初期値によるmuの事後分布をデータフレームに格納 posterior_mu_df <- tibble() for(k in 1:K) { # クラスタkのmuの事後分布を計算 tmp_density <- mvnfast::dmvn( X = mu_point_mat, mu = m_hat_kd[k, ], sigma = solve(beta_hat_k[k] * nu_hat_k[k] * w_hat_ddk[, , k]) ) # クラスタkの分布をデータフレームに格納 tmp_df <- tibble( mu_1 = mu_point_mat[, 1], mu_2 = mu_point_mat[, 2], density = tmp_density, cluster = as.factor(k) ) # 結果を結合 posterior_mu_df <- rbind(posterior_mu_df, tmp_df) }
事前分布のときと同様に計算します。ただし、超パラメータが$K$個になったので、多次元ガウス分布も$K$個になります。$K$回分の計算結果をデータフレームに格納していきます。
クラスタ$k$の事後分布の平均は$\hat{\mathbf{m}}$で、精度行列として$\hat{\beta}_k \mathbb{E}[\boldsymbol{\Lambda}_k]$を使います。事前分布のときと同様に、$\mathbb{E}[\boldsymbol{\Lambda}_k] = \hat{\nu}_k \hat{\mathbf{W}}_k$です。
作成したデータフレームは次のようになります。
# 確認 head(posterior_mu_df)
## # A tibble: 6 x 4 ## mu_1 mu_2 density cluster ## <dbl> <dbl> <dbl> <fct> ## 1 -51.6 -51.6 1.23e-200 1 ## 2 -51.2 -51.6 4.61e-199 1 ## 3 -50.7 -51.6 1.68e-197 1 ## 4 -50.3 -51.6 5.95e-196 1 ## 5 -49.8 -51.6 2.06e-194 1 ## 6 -49.3 -51.6 6.92e-193 1
$\boldsymbol{\mu}$の事後分布を作図します。
# 初期値によるmuの事後分布を作図 ggplot() + geom_contour(data = posterior_mu_df, aes(x = mu_1, y = mu_2, z = density, color = cluster)) + # 事後分布 geom_point(data = mu_df, aes(x = x_1, y = x_2), color = "red", shape = 4, size = 5) + # 真の値 labs(title = "Gaussian Distribution", subtitle = paste0("iter:", 0, ", N=", N), x = expression(mu[1]), y = expression(mu[2]))
m_hat_k, beta_hat_kとnu_hat_k, w_hat_ddkから求めた各パラメータの期待値を用いて、混合分布を計算します。
# 初期値による混合分布を計算 init_density <- 0 for(k in 1:K) { # クラスタkの分布の確率密度を計算 tmp_density <- mvnfast::dmvn( X = x_point_mat, mu = m_hat_kd[k, ], sigma = solve(nu_hat_k[k] * w_hat_ddk[, , k]) ) # K個の分布の加重平均を計算 init_density <- init_density + alpha_hat_k[k] / sum(alpha_hat_k) * tmp_density }
観測モデルのときと同様に、多次元ガウス混合分布の確率密度を計算します。ただし、パラメータ$\boldsymbol{\mu},\ \boldsymbol{\Lambda},\ \boldsymbol{\pi}$は未知なので、代わりに各パラメータの期待値$\mathbb{E}[\boldsymbol{\mu}],\ \mathbb{E}[\boldsymbol{\Lambda}],\ \mathbb{E}[\boldsymbol{\pi}]$を使うことにします。
クラスタ$k$における平均パラメータの期待値$\mathbb{E}[\boldsymbol{\mu}_k]$は、多次元ガウス分布の期待値(2.76)より
です。
また、混合比率パラメータの期待値は、先ほどと同様に$\mathbb{E}[\pi_k] = \frac{\hat{\alpha}_k}{\sum_{k'=1}^K \hat{\alpha}_{k'}}$です。
超パラメータの初期値から求めた分布を作図します。
# 初期値による分布をデータフレームに格納 init_df <- tibble( x_1 = x_point_mat[, 1], x_2 = x_point_mat[, 2], density = init_density ) # 初期値による分布を作図 ggplot() + geom_contour(data = model_df, aes(x = x_1, y = x_2, z = density, color = ..level..), alpha = 0.5, linetype = "dashed") + # 真の分布 geom_contour(data = init_df, aes(x = x_1, y = x_2, z = density, color = ..level..)) + # 初期値による分布 labs(title = "Gaussian Mixture Model", subtitle = paste0("iter:", 0, ", K=", K), x = expression(x[1]), y = expression(x[2]))
観測モデルのときと同様に作図します。
クラスタの初期値の散布図を作成します。
# 観測データとクラスタの初期値をデータフレームに格納 s_df <- tibble( x_n1 = x_nd[, 1], x_n2 = x_nd[, 2], cluster = which(t(s_nk) == 1, arr.ind = TRUE) %>% .[, "row"] %>% as.factor() ) # クラスタの初期値を作図 ggplot() + geom_contour(data = model_df, aes(x = x_1, y = x_2, z = density), color = "red", alpha = 0.5, linetype = "dashed") + # 真の分布 geom_contour_filled(data = init_df, aes(x = x_1, y = x_2, z = density, fill = ..level..), alpha = 0.6) + # 初期値による分布 geom_point(data = s_df, aes(x = x_n1, y = x_n2, color = cluster)) + # クラスタの初期値 labs(title = "Gaussian Mixture Model", subtitle = paste0("iter:", 0, ", N=", N, ", K=", K), x = expression(x[1]), y = expression(x[2]))
データの生成のときと同様にして作図します。
・崩壊型ギブスサンプリング
崩壊型ギブスサンプリングでは、潜在変数$\mathbf{s}_n$のサンプリングと、各パラメータ$\boldsymbol{\mu},\ \boldsymbol{\Lambda},\ \boldsymbol{\pi}$の事後分布$p(\boldsymbol{\mu} | \boldsymbol{\Lambda}, \mathbf{X}_{\backslash n}, \mathbf{S}_{\backslash n}),\ p(\boldsymbol{\Lambda} | \mathbf{X}_{\backslash n}, \mathbf{S}_{\backslash n}),\ p(\boldsymbol{\pi} | \mathbf{X}_{\backslash n}, \mathbf{S}_{\backslash n})$のパラメータ(超パラメータ)$\hat{\mathbf{m}}_{\backslash n},\ \hat{\boldsymbol{\beta}}_{\backslash n},\ \hat{\mathbf{W}}_{\backslash n},\ \hat{\boldsymbol{\nu}}_{\backslash n},\ \hat{\boldsymbol{\alpha}}_{\backslash n}$の更新を交互に行います。
まずは、超パラメータの初期値を用いて1番目の潜在変数$\mathbf{s}_1$の事後分布$p(\mathbf{s}_1 | \mathbf{X}, \mathbf{S}_{\backslash 1})$を計算します。得られた事後分布から$\mathbf{s}_1$サンプルします。次に、1番目を更新した$\mathbf{S}$を用いて超パラメータを更新します。さらに、更新した超パラメータを用いて$\mathbf{s}_2$の事後分布$p(\mathbf{s}_2 | \mathbf{X}, \mathbf{S}_{\backslash 2})$を計算して$\mathbf{s}_2$をサンプルし、超パラメータを更新します。これを$N$個全てのデータで行います。ここまでが1回の試行です。
さらに次の試行では、前回の最後に更新した超パラメータを用いて$\mathbf{s}_1$の事後分布の計算となります。この処理を指定した回数繰り返します。
・コード全体
スチューデントのt分布の確率密度の変数density_st_kに関して、添字を使って代入するため予め変数を用意しておく必要があります。
また、後で各パラメータの更新値や分布の推移を確認するために、それぞれ試行ごとにtrace_***に値を保存していきます。関心のあるパラメータを記録してください。
# 試行回数を指定 MaxIter <- 150 # 途中計算に用いる変数を作成 density_st_k <- rep(0, K) # 推移の確認用の受け皿を初期化 trace_s_in <- matrix(0, nrow = MaxIter + 1, ncol = N) trace_beta_ik <- matrix(0, nrow = MaxIter + 1, ncol = K) trace_m_ikd <- array(0, dim = c(MaxIter + 1, K, D)) trace_w_iddk <- array(0, dim = c(MaxIter + 1, D, D, K)) trace_nu_ik <- matrix(0, nrow = MaxIter + 1, ncol = K) trace_alpha_ik <- matrix(0, nrow = MaxIter + 1, ncol = K) # 初期値を記録 trace_s_in[1, ] <- which(t(s_nk) == 1, arr.ind = TRUE) %>% .[, "row"] trace_beta_ik[1, ] <- beta_hat_k trace_m_ikd[1, , ] <- m_hat_kd trace_w_iddk[1, , , ] <- w_hat_ddk trace_nu_ik[1, ] <- nu_hat_k trace_alpha_ik[1, ] <- alpha_hat_k # 崩壊型ギブスサンプリング for(i in 1:MaxIter) { for(n in 1:N) { # n番目のデータの潜在変数を初期化 s_nk[n, ] <- rep(0, K) # muの事後分布のパラメータを計算:式(4.128) beta_hat_k <- colSums(s_nk) + beta m_hat_kd <- t(t(x_nd) %*% s_nk + beta * m_d) / beta_hat_k # lambdaの事後分布のパラメータを計算:式(4.128) w_hat_ddk <- array(0, dim = c(D, D, K)) term_m_dd <- beta * matrix(m_d) %*% t(m_d) for(k in 1:K) { term_x_dd <- t(s_nk[, k] * x_nd) %*% x_nd term_m_hat_dd <- beta_hat_k[k] * matrix(m_hat_kd[k, ]) %*% t(m_hat_kd[k, ]) w_hat_ddk[, , k] <- solve(term_x_dd + term_m_dd - term_m_hat_dd + solve(w_dd)) } nu_hat_k <- colSums(s_nk) + nu # piの事後分布のパラメータを計算:式(4.73) alpha_hat_k <- colSums(s_nk) + alpha_k # スチューデントのt分布のパラメータを計算 nu_st_hat_k <- 1 - D + nu_hat_k mu_st_hat_kd <- m_hat_kd term_lmd_k <- nu_st_hat_k * beta_hat_k / (1 + beta_hat_k) lambda_st_hat_ddk <- w_hat_ddk * array(rep(term_lmd_k, each = D * D), dim = c(D, D, K)) # スチューデントのt分布の確率密度を計算:式(4.129) for(k in 1:K) { density_st_k[k] = mvnfast::dmvt( X = x_nd[n, ], mu = mu_st_hat_kd[k, ], sigma = solve(lambda_st_hat_ddk[, , k]), df = nu_st_hat_k[k] ) + 1e-7 } # カテゴリ分布のパラメータを計算:式(4.75) eta_k <- alpha_hat_k / sum(alpha_hat_k) # 潜在変数のサンプリング確率を計算:式(4.124) tmp_prob_k <- density_st_k * eta_k prob_s_k <- tmp_prob_k / sum(tmp_prob_k) # 正規化 # n番目のデータの潜在変数をサンプル s_nk[n, ] <- rmultinom(n = 1, size = 1, prob = prob_s_k) %>% as.vector() } # muの事後分布のパラメータを計算:式(4.99) beta_hat_k <- colSums(s_nk) + beta m_hat_kd <- t(t(x_nd) %*% s_nk + beta * m_d) / beta_hat_k # lambdaの事後分布のパラメータを計算:式(4.103) w_hat_ddk <- array(0, dim = c(D, D, K)) term_m_dd <- beta * matrix(m_d) %*% t(m_d) for(k in 1:K) { term_x_dd <- t(s_nk[, k] * x_nd) %*% x_nd term_m_hat_dd <- beta_hat_k[k] * matrix(m_hat_kd[k, ]) %*% t(m_hat_kd[k, ]) w_hat_ddk[, , k] <- solve(term_x_dd + term_m_dd - term_m_hat_dd + solve(w_dd)) } nu_hat_k <- colSums(s_nk) + nu # piの事後分布のパラメータを計算:式(4.45) alpha_hat_k <- colSums(s_nk) + alpha_k # i回目のパラメータを記録 trace_s_in[i + 1, ] <- which(t(s_nk) == 1, arr.ind = TRUE) %>% .[, "row"] trace_beta_ik[i + 1, ] <- beta_hat_k trace_m_ikd[i + 1, , ] <- m_hat_kd trace_w_iddk[i + 1, , , ] <- w_hat_ddk trace_nu_ik[i + 1, ] <- nu_hat_k trace_alpha_ik[i + 1, ] <- alpha_hat_k # 動作確認 #print(paste0(i, ' (', round(i / MaxIter * 100, 1), '%)')) }
・処理の解説
崩壊型ギブスサンプリングで行う処理を確認していきます。
・$n$番目のデータの除去
$n$番目の潜在変数$\mathbf{s}_n$を含めない超パラメータ$\mathbf{S}_{\backslash n}$を用いて、パラメータ$\hat{\mathbf{m}}_{\backslash n},\ \hat{\boldsymbol{\beta}}_{\backslash n},\ \hat{\mathbf{W}}_{\backslash n},\ \hat{\boldsymbol{\nu}}_{\backslash n},\ \hat{\boldsymbol{\alpha}}_{\backslash n}$を計算します。
n番目の潜在変数を0に置き換えます。
# n番目のデータの潜在変数を初期化 s_nk[n, ] <- rep(0, K)
s_nkのn行目を全て0にした上で、初期値のときと同じ計算すると、式(4.128)と式(4.73)
の計算となります。(本当は本にあるように$n$番目のデータに関して足し引きする処理をしたかったのですが、実装できませんでした。)\
・潜在変数のサンプリング
$\mathbf{s}_n$の事後分布(のパラメータ)を計算して、新たな$\mathbf{s}_n$をサンプルします。
$\mathbf{s}_n$の事後分布(サンプリング確率)は
を正規化することで求まります。サンプリングに用いるこの2つの分布を計算しましょう。
m_hat_kd, beta_hat_kとw_hat_kkd, nu_hat_kを用いて、式(4.124)の前の項のパラメータを計算します。
# スチューデントのt分布のパラメータを計算 nu_st_hat_k <- 1 - D + nu_hat_k mu_st_hat_kd <- m_hat_kd term_lmd_k <- nu_st_hat_k * beta_hat_k / (1 + beta_hat_k) lambda_st_hat_ddk <- w_hat_ddk * array(rep(term_lmd_k, each = D * D), dim = c(D, D, K))
$K$個のスチューデントのt分布のパラメータ$\hat{\boldsymbol{\nu}}^{(s)}_{\backslash n} = \{\hat{\nu}^{(s)}_{1 \backslash n}, \hat{\nu}^{(s)}_{2 \backslash n}, \cdots, \hat{\nu}^{(s)}_{K \backslash n}\}$、$\hat{\boldsymbol{\mu}}^{(s)}_{\backslash n} = \{\hat{\boldsymbol{\mu}}^{(s)}_{1 \backslash n}, \hat{\boldsymbol{\mu}}^{(s)}_{2 \backslash n}, \cdots, \hat{\boldsymbol{\mu}}^{(s)}_{K \backslash n}\}$、$\hat{\boldsymbol{\mu}}^{(s)}_{k \backslash n} = (\hat{\mu}^{(s)}_{k,1 \backslash n}, \hat{\mu}^{(s)}_{k,2 \backslash n}, \cdots, \hat{\mu}^{(s)}_{k,D \backslash n})$、$\hat{\boldsymbol{\Lambda}}^{(s)}_{\backslash n} = \{\hat{\boldsymbol{\Lambda}}^{(s)}_{1 \backslash n}, \hat{\boldsymbol{\Lambda}}^{(s)}_{2 \backslash n}, \cdots, \hat{\boldsymbol{\Lambda}}^{(s)}_{K \backslash n}\}$、$\hat{\boldsymbol{\Lambda}}^{(s)}_{k \backslash n} = (\hat{\lambda}^{(s)}_{k,1,1 \backslash n}, \cdots, \hat{\lambda}^{(s)}_{k,D,D \backslash n})$の各項を
で計算して、それぞれnu_st_hat_k, mu_st_hat_kd, lambda_st_hat_kddとします。
# 確認
nu_st_hat_k
mu_st_hat_kd
lambda_st_hat_ddk
## [1] 50 114 88 ## [,1] [,2] ## [1,] -17.542203 -12.62923 ## [2,] 3.943457 33.72711 ## [3,] 25.425326 -16.46922 ## , , 1 ## ## [,1] [,2] ## [1,] 0.004988647 0.001623865 ## [2,] 0.001623865 0.003768091 ## ## , , 2 ## ## [,1] [,2] ## [1,] 0.003889774 -0.001709859 ## [2,] -0.001709859 0.003740646 ## ## , , 3 ## ## [,1] [,2] ## [1,] 0.002766116 0.001080254 ## [2,] 0.001080254 0.003179440
nu_st_hat_k, mu_st_hat_kd, lambda_st_hat_kddを用いて、多次元スチューデントのt分布の確率密度を計算します。
# スチューデントのt分布の確率密度を計算:式(4.129) for(k in 1:K) { density_st_k[k] = mvnfast::dmvt( X = x_nd[n, ], mu = mu_st_hat_kd[k, ], sigma = solve(lambda_st_hat_ddk[, , k]), df = nu_st_hat_k[k] ) + 1e-7 }
クラスタごとに確率密度を計算して、density_st_kに代入していきます。多次元スチューデントのt分布の確率密度は、mvnfastパッケージのdmvt()で計算できます。データの引数Xにn番目のデータx_nd[n, ]、平均の引数muにmu_st_hat_kd[k, ]、スケール行列のsigmaにlambda_st_hat_ddk[, , k]の逆行列、自由度の引数dfにnu_st_hat_k[k]を指定します。
# 確認
density_st_k
## [1] 2.584541e-04 3.088012e-06 4.377623e-06
alpha_hat_kを用いて、式(4.124)の後の項のパラメータを計算します。
# カテゴリ分布のパラメータを計算:式(4.75) eta_k <- alpha_hat_k / sum(alpha_hat_k)
カテゴリ分布のパラメータ$\boldsymbol{\eta}_{\backslash n} = (\eta_{1 \backslash n}, \eta_{2 \backslash n}, \cdots, \eta_{K \backslash n})$の各項を
で計算して、eta_kとします。
# 確認 eta_k sum(eta_k)
## [1] 0.2000000 0.4509804 0.3490196 ## [1] 1
和が1になります。
式(4.124)の計算をして、結果を正規化し$\mathbf{s}_n$のサンプリング確率を求めます。
# 潜在変数のサンプリング確率を計算:式(4.124) tmp_prob_k <- density_st_k * eta_k prob_s_k <- tmp_prob_k / sum(tmp_prob_k) # 正規化
多次元スチューデントのt分布の確率密度とカテゴリ分布のパラメータの積を計算します。(式(4.75)の計算にて正規化を行っているため、tmp_prob_kの時点で和が1になります。なので、ここでの正規化の必要性はありません。)
# 確認 prob_s_k sum(prob_s_k)
## [1] 0.94652192 0.02550080 0.02797728 ## [1] 1
prob_s_kを用いて、新たな$\mathbf{s}_n$を生成します。
# n番目のデータの潜在変数をサンプル s_nk[n, ] <- rmultinom(n = 1, size = 1, prob = prob_s_k) %>% as.vector()
prob_s_kをパラメータとして持つカテゴリ分布(4.124)から$\mathbf{s}_n$をサンプルします。
出力は、全ての要素を0にしたs_nk[n, ]に出力を代入します。
これで$n$番目の潜在変数$\mathbf{s}_n$を更新できました。for()ループで$N$個全てのデータに関して同じ処理を行い全ての潜在変数$\mathbf{S}$を更新します。
次の$n + 1$番目のデータに関する処理では、始めに$\mathbf{S}_{\backslash n+1}$を用いて、超パラメータ$\hat{\mathbf{m}}_{\backslash n+1},\ \hat{\boldsymbol{\beta}}_{\backslash n+1},\ \hat{\mathbf{W}}_{\backslash n+1},\ \hat{\boldsymbol{\nu}}_{\backslash n+1},\ \hat{\boldsymbol{\alpha}}_{\backslash n+1}$を計算します。このとき、更新した$\mathbf{s}_n$を含め更新前の$\mathbf{s}_{n+1}$を除くことで超パラメータを更新しています。
最後に、初期値のときと同様にして、更新した$\mathbf{S}$を用いた超パラメータ$\hat{\mathbf{m}},\ \hat{\boldsymbol{\beta}},\ \hat{\mathbf{W}},\ \hat{\boldsymbol{\nu}},\ \hat{\boldsymbol{\alpha}}$を更新します。これは$N$個全てのデータを用いた超パラメータです。
以上が崩壊型ギブスサンプリングによる近似推論で行うの個々の処理です。これを指定した回数繰り返し行うことで、各パラメータの値を徐々に真の値に近付けていきます。
・推論結果の確認
推論結果を確認していきます。
・最後のパラメータの確認
超パラメータm_hat_k, beta_hat_kとnu_hat_k, w_hat_ddkの最後の更新値を用いて、$\boldsymbol{\mu}$の事後分布を計算します。
# muの事後分布をデータフレームに格納 posterior_mu_df <- tibble() for(k in 1:K) { # クラスタkのmuの事後分布を計算 tmp_density <- mvnfast::dmvn( X = mu_point_mat, mu = m_hat_kd[k, ], sigma = solve(beta_hat_k[k] * nu_hat_k[k] * w_hat_ddk[, , k]) ) # クラスタkの分布をデータフレームに格納 tmp_df <- tibble( mu_1 = mu_point_mat[, 1], mu_2 = mu_point_mat[, 2], density = tmp_density, cluster = as.factor(k) ) # 結果を結合 posterior_mu_df <- rbind(posterior_mu_df, tmp_df) } # muの事後分布を作図 ggplot() + geom_contour(data = posterior_mu_df, aes(x = mu_1, y = mu_2, z = density, color = cluster)) + # 事後分布 geom_point(data = mu_df, aes(x = x_1, y = x_2), color = "red", shape = 4, size = 5) + # 真の値 labs(title = "Gaussian Distribution:Collapsed Gibbs Sampling", subtitle = paste0("iter:", MaxIter, ", N=", N), x = expression(mu[1]), y = expression(mu[2]))
変数名やタイトル以外は「初期値の設定」のときと同じコードです。
$K$個のガウス分布がそれぞれ真の平均付近をピークとする分布を推定できています。
m_hat_k, beta_hat_kとnu_hat_k, w_hat_ddkから求めた各パラメータの期待値$\mathbb{E}[\boldsymbol{\mu}],\ \mathbb{E}[\boldsymbol{\Lambda}],\ \mathbb{E}[\boldsymbol{\pi}]$を用いて、混合分布を計算します。
# 最後に更新したパラメータの期待値による混合分布を計算 res_density <- 0 for(k in 1:K) { # クラスタkの分布の確率密度を計算 tmp_density <- mvnfast::dmvn( X = x_point_mat, mu = m_hat_kd[k, ], sigma = solve(nu_hat_k[k] * w_hat_ddk[, , k]) ) # K個の分布の加重平均を計算 res_density <- res_density + alpha_hat_k[k] / sum(alpha_hat_k) * tmp_density } # 最終的な分布をデータフレームに格納 res_df <- tibble( x_1 = x_point_mat[, 1], x_2 = x_point_mat[, 2], density = res_density ) # 最終的な分布を作図 ggplot() + geom_contour(data = model_df, aes(x = x_1, y = x_2, z = density, color = ..level..), alpha = 0.5, linetype = "dashed") + # 真の分布 geom_point(data = mu_df, aes(x = x_1, y = x_2), color = "red", shape = 4, size = 5) + # 真の平均 geom_contour(data = res_df, aes(x = x_1, y = x_2, z = density, color = ..level..)) + # 期待値による分布 geom_point(data = x_df, aes(x = x_n1, y = x_n2)) + # 観測データ labs(title = "Gaussian Mixture Model:Collapsed Gibbs Sampling", subtitle = paste0("iter:", MaxIter, ", N=", N, ", K=", K), x = expression(x[1]), y = expression(x[2]))
変数名やタイトル以外は「初期値の設定」のときと同じコードです。
真の分布に近い分布になっています。
最後にサンプルしたクラスタの散布図を作成します。
# 観測データと最後にサンプルしたクラスタ番号をデータフレームに格納 s_df <- tibble( x_n1 = x_nd[, 1], x_n2 = x_nd[, 2], cluster = which(t(s_nk) == 1, arr.ind = TRUE) %>% .[, "row"] %>% as.factor() ) # 最終的なクラスタを作図 ggplot() + geom_contour(data = model_df, aes(x = x_1, y = x_2, z = density), color = "red", alpha = 0.5, linetype = "dashed") + # 真の分布 geom_point(data = mu_df, aes(x = x_1, y = x_2), color = "red", shape = 4, size = 5) + # 真の平均 geom_contour_filled(data = res_df, aes(x = x_1, y = x_2, z = density, fill = ..level..), alpha = 0.6) + # 期待値による分布 geom_point(data = s_df, aes(x = x_n1, y = x_n2, color = cluster)) + # サンプルしたクラスタ labs(title = "Gaussian Mixture Model:Collapsed Gibbs Sampling", subtitle = paste0("iter:", MaxIter, ", N=", N, ", E[pi]=(", paste0(round(alpha_hat_k / sum(alpha_hat_k), 3), collapse = ", "), ")"), x = expression(x[1]), y = expression(x[2]))
タイトル以外は「初期値の設定」のときと同じコードです。
クラスタリングできていることが確認できます。ただし、クラスタの番号についてはランダムに決まるため真のクラスタの番号と一致しません。
・超パラメータの推移の確認
超パラメータ$\hat{\mathbf{m}},\ \hat{\boldsymbol{\beta}},\ \hat{\boldsymbol{\nu}}, \hat{\mathbf{W}},\ \hat{\boldsymbol{\alpha}}$の更新値の推移を確認します。
・コード(クリックで展開)
学習ごとの値を格納した各パラメータのマトリクスtrace_***をtidyr::pivot_longer()で縦長のデータフレームに変換して作図します。
# mの推移を作図 dplyr::as_tibble(trace_m_ikd) %>% # データフレームに変換 magrittr::set_names( paste0(rep(paste0("k=", 1:K), D), rep(paste0(", d=", 1:D), each = K)) ) %>% # 列名として次元情報を付与 cbind(iteration = 1:(MaxIter + 1)) %>% # 試行回数列を追加 tidyr::pivot_longer( cols = -iteration, # 変換しない列 names_to = "dim", # 現列名を格納する列名 values_to = "value" # 現要素を格納する列名 ) %>% # 縦持ちに変換 ggplot(aes(x = iteration, y = value, color = dim)) + geom_line() + labs(title = "Collapsed Gibbs Sampling", subtitle = expression(hat(m)))
# betaの推移を作図 dplyr::as_tibble(trace_beta_ik) %>% # データフレームに変換 cbind(iteration = 1:(MaxIter + 1)) %>% # 試行回数の列を追加 tidyr::pivot_longer( cols = -iteration, # 変換しない列 names_to = "cluster", # 現列名を格納する列名 names_prefix = "V", # 現列名の頭から取り除く文字列 names_ptypes = list(cluster = factor()), # 現列名を値とする際の型 values_to = "value" # 現セルを格納する列名 ) %>% # 縦持ちに変換 ggplot(aes(x = iteration, y = value, color = cluster)) + geom_line() + # 推定値 labs(title = "Collapsed Gibbs Sampling", subtitle = expression(hat(beta)))
# nuの推移を作図 dplyr::as_tibble(trace_nu_ik) %>% # データフレームに変換 cbind(iteration = 1:(MaxIter + 1)) %>% # 試行回数の列を追加 tidyr::pivot_longer( cols = -iteration, # 変換しない列 names_to = "cluster", # 現列名を格納する列名 names_prefix = "V", # 現列名の頭から取り除く文字列 names_ptypes = list(cluster = factor()), # 現列名を値とする際の型 values_to = "value" # 現セルを格納する列名 ) %>% # 縦持ちに変換 ggplot(aes(x = iteration, y = value, color = cluster)) + geom_line() + # 推定値 labs(title = "Collapsed Gibbs Sampling", subtitle = expression(hat(nu)))
# wの推移を作図 dplyr::as_tibble(trace_w_iddk) %>% # データフレームに変換 magrittr::set_names( paste0( rep(paste0("d=", 1:D), times = D * K), rep(rep(paste0(", d=", 1:D), each = D), times = K), rep(paste0(", k=", 1:K), each = D * D) ) ) %>% # 列名として次元情報を付与 cbind(iteration = 1:(MaxIter + 1)) %>% # 試行回数列を追加 tidyr::pivot_longer( cols = -iteration, # 変換しない列 names_to = "dim", # 現列名を格納する列名 values_to = "value" # 現要素を格納する列名 ) %>% # 縦持ちに変換 ggplot(aes(x = iteration, y = value, color = dim)) + geom_line(alpha = 0.5) + labs(title = "Collapsed Gibbs Sampling", subtitle = expression(hat(w)))
# alphaの推移を作図 dplyr::as_tibble(trace_alpha_ik) %>% # データフレームに変換 cbind(iteration = 1:(MaxIter + 1)) %>% # 試行回数の列を追加 tidyr::pivot_longer( cols = -iteration, # 変換しない列 names_to = "cluster", # 現列名を格納する列名 names_prefix = "V", # 現列名の頭から取り除く文字列 names_ptypes = list(cluster = factor()), # 現列名を値とする際の型 values_to = "value" # 現セルを格納する列名 ) %>% # 縦持ちに変換 ggplot(aes(x = iteration, y = value, color = cluster)) + geom_line() + # 推定値 labs(title = "Collapsed Gibbs Sampling", subtitle = expression(hat(alpha)))
・おまけ：アニメーションによる推移の確認
gganimateパッケージを利用して、各分布の推移のアニメーション(gif画像)を作成するためのコードです。
# 追加パッケージ library(gganimate)
・平均パラメータの事後分布の推移
・コード(クリックで展開)
$\boldsymbol{\mu}$の事後分布の推移を確認します。作図用のデータフレームを作成します。
# 作図用のデータフレームを作成 trace_posterior_df <- tibble() for(i in 1:(MaxIter + 1)) { for(k in 1:K) { # クラスタkのmuの事後分布を計算 tmp_density <- mvnfast::dmvn( X = mu_point_mat, mu = trace_m_ikd[i, k, ], sigma = solve(trace_beta_ik[i, k] * trace_nu_ik[i, k] * trace_w_iddk[i, , , k]) ) # クラスタkの分布をデータフレームに格納 tmp_df <- tibble( mu_1 = mu_point_mat[, 1], mu_2 = mu_point_mat[, 2], density = tmp_density, cluster = as.factor(k), iteration = as.factor(i - 1) ) # 結果を結合 trace_posterior_df <- rbind(trace_posterior_df, tmp_df) } # 動作確認 #print(paste0(i - 1, ' (', round((i - 1) / MaxIter * 100, 1), '%)')) }
各試行における計算結果を同じデータフレームに格納していく必要があります。処理に時間がかかる場合は、描画するフレーム数(1:(MaxIter + 1))や、mu_point_matの範囲や点の間隔を調整してください。
分布の推移をgif画像で出力します。
# muの事後分布を作図 trace_posterior_graph <- ggplot() + geom_contour(data = trace_posterior_df, aes(x = mu_1, y = mu_2, z = density, color = cluster)) + # 事後分布 geom_point(data = mu_df, aes(x = x_1, y = x_2), color = "red", shape = 4, size = 5) + # 真の値 gganimate::transition_manual(iteration) + # フレーム labs(title = "Gaussian Distribution:Collapsed Gibbs Sampling", subtitle = paste0("iter:{current_frame}", ", N=", N), x = expression(mu[1]), y = expression(mu[2])) # gif画像を作成 gganimate::animate(trace_posterior_graph, nframes = MaxIter + 1, fps = 10)
各フレームの順番を示す列をgganimate::transition_manual()に指定します。初期値を含むため、フレーム数はMaxIter + 1です。
・パラメータの期待値による分布とクラスタの推移
・コード(クリックで展開)
各パラメータの期待値を用いた分布とクラスタのサンプルの推移を確認します。同様に、作図用のデータフレームを作成します。
# 作図用のデータフレームを作成 trace_model_df <- tibble() trace_cluster_df <- tibble() for(i in 1:(MaxIter + 1)) { # i回目の混合分布を計算 res_density <- 0 for(k in 1:K) { # クラスタkの分布の確率密度を計算 tmp_density <- mvnfast::dmvn( X = x_point_mat, mu = trace_m_ikd[i, k, ], sigma = solve(trace_nu_ik[i, k] * trace_w_iddk[i, , , k]) ) # K個の分布の加重平均を計算 res_density <- res_density + trace_alpha_ik[i, k] / sum(trace_alpha_ik[i, ]) * tmp_density } # i回目の分布をデータフレームに格納 res_df <- tibble( x_1 = x_point_mat[, 1], x_2 = x_point_mat[, 2], density = res_density, iteration = as.factor(i - 1), label = paste0( "iter:", i - 1, ", N=", N, ", pi=(", paste0(round(trace_alpha_ik[i, ] / sum(trace_alpha_ik[i, ]), 3), collapse = ", "), ")" ) %>% as.factor() ) # 結果を結合 trace_model_df <- rbind(trace_model_df, res_df) # 観測データとi回目のクラスタ番号をデータフレームに格納 s_df <- tibble( x_n1 = x_nd[, 1], x_n2 = x_nd[, 2], cluster = as.factor(trace_s_in[i, ]), iteration = as.factor(i - 1), label = paste0( "iter:", i - 1, ", N=", N, ", pi=(", paste0(round(trace_alpha_ik[i, ] / sum(trace_alpha_ik[i, ]), 3), collapse = ", "), ")" ) %>% as.factor() ) # 結果を結合 trace_cluster_df <- rbind(trace_cluster_df, s_df) # 動作確認 #print(paste0(i - 1, ' (', round((i - 1) / MaxIter * 100, 1), '%)')) }
・パラメータの期待値による分布の推移
# 分布の推移を作図 trace_model_graph <- ggplot() + geom_contour(data = model_df, aes(x = x_1, y = x_2, z = density, color = ..level..), alpha = 0.5, linetype = "dashed") + # 真の分布 geom_point(data = mu_df, aes(x = x_1, y = x_2), color = "red", shape = 4, size = 5) + # 真の平均 geom_contour(data = trace_model_df, aes(x = x_1, y = x_2, z = density, color = ..level..)) + # 期待値による分布 geom_point(data = trace_cluster_df, aes(x = x_n1, y = x_n2)) + # 観測データ gganimate::transition_manual(iteration) + # フレーム labs(title = "Gaussian Mixture Model:Collapsed Gibbs Sampling", subtitle = paste0("iter:{current_frame}", ", N=", N, ", K=", K), x = expression(x[1]), y = expression(x[2])) # gif画像を作成 gganimate::animate(trace_model_graph, nframes = MaxIter + 1, fps = 10)
・クラスタの散布図の推移
# クラスタの推移を作図 trace_cluster_graph <- ggplot() + geom_contour(data = model_df, aes(x = x_1, y = x_2, z = density), color = "red", alpha = 0.5, linetype = "dashed") + # 真の分布 geom_point(data = mu_df, aes(x = x_1, y = x_2), color = "red", shape = 4, size = 5) + # 真の平均 geom_contour_filled(data = trace_model_df, aes(x = x_1, y = x_2, z = density), alpha = 0.6) + # 期待値による分布 geom_point(data = trace_cluster_df, aes(x = x_n1, y = x_n2, color = cluster)) + # サンプルしたクラスタ gganimate::transition_manual(label) + # フレーム labs(title = "Gaussian Mixture Model:Collapsed Gibbs Sampling", subtitle = "{current_frame}", x = expression(x[1]), y = expression(x[2])) # gif画像を作成 gganimate::animate(trace_cluster_graph, nframes = MaxIter + 1, fps = 10)
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
ついに実装できたあぁ。長かったぁ、お疲れ様です、お付き合いありがとうございました。5章はやりません。この章でこの本は修了とします。いやまだこの後Pythonの方を更新するんだけど。
【次節の内容】
直接は続きませんが、5.4節のトピックモデルの内容はこちらで詳しく書いてます。
