からっぽのしょこ

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

【R】10.1.3:一変数ガウス分布の変分推論【PRMLのノート】

はじめに

 『パターン認識と機械学習』の独学時のまとめです。一連の記事は「数式の行間埋め」または「R・Pythonでの実装」からアルゴリズムの理解を補助することを目的としています。本とあわせて読んでください。

 この記事は、10.1.3項の内容です。平均と精度が未知の1次元ガウス分布(正規分布)に対する変分推論(変分ベイズ)をR言語で実装します。

【数式読解編】

www.anarchive-beta.com

【他の節一覧】

www.anarchive-beta.com

【この節の内容】

・Rで実装

 人工データを用いて、簡単な変分推論を行ってみましょう。

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

# 10.1.3で利用するパッケージ
library(tidyverse)


・真の分布の設定

 まずは、真の分布を設定します。この項では、真の分布を平均$\mu$、精度$\tau$の1次元ガウス分布$\mathcal{N}(x | \mu, \tau^{-1})$とします。

 真の分布のパラメータを設定します。

# 真の平均パラメータを指定
mu_truth <- 5

# 真の精度パラメータを指定
tau_truth <- 0.5

 真の平均パラメータ$\mu$をmu_truthとして値を指定します。
 真の精度パラメータ$\tau$をtau_truthとして0より大きい値$\tau > 0$を指定します。精度は分散の逆数$\sigma^2 = \frac{1}{\tau}$です。値が大きいほど、散らばり具合が小さくなります。
 この例では、この2つのパラメータの値を分布推定するのが目的です。

 真の分布をグラフで確認しましょう。作図用の点を作成します。

# 作図用のxの値を作成
x_vec <- seq(
  mu_truth - 4 * sqrt(1 / tau_truth), 
  mu_truth + 4 * sqrt(1 / tau_truth), 
  length.out = 1000
)

 作図用に、ガウス分布に従う変数$x$がとり得る値をseq()で作成してx_vecとします。この例では、平均値を中心に標準偏差の4倍を範囲とします。length.out引数を使うと指定した要素数で等間隔に切り分けます。by引数を使うと切り分ける間隔を指定できます。処理が重い場合は、この値を調整してください。

 真の分布を計算します。

# 真の分布を計算
model_df <- tibble(
  x = x_vec, 
  density = dnorm(x = x, mean = mu_truth, sd = sqrt(1 / tau_truth))
)

 x_vecの値ごとに確率密度を計算します。1次元ガウス分布の確率密度は、dnorm()で計算できます。値の引数xx_vec、平均の引数meanmu_truth、標準偏差の引数sdtauの逆数の平方根sqrt(1 / tau_truth)を指定します。精度$\tau$は分散の逆数なので、その平方根が真の標準偏差$\sigma = \sqrt{\frac{1}{\tau}}$になります。

 作成したデータフレームは次のようになります。

# 確認
head(model_df)
## # A tibble: 6 x 2
##        x   density
##    <dbl>     <dbl>
## 1 -0.657 0.0000946
## 2 -0.646 0.0000977
## 3 -0.634 0.000101 
## 4 -0.623 0.000104 
## 5 -0.612 0.000108 
## 6 -0.600 0.000111


 真の分布を作図します。

# 真のを作図
ggplot(model_df, aes(x = x, y = density)) + 
  geom_line(color = "blue") + # 真の分布
  labs(title = "Gaussian Distribution", 
       subtitle = paste0("mu=", mu_truth, ", tau=", tau_truth))

f:id:anemptyarchive:20210602013125p:plain
真の分布:1次元ガウス分布


・観測データの生成

 続いて、設定した真の分布に従って観測データ$\mathcal{D} = \{x_1, x_2, \cdots, x_N\}$を生成します。

 ガウス分布に従う$N$個のデータをランダムに生成します。

# (観測)データ数を指定
N <- 50

# ガウス分布に従うデータを生成
x_n <- rnorm(n = N, mean = mu_truth, sd = sqrt(1 / tau_truth))

 生成するデータ数$N$をNとして値を指定します。

 ガウス分布に従う乱数は、rnorm()で生成できます。データ数の引数nNを指定します。他の引数についてはdnorm()と同じです。生成したN個のデータをx_nとします。

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

# 観測データをデータフレームに格納
data_df <- tibble(x_n = x_n)

