からっぽのしょこ

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

【R】3.4.1:多次元ガウス分布の学習と予測:平均が未知の場合【緑ベイズ入門のノート】

はじめに

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

 この記事は、項の内容です。「尤度関数を平均が未知の多次元ガウス分布(多変量正規分布)」、「事前分布を多次元ガウス分布」とした場合の「パラメータの事後分布」と「未観測値の予測分布」の計算をR言語で実装します。

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

【数式読解編】

www.anarchive-beta.com

【他の節の内容】

www.anarchive-beta.com

【この節の内容】

・Rでやってみよう

 人工的に生成したデータを用いて、ベイズ推論を行ってみましょう。

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

# 3.4.1項で利用パッケージ
library(tidyverse)
library(mvnfast)

 mvnfastパッケージは、多次元ガウス分布に関するパッケージです。多次元ガウス分布に従う乱数生成関数rmvn()と確率密度関数dmvn()を使います。

・観測モデルの構築

 まずは、観測モデルを設定します。この例では、尤度$p(\mathbf{X} | \boldsymbol{\mu})$を多次元ガウス分布$\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}, \boldsymbol{\Lambda}^{-1})$とします。

 尤度のパラメータを設定します。この実装例では2次元のグラフで表現するため、$D = 2$のときのみ動作します。

# 真のパラメータを指定
mu_truth_d <- c(25, 50)
sigma_dd <- matrix(c(20, 15, 15, 30), nrow = 2, ncol = 2)
lambda_dd <- solve(sigma_dd^2)

 平均パラメータ$\boldsymbol{\mu} = (\mu_1, \cdots, \mu_D)$をmu_truth_dとして値を指定します。この項では未知の値であり、この値を求めるのが目的です。
 分散共分散行列$\boldsymbol{\Sigma} = (\sigma_{1,1}^2, \cdots, \sigma_{D,D}^2)$ではなく、標準偏差$\sigma_{i,i}$と相関係数$\sigma_{i,j}$の行列$(\sigma_{1,1}, \cdots, \sigma_{D,D})$をsigma_ddとして値を指定します。これは作図時に標準偏差を利用するためです。matrix()のデフォルトの仕様上、$(\sigma_{1,1}, \sigma_{1,2}, \sigma_{2,1}, \sigma_{2,2})$の順番に値を指定します。
 sigma_ddから精度パラメータ(精度行列)$\boldsymbol{\Lambda} = \boldsymbol{\Sigma}^{-1}$を計算してlambda_ddとします。逆行列はsolve()で計算します。

 グラフ用の点を作成します。

# 作図用のxの点を作成
x_1_vec <- seq(
  mu_truth_d[1] - 4 * sigma_dd[1, 1], 
  mu_truth_d[1] + 4 * sigma_dd[1, 1], 
  length.out = 1000
)
x_2_vec <- seq(
  mu_truth_d[2] - 4 * sigma_dd[2, 2], 
  mu_truth_d[2] + 4 * sigma_dd[2, 2], 
  length.out = 1000
)
x_point_mat <- cbind(
  rep(x_1_vec, times = length(x_2_vec)), 
  rep(x_2_vec, each = length(x_1_vec))
)

 作図用に、ガウス分布に従う変数$x_{n,1}$がとり得る値(x軸の値)をseq()で作成してx_1_vecとします。この例では、平均値を中心に標準偏差の4倍を範囲とします。length.out引数を使うと指定した要素数で等間隔に切り分けます。by引数を使うと切り分ける間隔を指定できます。処理が重い場合は、この値を調整してください。
 2次元方向(y軸の値)$x_{n,2}$についても同様に作成してx_2_vecとします。
 x_1_vecx_2_vecの要素の全ての組み合わせを持つようにマトリクスを作成してx_point_matとします。これは確率密度を等高線図にする際に、格子状の点(2軸の全ての値が直交する点)を渡す必要があるためです。

 尤度の確率密度を計算して、作図用のデータフレームを作成します。

