はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。基本的な内容は「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」になります。「数式」と「プログラム」から理解するのが目標です。
この記事は、3.3.3項の内容です。「尤度関数を平均と精度が未知の1次元ガウス分布(正規分布)」、「事前分布をガウス・ガンマ分布」とした場合の「パラメータの事後分布」と「未観測値の予測分布」の計算をR言語で実装します。
省略してある内容等ありますので、本とあわせて読んでください。初学者な自分が理解できるレベルまで落として書き下していますので、分かる人にはかなりくどくなっています。同じような立場の人のお役に立てれば幸いです。
【数式読解編】
【他の節の内容】
【この節の内容】
・Rでやってみよう
人工的に生成したデータを用いて、ベイズ推論を行ってみましょう。
利用するパッケージを読み込みます。
# 3.3.2項で利用するパッケージ library(tidyverse)
・モデルの構築
まずは、モデルを設定します。
尤度(ガウス分布)のパラメータを設定します。
# 真のパラメータwを指定 mu_truth <- 25 lambda_truth <- 0.01
平均パラメータ$\mu$をmu_truth、精度パラメータ$\lambda$をlambda_truthとして値を指定します。精度は分散の逆数なので、値が大きいほど散らばり具合が小さくなり、逆数の平方根が標準偏差$\sigma = \sqrt{\lambda^{-1}}$になります。この2つが真のパラメータであり、この値を求めるのがここでの目的です。
尤度の確率密度を計算して、作図用のデータフレームを作成します。
# 作図用のxの値を設定 x_line <- seq( mu_truth - 4 * sqrt(1 / lambda_truth), mu_truth + 4 * sqrt(1 / lambda_truth), length.out = 1000 ) # 尤度を計算:式(2.64) model_df <- tibble( x = x_line, # x軸の値 ln_C_N = - 0.5 * (log(2 * pi) - log(lambda_truth)), # 正規化項(対数) density = exp(ln_C_N - 0.5 * lambda_truth * (x - mu_truth)^2) # 確率密度 #C_N = 1 / sqrt(2 * pi / lambda_truth), # 正規化項 #density = C_N * exp(- 0.5 * lambda_truth * (x - mu_truth)^2) # 確率密度 #density = dnorm(x = x, mean = mu_truth, sd = sqrt(1 / lambda_truth)) # 確率密度 )
作図用に、ガウス分布に従う変数$x_n$がとり得る値をseq()で作成してx_lineとします。この例では、平均値を中心に標準偏差の4倍を範囲とします。length.out引数を使うと指定した要素数で等間隔に切り分けます。by引数を使うと切り分ける間隔を指定できます。処理が重い場合は、この値を調整してください。
x_lineの値ごとに確率密度を計算します。1次元ガウス分布の確率密度は、対数をとった定義式
で計算して、最後にexp()をします。ガウス分布の確率密度関数dnorm()でも計算できます。
作成したデータフレームを確認しましょう。
# 確認 head(model_df)
## # A tibble: 6 x 3 ## x ln_C_N density ## <dbl> <dbl> <dbl> ## 1 -15 -3.22 0.0000134 ## 2 -14.9 -3.22 0.0000138 ## 3 -14.8 -3.22 0.0000143 ## 4 -14.8 -3.22 0.0000147 ## 5 -14.7 -3.22 0.0000152 ## 6 -14.6 -3.22 0.0000157
ggplot2パッケージを利用して作図するには、データフレームを渡す必要があります。
尤度を作図します。
# 尤度を作図 ggplot(model_df, aes(x = x, y = density)) + geom_line(color = "purple") + # 尤度 labs(title = "Gaussian Distribution", subtitle = paste0("mu=", mu_truth, ", lambda=", lambda_truth))
真のパラメータを求めることは、この真の分布を求めることを意味します。
・データの生成
続いて、構築したモデルに従って観測データ$\mathbf{X} = \{x_1, x_2, \cdots, x_N\}$を生成します。
ガウス分布に従う$N$個のデータをランダムに生成します。
# (観測)データ数を指定 N <- 50 # ガウス分布に従うデータを生成 x_n <- rnorm(n = N, mean = mu_truth, sd = sqrt(1 / lambda_truth))
生成するデータ数$N$をNとして値を指定します。
ガウス分布に従う乱数は、rnorm()で生成できます。試行回数の引数nにN、平均の引数meanにmu_truth、標準偏差の引数sdにsqrt(1 / lambda_truth)を指定します。生成したN個のデータをx_nとします。
観測したデータを確認しましょう。
# 確認 summary(x_n)
## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 1.354 19.451 24.297 24.630 29.685 49.632
観測データをヒストグラムでも確認します。
# 観測データのヒストグラムを作図 tibble(x = x_n) %>% ggplot(aes(x = x)) + geom_histogram(binwidth = 1) + # 観測データ labs(title = "Gaussian Distribution", subtitle = paste0("N=", N, ", mu=", mu_truth, ", sigma=", round(sqrt(1 / lambda_truth), 1)))
データ数が十分に大きいと、分布の形状が真の分布に近づきます。
・事前分布の設定
次に、尤度に対する共役事前分布を設定します。
$\mu$の事前分布(ガウス分布)のパラメータ(超パラメータ)を設定します。
# muの事前分布のパラメータを指定 m <- 0 beta <- 1
ガウス分布の平均パラメータ$m$と精度パラメータの係数$\beta$をそれぞれm, betaとして値を指定します。
また、$\lambda$の事前分布(ガンマ分布)のパラメータ(超パラメータ)も設定します。
# lambdaの事前分布のパラメータを指定 a <- 1 b <- 1
ガンマ分布のパラメータ$a,\ b$をそれぞれa, bとして正の実数値を指定します。
$\mu$の事前分布の計算用に、$\lambda$の事前分布から$\lambda$の期待値を求めます。
# lambdaの期待値を計算:式(2.59) E_lambda <- a / b
ガンマ分布の期待値は
で計算します。
設定したパラメータを使って、$\mu$の事前分布の確率密度を計算します。
# 作図用のmuの値を設定 mu_line <- seq(mu_truth - 50, mu_truth + 50, length.out = 1000) # muの事前分布を計算:式(2.64) prior_mu_df <- tibble( mu = mu_line, # x軸の値 ln_C_N = - 0.5 * (log(2 * pi) - log(beta * E_lambda)), # 正規化項(対数) density = exp(ln_C_N - 0.5 * beta * E_lambda * (mu - m)^2) # 確率密度 #density = dnorm(x = mu, mean = m, sd = sqrt(1 / beta / E_lambda)) # 確率密度 )
作図用に、ガウス分布に従う変数$\mu$がとり得る値をseq()で作成してmu_lineとします。この例では、真の値を中心に指定した範囲とします(本当は自動で良い感じの範囲になるようにしたかった)。
尤度のときと同様にして、確率密度を計算します。
計算結果は次のようになります。
# 確認 head(prior_mu_df)
## # A tibble: 6 x 3 ## mu ln_C_N density ## <dbl> <dbl> <dbl> ## 1 -25 -0.919 7.65e-137 ## 2 -24.9 -0.919 9.30e-136 ## 3 -24.8 -0.919 1.12e-134 ## 4 -24.7 -0.919 1.33e-133 ## 5 -24.6 -0.919 1.57e-132 ## 6 -24.5 -0.919 1.83e-131
$\mu$の事前分布を作図します。
# muの事前分布を作図 ggplot(prior_mu_df, aes(x = mu, y = density)) + geom_line(color = "purple") + # muの事前分布 labs(title = "Gaussian Distribution", subtitle = paste0("m=", m, ", beta=", beta, ', E_lambda=', E_lambda), x = expression(mu))
m, betaやa, bの値を変更することで、ガウス分布におけるパラメータと形状の関係を確認できます。
続いて、$\lambda$の事前分布の確率密度も計算します。
# 作図用のlambdaの値を設定 lambda_line <- seq(0, 4 * lambda_truth, length.out = 1000) # lambdaの事前分布を計算:式(2.56) prior_lambda_df <- tibble( lambda = lambda_line, # x軸の値 ln_C_Gam = a * log(b) - lgamma(a), # 正規化項(対数) density = exp(ln_C_Gam + (a - 1) * log(lambda) - b * lambda) # 確率密度 #density = dgamma(x = lambda, shape = a, rate = b) # 確率密度 )
作図用に、ガンマ分布に従う変数$\lambda$がとり得る値をseq()で作成してlambda_lineとします。この例では、0から真の値の4倍までを範囲とします。
lambda_lineの値ごとに確率密度を計算します。ガンマ分布の確率密度は、対数をとった定義式
で計算して、最後にexp()をします。ここで、$\Gamma(\cdot)$はガンマ関数です。
ガンマ関数の計算はgamma()で行えますが、値が大きくなると発散してしまします。そこで、対数をとったガンマ関数lgamma()で計算した後にexp()で戻します。
ガンマ分布の確率密度は、dgamma()でも計算できます。
計算結果は次のようになります。
# 確認 head(prior_lambda_df)
## # A tibble: 6 x 3 ## lambda ln_C_Gam density ## <dbl> <dbl> <dbl> ## 1 0 0 NaN ## 2 0.0000400 0 1.00 ## 3 0.0000801 0 1.00 ## 4 0.000120 0 1.00 ## 5 0.000160 0 1.00 ## 6 0.000200 0 1.00
$\lambda$の事前分布を作図します。
# lambdaの事前分布を作図 ggplot(prior_lambda_df, aes(x = lambda, y = density)) + geom_line(color = "purple") + # lambdaの事前分布 labs(title = "Gamma Distribution", subtitle = paste0("a=", a, ", b=", b))
(lambda_lineの範囲をもっと広くして)a, bの値を変更することで、ガウス分布におけるパラメータと形状の関係を確認できます。
・事後分布の計算
観測データ$\mathbf{X}$からパラメータ$\mu$と$\lambda$それぞれ事後分布を求めます(パラメータ$\mu,\ \lambda$を分布推定します)。
観測データx_nを用いて、$\mu$の事後分布(ガウス分布)のパラメータを計算します。
# muの事後分布のパラメータを計算:式(3.83) beta_hat <- N + beta m_hat <- (sum(x_n) + beta * m) / beta_hat
$\mu$の事後分布のパラメータは
で計算して、結果をbeta_hat, m_hatとします。
# 確認 beta_hat; m_hat
## [1] 51 ## [1] 24.14679
続いて、$\lambda$の事後分布(ガンマ分布)のパラメータも計算します。
# lambdaの事後分布のパラメータを計算:式(3.88) a_hat <- 0.5 * N + a b_hat <- 0.5 * (sum(x_n^2) + beta * m^2 - beta_hat * m_hat^2) + b
$\lambda$の事後分布のパラメータは
で計算して、結果をa_hat, b_hatとします。
# 確認 a_hat; b_hat
## [1] 26 ## [1] 2595.884
$\mu$の事後分布の計算用に、$\lambda$の事後分布から$\hat{\lambda}$の期待値を求めます。
# lambdaの期待値を計算:式(2.59) E_lambda_hat <- a_hat / b_hat # 確認 E_lambda_hat
## [1] 0.01001586
求めたパラメータを使って、$\mu$の事後分布の確率密度を計算します。
# muの事後分布を計算:式(2.64) posterior_mu_df <- tibble( mu = mu_line, # x軸の値 ln_C_N = - 0.5 * (log(2 * pi) - log(beta_hat * E_lambda_hat)), # 正規化項(対数) density = exp(ln_C_N - 0.5 * beta_hat * E_lambda_hat * (mu - m_hat)^2) # 確率密度 #density = dnorm(mu, mean = m_hat, sd = sqrt(1 / beta_hat / E_lambda_hat)) # 確率密度 )
更新した超パラメータm_hat, lambda_mu_hatを用いて、事前分布のときと同様にして計算します。
計算結果は次のようになります。
# 確認 head(posterior_mu_df)
## # A tibble: 6 x 3 ## mu ln_C_N density ## <dbl> <dbl> <dbl> ## 1 -25 -1.25 3.44e-269 ## 2 -24.9 -1.25 4.23e-268 ## 3 -24.8 -1.25 5.18e-267 ## 4 -24.7 -1.25 6.32e-266 ## 5 -24.6 -1.25 7.66e-265 ## 6 -24.5 -1.25 9.24e-264
$\mu$の事後分布を作図します。
# muの事後分布を作図 ggplot(posterior_mu_df, aes(x = mu, y = density)) + geom_line(color = "purple") + # muの事後分布 geom_vline(aes(xintercept = mu_truth), color = "red", linetype = "dashed") + # 真のmu labs(title = "Gaussian Distribution", subtitle = paste0("N=", N, ", m_hat=", round(m_hat, 1), ", beta_hat=", beta_hat, ", E_lambda=", round(E_lambda_hat, 3)), x = expression(mu))
パラメータ$\mu$の真の値付近をピークとする分布を推定できています。
続いて、$\lambda$の事後分布の確率密度も計算します。
# lambdaの事後分布を計算:式(2.56) posterior_lambda_df <- tibble( lambda = lambda_line, # x軸の値 ln_C_Gam = a_hat * log(b_hat) - lgamma(a_hat), # 正規化項(対数) density = exp(ln_C_Gam + (a_hat - 1) * log(lambda) - b_hat * lambda) # 確率密度 #density = dgamma(x = lambda, shape = a_hat, rate = b_hat) # 確率密度 )
更新した超パラメータa_hat, b_hatを用いて、事前分布のときと同様にして計算します。
計算結果は次のようになります。
# 確認 head(posterior_lambda_df)
## # A tibble: 6 x 3 ## lambda ln_C_Gam density ## <dbl> <dbl> <dbl> ## 1 0 146. 0. ## 2 0.0000400 146. 3.96e-47 ## 3 0.0000801 146. 1.20e-39 ## 4 0.000120 146. 2.73e-35 ## 5 0.000160 146. 3.27e-32 ## 6 0.000200 146. 7.79e-30
$\lambda$の事後分布を作図します。
# lambdaの事後分布を作図 ggplot(posterior_lambda_df, aes(x = lambda, y = density)) + geom_line(color = "purple") + # lambdaの事後分布 geom_vline(aes(xintercept = lambda_truth), color = "red", linetype = "dashed") + # 真のlambda labs(title = "Gamma Distribution", subtitle = paste0("N=", N, ", a_hat=", a_hat, ", b_hat=", round(b_hat, 1)), x = expression(lambda))
パラメータ$\lambda$の真の値付近をピークとする分布を推定できています。
・予測分布の計算
最後に、観測データ$\mathbf{X}$から未観測のデータ$x_{*}$の予測分布を求めます。
$\mu$の事後分布のパラメータm_hat, beta_hatと$\lambda$の事後分布のパラメータa_hat, b_hat、または観測データx_nと$\mu$の事前分布のパラメータm, beta、$\lambda$の事前分布のパラメータa, bを用いて、予測分布(スチューデントのt分布)のパラメータを計算します。
# 予測分布のパラメータを計算:式(3.95') mu_s_hat <- m_hat lambda_s_hat <- beta_hat * a_hat / (1 + beta_hat) / b_hat nu_s_hat <- 2 * a_hat # 予測分布のパラメータを計算:式(3.95') mu_s_hat <- (sum(x_n) + beta * m) / (N + beta) numer_lambda <- (N + beta) * (N / 2 + a) denom_lambda <- (N + 1 + beta) * ((sum(x_n^2) + beta * m^2 - beta_hat * m_hat^2) / 2 + b) lambda_s_hat <- numer_lambda / denom_lambda nu_s_hat <- N + 2 * a
予測分布のパラメータは
で計算して、結果をmu_s_hat, lambda_s_hat, nu_s_hatとします。$\mu_s$は平均、$\lambda_s$は何？、$\nu_s$は自由度です。
それぞれ上の式だと、事後分布のパラメータmu_hat, beta_hat、a_hat, b_hatを使って計算できます。下の式だと、観測データx_nと事前分布のパラメータmu, beta、a, bを使って計算できます(嘘ですね$\hat{\beta}$と$\hat{m}$が残ってますねでもめんどいよう…)。
# 確認 mu_s_hat; lambda_s_hat; nu_s_hat
## [1] 24.14679 ## [1] 0.009823245 ## [1] 52
$\mathbf{X}$から$\hat{\mu}_s,\ \hat{\lambda}_s,\ \hat{\nu}_s$を学習しているのが式からも分かります。
求めたパラメータを使って、予測分布を計算します。
# 予測分布を計算:式(3.76) predict_df <- tibble( x = x_line, # x軸の値 ln_C_St = lgamma(0.5 * (nu_s_hat + 1)) - lgamma(0.5 * nu_s_hat), # 正規化項(対数) ln_term1 = 0.5 * log(lambda_s_hat / pi / nu_s_hat), ln_term2 = - 0.5 * (nu_s_hat + 1) * log(1 + lambda_s_hat / nu_s_hat * (x - mu_s_hat)^2), density = exp(ln_C_St + ln_term1 + ln_term2) # 確率密度 )
x_lineの値ごとに確率密度を計算します。1次元スチューデントのt分布の確率密度は、対数をとった定義式
で計算して、最後にexp()をします。この$\pi$は円周率です。(t分布の確率密度関数dt()でもやりたかったけど$\lambda_s$用の引数が分からなかった。そもそも$\lambda_s$って何？)
計算結果は次のようになります。
# 確認 head(predict_df)
## # A tibble: 6 x 5 ## x ln_C_St ln_term1 ln_term2 density ## <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 -15 1.62 -4.86 -6.74 0.0000466 ## 2 -14.9 1.62 -4.86 -6.71 0.0000478 ## 3 -14.8 1.62 -4.86 -6.69 0.0000490 ## 4 -14.8 1.62 -4.86 -6.66 0.0000502 ## 5 -14.7 1.62 -4.86 -6.64 0.0000514 ## 6 -14.6 1.62 -4.86 -6.62 0.0000527
予測分布を尤度と重ねて作図します。
# 予測分布を作図 ggplot() + geom_line(data = predict_df, aes(x = x, y = density), color = "purple") + # 予測分布 geom_line(data = model_df, aes(x = x, y = density), color = "red", linetype = "dashed") + # 真の分布 labs(title = "Student's t Distribution", subtitle = paste0("N=", N, ", mu_s_hat=", round(mu_s_hat, 1), ", lambda_s_hat=", round(lambda_s_hat, 3), ", nu_s_hat=", nu_s_hat))
観測データが増えると、予測分布が真の分布に近づきます。
・おまけ：アニメーションで推移の確認
gganimateパッケージを利用して、事後分布と予測分布の推移をアニメーション(gif画像)で確認するためのコードです。
・コード(クリックで展開)
異なる点のみを簡単に解説します。
# 利用するパッケージ library(tidyverse) library(gganimate)
・モデルの設定
# 真のパラメータを指定 mu_truth <- 25 lambda_truth <- 0.01 # muの事前分布のパラメータを指定 m <- 0 beta <- 1 # lambdaの事前分布のパラメータを指定 a <- 1 b <- 1 # lambdaの期待値を計算:式(2.59) E_lambda <- a / b # 作図用のmuの値を設定 mu_line <- seq(mu_truth - 30, mu_truth + 30, length.out = 1000) # muの事前分布(ガウス分布)を計算:式(2.64) posterior_mu_df <- tibble( mu = mu_line, # x軸の値 density = dnorm(mu, mean = m, sd = sqrt(1 / beta / E_lambda)), # 確率密度 label = as.factor( paste0( "N=", 0, ", m=", round(m, 1), ", beta=", beta, ", E_lambda=", round(E_lambda, 3) ) ) # パラメータ ) # 作図用のlambdaの値を設定 lambda_line <- seq(0, 4 * lambda_truth, length.out = 1000) # lambdaの事前分布(ガンマ分布)を計算:式(2.56) posterior_lambda_df <- tibble( lambda = lambda_line, # x軸の値 density = dgamma(x = lambda, shape = a, rate = b), # 確率密度 label = as.factor(paste0("N=", 0, ", a=", a, ", b=", round(b, 1))) # パラメータ ) # 初期値による予測分布のパラメータを計算:式(3.95) mu_s <- m lambda_s <- beta * a / (1 + beta) / b nu_s <- 2 * a # 作図用のxの値を設定 x_line <- seq( mu_truth - 4 * sqrt(1 / lambda_truth), mu_truth + 4 * sqrt(1 / lambda_truth), length.out = 1000 ) # 初期値による予測分布(スチューデントのt分布)を計算:式(3.76) predict_df <- tibble( x = x_line, # x軸の値 ln_C_St = lgamma(0.5 * (nu_s + 1)) - lgamma(0.5 * nu_s), # 正規化項(対数) ln_term1 = 0.5 * log(lambda_s / pi / nu_s), ln_term2 = - 0.5 * (nu_s + 1) * log(1 + lambda_s / nu_s * (x - mu_s)^2), density = exp(ln_C_St + ln_term1 + ln_term2), # 確率密度 label = as.factor( paste0( "N=", 0, ", nu_s=", nu_s, ", mu_s=", round(mu_s, 1), ", lambda_s=", round(lambda_s, 3) ) ) # パラメータ )
各試行の結果を同じデータフレームに格納していく必要があります。$\mu$の事後分布をposterior_mu_df、$\lambda$の事後分布をposterior_lambda_df、予測分布をpredict_dfとして、初期値の結果を持つように作成しておきます。
・推論処理
# データ数(試行回数)を指定 N <- 100 # 観測データの受け皿を初期化 x_n <- rep(0, N) # ベイズ推論 for(n in 1:N){ # ガウス分布に従うデータを生成 x_n[n] <- rnorm(n = 1, mean = mu_truth, sd = sqrt(1 / lambda_truth)) # muの事後分布のパラメータを更新:式(3.83) beta_old <- beta beta <- 1 + beta m_old <- m m <- (x_n[n] + beta_old * m) / beta # lambdaの事後分布のパラメータを更新:式(3.88) a <- 1 / 2 + a b <- (x_n[n]^2 + beta_old * m_old^2 - beta * m^2) / 2 + b # lambdaの期待値を計算:式(2.59) E_lambda <- a / b # muの事後分布(ガウス分布)を計算:式(2.64) tmp_posterior_mu_df <- tibble( mu = mu_line, # x軸の値 density = dnorm(mu, mean = m, sd = sqrt(1 / beta / E_lambda)), # 確率密度 label = as.factor( paste0( "N=", n, ", m_hat=", round(m, 1), ", beta_hat=", beta, ", E_lambda_hat=", round(E_lambda, 3) ) ) # パラメータ ) # lambdaの事後分布(ガンマ分布)を計算:式(2.56) tmp_posterior_lambda_df <- tibble( lambda = lambda_line, # x軸の値 density = dgamma(x = lambda, shape = a, rate = b), # 確率密度 label = as.factor(paste0("N=", n, ", a_hat=", a, ", b_hat=", round(b, 1))) # パラメータ ) # 予測分布のパラメータを更新:式(3.95) mu_s <- m lambda_s <- beta * a / (1 + beta) / b nu_s <- 2 * a # 予測分布(スチューデントのt分布)を計算:式(3.76) tmp_predict_df <- tibble( x = x_line, # x軸の値 ln_C_St = lgamma(0.5 * (nu_s + 1)) - lgamma(0.5 * nu_s), # 正規化項(対数) ln_term1 = 0.5 * log(lambda_s / pi / nu_s), ln_term2 = - 0.5 * (nu_s + 1) * log(1 + lambda_s / nu_s * (x - mu_s)^2), density = exp(ln_C_St + ln_term1 + ln_term2), # 確率密度 label = as.factor( paste0( "N=", n, ", nu_s_hat=", nu_s, ", mu_s_hat=", round(mu_s, 1), ", lambda_s_hat=", round(lambda_s, 3) ) ) # パラメータ ) # 推論結果を結合 posterior_mu_df <- rbind(posterior_mu_df, tmp_posterior_mu_df) posterior_lambda_df <- rbind(posterior_lambda_df, tmp_posterior_lambda_df) predict_df <- rbind(predict_df, tmp_predict_df) }
観測された各データによってどのように学習する(分布が変化する)のかを確認するため、for()で1データずつ処理します。よって、データ数Nがイタレーション数になります。
超パラメータに関して、$\hat{m},\ \hat{\beta}$に対応するm_hat, beta_hatや$\hat{a},\ \hat{b}$に対応するa_hat, b_hatを新たに作るのではなく、m_hat, beta_hatとa, bをイタレーションごとに更新(上書き)していきます。
それに伴い、事後分布のパラメータの計算式(3.83)(3.88)の$\sum_{n=1}^N x_n$や$N$といった計算は、ループ処理によってN回繰り返しx_n[n]や1を加えることで行います。n回目のループ処理のときには、n-1回分のx_n[n]や1が既に各パラメータに加えられているわけです。
・$\mu$の事後分布の推移
# muの事後分布を作図 posterior_mu_graph <- ggplot(posterior_mu_df, aes(x = mu, y = density)) + geom_line(color = "purple") + # muの事後分布 geom_vline(aes(xintercept = mu_truth), color = "red", linetype = "dashed") + # 真のmu gganimate::transition_manual(label) + # フレーム labs(title = "Gaussian Distribution", subtitle = "{current_frame}", x = expression(mu)) # gif画像を出力 gganimate::animate(posterior_mu_graph, nframes = N + 1, fps = 10)
各フレームの順番を示す列(label)をgganimate::transition_manual()に指定します。分布の推移と共にパラメータの値を表示するようにlabel列を作成していますが、ややこしければas.factor(paste0("iter=", n))として試行回数だけ表示するだけでもそれっぽくなります。
初期値(事前分布)を含むため、フレーム数の引数nframesはN + 1です。
・$\lambda$の事後分布の推移
# lambdaの事後分布を作図 posterior_lambda_graph <- ggplot(posterior_lambda_df, aes(x = lambda, y = density)) + geom_line(color = "purple") + # lambdaの事後分布 geom_vline(aes(xintercept = lambda_truth), color = "red", linetype = "dashed") + # 真のlambda gganimate::transition_manual(label) + # フレーム labs(title = "Gamma Distribution", subtitle = "{current_frame}", x = expression(lambda)) # gif画像を出力 gganimate::animate(posterior_lambda_graph, nframes = N + 1, fps = 10)
・予測分布の推移
# 尤度を計算:式(2.64) model_df <- tibble( x = x_line, # x軸の値 density = dnorm(x = x, mean = mu_truth, sd = sqrt(1 / lambda_truth)) # 確率密度 ) # 予測分布を作図 predict_graph <- ggplot() + geom_line(data = predict_df, aes(x = x, y = density), color = "purple") + # 予測分布 geom_line(data = model_df, aes(x = x, y = density), color = "red", linetype = "dashed") + # 真の分布 gganimate::transition_manual(label) + # フレーム labs(title = "Student's t Distribution", subtitle = "{current_frame}") # gif画像を出力 gganimate::animate(predict_graph, nframes = N + 1, fps = 10)
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
- 2021.03.18：加筆修正の際に、数式読解編とRで実装編に記事を分割しました。
修正作業って大変だけど、勘違いしてる1年前の自分に教えてあげてるんだって思うと優しい気持ちになれる(？)
【次節の内容】