はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」によって、「数式」と「プログラム」から理解するのが目標です。
省略してある内容等ありますので、本とあわせて読んでください。
この記事では、「尤度関数をベルヌーイ分布」、「事前分布をベータ分布」とした場合の「パラメータの事後分布」と「未観測値の予測分布」の計算をR言語で実装します。
【数式読解編】
【他の節の内容】
【この節の内容】
3.2.1 ベルヌーイ分布の学習と予測
尤度関数をベルヌーイ分布、事前分布をベータ分布とするモデルに対するベイズ推論を実装します。人工的に生成したデータを用いて、ベルヌーイ分布のパラメータを推定し、また未観測データに対する予測分布を求めます。ベルヌーイ分布については「【R】ベルヌーイ分布の作図 - からっぽのしょこ」、ベータ分布については「【R】ベータ分布の作図 - からっぽのしょこ」を参照してください。
利用するパッケージを読み込みます。
# 利用パッケージ library(tidyverse) library(gganimate)
この記事では、基本的にパッケージ名::関数名()の記法を使うので、パッケージを読み込む必要はありません。ただし、作図コードがごちゃごちゃしないようにパッケージ名を省略しているため、ggplot2は読み込む必要があります。また、基本的にベースパイプ(ネイティブパイプ)演算子|>を使いますが、パイプ演算子%>%でないと処理できない部分があるため、magrittrも読み込む必要があります。
分布の変化をアニメーション(gif画像)で確認するのにgganimateパッケージを利用します。不要であれば省略してください。
ベイズ推論の実装
まずは、モデルを設定して、データを生成します。生成したデータを用いて、事後分布のパラメータを計算します。さらに、事後分布のパラメータを用いて、予測分布のパラメータを計算します。
生成分布の設定
データ生成分布(真の分布)として、ベルヌーイ分布$\mathrm{Bern}(x | \mu)$を設定します。
真の分布(ベルヌーイ分布)のパラメータ$\mu$を設定します。
# 真のパラメータを指定 mu_truth <- 0.25
成功確率($x_n = 1$となる確率)$0 \leq \mu \leq 1$をmu_truthとして値を指定します。mu_truthが真のパラメータであり、この値を求めるのがここでの目的です。
真の分布を計算して、作図用のデータフレームを作成します。
# xがとり得る値を作成 x_vec <- 0:1 # 真の分布を計算:式(2.16) model_df <- tibble::tibble( x = x_vec, # x軸の値 prob = c(1 - mu_truth, mu_truth) # 確率 ) model_df
## # A tibble: 2 × 2 ## x prob ## <int> <dbl> ## 1 0 0.75 ## 2 1 0.25
ベルヌーイ分布に従う確率変数がとり得る値$x \in \{0, 1\}$を作成して、値ごとに確率を計算します。
ベルヌーイ分布に従う確率は、二項分布の確率計算関数dbinom()で計算できますが、ここでは簡単に、失敗確率($x_n = 0$となる確率)$1 - \mu$と$\mu$をデータフレームに格納します。
真の分布のグラフを作成します。
# 真の分布を作図 ggplot() + geom_bar(data = model_df, mapping = aes(x = x, y = prob, fill = "model"), stat = "identity") + # 真の分布 scale_fill_manual(breaks = "model", values = "purple", labels = "true model", name = "") + # バーの色:(凡例表示用) scale_x_continuous(breaks = x_vec, minor_breaks = FALSE) + # x軸目盛 ylim(c(0, 1)) + # y軸の表示範囲 labs(title = "Bernoulli Distribution", subtitle = parse(text = paste0("mu==", mu_truth)), x = "x", y = "probability")
真のパラメータを求めることは、真の分布を求めることを意味します。
データの生成
続いて、設定した生成分布から観測データ$\mathbf{X} = \{x_1, x_2, \cdots, x_N\}$を生成します。
ベルヌーイ分布に従う$N$個のデータをランダムに生成します。
# (観測)データ数を指定 N <- 100 # ベルヌーイ分布に従うデータを生成 x_n <- rbinom(n = N, size = 1, prob = mu_truth) head(x_n)
## [1] 0 1 0 0 0 0
生成するデータ数$N$を指定します。
ベルヌーイ分布に従う乱数は、二項分布に従う乱数生成関数rbinom()の試行回数の引数sizeを1にすることで生成できます。データ数の引数nにN、成功確率の引数probにmu_truthを指定します。生成したN個のデータをx_nとします。
観測データ$\mathbf{X}$の度数を集計します。
# 観測データを集計 freq_df <- tibble::tibble( x = x_vec, # x軸の値 freq = c(N - sum(x_n), sum(x_n)) # 度数 ) freq_df
## # A tibble: 2 × 2 ## x freq ## <int> <dbl> ## 1 0 75 ## 2 1 25
ベルヌーイ分布の乱数$x_n$は0か1の値をとるので、1の要素数はx_nの総和で得られます。また、データ数はNなので、0の要素数はNと1の要素数の差で得られます。
観測データのヒストグラムを真の分布と重ねて確認します。
# 観測データのヒストグラムを作成 ggplot() + geom_bar(data = freq_df, mapping = aes(x = x, y = freq/N, fill = "data"), stat = "identity") + # 観測データ(相対度数) geom_bar(data = model_df, mapping = aes(x = x, y = prob, fill = "model", color = "model"), stat = "identity", 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 = "") + # 線の色:(凡例表示用) scale_x_continuous(breaks = x_vec, minor_breaks = FALSE) + # x軸目盛 ylim(c(0, 1)) + # y軸の表示範囲 labs(title = "Bernoulli Distribution", subtitle = parse(text = paste0("list(mu==", mu_truth, ", N==", N, "(", paste0(freq_df[["freq"]], collapse = ", "), "))")), x = "x", y = "relative frequency")
各値の度数(freq列)をデータ数Nで割って、相対度数をバーの高さとします。
データ数が十分に大きいと、ヒストグラムの形状が真の分布に近付きます。
事前分布の設定
次は、尤度に対する共役事前分布を設定します。ベルヌーイ分布のパラメータ$\mu$の事前分布としてベータ分布$\mathrm{Beta}(\mu | a, b)$を設定します。
事前分布(ベータ分布)のパラメータ(超パラメータ)$a, b$を設定します。
# 事前分布のパラメータを指定 a <- 1 b <- 1
$a > 0, b > 0$の値を指定します。
事前分布を計算します。
# グラフ用のmuの値を作成 mu_vec <- seq(0, 1, length.out = 501) # 事前分布を計算:式(2.41) prior_df <- tibble::tibble( mu = mu_vec, # x軸の値 density = dbeta(x = mu_vec, shape1 = a, shape2 = b) # 確率密度 ) prior_df
## # A tibble: 501 × 2 ## mu density ## <dbl> <dbl> ## 1 0 1 ## 2 0.002 1 ## 3 0.004 1 ## 4 0.006 1 ## 5 0.008 1 ## 6 0.01 1 ## 7 0.012 1 ## 8 0.014 1 ## 9 0.016 1 ## 10 0.018 1 ## # … with 491 more rows
ベータ分布の確率変数がとり得る値$0 \leq \mu \leq 1$を作成して、値ごとに確率密度を計算します。
ベータ分布の確率密度は、dbeta()で計算できます。確率変数の引数xにmu_vec、パラメータの引数shape1, shape2にa, bを指定します。
事前分布のグラフを作成します。
# 事前分布を作図 ggplot() + geom_line(data = prior_df, mapping = aes(x = mu, y = density, color = "prior"), size = 1) + # 事前分布 geom_vline(mapping = aes(xintercept = mu_truth, color = "param"), size = 1, linetype = "dashed", show.legend = FALSE) + # 真のパラメータ 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.8, 0.8), linetype = c(2, 1)))) + # 凡例の体裁:(凡例表示用) labs(title = "Beta Distribution", subtitle = parse(text = paste0("list(a==", a, ", b==", b, ")")), x = expression(mu), y = "density")
真のパラメータの値を赤色の破線で示します。
無情報事前分布として、全ての確率密度が等しい分布を設定しました。
事後分布の計算
観測データ$\mathbf{X}$からパラメータ$\mu$の事後分布を求めます(パラメータ$\mu$を分布推定します)。
観測データ$\mathbf{X}$を用いて、事後分布(ベータ分布)のパラメータ$\hat{a}, \hat{b}$を計算します。
# 事後分布のパラメータを計算:式(3.15) a_hat <- sum(x_n) + a b_hat <- N - sum(x_n) + b a_hat; b_hat
## [1] 26 ## [1] 76
事後分布のパラメータは、次の式で計算できます。
事前分布のパラメータ$a, b$に、それぞれ$x_n = 1, x_n = 0$の数を加えています。
事後分布を計算します。
# 事後分布を計算:式(2.41) posterior_df <- tibble::tibble( mu = mu_vec, # x軸の値 dens = dbeta(x = mu, shape1 = a_hat, shape2 = b_hat) # 確率密度 ) posterior_df
## # A tibble: 501 × 2 ## mu dens ## <dbl> <dbl> ## 1 0 0 ## 2 0.002 7.07e-43 ## 3 0.004 2.04e-35 ## 4 0.006 4.43e-31 ## 5 0.008 5.07e-28 ## 6 0.01 1.15e-25 ## 7 0.012 9.45e-24 ## 8 0.014 3.83e-22 ## 9 0.016 9.26e-21 ## 10 0.018 1.51e-19 ## # … with 491 more rows
更新した超パラメータa_hat, b_hatを使って、事前分布のときと同様にして処理します。
事後分布のグラフを作成します。
# 事後分布を作図 ggplot() + geom_line(data = posterior_df, mapping = aes(x = mu, y = dens, color = "posterior"), size = 1) + # 事後分布 geom_vline(mapping = aes(xintercept = mu_truth, color = "param"), size = 1, linetype = "dashed", show.legend = FALSE) + # 真のパラメータ 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.8, 0.8), linetype = c(2, 1)))) + # 凡例の体裁:(凡例表示用) labs(title = "Beta Distribution", subtitle = parse(text = paste0("list(N==", N, ", hat(a)==", a_hat, ", hat(b)==", b_hat, ")")), x = expression(mu), y = "density")
$\mu$の真の値付近をピークとする分布になってます。
予測分布の計算
最後に、観測データ$\mathbf{X}$または(観測データから求めた)事後分布のパラメータ$\hat{a}, \hat{b}$から未観測のデータ$x_{*}$の予測分布を求めます。
事後分布のパラメータ$\hat{a}, \hat{b}$、または観測データ$\mathbf{X}$と事前分布のパラメータ$a, b$を用いて、予測分布(ベルヌーイ分布)のパラメータを計算します。
# 予測分布のパラメータを計算:式(3.19') mu_star_hat <- a_hat / (a_hat + b_hat) #mu_star_hat <- (sum(x_n) + a) / (N + a + b) mu_star_hat
## [1] 0.254902
予測分布のパラメータは、次の式で計算できます。
1つ目の式だと、事後分布のパラメータa_hat, b_hatを使って計算できます。2つ目の式だと、観測データx_nと事前分布のパラメータa, bを使って計算できます。
$\hat{\mu}_{*}$は、$x_{*} = 1$となる確率(成功確率)を表し、$x_n = 1$の数が多いほど値が1に近付きます。
予測分布を計算します。
# 予測分布を計算:式(2.16) predict_df <- tibble::tibble( x = x_vec, # x軸の値 prob = c(1 - mu_star_hat, mu_star_hat) # 確率 ) predict_df
## # A tibble: 2 × 2 ## x prob ## <int> <dbl> ## 1 0 0.745 ## 2 1 0.255
予測分布のパラメータmu_star_hatを使って、真の分布のときと同様にして処理します。
予測分布のグラフを真の分布と重ねて作成します。
# 予測分布を作図 ggplot() + geom_bar(data = predict_df, mapping = aes(x = x, y = prob, fill = "predict"), stat = "identity") + # 予測分布 geom_bar(data = model_df, mapping = aes(x = x, y = prob, fill = "color", color = "model"), stat = "identity", size = 1, linetype = "dashed") + # 真の分布 scale_fill_manual(values = c(model = NA, predict ="purple"), na.value = NA, labels = c(model = "true model", predict = "predict"), name = "") + # バーの色:(凡例表示用) scale_color_manual(values = c(model = "red", predict ="purple"), labels = c(model = "true model", predict = "predict"), name = "") + # 線の色:(凡例表示用) scale_x_continuous(breaks = x_vec, minor_breaks = FALSE) + # x軸目盛 ylim(c(0, 1)) + # y軸の表示範囲 labs(title = "Bernoulli Distribution", subtitle = parse(text = paste0("list(N==", N, ", hat(mu)[s]==", round(mu_star_hat, 2), ")")), x = "x", y = "probability")
観測データが増えると、予測分布が真の分布に近付きます。
以上で、ベルヌーイ分布のベイズ推論を実装できました。
学習推移の可視化
gganimateパッケージを利用して、パラメータの推定値(更新値)の推移(分布の変化)をアニメーション(gif画像)で確認します。
モデルの設定
真の分布のパラメータ$\mu$と、事前分布のパラメータ(の初期値)$a, b$を設定します。
# 真のパラメータを指定 mu_truth <- 0.25 # 事前分布のパラメータを指定 a <- 1 b <- 1 # データ数(試行回数)を指定 N <- 300 # グラフ用のmuの値を作成 mu_vec <- seq(0, 1, length.out = 501) # グラフ用のxの値を作成 x_vec <- 0:1
実装時と同じく、パラメータを指定して、確率変数がとり得る値(x軸の値)を作成します。グラフが粗い場合や処理が重い場合は、mu_vecの要素数(値の間隔)を調整してください。
事後分布と予測分布のパラメータの計算(更新)処理について、2通りの方法を載せます。目的に応じて使い分けてください。
推論処理：for関数による処理
1つ目の処理方法では、for()を使って、1データずつ生成してパラメータの更新を繰り返し処理します。こちらの方が、前ステップで求めた事後分布(のパラメータ)を次ステップの事前分布(のパラメータ)として用いる逐次学習をイメージしやすいです。
・コード(クリックで展開)
超パラメータの初期値を使って、事前分布を計算します。
# 事前分布(ベータ分布)を計算:式(2.41) anime_posterior_df <- tibble::tibble( mu = mu_vec, # x軸の値 dens = dbeta(x = mu, shape1 = a, shape2 = b), # 確率密度 param = paste0("N=", 0, "=(0, 0), a=", a, ", b=", b) |> as.factor() # フレーム切替用ラベル ) anime_posterior_df
## # A tibble: 501 × 3 ## mu dens param ## <dbl> <dbl> <fct> ## 1 0 1 N=0=(0, 0), a=1, b=1 ## 2 0.002 1 N=0=(0, 0), a=1, b=1 ## 3 0.004 1 N=0=(0, 0), a=1, b=1 ## 4 0.006 1 N=0=(0, 0), a=1, b=1 ## 5 0.008 1 N=0=(0, 0), a=1, b=1 ## 6 0.01 1 N=0=(0, 0), a=1, b=1 ## 7 0.012 1 N=0=(0, 0), a=1, b=1 ## 8 0.014 1 N=0=(0, 0), a=1, b=1 ## 9 0.016 1 N=0=(0, 0), a=1, b=1 ## 10 0.018 1 N=0=(0, 0), a=1, b=1 ## # … with 491 more rows
現在のパラメータを文字列結合し因子型に変換して、フレーム切替用のラベル列とします。
また、事前分布のパラメータを使って、予測分布のパラメータと予測分布を計算します。
# 初期値による予測分布のパラメータを計算:式(3.19) mu_star <- a / (a + b) # 初期値による予測分布(ベルヌーイ分布)を計算:式(2.16) anime_predict_df <- tibble::tibble( x = x_vec, # x軸の値 prob = c(1 - mu_star, mu_star), # 対応する確率 param = paste0("N=", 0, "=(0, 0), mu_star=", round(mu_star, 2)) |> as.factor() # フレーム切替用ラベル ) anime_predict_df
## # A tibble: 2 × 3 ## x prob param ## <int> <dbl> <fct> ## 1 0 0.5 N=0=(0, 0), mu_star=0.5 ## 2 1 0.5 N=0=(0, 0), mu_star=0.5
同様に、フレーム切替用のラベル列を作成します。
パラメータの更新処理をN回繰り返します。
# 観測データの受け皿を作成 x_n <- rep(NA, times = N) # ベイズ推論 for(n in 1:N){ # ベルヌーイ分布に従うデータを生成 x_n[n] <- rbinom(n = 1, size = 1, prob = mu_truth) # 観測データを集計 freq_vec = c(n - sum(x_n[1:n]), sum(x_n[1:n])) # 事後分布のパラメータを更新:式(3.15) a <- x_n[n] + a b <- 1 - x_n[n] + b # 事後分布(ベータ分布)を計算:式(2.41) tmp_posterior_df <- tibble::tibble( mu = mu_vec, # x軸の値 dens = dbeta(x = mu, shape1 = a, shape2 = b), # 確率密度 param = paste0("N=", n, "=(", paste0(freq_vec, collapse = ", "), "), a=", a, ", b=", b) |> as.factor() # フレーム切替用ラベル ) # 予測分布のパラメーターを更新:式(3.19) mu_star <- a / (a + b) # 予測分布(ベルヌーイ分布)を計算:式(2.16) tmp_predict_df <- tibble::tibble( x = x_vec, # x軸の値 prob = c(1 - mu_star, mu_star), # 確率 param = paste0("N=", n, "=(", paste0(freq_vec, collapse = ", "), ")", ", mu_star=", round(mu_star, 3)) |> as.factor() # フレーム切替用ラベ ) # n回目の結果を結合 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.15)の$\sum_{n=1}^N$の計算は、ループ処理によってN回繰り返しx_n[n]や1-x_n[n]を加えることで行います。n回目のループ処理のときには、n-1回分のx_n[n]や1-x_n[n]が既にaとbに加えられています。
更新したパラメータを使って事後分布と予測分布を計算して、それぞれ計算結果をanime_***_dfに結合していきます。
結果を確認します。
# 確認 head(x_n); anime_posterior_df; anime_predict_df
## [1] 0 0 0 0 1 0 ## # A tibble: 75,651 × 3 ## mu dens param ## <dbl> <dbl> <fct> ## 1 0 1 N=0=(0, 0), a=1, b=1 ## 2 0.002 1 N=0=(0, 0), a=1, b=1 ## 3 0.004 1 N=0=(0, 0), a=1, b=1 ## 4 0.006 1 N=0=(0, 0), a=1, b=1 ## 5 0.008 1 N=0=(0, 0), a=1, b=1 ## 6 0.01 1 N=0=(0, 0), a=1, b=1 ## 7 0.012 1 N=0=(0, 0), a=1, b=1 ## 8 0.014 1 N=0=(0, 0), a=1, b=1 ## 9 0.016 1 N=0=(0, 0), a=1, b=1 ## 10 0.018 1 N=0=(0, 0), a=1, b=1 ## # … with 75,641 more rows ## # A tibble: 302 × 3 ## x prob param ## <int> <dbl> <fct> ## 1 0 0.5 N=0=(0, 0), mu_star=0.5 ## 2 1 0.5 N=0=(0, 0), mu_star=0.5 ## 3 0 0.667 N=1=(1, 0), mu_star=0.333 ## 4 1 0.333 N=1=(1, 0), mu_star=0.333 ## 5 0 0.75 N=2=(2, 0), mu_star=0.25 ## 6 1 0.25 N=2=(2, 0), mu_star=0.25 ## 7 0 0.8 N=3=(3, 0), mu_star=0.2 ## 8 1 0.2 N=3=(3, 0), mu_star=0.2 ## 9 0 0.833 N=4=(4, 0), mu_star=0.167 ## 10 1 0.167 N=4=(4, 0), mu_star=0.167 ## # … with 292 more rows
それぞれ「mu_vec, x_vecの要素数」掛ける「N+1」行のデータフレームになります。行数が増えすぎないように注意してくださいしてください。
推論処理：tidyverseパッケージによる処理
2つ目の処理方法では、tidyverseパッケージの関数を使って、一度の処理でN+1回分の計算をします。こちらの方が、処理時間が短いです。
・コード(クリックで展開)
ベルヌーイ分布に従う$N$個のデータを生成します。
# ベルヌーイ分布に従うデータを生成 x_n <- rbinom(n = N, size = 1, prob = mu_truth) head(x_n)
## [1] 0 0 0 0 1 0
ラベル用のテキストとして、観測データを集計します。
# 試行ごとに度数を集計:(ラベル用) freq_vec <- stringr::str_c( 0:N - cumsum(c(0, x_n)), cumsum(c(0, x_n)), sep = ", " ) head(freq_vec)
## [1] "0, 0" "1, 0" "2, 0" "3, 0" "4, 0" "4, 1"
試行ごとに、それまでの試行における0と1の数を求めて、str_c()でそれぞれ結合します。
事前分布を含めたN+1回分の事後分布を計算します。
# 試行ごとに事後分布(ベータ分布)を計算 anime_posterior_df <- tidyr::expand_grid( n = 0:N, # 試行回数 mu = mu_vec # x軸の値 ) |> # 全ての組み合わせを作成 dplyr::mutate( a = c(a, (cumsum(x_n) + a))[n+1], b = c(b, (1:N - cumsum(x_n) + b))[n+1] ) |> # 事後分布のパラメータを計算:式(3.15) dplyr::mutate( dens = dbeta(x = mu, shape1 = a, shape2 = b), # 確率密度 param = paste0("N=", n, "=(", freq_vec[n+1], "), a=", a, ", b=", b) %>% factor(., levels = unique(.)) # フレーム切替用ラベル ) # 事後分布を計算:式(2.41) anime_posterior_df
## # A tibble: 75,651 × 6 ## n mu a b dens param ## <int> <dbl> <dbl> <dbl> <dbl> <fct> ## 1 0 0 1 1 1 N=0=(0, 0), a=1, b=1 ## 2 0 0.002 1 1 1 N=0=(0, 0), a=1, b=1 ## 3 0 0.004 1 1 1 N=0=(0, 0), a=1, b=1 ## 4 0 0.006 1 1 1 N=0=(0, 0), a=1, b=1 ## 5 0 0.008 1 1 1 N=0=(0, 0), a=1, b=1 ## 6 0 0.01 1 1 1 N=0=(0, 0), a=1, b=1 ## 7 0 0.012 1 1 1 N=0=(0, 0), a=1, b=1 ## 8 0 0.014 1 1 1 N=0=(0, 0), a=1, b=1 ## 9 0 0.016 1 1 1 N=0=(0, 0), a=1, b=1 ## 10 0 0.018 1 1 1 N=0=(0, 0), a=1, b=1 ## # … with 75,641 more rows
試行回数(観測データ数)を表す0からNまでのN+1個の整数と、x軸の値mu_vecの全ての組み合わせをexpand_grid()で作成します。これにより、mu_vecの各要素をN+1フレーム分に複製できます。
事前分布とN回分の事後分布のパラメータa, b(のベクトル)を計算(作成)して、n列の値をインデックスとして使って各試行に対応する値を抽出します。
確率変数とパラメータの組み合わせごとに確率密度を計算して、パラメータごとにフレーム切替用のラベルを作成します。
同様に、初期値を含めたN+1回分の予測分布を計算します。
# 試行ごとに予測分布(ベルヌーイ分布)を計算 anime_predict_df <- tidyr::expand_grid( n = 0:N, # 試行回数 x = x_vec # x軸の値 ) |> # 全ての組み合わせを作成 dplyr::mutate( mu_star = c(a / (a + b), (cumsum(x_n) + a) / (1:N + a + b))[n+1] ) |> # 予測分布のパラメータを計算算:式(3.19) dplyr::mutate( prob = dbinom(x = x, size = 1, prob = mu_star), # 確率 param = paste0("N=", n, "=(", freq_vec[n+1], "), mu_star=", round(mu_star, 3)) %>% factor(., levels = unique(.)) # フレーム切替用ラベル ) # 予測分布を計算:式(2.16) anime_predict_df
## # A tibble: 302 × 5 ## n x mu_star prob param ## <int> <int> <dbl> <dbl> <fct> ## 1 0 0 0.5 0.5 N=0=(0, 0), mu_star=0.5 ## 2 0 1 0.5 0.5 N=0=(0, 0), mu_star=0.5 ## 3 1 0 0.333 0.667 N=1=(1, 0), mu_star=0.333 ## 4 1 1 0.333 0.333 N=1=(1, 0), mu_star=0.333 ## 5 2 0 0.25 0.75 N=2=(2, 0), mu_star=0.25 ## 6 2 1 0.25 0.25 N=2=(2, 0), mu_star=0.25 ## 7 3 0 0.2 0.8 N=3=(3, 0), mu_star=0.2 ## 8 3 1 0.2 0.2 N=3=(3, 0), mu_star=0.2 ## 9 4 0 0.167 0.833 N=4=(4, 0), mu_star=0.167 ## 10 4 1 0.167 0.167 N=4=(4, 0), mu_star=0.167 ## # … with 292 more rows
0からNまでの整数と、x軸の値x_vecの全ての組み合わせを作成します。
N+1回分の予測分布パラメータmu_star(のベクトル)を計算(作成)して、各試行に対応する値を抽出します。
確率変数とパラメータの組み合わせごとに確率を計算して、パラメータごとにフレーム切替用のラベルを作成します。
作図処理
事後分布と予測分布の推移をそれぞれアニメーションで可視化します。
・コード(クリックで展開)
観測データと対応するラベルをデータフレームに格納します。
# 観測データを格納 anime_data_df <- tibble::tibble( x = c(NA, x_n), param = unique(anime_posterior_df[["param"]]) # フレーム切替用ラベル ) anime_data_df
## # A tibble: 151 × 2 ## x param ## <int> <fct> ## 1 NA N=0=(0, 0), a=1, b=1 ## 2 0 N=1=(1, 0), a=1, b=2 ## 3 0 N=2=(2, 0), a=1, b=3 ## 4 0 N=3=(3, 0), a=1, b=4 ## 5 0 N=4=(4, 0), a=1, b=5 ## 6 1 N=5=(4, 1), a=2, b=5 ## 7 0 N=6=(5, 1), a=2, b=6 ## 8 1 N=7=(5, 2), a=3, b=6 ## 9 0 N=8=(6, 2), a=3, b=7 ## 10 0 N=9=(7, 2), a=3, b=8 ## # … with 141 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, 5), linetype = c(2, 1, 0), shape = c(NA, NA, 19)))) + # 凡例の体裁:(凡例表示用) labs(title = "Beta Distribution", subtitle = "{current_frame}", x = expression(mu), y = "density") # gif画像を作成 gganimate::animate(posterior_graph, nframes = N+1, fps = 10, width = 800, height = 600)
フレームの順番を示す列をtransition_manual()に指定して、animate()でgif画像を作成します。事前分布(初期値)を含むため、フレーム数はN+1です。
真の分布を複製して、対応するラベルとデータフレームに格納します。
# 生成分布をフレーム分に複製 anime_model_df <- tibble::tibble( x = rep(x_vec, times = N+1), # x軸の値 prob = rep(c(1 - mu_truth, mu_truth), times = N+1), # 確率 param = anime_predict_df[["param"]] # フレーム切替用ラベル ) anime_model_df
## # A tibble: 302 × 3 ## x prob param ## <int> <dbl> <fct> ## 1 0 0.75 N=0=(0, 0), mu_star=0.5 ## 2 1 0.25 N=0=(0, 0), mu_star=0.5 ## 3 0 0.75 N=1=(1, 0), mu_star=0.333 ## 4 1 0.25 N=1=(1, 0), mu_star=0.333 ## 5 0 0.75 N=2=(2, 0), mu_star=0.25 ## 6 1 0.25 N=2=(2, 0), mu_star=0.25 ## 7 0 0.75 N=3=(3, 0), mu_star=0.2 ## 8 1 0.25 N=3=(3, 0), mu_star=0.2 ## 9 0 0.75 N=4=(4, 0), mu_star=0.167 ## 10 1 0.25 N=4=(4, 0), mu_star=0.167 ## # … with 292 more rows
棒グラフの場合は、フレームごとにデータを用意する必要があります。
真の分布は全てのフレームで変化しないので、N+1個に複製します。
観測データと対応するラベルをデータフレームに格納します。
# 観測データを格納 anime_data_df <- tibble::tibble( x = c(NA, x_n), param = unique(anime_predict_df[["param"]]) # フレーム切替用ラベル ) anime_data_df
## # A tibble: 151 × 2 ## x param ## <int> <fct> ## 1 NA N=0=(0, 0), mu_star=0.5 ## 2 0 N=1=(1, 0), mu_star=0.333 ## 3 0 N=2=(2, 0), mu_star=0.25 ## 4 0 N=3=(3, 0), mu_star=0.2 ## 5 0 N=4=(4, 0), mu_star=0.167 ## 6 1 N=5=(4, 1), mu_star=0.286 ## 7 0 N=6=(5, 1), mu_star=0.25 ## 8 1 N=7=(5, 2), mu_star=0.333 ## 9 0 N=8=(6, 2), mu_star=0.3 ## 10 0 N=9=(7, 2), mu_star=0.273 ## # … with 141 more rows
こちらは、anime_predict_dfのラベル列を使います。
予測分布の推移のアニメーションを作成します。
# 予測分布のアニメーションを作図 predict_graph <- ggplot() + geom_bar(data = anime_predict_df, mapping = aes(x = x, y = prob, fill = "predict"), stat = "identity") + # 予測分布 geom_bar(data = anime_model_df, aes(x = x, y = prob, fill = "model", color = "model"), stat = "identity", size = 1, linetype = "dashed") + # 真の分布 geom_point(data = anime_data_df, aes(x = x, y = 0, color = "data"), size = 6) + # 観測データ gganimate::transition_manual(param) + # フレーム scale_fill_manual(values = c(model = NA, predict = "purple", data = "pink"), na.value = NA, labels = c(model = "true model", predict = "predict", data = "observation data"), name = "") + # バーの色:(凡例表示用) scale_color_manual(values = c(model = "red", predict = "purple", data = "pink"), labels = c(model = "true model", predict = "predict", data = "observation data"), name = "") + # 線の色:(凡例表示用) guides(fill = guide_legend(override.aes = list(fill = c(NA, "purple", NA))), color = guide_legend(override.aes = list(size = c(0.5, 0.5, 5), linetype = c(2, 0, 0), shape = c(NA, NA, 19)))) + # 凡例の体裁:(凡例表示用) scale_x_continuous(breaks = x_vec, minor_breaks = FALSE) + # x軸目盛 ylim(c(0, 1)) + # y軸の表示範囲 labs(title = "Bernoulli Distribution", subtitle = "{current_frame}", x = "x", y = "probability") # gif画像を作成 gganimate::animate(predict_graph, nframes = N+1, fps = 10, width = 600, height = 600)
データが増えるにつれて、真のパラメータ付近の確率密度が大きくなっていくのを確認できます。
予測分布が真の分布に近付いていくのを確認できます。
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
大幅加筆修正の際に記事を分割しました。
- 2022.07.09：加筆修正しました。
またまた書き直していきます。今回はグラフに載せる情報を増やすのが目的ですが、尤度と生成分布を勘違いしていたことに気付いたのでそれを早くどうにかしたいです。グラフを分かりやすくするために、作図コードが分かりにくくなってしまいました。その対応として、確率分布の計算・作図・乱数生成の記事を充実させたのでそちらも読んでみてください。
さすがにこのシリーズの修正はこれで最後にしたい。
【次節の内容】