# 確認
head(data_df)
## # A tibble: 6 x 1
##     x_n
##   <dbl>
## 1  4.46
## 2  5.56
## 3  3.76
## 4  3.48
## 5  4.76
## 6  6.16


 観測データのヒストグラムを作成します。

# 観測データのヒストグラムを作成
ggplot() + 
  geom_histogram(data = data_df, aes(x = x_n), binwidth = 0.5) + # 観測データ:(度数)
  #geom_histogram(data = data_df, aes(x = x_n, y = ..density..), binwidth = 1) + # 観測データ:(相対度数)
  #geom_line(data = model_df, aes(x = x, y = density), 
  #          color = "red", linetype = "dashed") + # 真の分布:(相対度数用)
  labs(title = "Gaussian Distribution", 
       subtitle = paste0("N=", N, ", mu=", mu_truth, ", tau=", tau_truth), 
       x = "x")

f:id:anemptyarchive:20210602013211p:plainf:id:anemptyarchive:20210602013214p:plain
観測データ:1次元ガウス分布

 左の図のy軸はデータ数です。カウントする間隔はbinwidthで調整できます。右の図のy軸は相対度数です。y = ..density..を指定すると相対度数のヒストグラムになります。

・事前分布の設定

 次に、尤度に対する共役事前分布を設定します。1次元ガウス分布の平均・精度パラメータ$\mu,\ \tau$の同時事前分布$p(\mu, \tau)$として、ガウス-ガンマ分布$\mathrm{NG}(\mu, \tau | \mu_0, \lambda_0, a_0, b_0)$を設定します。

 $\mu$と$\tau$の事前分布のパラメータ(超パラメータ)を指定します。

# muの事前分布のパラメータを指定
mu_0 <- 0
lambda_0 <- 1

# tauの事前分布のパラメータを指定
a_0 <- 1
b_0 <- 1

 $\mu$の事前分布$\mathcal{N}(\mu | \mu_0, (\lambda_0 \tau)^{-1})$の平均パラメータ$\mu_0$と精度パラメータの係数$\lambda_0$をそれぞれmu_0, lambda_0として$\lambda_0 > 0$の値を指定します。
 $\tau$の事前分布$\mathrm{Gam}(\tau | a_0, b_0)$のパラメータ$a_0,\ b_0$をa_0, b_0として$a_0 > 0,\ b_0 > 0$の値を指定します。

 $\mu,\ \tau$の事前分布をグラフで確認しましょう。作図用の点を作成します。

# 作図用のmuの値を作成
mu_vec <- seq(
  mu_truth - 4 * sqrt(1 / tau_truth), 
  mu_truth + 4 * sqrt(1 / tau_truth), 
  length.out = 500
)

# 作図用のtauの値を作成
tau_vec <- seq(0, tau_truth + 2, length.out = 500)

 1次元ガウス分布に従う変数$\mu$がとり得る値を作成してmu_vecとします。この例では、真の平均値を中心に真の標準偏差の4倍を範囲とします。
 ガンマ分布に従う変数$\tau$がとり得る値を作成してtau_vecとします。この例では、0からtau + 2を範囲とします。

 格子状の点を作成します。

# 作図用の点を作成
mu_tau_point_mat <- cbind(
  mu = rep(mu_vec, times = length(tau_vec)), 
  tau = rep(tau_vec, each = length(mu_vec))
)

 mu_vectau_vecの要素の全ての組み合わせを持つようにマトリクスを作成してmu_tau_point_matとします。これは、ガウス-ガンマ分布の確率密度を等高線図にする際に格子状の点(2軸の全ての値が直交する点)を渡す必要があるためです。

 作成した$\mu$と$\tau$の点は次のようになります。

# 確認
mu_tau_point_mat[1:5, ]
##              mu tau
## [1,] -0.6568542   0
## [2,] -0.6341815   0
## [3,] -0.6115087   0
## [4,] -0.5888360   0
## [5,] -0.5661632   0

 1列目が$\mu$、2列目が$\tau$の値です。

 $\mu,\ \tau$の事前分布を計算します。

# 同時事前分布を計算
prior_df <- tibble(
  mu = mu_tau_point_mat[, 1], 
  tau = mu_tau_point_mat[, 2], 
  N_dens = dnorm(x = mu, mean = mu_0, sd = sqrt(1 / lambda_0 / tau)), # muの事前分布
  Gam_dens = dgamma(x = tau, shape = a_0, rate = b_0), # tauの事前分布
  density = N_dens * Gam_dens
)

 ガウス-ガンマ分布の確率密度は、次のガウス分布とガンマ分布の確率密度の積です。

