からっぽのしょこ

読んだら書く!書いたら読む!同じ事は二度調べ(たく)ない

【R】4.4.2:ガウス混合モデルのギブスサンプリング【緑ベイズ入門のノート】

はじめに

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

 この記事は4.4.2項の内容です。ガウス混合モデルにおけるギブスサンプリングによる近似推論をRで実装します。

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

【数式読解編】

www.anarchive-beta.com

【前節の内容】

www.anarchive-beta.com

【他の節一覧】

www.anarchive-beta.com

【この節の内容】

・Rでやってみる

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

 利用するパッケージを読み込みます。

# 利用パッケージ
library(tidyverse)
library(mvnfast)
library(MCMCpack)

 mvnfastは、多次元ガウス分布に関するパッケージです。多次元ガウス分布に従う乱数の生成関数rmvn()、確率密度の計算関数dmvn()を使います。またディリクレ分布に従う乱数生成には、MCMCpackパッケージのrdirichlet()を使います。

・真の観測モデルの設定

 観測モデルのパラメータを設定します。作図に関しては、2次元のグラフで表現するため$D = 2$のときのみ動作します。

# 真の観測モデルのパラメータを指定
D <- 2 # (固定)
K <- 3
mu_true_kd <- matrix(
  c(0, 4, 
    -5, -5, 
    5, -2.5), nrow = K, ncol = D, byrow = TRUE
)
sigma2_true_ddk <- array(
  c(8, 0, 0, 8, 
    4, -2.5, -2.5, 4, 
    6.5, 4, 4, 6.5), dim = c(D, D, K)
)

# 確認
mu_true_kd
##      [,1] [,2]
## [1,]    0  4.0
## [2,]   -5 -5.0
## [3,]    5 -2.5
sigma2_true_ddk
## , , 1
## 
##      [,1] [,2]
## [1,]    8    0
## [2,]    0    8
## 
## , , 2
## 
##      [,1] [,2]
## [1,]  4.0 -2.5
## [2,] -2.5  4.0
## 
## , , 3
## 
##      [,1] [,2]
## [1,]  6.5  4.0
## [2,]  4.0  6.5

 真の観測モデルにおけるクラスタごとの平均パラメータ$\boldsymbol{\mu} = \{\boldsymbol{\mu_1}, \cdots, \boldsymbol{\mu}_K\}$、$\boldsymbol{\mu}_k = \{\mu_{k,1}, \cdots, \mu_{k,D}\}$を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_ddkとして、値を指定します。ただし$\boldsymbol{\Sigma}_k$は、正定値行列である必要があります。またmu_true_kdについては、値を指定しやすくするため引数にbyrow = TRUEを指定しておきます。
 この2つのパラメータの値を観測データから求めることが目的となります。

 続いて混合比率を設定します。

# 真の混合比率を指定
pi_true_k <- c(0.5, 0.2, 0.3)

 混合比率$\boldsymbol{\pi} = \{\pi_1, \cdots, \pi_K\}$をpi_true_kとして、値を指定します。ここで$\pi_k$は各データ$\mathbf{x}_n$のクラスタ(各潜在変数$\mathbf{s}_n$)が$k$となる確率であり、$0 \leq \pi_k \leq 1,\ \sum_{k=1}^K \pi_k = 1$の値をとります。このパラメータの値も求めます。

 以上が観測モデルの設定です。続いて、設定したガウス混合モデルに従ってデータを生成します。まずは各データが持つクラスタを割り当てます。

# 観測データの真のクラスタを生成
N <- 250
s_true_nk <- rmultinom(n = N, size = 1, prob = pi_true_k) %>% 
  t()

# 確認
s_true_nk[1:5, ]
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    1    0    0
## [4,]    0    0    1
## [5,]    1    0    0

 $\boldsymbol{\pi}$をパラメータとするカテゴリ分布に従い、潜在変数$\mathbf{S} = \{s_{1,1}, \cdots, s_{N,K}\}$を生成します。カテゴリ分布の乱数は、多項分布の乱数生成関数rmultinom()の引数size1を指定することで生成できます。確率値引数probpi_true_k、データ数引数nNを指定します。
 転置した出力は、各列がクラスタ1から$K$に対応していて、各行(各データ)に割り当てられたクラスタが1でそれ以外は0となります。これが(本来は観測できない)真のクラスタとなります。

 処理の都合上、割り当てられたクラスタを抽出しておきます。