# 尤度を計算:式(2.72)
model_df <- tibble(
  x_1 = x_point_mat[, 1], 
  x_2 = x_point_mat[, 2], 
  density = mvnfast::dmvn(
    X = x_point_mat, mu = mu_truth_d, sigma = solve(lambda_dd)
  )
)

 x_point_matの値(の組み合わせ)ごとに確率密度を計算します。多次元ガウス分布の確率密度は、mvnfast::dmvn()で計算できます。データの引数Xx_point_mat、平均の引数mumu_truth_d、分散共分散行列の引数sigmasolve(lambda_dd)を指定します。$\boldsymbol{\Sigma} = \boldsymbol{\Lambda}^{-1}$の計算をして精度行列から分散共分散行列に戻しています。sigma_dd^2を指定することもできます。

 作成したデータフレームを確認しましょう。

# 確認
head(model_df)
## # A tibble: 6 x 3
##     x_1   x_2       density
##   <dbl> <dbl>         <dbl>
## 1 -55     -70 0.00000000253
## 2 -54.8   -70 0.00000000259
## 3 -54.7   -70 0.00000000265
## 4 -54.5   -70 0.00000000271
## 5 -54.4   -70 0.00000000277
## 6 -54.2   -70 0.00000000284

 ggplot2パッケージを利用して作図するには、データフレームを渡す必要があります。

 尤度を作図します。

# 尤度を作図
ggplot(model_df, aes(x = x_1, y = x_2, z = density, color = ..level..)) + 
  geom_contour() + # 尤度
  labs(title = "Multivariate Gaussian Distribution", 
     subtitle = paste0("mu=(", paste(round(mu_truth_d, 1), collapse = ", "), ")", 
                       ", lambda=(", paste(round(lambda_dd, 5), collapse = ", "), ")"), 
     x = expression(x[1]), y = expression(x[2]), 
     color = "density")

f:id:anemptyarchive:20210410191446p:plain
尤度:多次元ガウス分布

 geom_contour()で、等高線グラフを描画します。格子状の点を渡す必要があります。

 真のパラメータを求めることは、この真の分布を求めることを意味します。

・データの生成

 続いて、構築したモデルに従って観測データ$\mathbf{X} = \{\mathbf{x}_1, \mathbf{x}_2, \cdots, \mathbf{x}_N\}$を生成します。

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

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

# 多次元ガウス分布に従うデータを生成
x_nd <- mvnfast::rmvn(n = N, mu = mu_truth_d, sigma = solve(lambda_dd))

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

 多次元ガウス分布に従う乱数は、mvnfast::rmvn()で生成できます。試行回数の引数nNを指定します。他の引数についてはmvnfast::dmvn()と同じです。生成したN個のデータをx_ndとします。

 観測したデータを確認しましょう。

# 確認
summary(x_nd)
##        V1               V2        
##  Min.   :-24.42   Min.   :-25.78  
##  1st Qu.: 16.90   1st Qu.: 30.59  
##  Median : 28.24   Median : 52.81  
##  Mean   : 28.07   Mean   : 52.14  
##  3rd Qu.: 37.76   3rd Qu.: 73.19  
##  Max.   : 74.33   Max.   :120.41


 観測データの散布図を尤度と重ねて確認します。

# 観測データのデータフレームを作成
x_df <- tibble(
  x_n1 = x_nd[, 1], 
  x_n2 = x_nd[, 2]
)

# 観測データの散布図を作図
ggplot() + 
  geom_point(data = x_df, aes(x = x_n1, y = x_n2)) + # 観測データ
  geom_contour(data = model_df, aes(x = x_1, y = x_2, z = density, color = ..level..)) + # 尤度
  labs(title = "Multivariate Gaussian Distribution", 
       subtitle = paste0("N=", N, ", mu=(", paste(mu_truth_d, collapse = ", "), ")", 
                         ", sigma=(", paste(sqrt(solve(lambda_dd)), collapse = ", "), ")"), 
      x = expression(x[1]), y = expression(x[2]), 
      color = "density")

