はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。基本的な内容は「数式の行間を読んでみた」とそれを「Rで組んでみた」になります。「数式」と「プログラム」から理解するのが目標です。
この記事は3.4.2項の内容です。尤度関数を多次元ガウス分布、事前分布をウィシャート分布とした場合の精度パラメータの事後分布を導出し、また学習した事後分布を用いた予測分布を導出します。またその推論過程をR言語で実装します。
省略してある内容等ありますので、本と併せて読んでください。初学者な自分が理解できるレベルまで落として書き下していますので、分かる人にはかなりくどくなっています。同じような立場の人のお役に立てれば幸いです。
【前節の内容】
【他の節一覧】
【この節の内容】
3.4.2 精度が未知の場合
・観測モデルの設定
多次元ガウス分布に従うと仮定する$N$個の観測データ$\mathbf{X} = \{\mathbf{x}_1, \mathbf{x}_2, \cdots, \mathbf{x}_N\}$を用いて精度パラメータ$\boldsymbol{\Lambda}$の事後分布を求めていく。
観測データ$\mathbf{X}$、平均パラメータ$\boldsymbol{\mu}$、未知の共分散行列$\boldsymbol{\Sigma}$、精度行列(パラメータ)$\boldsymbol{\Lambda}$をそれぞれ
としたとき、観測モデルは次のようになる。
また$\boldsymbol{\Lambda}$の事前分布$p(\boldsymbol{\Lambda})$を、ウィシャート分布
とする。$\nu$は自由度であり、$\nu > D - 1$の値をとる($\mathbf{W}$は何？)。
・事後分布の導出
観測データ$\mathbf{X}$によって学習した$\boldsymbol{\Lambda}$の事後分布$p(\boldsymbol{\Lambda} | \mathbf{X})$は、観測モデルに対してベイズの定理を用いて
となる。分母の$p(\mathbf{X})$は$\boldsymbol{\Lambda}$に影響しないため省略して、比例関係にのみ注目する。省略した部分については、最後に正規化することで対応できる。
次にこの分布の具体的な形状を明らかにしていく。対数をとって指数部分の計算を分かりやすくして、$\boldsymbol{\Lambda}$に関して整理すると
【途中式の途中式】
- それぞれ具体的な式に置き換える。
- 1つ目の因子をトレース$\mathrm{Tr}(\cdot)$を使って置き換える。
精度パラメータを
として、また計算を分かりやすくするため$(\mathbf{x}_n - \mu)^{\top} = (x_{n1} - \mu_1, x_{n2} - \mu_2, \cdots, x_{nD} - \mu_D) = (\tilde{x}_{n1}, \tilde{x}_{n2}, \cdots, \tilde{x}_{nD})$とおくと(他に適切な記号があれば教えて)、$(\mathbf{x}_n - \boldsymbol{\mu})^{\top} \boldsymbol{\Lambda} (\mathbf{x}_n - \boldsymbol{\mu})$は
となる。
この式を整理するために、$(\mathbf{x}_n - \boldsymbol{\mu}) (\mathbf{x}_n - \boldsymbol{\mu})^{\top} \boldsymbol{\Lambda}$を考える。
$(\mathbf{x}_n - \boldsymbol{\mu}) (\mathbf{x}_n - \boldsymbol{\mu})^{\top} \boldsymbol{\Lambda}$は、$D \times D$の正方行列になる。この行列の対角成分の和$\mathrm{Tr} \Bigl((\mathbf{x}_n - \boldsymbol{\mu}) (\mathbf{x}_n - \boldsymbol{\mu})^{\top} \boldsymbol{\Lambda}\Bigr)$は
となる。
よって
であることが分かる。
- 式(A.12)より、$\mathrm{Tr} \Bigl(\sum_{n=1}^N (\mathbf{x}_n - \mu) (\mathbf{x}_n - \mu)^{\top} \boldsymbol{\Lambda}\Bigr) + \mathrm{Tr}(\mathbf{W}^{-1} \boldsymbol{\Lambda})= \mathrm{Tr} \Bigl(\sum_{n=1}^N (\mathbf{x}_n - \mu) (\mathbf{x}_n - \mu)^{\top} \boldsymbol{\Lambda} + \mathbf{W}^{-1} \boldsymbol{\Lambda}\Bigr)$であることから、式を整理する。
となる。適宜$\boldsymbol{\Lambda}$に影響しない項を$\mathrm{const.}$にまとめている。
式の形から事後分布もウィシャート分布となることが分かる。そこで事後分布を次のようにおき
この式の対数をとり、$\boldsymbol{\Lambda}$に関して整理すると
となる。
従って式(3.114)との対応関係から、事後分布のパラメータは
と求められる。
ちなみに分散共分散行列$\boldsymbol{\Sigma}$の事後分布の場合、逆ウィシャート分布になる。
・予測分布の導出
続いて未観測のデータ$\mathbf{x}_{*}$に対する予測分布$p(\mathbf{x}_{*})$を求めていく。
1次元ガウス分布のときと同様に、パラメータ$\boldsymbol{\mu}$を周辺化(積分計算)して予測分布
を求めることは避け、ベイズの定理を用いて
予測分布を求める。
この式の両辺の対数をとり、$p(\mathbf{x}_{*})$に関して式を整理すると
で計算できることが分かる。
2つ目の項は、データが1つ($N = 1$)の事後分布(3.115)と捉えられることから、式(3.116)よりパラメータを
とおくと
となる。
この式を式(3.117)に代入して、$\mathbf{x}_{*}$に関して式を整理すると
【途中式の途中式】
- 具体的な式に置き換える。
- $\nu_{\mathbf{x}_{*}},\ \mathbf{W_{\mathbf{x}_{*}}}$にそれぞれ式(3.117)を代入する。このとき式(A.18)より、$|\mathbf{W}_{\mathbf{x}_{*}}| = |\mathbf{W}_{\mathbf{x}_{*}}^{-1}|^{-1}$である。
- 2行目の項は式(A.12)より、3行目の項は$\mathbf{W}^{-1} + (\mathbf{x}_{*} - \mu) (\mathbf{x}_{*} - \mu)^{\top} = \mathbf{W}^{-1} \{\mathbf{I}_D + \mathbf{W} (\mathbf{x}_{*} - \mu) (\mathbf{x}_{*} - \mu)^{\top}\}$とできるので式(A.17)より、それぞれ分解する。
- 式を整理する。
- 事後分布の導出時より、$(\mathbf{x}_n - \boldsymbol{\mu})^{\top} \boldsymbol{\Lambda} (\mathbf{x}_n - \boldsymbol{\mu}) = \mathrm{Tr} ((\mathbf{x}_n - \boldsymbol{\mu}) (\mathbf{x}_n - \boldsymbol{\mu})^{\top}\boldsymbol{\Lambda})$である。
- $\mathbf{W} (\mathbf{x}_{*} - \boldsymbol{\mu})$を1つの行列とみると、$(D \times 1)$、$(1 \times D)$の行列の積なので、式(A.19)の変形を行う。
- 行列式の定義より、スカラの行列式はスカラ$|a| = a$である。また式(A.2)より、転置する。
- $(\mathbf{x}_{*} - \boldsymbol{\mu})^{\top} \mathbf{W}^{\top} (\mathbf{x}_{*} - \boldsymbol{\mu})$はスカラなので転置しても影響を受けないため、$(\mathbf{x}_{*} - \boldsymbol{\mu})^{\top} \mathbf{W}^{\top} (\mathbf{x}_{*} - \boldsymbol{\mu}) = \{(\mathbf{x}_{*} - \boldsymbol{\mu})^{\top} \mathbf{W}^{\top} (\mathbf{x}_{*} - \boldsymbol{\mu})\}^{\top}$である。また式(A.2)と同様に、$(A B C)^{\top} = C^{\top} B^{\top} A^{\top}$である。
となる。適宜$\mathbf{x}_{*}$に影響しない項を$\mathrm{const.}$にまとめている。
式の形から予測分布は多次元のスチューデントのt分布となることが分かる。そこで次のようにおき
この式の対数をとり、$\mathbf{x}_{*}$に関して整理すると
となる。
従って式(3.120)との対応関係から、予測分布のパラメータは
と求まる。
事前分布のパラメータ$\nu,\ \mathbf{W}$を事後分布のパラメータ$\hat{\nu},\ \hat{\mathbf{W}}$に置き換えると、観測データによって学習を行った予測分布$p(\mathbf{x}_{*} | \mathbf{X})$のパラメータ
が得られる。
・Rでやってみよう
多次元ガウス分布に従いランダムに生成したデータを用いて、パラメータを推定してみましょう。
利用するパッケージを読み込みます。
# 利用パッケージ library(tidyverse) library(mvtnorm) library(mvnfast)
mvtnormパッケージは、多次元ガウス分布に関するパッケージです。多次元ガウス分布に従う乱数生成関数rmvnorm()と確率密度計算関数dmvnorm()を使います。
`mvnfastパッケージは、多次元スチューデントのt分布の確率密度関数dmvt()を使います。
・観測モデルの設定
観測モデルのパラメータを指定します。この例では2次元のグラフで表現するため、$D = 2$のときのみ動作します。
# データ数を指定 N <- 100 D <- 2 # (固定) # 観測モデルのパラメータを指定 mu_d <- c(25, 50) sigma_truth_dd <- matrix(c(50, 30, 30, 50), nrow = 2, ncol = 2) lambda_truth_dd <- solve(sigma_truth_dd^2)
観測モデルの平均パラメータ$\boldsymbol{\mu}$をmu_dとします。
分散共分散行列$\boldsymbol{\Sigma}$ではなく、分散共分散の平方根$(\sigma_{11}, \cdots, \sigma_{DD})$をsigma_truth_ddとします。これは作図時に標準偏差を利用するためです。matrix()のデフォルトの仕様上、$(\sigma_{11}, \sigma_{12}, \sigma_{21}, \sigma_{22})$の順番に値を指定します。
sigma_truth_ddの各要素を2乗したものが$\boldsymbol{\Sigma}$です。更にその逆行列が精度パラメータ$\boldsymbol{\Lambda}$でした。逆行列の計算はsolve()で行い、その計算結果をlambda_truth_ddとします。この値を推論するのがこの項の目的になります。
次に事前分布のパラメータを指定します。
# 事前分布のパラメータを指定 W_dd <- matrix(c(0.00005, 0, 0, 0.00005), nrow = 2, ncol = 2) nu <- 2
ウィシャート事前分布のパラメータ$\mathbf{W},\ \nu$をそれぞれW_dd、nuとします。
作図時に利用する格子状の点(データ)と平均の点のデータフレームを作成しておきます。
# 作図用の点を生成 x_vec <- seq(mu_d[1] - 2 * sigma_truth_dd[1, 1], mu_d[1] + 2 * sigma_truth_dd[1, 1], by = 0.5) y_vec <- seq(mu_d[2] - 2 * sigma_truth_dd[2, 2], mu_d[2] + 2 * sigma_truth_dd[2, 2], by = 0.5) point_df <- tibble( x = rep(x_vec, times = length(y_vec)), y = rep(y_vec, each = length(x_vec)) ) mu_df <- tibble( x = mu_d[1], y = mu_d[2] )
この例では、平均から標準偏差の2倍の範囲をグラフ化することにします。
設定した観測モデルのパラメータに従って乱数を生成します。
# 2次元ガウス分布に従うデータを生成 x_nd <- mvtnorm::rmvnorm(n = N, mean = mu_d, sigma = sigma_truth_dd^2) summary(x_nd)
## V1 V2 ## Min. :-95.30 Min. :-49.47 ## 1st Qu.: -9.20 1st Qu.: 22.05 ## Median : 29.89 Median : 62.35 ## Mean : 33.45 Mean : 61.04 ## 3rd Qu.: 75.12 3rd Qu.: 99.86 ## Max. :158.02 Max. :166.78
mvtnorm::rmvnorm()で多次元ガウス分布に従う乱数を生成し、観測データ$\mathbf{X}$(x_nd)とします。
引数nにはデータ数N、引数meanには平均パラメータmu_d、引数sigmaにはsimga_truth_ddの2乗を指定します。精度行列の逆行列solve(lambda_truth_dd)を引数sigmaに指定することもできます。
観測データを散布図で確認しましょう。
# 観測データのデータフレーム sample_df <- tibble( x = x_nd[, 1], y = x_nd[, 2] ) # 観測モデルのデータフレーム model_df <- tibble( xy = point_df, density = mvtnorm::dmvnorm(x = xy, mean = mu_d, sigma = sigma_truth_dd^2), # 確率密度 ) %>% dplyr::select(density) %>% cbind(point_df, .) # 観測データの散布図を作成 ggplot() + geom_point(data = sample_df, aes(x = x, y = y)) + # 観測データ geom_contour(data = model_df, aes(x, y, z = density, color = ..level..)) + # 観測モデル geom_point(data = mu_df, aes(x = x, y = y), color = "red", shape = 3, size = 5) + # 平均パラメータ labs(title = "Multivariate Gaussian Distribution", subtitle = paste0("N=", N, ", mu=(", paste(round(mu_d, 1), collapse = ", "), ")", ", sigma=(", paste(round(sigma_truth_dd, 1), collapse = ", "), ")"), x = expression(x[1]), y = expression(x[2]), color = "density") # ラベル
観測モデルの分布と重ねて描画します。多次元ガウス分布の確率密度は、mvtnorm::dmvnorm()で計算します。等高線グラフgeom_contour()には、格子状の点を渡す必要があります。
では事後分布を求めていきます。
・事後分布
観測データを使って、事後分布のパラメータを計算します。
# 事後分布のパラメータを計算 W_hat_dd <- solve( (t(x_nd) - mu_d) %*% t(t(x_nd) - mu_d) + solve(W_dd) ) nu_hat <- N + nu
事後分布のパラメータ$\hat{\mathbf{W}},\ \hat{\nu}$をそれぞれW_hat_dd、nu_hatとします。計算式(更新式)はそれぞれ次の式です。
ただし上の実装例では、効率よく処理するために転置して計算しています。この式の通りに計算するには、次のような処理となります。
sum_dot_x <- matrix(0, D, D) for(n in 1:N) { dot_x <- matrix(x_nd[n, ] - mu_d) %*% t(x_nd[n, ] - mu_d) sum_dot_x <- sum_dot_x + dot_x } W_hat_dd <- solve(sum_dot_x + solve(W_dd))
for()によってN回加算する処理が$\sum_{n=1}^N$の計算に対応します。
精度行列$\boldsymbol{\Lambda}$の事後分布$p(\boldsymbol{\Lambda} | \mathbf{X})$のパラメータを計算できました。これまでのように、このパラメータを用いて事後分布(ウィシャート分布)を計算できます。しかし精度パラメータの分布はイメージしにくいため、精度パラメータの期待値$\mathbb{E}[\hat{\boldsymbol{\Lambda}}]$を用いた多次元ガウス分布を可視化することで、観測モデルと比較しましょう。
ちなみにウィシャート分布の確率密度はMCMCpack::dwish()で計算できます。
# 精度パラメータの期待値を計算 lambda_E_dd <- nu_hat * W_hat_dd
ウィシャート分布の期待値(2.89)より
で計算できます。
求めたパラメータを使って、確率密度を計算します。
# 事後分布の期待値を用いた分布を計算 posterior_df <- tibble( xy = point_df, density = mvtnorm::dmvnorm(x = xy, mean = mu_d, sigma = solve(lambda_E_dd)) # 確率密度 ) %>% dplyr::select(density) %>% cbind(point_df, .) head(posterior_df)
## x y density ## 1 -75.0 -50 5.682318e-06 ## 2 -74.5 -50 5.742626e-06 ## 3 -74.0 -50 5.802989e-06 ## 4 -73.5 -50 5.863398e-06 ## 5 -73.0 -50 5.923840e-06 ## 6 -72.5 -50 5.984302e-06
作図用にデータフレームとしておきます。
多次元ガウス分布を作図します。
# 作図 ggplot() + geom_contour(data = posterior_df, aes(x, y, z = density, color = ..level..)) + # 精度の期待値を用いた分布 geom_contour(data = model_df, aes(x, y, z = density, color = ..level..), alpha = 0.5, linetype = "dashed") + # 観測モデル geom_point(data = mu_df, aes(x = x, y = y), color = "red", shape = 3, size = 5) + # 平均値 labs(title = "Multivariate Gaussian Distribution", subtitle = paste0("N=", N, ", mu=(", paste(round(mu_d, 1), collapse = ", "), ")", ", E_sigma_hat=(", paste(round(sqrt(solve(lambda_E_dd)), 1), collapse = ", "), ")"), x = expression(x[1]), y = expression(x[2]), color = "density") # ラベル
観測データの散布図のときと同様に、観測モデル(真の分布)と重ねて描画します。
続いて予測分布を求めていきます。
・予測分布
事後分布のパラメータを使って、予測分布のパラメータを計算します。
# 予測分布のパラメータを計算 mu_s_d <- mu_d lambda_s_hat_dd <- (1 - D + nu_hat) * W_hat_dd nu_s_hat <- 1 - D + nu_hat
予測分布(多次元スチューデントのt分布)のパラメータ$\boldsymbol{\mu}_s,\ \hat{\boldsymbol{\Lambda}}_s,\ \hat{\nu}_s$をそれぞれmu_s_d、lambda_s_hat_dd、nu_s_hatとします。計算式(更新式)はそれぞれ次の式です。
求めたパラメータを使って、予測分布(の確率密度)を計算します。
# 予測分布を計算 predict_df <- cbind( point_df, density = mvnfast::dmvt( X = as.matrix(point_df), mu = mu_s_d, sigma = solve(lambda_s_hat_dd), df = nu_s_hat ) # 確率密度 ) head(predict_df)
## x y density ## 1 -75.0 -50 5.801869e-06 ## 2 -74.5 -50 5.861283e-06 ## 3 -74.0 -50 5.920743e-06 ## 4 -73.5 -50 5.980238e-06 ## 5 -73.0 -50 6.039755e-06 ## 6 -72.5 -50 6.099284e-06
多次元スチューデントのt分布の確率密度は、mvnfast::dmvt()で計算できます。dmvt()の第1引数には、複数のデータをマトリクスで渡すことができます。また平均引数muにはmu_s_d、スケール引数sigmaにはlambda_s_hat_ddの逆行列、自由度引数dfにはnu_s_hatを指定します。
mvtnorm::dmvt()でも計算できるはずですが、ちょっとよく分からなかった…。
予測分布を作図します。
# 作図 ggplot() + geom_contour(data = predict_df, aes(x, y, z = density, color = ..level..)) + # 予測分布 geom_contour(data = model_df, aes(x, y, z = density, color = ..level..), alpha = 0.5, linetype = "dashed") + # 観測モデル geom_point(data = mu_df, aes(x = x, y = y), color = "red", shape = 3, size = 5) + # 平均パラメータ labs(title = "Multivariate Student's t Distribution", subtitle = paste0("N=", N, ", mu_s=(", paste(round(mu_s_d, 1), collapse = ", "), ")", ", lambda_s_hat=(", paste(round(lambda_s_hat_dd, 1), collapse = ", "), ")", ", df=", nu_s_hat), x = expression(x[1]), y = expression(x[2]), color = "density") # ラベル
これまでと同様に作図できます。
・おまけ
gganimateパッケージを利用して、パラメータの推定値の推移を確認するためのgif画像を作成するコードです。
・コード(クリックで展開)
# 追加パッケージ library(gganimate) # データ数を指定 N <- 100 D <- 2 # (固定) # 観測モデルのパラメータを指定 mu_d <- c(25, 50) sigma_truth_dd <- matrix(c(50, 30, 30, 50), nrow = 2, ncol = 2) lambda_truth_dd <- solve(sigma_truth_dd^2) # 事前分布のパラメータを指定 nu <- 2 W_dd <- matrix(c(0.00005, 0, 0, 0.00005), nrow = 2, ncol = 2) # 精度パラメータの期待値を計算 lambda_E_dd <- nu * W_dd # 作図用の点を生成 x_vec <- seq(mu_d[1] - 3 * sigma_truth_dd[1, 1], mu_d[1] + 3 * sigma_truth_dd[1, 1], by = 1) y_vec <- seq(mu_d[2] - 3 * sigma_truth_dd[2, 2], mu_d[2] + 3 * sigma_truth_dd[2, 2], by = 1) point_df <- tibble( x = rep(x_vec, times = length(y_vec)), y = rep(y_vec, each = length(x_vec)) ) mu_df <- tibble( x = mu_d[1], y = mu_d[2] ) # 観測モデルを計算 model_df <- tibble( xy = point_df, density = mvtnorm::dmvnorm(x = xy, mean = mu_d, sigma = sigma_truth_dd^2), # 確率密度 ) %>% dplyr::select(density) %>% cbind(point_df, .) # 事前分布の期待値を用いた分布を計算 posterior_df <- tidyr::tibble( xy = point_df, density = mvtnorm::dmvnorm(x = xy, mean = mu_d, sigma = solve(lambda_E_dd)), # 確率密度 iteration = 0 # 試行回数 ) %>% dplyr::select(density, iteration) %>% cbind(point_df, .) # 予測分布のパラメータを計算 mu_s_d <- mu_d lambda_s_dd <- (1 - D + nu) * W_dd nu_s <- 1 - D + nu # 初期値による予測分布を計算 predict_df <- cbind( point_df, density = mvnfast::dmvt( X = as.matrix(point_df), mu = mu_s_d, sigma = solve(lambda_s_dd), df = nu_s ), # 確率密度 iteration = 0 # 試行回数 ) # ベイズ推論 for(i in 1:N) { # 2次元ガウス分布に従うデータを生成 x_nd <- mvtnorm::rmvnorm(n = 1, mean = mu_d, sigma = sigma_truth_dd^2) # 観測データを記録 if(i > 1) { # 初回以外 # オブジェクトを結合 x_mat <- rbind(x_mat, x_nd) sample_df <- tibble( x = x_mat[, 1], y = x_mat[, 2], iteration = i ) %>% rbind(sample_df, .) } else if(i == 1){ # 初回 # オブジェクトを作成 x_mat <- x_nd sample_df <- tibble( x = x_mat[, 1], y = x_mat[, 2], iteration = i ) } # 事後分布のパラメータを更新 W_dd <- solve( (t(x_nd) - mu_d) %*% t(t(x_nd) - mu_d) + solve(W_dd) ) nu <- 1 + nu # 精度パラメータの期待値を計算 lambda_E_dd <- nu * W_dd # 事後分布の期待値を用いた分布を計算 tmp_posterior_df <- tidyr::tibble( xy = point_df, density = mvtnorm::dmvnorm(x = xy, mean = mu_d, sigma = solve(lambda_E_dd)), # 確率密度 iteration = i # 試行回数 ) %>% dplyr::select(density, iteration) %>% cbind(point_df, .) # 予測分布のパラメータを更新 mu_s_d <- mu_d lambda_s_dd <- (1 - D + nu) * W_dd nu_s <- 1 - D + nu # 予測分布を計算 tmp_predict_df <- cbind( point_df, density = mvnfast::dmvt( X = as.matrix(point_df), mu = mu_s_d, sigma = solve(lambda_s_dd), df = nu_s ), # 確率密度 iteration = i # 試行回数 ) # 推論結果を結合 posterior_df <- rbind(posterior_df, tmp_posterior_df) predict_df <- rbind(predict_df, tmp_predict_df) # 動作確認 print(i) } # 事後分布の期待値を用いた分布を作図 posterior_graph <- ggplot() + geom_contour(data = posterior_df, aes(x, y, z = density, color = ..level..)) + # 精度の期待値を用いた分布 geom_point(data = sample_df, aes(x = x, y = y)) + # 観測データ geom_contour(data = model_df, aes(x, y, z = density, color = ..level..), alpha = 0.5, linetype = "dashed") + # 観測モデル geom_point(data = mu_df, aes(x = x, y = y), color = "red", shape = 3, size = 5) + # 平均パラメータ transition_manual(iteration) + # フレーム labs(title = "Multivariate Gaussian Distribution", subtitle = "N={current_frame}", x = expression(x[1]), y = expression(x[2]), color = "density") # ラベル # gif画像を作成 animate(posterior_graph, nframes = N + 1, fps = 5) # 予測分布を作図 predict_graph <- ggplot() + geom_contour(data = predict_df, aes(x, y, z = density, color = ..level..)) + # 予測分布 geom_point(data = sample_df, aes(x = x, y = y)) + # 観測データ geom_contour(data = model_df, aes(x, y, z = density, color = ..level..), alpha = 0.5, linetype = "dashed") + # 観測モデル geom_point(data = mu_df, aes(x = x, y = y), color = "red", shape = 3, size = 5) + # 平均パラメータ transition_manual(iteration) + # フレーム labs(title = "Multivariate Student's t Distribution", subtitle = "N={current_frame}", x = expression(x[1]), y = expression(x[2]), color = "density") # ラベル # gif画像を作成 animate(predict_graph, nframes = N + 1, fps = 5)
異なる点のみを簡単に解説します。
各データによってどのように学習する(推定値が変化する)のかを確認するため、こちらはfor()で1データずつ処理します。よって観測データ数Nがイタレーション数になります。
一度の処理で事後分布のパラメータを計算するのではなく、事前分布のパラメータに対して繰り返し観測データの情報を与えることで更新(上書き)していきます。$\sum_{n=1}^N$の計算は、N回繰り返し計算することで実行されます。
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
やっぱり今はベイズ推論が楽しい！ところでスケール行列って何？
mvnfastパッケージにも多次元ガウス分布関連の関数がありますが、書き換えるのが面倒だったので止めました。
【次節の内容】続く