$$ \mathrm{NG}(\mu, \tau | \mu_0, \lambda_0, a_0, b_0) = \mathcal{N}(\mu | \mu_0, (\lambda_0 \tau)^{-1}) \mathrm{Gam}(\tau | a_0, b_0) $$

 真の分布のときと同様に、ガウス分布の確率密度をdnorm()で計算します。引数xには格子状の点に変換した$\mu$の値mu_tau_point_mat[, 1]を指定します。また、標準偏差は$\sqrt{\frac{1}{\lambda_0 \tau}}$です。$\tau$として格子状の点に変換した$\tau$の値mu_tau_point_mat[, 2]を使います。
 ガンマ分布の確率密度は、dgamma()で計算できます。値の引数xmu_tau_point_mat[, 2]、パラメータの引数shapea_0rateb_0を指定します。

 作成したデータフレームは次のようになります。

# 確認
head(prior_df)
## # A tibble: 6 x 5
##       mu   tau N_dens Gam_dens density
##    <dbl> <dbl>  <dbl>    <dbl>   <dbl>
## 1 -0.657     0      0        1       0
## 2 -0.634     0      0        1       0
## 3 -0.612     0      0        1       0
## 4 -0.589     0      0        1       0
## 5 -0.566     0      0        1       0
## 6 -0.543     0      0        1       0


 $\mu,\ \tau$の事前分布を作図します。

# 真のパラメータをデータフレームに格納
true_params_df <- tibble(mu = mu_truth, tau = tau_truth)

# 事前分布を作図
ggplot() + 
  geom_contour(data = prior_df, aes(x = mu, y = tau, z = density, color = ..level..)) + # 事前分布
  geom_point(data = true_params_df, aes(x = mu, y = tau), color = "red", shape = 4, size = 5) + # 真の値
  labs(title = "Gaussian-Gamma Distribution", 
       subtitle = paste0("mu_0=", mu_0, ", lambda_0=", lambda_0, 
                         ", a_0=", a_0, ", b_0=", b_0), 
       x = expression(mu), y = expression(tau))

f:id:anemptyarchive:20210602013332p:plain
事前分布:ガウス・ガンマ分布

 等高線図はgeom_contour()で描画します。また、真のパラメータをデータフレームに格納して、事前分布と重ねてプロットしています。

・真の事後分布の計算

 この項では、変分推論により事後分布の近似分布$q(\mu, \tau)$を求めます。ただし、この問題は解析的に求められます。解析的に求めた真の事後分布$p(\mu, \tau | \mathcal{D})$を確認します。事後分布もガウス-ガンマ分布$\mathrm{NG}(\mu, \tau | \hat{\mu}, \hat{\lambda}, \hat{a}, \hat{b})$になります。

 $\mu$の真の事後分布のパラメータを計算します。

# muの真の事後分布のパラメータを計算
lambda_hat <- lambda_0 + N
mu_hat <- (lambda_0 * mu_0 + sum(x_n)) / lambda_hat

 解析的に求めた$\mu$の事後分布の平均パラメータ$\hat{\mu}$、精度パラメータの係数$\hat{\lambda}$を次の式で計算します。

$$ \begin{aligned} \hat{\lambda} &= \lambda_0 + N \\ \hat{\mu} &= \frac{1}{\hat{\lambda}} \left( \lambda_0 \mu_0 + \sum_{n=1}^N x_n \right) \end{aligned} $$

 計算結果は次のようになります。

# 確認
mu_hat
lambda_hat
## [1] 4.976428
## [1] 51


 $\tau$の真の事後分布のパラメータを計算します。

# lambdaの真の事後分布のパラメータを計算
a_hat <- a_0 + 0.5 * N
b_hat <- b_0 + 0.5 * (sum(x_n^2) + lambda_0 * mu_0^2 - lambda_hat * mu_hat^2)

 解析的に求めた$\tau$の事後分布のパラメータ$\hat{a},\ \hat{b}$を次の式で計算します。

$$ \begin{aligned} \hat{a} &= a_0 + \frac{N}{2} \\ \hat{b} &= \frac{1}{2} \left( \sum_{n=1}^N x_n + \lambda_0 \mu_0^2 + \hat{\lambda} \hat{\mu}^2 \right) \end{aligned} $$

 計算結果は次のようになります。