f:id:anemptyarchive:20210410191512p:plain
観測データの散布図:多次元ガウス分布

 作図用に観測データをデータフレームに格納しておき、確率密度のグラフと重ねて表示します。

・事前分布の設定

 次に、尤度に対する共役事前分布を設定します。多次元ガウス分布の平均パラメータ$\boldsymbol{\mu}$に対する事前分布$p(\boldsymbol{\mu})$として、多次元ガウス分布$\mathcal{N}(\boldsymbol{\mu} | \mathbf{m}, \boldsymbol{\Lambda}_{\boldsymbol{\mu}}^{-1})$を設定します。

 $\boldsymbol{\mu}$の事前分布のパラメータ(超パラメータ)を設定します。

# muの事前分布のパラメータを指定
m_d <- c(0, 0)
sigma_mu_dd <- matrix(c(100, 0, 0, 100), nrow = 2, ncol = 2)
lambda_mu_dd <- solve(sigma_mu_dd^2)

 多次元ガウス分布の平均パラメータ$\mathbf{m} = (m_1, \cdots, m_D)$をm_d、精度パラメータ$\boldsymbol{\Lambda}_{\boldsymbol{\mu}} = (\lambda_{1,1}, \cdots, \lambda_{D,D})$をlambda_mu_ddとして、尤度と同様に値を指定します。

 $\boldsymbol{\mu}$の事前分布の確率密度を計算します。

# 作図用のmuの点を作成
mu_1_vec <- seq(mu_truth_d[1] - 100, mu_truth_d[1] + 100, length.out = 1000)
mu_2_vec <- seq(mu_truth_d[2] - 100, mu_truth_d[2] + 100, length.out = 1000)
mu_point_mat <- cbind(
  rep(mu_1_vec, times = length(mu_2_vec)), 
  rep(mu_2_vec, each = length(mu_1_vec))
)

# muの事前分布を計算:式(2.72)
prior_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(lambda_mu_dd)
  )
)

 尤度のときと同様に、多次元ガウス分布に従う変数$\boldsymbol{\mu}$がとり得る値を作成してmu_point_matとします。この例では、真の値を中心に指定した範囲とします(本当は自動で良い感じの範囲になるようにしたかった)。

 尤度のときと同様にして、確率密度を計算します。

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

# 確認
head(prior_df)
## # A tibble: 6 x 3
##    mu_1  mu_2   density
##   <dbl> <dbl>     <dbl>
## 1 -75     -50 0.0000106
## 2 -74.8   -50 0.0000106
## 3 -74.6   -50 0.0000106
## 4 -74.4   -50 0.0000106
## 5 -74.2   -50 0.0000107
## 6 -74.0   -50 0.0000107


 $\boldsymbol{\mu}$の事前分布を作図します。

# 作図用に真のmuのデータフレームを作成
mu_df <- tibble(
  mu_1 = mu_truth_d[1], 
  mu_2 = mu_truth_d[2]
)

# muの事前分布を作図
ggplot() + 
  geom_contour(data = prior_df, aes(x = mu_1, y = mu_2, z = density, color = ..level..)) + # muの事前分布
  geom_point(data = mu_df, aes(x = mu_1, y = mu_2), color = "red", shape = 4, size = 5) + # 真のmu
  labs(title = "Multivariate Gaussian Distribution", 
       subtitle = paste0("m=(", paste(round(m_d, 1), collapse = ", "), ")", 
                         ", lambda_mu=(", paste(round(lambda_mu_dd, 5), collapse = ", "), ")"), 
       x = expression(mu[1]), y = expression(mu[2]), 
       color = "density")

