からっぽのしょこ

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

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

はじめに

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

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

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

【数式読解編】

www.anarchive-beta.com

【他の節の内容】

www.anarchive-beta.com

【この節の内容】

・Rでやってみよう

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

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

# 3.4.2項で利用するパッケージ
library(tidyverse)
library(mvnfast)

 mvnfastパッケージは、多次元ガウス分布に関するパッケージです。多次元ガウス分布に従う乱数生成関数rmvn()と確率密度関数dmvn()を使います。また予測分布の計算において、多次元スチューデントのt分布の確率密度関数dmvt()も使います。

・観測モデルの設定

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

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

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

 平均パラメータ$\boldsymbol{\mu} = (\mu_1, \cdots, \mu_D)$をmu_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_truth_ddとして値を指定します。これは作図時に標準偏差を利用するためです。matrix()のデフォルトの仕様上、$(\sigma_{1,1}, \sigma_{1,2}, \sigma_{2,1}, \sigma_{2,2})$の順番に値を指定します。
 sigma_truth_ddから精度パラメータ(精度行列)$\boldsymbol{\Lambda} = \boldsymbol{\Sigma}^{-1}$を計算してlambda_truth_ddとします。逆行列はsolve()で計算します。この項では未知の値であり、この値を求めるのが目的です。

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

# 作図用のxの点を作成
x_1_vec <- seq(
  mu_d[1] - 4 * sigma_truth_dd[1, 1], 
  mu_d[1] + 4 * sigma_truth_dd[1, 1], 
  length.out = 1000
)
x_2_vec <- seq(
  mu_d[2] - 4 * sigma_truth_dd[2, 2], 
  mu_d[2] + 4 * sigma_truth_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次元方向$x_{n,2}$(y軸の値)についても同様に作成して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_d, sigma = solve(lambda_truth_dd)
  )
)

 x_point_matの値(の組み合わせ)ごとに確率密度を計算します。多次元ガウス分布の確率密度は、mvnfast::dmvn()で計算できます。データの引数Xx_point_mat、平均の引数mumu_d、分散共分散行列の引数sigmasolve(lambda_truth_dd)を指定します。$\boldsymbol{\Sigma} = \boldsymbol{\Lambda}^{-1}$の計算をして精度行列から分散共分散行列に戻しています。sigma_truth_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_d, 1), collapse = ", "), ")", 
                         ", lambda=(", paste(round(lambda_truth_dd, 5), collapse = ", "), ")"), 
       x = expression(x[1]), y = expression(x[2]), 
       color = "density")

f:id:anemptyarchive:20210410221307p: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_d, sigma = solve(lambda_truth_dd))

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

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

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

# 確認
summary(x_nd)
##        V1                V2         
##  Min.   :-21.275   Min.   :  1.202  
##  1st Qu.:  9.939   1st Qu.: 26.863  
##  Median : 19.116   Median : 46.124  
##  Mean   : 19.577   Mean   : 46.538  
##  3rd Qu.: 31.175   3rd Qu.: 61.092  
##  Max.   : 48.411   Max.   :101.771


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

# 観測データのデータフレームを作成
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_d, collapse = ", "), ")", 
                         ", sigma=(", paste(sqrt(solve(lambda_truth_dd)), collapse = ", "), ")"), 
       x = expression(x[1]), y = expression(x[2]), 
       color = "density")

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

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

・事前分布の設定

 次に、尤度に対する共役事前分布を設定します。ガウス分布の精度パラメータ$\boldsymbol{\Lambda}$に対する事前分布$p(\boldsymbol{\Lambda})$として、多次元ガウス分布$\mathcal{W}(\boldsymbol{\Lambda} | \nu, \mathbf{W})$を設定します。

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