# 確認
a_hat
b_hat
## [1] 26
## [1] 47.89597


 $\mu,\ \tau$の真の事後分布を作図します。

# 真の同時事後分布を計算
true_posterior_df <- tibble(
  mu = mu_tau_point_mat[, 1], 
  tau = mu_tau_point_mat[, 2], 
  N_dens = dnorm(x = mu, mean = mu_hat, sd = sqrt(1 / lambda_hat / tau)), # muの真の事後分布
  Gam_dens = dgamma(x = tau, shape = a_hat, rate = b_hat), # tauの真の事後分布
  density = N_dens * Gam_dens
)

# 真の事後分布を作図
ggplot() + 
  geom_contour(data = true_posterior_df, aes(x = mu, y = tau, z = density, color = ..level..)) + # 真の事後分布
  geom_point(data = true_params_df, aes(x = mu, y = tau), color = "red", shape = 4, size = 5) + # 真の値
  labs(title = "Gaussian-Gamma Distribution", 
       subtitle = paste0("N=", N, 
                         ", mu_hat=", round(mu_hat, 1), ", lambda_hat=", lambda_hat, 
                         ", a_hat=", a_hat, ", b_hat=", round(b_hat, 1)), 
       x = expression(mu), y = expression(tau))

f:id:anemptyarchive:20210602013437p:plain
真の事後分布:ガウス・ガンマ分布

 事前分布のときと同様に計算して作図します。

 各パラメータの真の事後分布も確認しておきましょう。

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

 $\mu$の真の事後分布を作図します。

# muの真の事後分布を作図
tibble(
  mu = mu_vec, 
  N_dens = dnorm(x = mu, mean = mu_hat, sd = sqrt(1 / (lambda_hat * a_hat / b_hat)))
) %>% 
  ggplot(aes(x = mu, y = N_dens)) + 
    geom_line() + # muの真の事後分布
    geom_vline(xintercept = mu_truth, color = "red", linetype = "dashed") + # muの真の値
    labs(title = "Gaussian Distribution", 
         subtitle = paste0("N=", N, 
                           ", mu_hat=", round(mu_hat, 1), 
                           ", lambda_hat=", lambda_hat, 
                           ", E[tau]=", round(a_hat / b_hat, 5)), 
         x = expression(mu), y = "density")

 重複していない$\mu$の値mu_vecの要素ごとに確率密度を計算します。ただし、$\tau$の代わりに、ガンマ分布の期待値(B.27)より、$\tau$の期待値$\mathbb{E}[\tau] = \frac{\hat{a}}{\hat{b}}$を用いて確率密度を計算します。
 真の値を垂直線geom_vline()で示しています。

 同様に、$\tau$の真の事後分布を作図します。

# tauの真の事後分布を作図
tibble(
  tau = tau_vec, 
  Gam_dens = dgamma(x = tau, shape = a_hat, rate = b_hat)
) %>% 
  ggplot(aes(x = tau, y = Gam_dens)) + 
    geom_line() + # tauの真の事後分布
    geom_vline(xintercept = tau_truth, color = "red", linetype = "dashed") + # tauの真の値
    labs(title = "Gamma Distribution", 
         subtitle = paste0("N=", N, ", a_hat=", a_hat, ", b_hat=", round(b_hat, 1)), 
         x = expression(tau), y = "density")

 こちらも重複していない$\tau$の値tau_vecの要素ごとに確率密度を計算します。


f:id:anemptyarchive:20210602013507p:plainf:id:anemptyarchive:20210602013511p:plain
各パラメータの真の事後分布:ガウス分布、ガンマ分布


 以上でモデルの構築ができたので、次は推論処理を行います。

・変分推論

 変分推論により$\mu,\ \tau$の近似事後分布$q(\mu, \tau)$を求めます。近似事後分布もガウス-ガンマ分布$\mathrm{NG}(\mu, \tau | \mu_N, \lambda_N, a_N, b_N)$になります。
 変分推論では、$\mu$の近似事後分布$q_{\mu}(\mu)$(のパラメータ$\mu_N,\ \lambda_N$)と$\tau$の近似事後分布$q_{\tau}(\tau)$(のパラメータ$a_N,\ b_N$)の更新を交互に行います。