f:id:anemptyarchive:20210410191605p:plain
事前分布:多次元ガウス分布

 $\boldsymbol{\mu}$の真の値をデータフレームに格納しておき、確率密度のグラフと重ねて表示します。

・事後分布の計算

 観測データ$\mathbf{X}$から平均パラメータ$\boldsymbol{\mu}$の事後分布$p(\boldsymbol{\mu} | \mathbf{X})$を求めます(平均パラメータ$\boldsymbol{\mu}$を分布推定します)。事後分布は多次元ガウス分布$\mathcal{N}(\boldsymbol{\mu} | \hat{\mathbf{m}}, \hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1})$になります。

 観測データx_ndを用いて、$\boldsymbol{\mu}$の事後分布のパラメータを計算します。

# muの事後分布のパラメータを計算:式(3.102),(3.103)
lambda_mu_hat_dd <- N * lambda_dd + lambda_mu_dd
m_hat_d <- solve(lambda_mu_hat_dd) %*% (lambda_dd %*% colSums(x_nd) + lambda_mu_dd %*% m_d) %>% 
  as.vector()

 事後分布のパラメータを

$$ \begin{align} \hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}} &= N \boldsymbol{\Lambda} + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \tag{3.102}\\ \hat{\mathbf{m}} &= \hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1} \left( \boldsymbol{\Lambda} \sum_{n=1}^N \mathbf{x}_n + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} \right) \tag{3.103} \end{align} $$

で計算して、結果をlambda_mu_hat_ddm_hat_dとします。

# 確認
m_hat_d; lambda_mu_hat_dd
## [1] 28.02528 52.02970
##             [,1]        [,2]
## [1,]  0.14555455 -0.03636364
## [2,] -0.03636364  0.06474646


 求めたパラメータを用いて、$\boldsymbol{\mu}$の事後分布の確率密度を計算します。

# muの事後分布を計算:式(2.72)
posterior_df <- tibble(
  mu_1 = mu_point_mat[, 1], 
  mu_2 = mu_point_mat[, 2], 
  density = mvnfast::dmvn(
    X = , mu_point_mat, mu = m_hat_d, sigma = solve(lambda_mu_hat_dd)
  )
)

 更新した超パラメータm_hat_d, lambda_mu_hat_ddを用いて、事前分布のときと同様にして計算します。

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

# 確認
head(posterior_df)
## # A tibble: 6 x 3
##    mu_1  mu_2   density
##   <dbl> <dbl>     <dbl>
## 1 -75     -50 2.09e-318
## 2 -74.8   -50 2.00e-317
## 3 -74.6   -50 1.90e-316
## 4 -74.4   -50 1.79e-315
## 5 -74.2   -50 1.68e-314
## 6 -74.0   -50 1.57e-313


 $\boldsymbol{\mu}$の事後分布を作図します。

# muの事後分布を作図
ggplot() + 
  geom_contour(data = posterior_df, aes(x = mu_1, y = mu_2, z = density, color = ..level..)) + # muの事後分布
  geom_point(data = mu_df, aes(x = mu_1, y = mu_2), color = "red", shape = 4, size = 5) + # 真のmu
  labs(title = "Multivariate Gaussian Distribution", 
       subtitle = paste0("N=", N, ", m_hat=(", paste(round(m_hat_d, 1), collapse = ", "), ")", 
                         ", sigma_mu_hat=(", paste(round(sqrt(solve(lambda_mu_hat_dd)), 1), collapse = ", "), ")"), 
       x = expression(mu[1]), y = expression(mu[2]), 
       color = "density")

f:id:anemptyarchive:20210410191638p:plain
事後分布:多次元ガウス分布

 $\boldsymbol{\mu}$の真の値付近をピークとする分布を推定できています。

 真の値を等高線の中心に捉えていないように見えるかもしれませんが、x軸とy軸の値に注目してください。事前分布と比べて範囲が狭くなっています。グラフの表示範囲を広げてみましょう。