# 各データのクラスタを抽出
res_s <- which(t(s_true_nk) == 1, arr.ind = TRUE)
s_true_n <- res_s[, "row"]

# 確認
s_true_n[1:5]
## [1] 1 2 1 3 1

 s_true_nkは、行がデータ、列がクラスタに対応しています。which()を使って行(データ)ごとに要素が1のインデックスを取得します。

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

# 観測データを生成
x_nd <- matrix(0, nrow = N, ncol = D)
for(n in 1:N) {
  k <- s_true_n[n] # クラスタを取得
  x_nd[n, ] = mvnfast::rmvn(n = 1, mu = mu_true_kd[k, ], sigma = sigma2_true_ddk[, , k])
}

# 確認
x_nd[1:5, ]
##            [,1]      [,2]
## [1,] -0.9709953  1.546786
## [2,] -6.3991780 -2.219533
## [3,]  3.6932798  2.003394
## [4,]  1.9321825 -5.118773
## [5,] -2.1357873  2.247528

 各クラスタのパラメータ$\boldsymbol{\mu}_k,\ \boldsymbol{\Sigma}_k$を持つ多次元ガウス分布に従い乱数を生成します。多次元ガウス分布の乱数は、mvnfast::rmvn()で生成できます。
 データごとにクラスタが異なるので、for()で1データずつs_true_nからクラスタの値を取り出し、そのクラスタに従ってデータ$\mathbf{x}_n$を生成します。平均引数mumu_true_kd[k, ]、分散共分散行列引数sigmasigma2_true_ddk[, , k]を指定します。また1データずつ生成するので、データ数引数n1です。

 グラフを作成して観測モデルと観測データを確認しましょう。作図用のデータフレームを作成します。

# 作図用の点を生成
x_line <- seq(-10, 10, by = 0.1)
point_df <- tibble(
  x1 = rep(x_line, times = length(x_line)), 
  x2 = rep(x_line, each = length(x_line))
)

# 作図用のデータフレームを作成
model_true_df <- tibble()
sample_df <- tibble()
for(k in 1:K) {
  # 真の観測モデルを計算
  tmp_model_df <- cbind(
    point_df, 
    density = mvnfast::dmvn(
      X = as.matrix(point_df), mu = mu_true_kd[k, ], sigma = sigma2_true_ddk[, , k]
    ), 
    cluster = as.factor(k)
  )
  model_true_df <- rbind(model_true_df, tmp_model_df)
  
  # 観測データのデータフレーム
  k_idx <- which(s_true_n == k)
  tmp_sample_df <- tibble(
    x1 = x_nd[k_idx, 1], 
    x2 = x_nd[k_idx, 2], 
    cluster = as.factor(k)
  )
  sample_df <- rbind(sample_df, tmp_sample_df)
}

# 確認
head(model_true_df)
##      x1  x2      density cluster
## 1 -10.0 -10 1.837732e-10       1
## 2  -9.9 -10 2.081122e-10       1
## 3  -9.8 -10 2.353803e-10       1
## 4  -9.7 -10 2.658886e-10       1
## 5  -9.6 -10 2.999760e-10       1
## 6  -9.5 -10 3.380107e-10       1
head(sample_df)
## # A tibble: 6 x 3
##       x1    x2 cluster
##    <dbl> <dbl> <fct>  
## 1 -0.971  1.55 1      
## 2  3.69   2.00 1      
## 3 -2.14   2.25 1      
## 4 -4.02   4.99 1      
## 5  4.20   9.33 1      
## 6 -0.526  5.14 1

 等高線グラフ用に格子状の点を用意する必要があります。seq()で描画範囲と間隔を決めて、x軸の点($x_{n,1}$がとり得る値)とy軸の点($x_{n,2}$がとり得る値)が直交する点(格子状の点)を作成します。
 各交点$(x_{n,1}, x_{n,2})$の確率密度のマトリクスをmvnfast::dmvn()の第1引数に渡すことで計算します。他の引数についてはmvnfast::rmvn()と同じです。これをクラスタごとに行います。
 観測データについても作図用にデータフレームを作成します。which(s_true_n == k)で各クラスタが割り当てられたデータのインデックスを取得しておき、そのインデックスを添字として用いてx_ndから各クラスタのデータを取り出します。

 観測モデル(多次元ガウス分布の確率密度)の等高線図と、観測データの散布図を重ねて作図します。

