はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。基本的な内容は「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」になります。「数式」と「プログラム」から理解するのが目標です。
この記事は、3.2.3項の内容です。尤度関数をポアソン分布、事前分布をガンマ分布とした場合のパラメータの事後分布と未観測値の予測分布の計算をR言語で実装します。
省略してある内容等ありますので、本とあせて読んでください。初学者な自分が理解できるレベルまで落として書き下していますので、分かる人にはかなりくどくなっています。同じような立場の人のお役に立てれば幸いです。
【数式読解編】
【他の節の内容】
【この節の内容】
・Rでやってみよう
人工的に生成したデータを用いて、ベイズ推論を行ってみましょう。
利用するパッケージを読み込みます。
# 3.2.3項で利用するパッケージ library(tidyverse)
・モデルの構築
まずは、モデルの設定を行います。
尤度(ポアソン分布)$p(\mathbf{X} | \lambda)$のパラメータ$\lambda$を設定します。
# 真のパラメータを指定 lambda_truth <- 4
$\lambda$をlambda_truthとして指定します。これが真のパラメータであり、この値を求めるのがここでの目的です。
尤度の確率を計算して、作図用のデータフレームを作成します。
# 作図用のxの値を設定 x_line <- seq(0, 4 * lambda_truth) # 尤度を計算:式(2.37) model_df <- tibble( x = x_line, # x軸の値 ln_C_poi = x * log(lambda_truth) - lgamma(x + 1), # 正規化項(対数) prob = exp(ln_C_poi - lambda_truth) # 確率 #prob = dpois(x = x, lambda = lambda_truth) # 確率 )
作図用に、ポアソン分布に従う変数$x_n$がとり得る非負の整数をseq()で作成します。この例では、0からlambda_truthの4倍を範囲とします。
x_lineの各要素の値となる確率は、ポアソン分布の定義式
や、ポアソン分布の確率計算関数dpois()で計算できます。
作成したデータフレームを確認しましょう。
# 確認 head(model_df)
## # A tibble: 6 x 3 ## x ln_C_poi prob ## <int> <dbl> <dbl> ## 1 0 0 0.0183 ## 2 1 1.39 0.0733 ## 3 2 2.08 0.147 ## 4 3 2.37 0.195 ## 5 4 2.37 0.195 ## 6 5 2.14 0.156
ggplot2パッケージを利用して作図するには、データフレームを渡す必要があります。
尤度を作図します。
# 尤度を作図 ggplot(model_df, aes(x = x, y = prob)) + geom_bar(stat = "identity", position = "dodge", fill = "purple") + # 尤度 labs(title = "Poisson Distribution", subtitle = paste0("lambda=", lambda_truth))
真のパラメータを求めることは、この真の分布を求めることを意味します。
・データの生成
続いて、構築した尤度に従って観測データ$\mathbf{X} = \{x_1, x_2, \cdots, x_N\}$を生成します。
ポアソン分布に従う$N$個のデータをランダムに生成します。
# (観測)データ数を指定 N <- 50 # ポアソン分布に従うデータを生成 x_n <- rpois(n = N ,lambda = lambda_truth)
生成するデータ数$N$をNとして、値を指定します。
ポアソン分布に従う乱数は、rpois()で生成できます。試行回数の引数nにN、パラメータの引数lambdaにlambda_truthを指定します。生成したN個のデータをx_nとします。
観測したデータ$\mathbf{X}$を確認しましょう。
# 確認 table(x_n)
## x_n ## 0 1 2 3 4 5 6 7 8 10 ## 1 1 10 9 13 4 6 3 2 1
$x_n$は、非負の整数になっています。
$\mathbf{X}$をヒストグラムでも確認します。
# 観測データのヒストグラムを作図 tibble(x = x_n) %>% ggplot(aes(x = x)) + geom_histogram(binwidth = 1) + # 観測データ labs(title = "Observation Data", subtitle = paste0("N=", N, ", lambda=", lambda_truth))
データ数が十分に大きいと、分布の形状が真の分布に近づきます。
・事前分布の設定
尤度に対する共役事前分布を設定します。
事前分布(ガンマ分布)のパラメータ(超パラメータ)を設定します。
# 事前分布のパラメータを指定 a <- 1 b <- 1
ガンマ分布のパラメータ$a,\ b$をそれぞれa, bとして、正の実数値を指定します。
事前分布の確率密度を計算します。
# 作図用のlambdaの値を設定 lambda_line <- seq(0, 2 * lambda_truth, by = 0.001) # 事前分布を計算:式(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$がとり得る0以上の値をseq()で用意します。by引数で間隔を指定できるので、グラフが粗かったり処理が重かったりする場合はこの値を調整してください。この例では、lambda_truthの2倍までとします。
lambda_lineの各要素の値に対して、確率密度を計算します。ガンマ分布の確率密度は、定義式
で計算します。ここで、$\Gamma(\cdot)$はガンマ関数であり、$\Gamma(x + 1) = x!$です。
ガンマ関数の計算は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.001 0 0.999 ## 3 0.002 0 0.998 ## 4 0.003 0 0.997 ## 5 0.004 0 0.996 ## 6 0.005 0 0.995
事前分布を作図します。
# 事前分布を作図 ggplot(prior_df, aes(x = lambda, y = density)) + geom_line(color = "purple") + # 事前分布 labs(title = "Gamma Distribution", subtitle = paste0("a=", a, ", b=", b), x = expression(lambda))
a, bの値を変更することで、ポアソン分布におけるパラメータと形状の関係を確認できます。
・事後分布の計算
観測データ$\mathbf{X}$からパラメータ$\lambda$の事後分布を求めます(パラメータ$\lambda$を分布推定します)。
観測データx_nを用いて事後分布のパラメータを計算します。
# 事後分布のパラメータを計算:式(3.38) a_hat <- sum(x_n) + a b_hat <- N + b
事後分布のパラメータは
で計算して、結果をa_hat, b_hatとします。
# 確認 a_hat; b_hat
## [1] 204 ## [1] 51
事前分布のパラメータ$a,\ b$に、それぞれ観測データの総和$\sum_{n=1}^N x_n$とデータ数$N$を加えています。
事後分布(ガンマ分布)の確率密度を計算します。
# 事後分布を計算:式(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 -77.1 0 ## 2 0.001 -77.1 0 ## 3 0.002 -77.1 0 ## 4 0.003 -77.1 0 ## 5 0.004 -77.1 0 ## 6 0.005 -77.1 0
事後分布を作図します。
# 事後分布を作図 ggplot(posterior_df, aes(x = lambda, y = density)) + geom_line(color = "purple") + # 事後分布 geom_vline(aes(xintercept = lambda_truth), color = "red", linetype = "dashed") + # 真のパラメータ labs(title = "Gamma Distribution", subtitle = paste0("N=", N, ", a_hat=", a_hat, ", b_hat=", b_hat), x = expression(lambda)) # ラベル
パラメータ$\lambda$の真の値付近をピークとする分布を推定できています。
・予測分布の計算
最後に、$\mathbf{X}$から未観測のデータ$x_{*}$の予測分布を求めます。
事後分布のパラメータa_hat, b_hat、または観測データx_nと事前分布のパラメータa, bを用いて予測分布(負の二項分布分布)のパラメータを計算します。
# 予測分布のパラメータを計算:式(3.44') r_hat <- a_hat p_hat <- 1 / (b_hat + 1) # 予測分布のパラメータを計算:式(3.44') r_hat <- sum(x_n) + a p_hat <- 1 / (N + 1 + b)
予測分布のパラメータは
で計算して、結果をそれぞれr_hat, p_hatとします。
上の式だと、事後分布のパラメータa_hat, b_hatを使って計算できます。下の式だと、観測データx_nと事前分布のパラメータa, bを使って計算できます。
# 確認 r_hat; p_hat
## [1] 204 ## [1] 0.01923077
$\hat{r},\ \hat{p}$は、$\mathbf{X}$から学習しているのが式からも分かります。
予測分布を計算します。
# 予測分布を計算:式(3.43) predict_df <- tibble( x = x_line, # x軸の値 ln_C_NB = lgamma(x + r_hat) - lgamma(x + 1) - lgamma(r_hat), # 正規化項(対数) prob = exp(ln_C_NB + r_hat * log(1 - p_hat) + x * log(p_hat)) # 確率 #prob = dnbinom(x = x, size = r_hat, prob = 1 - p_hat) # 確率 )
x_lineの各様の値に対して、確率を計算します。負の二項分布の確率は、定義式
または負の二項分布の確率計算関数dnbinom()で計算します。
計算結果は次のようになります。
# 確認 head(predict_df)
## # A tibble: 6 x 3 ## x ln_C_NB prob ## <int> <dbl> <dbl> ## 1 0 0 0.0190 ## 2 1 5.32 0.0747 ## 3 2 9.95 0.147 ## 4 3 14.2 0.194 ## 5 4 18.1 0.193 ## 6 5 21.9 0.155
予測分布を尤度と重ねて作図します。
# 予測分布を作図 ggplot() + geom_bar(data = predict_df, aes(x = x, y = prob), stat = "identity", position = "dodge", fill = "purple") + # 予測分布 geom_bar(data = model_df, aes(x = x, y = prob), stat = "identity", position = "dodge", alpha = 0, color = "red", linetype = "dashed") + # 真の分布 labs(title = "Negative Binomial Distribution", subtitle = paste0("N=", N, ", r_hat=", r_hat, ", p_hat=", round(p_hat, 3)))
観測データが増えると、予測分布が真の分布に近づきます。
・おまけ：推移の確認
gganimateパッケージを利用して、パラメータの推定値の推移のアニメーション(gif画像)を作成するためのコードです。
・コード(クリックで展開)
異なる点のみを簡単に解説します。
# 利用するパッケージ library(tidyverse) library(gganimate)
・モデルの設定
# 真のパラメータを指定 lambda_truth <- 4 # 事前分布のパラメータを指定 a <- 1 b <- 1 # 作図用のlambdaの値を設定 lambda_line <- seq(0, 2 * lambda_truth, by = 0.001) # 事前分布(ガンマ分布)を計算 posterior_df <- tibble( lambda = lambda_line, # x軸の値 density = dgamma(x = lambda, shape = a, rate = b), # 確率密度 label = as.factor(paste0("N=", 0, ", a=", a, ", b=", b)) # パラメータ ) # 初期値による予測分布のパラメータを計算 r <- a p <- 1 / (b + 1) # 作図用のxの値を設定 x_line <- seq(0, 4 * lambda_truth) # 初期値による予測分布(負の二項分布)を計算 predict_df <- tibble( x = x_line, # 作図用の値 prob = dnbinom(x = x, size = r, prob = 1 - p), # 確率 label = as.factor(paste0("N=", 0, ", r=", r, ", p=", round(p, 3))) # パラメータ )
試行ごとの結果を同じデータフレームに格納していく必要があります。事後分布をposterior_df、予測分布をpredict_dfとして、初期値の結果を持つように作成しておきます。
・推論処理
# データ数(試行回数)を指定 N <- 100 # 観測データの受け皿を初期化 x_n <- rep(0, N) # ベイズ推論 for(n in 1:N){ # ポアソン分布に従うデータを生成 x_n[n] <- rpois(n = 1 ,lambda = lambda_truth) # 事後分布のパラメータを更新:式(3.38) a <- sum(x_n[n] * 1) + a b <- 1 + b # 事後分布(ガンマ分布)を計算:式(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=", b)) # パラメータ ) # 予測分布のパラメータを更新:式(3.44) r <- a p <- 1 / (b + 1) # 予測分布(負の二項分布)を計算:式(3.43) tmp_predict_df <- tibble( x = x_line, # x軸の値 prob = dnbinom(x = x, size = r, prob = 1 - p), # 確率 label = as.factor(paste0("N=", n, ", r_hat=", r, ", p_hat=", round(p, 3))) # パラメータ ) # 推論結果を結合 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.38)の$\sum_{n=1}^N x_n$と$N$の計算は、ループ処理によってN回繰り返しx_n[n]と1を加えることで行います。n回目のループ処理のときには、n-1回分のx_n[n](bの場合は1)が既にaとbに加えられているわけです。
・事後分布の推移
# 事後分布を作図 posterior_graph <- ggplot(posterior_df, aes(x = lambda, y = density)) + geom_line(color = "purple") + # 事後分布 geom_vline(aes(xintercept = lambda_truth), color = "red", linetype = "dashed") + # 真のパラメータ gganimate::transition_manual(label) + # フレーム labs(title = "Gamma Distribution", subtitle = "{current_frame}", x = expression(lambda)) # gif画像を出力 gganimate::animate(posterior_graph, nframes = N + 1, fps = 10)
各フレームの順番を示す列をgganimate::transition_manual()に指定します。初期値(事前分布)を含むため、フレーム数はN + 1です。
・予測分布の推移
# 真のモデルをN+1フレーム分に複製 model_df <- tibble( x = rep(x_line, times = N + 1), # x軸の値 prob = rep(dpois(x = x_line, lambda = lambda_truth), times = N + 1), # 確率 label = predict_df[["label"]] # パラメータ ) # 予測分布を作図 predict_graph <- ggplot() + geom_bar(data = predict_df, aes(x = x, y = prob), stat = "identity", position = "dodge", fill = "purple") + # 予測分布 geom_bar(data = model_df, aes(x = x, y = prob), stat = "identity", position = "dodge", alpha = 0, color = "red", linetype = "dashed") + # 真のモデル gganimate::transition_manual(label) + # フレーム labs(title = "Negative Binomial Distribution", subtitle = "{current_frame}") # gif画像を出力 gganimate::animate(predict_graph, nframes = N + 1, fps = 10)
真の分布についても予測分布と同じフレーム数分用意する必要があります(たぶん？)。そのため、ラベルを対応させて値を複製しています。
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
加筆修正の際に記事を分割しました。もう少し修正のペースを上げたい。
【次節の内容】