# muの事後分布を作図
ggplot() + 
  geom_contour(data = posterior_df, aes(x = mu_1, y = mu_2, z = density, color = ..level..)) + # muの事後分布
  geom_point(data = mu_df, aes(x = mu_1, y = mu_2), color = "red", shape = 4, size = 5) + # 真のmu
  xlim(c(min(mu_point_mat[, 1]), max(mu_point_mat[, 1]))) + # x軸の表示範囲
  ylim(c(min(mu_point_mat[, 2]), max(mu_point_mat[, 2]))) + # y軸の表示範囲
  labs(title = "Multivariate Gaussian Distribution", 
       subtitle = paste0("N=", N, ", m_hat=(", paste(round(m_hat_d, 1), collapse = ", "), ")", 
                         ", sigma_mu_hat=(", paste(round(sqrt(solve(lambda_mu_hat_dd)), 1), collapse = ", "), ")"), 
       x = expression(mu[1]), y = expression(mu[2]), 
       color = "density")

f:id:anemptyarchive:20210410191657p:plain
事後分布:多次元ガウス分布

 先のとがった分布になっているのが分かります。

・予測分布の計算

 最後に、観測データ$\mathbf{X}$から未観測のデータ$\mathbf{x}_{*}$の予測分布$p(\mathbf{x}_{*} | \mathbf{X})$を求めます。予測分布は多次元ガウス分布$\mathcal{N}(\mathbf{x}_{*} | \hat{\boldsymbol{\mu}}_{*}, \hat{\boldsymbol{\Lambda}}_{*}^{-1})$になります。

 $\boldsymbol{\mu}$の事後分布のパラメータ、または観測データと$\boldsymbol{\mu}$の事前分布のパラメータを用いて、予測分布のパラメータを計算します。

# 予測分布のパラメータを計算:式(3.109'),(3.110')
lambda_star_hat_dd <- solve(solve(lambda_dd) + solve(lambda_mu_hat_dd))
mu_star_hat_d <- m_hat_d

 予測分布のパラメータを

$$ \begin{aligned} \hat{\boldsymbol{\Lambda}}_{*} &= ( \boldsymbol{\Lambda}^{-1} + \hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1} )^{-1} \\ &= \{ \boldsymbol{\Lambda}^{-1} + ( N \boldsymbol{\Lambda} + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} )^{-1} \}^{-1} \\ \hat{\boldsymbol{\mu}}_{*} &= \hat{\mathbf{m}} \\ &= \hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1} \left( \boldsymbol{\Lambda} \sum_{n=1}^N \mathbf{x}_n + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} \right) \end{aligned} $$

で計算して、結果をlambda_star_hat_ddm_star_hat_dとします。

 それぞれ上の式だと、事後分布のパラメータm_hat_d, lambda_mu_hat_ddで計算できます。下の式だと、観測データx_ndと事前分布のパラメータm_d, lambda_mu_ddで計算できます。

# 確認
mu_star_hat_d; lambda_star_hat_dd
## [1] 28.02528 52.02970
##               [,1]          [,2]
## [1,]  0.0028520883 -0.0007130125
## [2,] -0.0007130125  0.0012676161

 $\mathbf{X}$から$\hat{\boldsymbol{\mu}}_{*},\ \hat{\boldsymbol{\Lambda}}_{*}$を学習しているのが式からも分かります。

 求めたパラメータを用いて、予測分布の確率密度を計算します。

# 予測分布を計算:式(2.72)
predict_df <- tibble(
  x_1 = x_point_mat[, 1], 
  x_2 = x_point_mat[, 2], 
  density = mvnfast::dmvn(
    X = x_point_mat, mu = mu_star_hat_d, sigma = solve(lambda_star_hat_dd)
  )
)

 尤度のときと同様に、x_point_matの値ごとに確率密度を計算します。

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