# 真の観測モデルを作図
ggplot() + 
  geom_contour(data = model_true_df, aes(x1, x2, z = density, color = cluster)) + # 真の観測モデル
  geom_point(data = sample_df, aes(x1, x2, color = cluster)) + # 真の観測データ
  labs(title = "Gaussian Mixture Model", 
       subtitle = paste0('K=', K, ', N=', N), 
       x = expression(x[1]), y = expression(x[2]))

f:id:anemptyarchive:20201130190357p:plain
ガウス混合モデル

 等高線グラフはgeom_contour()、散布図はgeom_point()で描画します。

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

・パラメータの設定

 観測モデルのパラメータは本来知り得ないものです。そこで初期値を設定し、ギブスサンプリングにより真の値に近付けていきます。

# 観測モデルのパラメータの初期値を指定
mu_kd <- matrix(sample(seq(-10, 10, 0.1), size = K * D), nrow = K, ncol = D)
lambda_ddk <- array(
  c(0.01, 0, 0, 0.01, 
    0.01, 0, 0, 0.01, 
    0.01, 0, 0, 0.01), dim = c(D, D, K)
)

# 混合比率の初期値を指定
pi_k <- sample(seq(0, 1, by = 0.01), size = K)
pi_k <- pi_k / sum(pi_k) # 正規化

 $\boldsymbol{\mu},\ \boldsymbol{\Lambda},\ \boldsymbol{\pi}$をそれぞれmu_kdlambda_ddkpi_kとします。

 各事前分布のパラメータを設定します。

# 事前分布のパラメータを指定
beta <- 1
m_d <- rep(0, D)
nu <- D
w_dd <- diag(D) * 10
alpha_k <- rep(1, K)

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

 以上で推論に必要なモデルに関する設定は完了です。次はギブスサンプリングを実装します。

・ギブスサンプリング

 潜在変数(各データのクラスタ)$\mathbf{S}$とパラメータ$\boldsymbol{\mu},\ \boldsymbol{\Lambda},\ \boldsymbol{\pi}$のサンプリングを交互に行います。またサンプルされたパラメータを用いて、繰り返し事後分布のパラメータを計算(事前分布のパラメータを更新)します。

・コード全体

 本では縦ベクトルとしているところを、この例ではベクトルとして扱うため、転置などの処理が異なっている点に注意してください。

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

# 試行回数を指定
MaxIter <- 500

# 推移の確認用の受け皿
trace_s_in <- matrix(0, nrow = MaxIter, ncol = N)
trace_mu_ikd <- array(0, dim = c(MaxIter+1, K, D))
trace_lambda_iddk <- array(0, dim = c(MaxIter+1, D, D, K))
trace_mu_ikd[1, , ] <- mu_kd
trace_lambda_iddk[1, , , ] <- lambda_ddk

