からっぽのしょこ

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

【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)

# 真の分散共分散行列を指定
sigma2_truth_dd <- matrix(c(600, 90, 90, 400), nrow = 2, ncol = 2)

# 真の精度行列を計算
lambda_truth_dd <- solve(sigma2_truth_dd)

 平均パラメータ$\boldsymbol{\mu} = (\mu_1, \cdots, \mu_D)$をmu_dとして値を指定します。
 分散共分散行列$\boldsymbol{\Sigma} = (\sigma_{1,1}^2, \cdots, \sigma_{D,D}^2)$をsigma2_truth_ddとして値を指定します。作図時に標準偏差に変換して利用します。matrix()のデフォルトの仕様上、$(\sigma_{1,1}^2, \sigma_{2,1}^2, \sigma_{1,2}^2, \sigma_{2,2}^2)$の順番に値を指定します。$\sigma_{i,i}^2$は$i$次元方向の分散、$\sigma_{i,j}^2$は$i$次元と$j$次元方向の共分散です。
 分散共分散行列の逆行列を精度行列と呼びます。sigma2_truth_ddから精度パラメータ(精度行列)$\boldsymbol{\Lambda} = \boldsymbol{\Sigma}^{-1}$を計算してlambda_truth_ddとします。逆行列はsolve()で計算できます。この例では未知の値であり、この値を求めるのが目的です。

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

# 作図用のxのx軸の値を作成
x_1_vec <- seq(
  mu_d[1] - 3 * sqrt(sigma2_truth_dd[1, 1]), 
  mu_d[1] + 3 * sqrt(sigma2_truth_dd[1, 1]), 
  length.out = 500
)

# 作図用のxのy軸の値を作成
x_2_vec <- seq(
  mu_d[2] - 3 * sqrt(sigma2_truth_dd[2, 2]), 
  mu_d[2] + 3 * sqrt(sigma2_truth_dd[2, 2]), 
  length.out = 500
)

# 作図用の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_{n,1}, x_{n,2})$がとり得る点(2値の組み合わせ)を作成します。
 $x_{n,1}$がとり得る値(x軸の値)をseq()で作成してx_1_vecとします。この例では、平均値を中心に標準偏差の3倍を範囲とします。length.out引数を使うと指定した要素数で等間隔に切り分けます。by引数を使うと切り分ける間隔を指定できます。処理が重い場合は、この値を調整してください。$i$軸の標準偏差$\sigma_{i,i}$は、分散の平方根$\sigma_{i,i} = \sqrt{\sigma_{i,i}^2}$です。平方根はsqrt()で計算できます。同様に、$x_{n,2}$がとり得る値(y軸の値)も作成してx_2_vecとします。
 x_1_vecx_2_vecの要素の全ての組み合わせを持つようにマトリクスを作成してx_point_matとします。これは、確率密度を等高線図にする際に格子状の点(2軸の全ての値が直交する点)を渡す必要があるためです。

 作成した$\mathbf{x}$の点は次のようになります。

# 確認
x_point_mat[1:5, ]
##           [,1] [,2]
## [1,] -48.48469  -10
## [2,] -48.19016  -10
## [3,] -47.89564  -10
## [4,] -47.60111  -10
## [5,] -47.30658  -10

 1列目がx軸の値、2列目がy軸の値に対応しています。

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

# 尤度を計算:式(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、分散共分散行列の引数sigmalambda_truth_ddを逆行列$\boldsymbol{\Sigma} = \boldsymbol{\Lambda}^{-1}$に変換して指定します。勿論sigma2_truth_ddを指定すればいいのですが、ここでは式に合わせてlambda_truth_ddを使っています。

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

# 確認
head(model_df)
## # A tibble: 6 x 3
##     x_1   x_2     density
##   <dbl> <dbl>       <dbl>
## 1 -48.5   -10 0.000000165
## 2 -48.2   -10 0.000000170
## 3 -47.9   -10 0.000000175
## 4 -47.6   -10 0.000000181
## 5 -47.3   -10 0.000000186
## 6 -47.0   -10 0.000000192

 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:20210517204658p: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.   :-42.78   Min.   : 9.102  
##  1st Qu.: 10.07   1st Qu.:35.695  
##  Median : 28.44   Median :50.304  
##  Mean   : 27.51   Mean   :48.154  
##  3rd Qu.: 44.29   3rd Qu.:63.380  
##  Max.   : 75.96   Max.   :85.500


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

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

f:id:anemptyarchive:20210517204715p:plain
観測データ:多次元ガウス分布

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

・事前分布の設定

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

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

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

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

 これまでのように、このパラメータを用いて事前分布(の確率密度)を計算してグラフで確認したいところですが、精度行列の分布はイメージしにくいため、精度パラメータの期待値$\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 -48.5   -10 0.00000177
## 2 -48.2   -10 0.00000181
## 3 -47.9   -10 0.00000185
## 4 -47.6   -10 0.00000189
## 5 -47.3   -10 0.00000193
## 6 -47.0   -10 0.00000197


 $\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:20210517204739p: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

 $\boldsymbol{\Lambda}$の事後分布のパラメータを

$$ \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,]  3.381533e-05 -5.842352e-06
## [2,] -5.842352e-06  5.069963e-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.0017583973 -0.0003038023
## [2,] -0.0003038023  0.0026363805


 $\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 -48.5   -10 0.0000000976