# lambdaの事前分布のパラメータを指定
w_dd <- matrix(c(0.0005, 0, 0, 0.0005), nrow = 2, ncol = 2)
nu <- 2

 ウィシャート分布のパラメータ$\mathbf{W} = (w_{1,1}, \cdots, w_{D,D})$をW_dd、自由度$\nu$をnuとして値を指定します。$\mathbf{W}$は正定値行列、$\nu$は$D - 1$より大きい必要があります。

 これまでのように、このパラメータを用いて事前分布(の確率密度)を計算しグラフで確認したいところですが、精度行列の分布はイメージしにくいため、精度パラメータの期待値$\mathbb{E}[\boldsymbol{\Lambda}]$を用いた多次元ガウス分布を尤度と比較することにしましょう。ちなみに、ウィシャート分布の確率密度はMCMCpack::dwish()で計算できます(たぶん)。

# lambdaの期待値を計算:式(2.89)
E_lambda_dd <- nu * w_dd

# 確認
E_lambda_dd
##       [,1]  [,2]
## [1,] 0.001 0.000
## [2,] 0.000 0.001

 精度行列の期待値は、ウィシャート分布の期待値(2.89)より

$$ \mathbb{E}[\boldsymbol{\Lambda}] = \nu \mathbf{W} $$

で計算できます。

 $\boldsymbol{\Lambda}$の事前分布の期待値(精度パラメータの期待値)を用いて、多次元ガウス分布の確率密度を計算します。

# 事前分布の期待値を用いた分布を計算:式(2.72)
prior_df <- tibble(
  x_1 = x_point_mat[, 1], 
  x_2 = x_point_mat[, 2], 
  density = mvnfast::dmvn(
    X = x_point_mat, mu = mu_d, sigma = solve(E_lambda_dd)
  )
)

 mu_d, E_lambda_ddを用いて、尤度のときと同様にして計算します。

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

# 確認
head(prior_df)
## # A tibble: 6 x 3
##     x_1   x_2       density
##   <dbl> <dbl>         <dbl>
## 1 -55     -70 0.00000000484
## 2 -54.8   -70 0.00000000491
## 3 -54.7   -70 0.00000000497
## 4 -54.5   -70 0.00000000503
## 5 -54.4   -70 0.00000000510
## 6 -54.2   -70 0.00000000516


 $\boldsymbol{\Lambda}$の事前分布の期待値(精度パラメータの期待値)を用いた分布を作図します。

# 事前分布の期待値を用いた分布を作図
ggplot() + 
  geom_contour(data = prior_df, aes(x = x_1, y = x_2, z = density, color = ..level..)) + # lambdaの期待値を用いた分布
  geom_contour(data = model_df, aes(x = x_1, y = x_2, z = density, color = ..level..), 
               alpha = 0.5, linetype = "dashed") + # 尤度
  labs(title = "Multivariate Gaussian Distribution", 
       subtitle = paste0("nu=", nu, 
                         ", W=(", paste(round(w_dd, 5), collapse = ", "), ")"), 
       x = expression(x[1]), y = expression(x[2]), 
       color = "density")

f:id:anemptyarchive:20210410221356p:plain
事前分布の期待値を用いた分布:多次元ガウス分布

 尤度を破線で重ねて描画します。

・事後分布の計算

 観測データ$\mathbf{X}$から精度パラメータ$\boldsymbol{\Lambda}$の事後分布$p(\boldsymbol{\Lambda} | \mathbf{X})$を求めます(精度パラメータ$\boldsymbol{\Lambda}$を分布推定します)。事後分布はウィシャート分布$p(\boldsymbol{\Lambda} | \hat{\nu}, \hat{\mathbf{W}})$になります。

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

# lambdaの事後分布のパラメータを計算:式(3.116)
w_hat_dd <- solve(
  (t(x_nd) - mu_d) %*% t(t(x_nd) - mu_d) + solve(w_dd)
)
nu_hat <- N + nu

 事後分布のパラメータを

$$ \begin{aligned} \hat{\mathbf{W}} &= \sum_{n=1}^N (\mathbf{x}_n - \mu) (\mathbf{x}_n - \mu)^{\top} + \mathbf{W}^{-1} \\ \hat{\nu} &= N + \nu \end{aligned} \tag{3.116} $$