# ギブスサンプリング
for(i in 1:MaxIter) {
  
  # 初期化
  eta_nk <- matrix(0, nrow = N, ncol = K)
  s_nk <- matrix(0, nrow = N, ncol = K)
  beta_hat_k <- rep(0, K)
  m_hat_kd <- matrix(0, nrow = K, ncol = D)
  nu_hat_k <- rep(0, K)
  w_hat_ddk <- array(0, dim = c(D, D, K))
  alpha_hat_k <- rep(0, K)
  
  # 潜在変数変数のパラメータを計算:式(4.94)
  for(k in 1:K) {
    tmp_term_dn <- t(x_nd) - mu_kd[k, ]
    tmp_eta_n <- diag(
      t(tmp_term_dn) %*% lambda_ddk[, , k] %*% tmp_term_dn
    )
    tmp_eta <- 0.5 * log(det(lambda_ddk[, , k]) + 1e-7) + log(pi_k[k] + 1e-7)
    eta_nk[, k] <- exp(-0.5 * tmp_eta_n + tmp_eta)
  }
  eta_nk <- eta_nk / rowSums(eta_nk) # 正規化
  
  # 潜在変数をサンプル
  for(n in 1:N) {
    s_nk[n, ] <- rmultinom(n = 1, size = 1, prob = eta_nk[n, ]) %>% 
      as.vector()
  }
  
  # 観測モデルのパラメータをサンプル
  for(k in 1:K) {
    
    # muの事後分布のパラメータを計算:式(4.99)
    beta_hat_k[k] <- sum(s_nk[, k]) + beta
    m_hat_kd[k, ] <- (colSums(s_nk[, k] * x_nd) + beta * m_d) / beta_hat_k[k]
    
    # lambdaの事後分布のパラメータを計算:式(4.103)
    nu_hat_k[k] <- sum(s_nk[, k]) + nu
    tmp_w1_dd <- t(s_nk[, k] * x_nd) %*% x_nd
    tmp_w2_dd <- beta * matrix(m_d) %*% t(m_d)
    tmp_w3_dd <- beta_hat_k[k] * matrix(m_hat_kd[k, ]) %*% t(m_hat_kd[k, ])
    w_hat_ddk[, , k] <- solve(
      tmp_w1_dd + tmp_w2_dd - tmp_w3_dd + solve(w_dd)
    )
    
    # lambdaをサンプル:式(4.102)
    lambda_ddk[, , k] <- rWishart(n = 1, df = nu_hat_k[k], Sigma = w_hat_ddk[, , k])
    
    # muをサンプル:式(4.98)
    mu_kd[k, ] <- mvnfast::rmvn(
      n = 1, mu = m_hat_kd[k, ], sigma = solve(beta_hat_k[k] * lambda_ddk[, , k])
    ) %>% 
      as.vector()
  }
  
  # 混合比率のパラメータを計算:式(4.45)
  alpha_hat_k <- colSums(s_nk) + alpha_k
  
  # piをサンプル:式(4.44)
  pi_k <- MCMCpack::rdirichlet(n = 1, alpha = alpha_hat_k) %>% 
    as.vector()
  
  # 値を記録
  res_s <- which(t(s_nk) == 1, arr.ind = TRUE)
  trace_s_in[i, ] <- res_s[, "row"]
  trace_mu_ikd[i+1, , ] <- mu_kd
  trace_lambda_iddk[i+1, , , ] <- lambda_ddk
  
  # 動作確認
  print(paste0(i, ' (', round(i / MaxIter * 100, 1), '%)'))
}

 trace_***は、パラメータや分布の推移を確認するためのオブジェクトです。推論自体には不要です。


・処理の解説

 始めに、パラメータ$\boldsymbol{\eta}_n = \{\eta_{n,1}, \cdots, \eta_{n,K}\}$を持つカテゴリ分布に従い各潜在変数$\mathbf{s}_n$をサンプルします(各データにクラスタを割り当てます)。

 まずは各データ$\mathbf{x}_n$のクラスタが$k$である確率を表すパラメータ$\eta_{n,k}$を計算します。

# 初期化
eta_nk <- matrix(0, nrow = N, ncol = K)

# 潜在変数変数のパラメータを計算:式(4.94)
for(k in 1:K) {
  tmp_term_dn <- t(x_nd) - mu_kd[k, ]
  tmp_eta_n <- diag(
    t(tmp_term_dn) %*% lambda_ddk[, , k] %*% tmp_term_dn
  )
  tmp_eta <- 0.5 * log(det(lambda_ddk[, , k]) + 1e-7) + log(pi_k[k] + 1e-7)
  eta_nk[, k] <- exp(-0.5 * tmp_eta_n + tmp_eta)
}

# 確認
round(eta_nk[1:5, ], 4)
##        [,1]   [,2]   [,3]
## [1,] 0.0000 0.0391 0.0000
## [2,] 0.0000 0.0006 0.0191
## [3,] 0.0000 0.0265 0.0000
## [4,] 0.0471 0.0020 0.0000
## [5,] 0.0000 0.0331 0.0000

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