## 2 -48.2   -10 0.000000101 
## 3 -47.9   -10 0.000000104 
## 4 -47.6   -10 0.000000108 
## 5 -47.3   -10 0.000000111 
## 6 -47.0   -10 0.000000115


 $\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:20210517204917p: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}$の事後分布のパラメータを用いて、予測分布のパラメータを計算します。

# 次元数を取得
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}} \\ \hat{\nu}_s &= 1 - D + \hat{\nu} \end{aligned} $$

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

# 確認
mu_s_d
lambda_s_hat_dd
nu_s_hat
## [1] 25 50
##              [,1]         [,2]
## [1,]  0.001724582 -0.000297960
## [2,] -0.000297960  0.002585681
## [1] 51

 $\mathbf{X}$から$\hat{\boldsymbol{\nu}}_{*},\ \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 -48.5   -10 0.000000241
## 2 -48.2   -10 0.000000248
## 3 -47.9   -10 0.000000254
## 4 -47.6   -10 0.000000260
## 5 -47.3   -10 0.000000267
## 6 -47.0   -10 0.000000274


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

# 予測分布を作図
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, 
                         ", lambda_s_hat=(", paste(round(lambda_s_hat_dd, 5), collapse = ", "), ")", 
                         ", nu_s_hat=", nu_s_hat), 
       x = expression(x[1]), y = expression(x[2]), 
       color = "density")

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

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

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

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

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

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

・ライブラリの読み込み

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


・モデルの設定

# (既知の)平均パラメータを指定
mu_d <- c(25, 50)

# 真の分散共分散行列を指定
sigma2_truth_dd <- matrix(c(600, 90, 90, 400), nrow = 2, ncol = 2)

# 真の精度行列を計算
lambda_truth_dd <- solve(sigma2_truth_dd)


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

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

 事前分布のパラメータを用いて、予測分布のパラメータを計算しておきます。

・作図用のデータフレームの初期化

# 作図用のxのx軸の値を作成
x_1_vec <- seq(
  mu_d[1] - 3 * sqrt(sigma2_truth_dd[1, 1]), 
  mu_d[1] + 3 * sqrt(sigma2_truth_dd[1, 1]), 
  length.out = 250
)

# 作図用のxのy軸の値を作成
x_2_vec <- seq(
  mu_d[2] - 3 * sqrt(sigma2_truth_dd[2, 2]), 
  mu_d[2] + 3 * sqrt(sigma2_truth_dd[2, 2]), 
  length.out = 250
)

# 作図用のxの点を作成
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(nu * w_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, 
      ", lambda_s=(", paste(round(lambda_s_dd, 5), collapse = ", "), ")", 
      ", nu_s=", nu_s
    )
  )
)

 各試行の結果を同じデータフレームに格納していく必要があります。事後分布をposterior_df、予測分布をpredict_dfとして、初期値の結果を持つように作成しておきます。分布の推移と共にパラメータの値を表示するために、label列にパラメータの値を記録しておきます。ややこしければas.factor(paste0("iter=", n))として試行回数だけ表示するだけでもそれっぽくなります。

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

# 確認
head(posterior_df)
head(predict_df)
## # A tibble: 6 x 4
##     x_1   x_2    density label                            
##   <dbl> <dbl>      <dbl> <fct>                            
## 1 -48.5   -10 0.00000177 N=0, nu=2, W=(5e-04, 0, 0, 5e-04)
## 2 -47.9   -10 0.00000185 N=0, nu=2, W=(5e-04, 0, 0, 5e-04)
## 3 -47.3   -10 0.00000193 N=0, nu=2, W=(5e-04, 0, 0, 5e-04)
## 4 -46.7   -10 0.00000201 N=0, nu=2, W=(5e-04, 0, 0, 5e-04)
## 5 -46.1   -10 0.00000210 N=0, nu=2, W=(5e-04, 0, 0, 5e-04)
## 6 -45.5   -10 0.00000219 N=0, nu=2, W=(5e-04, 0, 0, 5e-04)
## # A tibble: 6 x 4
##     x_1   x_2    density label                                     
##   <dbl> <dbl>      <dbl> <fct>                                     
## 1 -48.5   -10 0.00000617 N=0, lambda_s=(5e-04, 0, 0, 5e-04), nu_s=1
## 2 -47.9   -10 0.00000624 N=0, lambda_s=(5e-04, 0, 0, 5e-04), nu_s=1
## 3 -47.3   -10 0.00000632 N=0, lambda_s=(5e-04, 0, 0, 5e-04), nu_s=1
## 4 -46.7   -10 0.00000639 N=0, lambda_s=(5e-04, 0, 0, 5e-04), nu_s=1
## 5 -46.1   -10 0.00000647 N=0, lambda_s=(5e-04, 0, 0, 5e-04), nu_s=1
## 6 -45.5   -10 0.00000655 N=0, lambda_s=(5e-04, 0, 0, 5e-04), nu_s=1


・推論処理

# データ数(試行回数)を指定
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
  
  # 事後分布の期待値を用いた分布を更新:式(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(nu * w_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, 
        ", lambda_s_hat=(", paste(round(lambda_s_dd, 5), collapse = ", "), ")", 
        ", nu_s_hat=", nu_s
      )
    )
  )
  
  # 推論結果を結合
  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()に指定します。

 初期値(事前分布)を含むため、フレーム数の引数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:20210517205031g:plain
精度パラメータの事後分布の期待値による分布の推移:多次元ガウス分布

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

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

参考文献

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

おわりに

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

【次節の内容】

www.anarchive-beta.com


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

 mvtnormパッケージからmvnfastパッケージに置き替えました。