# 確認
head(predict_df)
## # A tibble: 6 x 3
##     x_1   x_2       density
##   <dbl> <dbl>         <dbl>
## 1 -55     -70 0.00000000165
## 2 -54.8   -70 0.00000000169
## 3 -54.7   -70 0.00000000173
## 4 -54.5   -70 0.00000000177
## 5 -54.4   -70 0.00000000181
## 6 -54.2   -70 0.00000000186


 予測分布を尤度と重ねて作図します。

# 予測分布を作図
ggplot() + 
  geom_contour(data = predict_df, aes(x = x_1, y = x_2, z = density, color = ..level..)) + # 予測分布
  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 = mu_1, y = mu_2), color = "red", shape = 4, size = 5) + # 真のmu
  #geom_point(data = x_df, aes(x = x_n1, y = x_n2)) + # 観測データ
  labs(title = "Multivariate Gaussian Distribution", 
       subtitle = paste0("N=", N, ", mu_star_hat=(", paste(round(mu_star_hat_d, 1), collapse = ", "), ")", 
                         ", sigma_star_hat=(", paste(round(sqrt(solve(lambda_star_hat_dd)), 1), collapse = ", "), ")"), 
       x = expression(x[1]), y = expression(x[2]), 
       color = "density")

f:id:anemptyarchive:20210410191716p:plain
予測分布:多次元ガウス分布

 観測データが増えると、予測分布が真の分布に近づきます。(なぜか精度から標準偏差$\sigma_{d,d}$にしていますが深い意味はありません。)

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

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

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

 異なる点のみを簡単に解説します。

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


・モデルの設定

# 真のパラメータを指定
mu_truth_d <- c(25, 50)
sigma_dd <- matrix(c(20, 15, 15, 30), nrow = 2, ncol = 2)
lambda_dd <- solve(sigma_dd^2)

# muの事前分布のパラメータを指定
m_d <- c(0, 0)
sigma_mu_dd <- matrix(c(100, 0, 0, 100), nrow = 2, ncol = 2)
lambda_mu_dd <- solve(sigma_mu_dd^2)

# 作図用のmuの点を作成
mu_1_vec <- seq(mu_truth_d[1] - 100, mu_truth_d[1] + 100, length.out = 500)
mu_2_vec <- seq(mu_truth_d[2] - 100, mu_truth_d[2] + 100, length.out = 500)
mu_point_mat <- cbind(
  rep(mu_1_vec, times = length(mu_2_vec)), 
  rep(mu_2_vec, each = length(mu_1_vec))
)

# muの事前分布を計算:式(2.72)
posterior_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(lambda_mu_dd)
  ), 
  label = as.factor(
    paste0(
      "N=", 0, 
      ", m=(", paste(round(m_d, 1), collapse = ", "), ")", 
      ", lambda_mu=(", paste(round(lambda_mu_dd, 5), collapse = ", "), ")"
    )
  )
)

# 初期値による予測分布のパラメータを計算:式(3.109),(3.110)
lambda_star_dd <- solve(solve(lambda_dd) + solve(lambda_mu_dd))
mu_star_d <- m_d

# 作図用のxの点を作成
x_1_vec <- seq(mu_truth_d[1] - 4 * sigma_dd[1, 1], mu_truth_d[1] + 4 * sigma_dd[1, 1], length.out = 500)
x_2_vec <- seq(mu_truth_d[2] - 4 * sigma_dd[2, 2], mu_truth_d[2] + 4 * sigma_dd[2, 2], length.out = 500)
x_point_mat <- cbind(
  rep(x_1_vec, times = length(x_2_vec)), 
  rep(x_2_vec, each = length(x_1_vec))
)