$$ \eta_{n,k} \propto \exp \left\{ - \frac{1}{2} (\mathbf{x}_n - \boldsymbol{\mu}_k)^{\top} \boldsymbol{\Lambda}_k (\mathbf{x}_n - \boldsymbol{\mu}_k) + \frac{1}{2} \ln |\boldsymbol{\Lambda}_k| + \ln \pi_k \right\} \tag{4.94} $$

 ただし全てのデータ($n = 1, \cdots, N$)を一度に処理するために、$(\mathbf{x}_n - \boldsymbol{\mu}_k)^{\top} \boldsymbol{\Lambda}_k (\mathbf{x}_n - \boldsymbol{\mu}_k)$の計算を1から$N$の全ての組み合わせで行います(for()内の3行目)。この計算結果は、$N \times N$のマトリクスになります。これは例えば、1行$N$列目の要素は$(\mathbf{x}_1 - \boldsymbol{\mu}_k)^{\top} \boldsymbol{\Lambda}_k (\mathbf{x}_N - \boldsymbol{\mu}_k)$の計算結果です。これは不要なので、diag()で対角成分(同じデータによる計算結果)のみ取り出します。
 またこの計算において$\ln 0$とならないように、微小な値1e-7($10^{-7} = 0.0000001$)を加えています。

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

# 正規化
eta_nk <- eta_nk / rowSums(eta_nk)

# 確認
round(eta_nk[1:5, ], 4)
##        [,1]   [,2]   [,3]
## [1,] 0.0000 0.9998 0.0002
## [2,] 0.0000 0.0288 0.9712
## [3,] 0.0008 0.9992 0.0000
## [4,] 0.9590 0.0408 0.0002
## [5,] 0.0000 0.9995 0.0005
rowSums(eta_nk[1:5, ])
## [1] 1 1 1 1 1

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

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


 パラメータが得られたので、カテゴリ分布に従い潜在変数$\mathbf{S}$を生成します。

# 初期化
s_nk <- matrix(0, nrow = N, ncol = K)

# 潜在変数をサンプル
for(n in 1:N) {
  s_nk[n, ] <- rmultinom(n = 1, size = 1, prob = eta_nk[n, ]) %>% 
    as.vector()
}

# 確認
s_nk[1:5, ]
##      [,1] [,2] [,3]
## [1,]    0    1    0
## [2,]    0    0    1
## [3,]    0    1    0
## [4,]    1    0    0
## [5,]    0    1    0

 真のクラスタのときと同様に、潜在変数をサンプルします。ただしデータによってどのクラスタが割り当てられるのかの確率が異なるため、for()で1データずつ生成します。

 続いて、それぞれ仮定した分布に従いパラメータ$\boldsymbol{\mu},\ \boldsymbol{\Lambda},\ \boldsymbol{\pi}$をサンプルします。ここからの内容は、for()でクラスタごとに行います。

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

# クラスタを指定
k <- 1

# 初期化
beta_hat_k <- rep(0, K)
m_hat_kd <- matrix(0, nrow = K, ncol = D)

# muの事後分布のパラメータを計算:式(4.99)
beta_hat_k[k] <- sum(s_nk[, k]) + beta
m_hat_kd[k, ] <- (colSums(s_nk[, k] * x_nd) + beta * m_d) / beta_hat_k[k]

# 確認
beta_hat_k
## [1] 69  0  0
m_hat_kd
##          [,1]      [,2]
## [1,] 3.891653 -3.188316
## [2,] 0.000000  0.000000
## [3,] 0.000000  0.000000

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

$$ \begin{align} \hat{\beta}_k &= \sum_{n=1}^N s_{n,k} + \beta \\ \hat{\mathbf{m}}_k &= \frac{ \sum_{n=1}^N s_{n,k} \mathbf{x}_n + \beta \mathbf{m} }{ \hat{\beta}_k } \tag{4.99} \end{align} $$

 計算結果をそれぞれbeta_hat_km_hat_kdとします。

 ちなみに$\sum_{n=1}^N s_{n,k} \mathbf{x}_n$の計算を確認してみると

round(x_nd[1:5, ], 2)
##       [,1]  [,2]
## [1,] -0.97  1.55
## [2,] -6.40 -2.22
## [3,]  3.69  2.00
## [4,]  1.93 -5.12
## [5,] -2.14  2.25
s_nk[1:5, k]
## [1] 0 0 0 1 0
round(s_nk[, k] * x_nd, 2)[1:5, ]
##      [,1]  [,2]
## [1,] 0.00  0.00
## [2,] 0.00  0.00
## [3,] 0.00  0.00
## [4,] 1.93 -5.12
## [5,] 0.00  0.00

対象としているクラスタのデータには1を掛けてそのままの値で、それ以外のクラスタのデータには0を掛けて値が0になっています。つまり対象となるクラスタのデータのみの和をとっていることが分かります。

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