・コード全体

 事前分布のパラメータの値を初期値として持つように近似事後分布のパラメータを作成しておきます。
 また、後で各パラメータの更新値や分布の推移を確認するために、それぞれ試行ごとにtrace_***に値を保存していきます。ただし、この例はすぐに収束してしまうので、$q_{\mu}(\mu)$と$q_{\tau}(\tau)$をそれぞれ更新する度に値を記録することにします。

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

# 初期値を代入
mu_N <- mu_0
lambda_N <- lambda_0 * a_0 / b_0
a_N <- a_0
b_N <- b_0

# 推移の確認用の受け皿を作成
trace_mu_i     <- rep(0, MaxIter * 2 + 1)
trace_lambda_i <- rep(0, MaxIter * 2 + 1)
trace_a_i      <- rep(0, MaxIter * 2 + 1)
trace_b_i      <- rep(0, MaxIter * 2 + 1)

# 初期値を記録
trace_mu_i[1]     <- mu_0
trace_lambda_i[1] <- lambda_0
trace_a_i[1]      <- a_0
trace_b_i[1]      <- b_0

# 変分推論
for(i in 1:MaxIter) {
  
  # muの近似事後分布のパラメータを計算:式(10.26)(10.27)
  mu_N <- (lambda_0 * mu_0 + sum(x_n)) / (lambda_0 + N)
  lambda_N <- (lambda_0 + N) * a_N / b_N
  
  # i回目のmuの近似事後分布の更新後の結果を記録
  trace_mu_i[i * 2]     <- mu_N
  trace_lambda_i[i * 2] <- lambda_N
  trace_a_i[i * 2]      <- a_N
  trace_b_i[i * 2]      <- b_N
  
  # lambdaの近似事後分布のパラメータを計算:式(10.29)(10.30)
  a_N <- a_0 + 0.5 * (N + 1)
  term_x <- 0.5 * (lambda_0 * mu_0^2 + sum(x_n^2))
  term_mu2 <- 0.5 * (lambda_0 + N) * (mu_N^2 + 1 / lambda_N)
  term_mu <- (lambda_0 * mu_0 + sum(x_n)) * mu_N
  b_N <- b_0 + term_x + term_mu2 - term_mu
  
  # i回目のlambdaの近似事後分布の更新後の結果を記録
  trace_mu_i[i * 2 + 1]     <- mu_N
  trace_lambda_i[i * 2 + 1] <- lambda_N
  trace_a_i[i * 2 + 1]      <- a_N
  trace_b_i[i * 2 + 1]      <- b_N
  
  # 動作確認
  #print(paste0(i, ' (', round(i / MaxIter * 100, 1), '%)'))
}


・処理の解説

 続いて、変分推論で行う処理を確認していきます。

・平均パラメータの近似事後分布の更新

 先に、$q_{\tau}(\tau)$を固定した下での$\mu$の近似事後分布の最適解$q_{\mu}^{*}(\mu)$を計算します。

 $\mu$の近似事後分布のパラメータを計算します。 

# muの近似事後分布のパラメータを計算:式(10.26)(10.27)
mu_N <- (lambda_0 * mu_0 + sum(x_n)) / (lambda_0 + N)
lambda_N <- (lambda_0 + N) * a_N / b_N

 変分推論により求めた$\mu$の事後分布の平均パラメータ$\mu_N$、精度パラメータ$\lambda_N$を次の式で計算します。

$$ \begin{align} \mu_N &= \frac{ \lambda_0 \mu_0 + \sum_{n=1}^N x_n }{ \lambda_0 + N } \tag{10.26}\\ \lambda_N &= (\lambda_0 + N) \mathbb{E}_{q_{\tau}(\tau)}[\tau] \tag{10.27}\\ \mathbb{E}_{q_{\tau}^{*}(\tau)}[\tau] &= \frac{a_N}{b_N} \end{align} $$

 $\lambda_N$は、これまでのように精度パラメータの係数ではなく、精度パラメータそのものであることに注意してください。

 計算結果は次のようになります。

# 確認
mu_N
lambda_N
## [1] 4.976428
## [1] 28.19653


・精度パラメータの近似事後分布の更新

 続いて、$q_{\mu}(\mu)$を固定した下での$\tau$の近似事後分布の最適解$q_{\tau}^{*}(\tau)$を計算します。

 $\tau$の近似事後分布のパラメータを計算します。 

