はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。基本的な内容は「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」になります。「数式」と「プログラム」から理解するのが目標です。
この記事は、3.4.3項の内容です。「尤度関数を平均と精度が未知の多次元ガウス分布(多変量正規分布)」、「事前分布をガウス・ウィシャート分布」とした場合の「パラメータの事後分布」と「未観測値の予測分布」の計算をR言語で実装します。
省略してある内容等ありますので、本とあわせて読んでください。初学者な自分が理解できるレベルまで落として書き下していますので、分かる人にはかなりくどくなっています。同じような立場の人のお役に立てれば幸いです。
【数式読解編】
【他の節の内容】
【この節の内容】
・Rでやってみよう
人工的に生成したデータを用いて、ベイズ推論を行ってみましょう。
利用するパッケージを読み込みます。
# 3.4.3項で利用するパッケージ library(tidyverse) library(mvnfast)
mvnfastパッケージは、多次元ガウス分布に関するパッケージです。多次元ガウス分布に従う乱数生成関数rmvn()と確率密度関数dmvn()を使います。また、予測分布の計算において多次元スチューデントのt分布の確率密度関数dmvt()も使います。
・観測モデルの構築
まずは、観測モデルを設定します。この例では、尤度$p(\mathbf{X} | \boldsymbol{\mu}, \boldsymbol{\Lambda})$をデータごとに独立した多次元ガウス分布$\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}, \boldsymbol{\Lambda}^{-1})$とします。
尤度のパラメータを設定します。この実装例では2次元のグラフで表現するため、$D = 2$のときのみ動作します。
# 真の平均パラメータを指定 mu_truth_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_truth_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()で計算できます。
この例ではmu_truth_dとlambda_truth_ddの両方が未知の値であり、この値を求めるのが目的です。
グラフ用の点を作成します。
# 作図用のxのx軸の値を作成 x_1_vec <- seq( mu_truth_d[1] - 3 * sqrt(sigma2_truth_dd[1, 1]), mu_truth_d[1] + 3 * sqrt(sigma2_truth_dd[1, 1]), length.out = 500 ) # 作図用のxのy軸の値を作成 x_2_vec <- seq( mu_truth_d[2] - 3 * sqrt(sigma2_truth_dd[2, 2]), mu_truth_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_vecとx_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_truth_d, sigma = solve(lambda_truth_dd) ) )
x_point_matの値(の組み合わせ)ごとに確率密度を計算します。多次元ガウス分布の確率密度は、mvnfast::dmvn()で計算できます。データの引数Xにx_point_mat、平均の引数muにmu_truth_d、分散共分散行列の引数sigmaにlambda_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_truth_d, 1), collapse = ", "), ")", ", lambda=(", paste(round(lambda_truth_dd, 5), collapse = ", "), ")"), x = expression(x[1]), y = expression(x[2]), color = "density")
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_truth_dd))
生成するデータ数$N$をNとして値を指定します。
多次元ガウス分布に従う乱数は、mvnfast::rmvn()で生成できます。試行回数の引数nにNを指定します。他の引数についてはmvnfast::dmvn()と同じです。生成したN個のデータをx_ndとします。
観測したデータを確認しましょう。
# 確認 summary(x_nd)
## V1 V2 ## Min. :-27.140 Min. : 9.207 ## 1st Qu.: 8.703 1st Qu.:40.316 ## Median : 28.062 Median :55.253 ## Mean : 28.088 Mean :54.263 ## 3rd Qu.: 41.454 3rd Qu.:70.668 ## Max. : 91.821 Max. :89.581
観測データの散布図を尤度と重ねて確認します。
# 観測データのデータフレームを作成 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 = ", "), ")", ", lambda=(", paste(round(lambda_truth_dd, 5), collapse = ", "), ")"), x = expression(x[1]), y = expression(x[2]), color = "density")
作図用に観測データをデータフレームに格納しておき、確率密度のグラフと重ねて表示します。
・事前分布の設定
次に、尤度に対する共役事前分布を設定します。多次元ガウス分布の平均パラメータ$\boldsymbol{\mu}$と精度パラメータ$\boldsymbol{\Lambda}$に対する事前分布$p(\boldsymbol{\mu}, \boldsymbol{\Lambda})$として、ガウス・ウィシャート分布$\mathrm{NW}(\boldsymbol{\mu}, \boldsymbol{\Lambda} | \mathbf{m}, \beta, \nu, \mathbf{W})$を設定します。
$\boldsymbol{\mu}$の事前分布のパラメータ(超パラメータ)を設定します。
# muの事前分布のパラメータを指定 m_d <- c(0, 0) beta <- 1 # lambdaの事前分布のパラメータを指定 nu <- 2 w_dd <- matrix(c(0.0005, 0, 0, 0.0005), nrow = 2, ncol = 2)
多次元ガウス分布の平均パラメータ$\mathbf{m} = (m_1, \cdots, m_D)$をm_d、精度パラメータの係数$\beta$をbetaとして値を指定します。
ウィシャート分布の自由度$\nu$をnu、パラメータ$\mathbf{W} = (w_{1,1}, \cdots, w_{D,D})$をw_ddとして値を指定します。$\nu$は$D - 1$より大きい値、$\mathbf{W}$は正定値行列である必要があります。
事前分布のパラメータを設定できたので、これまでのように事前分布をグラフで確認しましょう。しかし、この例では精度パラメータ$\boldsymbol{\Lambda}$が未知の値なので、$\boldsymbol{\mu}$の事前分布$\mathcal{N}(\boldsymbol{\mu} | \mathbf{m}, (\beta \boldsymbol{\Lambda})^{-1})$の精度パラメータ$\beta \boldsymbol{\Lambda}$が分かりません。そこで、$\boldsymbol{\Lambda}$の事前分布$\mathcal{W}(\boldsymbol{\Lambda} | \nu, \mathbf{W})$から精度パラメータの期待値$\mathbb{E}[\boldsymbol{\Lambda}]$を求めて、$\boldsymbol{\Lambda}$の代わりに利用することにします。
# lambdaの期待値を計算:式(2.89) E_lambda_dd <- nu * w_dd
精度パラメータの期待値は、ウィシャート分布の期待値(2.89)より
で計算できます。
作図用の$\boldsymbol{\mu}$の点を作成します。
# muの事前分布の分散共分散行列を計算 E_sigma2_mu_dd <- solve(beta * E_lambda_dd) # 作図用のmuのx軸の値を作成 mu_1_vec <- seq( mu_truth_d[1] - 2 * sqrt(E_sigma2_mu_dd[1, 1]), mu_truth_d[1] + 2 * sqrt(E_sigma2_mu_dd[1, 1]), length.out = 500 ) # 作図用のmuのy軸の値を作成 mu_2_vec <- seq( mu_truth_d[2] - 2 * sqrt(E_sigma2_mu_dd[2, 2]), mu_truth_d[2] + 2 * sqrt(E_sigma2_mu_dd[2, 2]), length.out = 500 ) # 作図用のmuの点を作成 mu_point_mat <- cbind( rep(mu_1_vec, times = length(mu_2_vec)), rep(mu_2_vec, each = length(mu_1_vec)) )
$\mathbf{x}$の点のときと同様に、2次元ガウス分布に従う変数$\boldsymbol{\mu} = (\mu_1, \mu_2)$がとり得る値を作成してmu_point_matとします。この例では、真の値を中心に$\boldsymbol{\mu}$の事前分布の標準偏差の2倍を範囲とします。$(\beta \mathbb{E}[\boldsymbol{\Lambda}])^{-1}$を$\boldsymbol{\mu}$の事前分布の分散共分散行列E_sigma2_mu_ddとして標準偏差を計算しています。
$\boldsymbol{\mu}$の事前分布の確率密度を計算します。
# muの事前分布を計算:式(2.72) prior_mu_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(beta * E_lambda_dd) ) )
m_dとbeta, E_lambda_ddを用いて、尤度のときと同様にしてmu_point_matの値ごとに確率密度を計算します。
計算結果は次のようになります。
# 確認 head(prior_mu_df)
## # A tibble: 6 x 3 ## mu_1 mu_2 density ## <dbl> <dbl> <dbl> ## 1 -38.2 -13.2 0.0000702 ## 2 -38.0 -13.2 0.0000708 ## 3 -37.7 -13.2 0.0000715 ## 4 -37.5 -13.2 0.0000722 ## 5 -37.2 -13.2 0.0000729 ## 6 -37.0 -13.2 0.0000736
$\boldsymbol{\mu}$の事前分布を作図します。
# 作図用に真のmuのデータフレームを作成 mu_df <- tibble( mu_1 = mu_truth_d[1], mu_2 = mu_truth_d[2] ) # muの事前分布を作図 ggplot() + geom_contour(data = prior_mu_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 = ", "), ")", ", beta=", beta, ", E[lambda]=(", paste(round(E_lambda_dd, 5), collapse = ", "), ")"), x = expression(mu[1]), y = expression(mu[2]), color = "density")
$\boldsymbol{\mu}$の真の値をデータフレームに格納しておき、確率密度のグラフと重ねて表示します。
$\boldsymbol{\Lambda}$についても事前分布(の確率密度)を計算してグラフで確認したいところですが、精度行列の分布はイメージしにくいため、平均パラメータの期待値$\mathbb{E}[\boldsymbol{\mu}]$と精度パラメータの期待値$\mathbb{E}[\boldsymbol{\Lambda}]$を用いた多次元ガウス分布を尤度と比較することにしましょう。ちなみに、ウィシャート分布の確率密度はMCMCpack::dwish()で計算できます(たぶん)。
# muの期待値を計算:式(2.76) E_mu_d <- m_d
平均パラメータの期待値は、ガウス分布の期待値(3.76)より
です。
$\boldsymbol{\mu}$の事前分布の期待値(平均パラメータの期待値)と$\boldsymbol{\Lambda}$の事前分布の期待値(精度パラメータの期待値)を用いて、多次元ガウス分布の確率密度を計算します。
# 事前分布の期待値を用いた分布を計算:式(2.72) prior_E_df <- tibble( x_1 = x_point_mat[, 1], x_2 = x_point_mat[, 2], density = mvnfast::dmvn( X = x_point_mat, mu = E_mu_d, sigma = solve(E_lambda_dd) ) )
E_mu_d, E_lambda_ddを用いて、尤度のときと同様にして計算します。
計算結果は次のようになります。
# 確認 head(prior_E_df)
## # A tibble: 6 x 3 ## x_1 x_2 density ## <dbl> <dbl> <dbl> ## 1 -48.5 -10 0.0000467 ## 2 -48.2 -10 0.0000474 ## 3 -47.9 -10 0.0000481 ## 4 -47.6 -10 0.0000488 ## 5 -47.3 -10 0.0000494 ## 6 -47.0 -10 0.0000501
$\boldsymbol{\mu}$と$\boldsymbol{\Lambda}$の事前分布の期待値(平均と精度パラメータの期待値)を用いた分布を作図します。
# 事前分布の期待値を用いた分布を作図 ggplot() + geom_contour(data = prior_E_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 labs(title = "Multivariate Gaussian Distribution", subtitle = paste0("m=(", paste(m_d, collapse = ", "), ")", ", nu=", nu, ", W=(", paste(round(w_dd, 5), collapse = ", "), ")"), x = expression(x[1]), y = expression(x[2]), color = "density")
尤度を破線で重ねて描画します。
・事後分布の計算
観測データ$\mathbf{X}$からパラメータ$\boldsymbol{\mu},\ \boldsymbol{\Lambda}$の事後分布$p(\boldsymbol{\mu}, \boldsymbol{\Lambda} | \mathbf{X})$を求めます(パラメータ$\boldsymbol{\mu},\ \boldsymbol{\Lambda}$を分布推定します)。事後分布はガウス・ウィシャート分布$\mathrm{NW}(\boldsymbol{\mu} | \hat{\mathbf{m}}, \hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1}, \hat{\nu}, \hat{\mathbf{W}})$になります。
観測データx_ndを用いて、$\boldsymbol{\mu}$の事後分布のパラメータを計算します。
# muの事後分布のパラメータを計算:式(3.129) beta_hat <- N + beta m_hat_d <- (colSums(x_nd) + beta * m_d) / beta_hat
$\boldsymbol{\mu}$の事後分布のパラメータを
で計算して、結果をbeta_hat、m_hat_dとします。
# 確認
beta_hat
m_hat_d
## [1] 51 ## [1] 27.53769 53.19898
$\boldsymbol{\Lambda}$の事後分布のパラメータも計算します。
# lambdaの事後分布のパラメータを計算:式(3.133) term_x <- t(x_nd) %*% as.matrix(x_nd) term_m <- beta * as.matrix(m_d) %*% t(m_d) term_m_hat <- beta_hat * as.matrix(m_hat_d) %*% t(m_hat_d) w_hat_dd <- solve( term_x + term_m - term_m_hat + solve(w_dd) ) nu_hat <- N + nu
$\boldsymbol{\Lambda}$の事後分布のパラメータを
で計算して、結果をw_hat_dd、nu_hatとします。ただし、$\sum_{n=1}^N \mathbf{x}_n \mathbf{x}_n^{\top}$の計算を効率よく処理するために、$\mathbf{X}^{\top} \mathbf{X}$で計算しています。
# 確認
nu_hat
w_hat_dd
## [1] 52 ## [,1] [,2] ## [1,] 2.861227e-05 -9.771505e-06 ## [2,] -9.771505e-06 4.348895e-05
事前分布のときと同様に、求めたパラメータを用いて、学習後の精度パラメータの期待値$\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.0014878378 -0.0005081183 ## [2,] -0.0005081183 0.0022614254
$\boldsymbol{\mu}$の事後分布の確率密度を計算します。
# muの事後分布を計算:式(2.72) posterior_mu_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(beta_hat * E_lambda_hat_dd) ) )
更新した超パラメータm_hat_d, lambda_lambda_hat_ddを用いて、事前分布のときと同様にして計算します。
計算結果は次のようになります。
# 確認 head(posterior_mu_df)
## # A tibble: 6 x 3 ## mu_1 mu_2 density ## <dbl> <dbl> <dbl> ## 1 -38.2 -13.2 3.00e-135 ## 2 -38.0 -13.2 6.85e-135 ## 3 -37.7 -13.2 1.56e-134 ## 4 -37.5 -13.2 3.53e-134 ## 5 -37.2 -13.2 7.94e-134 ## 6 -37.0 -13.2 1.78e-133
$\boldsymbol{\mu}$の事後分布を作図します。
# muの事後分布を作図 ggplot() + geom_contour(data = posterior_mu_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 = ", "), ")", ", beta_hat=", beta_hat, ", E[lambda]=(", paste(round(E_lambda_hat_dd, 5), collapse = ", "), ")"), x = expression(mu[1]), y = expression(mu[2]), color = "density")
$\boldsymbol{\mu}$の真の値付近をピークとする分布を推定できています。
x軸とy軸の値に注目すると、事前分布と比べて範囲が狭くなっています。グラフの表示範囲を広げてみましょう。
# muの事後分布を作図 ggplot() + geom_contour(data = posterior_mu_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 = ", "), ")", ", beta_hat=", beta_hat, ", E[lambda]=(", paste(round(E_lambda_hat_dd, 5), collapse = ", "), ")"), x = expression(mu[1]), y = expression(mu[2]), color = "density")
先のとがった分布になっているのが分かります。
$\boldsymbol{\Lambda}$の事後分布についても事前分布のときと同様に、パラメータの期待値を用いた分布を確認しましょう。学習後の平均パラメータの期待値$\mathbb{E}[\hat{\boldsymbol{\mu}}]$を計算します。
# muの期待値を計算:式(2.76) E_mu_hat_d <- m_hat_d
$\boldsymbol{\mu}$と$\boldsymbol{\Lambda}$の事後分布の期待値(平均と精度パラメータの期待値)を用いて、多次元ガウス分布の確率密度を計算します。
# 事後分布の期待値を用いた分布を計算:式(2.72) posterior_E_df <- tibble( x_1 = x_point_mat[, 1], x_2 = x_point_mat[, 2], density = mvnfast::dmvn( X = x_point_mat, mu = E_mu_hat_d, sigma = solve(E_lambda_hat_dd) ) )
E_mu_hat_d, E_lambda_hat_ddを用いて、尤度のときと同様にして計算します。
計算結果は次のようになります。
# 確認 head(posterior_E_df)
## # A tibble: 6 x 3 ## x_1 x_2 density ## <dbl> <dbl> <dbl> ## 1 -48.5 -10 0.000000478 ## 2 -48.2 -10 0.000000490 ## 3 -47.9 -10 0.000000501 ## 4 -47.6 -10 0.000000513 ## 5 -47.3 -10 0.000000526 ## 6 -47.0 -10 0.000000538
$\boldsymbol{\mu}$と$\boldsymbol{\Lambda}$の事後分布の期待値(平均と精度パラメータの期待値)を用いた分布を作図します。
# 事後分布の期待値を用いた分布を作図 ggplot() + geom_contour(data = posterior_E_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 labs(title = "Multivariate Gaussian Distribution", subtitle = paste0("N=", N, ", m_hat=(", paste(round(m_hat_d, 1), collapse = ", "), ")", ", nu_hat=", nu_hat, ", W_hat=(", paste(round(w_hat_dd, 5), collapse = ", "), ")"), x = expression(x[1]), y = expression(x[2]), color = "density")
分布の形状が近づいているのを確認できます。
・予測分布の計算
最後に、観測データ$\mathbf{X}$から未観測のデータ$\mathbf{x}_{*}$の予測分布$p(\mathbf{x}_{*} | \mathbf{X})$を求めます。予測分布はスチューデントのt分布$p(\mathbf{x}_{*} | \hat{\boldsymbol{\mu}}_s, \hat{\boldsymbol{\Lambda}}_s, \hat{\nu}_s)$になります。
事後分布のパラメータを用いて、予測分布のパラメータを計算します。
# 次元数を取得 D <- length(m_d) # 予測分布のパラメータを計算:式(3.140') mu_s_hat_d <- m_hat_d lambda_s_hat_dd <- (1 - D + nu_hat) * beta_hat / (1 + beta_hat) * w_hat_dd nu_s_hat <- 1 - D + nu_hat
予測分布のパラメータを
で計算して、結果をmu_s_hat_d、lambda_s_hat_dd、nu_s_hatとします。
# 確認
mu_s_hat_d
lambda_s_hat_dd
nu_s_hat
## [1] 27.53769 53.19898 ## [,1] [,2] ## [1,] 0.0014311635 -0.0004887632 ## [2,] -0.0004887632 0.0021752838 ## [1] 51
$\mathbf{X}$から$\hat{\boldsymbol{\mu}}_s,\ \hat{\boldsymbol{\Lambda}}_s,\ \hat{\nu}_s$を学習しているのが式からも分かります。
求めたパラメータを用いて、予測分布の確率密度を計算します。
# 予測分布を計算:式(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_hat_d, sigma = solve(lambda_s_hat_dd), df = nu_s_hat ) )
x_point_matの値ごとに確率密度を計算します。多次元のスチューデントのt分布の確率密度は、mvnfast::dmvt()で計算できます。データの引数Xにはx_point_mat、平均の引数muにmu_s_hat_d、スケール行列(?)の引数sigmaにlambda_s_hat_ddの逆行列、自由度の引数dfにnu_s_hatを指定します。
計算結果は次のようになります。
# 確認 head(predict_df)
## # A tibble: 6 x 3 ## x_1 x_2 density ## <dbl> <dbl> <dbl> ## 1 -48.5 -10 0.000000894 ## 2 -48.2 -10 0.000000911 ## 3 -47.9 -10 0.000000929 ## 4 -47.6 -10 0.000000946 ## 5 -47.3 -10 0.000000965 ## 6 -47.0 -10 0.000000983
予測分布を尤度と重ねて作図します。
# 予測分布を作図 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 Student's t Distribution", subtitle = paste0("N=", N, ", mu_s_hat=(", paste(round(mu_s_hat_d, 1), collapse = ", "), ")", ", 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")
観測データが増えると、予測分布が真の分布に近づきます。
・おまけ：アニメーションで推移の確認
gganimateパッケージを利用して、事後分布と予測分布の推移をアニメーション(gif画像)で確認するためのコードです。
・コード(クリックで展開)
異なる点のみを簡単に解説します。
・ライブラリの読み込み
# 3.4.3項で利用するパッケージ library(tidyverse) library(mvnfast) library(gganimate)
・モデルの設定
# 真の平均パラメータを指定 mu_truth_d <- c(25, 50) # 真の分散共分散行列を指定 sigma2_truth_dd <- matrix(c(600, 90, 90, 400), nrow = 2, ncol = 2) # 真の精度行列を計算 lambda_truth_dd <- solve(sigma2_truth_dd) # muの事前分布のパラメータを指定 m_d <- c(0, 0) beta <- 1 # lambdaの事前分布のパラメータを指定 nu <- 2 w_dd <- matrix(c(0.0005, 0, 0, 0.0005), nrow = 2, ncol = 2) # 初期値による予測分布のパラメータを計算:式(3.140) mu_s_d <- m_d lambda_s_dd <- (nu - 1) * beta / (1 + beta) * w_dd nu_s <- nu - 1
事前分布のパラメータを用いて、予測分布のパラメータを計算しておきます。
・作図用のデータフレームの初期化
# 作図用のmuのx軸の値を作成 mu_1_vec <- seq( mu_truth_d[1] - 1 * sqrt(solve(beta * nu * w_dd)[1, 1]), mu_truth_d[1] + 1 * sqrt(solve(beta * nu * w_dd)[1, 1]), length.out = 250 ) # 作図用のmuのy軸の値を作成 mu_2_vec <- seq( mu_truth_d[2] - 1 * sqrt(solve(beta * nu * w_dd)[2, 2]), mu_truth_d[2] + 1 * sqrt(solve(beta * nu * w_dd)[2, 2]), length.out = 250 ) # 作図用のmuの点を作成 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_mu_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(beta * nu * w_dd) ), label = as.factor( paste0( "N=", 0, ", m=(", paste(round(m_d, 1), collapse = ", "), ")", ", beta=", beta, ", E[lambda]=(", paste0(round(nu * w_dd, 5), collapse = ", "), ")" ) ) ) # 作図用のxのx軸の値を作成 x_1_vec <- seq( mu_truth_d[1] - 3 * sqrt(sigma2_truth_dd[1, 1]), mu_truth_d[1] + 3 * sqrt(sigma2_truth_dd[1, 1]), length.out = 250 ) # 作図用のxのy軸の値を作成 x_2_vec <- seq( mu_truth_d[2] - 3 * sqrt(sigma2_truth_dd[2, 2]), mu_truth_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)) ) # 初期値による予測分布を計算:式(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, ", mu_s=(", paste(round(mu_s_d, 1), collapse = ", "), ")", ", lambda_s=(", paste(round(lambda_s_dd, 5), collapse = ", "), ")", ", nu_s=", nu_s ) ) )
各試行の結果を同じデータフレームに格納していく必要があります。事後分布をposterior_mu_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.0000194 N=0, mu_s=(0, 0), lambda_s=(0.00025, 0, 0, 0.00025), nu~ ## 2 -47.9 -10 0.0000197 N=0, mu_s=(0, 0), lambda_s=(0.00025, 0, 0, 0.00025), nu~ ## 3 -47.3 -10 0.0000200 N=0, mu_s=(0, 0), lambda_s=(0.00025, 0, 0, 0.00025), nu~ ## 4 -46.7 -10 0.0000202 N=0, mu_s=(0, 0), lambda_s=(0.00025, 0, 0, 0.00025), nu~ ## 5 -46.1 -10 0.0000205 N=0, mu_s=(0, 0), lambda_s=(0.00025, 0, 0, 0.00025), nu~ ## 6 -45.5 -10 0.0000208 N=0, mu_s=(0, 0), lambda_s=(0.00025, 0, 0, 0.00025), nu~
・推論処理
# データ数(試行回数)を指定 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_truth_dd)) %>% as.vector() # muの事後分布のパラメータを更新:式(3.129) old_beta <- beta old_m_d <- m_d beta <- 1 + beta m_d <- (x_nd[n, ] + old_beta * m_d) / beta # lambdaの事後分布のパラメータを更新:式(3.133) term_x <- x_nd[n, ] %*% t(x_nd[n, ]) term_m <- old_beta * old_m_d %*% t(old_m_d) term_m_hat <- beta * m_d %*% t(m_d) w_dd <- solve( term_x + term_m - term_m_hat + solve(w_dd) ) nu <- 1 + nu # muの事前分布を計算:式(2.72) tmp_posterior_mu_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(beta * nu * w_dd) ), label = as.factor( paste0( "N=", n, ", m_hat=(", paste(round(m_d, 1), collapse = ", "), ")", ", beta_hat=", beta, ", E[lambda]=(", paste0(round(nu * w_dd, 5), collapse = ", "), ")" ) ) ) # 予測分布のパラメータを更新:式(1.40) mu_s_d <- m_d lambda_s_dd <- (nu - 1) * beta / (1 + beta) * 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, ", mu_s_hat=(", paste(round(mu_s_d, 1), collapse = ", "), ")", ", lambda_s_hat=(", paste(round(lambda_s_dd, 5), collapse = ", "), ")", ", nu_s=", nu_s ) ) ) # 推論結果を結合 posterior_mu_df <- rbind(posterior_mu_df, tmp_posterior_mu_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回繰り返し1とx_nd[n, ]を加えることで行います。n回目のループ処理のときには、n-1回分の1とx_nd[n, ]が既にm_dとlambda_mu_ddに加えられているわけです。
ただし、事後分布のパラメータの計算において、更新前(1ステップ前)のパラメータを使うため、それぞれold_***として値を一時的に保存しておきます。
・平均パラメータの事後分布の推移
# 作図用に真のmuのデータフレームを作成 mu_df <- tibble( mu_1 = mu_truth_d[1], mu_2 = mu_truth_d[2] ) # muの事後分布を作図 posterior_graph <- ggplot() + geom_contour(data = posterior_mu_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 transition_manual(label) + # フレーム labs(title = "Multivariate Gaussian Distribution", subtitle = "{current_frame}", x = expression(mu[1]), y = expression(mu[2]), color = "density") # gif画像を作成 animate(posterior_graph, nframes = N + 1, fps = 10)
各フレームの順番を示す列(label)をgganimate::transition_manual()に指定します。
初期値(事前分布)を含むため、フレーム数の引数nframesはN + 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_truth_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_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..), alpha = 0.5, linetype = "dashed") + # 尤度 geom_point(data = mu_df, aes(x = mu_1, y = mu_2), color = "red", shape = 4, size = 5) + # 真のmu 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と一致している必要があります。
(始めの方は精度が低い(分散が大きい)ため確率密度が低く等高線が表示されていません。)
新たなデータによって平均(分布の中心)と精度(分布の広がり)が推移しているのを確認できます。
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
データに従って変形する予測分布を見てると落ち着く。ところで、事後分布のアニメーションが何か少し変な気が、予測分布はアニメーションの方だとうまくいってるけどメインの方が何か変な気がする(データ数が前2つよりも多めに必要なんだと思う)。
【次節の内容】
- 2021/04/12：加筆修正しました。その際に数式読解編と記事を分割しました。