で計算して、結果をw_hat_ddnu_hatとします。ただし、$\sum_{n=1}^N (\mathbf{x}_n - \mu) (\mathbf{x}_n - \mu)^{\top}$の計算を効率よく処理するために、$\tilde{x}_{n,d} = x_{n,d} - \mu_d$として$\tilde{\mathbf{X}}^{\top} \tilde{\mathbf{X}}$で計算しています。

# 確認
w_hat_dd; nu_hat
##               [,1]          [,2]
## [1,]  6.602907e-05 -1.767267e-05
## [2,] -1.767267e-05  3.246292e-05
## [1] 52


 事前分布のときと同様に、求めたパラメータを用いて、学習後の精度パラメータ$\hat{\boldsymbol{\Lambda}}$の期待値$\mathbb{E}[\hat{\boldsymbol{\Lambda}}]$を計算します。

# lambdaの期待値を計算:式(2.89)
E_lambda_hat_dd <- nu_hat * w_hat_dd

# 確認
E_lambda_hat_dd
##              [,1]         [,2]
## [1,]  0.003433511 -0.000918979
## [2,] -0.000918979  0.001688072


 $\boldsymbol{\Lambda}$の事後分布の期待値(精度パラメータの期待値)を用いた多次元ガウス分布の確率密度を計算します。

# 事後分布の期待値を用いた分布を計算:式(2.72)
posterior_df <- tibble(
  x_1 = x_point_mat[, 1], 
  x_2 = x_point_mat[, 2], 
  density = mvnfast::dmvn(
    X = x_point_mat, mu = mu_d, sigma = solve(E_lambda_hat_dd)
  )
)

 mu_d, E_lambda_hat_ddを用いて、尤度のときと同様にして計算します。

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

# 確認
head(posterior_df)
## # A tibble: 6 x 3
##     x_1   x_2  density
##   <dbl> <dbl>    <dbl>
## 1 -55     -70 2.14e-10
## 2 -54.8   -70 2.20e-10
## 3 -54.7   -70 2.26e-10
## 4 -54.5   -70 2.32e-10
## 5 -54.4   -70 2.38e-10
## 6 -54.2   -70 2.44e-10


 $\boldsymbol{\Lambda}$の事後分布の期待値(精度パラメータの期待値)を用いた分布を作図します。

# 事後分布の期待値を用いた分布を作図
ggplot() + 
  geom_contour(data = posterior_df, aes(x = x_1, y = x_2, z = density, color = ..level..)) + # lambdaの期待値を用いた分布
  geom_contour(data = model_df, aes(x = x_1, y = x_2, z = density, color = ..level..), 
               alpha = 0.5, linetype = "dashed") + # 尤度
  labs(title = "Multivariate Gaussian Distribution", 
       subtitle = paste0("N=", N, ", nu_hat=", nu_hat, 
                         ", W_hat=(", paste(round(w_hat_dd, 5), collapse = ", "), ")"), 
       x = expression(x[1]), y = expression(x[2]), 
       color = "density")

f:id:anemptyarchive:20210410221428p:plain
事後分布の期待値を用いた分布:多次元ガウス分布

 分布の形状が尤度に近づいているのを確認できます。

・予測分布の計算

 最後に、観測データ$\mathbf{X}$から未観測のデータ$\mathbf{x}_{*}$の予測分布$p(\mathbf{x}_{*} | \mathbf{X})$を求めます。予測分布は多次元スチューデントのt分布$p(\mathbf{x}_{*} | \boldsymbol{\mu}_s, \hat{\boldsymbol{\Lambda}}_s, \hat{\nu}_s)$になります。

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

# 次元数を取得
D <- length(mu_d)

# 予測分布のパラメータを計算:式(3.124)
mu_s_d <- mu_d
lambda_s_hat_dd <- (1 - D + nu_hat) * w_hat_dd
nu_s_hat <- 1 - D + nu_hat

 予測分布のパラメータを