# 予測分布を計算:式(2.72)
predict_df <- tibble(
  x_1 = x_point_mat[, 1], 
  x_2 = x_point_mat[, 2], 
  density = mvnfast::dmvn(
    X = x_point_mat, mu = mu_star_d, sigma = solve(lambda_star_dd)
  ), 
  label = as.factor(
    paste0(
      "N=", 0, 
      ", mu_star=(", paste(round(mu_star_d, 1), collapse = ", "), ")", 
      ", lambda_star=(", paste(round(lambda_star_dd, 5), collapse = ", "), ")"
    )
  )
)

 各試行の結果を同じデータフレームに格納していく必要があります。事後分布をposterior_df、予測分布をpredict_dfとして、初期値の結果を持つように作成しておきます。

・推論処理

# データ数(試行回数)を指定
N <- 100

# 観測データの受け皿を初期化
x_nd <- matrix(0, nrow = N, ncol = 2)

# ベイズ推論
for(n in 1:N) {
  
  # 多次元ガウス分布に従うデータを生成
  x_nd[n, ] <- mvnfast::rmvn(n = 1, mu = mu_truth_d, sigma = solve(lambda_dd)) %>% 
    as.vector()
  
  # 事後分布のパラメータを更新:式(3.102),(3.103)
  old_lambda_mu_dd <- lambda_mu_dd
  lambda_mu_dd <- lambda_dd + lambda_mu_dd
  m_d <- solve(lambda_mu_dd) %*% (lambda_dd %*% x_nd[n, ] + old_lambda_mu_dd %*% m_d) %>% 
    as.vector()
  
  # muの事後分布を計算:式(2.72)
  tmp_posterior_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(lambda_mu_dd)
    ), 
    label = as.factor(
      paste0(
        "N=", n, 
        ", m_hat=(", paste(round(m_d, 1), collapse = ", "), ")", 
        ", lambda_mu_hat=(", paste(round(lambda_mu_dd, 5), collapse = ", "), ")"
      )
    )
  )
  
  # 予測分布のパラメータを更新:式(3.109),(3.110)
  lambda_star_dd <- solve(solve(lambda_dd) + solve(lambda_mu_dd))
  mu_star_d <- m_d
  
  # 予測分布を計算:式(2.72)
  tmp_predict_df <- tibble(
    x_1 = x_point_mat[, 1], 
    x_2 = x_point_mat[, 2], 
    density = mvnfast::dmvn(
      X = x_point_mat, mu = mu_star_d, sigma = solve(lambda_star_dd)
    ), 
    label = as.factor(
      paste0(
        "N=", n, 
        ", mu_star_hat=(", paste(round(mu_star_d, 1), collapse = ", "), ")", 
        ", lambda_star_hat=(", paste(round(lambda_star_dd, 5), collapse = ", "), ")"
      )
    )
  )
  
  # 推論結果を結合
  posterior_df <- rbind(posterior_df, tmp_posterior_df)
  predict_df <- rbind(predict_df, tmp_predict_df)
  
  # 動作確認
  #print(paste0("n=", n, " (", round(n / N * 100, 1), "%)"))
}

 観測された各データによってどのように学習する(分布が変化する)のかを確認するため、for()で1データずつ処理します。よって、データ数Nがイタレーション数になります。

 一度の処理で事後分布のパラメータを計算するのではなく、事前分布(1ステップ前の事後分布)に対して繰り返し観測データの情報を与えることでパラメータを更新(上書き)していきます。
 それに伴い、事後分布のパラメータの計算式(3.102-103)の$N$と$\sum_{n=1}^N \mathbf{x}_n$の計算は、ループ処理によって$N$回繰り返し1x_nd[n, ]を加えることで行います。$n$回目のループ処理のときには、$n-1$回分の1x_nd[n, ]が既にm_dlambda_mu_ddに加えられているわけです。
 ただし、事後分布のパラメータの計算において、更新前(1ステップ前)のパラメータを使うため、old_lambda_mu_ddとして値を一時的に保存しておきます。

・事後分布の推移