# 初期化
nu_hat_k <- rep(0, K)
w_hat_ddk <- array(0, dim = c(D, D, K))

# lambdaの事後分布のパラメータを計算:式(4.103)
nu_hat_k[k] <- sum(s_nk[, k]) + nu
tmp_w1_dd <- t(s_nk[, k] * x_nd) %*% x_nd
tmp_w2_dd <- beta * matrix(m_d) %*% t(m_d)
tmp_w3_dd <- beta_hat_k[k] * matrix(m_hat_kd[k, ]) %*% t(m_hat_kd[k, ])
w_hat_ddk[, , k] <- solve(
  tmp_w1_dd + tmp_w2_dd - tmp_w3_dd + solve(w_dd)
)

# 確認
nu_hat_k
## [1] 70  0  0
round(w_hat_ddk, 3)
## , , 1
## 
##        [,1]   [,2]
## [1,]  0.004 -0.004
## [2,] -0.004  0.006
## 
## , , 2
## 
##      [,1] [,2]
## [1,]    0    0
## [2,]    0    0
## 
## , , 3
## 
##      [,1] [,2]
## [1,]    0    0
## [2,]    0    0

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

$$ \begin{align} \hat{\mathbf{W}}_k^{-1} &= \sum_{n=1}^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} \\ \hat{\nu}_k &= \sum_{n=1}^N s_{n,k} + \nu \tag{4.103} \end{align} $$

 計算結果をそれぞれnu_hat_kw_hat_ddkとします。

 パラメータが得られたので、ウィシャート分布に従い$\boldsymbol{\lambda}_k$をサンプルします。

# lambdaをサンプル:式(4.102)
lambda_ddk[, , k] <- rWishart(n = 1, df = nu_hat_k[k], Sigma = w_hat_ddk[, , k])

# 確認
lambda_ddk[, , k]
##            [,1]       [,2]
## [1,]  0.2204358 -0.2406156
## [2,] -0.2406156  0.4220070

 rWishart()で$\boldsymbol{\lambda}_k$を生成(サンプル)します。クラスタごとに生成するので、n引数に1を指定します。また自由度引数dfnu_hat_k[k]、スケール行列引数Sigmaw_hat_ddk[, , k]を指定します。

 多次元ガウス分布に従い$\boldsymbol{\mu}_k$もサンプルします。

# muをサンプル:式(4.98)
mu_kd[k, ] <- mvnfast::rmvn(
  n = 1, mu = m_hat_kd[k, ], sigma = solve(beta_hat_k[k] * lambda_ddk[, , k])
) %>% 
  as.vector()

# 確認
mu_kd[k, ]
## [1]  3.769087 -3.291953

 mvnfast::rmvn()で$\boldsymbol{\mu}_k$を生成(サンプル)します。こちらもクラスタごとに生成するので、n引数に1を指定します。また平均引数mum_hat_kd[k]、分散共分散行列引数sigmaに先ほどサンプルしたlambda_ddk[, , k]と係数beta_hat_k[k]との積をsolve()で逆行列に変換して指定します。
 サイズが1でもマトリクスで出力されるため、as.vector()でベクトルに変換してから代入しています。

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

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

# 初期化
alpha_hat_k <- rep(0, K)

# 混合比率のパラメータを計算:式(4.45)
alpha_hat_k <- colSums(s_nk) + alpha_k

# 確認
alpha_hat_k
## [1]  69 138  46

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

$$ \hat{\alpha}_k = \sum_{n=1}^N s_{n,k} + \alpha_k \tag{4.45} $$

 計算結果をalpha_hat_kとします。

 パラメータが得られたので、ディリクレ分布に従い$\boldsymbol{\pi}$をサンプルします。

# piをサンプル:式(4.44)
pi_k <- MCMCpack::rdirichlet(n = 1, alpha = alpha_hat_k) %>% 
  as.vector()

# 確認
pi_k
## [1] 0.2299443 0.5912498 0.1788058

 MCMCpack::rdirichlet()で$\boldsymbol{\pi}$を生成(サンプル)します。これは全データに対して共通するパラメータなので、size引数に1を指定します。またパラメータ引数alphaalpha_hat_kを指定します。
 こちらもマトリクスで出力されるため、as.vector()でベクトルに変換してから代入しています。

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

