はじめに
『パターン認識と機械学習』の独学時のまとめです。一連の記事は「数式の行間埋め」または「R・Pythonでの実装」からアルゴリズムの理解を補助することを目的としています。本とあわせて読んでください。
この記事は、9.2.2項の内容です。多次元混合ガウス分布(多変量混合正規分布)に対するEMアルゴリズムによる最尤推定をR言語で実装します。
【数式読解編】
【他の節一覧】
【この節の内容】
・Rでやってみる
人工データを用いて、EMアルゴリズムによる最尤推定を行ってみましょう。
利用するパッケージを読み込みます。
# 9.2.2項で利用するパッケージ library(tidyverse) library(mvnfast)
mvnfastは、多次元ガウス分布に関するパッケージです。多次元ガウス分布に従う乱数生成関数rmvn()と確率密度関数dmvn()を使います。
・真の分布の設定
まずは、データを生成する真の分布を設定します。この節では、$K$個の$D$次元ガウス分布からなる混合ガウス分布$\sum_{k=1}^K \pi_k \mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Lambda}_k^{-1})$とします。また、各データに対するクラスタの割り当てをカテゴリ分布$\mathrm{Cat}(\mathbf{z}_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_{k1}, \mu_{k2}, \cdots, \mu_{kD})$をmu_truth_kd、共分散行列パラメータ$\boldsymbol{\Sigma} = \{\boldsymbol{\Sigma}_1, \boldsymbol{\Sigma}_2, \cdots, \boldsymbol{\Sigma}_K\}$、$\boldsymbol{\Sigma}_k = (\sigma_{k11}^2, \cdots, \sigma_{kDD}^2)$をsigma2_truth_ddkとして、値を指定します。$\boldsymbol{\Sigma}_k$は正定値行列です。
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$となる($z_{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_{n1}, x_{n2})$がとり得る点(値の組み合わせ)を作成します。$x_{n1}$がとり得る値(x軸の値)をseq()で作成してx_1_vecとします。この例では、$K$個のクラスタそれぞれで平均値から標準偏差の3倍を引いた値と足した値を計算して、その最小値から最大値までを範囲とします。length.out引数を使うと指定した要素数で等間隔に切り分けます。by引数を使うと切り分ける間隔を指定できます。処理が重い場合は、この値を調整してください。同様に、$x_{n2}$がとり得る値(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_dens <- 0 for(k in 1:K) { # クラスタkの分布の確率密度を計算 tmp_dens <- mvnfast::dmvn( X = x_point_mat, mu = mu_truth_kd[k, ], sigma = sigma2_truth_ddk[, , k] ) # K個の分布を線形結合 model_dens <- model_dens + pi_truth_k[k] * tmp_dens }
$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_densに加算していきます。
計算結果をデータフレームに格納します。
# 真の分布をデータフレームに格納 model_df <- tibble( x_1 = x_point_mat[, 1], x_2 = x_point_mat[, 2], density = model_dens ) # 確認 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{Z} = \{\mathbf{z}_1, \mathbf{z}_2, \cdots, \mathbf{z}_N\}$、$\mathbf{z}_n = (z_{n1}, z_{n2}, \cdots, z_{nK})$を生成し、各データにクラスタを割り当てます。次に、割り当てられたクラスタに従い観測データ$\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{Z}$を生成します。
# (観測)データ数を指定 N <- 250 # 潜在変数を生成 z_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$次元ベクトルをz_truth_nkとします。ただし、クラスタが行でデータが列方向に並んで出力されるので、t()で転置しておきます。
結果は次のようになります。
# 確認 z_truth_nk[1:5, ]
## [,1] [,2] [,3] ## [1,] 1 0 0 ## [2,] 0 0 1 ## [3,] 0 1 0 ## [4,] 1 0 0 ## [5,] 0 0 1
これが本来は観測できない潜在変数であり、真のクラスタです。各行がデータを表し、値が1の列番号が割り当てられたクラスタ番号を表します。
z_truth_nkから各データのクラスタ番号を抽出します。
# クラスタ番号を取得 z_truth_n <- which(t(z_truth_nk) == 1, arr.ind = TRUE) %>% .[, "row"] # 確認 z_truth_n[1:5]
## [1] 1 3 2 1 3
which()を使って行(データ)ごとに要素が1の列番号を抽出します。
各データに割り当てられたクラスタに従い$N$個のデータを生成します。
# (観測)データを生成 x_nd <- matrix(0, nrow = N, ncol = D) for(n in 1:N) { # n番目のデータのクラスタ番号を取得 k <- z_truth_n[n] # n番目のデータを生成 x_nd[n, ] = mvnfast::rmvn(n = 1, mu = mu_truth_kd[k, ], sigma = sigma2_truth_ddk[, , k]) }
各データ$\mathbf{x}_n$に割り当てられたクラスタ($z_{nk} = 1$である$k$)のパラメータ$\boldsymbol{\mu}_k,\ \boldsymbol{\Sigma}_k$を持つ多次元ガウス分布に従いデータ(値)を生成します。多次元ガウス分布の乱数は、mvnfastパッケージのrmvn()で生成できます。
データごとにクラスタが異なるので、for()で1データずつ処理します。$n$番目のデータのクラスタ番号z_true_n[n]を取り出してkとします。そのクラスタに従ってデータ$\mathbf{x}_n$を生成します。平均の引数muにmu_true_kd[k, ]、共分散行列の引数sigmaにsigma2_true_ddk[, , k]を指定します。また1データずつ生成するので、データ数の引数nは1です。
観測データもグラフで確認しましょう。作図用のデータフレームを作成します。
# 観測データと真のクラスタ番号をデータフレームに格納 data_df <- tibble( x_n1 = x_nd[, 1], x_n2 = x_nd[, 2], cluster = as.factor(z_truth_n) ) # 確認 head(data_df)
## # A tibble: 6 x 3 ## x_n1 x_n2 cluster ## <dbl> <dbl> <fct> ## 1 -9.18 26.3 1 ## 2 26.9 -30.9 3 ## 3 -21.8 -30.4 2 ## 4 17.5 49.8 1 ## 5 -1.04 -17.0 3 ## 6 13.9 37.3 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 = data_df, aes(x = x_n1, y = x_n2, color = cluster)) + # 真のクラスタ labs(title = "Mixture of Gaussians", subtitle = paste0('N=', N, ", pi=(", paste0(pi_truth_k, collapse = ", "), ")"), x = expression(x[1]), y = expression(x[2]))
散布図はgeom_point()で描画します。クラスタごとに色分けしています。観測データから各データのクラスタも推定します。
・初期値の設定
パラメータ$\boldsymbol{\mu},\ \boldsymbol{\Lambda},\ \boldsymbol{\pi}$の初期値を設定します。
ランダムに値を生成してパラメータを初期化します。
# 平均パラメータの初期値を生成 mu_kd <- matrix(0, nrow = K, ncol = D) for(d in 1:D) { mu_kd[, d] <- runif(n = K, min = min(x_nd[, d]), max = max(x_nd[, d])) } # 共分散行列の初期値を指定 sigma2_ddk <- (diag(D) * 1000) %>% rep(times = K) %>% array(dim = c(D, D, K)) # 混合係数の初期値を生成 pi_k <- runif(n = K, min = 0, max = 1) pi_k <- pi_k / sum(pi_k)
$K$個の平均パラメータ$\boldsymbol{\mu}$をmu_kdとして値を生成します。この例では、次元ごとに観測データx_ndの最小値から最大値を範囲とする一様乱数を設定します。一様乱数はrunif()で生成できます。
$K$個の共分散行列$\boldsymbol{\mu}$をsigma2_ddkとして値を指定します。(ランダムに共分散行列を作成するのは少し手間なので)この例では、全てのクラスタを分散1000、共分散0とします。$D \times D$の単位行列をdiag(D)で作成して1000を掛けます。それをrep()で複製して、array()で3次元配列に整形します。
混合係数$\boldsymbol{\pi}$をpi_kとして値を生成します。この例では、0以上の一様乱数を$K$個生成して、総和で割ることで$\sum_{k=1}^K \pi_k = 1,\ \pi_k \geq 0$となるように正規化します。
パラメータmu_kd, sigma2_ddkとpi_kの初期値を用いて、混合分布を計算します。
# 初期値による混合分布を計算 init_dens <- 0 for(k in 1:K) { # クラスタkの分布の確率密度を計算 tmp_dens <- mvnfast::dmvn( X = x_point_mat, mu = mu_kd[k, ], sigma = sigma2_ddk[, , k] ) # K個の分布を線形結合 init_dens <- init_dens + pi_k[k] * tmp_dens } # 初期値による分布をデータフレームに格納 init_df <- tibble( x_1 = x_point_mat[, 1], x_2 = x_point_mat[, 2], density = init_dens ) # 初期値による分布を作図 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 = "Mixture of Gaussians", subtitle = paste0("iter:", 0, ", K=", K), x = expression(x[1]), y = expression(x[2]))
真の分布のときと同様に計算して作図します。更新を繰り返すことで真の分布に近付けていきます。
初期値による対数尤度を計算します。
# 初期値による対数尤度を計算:式(9.14) term_dens_nk <- matrix(0, nrow = N, ncol = K) for(k in 1:K) { # クラスタkの混合分布の確率密度を計算:式(9.14)の波括弧 tmp_N_dens_n <- mvnfast::dmvn( X = x_nd, mu = mu_kd[k, ], sigma = sigma2_ddk[, , k] ) term_dens_nk[, k] <- pi_k[k] * tmp_N_dens_n } L <- sum(log(rowSums(term_dens_nk)))
次の式で対数尤度を計算します。
計算結果は次のようになります。
# 確認
L
## [1] -2434.502
以上でモデルが整いました。次は、推論を行います。
・最尤推定
最尤推定では、EステップとMステップを交互に繰り返します。Eステップでは、負担率$\gamma(\mathbf{Z})$を更新します。Mステップでは、対数尤度が最大化されるようにパラメータ$\boldsymbol{\mu},\ \boldsymbol{\Sigma},\ \boldsymbol{\pi}$を更新します。
まずは、パラメータの初期値を用いて負担率を更新します。続いて、更新した負担率を用いて各パラメータを更新します。ここまでが1回の試行です。
次の試行では、前回更新したパラメータを用いて負担率を更新となります。この処理を指定した回数繰り返します。
・コード全体
負担率の計算時に利用するtmp_gamma_nkは、添字を使って代入するため予め変数を用意しておく必要があります。
また、後で各パラメータの更新値や分布の推移を確認するために、それぞれ試行ごとにtrace_***に値を保存していきます。関心のあるパラメータを記録してください。
# 試行回数を指定 MaxIter <- 100 # 計算に用いる中間変数を作成 tmp_gamma_nk <- matrix(0, nrow = N, ncol = K) # 推移の確認用の受け皿を作成 trace_L_i <- rep(0, MaxIter + 1) trace_gamma_ink <- array(0, dim = c(MaxIter + 1, N, K)) trace_mu_ikd <- array(0, dim = c(MaxIter + 1, K, D)) trace_sigma2_iddk <- array(0, dim = c(MaxIter + 1, D, D, K)) trace_pi_ik <- matrix(0, nrow = MaxIter + 1, ncol = K) # 初期値を記録 trace_L_i[1] <- L trace_gamma_ink[1, , ] <- NA trace_mu_ikd[1, , ] <- mu_kd trace_sigma2_iddk[1, , , ] <- sigma2_ddk trace_pi_ik[1, ] <- pi_k # 最尤推定 for(i in 1:MaxIter) { # 負担率を計算:式(9.13) for(k in 1:K) { # 正規化前の負担率を計算:式(9.13)の分子 tmp_N_dens <- mvnfast::dmvn( X = x_nd, mu = mu_kd[k, ], sigma = sigma2_ddk[, , k] ) tmp_gamma_nk[, k] <- pi_k[k] * tmp_N_dens } gamma_nk <- tmp_gamma_nk / rowSums(tmp_gamma_nk) # 正規化 # 各クラスタのデータ数の期待値を計算:(9.18) N_k <- colSums(gamma_nk) for(k in 1:K) { # 平均パラメータの最尤解を計算:式(9.17) mu_kd[k, ] <- (t(gamma_nk[, k]) %*% x_nd / N_k[k]) %>% as.vector() # 共分散行列の最尤解を計算:(9.19) term_x_dn <- t(x_nd) - mu_kd[k, ] sigma2_ddk[, , k] <- term_x_dn %*% (gamma_nk[, k] * t(term_x_dn)) / N_k[k] } # 混合係数の最尤解を計算:式(9.22) pi_k <- N_k / N # 対数尤度を計算:式(9.14) for(k in 1:K) { # クラスタkの混合分布の確率密度を計算:式(9.14)の波括弧 tmp_N_dens_n <- mvnfast::dmvn( X = x_nd, mu = mu_kd[k, ], sigma = sigma2_ddk[, , k] ) term_dens_nk[, k] <- pi_k[k] * tmp_N_dens_n } L <- sum(log(rowSums(term_dens_nk))) # i回目の結果を記録 trace_L_i[i + 1] <- L trace_gamma_ink[i + 1, , ] <- gamma_nk trace_mu_ikd[i + 1, , ] <- mu_kd trace_sigma2_iddk[i + 1, , , ] <- sigma2_ddk trace_pi_ik[i + 1, ] <- pi_k # 動作確認 #print(paste0(i, ' (', round(i / MaxIter * 100, 1), '%)')) }
・処理の解説
続いて、EMアルゴリズムで行う処理を確認していきます。
・Eステップ
Eステップでは、負担率$\gamma(\mathbf{Z})$を更新(最尤解を計算)します。
前回のパラメータを用いて、負担率を計算します。
# 負担率を計算:式(9.13) for(k in 1:K) { # 正規化前の負担率を計算:式(9.13)の分子 tmp_N_dens <- mvnfast::dmvn( X = x_nd, mu = mu_kd[k, ], sigma = sigma2_ddk[, , k] ) tmp_gamma_nk[, k] <- pi_k[k] * tmp_N_dens } gamma_nk <- tmp_gamma_nk / rowSums(tmp_gamma_nk) # 正規化
次の式で負担率を計算します。
処理上は、全てのクラスタで分子の計算をしてtmp_gamma_nkとします。tmp_gamma_nkを行(データ)ごとに和をとったもので割る(正規化する)ことで、式(9.13)の計算になります。
計算結果は次のようになります。
# 確認 round(gamma_nk[1:5, ], 5) rowSums(gamma_nk[1:5, ])
## [,1] [,2] [,3] ## [1,] 0.00610 0.98909 0.00481 ## [2,] 0.00024 0.00005 0.99971 ## [3,] 0.99438 0.00071 0.00491 ## [4,] 0.00000 0.99899 0.00101 ## [5,] 0.69809 0.01484 0.28708 ## [1] 1 1 1 1 1
行(データ)ごとの和が1になります。
各クラスタが割り当てられたデータ数の期待値を計算します。
# 各クラスタのデータ数の期待値を計算:(9.18) N_k <- colSums(gamma_nk)
次の式を計算してN_kとします。
これは次のパラメータの更新に使います。
・Mステップ
Mステップでは、対数尤度を最大にするパラメータ$\boldsymbol{\pi},\ \boldsymbol{\Sigma},\ \boldsymbol{\pi}$を更新(最尤解を計算)します。
平均パラメータを更新します。
# 平均パラメータの最尤解を計算:式(9.17) mu_kd[k, ] <- (t(gamma_nk[, k]) %*% x_nd / N_k[k]) %>% as.vector()
$\boldsymbol{\mu}_k$の最尤解をクラスタごとに次の式で更新します。
# 確認
mu_kd
## [,1] [,2] ## [1,] -20.959732 -9.624021 ## [2,] 6.669279 36.245674 ## [3,] 29.332196 -18.487717
共分散行列を更新します。
# 共分散行列の最尤解を計算:(9.19) term_x_dn <- t(x_nd) - mu_kd[k, ] sigma2_ddk[, , k] <- term_x_dn %*% (gamma_nk[, k] * t(term_x_dn)) / N_k[k]
$\boldsymbol{\Sigma}_k$の最尤解をクラスタごとに次の式で更新します。
ただし、$\sum_{n=1}^N$の計算を効率よく処理するために、$\tilde{\mathbf{X}} = \mathbf{x}_n - \boldsymbol{\mu}_k$として、$\tilde{\mathbf{X}}^{\top} \tilde{\mathbf{X}}$で計算しています。
# 確認
sigma2_ddk
## , , 1 ## ## [,1] [,2] ## [1,] 128.33657 -44.56783 ## [2,] -44.56783 168.79502 ## ## , , 2 ## ## [,1] [,2] ## [1,] 287.43758 57.50498 ## [2,] 57.50498 271.85665 ## ## , , 3 ## ## [,1] [,2] ## [1,] 226.17242 -34.17174 ## [2,] -34.17174 329.27490
混合係数を更新します。
# 混合係数の最尤解を計算:式(9.22) pi_k <- N_k / N
$\boldsymbol{\pi}$の最尤解を次の式で計算します。
# 確認 pi_k sum(pi_k)
## [1] 0.2509721 0.4184309 0.3305970 ## [1] 1
これでパラメータ$\boldsymbol{\mu},\ \boldsymbol{\Sigma},\ \boldsymbol{\pi}$を更新できました。更新したパラメータを用いて、対数尤度を計算して記録しておきます。計算方法は「初期値の設定」のときと同じです。
以上がEMアルゴリズムで行う個々の処理です。これを繰り返し行うことで、徐々にパラメータを真の値に近付けていきます。
・推論結果の確認
推論結果を確認していきます。
・最後のパラメータの確認
パラメータmu_kd, sigma2_ddkとpi_kの最後の更新値を用いて、混合分布を計算します。
# 最後の更新値による混合分布を計算 res_dens <- 0 for(k in 1:K) { # クラスタkの分布の確率密度を計算 tmp_dens <- mvnfast::dmvn( X = x_point_mat, mu = mu_kd[k, ], sigma = sigma2_ddk[, , k] ) # K個の分布を線形結合 res_dens <- res_dens + pi_k[k] * tmp_dens } # 最終的な分布をデータフレームに格納 res_df <- tibble( x_1 = x_point_mat[, 1], x_2 = x_point_mat[, 2], density = res_dens ) # 真の平均をデータフレームに格納 mu_df <- tibble( x_1 = mu_truth_kd[, 1], x_2 = mu_truth_kd[, 2] ) # 最終的な分布を作図 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..)) + # 推定値による分布 labs(title = "Mixture of Gaussians:Maximum Likelihood", subtitle = paste0("iter:", MaxIter, ", L=", round(L, 1), ", N=", N, ", pi=(", paste0(round(pi_k, 3), collapse = ", "), ")"), x = expression(x[1]), y = expression(x[2]))
「初期値の設定」のときと同様にして作図できます。ここではさらに、目安として真の平均値をプロットしています。
各データに割り当てられたクラスタを散布図で確認します。
# 観測データと負担率が最大のクラスタ番号をデータフレームに格納 gamma_df <- tibble( x_n1 = x_nd[, 1], x_n2 = x_nd[, 2], cluster = as.factor(max.col(gamma_nk)), # 負担率が最大のクラスタ番号 prob = gamma_nk[cbind(1:N, max.col(gamma_nk))] # 負担率の最大値 )
各データにおいて負担率(確率)が最大のクラスタをそのデータのクラスタとします。max.col()を使って、gamma_nkの行(データ)ごとに値が最大の列番号(クラスタ番号)を抽出します。これをcluster列とします。また、その値をgamma_nk[cbind(1:N, max.col(gamma_nk))]で取り出してprob列とします。添字にマトリクスを指定することで抽出できます。
作成したデータフレームは次のようになります。
# 確認 head(gamma_df)
## # A tibble: 6 x 4 ## x_n1 x_n2 cluster prob ## <dbl> <dbl> <fct> <dbl> ## 1 -9.18 26.3 2 0.989 ## 2 26.9 -30.9 3 1.00 ## 3 -21.8 -30.4 1 0.994 ## 4 17.5 49.8 2 0.999 ## 5 -1.04 -17.0 1 0.698 ## 6 13.9 37.3 2 0.994
x_n1列とx_n2列が観測データのx軸とy軸の値です。そのデータx_nd[n, ]に割り当てられた(とみなした)クラスタ番号がcluster列で、そのクラスタとなる確率(事後確率)がprob列です。
クラスタの散布図を作成します。
# 最終的なクラスタを作図 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 = gamma_df, aes(x = x_n1, y = x_n2, color = cluster), alpha = gamma_df[["prob"]]) + # 負担率によるクラスタ labs(title = "Mixture of Gaussians:Maximum Likelihood", subtitle = paste0("iter:", MaxIter, ", L=", round(L, 1), ", N=", N, ", pi=(", paste0(round(pi_k, 3), collapse = ", "), ")"), x = expression(x[1]), y = expression(x[2]))
データの生成のときと同様にして作図します。ただしここでは、geom_point()の透過度の引数alphaにgamma_df[["prob"]]を指定することで、各データに割り当てられたクラスタとなる確率に従って濃淡を付けています。
対数尤度の推移を作図します。
# 対数尤度の推移をデータフレームに格納 L_df <- tibble( iteration = 0:MaxIter, L = trace_L_i ) # 対数尤度の推移を作図 ggplot(L_df, aes(x = iteration, y = L)) + geom_line() + labs(title = "Maximum Likelihood", subtitle = "Log Likelihood")
更新を繰り返すごとに対数尤度が増加しているのを確認できます。
・パラメータの推移の確認
パラメータの更新値の推移を確認します。
・コード(クリックで展開)
学習ごとの値を格納した各パラメータの配列trace_***をtidyr::pivot_longer()で縦長のデータフレームに変換して作図します。
# 平均パラメータの推移を作図 dplyr::as_tibble(trace_mu_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() + geom_hline(yintercept = as.vector(mu_truth_kd), color = "red", linetype = "dashed") + # 真の値 labs(title = "Maximum Likelihood", subtitle = expression(mu))
# 共分散行列の推移を作図 dplyr::as_tibble(trace_sigma2_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) + geom_hline(yintercept = as.vector(sigma2_truth_ddk), color = "red", linetype = "dashed") + # 真の値 labs(title = "Maximum Likelihood", subtitle = expression(Sigma))
# 混合係数の推移を作図 dplyr::as_tibble(trace_pi_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() + # 推定値 geom_hline(yintercept = pi_truth_k, color = "red", linetype = "dashed") + # 真の値 labs(title = "Maximum Likelihood", subtitle = expression(pi))
更新を繰り返すに従って真の値に近付いているのを確認できます。
・おまけ：アニメーションによる推移の確認
gganimateパッケージを利用して、各分布の推移のアニメーション(gif画像)を作成するためのコードです。
・コード(クリックで展開)
# 追加パッケージ library(gganimate)
各パラメータのを用いた分布とクラスタの推移を確認します。作図用のデータフレームを作成します。
# 作図用のデータフレームを作成 trace_model_df <- tibble() trace_cluster_df <- tibble() for(i in 1:(MaxIter + 1)) { # i回目の混合分布を計算 res_dens <- 0 for(k in 1:K) { # クラスタkの分布の確率密度を計算 tmp_dens <- mvnfast::dmvn( X = x_point_mat, mu = trace_mu_ikd[i, k, ], sigma = trace_sigma2_iddk[i, , , k] ) # K個の分布の加重平均を計算 res_dens <- res_dens + trace_pi_ik[i, k] * tmp_dens } # i回目の分布をデータフレームに格納 res_df <- tibble( x_1 = x_point_mat[, 1], x_2 = x_point_mat[, 2], density = res_dens, label = paste0( "iter:", i - 1, ", L=", round(trace_L_i[i], 1), ", N=", N, ", pi=(", paste0(round(trace_pi_ik[i, ], 3), collapse = ", "), ")" ) %>% as.factor() ) # 結果を結合 trace_model_df <- rbind(trace_model_df, res_df) # 観測データとi回目のクラスタ番号をデータフレームに格納 tmp_gamma_nk <- trace_gamma_ink[i, , ] # i回目の結果 gamma_df <- tibble( x_n1 = x_nd[, 1], x_n2 = x_nd[, 2], cluster = as.factor(max.col(tmp_gamma_nk)), # 負担率が最大のクラスタ番号 prob = tmp_gamma_nk[cbind(1:N, max.col(tmp_gamma_nk))], # 負担率の最大値 label = paste0( "iter:", i - 1, ", L=", round(trace_L_i[i], 1), ", N=", N, ", pi=(", paste0(round(trace_pi_ik[i, ], 3), collapse = ", "), ")" ) %>% as.factor() ) # 結果を結合 trace_cluster_df <- rbind(trace_cluster_df, gamma_df) # 動作確認 #print(paste0(i - 1, ' (', round((i - 1) / MaxIter * 100, 1), '%)')) }
各試行における計算結果を同じデータフレームに格納していく必要があります。処理に時間がかかる場合は、描画するフレーム数(1:(MaxIter + 1))や、x_point_matの範囲や間隔を調整してください。
作成したデータフレームは次のようになります。
# 確認 head(trace_model_df) head(trace_cluster_df)
## # A tibble: 6 x 4 ## x_1 x_2 density label ## <dbl> <dbl> <dbl> <fct> ## 1 -53.5 -67.4 0.00000637 iter:0, L=-2434.5, N=250, pi=(0.66, 0.059, 0.281) ## 2 -53.1 -67.4 0.00000647 iter:0, L=-2434.5, N=250, pi=(0.66, 0.059, 0.281) ## 3 -52.7 -67.4 0.00000657 iter:0, L=-2434.5, N=250, pi=(0.66, 0.059, 0.281) ## 4 -52.3 -67.4 0.00000668 iter:0, L=-2434.5, N=250, pi=(0.66, 0.059, 0.281) ## 5 -51.8 -67.4 0.00000679 iter:0, L=-2434.5, N=250, pi=(0.66, 0.059, 0.281) ## 6 -51.4 -67.4 0.00000689 iter:0, L=-2434.5, N=250, pi=(0.66, 0.059, 0.281) ## # A tibble: 6 x 5 ## x_n1 x_n2 cluster prob label ## <dbl> <dbl> <fct> <dbl> <fct> ## 1 -9.18 26.3 <NA> NA iter:0, L=-2434.5, N=250, pi=(0.66, 0.059, 0.281) ## 2 26.9 -30.9 <NA> NA iter:0, L=-2434.5, N=250, pi=(0.66, 0.059, 0.281) ## 3 -21.8 -30.4 <NA> NA iter:0, L=-2434.5, N=250, pi=(0.66, 0.059, 0.281) ## 4 17.5 49.8 <NA> NA iter:0, L=-2434.5, N=250, pi=(0.66, 0.059, 0.281) ## 5 -1.04 -17.0 <NA> NA iter:0, L=-2434.5, N=250, pi=(0.66, 0.059, 0.281) ## 6 13.9 37.3 <NA> NA iter:0, L=-2434.5, N=250, pi=(0.66, 0.059, 0.281)
分布の推移をgif画像で出力します。
# 分布の推移を作図 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(label) + # フレーム labs(title = "Mixture of Gaussians:Maximum Likelihood", subtitle = paste0("{current_frame}"), x = expression(x[1]), y = expression(x[2])) # gif画像を作成 gganimate::animate(trace_model_graph, nframes = MaxIter + 1, fps = 10)
各フレームの順番を示す列をgganimate::transition_manual()に指定します。初期値を含むため、フレーム数はMaxIter + 1です。
クラスタの推移をgif画像で出力します。
# クラスタの推移を作図 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 = "Mixture of Gaussians:Maximum Likelihood", subtitle = "{current_frame}", x = expression(x[1]), y = expression(x[2])) # gif画像を作成 gganimate::animate(trace_cluster_graph, nframes = MaxIter + 1, fps = 10)
参考文献
- C.M.ビショップ著,元田 浩・他訳『パターン認識と機械学習 上下』,丸善出版,2012年.
おわりに
やってることベイズ推論入門のときと大差ないな。
【関連する内容】
