はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。基本的な内容は「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」になります。「数式」と「プログラム」から理解するのが目標です。
この記事は、3.2.1項の内容です。尤度関数をベルヌーイ分布、事前分布をベータ分布とした場合のパラメータの事後分布と未観測値の予測分布の計算をR言語で実装します。
省略してある内容等ありますので、本とあせて読んでください。初学者な自分が理解できるレベルまで落として書き下していますので、分かる人にはかなりくどくなっています。同じような立場の人のお役に立てれば幸いです。
【数式読解編】
【他の節の内容】
【この節の内容】
3.2.1 ベルヌーイ分布の学習と予測
人工的に生成したデータを用いて、ベイズ推論を行ってみましょう。
利用するパッケージを読み込みます。
# 3.2.1項で利用するパッケージ library(tidyverse)
・モデルの構築
まずは、モデルの設定を行います。
尤度(ベルヌーイ分布)$p(x_n | \mu)$のパラメータ$\mu$を設定します。
# 真のパラメータを指定 mu_truth <- 0.25
$x_n = 1$となる確率$\mu$をmu_truthとして、$0 \leq \mu \leq 1$の値を指定します。これが真のパラメータであり、この値を求めるのがここでの目的です。
尤度の確率を計算して、作図用のデータフレームを作成します。
# 尤度(ベルヌーイ分布)を計算 model_df <- tibble( x = c(0, 1), # xがとり得る値 prob = c(1 - mu_truth, mu_truth) # 対応する確率 #prob = mu_truth^x * (1 - mu_truth)^(1 - x) # 確率 #prob = dbinom(x = x, size = 1, prob = mu_truth) # 確率 )
ここでは簡単に、$x_n = 0$となる確率$1 - \mu$と合わせて、データフレームに格納します。
ベルヌーイ分布の定義式
や、二項分布の確率計算関数dbinom()を使っても計算できます。ただし$a^0 = 1$なので、$x = 0$のとき$\mu^0 (1 - \mu)^1 = 1 - \mu$、$x = 1$のとき$\mu^1 (1 - \mu)^0 = \mu$になります。よって簡単に$1 - \mu$だけを計算しています。
作成したデータフレームを確認しましょう。
# 確認 head(model_df)
## # A tibble: 2 x 2 ## x prob ## <dbl> <dbl> ## 1 0 0.75 ## 2 1 0.25
ggplot2パッケージを利用して作図するには、データフレームを渡す必要があります。
尤度をプロットします。
# 尤度を作図 ggplot(model_df, aes(x = x, y = prob)) + geom_bar(stat = "identity", position = "dodge", fill = "purple") + # 棒グラフ ylim(c(0, 1)) + # y軸の表示範囲 labs(title = "Bernoulli Distribution", subtitle = paste0("mu=", mu_truth))
真のパラメータを求めることは、この真の分布を求めることを意味します。
・データの生成
続いて、構築した尤度に従って観測データ$\mathbf{X} = \{x_1, x_2, \cdots, x_N\}$を生成します。
ベルヌーイ分布に従う$N$個のデータをランダムに生成します。
# データ数を指定 N <- 50 # (観測)データを生成 x_n <- rbinom(n = N, size = 1, prob = mu_truth)
生成するデータ数$N$をNとして、値を指定します。
ベルヌーイ分布に従う乱数は、二項分布に従う乱数生成関数rbinom()のsize引数を1にすることで生成できます。また、試行回数の引数nにN、確率の引数probにmu_truthを指定します。生成したN個のデータをx_nとします。
観測したデータ$\mathbf{X}$を確認しましょう。
# 確認 table(x_n)
## x_n ## 0 1 ## 38 12
$x_n = 1$の数はsum(x_n)で、$x_n = 0$の数はN - sum(x_n)でも得られます。
$\mathbf{X}$をヒストグラムでも確認します。
# 観測データのヒストグラムを作図 tibble(x = x_n) %>% ggplot(aes(x = x)) + geom_histogram(bins = 2) + # ヒストグラム labs(title = "Observation Data", subtitle = paste0("N=", N, ", mu=", mu_truth))
データ数が十分に大きいと、分布の形状が真の分布に近づきます。
・事前分布の設定
尤度に対する共役事前分布を設定します。
事前分布(ベータ分布)のパラメータ(超パラメータ)を設定します。
# 事前分布のパラメータを指定 a <- 1 b <- 1
ベータ分布のパラメータ$a,\ b$をそれぞれa, bとして、$a > 0,\ b > 0$の値を指定します。
事前分布の確率密度を計算します。
# 事前分布(ベータ分布)を計算 prior_df <- tibble( mu = seq(0, 1, by = 0.001), # muがとり得る値 ln_C_beta = lgamma(a + b) - lgamma(a) - lgamma(b), # 正規化項(対数) density = exp(ln_C_beta) * mu^(a - 1) * (1 - mu)^(b - 1) # 確率密度 #density = dbeta(x = mu, shape1 = a, shape2 = b) # 確率密度 )
seq(0, 1)で、$\mu$がとり得る0から1までの値を用意します。by引数で間隔を指定できるので、グラフが粗かったり処理が重かったりする場合はこの値を調整してください。
mu列の各要素に対して、確率密度を計算します。ベータ分布の確率密度は、定義式
で計算します。ここで、$\Gamma(\cdot)$はガンマ関数です。
ガンマ関数の計算はgamma()で行えますが、値が大きくなると発散してしまします。そこで、対数をとったガンマ関数lgamma()で計算した後にexp()を使います。
ベータ分布の確率密度は、dbeta()でも計算できます。
計算結果は次のようになります。
# 確認 head(prior_df)
## # A tibble: 6 x 3 ## mu ln_C_beta density ## <dbl> <dbl> <dbl> ## 1 0 0 1 ## 2 0.001 0 1 ## 3 0.002 0 1 ## 4 0.003 0 1 ## 5 0.004 0 1 ## 6 0.005 0 1
事前分布を描画します。
# 事前分布を作図 ggplot(prior_df, aes(x = mu, y = density)) + geom_line(color = "purple") + # 折れ線グラフ labs(title = "Beta Distribution", subtitle = paste0("a=", a, ", b=", b), x = expression(mu))
a, bの値を変更することで、ベータ分布におけるパラメータと形状の関係を確認できます。
・事後分布の計算
観測データ$\mathbf{X}$からパラメータ$\mu$の事後分布を求めます(パラメータ$\mu$を分布推定します)。
観測データx_nを用いて事後分布のパラメータを計算します。
# 事後分布のパラメータを計算 a_hat <- sum(x_n) + a b_hat <- N - sum(x_n) + b
事後分布のパラメータは
で計算して、結果をa_hat, b_hatとします。
# 確認 a_hat; b_hat
## [1] 13 ## [1] 39
事前分布のパラメータ$a,\ b$に、それぞれ$x_n = 1,\ x_n = 0$の数を加えています。
事後分布(ベータ分布)の確率密度を計算します。
# 事後分布(ベータ分布)を計算 posterior_df <- tibble( mu = seq(0, 1, by = 0.001), # 作図用の値 ln_C_beta = lgamma(a_hat + b_hat) - lgamma(a_hat) - lgamma(b_hat), # 正規化項(対数) density = exp(ln_C_beta) * mu^(a_hat - 1) * (1 - mu)^(b_hat - 1) # 確率密度 #density = dbeta(x = mu, shape1 = a_hat, shape2 = b_hat) # 確率密度 )
更新した超パラメータa_hat, b_hatを用いて、事前分布のときと同様にして計算します。
計算結果は次のようになります。
# 確認 head(posterior_df)
## # A tibble: 6 x 3 ## mu ln_C_beta density ## <dbl> <dbl> <dbl> ## 1 0 29.5 0. ## 2 0.001 29.5 5.96e-24 ## 3 0.002 29.5 2.35e-20 ## 4 0.003 29.5 2.94e-18 ## 5 0.004 29.5 8.92e-17 ## 6 0.005 29.5 1.25e-15
事後分布を描画します。
# 事後分布を作図 ggplot(posterior_df, aes(x = mu, y = density)) + geom_line(color = "purple") + # 事後分布 geom_vline(aes(xintercept = mu_truth), color = "red", linetype = "dashed") + # 真のパラメータ labs(title = "Beta Distribution", subtitle = paste0("N=", N, ", a_hat=", a_hat, ", b_hat=", b_hat), x = expression(mu))
パラメータ$\mu$の真の値付近をピークとする分布を推定できています。
・予測分布の計算
最後に、$\mathbf{X}$から未観測のデータ$x_{*}$の予測分布を求めます。
事後分布のパラメータa_hat, b_hat、または観測データx_nと事前分布のパラメータa, bを用いて予測分布(ベルヌーイ分布)のパラメータを計算します。
# 予測分布のパラメータを計算 mu_hat_star <- a_hat / (a_hat + b_hat) mu_hat_star <- (sum(x_n) + a) / (N + a + b)
予測分布のパラメータの計算式
の結果をmu_hat_starとします。
上の式だと、事後分布のパラメータa_hat, b_hatを使って計算できます。下の式だと、観測データx_nと事前分布のパラメータa, bを使って計算できます。
# 確認
mu_hat_star
## [1] 0.25
$\hat{\mu}_{*}$は、$x_{*} = 1$となる確率を表し、$\mathbf{X}$から学習しているのが式からも分かります。
予測分布を計算します。
# 予測分布(ベルヌーイ分布)を計算 predict_df <- tibble( x = c(0, 1), # xがとり得る値 prob = c(1 - mu_hat_star, mu_hat_star) # 対応する確率 #prob = mu_hat_star^x * (1 - mu_hat_star)^(1 - x) # 確率 )
尤度のときと同様に処理できます。
計算結果は次のようになります。
# 確認 head(predict_df)
## # A tibble: 2 x 2 ## x prob ## <dbl> <dbl> ## 1 0 0.75 ## 2 1 0.25
予測分布を尤度と重ねて描画します。
# 予測分布を作図 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") + # 真の分布 ylim(c(0, 1)) + # y軸の表示範囲 labs(title = "Beta Distribution", subtitle = paste0("N=", N, ", mu_hat_star=", round(mu_hat_star, 2)))
観測データが増えると、予測分布が真の分布に近づきます(ここまでぴったりなのはたまたまです)。
おまけ：推移の確認
gganimateパッケージを利用して、パラメータの推定値の推移のアニメーション(gif画像)を作成するためのコードです。
・コード(クリックで展開)
# 利用するパッケージ library(tidyverse) library(gganimate)
異なる点のみを簡単に解説します。
### モデルの設定 # 真のパラメータを指定 mu_truth <- 0.25 # 事前分布のパラメータを指定 a <- 1 b <- 1 # 事前分布(ベータ分布)を計算 posterior_df <- tibble( mu = seq(0, 1, by = 0.001), # muがとり得る値 density = dbeta(x = mu, shape1 = a, shape2 = b), # 確率密度 label = as.factor(paste("N=", 0, ", a=", a, ", b", b)) # 試行回数とパラメータのラベル ) # 初期値による予測分布のパラメーターを計算 mu_star <- a / (a + b) # 初期値による予測分布(ベルヌーイ分布)を計算 predict_df <- tibble( x = c(0, 1), # xがとり得る値 prob = c(1 - mu_star, mu_star), # 対応する確率 label = as.factor(paste0("N=", 0, ", mu_hat_star=", round(mu_star, 2))) # 試行回数とパラメータのラベル )
試行ごとの結果を同じデータフレームに格納していく必要があります。事後分布をposterior_df、予測分布をpredict_dfとして、初期値の結果を持つように作成しておきます。
### 推論 # データ数(試行回数)を指定 N <- 100 # 推論処理 x_n <- rep(0, N) # 受け皿を初期化 for(n in 1:N){ # ベルヌーイ分布に従うデータを生成 x_n[n] <- rbinom(n = 1, size = 1, prob = mu_truth) # 事後分布のパラメータを更新 a <- x_n[n] + a b <- 1 - x_n[n] + b # 事後分布を計算 tmp_posterior_df <- tibble( mu = seq(0, 1, by = 0.001), # muがとり得る値 density = dbeta(x = mu, shape1 = a, shape2 = b), # 確率密度 label = as.factor(paste0("N=", n, ", a_hat=", a, ", b_hat=", b)) # 試行回数とパラメータのラベル ) # 予測分布のパラメーターを更新 mu_star <- a / (a + b) # 予測分布を計算 tmp_predict_df <- tibble( x = c(0, 1), # xがとり得る値 prob = c(1 - mu_star, mu_star), # 対応する確率 label = as.factor(paste0("N=", n, ", mu_hat_star=", round(mu_star, 2))) # 試行回数とパラメータのラベル ) # 推論結果を結合 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.15)の$\sum_{n=1}^N$の計算は、forループによってN回繰り返しx_n[n]を加えることで行います。n回目のループ処理のときには、n-1回分のx_n[n](bの場合は1)が既にaとbに加えられているわけです。
・事後分布の推移
### 作図 # 事後分布を作図 posterior_graph <- ggplot(posterior_df, aes(x = mu, y = density)) + geom_line(color = "purple") + # 事後分布 geom_vline(aes(xintercept = mu_truth), color = "red", linetype = "dashed") + # 真のパラメータ gganimate::transition_manual(label) + # フレーム labs(title = "Beta Distribution", subtitle = "{current_frame}", x = expression(mu)) # gif画像を出力 gganimate::animate(posterior_graph, nframes = (N + 1), fps = 10)
各フレームの順番を示す列をgganimate::transition_manual()に指定します。初期値(事前分布)を含むため、フレーム数はN + 1です。
・予測分布の推移
# Nフレーム分の真のモデル(ベルヌーイ分布)を計算 label_vec <- unique(predict_df[["label"]]) model_df <- tibble() for(n in 1:(N + 1)) { # n番目のフレーム用に計算 tmp_df <- tibble( x = c(0, 1), prob = c(1 - mu_truth, mu_truth), label = label_vec[n] ) # 結果を結合 model_df <- rbind(model_df, tmp_df) } # 予測分布を作図 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 = "Bernoulli Distribution", subtitle = "{current_frame}") # gif画像を出力 gganimate::animate(predict_graph, nframes = (N + 1), fps = 10)
真の分布についても予測分布と同じフレーム数分用意する必要があります(たぶん？)。そのため、ラベルを対応させて値を複製しています。
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
大幅加筆修正の際に記事を分割しました。
【次節の内容】