・結果の確認

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

# 作図用のデータフレームを作成
model_df <- tibble()
sample_df <- tibble()
for(k in 1:K) {
  # 近似事後分布を計算
  tmp_model_df <- cbind(
    point_df, 
    density = mvnfast::dmvn(
      as.matrix(point_df), mu = mu_kd[k, ], sigma = solve(lambda_ddk[, , k])
    ), 
    cluster = as.factor(k)
  )
  model_df <- rbind(model_df, tmp_model_df)
  
  # 観測データのクラスタを抽出
  k_idx <- which(s_nk[, k] == 1)
  tmp_sample_df <- tibble(
    x1 = x_nd[k_idx, 1], 
    x2 = x_nd[k_idx, 2], 
    cluster = as.factor(k)
  )
  sample_df <- rbind(sample_df, tmp_sample_df)
}

# 確認
head(model_df)
##      x1  x2      density cluster
## 1 -10.0 -10 3.381271e-07       1
## 2  -9.9 -10 4.044903e-07       1
## 3  -9.8 -10 4.825300e-07       1
## 4  -9.7 -10 5.740221e-07       1
## 5  -9.6 -10 6.809590e-07       1
## 6  -9.5 -10 8.055665e-07       1


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

# 近似事後分布を作図
ggplot() + 
  geom_contour(data = model_df, aes(x1, x2, z = density, color = cluster)) + # 近似事後分布
  geom_contour(data = model_true_df, aes(x1, x2, z = density, color = cluster), 
               linetype = "dotted", alpha = 0.6) + # 真の観測モデル
  geom_point(data = sample_df, aes(x1, x2, color = cluster)) + # 観測データ
  labs(title = "Gibbs Sampling", 
       subtitle = paste0('K=', K, ', N=', N, ', iter:', MaxIter), 
       x = expression(x[1]), y = expression(x[2]))

f:id:anemptyarchive:20201130190437p:plain
ギブスサンプリングによるガウス混合モデルの近似事後分布

 うまく推定できていることが確認できます。ただしクラスタの順番についてはランダムに決まります。

 実際には、始めの頃のサンプルは初期値に依存しているため破棄する期間(burn-in period)を設ける必要がありますが、この資料ではアルゴリズムの理解を焦点にしているためこれ以上は扱いません。

 サンプルされた$\boldsymbol{\mu},\ \boldsymbol{\Lambda}$の値の推移を確認しましょう。

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

 各クラスタと次元の値をそれぞれ取り出してデータフレームにまとめます。

# 作図用のデータフレームを作成
trace_mu_df <- tibble()
trace_lambda_df <- tibble()
for(k in 1:K) {
  for(d1 in 1:D) {
    # muの値を取得
    tmp_mu_df <- tibble(
      iteration = seq(0, MaxIter), 
      value = trace_mu_ikd[, k, d1], 
      label = as.factor(
        paste0("k=", k, ", d=", d1)
      )
    )
    trace_mu_df <- rbind(trace_mu_df, tmp_mu_df)
    
    for(d2 in 1:D) {
      # lambdaの値を取得
      tmp_lambda_df <- tibble(
        iteration = seq(0, MaxIter), 
        value = trace_lambda_iddk[, d1, d2, k], 
        label = as.factor(
          paste0("k=", k, ", d=", d1, ", d'=", d2)
        )
      )
      trace_lambda_df <- rbind(trace_lambda_df, tmp_lambda_df)
    }
  }
}

# 確認
head(trace_mu_df)
## # A tibble: 6 x 3
##   iteration value label   
##       <int> <dbl> <fct>   
## 1         0  7.6  k=1, d=1
## 2         1  1.44 k=1, d=1
## 3         2  1.66 k=1, d=1
## 4         3  2.60 k=1, d=1
## 5         4  2.52 k=1, d=1
## 6         5  2.20 k=1, d=1
head(trace_lambda_df)
## # A tibble: 6 x 3
##   iteration  value label         
##       <int>  <dbl> <fct>         
## 1         0 0.01   k=1, d=1, d'=1
## 2         1 0.0521 k=1, d=1, d'=1
## 3         2 0.0546 k=1, d=1, d'=1
## 4         3 0.0684 k=1, d=1, d'=1
## 5         4 0.0669 k=1, d=1, d'=1
## 6         5 0.0802 k=1, d=1, d'=1


 それぞれ折れ線グラフで可視化します。