$$ \begin{aligned} \boldsymbol{\mu}_s &= \boldsymbol{\mu} \\ \hat{\boldsymbol{\Lambda}}_s &= (1 - D + \hat{\nu}) \hat{\mathbf{W}} \\ &= (1 - D + N + \nu) \left\{ \sum_{n=1}^N (\mathbf{x}_n - \mu) (\mathbf{x}_n - \mu)^{\top} + \mathbf{W}^{-1} \right\} \\ \hat{\nu}_s &= 1 - D + \hat{\nu} \\ &= 1 - D + N + \nu \end{aligned} $$

で計算して、結果をmu_s_d, lambda_s_hat_ddm_s_hat_dとします。

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

# 確認
mu_s_d; lambda_s_hat_dd; nu_s_hat
## [1] 25 50
##               [,1]          [,2]
## [1,]  0.0033674824 -0.0009013063
## [2,] -0.0009013063  0.0016556091
## [1] 51

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

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

# 予測分布を計算:式(3.121)
predict_df <- tibble(
  x_1 = x_point_mat[, 1], 
  x_2 = x_point_mat[, 2], 
  density = mvnfast::dmvt(
    X = x_point_mat, mu = mu_s_d, sigma = solve(lambda_s_hat_dd), df = nu_s_hat
  )
)

 尤度のときと同様に、x_point_matの値ごとに確率密度を計算します。多次元のスチューデントのt分布の確率密度は、mvnfast::dmvt()で計算できます。データの引数Xにはx_point_mat、平均の引数mumu_s_d、スケール行列?の引数sigmalambda_s_hat_ddの逆行列、自由度の引数dfnu_s_hatを指定します。

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

# 確認
head(predict_df)
## # A tibble: 6 x 3
##     x_1   x_2       density
##   <dbl> <dbl>         <dbl>
## 1 -55     -70 0.00000000310
## 2 -54.8   -70 0.00000000316
## 3 -54.7   -70 0.00000000321
## 4 -54.5   -70 0.00000000327
## 5 -54.4   -70 0.00000000332
## 6 -54.2   -70 0.00000000338


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

# 予測分布を作図
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)) + # 観測データ
  labs(title = "Multivariate Student's t Distribution", 
       subtitle = paste0("N=", N, ", nu_s_hat=", nu_s_hat, 
                         ", lambda_s_hat=(", paste(round(lambda_s_hat_dd, 5), collapse = ", "), ")"), 
       x = expression(x[1]), y = expression(x[2]), 
       color = "density")

f:id:anemptyarchive:20210410221457p:plain
予測分布:多次元スチューデントのt分布

 観測データが増えると、予測分布が真の分布に近づきます。

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

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

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

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

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


・モデルの設定

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

# 事前分布のパラメータを指定
nu <- 2
w_dd <- matrix(c(0.0005, 0, 0, 0.0005), nrow = 2, ncol = 2)

# 初期値によるlambdaの期待値を計算:式(2.89)
E_lambda_dd <- nu * w_dd

# 初期値による予測分布のパラメータを計算:式(3.124)
mu_s_d <- mu_d
lambda_s_dd <- (nu - 1) * w_dd
nu_s <- nu - 1


# 作図用のxの点を作成
x_1_vec <- seq(mu_d[1] - 4 * sigma_truth_dd[1, 1], mu_d[1] + 4 * sigma_truth_dd[1, 1], length.out = 500)
x_2_vec <- seq(mu_d[2] - 4 * sigma_truth_dd[2, 2], mu_d[2] + 4 * sigma_truth_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)
posterior_df <- tibble(
  x_1 = x_point_mat[, 1], 
  x_2 = x_point_mat[, 2], 
  density = mvnfast::dmvn(
    X = x_point_mat, mu = mu_d, sigma = solve(E_lambda_dd)
  ), 
  label = as.factor(
    paste0(
      "N=", 0, ", nu=", nu, 
      ", W=(", paste(round(w_dd, 5), collapse = ", "), ")"
    )
  )
)

