はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」によって、「数式」と「プログラム」から理解するのが目標です。
省略してある内容等ありますので、本とあわせて読んでください。
この記事では、「尤度関数を1次元ガウス分布」、「事前分布をガウス・ガンマ分布」とした場合の「パラメータの事後分布」と「未観測値の予測分布」の計算をR言語で実装します。
【数式読解編】
【前節の内容】
【他の節の内容】
【この節の内容】
3.3.3 1次元ガウス分布の学習と予測：平均・精度が未知の場合
尤度関数を1次元ガウス分布、事前分布をガウス・ガンマ分布とするモデルに対するベイズ推論を実装します。人工的に生成したデータを用いて、ガウス分布の平均パラメータを推定し、また未観測データに対する予測分布を求めます。
ガウス分布については「1次元ガウス分布の定義式 - からっぽのしょこ」、ガウス・ガンマ分布については「ガウス-ガンマ分布の定義式 - からっぽのしょこ」、t分布については「1次元スチューデントのt分布の定義式 - からっぽのしょこ」を参照してください。
利用するパッケージを読み込みます。
# 利用パッケージ library(tidyverse) library(gganimate) library(LaplacesDemon)
この記事では、基本的にパッケージ名::関数名()の記法を使うので、パッケージを読み込む必要はありません。ただし、作図コードがごちゃごちゃしないようにパッケージ名を省略しているため、ggplot2は読み込む必要があります。
magrittrパッケージのパイプ演算子%>%ではなく、ネイティブパイプ(ベースパイプ)演算子|>を使っています。%>%に置き換えても処理できます。
学習の推移をアニメーション(gif画像)で確認するのにgganimateパッケージを利用します。不要であれば省略してください。
ベイズ推論の実装
まずは、モデルを設定して、データを生成します。生成したデータを用いて、事後分布のパラメータを計算します。さらに、事後分布のパラメータを用いて、予測分布のパラメータを計算します。
生成分布の設定
データ生成分布(真の分布)として、ガウス分布$\mathcal{N}(x | \mu, \lambda^{-1})$を設定します。
ガウス分布のグラフ作成については「【R】1次元ガウス分布の作図 - からっぽのしょこ」を参照してください。
真の分布(ガウス分布)のパラメータ$\mu, \lambda$を設定します。
# 真の平均パラメータを指定 mu_truth <- 25 # 真の精度パラメータを指定 lambda_truth <- 0.01 sqrt(1 / lambda_truth) # 標準偏差
## [1] 10
平均パラメータ(実数)$\mu$をmu_truth、精度パラメータ(正の実数)$\lambda$をlambda_truthとして値を指定します。精度は分散の逆数なので、値が大きいほど散らばり具合が小さくなり、逆数の平方根が標準偏差$\sigma = \frac{1}{\sqrt{\lambda}}$です。この2つが真のパラメータであり、この値を求めるのがここでの目的です。
真の分布を計算して、作図用のデータフレームを作成します。
# グラフ用のxの値を作成 x_vec <- seq( mu_truth - 1/sqrt(lambda_truth) * 4, mu_truth + 1/sqrt(lambda_truth) * 4, length.out = 501 ) # 真の分布を計算:式(2.64) model_df <- tibble::tibble( x = x_vec, # 確率変数 dens = dnorm(x = x_vec, mean = mu_truth, sd = 1/sqrt(lambda_truth)) # 確率密度 ) model_df
## # A tibble: 501 × 2 ## x dens ## <dbl> <dbl> ## 1 -15 0.0000134 ## 2 -14.8 0.0000143 ## 3 -14.7 0.0000152 ## 4 -14.5 0.0000162 ## 5 -14.4 0.0000173 ## 6 -14.2 0.0000184 ## 7 -14.0 0.0000196 ## 8 -13.9 0.0000208 ## 9 -13.7 0.0000221 ## 10 -13.6 0.0000236 ## # … with 491 more rows
ガウス分布の変数がとり得る値(実数)$x$作成して、値ごとに確率密度を計算します。この例では、平均を中心に標準偏差の4倍を範囲とします。
ガウス分布の確率密度は、dnorm()で計算できます。確率変数の引数xにx_vec、平均の引数meanにmu_truth、標準偏差の引数sdにlambda_truthの平方根の逆数を指定します。
真の分布のグラフを作成します。
# 真の分布を作図 ggplot() + geom_line(data = model_df, mapping = aes(x = x, y = dens, color = "model"), size = 1) + # 真の分布 scale_color_manual(breaks = "model", values = "purple", labels = "true model", name = "") + # 線の色:(凡例表示用) labs(title = "Gaussian Distribution", subtitle = parse(text = paste0("list(mu==", mu_truth, ", lambda==", lambda_truth, ")")), x = "x", y = "density")
真のパラメータを求めることは、真の分布を求めることを意味します。
データの生成
続いて、構築したモデルに従って観測データ$\mathbf{X} = \{x_1, x_2, \cdots, x_N\}$を生成します。
ガウス分布の乱数生成については「【R】1次元ガウス分布の乱数生成 - からっぽのしょこ」を参照してください。
ガウス分布に従う$N$個のデータをランダムに生成します。
# (観測)データ数を指定 N <- 50 # ガウス分布に従うデータを生成 x_n <- rnorm(n = N, mean = mu_truth, sd = 1/sqrt(lambda_truth)) head(x_n)
## [1] 22.67110 22.96440 22.71595 26.66057 39.23895 33.45446
生成するデータ数$N$を指定します。
ガウス分布に従う乱数は、rnorm()で生成できます。試行回数の引数nにN、平均の引数meanにmu_truth、標準偏差の引数sdにlambda_truthの平方根の逆数を指定します。生成したN個のデータをx_nとします。
観測データのヒストグラムを真の分布と重ねて確認します。
# 観測データを格納 data_df <- tibble::tibble(x = x_n) # 観測データのヒストグラムを作成 ggplot() + geom_histogram(data = data_df, mapping = aes(x = x, y = ..density.., fill = "data"), bins = 30) + # 観測データ:(密度) geom_line(data = model_df, mapping = aes(x = x, y = dens, color = "model"), size = 1, linetype = "dashed") + # 真の分布 scale_fill_manual(values = c(model = NA, data = "pink"), na.value = NA, labels = c(model = "true model", data = "observation data"), name = "") + # バーの色:(凡例表示用) scale_color_manual(values = c(model = "red", data = "pink"), labels = c(model = "true model", data = "observation data"), name = "") + # 線の色:(凡例表示用) guides(color = guide_legend(override.aes = list(size = c(0.5, 0.5), linetype = c("dashed", "blank")))) + # 凡例の体裁:(凡例表示用) labs(title = "Gaussian Distribution", subtitle = parse(text = paste0("list(mu==", mu_truth, ", lambda==", lambda_truth, ", N==", N, ")")), x = "x", y = "density")
観測データx_nをデータフレームに格納して、geom_histogram()でヒストグラムを描画します。y軸の引数yに..density..を指定すると、度数を密度に変換して描画します。
データ数が十分に大きいと、分布の形状が真の分布に近付きます。
事前分布の設定
次は、尤度に対する共役事前分布を設定します。ガウス分布の平均パラメータ$\mu$と精度パラメータ$\lambda$の同時事前分布としてガウス・ガンマ分布$\mathrm{NG}(\nu, \lambda | m, \beta, a, b)$を設定します。
ガンマ分布のグラフ作成については「【R】ガウス-ガンマ分布の作図 - からっぽのしょこ」を参照してください。
$\mu$の事前分布(ガウス分布)のパラメータ(超パラメータ)$m, \beta$と、$\lambda$の事前分布(ガンマ分布)のパラメータ(超パラメータ)$a, b$を設定します。
# μの事前分布のパラメータを指定 m <- 0 beta <- 1 # λの事前分布のパラメータを指定 a <- 1 b <- 1
平均パラメータ$m$と精度パラメータの係数$\beta > 0$、形状パラメータ$a > 0$と尺度パラメータ$b > 0$を指定します。
ガウス分布の確率変数がとり得る値$\mu$と、ガンマ分布の確率変数がとり得る値$\lambda$を作成します。
# グラフ用のμの値を作成 mu_vec <- seq(mu_truth - 40, mu_truth + 40, length.out = 201) # グラフ用のλの値を作成 lambda_vec <- seq(0, lambda_truth * 4, length.out = 201) head(mu_vec); head(lambda_vec)
## [1] -15.0 -14.6 -14.2 -13.8 -13.4 -13.0 ## [1] 0e+00 2e-04 4e-04 6e-04 8e-04 1e-03
この例では、$\mu$を真の値を中心に40の範囲、$\lambda$を0Pから真の値の4`倍を範囲とします。
ガウス・ガンマ分布の確率変数がとり得る値(点)$(\mu, \lambda)$を作成します。
# グラフ用のμとλの点を作成 mu_lambda_mat <- tidyr::expand_grid(mu = mu_vec, lambda = lambda_vec) |> # 格子点を作成 as.matrix() # マトリクスに変換 head(mu_lambda_mat)
## mu lambda ## [1,] -15 0e+00 ## [2,] -15 2e-04 ## [3,] -15 4e-04 ## [4,] -15 6e-04 ## [5,] -15 8e-04 ## [6,] -15 1e-03
$\mu$の値mu_vevと$\lambda$の値lambda_vecの要素の全ての組み合わせ($\mu$と$\lambda$の格子点)をexpand.grid()で作成して、マトリクスに変換します。
$\mu$と$\lambda$の事前分布を計算します。
# μとλの(同時)事前分布を計算 prior_df <- tidyr::tibble( mu = mu_lambda_mat[, 1], # 確率変数μ lambda = mu_lambda_mat[, 2], # 確率変数λ N_dens = dnorm(x = mu, mean = m, sd = 1/sqrt(beta*lambda)), # μの確率密度 Gam_dens = dgamma(x = lambda, shape = a, rate = b), # λの確率密度 density = N_dens * Gam_dens # 確率密度 ) prior_df
## # A tibble: 40,401 × 5 ## mu lambda N_dens Gam_dens density ## <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 -15 0 0 1 0 ## 2 -15 0.0002 0.00552 1.00 0.00552 ## 3 -15 0.0004 0.00763 1.00 0.00762 ## 4 -15 0.0006 0.00913 0.999 0.00913 ## 5 -15 0.0008 0.0103 0.999 0.0103 ## 6 -15 0.001 0.0113 0.999 0.0113 ## 7 -15 0.0012 0.0121 0.999 0.0121 ## 8 -15 0.0014 0.0128 0.999 0.0127 ## 9 -15 0.0016 0.0133 0.998 0.0133 ## 10 -15 0.0018 0.0138 0.998 0.0138 ## # … with 40,391 more rows
$\mu$と$\lambda$の組み合わせ(点$(\mu, \lambda)$)ごとに確率密度を計算します。
dnorm()の確率変数の引数xにmu_valsの値、平均の引数meanにm、標準偏差の引数sdにbetaとlambda_valsの値の積の平方根の逆数を指定して、ガウス分布の確率密度を計算します。
dgamma()の確率変数の引数xにlambda_valsの値、形状の引数shapeにa、尺度の引数rateにbを指定して、ガンマ分布の確率密度を計算します。
2つの確率密度の積がガウス・ガンマ分布の確率密度です。
$\mu$と$\lambda$の事前分布のグラフを作成します。
# μとλの(同時)事前分布を作図:等高線図 ggplot() + geom_contour(data = prior_df, aes(x = mu, y = lambda, z = density, color = ..level.., alpha = "prior")) + # μとλの事前分布 geom_point(mapping = aes(x = mu_truth, y = lambda_truth, alpha = "data"), color = "red", size = 6, shape = 4) + # 真のパラメータ scale_alpha_manual(values = c(param = 1, prior = 1), labels = c(param = "true parameter", prior = "prior"), name = "") + # (凡例表示用の黒魔術) guides(alpha = guide_legend(override.aes = list(shape = c(4, NA), linetype = c("blank", "solid")))) + # (凡例表示用の黒魔術) labs(title = "Gaussian-Gamma Distribution", subtitle = parse(text = paste0("list(m==", m, ", beta==", beta, ", a==", a, ", b==", b, ")")), color = "density", x = expression(mu), y = expression(lambda))
# μとλの(同時)事前分布を作図:塗りつぶし等高線図 ggplot() + geom_contour_filled(data = prior_df, aes(x = mu, y = lambda, z = density, fill = ..level..), alpha = 0.8, size = 0) + # μとλの事前分布 geom_point(mapping = aes(x = mu_truth, y = lambda_truth, color = "param"), size = 6, shape = 4) + # 真のパラメータ scale_color_manual(breaks = "param", values = "red", labels = "true parameter", name = "") + # (凡例表示用の黒魔術) labs(title = "Gaussian-Gamma Distribution", subtitle = parse(text = paste0("list(m==", m, ", beta==", beta, ", a==", a, ", b==", b, ")")), fill = "density", x = expression(mu), y = expression(lambda))
geom_contour()またはgeom_contour_filled()で確率密度の等高線図を描画します。
真のパラメータの値を赤色のバツ印で示します。
事後分布の計算
観測データ$\mathbf{X}$から平均パラメータ$\mu$と精度パラメータ$\lambda$の同時事後分布$\mathrm{NG}(\mu, \lambda | \hat{m}, \hat{\beta}, \hat{a}, \hat{b})$を求めます(パラメータ$\mu, \lambda$を分布推定します)。
観測データ$\mathbf{X}$を用いて、$\mu$の事後分布(ガウス分布)のパラメータ$\hat{m}, \hat{\beta}$を計算します。
# μの事後分布のパラメータを計算:式(3.83) beta_hat <- N + beta m_hat <- (sum(x_n) + beta * m) / beta_hat m_hat; beta_hat
## [1] 24.27453 ## [1] 51
$\mu$の事後分布のパラメータは、次の式で計算できます。
続いて、$\lambda$の事後分布(ガンマ分布)のパラメータ$\hat{a}, \hat{b}$を計算します。
# λの事後分布のパラメータを計算:式(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 a_hat; b_hat
## [1] 26 ## [1] 2555.685
$\lambda$の事後分布のパラメータは、次の式で計算できます。
$\mu$と$\lambda$の事後分布(ガウス・ガンマ分布)します。
# μとλの(同時)事後分布を計算 posterior_df <- tidyr::tibble( mu = mu_lambda_mat[, 1], # 確率変数μ lambda = mu_lambda_mat[, 2], # 確率変数λ N_dens = dnorm(x = mu, mean = m_hat, sd = 1/sqrt(beta_hat*lambda)), # μの確率密度 Gam_dens = dgamma(x = lambda, shape = a_hat, rate = b_hat), # λの確率密度 density = N_dens * Gam_dens # 確率密度 ) posterior_df
## # A tibble: 40,401 × 5 ## mu lambda N_dens Gam_dens density ## <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 -15 0 0 0 0 ## 2 -15 0.0002 1.54e- 5 5.11e-30 7.89e-35 ## 3 -15 0.0004 8.37e- 9 1.03e-22 8.61e-31 ## 4 -15 0.0006 3.93e-12 1.56e-18 6.12e-30 ## 5 -15 0.0008 1.74e-15 1.24e-15 2.16e-30 ## 6 -15 0.001 7.45e-19 1.97e-13 1.47e-31 ## 7 -15 0.0012 3.13e-22 1.13e-11 3.53e-33 ## 8 -15 0.0014 1.30e-25 3.19e-10 4.13e-35 ## 9 -15 0.0016 5.31e-29 5.39e- 9 2.86e-37 ## 10 -15 0.0018 2.16e-32 6.14e- 8 1.33e-39 ## # … with 40,391 more rows
更新した超パラメータm_hat, beta_hat, a_hat, b_hatを用いて、事前分布のときと同様にして計算します。
$\mu$と$\lambda$の事後分布のグラフを作成します。
# パラメータラベル用の文字列を作成 param_text <- paste0( "list(hat(m)==", round(m_hat, 2), ", hat(beta)==", beta_hat, ", hat(a)==", a_hat, ", hat(b)==", round(b_hat, 1), ")" ) # μとλの(同時)事後分布を作図:等高線図 ggplot() + geom_contour(data = posterior_df, aes(x = mu, y = lambda, z = density, color = ..level.., alpha = "posterior")) + # μとλの事後分布 geom_point(mapping = aes(x = mu_truth, y = lambda_truth, alpha = "param"), color = "red", size = 6, shape = 4) + # 真のパラメータ scale_alpha_manual(values = c(param = 1, posterior = 1), labels = c(param = "true parameter", posterior = "posterior"), name = "") + # (凡例表示用の黒魔術) guides(alpha = guide_legend(override.aes = list(shape = c(4, NA), linetype = c("blank", "solid")))) + # (凡例表示用の黒魔術) labs(title = "Gaussian-Gamma Distribution", subtitle = parse(text = param_text), color = "density", x = expression(mu), y = expression(lambda))
# μとλの(同時)事後分布を作図:塗りつぶし等高線図 ggplot() + geom_contour_filled(data = posterior_df, aes(x = mu, y = lambda, z = density, fill = ..level..), alpha = 0.8) + # μとλの事後分布 geom_point(mapping = aes(x = mu_truth, y = lambda_truth, color = "param"), size = 6, shape = 4) + # 真のパラメータ scale_color_manual(breaks = "param", values = "red", labels = "true parameter", name = "") + # (凡例表示用の黒魔術) labs(title = "Gaussian-Gamma Distribution", subtitle = parse(text = param_text), fill = "density", x = expression(mu), y = expression(lambda))
$\mu, \lambda$の真の値付近をピークとする分布を推定できています。
予測分布の計算
最後に、観測データ$\mathbf{X}$または(観測データから求めた)事後分布のパラメータ$\hat{a}, \hat{b}$から未観測のデータ$x_{*}$の予測分布$\mathrm{St}(x_{*} | \hat{\mu}_s, \hat{\lambda}_s, \hat{\nu}_s)$を求めます。
t分布のグラフ作成については「【R】1次元スチューデントのt分布の作図 - からっぽのしょこ」を参照してください。
$\mu, \lambda$の事後分布のパラメータ$\hat{m}, \hat{\beta}, \hat{a}, \hat{b}$、または観測データ$\mathbf{X}$と$\mu, \lambda$の事前分布のパラメータ$m, \beta, a, b$を用いて、予測分布(スチューデントのt分布)のパラメータ$\hat{\mu}_s, \hat{\lambda}_s, \hat{\nu}_s$を計算します。
# 予測分布のパラメータを計算:式(3.95') mu_st_hat <- m_hat lambda_st_hat <- beta_hat * a_hat / (1 + beta_hat) / b_hat nu_st_hat <- 2 * a_hat #mu_st_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_st_hat <- numer_lambda / denom_lambda #nu_st_hat <- N + 2 * a mu_st_hat; lambda_st_hat; nu_st_hat
## [1] 24.27453 ## [1] 0.009977756 ## [1] 52
予測分布のパラメータは、次の式で計算できます。
それぞれの1つ目の式だと、事後分布のパラメータmu_hat, beta_hat, a_hat, b_hat、2つ目の式だと、観測データx_nと事前分布のパラメータmu, beta, a, bを使って計算できます(嘘ですね$\hat{\beta}$と$\hat{m}$が残ってますねでもめんどい…)。
予測分布を計算します。
# 予測分布を計算:式(3.76) predict_df <- tibble::tibble( x = x_vec, # 確率変数 dens = LaplacesDemon::dst(x = x_vec, mu = mu_st_hat, sigma = 1/sqrt(lambda_st_hat), nu = nu_st_hat) # 確率密度 ) predict_df
## # A tibble: 501 × 2 ## x dens ## <dbl> <dbl> ## 1 -15 0.0000412 ## 2 -14.8 0.0000432 ## 3 -14.7 0.0000454 ## 4 -14.5 0.0000477 ## 5 -14.4 0.0000501 ## 6 -14.2 0.0000526 ## 7 -14.0 0.0000552 ## 8 -13.9 0.0000580 ## 9 -13.7 0.0000609 ## 10 -13.6 0.0000639 ## # … with 491 more rows
スチューデントのt分布の確率密度は、LaplacesDemonパッケージのdst()で計算できます。確率変数の引数xにx_vec、位置パラメータの引数muにmu_st_hat、スケールパラメータの引数sigmaにlambda_st_hatの平方根の逆数、形状パラメータ(自由度)の引数nuにnu_st_hatを指定します。
予測分布のグラフを真の分布と重ねて作成します。
# パラメータラベル用の文字列を作成 param_text <- paste0( "list(N==", N, ", hat(mu)[s]==", round(mu_st_hat, 2), ", hat(lambda)[s]==", round(lambda_st_hat, 5), ", hat(nu)[s]==", nu_st_hat, ")" ) # 予測分布を作図 ggplot() + geom_line(data = model_df, mapping = aes(x = x, y = dens, color = "model"), size = 1, linetype = "dashed") + # 真の分布 geom_line(data = predict_df, mapping = aes(x = x, y = dens, color = "predict"), size = 1) + # 予測分布 scale_color_manual(values = c(model = "red", predict = "purple"), labels = c(model = "true model", predict = "predict"), name = "") + # 線の色:(凡例表示用) guides(color = guide_legend(override.aes = list(size = c(0.5, 0.5), linetype = c("dashed", "solid")))) + # 凡例の体裁:(凡例表示用) labs(title = "Student's t Distribution", subtitle = parse(text = param_text), x = "x", y = "density")
観測データが増えると、予測分布が真の分布に近付きます。
以上で、ガウス分布のベイズ推論を実装できました。
学習推移の可視化
gganimateパッケージを利用して、パラメータの推定値(更新値)の推移(分布の変化)をアニメーション(gif画像)で確認します。
モデルの設定
真の分布のパラメータ$\mu, \lambda$と、$\mu$と$\lambda$の事前分布のパラメータ(の初期値)$m, \beta, a, b$を設定します。
# 真の平均パラメータを指定 mu_truth <- 25 # 真の精度パラメータを指定 lambda_truth <- 0.01 # μの事前分布のパラメータを指定 m <- 0 beta <- 1 # λの事前分布のパラメータを指定 a <- 1 b <- 1 # データ数(試行回数)を指定 N <- 100 # グラフ用のμの値を作成 mu_vec <- seq(mu_truth - 30, mu_truth + 30, length.out = 201) # グラフ用のλの値を作成 lambda_vec <- seq(0, lambda_truth * 4, length.out = 201) # グラフ用のμとλの点を作成 mu_lambda_mat <- tidyr::expand_grid(mu = mu_vec, lambda = lambda_vec) |> # 格子点を作成 as.matrix() # マトリクスに変換 # グラフ用のxの値を作成 x_vec <- seq( mu_truth - 1/sqrt(lambda_truth) * 4, mu_truth + 1/sqrt(lambda_truth) * 4, length.out = 251 )
実装時と同じく、パラメータを指定して、確率変数がとり得る値(x軸とy軸の値)を作成します。グラフが粗い場合や処理が重い場合は、mu_vec, lambda_vec, x_vecの要素数(by引数の値)または値の間隔(length.out引数の値)を調整してください。
事後分布と予測分布のパラメータの計算(更新)処理について、2通りの方法を載せます。目的に応じて使い分けてください。
推論処理：for関数による処理
1つ目の処理方法では、for()を使って、1データずつ生成してパラメータの更新を繰り返し処理します。こちらの方が、前ステップで求めた事後分布(のパラメータ)を次ステップの事前分布(のパラメータ)として用いる逐次学習をイメージしやすいです。
・コード(クリックで展開)
超パラメータの初期値を使って、事前分布を計算します。
# μとλの事後分布(ガウス・ガンマ分布)を計算 anime_posterior_df <- tidyr::tibble( mu = mu_lambda_mat[, 1], # 確率変数μ lambda = mu_lambda_mat[, 2], # 確率変数λ N_dens = dnorm(x = mu, mean = m, sd = 1/sqrt(beta*lambda)), # μの確率密度 Gam_dens = dgamma(x = lambda, shape = a, rate = b), # λの確率密度 density = N_dens * Gam_dens, # 確率密度 param = paste0("N=", 0, ", m=", m, ", beta=", beta, ", a=", a, ", b=", b) |> as.factor() # フレーム切替用ラベル ) anime_posterior_df
## # A tibble: 40,401 × 6 ## mu lambda N_dens Gam_dens density param ## <dbl> <dbl> <dbl> <dbl> <dbl> <fct> ## 1 -5 0 0 1 0 N=0, m=0, beta=1, a=1, b=1 ## 2 -5 0.0002 0.00563 1.00 0.00563 N=0, m=0, beta=1, a=1, b=1 ## 3 -5 0.0004 0.00794 1.00 0.00794 N=0, m=0, beta=1, a=1, b=1 ## 4 -5 0.0006 0.00970 0.999 0.00969 N=0, m=0, beta=1, a=1, b=1 ## 5 -5 0.0008 0.0112 0.999 0.0112 N=0, m=0, beta=1, a=1, b=1 ## 6 -5 0.001 0.0125 0.999 0.0124 N=0, m=0, beta=1, a=1, b=1 ## 7 -5 0.0012 0.0136 0.999 0.0136 N=0, m=0, beta=1, a=1, b=1 ## 8 -5 0.0014 0.0147 0.999 0.0146 N=0, m=0, beta=1, a=1, b=1 ## 9 -5 0.0016 0.0156 0.998 0.0156 N=0, m=0, beta=1, a=1, b=1 ## 10 -5 0.0018 0.0165 0.998 0.0165 N=0, m=0, beta=1, a=1, b=1 ## # … with 40,391 more rows
現在のパラメータを文字列結合し因子型に変換して、フレーム切替用のラベル列とします。
また、事前分布のパラメータを使って、予測分布のパラメータと予測分布を計算します。
# 初期値による予測分布のパラメータを計算:式(3.95) mu_st <- m lambda_st <- beta * a / (1 + beta) / b nu_st <- 2 * a # 初期値による予測分布(スチューデントのt分布)を計算:式(3.76) anime_predict_df <- tibble::tibble( x = x_vec, # 確率変数 dens = LaplacesDemon::dst(x = x_vec, mu = mu_st, sigma = 1/sqrt(lambda_st), nu = nu_st), # 確率密度 param = paste0("N=", 0, ", mu_s=",mu_st, ", lambda_s=", lambda_st, ", nu_s=", nu_st) |> as.factor() # フレーム切替用ラベル ) anime_predict_df
## # A tibble: 251 × 3 ## x dens param ## <dbl> <dbl> <fct> ## 1 -15 0.000577 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 2 -14.7 0.000615 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 3 -14.4 0.000656 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 4 -14.0 0.000701 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 5 -13.7 0.000750 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 6 -13.4 0.000804 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 7 -13.1 0.000863 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 8 -12.8 0.000928 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 9 -12.4 0.00100 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 10 -12.1 0.00108 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## # … with 241 more rows
同様に、フレーム切替用のラベル列を作成します。
パラメータの更新処理をN回繰り返します。
# 観測データの受け皿を作成 x_n <- rep(NA, times = N) # ベイズ推論 for(n in 1:N){ # ガウス分布に従うデータを生成 x_n[n] <- rnorm(n = 1, mean = mu_truth, sd = 1/sqrt(lambda_truth)) # μの事後分布のパラメータを更新:式(3.83) beta_old <- beta m_old <- m beta <- 1 + beta m <- (x_n[n] + beta_old * m) / beta # λの事後分布のパラメータを更新:式(3.88) a <- 0.5 + a b <- 0.5 * (x_n[n]^2 + beta_old * m_old^2 - beta * m^2) + b # μとλの事後分布(ガウス・ガンマ分布)を計算 tmp_posterior_df <- tidyr::tibble( mu = mu_lambda_mat[, 1], # 確率変数μ lambda = mu_lambda_mat[, 2], # 確率変数λ N_dens = dnorm(x = mu, mean = m, sd = 1/sqrt(beta*lambda)), # μの確率密度 Gam_dens = dgamma(x = lambda, shape = a, rate = b), # λの確率密度 density = N_dens * Gam_dens, # 確率密度 param = paste0("N=", n, ", m=", round(m, 2), ", beta=", beta, ", a=", a, ", b=", round(b, 1)) |> as.factor() # フレーム切替用ラベル ) # 予測分布のパラメータを更新:式(3.95) mu_st <- m lambda_st <- beta * a / (1 + beta) / b nu_st <- 2 * a # 予測分布(スチューデントのt分布)を計算:式(3.76) tmp_predict_df <- tibble::tibble( x = x_vec, # 確率変数 dens = LaplacesDemon::dst(x = x_vec, mu = mu_st, sigma = 1/sqrt(lambda_st), nu = nu_st), # 確率密度 param = paste0( "N=", n, ", mu_s=", round(mu_st, 2), ", lambda_s=", round(lambda_st, 5), ", nu_s=", nu_st ) |> as.factor() # フレーム切替用ラベル ) # 推論結果を結合 anime_posterior_df <- rbind(anime_posterior_df, tmp_posterior_df) anime_predict_df <- rbind(anime_predict_df, tmp_predict_df) }
超パラメータに関して、$\hat{m}, \hat{\beta}$に対応するm_hat, beta_hatや$\hat{a}, \hat{b}$に対応するa_hat, b_hatを新たに作るのではなく、m, betaと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が既に各パラメータに加えられています。
更新したパラメータを使って事後分布と予測分布を計算して、それぞれ計算結果をanime_***_dfに結合していきます。
結果を確認します。
# 確認 head(x_n); anime_posterior_df; anime_predict_df
## [1] 30.5292221 36.0197926 37.6726401 -0.2926592 18.7791670 18.7824555 ## # A tibble: 4,080,501 × 6 ## mu lambda N_dens Gam_dens density param ## <dbl> <dbl> <dbl> <dbl> <dbl> <fct> ## 1 -5 0 0 1 0 N=0, m=0, beta=1, a=1, b=1 ## 2 -5 0.0002 0.00563 1.00 0.00563 N=0, m=0, beta=1, a=1, b=1 ## 3 -5 0.0004 0.00794 1.00 0.00794 N=0, m=0, beta=1, a=1, b=1 ## 4 -5 0.0006 0.00970 0.999 0.00969 N=0, m=0, beta=1, a=1, b=1 ## 5 -5 0.0008 0.0112 0.999 0.0112 N=0, m=0, beta=1, a=1, b=1 ## 6 -5 0.001 0.0125 0.999 0.0124 N=0, m=0, beta=1, a=1, b=1 ## 7 -5 0.0012 0.0136 0.999 0.0136 N=0, m=0, beta=1, a=1, b=1 ## 8 -5 0.0014 0.0147 0.999 0.0146 N=0, m=0, beta=1, a=1, b=1 ## 9 -5 0.0016 0.0156 0.998 0.0156 N=0, m=0, beta=1, a=1, b=1 ## 10 -5 0.0018 0.0165 0.998 0.0165 N=0, m=0, beta=1, a=1, b=1 ## # … with 4,080,491 more rows ## # A tibble: 25,351 × 3 ## x dens param ## <dbl> <dbl> <fct> ## 1 -15 0.000577 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 2 -14.7 0.000615 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 3 -14.4 0.000656 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 4 -14.0 0.000701 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 5 -13.7 0.000750 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 6 -13.4 0.000804 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 7 -13.1 0.000863 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 8 -12.8 0.000928 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 9 -12.4 0.00100 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 10 -12.1 0.00108 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## # … with 25,341 more rows
それぞれ「lambda_vec, x_vecの要素数」掛ける「N+1」行のデータフレームになります。行数が増えすぎないように注意してください。
推論処理：tidyverseパッケージによる処理
2つ目の処理方法では、tidyverseパッケージの関数を使って、一度の処理でN+1回分の計算をします。こちらの方が、処理時間が短いです。
・コード(クリックで展開)
ガウス分布に従う$N$個のデータを生成します。
# ガウス分布に従うデータを生成 x_n <- rnorm(n = N, mean = mu_truth, sd = 1/sqrt(lambda_truth)) head(x_n)
## [1] 30.5292221 36.0197926 37.6726401 -0.2926592 18.7791670 18.7824555
事前分布を含めたN+1回分の事後分布と予測分布のパラメータを計算します。
# 試行ごとに事後分布と予測分布のパラメータを計算 anime_param_df <- tibble::tibble( # 計算用の値 n = 0:N, sum_x = cumsum(c(0, x_n)), sum_x2 = cumsum(c(0, x_n)^2), # μの事後分布のパラメータ:式(3.83) beta_hat = n + beta, m_hat = (sum_x + beta*m) / (n + beta), # λの事後分布のパラメータ:式(3.88) a_hat = 0.5*n + a, b_hat = 0.5 * (sum_x2 + beta*m^2 - beta_hat*m_hat^2) + b, # 予測分布のパラメータ:式(3.95) mu_st_hat = m_hat, lambda_st_hat = beta_hat * a_hat / (1 + beta_hat) / b_hat, nu_st_hat = 2 * a_hat ) |> # パラメータを計算 dplyr::select(!c(sum_x, sum_x2)) # 不要な列を削除 anime_param_df
## # A tibble: 101 × 8 ## n beta_hat m_hat a_hat b_hat mu_st_hat lambda_st_hat nu_st_hat ## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 0 1 0 1 1 0 0.5 2 ## 2 1 2 15.3 1.5 234. 15.3 0.00427 3 ## 3 2 3 22.2 2 378. 22.2 0.00397 4 ## 4 3 4 26.1 2.5 468. 26.1 0.00428 5 ## 5 4 5 20.8 3 745. 20.8 0.00335 6 ## 6 5 6 20.5 3.5 747. 20.5 0.00402 7 ## 7 6 7 20.2 4 748. 20.2 0.00468 8 ## 8 7 8 20.8 4.5 757. 20.8 0.00528 9 ## 9 8 9 21.7 5 790. 21.7 0.00569 10 ## 10 9 10 22.7 5.5 834. 22.7 0.00600 11 ## # … with 91 more rows
初期値に対応する0を含む観測データ数(試行回数)の列n(0からNの整数)を作成します。
観測データの累積和$\sum_{i=0}^n x_i^2$の列sum_xと2乗の累積和$\sum_{i=0}^n x_i$の列sum_x2を作成します。累積和はcumsum()で計算できます。
この3つの列と初期値のオブジェクトを使って、事後分布と予測分布のパラメータの計算します。
事前分布を含めたN+1回分の事後分布を計算します。
# 試行ごとにμとλの(同時)事後分布を計算 anime_posterior_df <- tidyr::expand_grid( anime_param_df |> dplyr::select(n, m = m_hat, beta = beta_hat, a = a_hat, b = b_hat), # 利用するパラメータを取得 mu = mu_vec, # 確率変数μ lambda = lambda_vec # 確率変数λ ) |> # パラメータごとに確率変数(格子点)を複製 dplyr::mutate( N_dens = dnorm(x = mu, mean = m, sd = 1/sqrt(beta*lambda)), # μの確率密度 Gam_dens = dgamma(x = lambda, shape = a, rate = b), # λの確率密度 density = N_dens * Gam_dens, # 確率密度 param = paste0("N=", n, ", m=", round(m, 2), ", beta=", beta, ", a=", a, ", b=", round(b, 1)) |> (\(.){factor(., levels = unique(.))})() # フレーム切替用ラベル ) anime_posterior_df
## # A tibble: 4,080,501 × 11 ## n m beta a b mu lambda N_dens Gam_dens density param ## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fct> ## 1 0 0 1 1 1 -5 0 0 1 0 N=0, m=0… ## 2 0 0 1 1 1 -5 0.0002 0.00563 1.00 0.00563 N=0, m=0… ## 3 0 0 1 1 1 -5 0.0004 0.00794 1.00 0.00794 N=0, m=0… ## 4 0 0 1 1 1 -5 0.0006 0.00970 0.999 0.00969 N=0, m=0… ## 5 0 0 1 1 1 -5 0.0008 0.0112 0.999 0.0112 N=0, m=0… ## 6 0 0 1 1 1 -5 0.001 0.0125 0.999 0.0124 N=0, m=0… ## 7 0 0 1 1 1 -5 0.0012 0.0136 0.999 0.0136 N=0, m=0… ## 8 0 0 1 1 1 -5 0.0014 0.0147 0.999 0.0146 N=0, m=0… ## 9 0 0 1 1 1 -5 0.0016 0.0156 0.998 0.0156 N=0, m=0… ## 10 0 0 1 1 1 -5 0.0018 0.0165 0.998 0.0165 N=0, m=0… ## # … with 4,080,491 more rows
anime_param_dfから利用する列を取り出して、x軸の値mu_vecとy軸の値lambda_vecとの全ての組み合わせをexpand_grid()で作成します。これにより、mu_vec, lambda_vecの格子点をN+1フレーム分に複製できます。
確率変数(格子点)とパラメータの組み合わせごとに確率密度を計算して、パラメータごとにフレーム切替用のラベルを作成します。
因子型への変換処理では、無名関数function()の省略記法\()を使って、(\(引数){引数を使った具体的な処理})()としています。直前のパイプ演算子を%>%にすると、行全体(\(引数){処理})()を{}の中の処理(この例だとfactor(., levels = unique(.)))に置き換えられます(置き換えられるように引数名を.にしています)。
同様に、初期値を含めたN+1回分の予測分布を計算します。
# 試行ごとに予測分布(スチューデントのt分布)を計算 anime_predict_df <- tidyr::expand_grid( anime_param_df |> dplyr::select(n, mu_st = mu_st_hat, lambda_st = lambda_st_hat, nu_st = nu_st_hat), # 利用するパラメータを取得 x = x_vec # 確率変数 ) |> # パラメータごとに確率変数を複製 dplyr::mutate( dens = LaplacesDemon::dst(x = x_vec, mu = mu_st, sigma = 1/sqrt(lambda_st), nu = nu_st), # 確率密度 param = paste0("N=", n, ", mu_s=", round(mu_st, 2), ", lambda_s=", round(lambda_st, 5), ", nu_s=", nu_st) |> (\(.){factor(., levels = unique(.))})() # フレーム切替用ラベル ) # 事後分布を計算:式(3.76) anime_predict_df
## # A tibble: 25,351 × 7 ## n mu_st lambda_st nu_st x dens param ## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <fct> ## 1 0 0 0.5 2 -15 0.000577 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 2 0 0 0.5 2 -14.7 0.000615 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 3 0 0 0.5 2 -14.4 0.000656 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 4 0 0 0.5 2 -14.0 0.000701 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 5 0 0 0.5 2 -13.7 0.000750 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 6 0 0 0.5 2 -13.4 0.000804 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 7 0 0 0.5 2 -13.1 0.000863 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 8 0 0 0.5 2 -12.8 0.000928 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 9 0 0 0.5 2 -12.4 0.00100 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 10 0 0 0.5 2 -12.1 0.00108 N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## # … with 25,341 more rows
anime_param_dfから利用する列を取り出して、x軸の値x_vecとの全ての組み合わせをexpand_grid()で作成します。これにより、x_vecの各要素をN+1フレーム分に複製できます。
確率変数とパラメータの組み合わせごとに確率を計算して、パラメータごとにフレーム切替用のラベルを作成します。
作図処理
事後分布と予測分布の推移をそれぞれアニメーションで可視化します。
・コード(クリックで展開)
観測データと対応するラベルをデータフレームに格納します。
# 観測データを格納 anime_data_df <- tibble::tibble( x = c(NA, x_n), scaled_x = c(NA, 1/(x_n - mu_truth)^2), # 偏差の2乗の逆数に変換 param = unique(anime_posterior_df[["param"]]) # フレーム切替用ラベル ) anime_data_df
## # A tibble: 101 × 3 ## x scaled_x param ## <dbl> <dbl> <fct> ## 1 NA NA N=0, m=0, beta=1, a=1, b=1 ## 2 30.5 0.0327 N=1, m=15.26, beta=2, a=1.5, b=234 ## 3 36.0 0.00823 N=2, m=22.18, beta=3, a=2, b=377.6 ## 4 37.7 0.00623 N=3, m=26.06, beta=4, a=2.5, b=467.6 ## 5 -0.293 0.00156 N=4, m=20.79, beta=5, a=3, b=745.3 ## 6 18.8 0.0258 N=5, m=20.45, beta=6, a=3.5, b=746.9 ## 7 18.8 0.0259 N=6, m=20.21, beta=7, a=4, b=748.1 ## 8 24.8 24.5 N=7, m=20.79, beta=8, a=4.5, b=757.3 ## 9 29.4 0.0514 N=8, m=21.74, beta=9, a=5, b=790.4 ## 10 31.6 0.0231 N=9, m=22.73, beta=10, a=5.5, b=834 ## # … with 91 more rows
事前分布には観測データが影響しないので欠損値NAとして、anime_posterior_dfのラベル列と合わせて格納します。観測データと精度パラメータを対応させるため、各データから平均$\mu$を引き、2乗の逆数にして格納します。
参考として、各試行における観測データを描画するのに使います。
事後分布の推移のアニメーションを作成します。
# μとλの(同時)事後分布のアニメーションを作図:等高線図 posterior_graph <- ggplot() + geom_contour(data = anime_posterior_df, mapping = aes(x = mu, y = lambda, z = density, color = ..level.., alpha = "posterior")) + # μとλの事後分布 geom_point(mapping = aes(x = mu_truth, y = lambda_truth, alpha = "param"), color = "red", size = 6, shape = 4) + # 真のパラメータ geom_point(data = anime_data_df, mapping = aes(x = x, y = scaled_x, alpha = "data"), color = "pink", size = 6) + # 観測データ gganimate::transition_manual(param) + # フレーム scale_alpha_manual(breaks = c("posterior", "param", "data"), values = c(1, 1, 1), labels = c("posterior", "true parameter", "observation data"), name = "") + # (凡例表示用の黒魔術) guides(alpha = guide_legend(override.aes = list(alpha = c(1, 1, 1), shape = c(NA, 4, 16), linetype = c("solid", "blank", "blank")))) + # 凡例の体裁:(凡例表示用) coord_cartesian(xlim = c(min(mu_vec), max(mu_vec)), ylim = c(0, max(lambda_vec))) + labs(title = "Gaussian-Gamma Distribution", subtitle = "{current_frame}", color = "density", x = expression(mu), y = expression(lambda)) # ラベル # gif画像を作図 gganimate::animate(posterior_graph, nframes = N+1+10, end_pause=10, fps = 10, width = 800, height = 600)
# μとλの(同時)事後分布のアニメーションを作図:塗りつぶし等高線図 posterior_graph <- ggplot() + geom_contour_filled(data = anime_posterior_df, mapping = aes(x = mu, y = lambda, z = density, fill = ..level..), alpha = 0.8) + # μとλの事後分布 geom_point(mapping = aes(x = mu_truth, y = lambda_truth, color = "param"), size = 6, shape = 4) + # 真のパラメータ geom_point(data = anime_data_df, mapping = aes(x = x, y = scaled_x, color = "data"), size = 6) + # 観測データ gganimate::transition_manual(param) + # フレーム scale_color_manual(breaks = c("param", "data"), values = c("red", "pink"), labels = c("true parameter", "observation data"), name = "") + # (凡例表示用の黒魔術) guides(color = guide_legend(override.aes = list(shape = c(4, 16)))) + # 凡例の体裁:(凡例表示用) coord_cartesian(xlim = c(min(mu_vec), max(mu_vec)), ylim = c(0, max(lambda_vec))) + labs(title = "Gaussian-Gamma Distribution", subtitle = "{current_frame}", fill = "density", x = expression(mu), y = expression(lambda)) # ラベル # gif画像を作図 gganimate::animate(posterior_graph, nframes = N+1+10, end_pause=10, fps = 10, width = 800, height = 600)
フレームの順番を示す列をtransition_manual()に指定して、animate()でgif画像を作成します。事前分布(初期値)を含むため、フレーム数はN+1です。
観測データと対応するラベルをデータフレームに格納します。
# 観測データを格納 anime_data_df <- tibble::tibble( x = c(NA, x_n), param = unique(anime_predict_df[["param"]]) # フレーム切替用ラベル ) anime_data_df
## # A tibble: 101 × 2 ## x param ## <dbl> <fct> ## 1 NA N=0, mu_s=0, lambda_s=0.5, nu_s=2 ## 2 30.5 N=1, mu_s=15.26, lambda_s=0.00427, nu_s=3 ## 3 36.0 N=2, mu_s=22.18, lambda_s=0.00397, nu_s=4 ## 4 37.7 N=3, mu_s=26.06, lambda_s=0.00428, nu_s=5 ## 5 -0.293 N=4, mu_s=20.79, lambda_s=0.00335, nu_s=6 ## 6 18.8 N=5, mu_s=20.45, lambda_s=0.00402, nu_s=7 ## 7 18.8 N=6, mu_s=20.21, lambda_s=0.00468, nu_s=8 ## 8 24.8 N=7, mu_s=20.79, lambda_s=0.00528, nu_s=9 ## 9 29.4 N=8, mu_s=21.74, lambda_s=0.00569, nu_s=10 ## 10 31.6 N=9, mu_s=22.73, lambda_s=0.006, nu_s=11 ## # … with 91 more rows
こちらは、anime_predict_dfのラベル列を使います。
また、各試行までの観測データをデータフレームに格納します。
# 観測データを複製して格納 anime_alldata_df <- tidyr::expand_grid( frame = 1:N, # フレーム番号 n = 1:N # 試行回数 ) |> # 全ての組み合わせを作成 dplyr::filter(n < frame) |> # フレーム番号以前のデータを抽出 dplyr::mutate( x = x_n[n], param = unique(anime_predict_df[["param"]])[frame+1] ) # 対応するデータとラベルを抽出 anime_alldata_df
## # A tibble: 4,950 × 4 ## frame n x param ## <int> <int> <dbl> <fct> ## 1 2 1 30.5 N=2, mu_s=22.18, lambda_s=0.00397, nu_s=4 ## 2 3 1 30.5 N=3, mu_s=26.06, lambda_s=0.00428, nu_s=5 ## 3 3 2 36.0 N=3, mu_s=26.06, lambda_s=0.00428, nu_s=5 ## 4 4 1 30.5 N=4, mu_s=20.79, lambda_s=0.00335, nu_s=6 ## 5 4 2 36.0 N=4, mu_s=20.79, lambda_s=0.00335, nu_s=6 ## 6 4 3 37.7 N=4, mu_s=20.79, lambda_s=0.00335, nu_s=6 ## 7 5 1 30.5 N=5, mu_s=20.45, lambda_s=0.00402, nu_s=7 ## 8 5 2 36.0 N=5, mu_s=20.45, lambda_s=0.00402, nu_s=7 ## 9 5 3 37.7 N=5, mu_s=20.45, lambda_s=0.00402, nu_s=7 ## 10 5 4 -0.293 N=5, mu_s=20.45, lambda_s=0.00402, nu_s=7 ## # … with 4,940 more rows
フレーム番号と試行回数(データ番号)の全ての組み合わせを作成して、フレーム番号未満のデータを抽出します。フレーム切替用のラベルは、フレーム番号を使って抽出します。
観測データの無い初回(0回目)と1回目は描画しないので、フレーム番号が2からになります。
予測分布の推移のアニメーションを作成します。
# 真の分布を計算:式(2.64) model_df <- tibble::tibble( x = x_vec, # 確率変数 dens = dnorm(x = x_vec, mean = mu_truth, sd = 1/sqrt(lambda_truth)) # 確率密度 ) # 予測分布のアニメーションを作図 predict_graph <- ggplot() + geom_line(data = model_df, mapping = aes(x = x, y = dens, color = "model"), size = 1, linetype = "dashed") + # 真の分布 geom_line(data = anime_predict_df, mapping = aes(x = x, y = dens, color = "predict"), size = 1) + # 予測分布 geom_point(data = anime_alldata_df, mapping = aes(x = x, y = 0), color = "pink", alpha = 0.5, size = 3) + # 過去の観測データ geom_point(data = anime_data_df, mapping = aes(x = x, y = 0, color = "data"), size = 6) + # 観測データ gganimate::transition_manual(param) + # フレーム scale_color_manual(breaks = c("model", "predict", "data"), values = c("red", "purple", "pink"), labels = c("true model", "predict", "observation data"), name = "") + # 線の色:(凡例表示用) guides(color = guide_legend(override.aes = list(size = c(0.5, 0.5, 3), linetype = c("dashed", "solid", "blank"), shape = c(NA, NA, 19)))) + # 凡例の体裁:(凡例表示用) coord_cartesian(ylim = c(0, max(model_df[["dens"]])*2)) + # 軸の表示範囲 labs(title = "Student's t Distribution", subtitle = "{current_frame}", x = "x", y = "density") # gif画像を作図 gganimate::animate(predict_graph, nframes = N+1+10, end_pause=10, fps = 10, width = 800, height = 600)
データが増えるにつれて、真のパラメータ付近の確率密度が大きくなっていくのを確認できます。
予測分布が真の分布に近付いていくのを確認できます。
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
- 2021.03.18：加筆修正の際に、数式読解編とRで実装編に記事を分割しました。
修正作業って大変だけど、勘違いしてる1年前の自分に教えてあげてるんだって思うと優しい気持ちになれる(？)
- 2022.09.15：加筆修正しました。
ガウス・ガンマ分布のグラフを描けるようになりました。
【次節の内容】