# muの推移を確認
ggplot(trace_mu_df, aes(x = iteration, y = value, color = label)) + 
  geom_line() + 
  labs(title = expression(bold(mu)))
# lambdaの推移を確認
ggplot(trace_lambda_df, aes(x = iteration, y = value, color = label)) + 
  geom_line() + 
  labs(title = expression(bolditalic(Lambda)))

f:id:anemptyarchive:20201130190509p:plain
$\boldsymbol{\mu}$の推移

f:id:anemptyarchive:20201130190540p:plain
$\boldsymbol{\Lambda}$の推移


 以上がギブスサンプリングによる近似推論の処理です。

・おまけ

 gganimateパッケージを利用して、gif画像を作成するコードです。

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

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

# 追加パッケージ
library(gganimate)
MaxIter <- 150
# 作図用のデータフレームを作成
model_df <- tibble()
sample_df <- tibble()
for(i in 1:(MaxIter + 1)) {
  for(k in 1:K) {
    # 観測モデルを計算
    tmp_model_df <- cbind(
      point_df, 
      density = mvnfast::dmvn(
        as.matrix(point_df), mu = trace_mu_ikd[i, k, ], sigma = solve(trace_lambda_iddk[i, , , k])
      ), 
      cluster = as.factor(k), 
      iteration = as.factor(i-1)
    )
    model_df <- rbind(model_df, tmp_model_df)
    
    # 観測データのデータフレーム
    if(i > 1) { # 初期値以外のとき
      k_idx <- which(trace_s_in[i - 1, ] == k)
      tmp_sample_df <- tibble(
        x1 = x_nd[k_idx, 1], 
        x2 = x_nd[k_idx, 2], 
        cluster = as.factor(k), 
        iteration = as.factor(i - 1)
      )
      sample_df <- rbind(sample_df, tmp_sample_df)
    }
  }
  
  if(i == 1) { # 初期値のとき
    tmp_sample_df <- tibble(
      x1 = x_nd[, 1], 
      x2 = x_nd[, 2], 
      cluster = NA, 
      iteration = as.factor(i - 1)
    )
    sample_df <- rbind(sample_df, tmp_sample_df)
  }
  
  # 動作確認
  #print(paste0(i - 1, ' (', round((i - 1) / MaxIter * 100, 1), '%)'))
}

# 近似事後分布を作図
trace_graphe <- ggplot() + 
  geom_contour(data = model_df, aes(x1, x2, z = density, color = cluster)) + # 近似事後分布
  geom_contour(data = model_true_df, aes(x1, x2, z = density, color = cluster), 
               linetype = "dotted", alpha = 0.6) + # 真の観測モデル
  geom_point(data = sample_df, aes(x1, x2, color = cluster)) + # 観測データ
  transition_manual(iteration) + # フレーム
  labs(title = "Gibbs Sampling", 
       subtitle = paste0('K=', K, ', N=', N, ', iter:{current_frame}'), 
       x = expression(x[1]), y = expression(x[2]))

# gif画像を作成
animate(trace_graphe, nframes = MaxIter + 1, fps = 10)

 先ほど同様に、イタレーションごとに近似事後分布の確率密度を計算して、データフレームに結合していきます。gganimate::transition_manual()に、フレームの切り替えの基準となる列を指定します。全てのフレーム(イタレーション)のグラフがtrace_grapheに格納されるので、gganimate::animate()でgif画像として出力します。
 処理に時間がかかる場合は、描画範囲や点の間隔(x_line)を調整してください。

f:id:anemptyarchive:20201130190634g:plain
ギブスサンプリングによる近似事後分布の推移

 (作図処理に時間がかかるので150回までで切っています。)

参考文献

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

おわりに

 ブログ開設2年目の最終日です。170ちょい目の記事です。これからもよろしくお願いします。

 2020年11月30日はモーニング娘。'20の加賀楓さん21歳のお誕生日です!おめでとうございます!

 曲後半でセンターしてるショートボブの方です!
 加賀温泉郷の観光大使もされているので、そちらも是非行きたいっ!!!

kaineya.jp

 来年こそは、加賀!温泉♨紅葉狩り🍁

【次節の内容】

www.anarchive-beta.com