# 予測分布を計算:式(3.121)
predict_df <- tibble(
  x_1 = x_point_mat[, 1], 
  x_2 = x_point_mat[, 2], 
  density = mvnfast::dmvt(
    X = x_point_mat, mu = mu_s_d, sigma = solve(lambda_s_dd), df = nu_s
  ), 
  label = as.factor(
    paste0(
      "N=", 0, ", nu_s=", nu_s, 
      ", lambda_s=(", paste(round(lambda_s_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_d, sigma = solve(lambda_truth_dd)) %>% 
    as.vector()
  
  # 事後分布のパラメータを更新:式(3.116)
  w_dd <- solve(
    (x_nd[n, ] - mu_d) %*% t(x_nd[n, ] - mu_d) + solve(w_dd)
  )
  nu <- 1 + nu
  
  # lambdaの期待値を計算:式(2.89)
  E_lambda_dd <- nu * w_dd
  
  # 事後分布の期待値を用いた分布を更新:式(2.72)
  tmp_posterior_df <- tibble(
    x_1 = x_point_mat[, 1], 
    x_2 = x_point_mat[, 2], 
    density = mvnfast::dmvn(
      X = x_point_mat, mu = mu_d, sigma = solve(E_lambda_dd)
    ), 
    label = as.factor(
      paste0(
        "N=", n, ", nu_hat=", nu, 
        ", W_hat=(", paste(round(w_dd, 5), collapse = ", "), ")"
      )
    )
  )
  
  # 予測分布のパラメータを更新:式(3.124)
  #mu_s_d <- mu_d
  lambda_s_dd <- (nu - 1) * w_dd
  nu_s <- nu - 1
  
  # 予測分布を計算:式(3.121)
  tmp_predict_df <- tibble(
    x_1 = x_point_mat[, 1], 
    x_2 = x_point_mat[, 2], 
    density = mvnfast::dmvt(
      X = x_point_mat, mu = mu_s_d, sigma = solve(lambda_s_dd), df = nu_s
    ), 
    label = as.factor(
      paste0(
        "N=", n, ", nu_s_hat=", nu_s, 
        ", lambda_s_hat=(", paste(round(lambda_s_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.116)の$\sum_{n=1}^N \mathbf{x}_n$と$N$の計算は、ループ処理によって$N$回繰り返しx_nd[n, ]1を加えることで行います。$n$回目のループ処理のときには、$n-1$回分のx_nd[n, ]1が既にw_ddnuに加えられているわけです。

・事後分布の推移

# 尤度を計算:式(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_d, sigma = solve(lambda_truth_dd)
  )
)


# 作図用に観測データのデータフレームを作成
label_vec <- unique(posterior_df[["label"]]) # lambdaの期待値を用いた分布のラベルを抽出
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)
}

# 事後分布の期待値を用いた分布を作図
posterior_graph <- ggplot() + 
  geom_contour(data = posterior_df, aes(x = x_1, y = x_2, z = density, color = ..level..)) + # lambdaの期待値を用いた分布
  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)) + # 観測データ
  transition_manual(label) + # フレーム
  labs(title = "Multivariate Gaussian Distribution", 
       subtitle = "{current_frame}", 
       x = expression(x[1]), y = expression(x[2]), 
       color = "density")

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

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

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

・予測分布の推移

# 作図用に観測データのデータフレームを作成
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)) + # 観測データ
  transition_manual(label) + # フレーム
  labs(title = "Multivariate Student's t Distribution", 
       subtitle = "{current_frame}", 
       x = expression(x[1]), y = expression(x[2]), 
       color = "density")

# gif画像を作成
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:20210410221551g:plain
事後分布の期待値を用いた分布の推移:多次元ガウス分布

f:id:anemptyarchive:20210410221609g:plain
予測分布の推移:多次元スチューデントのt分布

 新たなデータによって精度(分布の広がり)が推移しているのを確認できます。

参考文献

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

おわりに

 mvnfastパッケージにも多次元ガウス分布関連の関数がありますが、書き換えるのが面倒だったので止めました。mvtnormパッケージからmvnfastパッケージに置き替えました。

【次節の内容】

明日更新予定!


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