はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」によって、「数式」と「プログラム」から理解するのが目標です。
省略してある内容等ありますので、本とあわせて読んでください。
この記事では、「尤度関数を1次元ガウス分布」、「事前分布を1次元ガウス分布」とした場合の「パラメータの事後分布」と「未観測値の予測分布」の計算をR言語で実装します。
【数式読解編】
【前の節の内容】
【他の節の内容】
【この節の内容】
3.3.2 1次元ガウス分布の学習と予測：精度が未知の場合
尤度関数を1次元ガウス分布、事前分布をガンマ分布とするモデルに対するベイズ推論を実装します。人工的に生成したデータを用いて、ガウス分布の平均パラメータを推定し、また未観測データに対する予測分布を求めます。
ガウス分布については「1次元ガウス分布の定義式 - からっぽのしょこ」、t分布については「1次元スチューデントのt分布の定義式 - からっぽのしょこ」を参照してください。ガンマ分布についてはいつか書きます。
利用するパッケージを読み込みます。
# 利用パッケージ library(tidyverse) library(gganimate) library(LaplacesDemon)
この記事では、基本的にパッケージ名::関数名()の記法を使うので、パッケージを読み込む必要はありません。ただし、作図コードがごちゃごちゃしないようにパッケージ名を省略しているため、ggplot2は読み込む必要があります。
magrittrパッケージのパイプ演算子%>%ではなく、ネイティブパイプ(ベースパイプ)演算子|>を使っています。%>%に置き換えても処理できます。
学習の推移をアニメーション(gif画像)で確認するのにgganimateパッケージを利用します。不要であれば省略してください。
スチューデントのt分布の計算にLaplacesDemonパッケージを利用します。
ベイズ推論の実装
まずは、モデルを設定して、データを生成します。生成したデータを用いて、事後分布のパラメータを計算します。さらに、事後分布のパラメータを用いて、予測分布のパラメータを計算します。
生成分布の設定
データ生成分布(真の分布)として、ガウス分布$\mathcal{N}(x | \mu, \lambda^{-1})$を設定します。
ガウス分布のグラフ作成については「【R】1次元ガウス分布の作図 - からっぽのしょこ」を参照してください。
真の分布(ガウス分布)のパラメータ$\mu, \lambda$を設定します。
# 既知の平均パラメータを指定 mu <- 25 # 真の精度パラメータを指定 lambda_truth <- 0.01 sqrt(1 / lambda_truth) # 標準偏差
## [1] 10
平均パラメータ(実数)$\mu$をmuとして値を指定します。muは与えられている値として扱います。
精度パラメータ(正の実数)$\lambda$をlambda_truthとして値を指定します。精度は分散の逆数なので、値が大きいほど散らばり具合が小さくなり、逆数の平方根が標準偏差$\sigma = \frac{1}{\sqrt{\lambda}}$です。lambda_truthが真のパラメータであり、この値を求めるのがここでの目的です。
真の分布を計算して、作図用のデータフレームを作成します。
# グラフ用のxの値を作成 x_vec <- seq( mu - 1/sqrt(lambda_truth) * 4, mu + 1/sqrt(lambda_truth) * 4, length.out = 501 ) # 真の分布を計算:式(2.64) model_df <- tibble::tibble( x = x_vec, # 確率変数 dens = dnorm(x = x_vec, mean = mu, 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、標準偏差の引数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, ", 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, sd = 1/sqrt(lambda_truth)) head(x_n)
## [1] 19.37269 14.16877 18.28887 10.79315 11.11893 38.08840
生成するデータ数$N$を指定します。
ガウス分布に従う乱数は、rnorm()で生成できます。試行回数の引数nにN、平均の引数meanにmu、標準偏差の引数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, ", lambda==", lambda_truth, ", N==", N, ")")), x = "x", y = "density")
観測データx_nをデータフレームに格納して、geom_histogram()でヒストグラムを描画します。y軸の引数yに..density..を指定すると、度数を密度に変換して描画します。
データ数が十分に大きいと、分布の形状が真の分布に近付きます。
事前分布の設定
次は、尤度に対する共役事前分布を設定します。ガウス分布の精度パラメータ$\lambda$の事前分布としてガンマ分布$\mathrm{Gam}(\lambda | a, b)$を設定します。
ガンマ分布のグラフ作成については「【R】ガンマ分布の作図 - からっぽのしょこ」を参照してください。
$\lambda$の事前分布(ガンマ分布)のパラメータ(超パラメータ)$a, b$を設定します。
# λの事前分布のパラメータを指定 a <- 1 b <- 1
形状パラメータ$a > 0$と尺度パラメータ$b > 0$を指定します。
$\lambda$の事前分布を計算します。
# グラフ用のλの値を作成 lambda_vec <- seq(0, lambda_truth * 4, length.out = 501) # λの事前分布を計算:式(2.56) prior_df <- tibble::tibble( lambda = lambda_vec, # 確率変数 dens = dgamma(x = lambda_vec, shape = a, rate = b) # 確率密度 ) prior_df
## # A tibble: 501 × 2 ## lambda dens ## <dbl> <dbl> ## 1 0 1 ## 2 0.00008 1.00 ## 3 0.00016 1.00 ## 4 0.00024 1.00 ## 5 0.00032 1.00 ## 6 0.0004 1.00 ## 7 0.00048 1.00 ## 8 0.00056 0.999 ## 9 0.00064 0.999 ## 10 0.00072 0.999 ## # … with 491 more rows
ガンマ分布の確率変数がとり得る値$\lambda$を作成して、値ごとに確率密度を計算します。この例では、0から真の値の4倍までを範囲とします。
ガンマ分布の確率密度は、dgamma()で計算します。確率変数の引数xにlambda_vec、形状の引数shapeにa、尺度の引数rateにbを指定します。
$\lambda$の事前分布のグラフを作成します。
# λの事前分布を作図 ggplot() + geom_vline(mapping = aes(xintercept = lambda_truth, color = "param"), size = 1, linetype = "dashed", show.legend = FALSE) + # 真のパラメータ geom_line(data = prior_df, mapping = aes(x = lambda, y = dens, color = "prior"), size = 1) + # λの事前分布 scale_color_manual(values = c(param = "red", prior = "purple"), labels = c(param = "true parameter", prior = "prior"), name = "") + # 線の色:(凡例表示用) guides(color = guide_legend(override.aes = list(size = c(0.5, 0.5), linetype = c("dashed", "solid")))) + # 凡例の体裁:(凡例表示用) labs(title = "Gamma Distribution", subtitle = paste0("a=", a, ", b=", b), x = expression(lambda), y = "density")
真のパラメータの値を赤色の破線で示します。
事後分布の計算
観測データ$\mathbf{X}$から精度パラメータ$\lambda$の事後分布$\mathrm{Gam}(\lambda | \hat{a}, \hat{b})$を求めます(精度パラメータ$\lambda$を分布推定します)。
観測データ$\mathbf{X}$を用いて、$\lambda$の事後分布(ガンマ分布)のパラメータ$\hat{a}, \hat{b}$を計算します。
# λの事後分布のパラメータを計算:式(3.69) a_hat <- 0.5 * N + a b_hat <- 0.5 * sum((x_n - mu)^2) + b a_hat; b_hat
## [1] 26 ## [1] 2394.563
$\lambda$の事後分布のパラメータは、次の式で計算できます。
$\lambda$の事後分布を計算します。
# λの事後分布を計算:式(2.56) posterior_df <-tibble:: tibble( lambda = lambda_vec, # 確率変数 dens = dgamma(x = lambda_vec, shape = a_hat, rate = b_hat) # 確率密度 ) posterior_df
## # A tibble: 501 × 2 ## lambda dens ## <dbl> <dbl> ## 1 0 0 ## 2 0.00008 1.46e-40 ## 3 0.00016 4.04e-33 ## 4 0.00024 8.41e-29 ## 5 0.00032 9.23e-26 ## 6 0.0004 2.02e-23 ## 7 0.00048 1.59e-21 ## 8 0.00056 6.19e-20 ## 9 0.00064 1.44e-18 ## 10 0.00072 2.26e-17 ## # … with 491 more rows
更新した超パラメータa_hat, b_hatを用いて、事前分布のときと同様にして計算します。
$\lambda$の事後分布のグラフを作成します。
# λの事後分布を作図 ggplot() + geom_vline(aes(xintercept = lambda_truth, color = "param"), size = 1, linetype = "dashed", show.legend = FALSE) + # 真のパラメータ geom_line(data = posterior_df, mapping = aes(x = lambda, y = dens, color = "posterior"), size = 1) + # λの事後分布 scale_color_manual(values = c(param = "red", posterior = "purple"), labels = c(param = "true parameter", posterior = "posterior"), name = "") + # 線の色:(凡例表示用) guides(color = guide_legend(override.aes = list(size = c(0.5, 0.5), linetype = c("dashed", "solid")))) + # 凡例の体裁:(凡例表示用) labs(title = "Gamma Distribution", subtitle = parse( text = paste0("list(N==", N, ", hat(a)==", a_hat, ", hat(b)==", round(b_hat, 1), ")") ), x = expression(lambda), y = "density")
$\lambda$の真の値付近をピークとする分布を推定できています。
予測分布の計算
最後に、観測データ$\mathbf{X}$または(観測データから求めた)事後分布のパラメータ$\hat{a}, \hat{b}$から未観測のデータ$x_{*}$の予測分布$\mathrm{St}(x_{*} | \mu_s, \hat{\lambda}_s, \hat{\nu}_s)$を求めます。
t分布のグラフ作成については「【R】1次元スチューデントのt分布の作図 - からっぽのしょこ」を参照してください。
$\lambda$の事後分布のパラメータ$\hat{a}, \hat{b}$、または観測データ$\mathbf{X}$と$\lambda$の事前分布のパラメータ$a, b$を用いて、予測分布(スチューデントのt分布)のパラメータ$\mu_s, \hat{\lambda}_s, \hat{\nu}_s$を計算します。
# 予測分布のパラメータを計算:式(3.79) mu_st <- mu lambda_st_hat <- a_hat / b_hat nu_st_hat <- 2 * a_hat #lambda_st_hat <- (N + 2 * a) / (sum((x_n - mu)^2) + 2 * b) #nu_st_hat <- N + 2 * a mu_st; lambda_st_hat; nu_st_hat
## [1] 25 ## [1] 0.01085793 ## [1] 52
予測分布のパラメータは、次の式で計算できます。
それぞれ1つ目の式だと事後分布のパラメータa_hat, b_hat、2つ目の式だと観測データx_nと事前分布のパラメータa, bを使って計算できます。
予測分布を計算します。
# 予測分布を計算:式(3.76) predict_df <- tibble::tibble( x = x_vec, # 確率変数 dens = LaplacesDemon::dst(x = x_vec, mu = mu_st, sigma = 1/sqrt(lambda_st_hat), nu = nu_st_hat) # 確率密度 ) predict_df
## # A tibble: 501 × 2 ## x dens ## <dbl> <dbl> ## 1 -15 0.0000199 ## 2 -14.8 0.0000210 ## 3 -14.7 0.0000221 ## 4 -14.5 0.0000233 ## 5 -14.4 0.0000246 ## 6 -14.2 0.0000259 ## 7 -14.0 0.0000273 ## 8 -13.9 0.0000288 ## 9 -13.7 0.0000304 ## 10 -13.6 0.0000320 ## # … with 491 more rows
スチューデントのt分布の確率密度は、LaplacesDemonパッケージのdst()で計算できます。確率変数の引数xにx_vec、位置パラメータの引数muにmu_st、スケールパラメータの引数sigmaにlambda_st_hatの平方根の逆数、形状パラメータ(自由度)の引数nuに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 = paste0( "list(N==", N, ", mu[s]==", round(mu_st, 2), ", hat(lambda)[s]==", round(lambda_st_hat, 5), ", hat(nu)[s]==", nu_st_hat, ")" ) ), x = "x", y = "density")
観測データが増えると、予測分布が真の分布に近付きます。
以上で、ガウス分布のベイズ推論を実装できました。
学習推移の可視化
gganimateパッケージを利用して、パラメータの推定値(更新値)の推移(分布の変化)をアニメーション(gif画像)で確認します。
モデルの設定
真の分布のパラメータ$\mu, \lambda$と、$\lambda$の事前分布のパラメータ(の初期値)$a, b$を設定します。
# 既知のパラメータを指定 mu <- 25 # 真の精度パラメータを指定 lambda_truth <- 0.01 # λの事前分布のパラメータを指定 a <- 1 b <- 1 # データ数(試行回数)を指定 N <- 100 # グラフ用のλの値を作成 lambda_vec <- seq(0, lambda_truth * 5, length.out = 501) # グラフ用のxの値を作成 x_vec <- seq( mu - 1/sqrt(lambda_truth) * 4, mu + 1/sqrt(lambda_truth) * 4, length.out = 501 )
実装時と同じく、パラメータを指定して、確率変数がとり得る値(x軸の値)を作成します。グラフが粗い場合や処理が重い場合は、mu_vec, x_vecの要素数(by引数の値)または値の間隔(length.out引数の値)を調整してください。
事後分布と予測分布のパラメータの計算(更新)処理について、2通りの方法を載せます。目的に応じて使い分けてください。
推論処理：for関数による処理
1つ目の処理方法では、for()を使って、1データずつ生成してパラメータの更新を繰り返し処理します。こちらの方が、前ステップで求めた事後分布(のパラメータ)を次ステップの事前分布(のパラメータ)として用いる逐次学習をイメージしやすいです。
・コード(クリックで展開)
超パラメータの初期値を使って、事前分布を計算します。
# λの事前分布(ガンマ分布)を計算:式(2.56) anime_posterior_df <- tibble::tibble( lambda = lambda_vec, # 確率変数 dens = dgamma(x = lambda_vec, shape = a, rate = b), # 確率密度 param = paste0("N=", 0, ", a=", a, ", b=", b) |> as.factor() # フレーム切替用ラベル ) anime_posterior_df
## # A tibble: 501 × 3 ## lambda dens param ## <dbl> <dbl> <fct> ## 1 0 1 N=0, a=1, b=1 ## 2 0.0001 1.00 N=0, a=1, b=1 ## 3 0.0002 1.00 N=0, a=1, b=1 ## 4 0.0003 1.00 N=0, a=1, b=1 ## 5 0.0004 1.00 N=0, a=1, b=1 ## 6 0.0005 1.00 N=0, a=1, b=1 ## 7 0.0006 0.999 N=0, a=1, b=1 ## 8 0.0007 0.999 N=0, a=1, b=1 ## 9 0.0008 0.999 N=0, a=1, b=1 ## 10 0.0009 0.999 N=0, a=1, b=1 ## # … with 491 more rows
現在のパラメータを文字列結合し因子型に変換して、フレーム切替用のラベル列とします。
また、事前分布のパラメータを使って、予測分布のパラメータと予測分布を計算します。
# 初期値による予測分布のパラメータを計算:式(3.79) mu_st <- mu lambda_st <- a / 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=", round(lambda_st, 5), ", nu_s=", round(nu_st, 2)) |> as.factor() # フレーム切替用ラベル ) anime_predict_df
## # A tibble: 501 × 3 ## x dens param ## <dbl> <dbl> <fct> ## 1 -15 0.0000156 N=0, mu_s=25, lambda_s=1, nu_s=2 ## 2 -14.8 0.0000158 N=0, mu_s=25, lambda_s=1, nu_s=2 ## 3 -14.7 0.0000160 N=0, mu_s=25, lambda_s=1, nu_s=2 ## 4 -14.5 0.0000162 N=0, mu_s=25, lambda_s=1, nu_s=2 ## 5 -14.4 0.0000164 N=0, mu_s=25, lambda_s=1, nu_s=2 ## 6 -14.2 0.0000166 N=0, mu_s=25, lambda_s=1, nu_s=2 ## 7 -14.0 0.0000168 N=0, mu_s=25, lambda_s=1, nu_s=2 ## 8 -13.9 0.0000170 N=0, mu_s=25, lambda_s=1, nu_s=2 ## 9 -13.7 0.0000172 N=0, mu_s=25, lambda_s=1, nu_s=2 ## 10 -13.6 0.0000174 N=0, mu_s=25, lambda_s=1, nu_s=2 ## # … with 491 more rows
同様に、フレーム切替用のラベル列を作成します。
パラメータの更新処理をN回繰り返します。
# 観測データの受け皿を作成 x_n <- rep(NA, times = N) # ベイズ推論 for(n in 1:N){ # ガウス分布に従うデータを生成 x_n[n] <- rnorm(n = 1, mean = mu, sd = 1/sqrt(lambda_truth)) # λの事後分布のパラメータを更新:式(3.69) a <- 1 / 2 + a b <- 0.5 * (x_n[n] - mu)^2 + b # λの事後分布(ガンマ分布)を計算:式(2.56) tmp_posterior_df <- tibble::tibble( lambda = lambda_vec, # 確率変数 dens = dgamma(x = lambda, shape = a, rate = b), # 確率密度 param = paste0("N=", n, ", a=", a, ", b=", round(b, 1)) |> as.factor() # フレーム切替用ラベル ) # 予測分布のパラメータを更新:式(3.79) mu_st <- mu lambda_st <- a / 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=", mu_st, ", lambda_s=", round(lambda_st, 5), ", nu_s=", round(nu_st, 2)) |> as.factor() # フレーム切替用ラベル ) # 推論結果を結合 anime_posterior_df <- rbind(anime_posterior_df, tmp_posterior_df) anime_predict_df <- rbind(anime_predict_df, tmp_predict_df) }
超パラメータに関して、$\hat{a}, \hat{b}$に対応するa_hat, b_hatを新たに作るのではなく、a, bを繰り返し更新(上書き)していきます。これにより、事後分布のパラメータの計算式(3.69)の$\frac{N}{2}$と$\sum_{n=1}^N (x_n - \mu)^2$の計算は、ループ処理によってN回繰り返し0.5と(x_n[n] - mu)^2を加えることで行います。n回目のループ処理のときには、n-1回分の0.5と(x_n[n] - mu)^2が既にaとbに加えられています。
更新したパラメータを使って事後分布と予測分布を計算して、それぞれ計算結果をanime_***_dfに結合していきます。
結果を確認します。
# 確認 head(x_n); anime_posterior_df; anime_predict_df
## [1] 8.914162 19.838103 25.662442 15.688073 24.953831 27.057376 ## # A tibble: 50,601 × 3 ## lambda dens param ## <dbl> <dbl> <fct> ## 1 0 1 N=0, a=1, b=1 ## 2 0.0001 1.00 N=0, a=1, b=1 ## 3 0.0002 1.00 N=0, a=1, b=1 ## 4 0.0003 1.00 N=0, a=1, b=1 ## 5 0.0004 1.00 N=0, a=1, b=1 ## 6 0.0005 1.00 N=0, a=1, b=1 ## 7 0.0006 0.999 N=0, a=1, b=1 ## 8 0.0007 0.999 N=0, a=1, b=1 ## 9 0.0008 0.999 N=0, a=1, b=1 ## 10 0.0009 0.999 N=0, a=1, b=1 ## # … with 50,591 more rows ## # A tibble: 50,601 × 3 ## x dens param ## <dbl> <dbl> <fct> ## 1 -15 0.0000156 N=0, mu_s=25, lambda_s=1, nu_s=2 ## 2 -14.8 0.0000158 N=0, mu_s=25, lambda_s=1, nu_s=2 ## 3 -14.7 0.0000160 N=0, mu_s=25, lambda_s=1, nu_s=2 ## 4 -14.5 0.0000162 N=0, mu_s=25, lambda_s=1, nu_s=2 ## 5 -14.4 0.0000164 N=0, mu_s=25, lambda_s=1, nu_s=2 ## 6 -14.2 0.0000166 N=0, mu_s=25, lambda_s=1, nu_s=2 ## 7 -14.0 0.0000168 N=0, mu_s=25, lambda_s=1, nu_s=2 ## 8 -13.9 0.0000170 N=0, mu_s=25, lambda_s=1, nu_s=2 ## 9 -13.7 0.0000172 N=0, mu_s=25, lambda_s=1, nu_s=2 ## 10 -13.6 0.0000174 N=0, mu_s=25, lambda_s=1, nu_s=2 ## # … with 50,591 more rows
それぞれ「lambda_vec, x_vecの要素数」掛ける「N+1」行のデータフレームになります。行数が増えすぎないように注意してください。
推論処理：tidyverseパッケージによる処理
2つ目の処理方法では、tidyverseパッケージの関数を使って、一度の処理でN+1回分の計算をします。こちらの方が、処理時間が短いです。
・コード(クリックで展開)
ガウス分布に従う$N$個のデータを生成します。
# ガウス分布に従うデータを生成 x_n <- rnorm(n = N, mean = mu, sd = 1/sqrt(lambda_truth)) head(x_n)
## [1] 8.914162 19.838103 25.662442 15.688073 24.953831 27.057376
事前分布を含めたN+1回分の事後分布を計算します。
# 試行ごとに事後分布(ガンマ分布)を計算 anime_posterior_df <- tidyr::expand_grid( n = 0:N, # 試行回数 lambda = lambda_vec # 確率変数 ) |> # 全ての組み合わせを作成 dplyr::mutate( a = 0.5*n + a, b = c(b, 0.5*cumsum((x_n - mu)^2) + b)[n+1] ) |> # 事後分布のパラメータを計算:式(3.69) dplyr::mutate( dens = dgamma(x = lambda, shape = a, rate = b), # 確率密度 param = paste0("N=", n, ", a=", a, ", b=", round(b, 1)) |> (\(.){factor(., levels = unique(.))})() # フレーム切替用ラベル ) # 事後分布を計算:式(2.56) anime_posterior_df
## # A tibble: 50,601 × 6 ## n lambda a b dens param ## <int> <dbl> <dbl> <dbl> <dbl> <fct> ## 1 0 0 1 1 1 N=0, a=1, b=1 ## 2 0 0.0001 1 1 1.00 N=0, a=1, b=1 ## 3 0 0.0002 1 1 1.00 N=0, a=1, b=1 ## 4 0 0.0003 1 1 1.00 N=0, a=1, b=1 ## 5 0 0.0004 1 1 1.00 N=0, a=1, b=1 ## 6 0 0.0005 1 1 1.00 N=0, a=1, b=1 ## 7 0 0.0006 1 1 0.999 N=0, a=1, b=1 ## 8 0 0.0007 1 1 0.999 N=0, a=1, b=1 ## 9 0 0.0008 1 1 0.999 N=0, a=1, b=1 ## 10 0 0.0009 1 1 0.999 N=0, a=1, b=1 ## # … with 50,591 more rows
試行回数(観測データ数)を表す0からNまでのN+1個の整数と、x軸の値lambda_vecの全ての組み合わせをexpand_grid()で作成します。これにより、lambda_vecの各要素をN+1フレーム分に複製できます。
n列を使ってaを計算します。事前分布とN回分の事後分布のパラメータb(のベクトル)を計算(作成)して、n列の値をインデックスとして使って各試行に対応する値を抽出します。
確率変数とパラメータの組み合わせごとに確率密度を計算して、パラメータごとにフレーム切替用のラベルを作成します。
因子型への変換処理では、無名関数function()の省略記法\()を使って、(\(引数){引数を使った具体的な処理})()としています。直前のパイプ演算子を%>%にすると、行全体(\(引数){処理})()を{}の中の処理(この例だとfactor(., levels = unique(.)))に置き換えられます(置き換えられるように引数名を.にしています)。
同様に、初期値を含めたN+1回分の予測分布を計算します。
# 試行ごとに予測分布(スチューデントのt分布)を計算 anime_predict_df <- tidyr::expand_grid( n = 0:N, # 試行回数 x = x_vec # 確率変数 ) |> # 全ての組み合わせを作成 dplyr::mutate( mu_st = mu, lambda_st = c(a/b, (1:N + 2*a) / (cumsum((x_n - mu)^2) + 2*b))[n+1], nu_st = n + 2*a ) |> # 予測分布のパラメータを計算:式(3.79) 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=", mu_st, ", lambda_s=", round(lambda_st, 5), ", nu_s=", round(nu_st, 2)) |> (\(.){factor(., levels = unique(.))})() # フレーム切替用ラベル ) # 事後分布を計算:式(3.76) anime_predict_df
## # A tibble: 50,601 × 7 ## n x mu_st lambda_st nu_st dens param ## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <fct> ## 1 0 -15 25 1 2 0.0000156 N=0, mu_s=25, lambda_s=1, nu_s=2 ## 2 0 -14.8 25 1 2 0.0000158 N=0, mu_s=25, lambda_s=1, nu_s=2 ## 3 0 -14.7 25 1 2 0.0000160 N=0, mu_s=25, lambda_s=1, nu_s=2 ## 4 0 -14.5 25 1 2 0.0000162 N=0, mu_s=25, lambda_s=1, nu_s=2 ## 5 0 -14.4 25 1 2 0.0000164 N=0, mu_s=25, lambda_s=1, nu_s=2 ## 6 0 -14.2 25 1 2 0.0000166 N=0, mu_s=25, lambda_s=1, nu_s=2 ## 7 0 -14.0 25 1 2 0.0000168 N=0, mu_s=25, lambda_s=1, nu_s=2 ## 8 0 -13.9 25 1 2 0.0000170 N=0, mu_s=25, lambda_s=1, nu_s=2 ## 9 0 -13.7 25 1 2 0.0000172 N=0, mu_s=25, lambda_s=1, nu_s=2 ## 10 0 -13.6 25 1 2 0.0000174 N=0, mu_s=25, lambda_s=1, nu_s=2 ## # … with 50,591 more rows
0からNまでの整数と、x軸の値x_vecの全ての組み合わせを作成します。
mu_stはmuです。N+1回分の予測分布パラメータlambda_st(のベクトル)を計算(作成)して、各試行に対応する値を抽出します。nu_stについてはn列を使って計算できます。
確率変数とパラメータの組み合わせごとに確率を計算して、パラメータごとにフレーム切替用のラベルを作成します。
作図処理
事後分布と予測分布の推移をそれぞれアニメーションで可視化します。
・コード(クリックで展開)
観測データと対応するラベルをデータフレームに格納します。
# 観測データを格納 anime_data_df <- tibble::tibble( scaled_x = c(NA, 1/(x_n - mu)^2), # 偏差の2乗の逆数に変換 param = unique(anime_posterior_df[["param"]]) # フレーム切替用ラベル ) anime_data_df
## # A tibble: 101 × 2 ## scaled_x param ## <dbl> <fct> ## 1 NA N=0, a=1, b=1 ## 2 0.00386 N=1, a=1.5, b=130.4 ## 3 0.0375 N=2, a=2, b=143.7 ## 4 2.28 N=3, a=2.5, b=143.9 ## 5 0.0115 N=4, a=3, b=187.3 ## 6 469. N=5, a=3.5, b=187.3 ## 7 0.236 N=6, a=4, b=189.4 ## 8 0.0282 N=7, a=4.5, b=207.1 ## 9 0.00262 N=8, a=5, b=398.3 ## 10 0.00353 N=9, a=5.5, b=539.8 ## # … with 91 more rows
事前分布には観測データが影響しないので欠損値NAとして、anime_posterior_dfのラベル列と合わせて格納します。観測データと精度パラメータを対応させるため、各データから平均$\mu$を引き、2乗の逆数にして格納します。
参考として、各試行における観測データを描画するのに使います。
事後分布の推移のアニメーションを作成します。
# λの事後分布のアニメーションを作図 posterior_graph <- ggplot() + geom_vline(mapping = aes(xintercept = lambda_truth, color = "param"), size = 1, linetype = "dashed", show.legend = FALSE) + # 真のパラメータ geom_line(data = anime_posterior_df, mapping = aes(x = lambda, y = dens, color = "posterior"), size = 1) + # λの事後分布 geom_point(data = anime_data_df, mapping = aes(x = scaled_x, y = 0, color = "data"), size = 6) + # 観測データ gganimate::transition_manual(param) + # フレーム scale_color_manual(breaks = c("param", "posterior", "data"), values = c("red", "purple", "pink"), labels = c("true parameter", "posterior", "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(xlim = c(0, max(lambda_vec))) + # 軸の表示範囲 labs(title = "Gamma Distribution", subtitle = "{current_frame}", x = expression(lambda), y = "density") # 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=25, lambda_s=1, nu_s=2 ## 2 8.91 N=1, mu_s=25, lambda_s=0.01151, nu_s=3 ## 3 19.8 N=2, mu_s=25, lambda_s=0.01392, nu_s=4 ## 4 25.7 N=3, mu_s=25, lambda_s=0.01737, nu_s=5 ## 5 15.7 N=4, mu_s=25, lambda_s=0.01602, nu_s=6 ## 6 25.0 N=5, mu_s=25, lambda_s=0.01869, nu_s=7 ## 7 27.1 N=6, mu_s=25, lambda_s=0.02112, nu_s=8 ## 8 19.0 N=7, mu_s=25, lambda_s=0.02173, nu_s=9 ## 9 44.6 N=8, mu_s=25, lambda_s=0.01255, nu_s=10 ## 10 41.8 N=9, mu_s=25, lambda_s=0.01019, 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 8.91 N=2, mu_s=25, lambda_s=0.01392, nu_s=4 ## 2 3 1 8.91 N=3, mu_s=25, lambda_s=0.01737, nu_s=5 ## 3 3 2 19.8 N=3, mu_s=25, lambda_s=0.01737, nu_s=5 ## 4 4 1 8.91 N=4, mu_s=25, lambda_s=0.01602, nu_s=6 ## 5 4 2 19.8 N=4, mu_s=25, lambda_s=0.01602, nu_s=6 ## 6 4 3 25.7 N=4, mu_s=25, lambda_s=0.01602, nu_s=6 ## 7 5 1 8.91 N=5, mu_s=25, lambda_s=0.01869, nu_s=7 ## 8 5 2 19.8 N=5, mu_s=25, lambda_s=0.01869, nu_s=7 ## 9 5 3 25.7 N=5, mu_s=25, lambda_s=0.01869, nu_s=7 ## 10 5 4 15.7 N=5, mu_s=25, lambda_s=0.01869, nu_s=7 ## # … with 4,940 more rows
フレーム番号と試行回数(データ番号)の全ての組み合わせを作成して、フレーム番号未満のデータを抽出します。フレーム切替用のラベルは、フレーム番号を使って抽出します。
観測データの無い初回(0回目)と1回目は描画しないので、フレーム番号が2からになります。
予測分布の推移のアニメーションを作成します。
# 真の分布を計算 model_df <- tibble::tibble( x = x_vec, # 確率変数 dens = dnorm(x, mean = mu, 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.17：加筆修正の際に、数式読解編とRで実装編に記事を分割しました。
ガウスなのかガンマなのか字面が似ていてこんがらがってくる。
- 2022.08.08：加筆修正しました。
【次の節の内容】