# 作図用に真のmuのデータフレームを作成
mu_df <- tibble(
  mu_1 = mu_truth_d[1], 
  mu_2 = mu_truth_d[2]
)

# 事後分布を作図
posterior_graph <- ggplot() + 
  geom_contour(data = posterior_df, aes(x = mu_1, y = mu_2, z = density, color = ..level..)) + # muの事後分布
  geom_point(data = mu_df, aes(x = mu_1, y = mu_2), color = "red", shape = 4, size = 5) + # 真のmu
  gganimate::transition_manual(label) + # フレーム
  labs(title = "Multivariate Gaussian Distribution", 
       subtitle = "{current_frame}", 
       x = expression(mu[1]), y = expression(mu[2]), 
       color = "density")

# gif画像を作成
gganimate::animate(posterior_graph, nframes = N + 1, fps = 10)

 各フレームの順番を示す列(label)をgganimate::transition_manual()に指定します。分布の推移と共にパラメータの値を表示するようにlabel列を作成していますが、ややこしければas.factor(paste0("iter=", n))として試行回数だけ表示するだけでもそれっぽくなります。

 初期値(事前分布)を含むため、フレーム数の引数nframesN + 1です。

・予測分布の推移

# 尤度を計算:式(2.72)
model_df <- tibble(
  x_1 = x_point_mat[, 1], 
  x_2 = x_point_mat[, 2], 
  density = mvnfast::dmvn(
    X = x_point_mat, mu = mu_truth_d, sigma = solve(lambda_dd)
  )
)

# 作図用に観測データのデータフレームを作成
label_vec <- unique(predict_df[["label"]]) # ラベルを抽出
x_df <- tibble(x_n1 = NA, x_n2 = NA, label = label_vec[1]) # 初期値用
for(n in 1:N) {
  # n個目までのデータフレームを作成
  tmp_x_df <- tibble(
    x_n1 = x_nd[1:n, 1], 
    x_n2 = x_nd[1:n, 2], 
    label = label_vec[n + 1]
  )
  
  # 結合
  x_df <- rbind(x_df, tmp_x_df)
}

# 予測分布を作図
predict_graph <- ggplot() + 
  geom_contour(data = predict_df, aes(x = x_1, y = x_2, z = density, color = ..level..)) + # 予測分布
  geom_contour(data = model_df, aes(x = x_1, y = x_2, z = density, color = ..level..), 
               alpha = 0.5, linetype = "dashed") + # 真の分布
  geom_point(data = x_df, aes(x = x_n1, y = x_n2)) + # 観測データ
  gganimate::transition_manual(label) + # フレーム
  labs(title = "Multivariate Gaussian Distribution", 
       subtitle = "{current_frame}", 
       x = expression(x[1]), y = expression(x[2]), 
       color = "density")

# gif画像を作成
gganimate::animate(predict_graph, nframes = N + 1, fps = 10)

 n番目のフレームのグラフにおいて、n回目までに生成(観測)されたデータ$\mathbf{x}_1, \cdots, \mathbf{x}_n$を表示するために、Nフレーム分までを順番に取り出してx_dfに格納していきます。このときlabel列がpredict_dfと一致している必要があります。


f:id:anemptyarchive:20210410233413g:plain
事後分布の推移:多次元ガウス分布

 (始めの方は精度が小さい(分散が大きい)せいで確率密度が低くなっているため等高線が表示されていません。)

f:id:anemptyarchive:20210410233545g:plain
予測分布の推移:多次元ガウス分布

 新たなデータによって平均(分布の中心)が推移しているのを確認できます。

参考文献

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

おわりに

 その時々では理解したつもりで記事にしているのですが、勉強を続けていると勘違いしていたことに気付くことも多々あります。適宜修正していきたいですけど結構大変。独学あるある?

【次節の内容】

明日更新する!

  • 2021/04/10:加筆修正しました。その際に数式読解編と記事を分割しました。