はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。基本的な内容は「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」になります。「数式」と「プログラム」から理解するのが目標です。
この記事は、3.3.2項の内容です。「尤度関数を精度が未知の1次元ガウス分布(正規分布)」、「事前分布をガンマ分布」とした場合の「パラメータの事後分布」と「未観測値の予測分布」の計算をR言語で実装します。
省略してある内容等ありますので、本とあわせて読んでください。初学者な自分が理解できるレベルまで落として書き下していますので、分かる人にはかなりくどくなっています。同じような立場の人のお役に立てれば幸いです。
【数式読解編】
【他の節の内容】
【この節の内容】
・Rでやってみよう
人工的に生成したデータを用いて、ベイズ推論を行ってみましょう。
利用するパッケージを読み込みます。
# 3.3.2項で利用するパッケージ library(tidyverse)
・モデルの構築
まずは、モデルを設定します。
尤度(ガウス分布)のパラメータを設定します。
# 真のパラメータwを指定 mu <- 25 lambda_truth <- 0.01
平均パラメータ$\mu$をmuとして値を指定します。これは与えられている値として使います。
精度パラメータ$\lambda$をlambda_truthとして値を指定します。精度は分散の逆数なので、値が大きいほど散らばり具合が小さくなり、逆数の平方根が標準偏差$\sigma = \sqrt{\lambda^{-1}}$になります。これが真のパラメータであり、この値を求めるのがここでの目的です。
尤度の確率密度を計算して、作図用のデータフレームを作成します。
# 作図用のxの値を設定 x_line <- seq( mu - 4 * sqrt(1 / lambda_truth), mu + 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)^2) # 確率密度 #density = dnorm(x, mean = mu, 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=", round(mu, 2), ", lambda=",lambda_truth))
真のパラメータを求めることは、この真の分布を求めることを意味します。
・データの生成
続いて、構築したモデルに従って観測データ$\mathbf{X} = \{x_1, x_2, \cdots, x_N\}$を生成します。
ガウス分布に従う$N$個のデータをランダムに生成します。
# (観測)データ数を指定 N <- 50 # ガウス分布に従うデータを生成 x_n <- rnorm(n = N, mean = mu, sd = sqrt(1 / lambda_truth))
生成するデータ数$N$をNとして値を指定します。
ガウス分布に従う乱数は、rnorm()で生成できます。試行回数の引数nにN、平均の引数meanにmu、標準偏差の引数sdにsqrt(1 / lambda_truth)を指定します。生成したN個のデータをx_nとします。
観測したデータを確認しましょう。
# 確認 summary(x_n)
## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 5.29 16.70 26.99 24.84 30.37 50.08
観測データをヒストグラムでも確認します。
# 観測データのヒストグラムを作図 tibble(x = x_n) %>% ggplot(aes(x = x)) + geom_histogram(binwidth = 1) + # 観測データ labs(title = "Gaussian Distribution", subtitle = paste0("N=", N, ", mu=", mu, ", sigma=", round(sqrt(1 / lambda_truth), 1)))
データ数が十分に大きいと、分布の形状が真の分布に近づきます。
・事前分布の設定
次に、尤度に対する共役事前分布を設定します。
$\lambda$の事前分布(ガンマ分布)のパラメータ(超パラメータ)を設定します。
# lambdaの事前分布のパラメータを指定 a <- 1 b <- 1
ガンマ分布のパラメータ$a,\ b$をそれぞれa, bとして正の実数値を指定します。
$\lambda$の事前分布の確率密度を計算します。
# 作図用のlambdaの値を設定 lambda_line <- seq(0, 4 * lambda_truth, length.out = 1000) # lambdaの事前分布を計算:式(2.56) prior_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_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
log(0)が-Infになります。
$\lambda$の事前分布を作図します。
# lambdaの事前分布を作図 ggplot(prior_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}$からパラメータ$\lambda$の事後分布を求めます(パラメータ$\lambda$を分布推定します)。
観測データx_nを用いて、$\lambda$の事後分布(ガンマ分布)のパラメータを計算します。
# lambdaの事後分布のパラメータを計算:式(3.69) a_hat <- 0.5 * N + a b_hat <- 0.5 * sum((x_n - mu)^2) + b
$\lambda$の事後分布のパラメータは
で計算して、結果をa_hat, b_hatとします。
# 確認 a_hat; b_hat
## [1] 26 ## [1] 2275.808
求めたパラメータを使って、$\lambda$の事後分布の確率密度を計算します。
# lambdaの事後分布を計算:式(2.56) posterior_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_df)
## # A tibble: 6 x 3 ## lambda ln_C_Gam density ## <dbl> <dbl> <dbl> ## 1 0 143. 0. ## 2 0.0000400 143. 1.31e-48 ## 3 0.0000801 143. 4.02e-41 ## 4 0.000120 143. 9.26e-37 ## 5 0.000160 143. 1.12e-33 ## 6 0.000200 143. 2.71e-31
$\lambda$の事後分布を作図します。
# lambdaの事後分布を作図 ggplot(posterior_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_{*}$の予測分布を求めます。
$\lambda$の事後分布のパラメータa_hat, b_hat、または観測データx_nと$\lambda$の事前分布のパラメータa, bを用いて、予測分布(スチューデントのt分布)のパラメータを計算します。
# 予測分布のパラメータを計算:式(3.79) mu_s <- mu lambda_s_hat <- a_hat / b_hat nu_s_hat <- 2 * a_hat # 予測分布のパラメータを計算:式(3.79) lambda_s_hat <- (N + 2 * a) / (sum((x_n - mu)^2) + 2 * b) nu_s_hat <- N + 2 * a
予測分布のパラメータは
で計算して、結果をmu_s, lambda_s_hat, nu_s_hatとします。$\mu_s$は平均、$\lambda_s$は何？、$\nu_s$は自由度です。
それぞれ上の式だと、事後分布のパラメータa_hat, b_hatを使って計算できます。下の式だと、観測データx_nと事前分布のパラメータa, bを使って計算できます。
# 確認 mu_s; lambda_s_hat; nu_s_hat
## [1] 25 ## [1] 0.01142451 ## [1] 52
$\mathbf{X}$から$\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)^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.78 -7.98 0.0000145 ## 2 -14.9 1.62 -4.78 -7.96 0.0000149 ## 3 -14.8 1.62 -4.78 -7.93 0.0000153 ## 4 -14.8 1.62 -4.78 -7.90 0.0000157 ## 5 -14.7 1.62 -4.78 -7.87 0.0000162 ## 6 -14.6 1.62 -4.78 -7.85 0.0000166
予測分布を尤度と重ねて作図します。
# 予測分布を作図 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=", mu_s, ", lambda_s_hat=", round(lambda_s_hat, 3), ", nu_s_hat=", nu_s_hat))
観測データが増えると、予測分布が真の分布に近づきます。
・おまけ：アニメーションで推移の確認
gganimateパッケージを利用して、事後分布と予測分布の推移をアニメーション(gif画像)で確認するためのコードです。
・コード(クリックで展開)
異なる点のみを簡単に解説します。
# 利用するパッケージ library(tidyverse) library(gganimate)
・モデルの設定
# 真のパラメータを指定 mu <- 25 lambda_truth <- 0.01 # lambdaの事前分布のパラメータを指定 a <- 1 b <- 1 # 作図用のlambdaの値を設定 lambda_line <- seq(0, 4 * lambda_truth, length.out = 1000) # lambdaの事前分布(ガンマ分布)を計算:式(2.56) posterior_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.79) mu_s <- mu lambda_s <- a / b nu_s <- 2 * a # 作図用のxの値を設定 x_line <- seq( mu - 4 * sqrt(1 / lambda_truth), mu + 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, ", mu_s=", mu_s, ", lambda_s=", round(lambda_s, 3), ", nu_s", nu_s) ) # パラメータ )
各試行の結果を同じデータフレームに格納していく必要があります。事後分布をposterior_df、予測分布をpredict_dfとして、初期値の結果を持つように作成しておきます。
・推論処理
# データ数(試行回数)を指定 N <- 100 # 観測データの受け皿を初期化 x_n <- rep(0, N) # ベイズ推論 for(n in 1:N){ # ガウス分布に従うデータを生成 x_n[n] <- rnorm(n = 1, mean = mu, sd = sqrt(1 / lambda_truth)) # lambdaの事後分布のパラメータを更新:式(3.69) a <- 1 / 2 + a b <- 0.5 * (x_n[n] - mu)^2 + b # lambdaの事後分布(ガンマ分布)を計算:式(2.56) tmp_posterior_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.79) mu_s <- mu lambda_s <- a / 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, ", mu_s=", mu_s, ", lambda_s_hat=", round(lambda_s, 3), ", nu_s_hat=", nu_s) ) # パラメータ ) # 推論結果を結合 posterior_df <- rbind(posterior_df, tmp_posterior_df) predict_df <- rbind(predict_df, tmp_predict_df) }
観測された各データによってどのように学習する(分布が変化する)のかを確認するため、for()で1データずつ処理します。よって、データ数Nがイタレーション数になります。
超パラメータに関して、$\hat{a},\ \hat{b}$に対応するa_hat, b_hatを新たに作るのではなく、a, bをイタレーションごとに更新(上書き)していきます。
それに伴い、事後分布のパラメータの計算式(3.69)の$\sum_{n=1}^N (x_n - \mu)^2$と$\frac{N}{2}$の計算は、ループ処理によってN回繰り返し(x_n[n] - mu)^2と0.5を加えることで行います。n回目のループ処理のときには、n-1回分の(x_n[n] - mu)^2と0.5が既にaとbに加えられているわけです。
・事後分布の推移
# lambdaの事後分布を作図 posterior_graph <- ggplot(posterior_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_graph, nframes = N + 1, fps = 10)
各フレームの順番を示す列(label)をgganimate::transition_manual()に指定します。分布の推移と共にパラメータの値を表示するようにlabel列を作成していますが、ややこしければas.factor(paste0("iter=", n))として試行回数だけ表示するだけでもそれっぽくなります。
初期値(事前分布)を含むため、フレーム数の引数nframesはN + 1です。
・予測分布の推移
# 尤度を計算:式(2.64) model_df <- tibble( x = x_line, # x軸の値 density = dnorm(x, mean = mu, 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.17：加筆修正の際に、数式読解編とRで実装編に記事を分割しました。
ガウスなのかガンマなのか字面が似ていてこんがらがってくる。
【次節の内容】