# lambdaの近似事後分布のパラメータを計算:式(10.29)(10.30)
a_N <- a_0 + 0.5 * (N + 1)
term_x <- 0.5 * (lambda_0 * mu_0^2 + sum(x_n^2))
term_mu2 <- 0.5 * (lambda_0 + N) * (mu_N^2 + 1 / lambda_N)
term_mu <- (lambda_0 * mu_0 + sum(x_n)) * mu_N
b_N <- b_0 + term_x + term_mu2 - term_mu

 変分推論により求めた$\tau$の事後分布のパラメータを$a_N,\ b_N$を次の式で計算します。

$$ \begin{align} a_N &= a_0 + \frac{N + 1}{2} \tag{10.29}\\ b_N &= b_0 + \frac{1}{2} \left( \lambda_0 \mu_0^2 + \sum_{n=1}^N x_n^2 \right) + \frac{\lambda_0 + N}{2} \mathbb{E}_{g_{\mu}(\mu)} [\mu^2] - \left( \lambda_0 \mu_0 + \sum_{n=1}^N x_n \right) \mathbb{E}_{g_{\mu}(\mu)} [\mu] \tag{10.30'}\\ \mathbb{E}_{g_{\mu}^{*}(\mu)} [\mu] &= \mu_N \\ \mathbb{E}_{g_{\mu}^{*}(\mu)} [\mu^2] &= \mu_N^2 + \frac{1}{\lambda_N} \end{align} $$

 計算結果は次のようになります。

# 確認
a_N
b_N
## [1] 26.5
## [1] 47.93144


 これで超パラメータ$\mu_N,\ \lambda_N$と$a_N,\ b_N$を更新できました。
 以上が変分推論で行う個々の処理です。繰り返し行うことで、徐々に真の事後分布に近付けていきます。閾値に達するまでループすることもできますが、この例では簡単に、試行回数を設定してループすることにします。

・推論結果の確認

 推論結果を確認していきます。

・近似事後分布の確認

 mu_N, lambda_Na_N, b_Nの最後の更新値を用いて、$\mu,\ \tau$の近似事後分布を計算します。

# 同時近似事後分布を計算
posterior_df <- tibble(
  mu = mu_tau_point_mat[, 1], 
  tau = mu_tau_point_mat[, 2], 
  N_dens = dnorm(x = mu, mean = mu_N, sd = sqrt(1 / (lambda_N * b_N / a_N * tau))), # muの近似事後分布
  Gam_dens = dgamma(x = tau, shape = a_N, rate = b_N), # tauの近似事後分布
  density = N_dens * Gam_dens
)

 事前分布のときと同様に計算します。ただし、$\mu$の近似事後分布の精度パラメータ$\lambda_N$には、既に$\tau$の情報が期待値$\mathbb{E}[\tau]$として含まれています。そこで、$\mathbb{E}[\tau] = \frac{a_N}{b_N}$の逆数を掛けて打ち消してから、$\tau$の値mu_tau_point_mat[, 2]を掛けて確率密度を計算します。(本当にこんなことをしていいのかよく分かっていません。やるにしてももう少し納得できる理由が欲しい。)

 $\mu,\ \tau$の同時近似事後分布を作図します。

# 近似事後分布を作図
ggplot() + 
  geom_contour(data = posterior_df, aes(x = mu, y = tau, z = density, color = ..level..)) + # 近似事後分布
  geom_contour(data = true_posterior_df, aes(x = mu, y = tau, z = density, color = ..level..), 
               alpha = 0.5, linetype = "dashed") + # 真の事後分布
  geom_point(data = true_params_df, aes(x = mu, y = tau), color = "red", shape = 4, size = 5) + # 真の値
  labs(title = "Gaussian-Gamma Distribution", 
       subtitle = paste0("N=", N, 
                         ", mu_N=", round(mu_N, 1), ", lambda_N=", round(lambda_N, 5), 
                         ", a_N=", a_N, ", b_N=", round(b_N, 1)), 
       x = expression(mu), y = expression(tau))

f:id:anemptyarchive:20210602013713p:plain
近似事後分布:ガウス・ガンマ分布

 事前分布のときと同様に作図できます。真の事後分布と重ねて描画しています。
 真の事後分布を近似できています。

 各パラメータの近似事後分布も確認しましょう。

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

 $\mu$の近似事後分布を作図します。

# muの近似事後分布を作図
tibble(
  mu = mu_vec, 
  N_dens = dnorm(x = mu, mean = mu_N, sd = sqrt(1 / lambda_N))
) %>% 
  ggplot(aes(x = mu, y = N_dens)) + 
    geom_line() + # muの近似事後分布
    geom_vline(xintercept = mu_truth, color = "red", linetype = "dashed") + # muの真の値
    labs(title = "Gaussian Distribution", 
         subtitle = paste0("N=", N, 
                           ", mu_N=", round(mu_N, 1), 
                           ", E[tau]=", round(a_N / b_N, 5)), 
        x = expression(mu), y = "density")

 こちらは、lambda_Nをそのまま精度パラメータとして用います。

 $\tau$の近似事後分布を作図します。

# tauの近似事後分布を作図
tibble(
  tau = tau_vec, 
  Gam_dens = dgamma(x = tau, shape = a_N, rate = b_N)
) %>% 
  ggplot(aes(x = tau, y = Gam_dens)) + 
    geom_line() + # tauの近似事後分布
    geom_vline(xintercept = tau_truth, color = "red", linetype = "dashed") + # tauの真の値
    labs(title = "Gamma Distribution", 
         subtitle = paste0("N=", N, ", a_N=", a_N, ", b_N=", round(b_N, 1)), 
         x = expression(tau), y = "density")

 真の事後分布のときと同様に作図できます。


f:id:anemptyarchive:20210602013733p:plainf:id:anemptyarchive:20210602013736p:plain
各パラメータの近似事後分布:ガウス分布、ガンマ分布


・超パラメータの推移の確認

 超パラメータの更新値の推移を確認します。

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

 学習ごとの値を格納した各パラメータの配列trace_***tidyr::pivot_longer()で縦長のデータフレームに変換して作図します。

# mu_Nの推移を作図
tibble(
  iteration = 0:(MaxIter * 2) * 0.5, 
  value = trace_mu_i
) %>% 
  ggplot(aes(x = iteration, y = value)) + 
    geom_line() + 
    labs(title = "Variational Inference", 
         subtitle = expression(mu[N]))
# lambda_Nの推移を作図
tibble(
  iteration = 0:(MaxIter * 2) * 0.5, 
  value = trace_lambda_i
) %>% 
  ggplot(aes(x = iteration, y = value)) + 
    geom_line() + 
    labs(title = "Variational Inference", 
         subtitle = expression(lambda[N]))
# a_Nの推移を作図
tibble(
  iteration = 0:(MaxIter * 2) * 0.5, 
  value = trace_a_i
) %>% 
  ggplot(aes(x = iteration, y = value)) + 
    geom_line() + 
    labs(title = "Variational Inference", 
         subtitle = expression(a[N]))
# b_Nの推移を作図
tibble(
  iteration = 0:(MaxIter * 2) * 0.5, 
  value = trace_b_i
) %>% 
  ggplot(aes(x = iteration, y = value)) + 
    geom_line() + 
    labs(title = "Variational Inference", 
         subtitle = expression(b[N]))


f:id:anemptyarchive:20210602013814p:plainf:id:anemptyarchive:20210602013818p:plain
平均パラメータンの近似事後分布のパラメータの推移

f:id:anemptyarchive:20210602013840p:plainf:id:anemptyarchive:20210602013844p:plain
精度パラメータの近似事後分布のパラメータの推移

 平均パラメータに関しては1回目の試行で、精度パラメータに関しては2回目の試行で収束しているのが確認できます。

・おまけ:アニメーションによる推移の確認

 gganimateパッケージを利用して、分布の推移のアニメーション(gif画像)を作成するためのコードです。

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

# 追加パッケージ
library(gganimate)


 事後分布の推移を確認します。作図用のデータフレームを作成します。

# 作図用のデータフレームを作成
trace_posterior_df <- tibble()
for(i in 1:(MaxIter * 2 + 1)) {
  
  # i回目の近似事後分布を計算
  tmp_posterior_df <- tibble(
    mu = mu_tau_point_mat[, 1], 
    tau = mu_tau_point_mat[, 2], 
    N_dens = dnorm(
      x = mu, 
      mean = trace_mu_i[i], 
      sd = sqrt(1 / (trace_lambda_i[i] * trace_b_i[i] / trace_a_i[i] * tau))
    ), # muの近似事後分布
    Gam_dens = dgamma(x = tau, shape = trace_a_i[i], rate = trace_b_i[i]), # tauの近似事後分布
    density = N_dens * Gam_dens, 
    label = as.factor(
      paste0("iter:", (i - 1) * 0.5, ", N=", N, 
             ", mu_N=", round(trace_mu_i[i], 1), ", lambda_N=", round(trace_lambda_i[i], 5), 
             ", a_N=", trace_a_i[i], ", b_N=", round(trace_b_i[i], 1))
    ) # フレーム切替用のラベル
  )
  
  # 結合
  trace_posterior_df <- rbind(trace_posterior_df, tmp_posterior_df)
  
  # 動作確認
  #print(paste0((i - 1) * 0.5, ' (', round((i - 1) * 0.5 / MaxIter * 100, 1), '%)'))
}

 各試行における計算結果を同じデータフレームに格納していく必要があります。処理に時間がかかる場合は、mu_tau_point_matの範囲や間隔を調整してください。

 作成したデータフレームは次のようになります。

# 確認
head(trace_posterior_df)
## # A tibble: 6 x 6
##       mu   tau N_dens Gam_dens density label                                    
##    <dbl> <dbl>  <dbl>    <dbl>   <dbl> <fct>                                    
## 1 -0.657     0      0        1       0 iter:0, N=50, mu_N=0, lambda_N=1, a_N=1,~
## 2 -0.634     0      0        1       0 iter:0, N=50, mu_N=0, lambda_N=1, a_N=1,~
## 3 -0.612     0      0        1       0 iter:0, N=50, mu_N=0, lambda_N=1, a_N=1,~
## 4 -0.589     0      0        1       0 iter:0, N=50, mu_N=0, lambda_N=1, a_N=1,~
## 5 -0.566     0      0        1       0 iter:0, N=50, mu_N=0, lambda_N=1, a_N=1,~
## 6 -0.543     0      0        1       0 iter:0, N=50, mu_N=0, lambda_N=1, a_N=1,~

 フレーム切替用のラベルの列labelをパラメータの値を使って作成していますが、ややこしければas.factor(paste0("iter:", (i - 1) * 0.5))だけにして試行回数を表示するだけでもそれっぽくなります。

 分布の推移をgif画像で出力します。

# 近似事後分布を作図
trace_graph <- ggplot() + 
  geom_contour(data = trace_posterior_df, aes(x = mu, y = tau, z = density, color = ..level..)) + # 近似事後分布
  geom_contour(data = true_posterior_df, aes(x = mu, y = tau, z = density, color = ..level..), 
               alpha = 0.5, linetype = "dashed") + # 真の事後分布
  geom_point(data = true_params_df, aes(x = mu, y = tau), color = "red", shape = 4, size = 5) + # 真の値
  gganimate::transition_manual(label) + # フレーム
  labs(title = "Gaussian-Gamma Distribution", 
       subtitle = "{current_frame}", 
       x = expression(mu), y = expression(tau))

# gif画像を作成
gganimate::animate(trace_graph, nframes = MaxIter * 2 + 1, fps = 5)

 各フレームの順番を示す列をgganimate::transition_manual()に指定します。初期値を含み各試行で2回記録しているため、フレーム数はMaxIter * 2 + 1です。


f:id:anemptyarchive:20210602013914g:plain
近似事後分布の推移:ガウス・ガンマ分布

 初期値(事前分布)は分散が大きく確率密度が小さいため表示されていません。平均は初期値が0なので、0.5回目の試行で収束しています。精度は2回目の試行で収束しているのが分かります。

参考文献

  • C.M.ビショップ著,元田 浩・他訳『パターン認識と機械学習 上下』,丸善出版,2012年.

おわりに

 本では無情報事前分布として全ての超パラメータを0にすると書いてあるけど0はダメなのでは?それに初期値が0だと図10.4の(a)のグラフにはならないと思うし。それとは別に、(あんまり努力はしてないけど)図10.4を再現できなかった。今回は色々ちょっと自信ない。

 そんなことよりこのブログを投稿する数時間前に、元つばきファクトリーの小片リサさんの活動再開が発表されました!ソロデビューおめでとう!!

 同時に公開されたカバー動画です♪何はともあれ良かった良かった。

【関連する内容】

 解析的に解く(真の事後分布を求める)はこちらの記事でやりました。

www.anarchive-beta.com