はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」によって、「数式」と「プログラム」から理解するのが目標です。
省略してある内容等ありますので、本とあわせて読んでください。
この記事では、「尤度関数を1次元ガウス分布」、「事前分布を1次元ガウス分布」とした場合の「パラメータの事後分布」と「未観測値の予測分布」の計算をR言語で実装します。
【数式読解編】
【他の節の内容】
【この節の内容】
3.3.1 1次元ガウス分布の学習と予測：平均が未知の場合
尤度関数を1次元ガウス分布(Gaussian Distribution)・一変量正規分布(Normal Distribution)、事前分布を1次元ガウス分布とするモデルに対するベイズ推論を実装します。人工的に生成したデータを用いて、ガウス分布の平均パラメータを推定し、また未観測データに対する予測分布を求めます。
ガウス分布については「1次元ガウス分布の定義式の確認 - からっぽのしょこ」を参照してください。
利用するパッケージを読み込みます。
# 利用パッケージ library(tidyverse) library(gganimate)
この記事では、基本的にパッケージ名::関数名()の記法を使うので、パッケージを読み込む必要はありません。ただし、作図コードがごちゃごちゃしないようにパッケージ名を省略しているため、ggplot2は読み込む必要があります。
magrittrパッケージのパイプ演算子%>%ではなく、ネイティブパイプ(ベースパイプ)演算子|>を使っています。%>%に置き換えても処理できます。
学習の推移をアニメーション(gif画像)で確認するのにgganimateパッケージを利用します。不要であれば省略してください。
ベイズ推論の実装
まずは、モデルを設定して、データを生成します。生成したデータを用いて、事後分布のパラメータを計算します。さらに、事後分布のパラメータを用いて、予測分布のパラメータを計算します。
生成分布の設定
データ生成分布(真の分布)として、ガウス分布$\mathcal{N}(x | \mu, \lambda^{-1})$を設定します。
ガウス分布のグラフ作成については「【R】1次元ガウス分布の作図 - からっぽのしょこ」を参照してください。
真の分布(ガウス分布)のパラメータ$\mu, \lambda$を設定します。
# 真の平均パラメータを指定 mu_truth <- 25 # 既知の精度パラメータを指定 lambda <- 0.01 sqrt(1 / lambda) # 標準偏差
## [1] 10
平均パラメータ(実数)$\mu$をmu_truthとして値を指定します。mu_truthが真のパラメータであり、この値を求めるのがここでの目的です。
精度パラメータ(正の実数)$\lambda$をlambdaとして値を指定します。こちらは与えられている値として使います。精度は分散の逆数なので、値が大きいほど散らばり具合が小さくなり、逆数の平方根が標準偏差$\sigma = \frac{1}{\sqrt{\lambda}}$です。
真の分布を計算して、作図用のデータフレームを作成します。
# グラフ用のxの値を作成 x_vec <- seq( mu_truth - 1/sqrt(lambda) * 4, mu_truth + 1/sqrt(lambda) * 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)) # 確率密度 ) 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の平方根の逆数を指定します。
真の分布のグラフを作成します。
# 真の分布を作図 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, ")")), 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)) head(x_n)
## [1] 23.42248 30.08129 21.61894 28.09246 17.31935 26.12295
生成するデータ数$N$を指定します。
ガウス分布に従う乱数は、rnorm()で生成できます。試行回数の引数nにN、平均の引数meanにmu_truth、標準偏差の引数sdにlambdaの平方根の逆数を指定します。生成したN個のデータをx_nとします。
観測データのヒストグラムを真の分布と重ねて確認します。
# 観測データを格納 data_df <- tibble::tibble(x = x_n) # 観測データのヒストグラムを作成 ggplot() + geom_histogram(data = data_df, 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, ", N==", N, ")")), x = "x", y = "density")
観測データx_nをデータフレームに格納して、geom_histogram()でヒストグラムを描画します。y軸の引数yに..density..を指定すると、度数を密度に変換して描画します。
データ数が十分に大きいと、分布の形状が真の分布に近付きます。
事前分布の設定
次は、尤度に対する共役事前分布を設定します。ガウス分布の平均パラメータ$\mu$の事前分布としてガウス分布$\mathcal{N}(\mu | m, \lambda_{\mu})$を設定します。
$\mu$の事前分布(ガウス分布)のパラメータ(超パラメータ)$m, \lambda_{\mu}$を設定します。
# μの事前分布の平均パラメータを指定 m <- 0 # μの事前分布の精度パラメータを指定 lambda_mu <- 0.0016 sqrt(1 / lambda_mu) # 標準偏差
## [1] 25
平均パラメータ(実数)$m$と精度パラメータ(正の実数)$\lambda_{\mu}$の値を指定します。
$\mu$の事前分布を計算します。
# グラフ用のμの値を作成 mu_vec <- seq( mu_truth - 1/sqrt(lambda_mu) * 2, mu_truth + 1/sqrt(lambda_mu) * 2, length.out = 501 ) # μの事前分布を計算:式(2.64) prior_df <- tibble::tibble( mu = mu_vec, # 確率変数 dens = dnorm(x = mu_vec, mean = m, sd = 1/sqrt(lambda_mu)) # 確率密度 ) prior_df
## # A tibble: 501 × 2 ## mu dens ## <dbl> <dbl> ## 1 -25 0.00968 ## 2 -24.8 0.00976 ## 3 -24.6 0.00983 ## 4 -24.4 0.00991 ## 5 -24.2 0.00999 ## 6 -24 0.0101 ## 7 -23.8 0.0101 ## 8 -23.6 0.0102 ## 9 -23.4 0.0103 ## 10 -23.2 0.0104 ## # … with 491 more rows
真の分布と同様に、ガウス分布の確率変数がとり得る値$\mu$を作成して、値ごとに確率密度を計算します。この例では、真の値を中心に標準偏差の2倍を範囲とします。
$\mu$の事前分布のグラフを作成します。
# μの事前分布を作図 ggplot() + geom_vline(mapping = aes(xintercept = mu_truth, color = "param"), size = 1, linetype = "dashed", show.legend = FALSE) + # 真のパラメータ geom_line(data = prior_df, mapping = aes(x = mu, 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 = "Gaussian Distribution", subtitle = parse(text = paste0("list(m==", m, ", lambda[mu]==", lambda_mu, ")")), x = expression(mu), y = "density")
真のパラメータの値を赤色の破線で示します。
事後分布の計算
観測データ$\mathbf{X}$から平均パラメータ$\mu$の事後分布$\mathcal{N}(\mu | \hat{m}, \hat{\lambda}_{\mu})$を求めます(平均パラメータ$\mu$を分布推定します)。
観測データ$\mathbf{X}$を用いて、$\mu$の事後分布(ガウス分布)のパラメータ$\hat{m}, \hat{\lambda}_{\mu}$を計算します。
# μの事後分布のパラメータを計算:式(3.53),(3.54) lambda_mu_hat <- N * lambda + lambda_mu m_hat <- (lambda * sum(x_n) + lambda_mu * m) / lambda_mu_hat m_hat; lambda_mu_hat
## [1] 25.4732 ## [1] 0.5016
$\mu$の事後分布のパラメータは、次の式で計算できます。
$\mu$の事後分布を計算します。
# μの事後分布を計算:式(2.64) posterior_df <- tibble::tibble( mu = mu_vec, # 確率変数 dens = dnorm(x = mu_vec, mean = m_hat, sd = 1/sqrt(lambda_mu_hat)) # 確率密度 ) posterior_df
## # A tibble: 501 × 2 ## mu dens ## <dbl> <dbl> ## 1 -25 9.33e-279 ## 2 -24.8 1.46e-276 ## 3 -24.6 2.24e-274 ## 4 -24.4 3.37e-272 ## 5 -24.2 4.97e-270 ## 6 -24 7.18e-268 ## 7 -23.8 1.02e-265 ## 8 -23.6 1.41e-263 ## 9 -23.4 1.92e-261 ## 10 -23.2 2.56e-259 ## # … with 491 more rows
更新した超パラメータm_hat, lambda_mu_hatを使って、事前分布のときと同様にして処理します。
$\mu$の事後分布のグラフを作図します。
# μの事後分布を作図 ggplot() + geom_vline(aes(xintercept = mu_truth, color = "param"), size = 1, linetype = "dashed", show.legend = FALSE) + # 真のパラメータ geom_line(data = posterior_df, mapping = aes(x = mu, 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 = "Gaussian Distribution", subtitle = parse( text = paste0("list(N==", N, ", hat(m)==", round(m_hat, 2), ", hat(lambda)[mu]==", round(lambda_mu_hat, 5), ")") ), x = expression(mu), y = "density")
$\mu$の真の値付近をピークとする分布になっています。
予測分布の計算
最後に、観測データ$\mathbf{X}$または(観測データから求めた)事後分布のパラメータ$\hat{m}, \hat{\lambda}_{\mu}$から未観測のデータ$x_{*}$の予測分布$\mathcal{N}(x_{*} | \hat{\mu}_{*}, \hat{\lambda}_{*})$を求めます。
$\mu$の事後分布のパラメータ$\hat{m}, \hat{\lambda}_{\mu}$、または観測データ$\mathbf{X}$と$\mu$の事前分布のパラメータ$m, \lambda_{\mu}$を用いて、予測分布(ガウス分布)のパラメータ$\hat{\mu}_{*}, \hat{\lambda}_{*}$を計算します。
# 予測分布のパラメータを計算:式(3.62') lambda_star_hat <- lambda * lambda_mu_hat / (lambda + lambda_mu_hat) mu_star_hat <- m_hat #lambda_star_hat <- (N * lambda + lambda_mu) * lambda / ((N + 1) * lambda + lambda_mu) #mu_star_hat <- (lambda * sum(x_n) + lambda_mu * m) / (N * lambda + lambda_mu) mu_star_hat; lambda_star_hat
## [1] 25.4732 ## [1] 0.009804535
予測分布のパラメータは、次の式で計算できます。
それぞれ1つ目の式だと事後分布のパラメータm_hat, lambda_mu_hat、2つ目の式だと、観測データx_nと事前分布のパラメータm, lambda_muを使って計算できます。
予測分布を計算します。
# 予測分布を計算:式(2.64) predict_df <- tibble::tibble( x = x_vec, # 確率変数 dens = dnorm(x = x_vec, mean = mu_star_hat, sd = 1/sqrt(lambda_star_hat)) # 確率密度 ) predict_df
## # A tibble: 501 × 2 ## x dens ## <dbl> <dbl> ## 1 -15 0.0000129 ## 2 -14.8 0.0000137 ## 3 -14.7 0.0000146 ## 4 -14.5 0.0000155 ## 5 -14.4 0.0000165 ## 6 -14.2 0.0000176 ## 7 -14.0 0.0000187 ## 8 -13.9 0.0000199 ## 9 -13.7 0.0000212 ## 10 -13.6 0.0000225 ## # … with 491 more rows
真の分布のときと同様にして処理します。
予測分布のグラフを真の分布と重ねて作成します。
# 予測分布を作図 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 = "Gaussian Distribution", subtitle = parse( text = paste0("list(N==", N, ", hat(mu)[s]==", round(mu_star_hat, 2), ", hat(lambda)[s]==", round(lambda_star_hat, 5), ")") ), x = "x", y = "density")
観測データが増えると、予測分布が真の分布に近付きます。
以上で、ガウス分布のベイズ推論を実装できました。
学習推移の可視化
gganimateパッケージを利用して、パラメータの推定値(更新値)の推移(分布の変化)をアニメーション(gif画像)で確認します。
モデルの設定
真の分布のパラメータ$\mu, \lambda$と、$\mu$の事前分布のパラメータ(の初期値)$m, \lambda_{\mu}$を設定します。
# 真の平均パラメータを指定 mu_truth <- 25 # 既知の精度パラメータを指定 lambda <- 0.01 # μの事前分布の平均パラメータを指定 m <- 0 # μの事前分布の精度パラメータを指定 lambda_mu <- 0.0016 # データ数(試行回数)を指定 N <- 100 # グラフ用のμの値を作成 mu_vec <- seq( mu_truth - 1/sqrt(lambda_mu) * 2, mu_truth + 1/sqrt(lambda_mu) * 2, length.out = 501 ) # グラフ用のxの値を作成 x_vec <- seq( mu_truth - 1/sqrt(lambda) * 4, mu_truth + 1/sqrt(lambda) * 4, length.out = 501 )
実装時と同じく、パラメータを指定して、確率変数がとり得る値(x軸の値)を作成します。グラフが粗い場合や処理が重い場合は、mu_vec, x_vecの要素数(by引数の値)または値の間隔(length.out引数の値)を調整してください。
事後分布と予測分布のパラメータの計算(更新)処理について、2通りの方法を載せます。目的に応じて使い分けてください。
推論処理：for関数による処理
1つ目の処理方法では、for()を使って、1データずつ生成してパラメータの更新を繰り返し処理します。こちらの方が、前ステップで求めた事後分布(のパラメータ)を次ステップの事前分布(のパラメータ)として用いる逐次学習をイメージしやすいです。
・コード(クリックで展開)
超パラメータの初期値を使って、事前分布を計算します。
# μの事前分布(ガウス分布)を計算:式(2.64) anime_posterior_df <- tibble::tibble( mu = mu_vec, # 確率変数 dens = dnorm(x = mu_vec, mean = m, sd = 1/sqrt(lambda_mu)), # 確率密度 param = paste0("N=", 0, ", m=", m, ", lambda_mu=", lambda_mu) |> as.factor() # フレーム切替用ラベル ) anime_posterior_df
## # A tibble: 501 × 3 ## mu dens param ## <dbl> <dbl> <fct> ## 1 -25 0.00968 N=0, m=0, lambda_mu=0.0016 ## 2 -24.8 0.00976 N=0, m=0, lambda_mu=0.0016 ## 3 -24.6 0.00983 N=0, m=0, lambda_mu=0.0016 ## 4 -24.4 0.00991 N=0, m=0, lambda_mu=0.0016 ## 5 -24.2 0.00999 N=0, m=0, lambda_mu=0.0016 ## 6 -24 0.0101 N=0, m=0, lambda_mu=0.0016 ## 7 -23.8 0.0101 N=0, m=0, lambda_mu=0.0016 ## 8 -23.6 0.0102 N=0, m=0, lambda_mu=0.0016 ## 9 -23.4 0.0103 N=0, m=0, lambda_mu=0.0016 ## 10 -23.2 0.0104 N=0, m=0, lambda_mu=0.0016 ## # … with 491 more rows
現在のパラメータを文字列結合し因子型に変換して、フレーム切替用のラベル列とします。
また、事前分布のパラメータを使って、予測分布のパラメータと予測分布を計算します。
# 初期値による予測分布のパラメータを計算:式(3.62) lambda_star <- lambda * lambda_mu / (lambda + lambda_mu) mu_star <- m # 初期値による予測分布(ガウス分布)を計算:式(2.64) anime_predict_df <- tibble::tibble( x = x_vec, # 確率変数 dens = dnorm(x = x_vec, mean = mu_star, sd = 1/sqrt(lambda_star)), # 確率密度 param = paste0("N=", 0, ", mu_star=", mu_star, ", lambda_star=", round(lambda_star, 5)) |> as.factor() # フレーム切替用ラベル ) anime_predict_df
## # A tibble: 501 × 3 ## x dens param ## <dbl> <dbl> <fct> ## 1 -15 0.0127 N=0, mu_star=0, lambda_star=0.00138 ## 2 -14.8 0.0127 N=0, mu_star=0, lambda_star=0.00138 ## 3 -14.7 0.0128 N=0, mu_star=0, lambda_star=0.00138 ## 4 -14.5 0.0128 N=0, mu_star=0, lambda_star=0.00138 ## 5 -14.4 0.0129 N=0, mu_star=0, lambda_star=0.00138 ## 6 -14.2 0.0129 N=0, mu_star=0, lambda_star=0.00138 ## 7 -14.0 0.0129 N=0, mu_star=0, lambda_star=0.00138 ## 8 -13.9 0.0130 N=0, mu_star=0, lambda_star=0.00138 ## 9 -13.7 0.0130 N=0, mu_star=0, lambda_star=0.00138 ## 10 -13.6 0.0131 N=0, mu_star=0, lambda_star=0.00138 ## # … with 491 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)) # μの事後分布のパラメータを更新:式(3.53),(3.54) lambda_mu_old <- lambda_mu lambda_mu <- lambda + lambda_mu m <- (lambda * x_n[n] + lambda_mu_old * m) / lambda_mu # μの事後分布(ガウス分布)を計算:式(2.64) tmp_posterior_df <- tibble::tibble( mu = mu_vec, # 確率変数 dens = dnorm(x = mu_vec, mean = m, sd = 1/sqrt(lambda_mu)), # 確率密度 param = paste0("N=", n, ", m=", round(m, 2), ", lambda_mu=", round(lambda_mu, 5)) |> as.factor() # フレーム切替用ラベル ) # 予測分布のパラメータを更新:式(3.62) lambda_star <- lambda * lambda_mu / (lambda + lambda_mu) mu_star <- m # 予測分布(ガウス分布)を計算:式(2.64) tmp_predict_df <- tibble::tibble( x = x_vec, # 確率変数 dens = dnorm(x = x_vec, mean = mu_star, sd = 1/sqrt(lambda_star)), # 確率密度 param = paste0("N=", n, ", mu_star=", round(mu_star, 2), ", lambda_star=", round(lambda_star, 5)) |> 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{\lambda}_{\mu}$に対応するm_hat, lambda_mu_hatを新たに作るのではなく、m, lambda_muを繰り返し更新(上書き)していきます。これにより、事後分布のパラメータの計算式(3.53,3.54)の$N \lambda$と$\sum_{n=1}^N \lambda x_n$の計算は、ループ処理によってN回繰り返しlambdaとlambda*x_n[n]を加えることで行います。n回目のループ処理のときには、n-1回分のlambdaとlambda*x_n[n]が既にmとlambda_muに加えられています。
更新したパラメータを使って事後分布と予測分布を計算して、それぞれ計算結果をanime_***_dfに結合していきます。
結果を確認します。
# 確認 head(x_n); anime_posterior_df; anime_predict_df
## [1] 9.827869 27.355570 22.516882 37.941007 17.475327 19.148907 ## # A tibble: 50,601 × 3 ## mu dens param ## <dbl> <dbl> <fct> ## 1 -25 0.00968 N=0, m=0, lambda_mu=0.0016 ## 2 -24.8 0.00976 N=0, m=0, lambda_mu=0.0016 ## 3 -24.6 0.00983 N=0, m=0, lambda_mu=0.0016 ## 4 -24.4 0.00991 N=0, m=0, lambda_mu=0.0016 ## 5 -24.2 0.00999 N=0, m=0, lambda_mu=0.0016 ## 6 -24 0.0101 N=0, m=0, lambda_mu=0.0016 ## 7 -23.8 0.0101 N=0, m=0, lambda_mu=0.0016 ## 8 -23.6 0.0102 N=0, m=0, lambda_mu=0.0016 ## 9 -23.4 0.0103 N=0, m=0, lambda_mu=0.0016 ## 10 -23.2 0.0104 N=0, m=0, lambda_mu=0.0016 ## # … with 50,591 more rows ## # A tibble: 50,601 × 3 ## x dens param ## <dbl> <dbl> <fct> ## 1 -15 0.0127 N=0, mu_star=0, lambda_star=0.00138 ## 2 -14.8 0.0127 N=0, mu_star=0, lambda_star=0.00138 ## 3 -14.7 0.0128 N=0, mu_star=0, lambda_star=0.00138 ## 4 -14.5 0.0128 N=0, mu_star=0, lambda_star=0.00138 ## 5 -14.4 0.0129 N=0, mu_star=0, lambda_star=0.00138 ## 6 -14.2 0.0129 N=0, mu_star=0, lambda_star=0.00138 ## 7 -14.0 0.0129 N=0, mu_star=0, lambda_star=0.00138 ## 8 -13.9 0.0130 N=0, mu_star=0, lambda_star=0.00138 ## 9 -13.7 0.0130 N=0, mu_star=0, lambda_star=0.00138 ## 10 -13.6 0.0131 N=0, mu_star=0, lambda_star=0.00138 ## # … with 50,591 more rows
それぞれ「mu_vec, x_vecの要素数」掛ける「N+1」行のデータフレームになります。行数が増えすぎないように注意してください。
推論処理：tidyverseパッケージによる処理
2つ目の処理方法では、tidyverseパッケージの関数を使って、一度の処理でN+1回分の計算をします。こちらの方が、処理時間が短いです。
・コード(クリックで展開)
ガウス分布に従う$N$個のデータを生成します。
# ガウス分布に従うデータを生成 x_n <- rnorm(n = N, mean = mu_truth, sd = 1/sqrt(lambda)) head(x_n)
## [1] 9.827869 27.355570 22.516882 37.941007 17.475327 19.148907
事前分布を含めたN+1回分の事後分布を計算します。
# 試行ごとに事後分布(ガウス分布)を計算 anime_posterior_df <- tidyr::expand_grid( n = 0:N, # 試行回数 mu = mu_vec # 確率変数 ) |> # 全ての組み合わせを作成 dplyr::mutate( m = c(m, (lambda*cumsum(x_n) + lambda_mu*m) / (1:N*lambda + lambda_mu))[n+1], lambda_mu = n*lambda + lambda_mu ) |> # 事後分布のパラメータを計算:式(3.53),(3.54) dplyr::mutate( dens = dnorm(x = mu, mean = m, sd = 1/sqrt(lambda_mu)), # 確率密度 param = paste0("N=", n, ", m=", round(m, 2), ", lambda_mu=", round(lambda_mu, 5)) |> (\(.){factor(., levels = unique(.))})() # フレーム切替用ラベル ) # 事後分布を計算:式(2.64) anime_posterior_df
## # A tibble: 50,601 × 6 ## n mu m lambda_mu dens param ## <int> <dbl> <dbl> <dbl> <dbl> <fct> ## 1 0 -25 0 0.0016 0.00968 N=0, m=0, lambda_mu=0.0016 ## 2 0 -24.8 0 0.0016 0.00976 N=0, m=0, lambda_mu=0.0016 ## 3 0 -24.6 0 0.0016 0.00983 N=0, m=0, lambda_mu=0.0016 ## 4 0 -24.4 0 0.0016 0.00991 N=0, m=0, lambda_mu=0.0016 ## 5 0 -24.2 0 0.0016 0.00999 N=0, m=0, lambda_mu=0.0016 ## 6 0 -24 0 0.0016 0.0101 N=0, m=0, lambda_mu=0.0016 ## 7 0 -23.8 0 0.0016 0.0101 N=0, m=0, lambda_mu=0.0016 ## 8 0 -23.6 0 0.0016 0.0102 N=0, m=0, lambda_mu=0.0016 ## 9 0 -23.4 0 0.0016 0.0103 N=0, m=0, lambda_mu=0.0016 ## 10 0 -23.2 0 0.0016 0.0104 N=0, m=0, lambda_mu=0.0016 ## # … with 50,591 more rows
試行回数(観測データ数)を表す0からNまでのN+1個の整数と、x軸の値mu_vecの全ての組み合わせをexpand_grid()で作成します。これにより、mu_vecの各要素をN+1フレーム分に複製できます。
事前分布とN回分の事後分布のパラメータm(のベクトル)を計算(作成)して、n列の値をインデックスとして使って各試行に対応する値を抽出します。lambda_muについてはn列を使って計算できます。
確率変数とパラメータの組み合わせごとに確率密度を計算して、パラメータごとにフレーム切替用のラベルを作成します。
因子型への変換処理では、無名関数function()の省略記法\()を使って、(\(引数){引数を使った具体的な処理})()としています。直前のパイプ演算子を%>%にすると、行全体(\(引数){処理})()を{}の中の処理(この例だとfactor(., levels = unique(.)))に置き換えられます(置き換えられるように引数名を.にしています)。
同様に、初期値を含めたN+1回分の予測分布を計算します。
# 試行ごとに予測分布(ガウス分布)を計算 anime_predict_df <- tidyr::expand_grid( n = 0:N, # 試行回数 x = x_vec # 確率変数 ) |> # 全ての組み合わせを作成 dplyr::mutate( mu_star = c(m, (lambda*cumsum(x_n) + lambda_mu*m) / (1:N*lambda + lambda_mu))[n+1], lambda_star = lambda*(n * lambda + lambda_mu) / ((n + 1)*lambda + lambda_mu) ) |> # 予測分布のパラメータを計算:式(3.54),(3.62') dplyr::mutate( dens = dnorm(x = x, mean = mu_star, sd = 1/sqrt(lambda_star)), # 確率密度 param = paste0("N=", n, ", mu_star=", round(mu_star, 2), ", lambda_star=", round(lambda_star, 5)) |> (\(.){factor(., levels = unique(.))})() # フレーム切替用ラベル ) # 事後分布を計算:式(2.64) anime_predict_df
## # A tibble: 50,601 × 6 ## n x mu_star lambda_star dens param ## <int> <dbl> <dbl> <dbl> <dbl> <fct> ## 1 0 -15 0 0.00138 0.0127 N=0, mu_star=0, lambda_star=0.00138 ## 2 0 -14.8 0 0.00138 0.0127 N=0, mu_star=0, lambda_star=0.00138 ## 3 0 -14.7 0 0.00138 0.0128 N=0, mu_star=0, lambda_star=0.00138 ## 4 0 -14.5 0 0.00138 0.0128 N=0, mu_star=0, lambda_star=0.00138 ## 5 0 -14.4 0 0.00138 0.0129 N=0, mu_star=0, lambda_star=0.00138 ## 6 0 -14.2 0 0.00138 0.0129 N=0, mu_star=0, lambda_star=0.00138 ## 7 0 -14.0 0 0.00138 0.0129 N=0, mu_star=0, lambda_star=0.00138 ## 8 0 -13.9 0 0.00138 0.0130 N=0, mu_star=0, lambda_star=0.00138 ## 9 0 -13.7 0 0.00138 0.0130 N=0, mu_star=0, lambda_star=0.00138 ## 10 0 -13.6 0 0.00138 0.0131 N=0, mu_star=0, lambda_star=0.00138 ## # … with 50,591 more rows
0からNまでの整数と、x軸の値x_vecの全ての組み合わせを作成します。
N+1回分の予測分布パラメータmu_star(のベクトル)を計算(作成)して、各試行に対応する値を抽出します。lambda_starについてはn列を使って計算できます。
確率変数とパラメータの組み合わせごとに確率を計算して、パラメータごとにフレーム切替用のラベルを作成します。
作図処理
事後分布と予測分布の推移をそれぞれアニメーションで可視化します。
・コード(クリックで展開)
観測データと対応するラベルをデータフレームに格納します。
# 観測データを格納 anime_data_df <- tibble::tibble( x = c(NA, x_n), param = unique(anime_posterior_df[["param"]]) # フレーム切替用ラベル ) anime_data_df
## # A tibble: 101 × 2 ## x param ## <dbl> <fct> ## 1 NA N=0, m=0, lambda_mu=0.0016 ## 2 9.83 N=1, m=8.47, lambda_mu=0.0116 ## 3 27.4 N=2, m=17.21, lambda_mu=0.0216 ## 4 22.5 N=3, m=18.89, lambda_mu=0.0316 ## 5 37.9 N=4, m=23.47, lambda_mu=0.0416 ## 6 17.5 N=5, m=22.31, lambda_mu=0.0516 ## 7 19.1 N=6, m=21.8, lambda_mu=0.0616 ## 8 17.7 N=7, m=21.22, lambda_mu=0.0716 ## 9 20.1 N=8, m=21.09, lambda_mu=0.0816 ## 10 37.5 N=9, m=22.89, lambda_mu=0.0916 ## # … with 91 more rows
事前分布には観測データが影響しないので欠損値NAとして、anime_posterior_dfのラベル列と合わせて格納します。
参考として、各試行における観測データを描画するのに使います。
事後分布の推移のアニメーションを作成します。
# μの事後分布のアニメーションを作図 posterior_graph <- ggplot() + geom_vline(mapping = aes(xintercept = mu_truth, color = "param"), size = 1, linetype = "dashed", show.legend = FALSE) + # 真のパラメータ geom_line(data = anime_posterior_df, mapping = aes(x = mu, y = dens, color = "posterior"), size = 1) + # μの事後分布 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("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)))) + # 凡例の体裁:(凡例表示用) labs(title = "Gaussian Distribution", subtitle = "{current_frame}", x = expression(mu), 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_predict_dfのラベル列を使って観測データを格納します。
# 観測データを格納 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_star=0, lambda_star=0.00138 ## 2 9.83 N=1, mu_star=8.47, lambda_star=0.00537 ## 3 27.4 N=2, mu_star=17.21, lambda_star=0.00684 ## 4 22.5 N=3, mu_star=18.89, lambda_star=0.0076 ## 5 37.9 N=4, mu_star=23.47, lambda_star=0.00806 ## 6 17.5 N=5, mu_star=22.31, lambda_star=0.00838 ## 7 19.1 N=6, mu_star=21.8, lambda_star=0.0086 ## 8 17.7 N=7, mu_star=21.22, lambda_star=0.00877 ## 9 20.1 N=8, mu_star=21.09, lambda_star=0.00891 ## 10 37.5 N=9, mu_star=22.89, lambda_star=0.00902 ## # … with 91 more rows
各試行における観測データを描画するのに使います。
また、各試行までの観測データを格納します。
# 過去の観測データを複製 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 9.83 N=2, mu_star=17.21, lambda_star=0.00684 ## 2 3 1 9.83 N=3, mu_star=18.89, lambda_star=0.0076 ## 3 3 2 27.4 N=3, mu_star=18.89, lambda_star=0.0076 ## 4 4 1 9.83 N=4, mu_star=23.47, lambda_star=0.00806 ## 5 4 2 27.4 N=4, mu_star=23.47, lambda_star=0.00806 ## 6 4 3 22.5 N=4, mu_star=23.47, lambda_star=0.00806 ## 7 5 1 9.83 N=5, mu_star=22.31, lambda_star=0.00838 ## 8 5 2 27.4 N=5, mu_star=22.31, lambda_star=0.00838 ## 9 5 3 22.5 N=5, mu_star=22.31, lambda_star=0.00838 ## 10 5 4 37.9 N=5, mu_star=22.31, lambda_star=0.00838 ## # … with 4,940 more rows
フレーム番号と試行回数(データ番号)の全ての組み合わせを作成して、フレーム番号未満のデータを抽出します。フレーム切替用のラベルは、フレーム番号を使って抽出します。
各試行以前のデータを描画するのに使うので、フレーム番号が2からになります。
予測分布の推移のアニメーションを作成します。
# 真の分布を計算 model_df <- tibble::tibble( x = x_vec, # 確率変数 dens = dnorm(x, mean = mu_truth, sd = 1/sqrt(lambda)) # 確率密度 ) # 予測分布のアニメーションを作図 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)))) + # 凡例の体裁:(凡例表示用) labs(title = "Gaussian 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.16：加筆修正の際に、数式読解編とRで実装編に記事を分割しました。
推論結果だけでなく何がどうなったかを明確にするためにそれなりに加筆しました。それによって自分の中ではかなり整理できたのですが、その結果文量が多くなってしまって読みにくくなった気がしなくもない。
- 2022.08.08：加筆修正しました。
【次節の内容】