はじめに

　『ベイズ推論による機械学習入門』(MLSシリーズ)の独学時のノートです。各種のモデルやアルゴリズムについて「数式・プログラム・図」を用いて解説します。
　本の補助として読んでください。

　この記事では、ベルヌーイ分布に対するベイズ推論をR言語でスクラッチ実装します。

【前節の内容】

【他の節の内容】

【この節の内容】

3.2.1 ベルヌーイ分布のベイズ推論の実装

　ベルヌーイモデル(Bernoulli model)に対するベイズ推論(Bayesian inference)を実装して、人工的に生成したデータを用いて、パラメータの学習と未知変数の予測を行う。ベルヌーイモデルでは、尤度関数をベルヌーイ分布(Bernoulli distribution)、事前分布をベータ分布(Beta distribution)とする。この記事では、Rを利用して実装する。
　ベルヌーイモデルについては「3.2.1：ベルヌーイモデルの生成モデルの導出【緑ベイズ入門のノート】 - からっぽのしょこ」、ベイズ推論については「3.2.1：ベルヌーイ分布のベイズ推論の導出【緑ベイズ入門のノート】 - からっぽのしょこ」、Pythonを利用する場合は「【Python】3.2.1：ベルヌーイ分布のベイズ推論の実装【緑ベイズ入門のノート】 - からっぽのしょこ」を参照のこと。

　利用するパッケージを読み込む。

# 利用パッケージ
library(tidyverse)

　この記事では、基本的に パッケージ名::関数名() の記法を使うので、パッケージの読み込みは不要である。ただし、作図コードについては(ごちゃごちゃしないように)パッケージ名を省略するため、ggplot2 を読み込む必要がある。
　また、ネイティブパイプ演算子 |> を使う。(基本的には)magrittrパッケージのパイプ演算子 %>% に置き換えられるが、その場合は magrittr を読み込む必要がある。

ベイズ推論の実装

　まずは、ベルヌーイ分布に対するベイズ推論における一連の処理を確認する。
　生成モデル(ベルヌーイモデル)を設定して、モデルに従うデータ(トイデータ)を生成する。続いて、生成した観測データを用いて、事後分布の計算(パラメータの推定)を行う。さらに、事後分布のパラメータ(または観測データ)を用いて、予測分布の計算(未観測データの予測)を行う。

生成分布の設定

　データの生成分布(真の分布・ベルヌーイ分布) $p(x_n \mid \mu_{\mathrm{truth}}) = \mathrm{Bern}(x_n \mid \mu_{\mathrm{truth}})$ のパラメータ(真のパラメータ) $\mu_{\mathrm{truth}}$ を設定する。
　ベルヌーイ分布については「【R】ベルヌーイ分布の作図 - からっぽのしょこ」を参照のこと。

　生成分布のパラメータ $\mu_{\mathrm{truth}}$ を設定する。

# 真のパラメータを指定
mu_truth <- 0.25
mu_truth

[1] 0.25

　ベルヌーイ分布の成功確率パラメータ(０から1の値) $0 \leq \mu_{\mathrm{truth}} \leq 1$ を指定する。
　 $\mu_{\mathrm{truth}}$ がベルヌーイモデルにおける真のパラメータであり、この値を求めるのがここでの目的(学習)である。

　生成分布の確率変数 $x$ の作図範囲を設定する。

# x軸の範囲を設定
x_min <- 0 # (固定)
x_max <- 1 # (固定)

# x軸の値を作成
x_vec <- seq(from = x_min, to = x_max, by = 1)
x_min; x_max; head(x_vec)

[1] 0
[1] 1
[1] 0 1

　0か1の2値変数 $x \in \{0, 1\}$ がとり得る範囲を設定する。

　生成分布の確率を計算する。

# 生成分布の確率を計算:式(2.16)
model_df <- tibble::tibble(
  x    = x_vec, # 確率変数
  mu   = mu_truth, # 成功確率パラメータ
  prob = c(1-mu_truth, mu_truth) # 確率
)
model_df

# A tibble: 2 × 3
      x    mu  prob
  <dbl> <dbl> <dbl>
1     0  0.25  0.75
2     1  0.25  0.25

　 $x$ の値ごとに、ベルヌーイ分布に従う確率 $\mathrm{Bern}(x \mid \mu_{\mathrm{truth}})$ を計算する。
　この例では、失敗確率 $1-\mu$ と成功確率 $\mu$ をデータフレームに格納している。

　生成分布のグラフを作成する。

# 生成分布のラベルを作成
model_param_lbl <- paste0(
  "mu[truth] == ", round(mu_truth, digits = 2)
) |> 
  parse(text = _)

# 生成分布を作図
ggplot() + 
  geom_bar(
    data    = model_df, 
    mapping = aes(x = x, y = prob, fill = "model"), 
    stat = "identity", position = "identity"
  ) + # 真の分布
  scale_x_continuous(breaks = x_vec) + # x軸目盛
  scale_fill_manual(
    breaks = "model", 
    values = "red", 
    labels = "true model", 
    name   = ""
  ) + # (凡例表示用)
  theme(
    legend.title           = element_text(size = 0), 
    legend.position        = "inside", 
    legend.position.inside = c(1, 1), 
    legend.justification   = c(1, 1), 
    legend.background      = element_rect(fill = alpha("white", alpha = 0.8))
  ) + 
  labs(
    title = "Bernoulli distribution", 
    subtitle = model_param_lbl, 
    x = expression(x), 
    y = "probability"
  )

　真の分布(ベルヌーイ分布)を赤色の棒グラフで示す。

　真のパラメータ $\mu_{\mathrm{truth}}$ を求めることは、真の分布 $\mathrm{Bern}(x_n \mid \mu_{\mathrm{truth}})$ を求めることを意味する。

データの生成

　設定した生成分布(ベルヌーイ分布) $p(x_n \mid \mu_{\mathrm{truth}}) = \mathrm{Bern}(x_n \mid \mu_{\mathrm{truth}})$ に従うデータ(観測データ) $\mathbf{X}$ を作成する。
　ベルヌーイモデルのデータ生成については「【R】ベルヌーイ分布の乱数生成 - からっぽのしょこ」を参照のこと。

　生成分布からデータ $\mathbf{X}$ を生成する。

# データ数を指定
N <- 100

# 観測データを生成
x_n <- rbinom(n = N, size = 1, prob = mu_truth)
head(x_n)

[1] 1 0 1 0 0 1

　データ数 $N$ を指定して、ベルヌーイ分布に従う乱数 $\mathbf{X} = \{x_1, \cdots, x_N\}$ を生成する。
　二項分布の乱数生成関数 rbinom() のサンプルサイズの引数 n に $N$ 、試行回数の引数 size に 1、成功確率の引数 prob に $\mu_{\mathrm{truth}}$ を指定する。

　観測データ $\mathbf{X}$ を集計する。

# 観測データを集計
obs_df <- tibble::tibble(
  x        = x_vec,                   # 観測値
  freq     = c(N-sum(x_n), sum(x_n)), # 度数
  rel_freq = freq / N                 # 相対度数
)
obs_df

# A tibble: 2 × 3
      x  freq rel_freq
  <dbl> <dbl>    <dbl>
1     0    70      0.7
2     1    30      0.3

　 $x$ の値ごとに、観測データ x_n に含まれる要素数をカウントして、度数 $N_x$ と相対度数 $\frac{N_x}{N}$ を求める。
　この例では、失敗回数 $N_0 = N - N_1$ と成功回数 $N_1 = \sum_{n=1}^N x_n$ をデータフレームに格納している。

　観測データ $\mathbf{X}$ のグラフを作成する。

# 観測データのラベルを作成
obs_param_lbl <- paste0(
  "list(", 
  "N == ",         N, ", ", 
  "mu[truth] == ", round(mu_truth, digits = 2), 
  ")"
) |> 
  parse(text = _)

# 観測データを作図
ggplot() + 
  geom_bar(
    data    = model_df, 
    mapping = aes(x = x, y = prob, linetype = "model"), 
    stat = "identity", position = "identity",
    fill = NA, color = "red", linewidth = 1
  ) + # 生成分布
  geom_bar(
    data    = obs_df, 
    mapping = aes(x = x, y = rel_freq, linetype = "data"), 
    stat = "identity", position = "identity", 
    fill = "hotpink", color = NA, alpha = 0.5
  ) + # 観測データ
  scale_y_continuous(
    sec.axis = sec_axis(transform = ~ .*N, name = "frequency")
  ) + # 度数軸目盛
  scale_linetype_manual(
    breaks = c("model", "data"), 
    values = c("dashed", "blank"), 
    labels = c("true model\n(generation distribution)", "observation data"), 
    name = ""
  ) + # (凡例表示用)
  guides(
    linetype = guide_legend(override.aes = list(linewidth = 0.5))
  ) + 
  theme(
    legend.title           = element_text(size = 0), 
    legend.position        = "inside", 
    legend.position.inside = c(1, 1), 
    legend.justification   = c(1, 1), 
    legend.background      = element_rect(fill = alpha("white", alpha = 0.8))
  ) + 
  labs(
    title = "Bernoulli distribution", 
    subtitle = obs_param_lbl, 
    x = expression(x), 
    y = "probability, relative frequency"
  )

　 $x$ の値ごとの生成確率(生成分布)を赤色(白抜き)の棒グラフ(破線)、観測データ $\mathbf{X}$ の度数 $N_x$ (相対度数 $\frac{N_x}{N}$ )を桃色(塗りつぶし)の棒グラフ(実線)で示す。

　データ数 $N$ が十分に大きいと、観測データ $\mathbf{X}$ のヒストグラムの形状が生成分布 $\mathrm{Bern}(x_n \mid \mu_{\mathrm{truth}})$ に近付く。

事前分布の設定

　パラメータ $\mu$ の事前分布(ベータ分布) $p(\mu \mid a, b) = \mathrm{Beta}(\mu \mid a, b)$ のパラメータ(超パラメータ) $a, b$ を設定する。
　ベータ分布については「【R】ベータ分布の作図 - からっぽのしょこ」を参照のこと。

　事前分布のパラメータ $a, b$ を設定する。

# 事前分布のパラメータを指定
a <- 1
b <- 1
a; b

[1] 1
[1] 1

　ベータ分布の形状パラメータ(正の値) $a \gt 0, b \gt 0$ を指定する。

　事前分布の確率変数 $\mu$ の作図範囲を設定する。

# μ軸の範囲を設定
mu_min <- 0
mu_max <- 1

# μ軸の値を作成
mu_vec <- seq(from = mu_min, to = mu_max, length.out = 1001)
mu_min; mu_max; head(mu_vec)

[1] 0
[1] 1
[1] 0.000 0.001 0.002 0.003 0.004 0.005

　0から1の値の変数 $\mu \in (0, 1)$ がとり得る範囲を設定する。

　事前分布の確率密度を計算する。

# 事前分布の確率密度を計算:式(2.41)
prior_df <- tibble::tibble(
  mu   = mu_vec, # 確率変数
  a    = a, # 形状パラメータ
  b    = b, # 形状パラメータ
  dens = dbeta(x = mu, shape1 = a, shape2 = b) # 確率密度
)
prior_df

# A tibble: 1,001 × 4
      mu     a     b  dens
   <dbl> <dbl> <dbl> <dbl>
 1 0         1     1     1
 2 0.001     1     1     1
 3 0.002     1     1     1
 4 0.003     1     1     1
 5 0.004     1     1     1
 6 0.005     1     1     1
 7 0.006     1     1     1
 8 0.007     1     1     1
 9 0.008     1     1     1
10 0.009     1     1     1
# ℹ 991 more rows

　 $\mu$ の値ごとに、ベータ分布に従う確率密度 $\mathrm{Bern}(\mu \mid a, b)$ を計算する。
　ベータ分布の確率密度関数 dbeta() の確率変数の引数 x に $\mu$ 、形状パラメータの引数 shape1, shape2 に $a, b$ を指定する。

　事前分布のグラフを作成する。

# 事前分布のラベルを作成
prior_param_lbl <- paste0(
  "list(", 
  "mu[truth] == ", round(mu_truth, digits = 2), ", ", 
  "a == ",         round(a,        digits = 1), ", ", 
  "b == ",         round(b,        digits = 1),  
  ")"
) |> 
  parse(text = _)

# 事前分布を作図
ggplot() + 
  geom_vline(
    mapping = aes(xintercept = mu_truth, linetype = "model"), 
    color = "red", linewidth = 1
  ) + # 真のパラメータ
  geom_line(
    data    = prior_df, 
    mapping = aes(x = mu, y = dens, linetype = "prior"), 
    color = "purple", linewidth = 1
  ) + # 事前分布
  scale_x_continuous(
    sec.axis = sec_axis(
      transform = ~ ., 
      breaks    = mu_truth, 
      labels    = expression(mu[truth])
    ) # パラメータラベル
  ) + 
  scale_linetype_manual(
    breaks = c("model", "prior"), 
    values = c("dashed", "solid"), 
    labels = c("true parameter", "prior distribution"), 
    name = ""
  ) + # (凡例表示用)
  guides(
    linetype = guide_legend(override.aes = list(linewidth = 0.5))
  ) + 
  theme(
    legend.title           = element_text(size = 0), 
    legend.position        = "inside", 
    legend.position.inside = c(1, 1), 
    legend.justification   = c(1, 1), 
    legend.background      = element_rect(fill = alpha("white", alpha = 0.8))
  ) + 
  labs(
    title = "Beta distribution", 
    subtitle = prior_param_lbl, 
    x = expression(mu), 
    y = "density"
  )

　真のパラメータを赤色の垂直線(破線)、事前分布(ベータ分布)を紫色の曲線(実線)で示す。

　真のパラメータ(真の分布のパラメータ) $\mu_{\mathrm{truth}}$ と、パラメータ $\mu$ の事前分布 $\mathrm{Beta}(\mu \mid a, b)$ の位置関係を図で確認する。

　この例では、無情報事前分布として、全ての確率密度が等しい分布を設定している。

事後分布の計算

　観測データ $\mathbf{X}$ から、パラメータ $\mu$ の事後分布(ベータ分布) $p(\mu \mid \mathbf{X}, a, b) = \mathrm{Beta}(\mu \mid \hat{a}, \hat{b})$ のパラメータ(超パラメータ) $\hat{a}, \hat{b}$ を求める(真のパラメータ $\mu_{\mathrm{truth}}$ を分布推定する)。

　事後分布のパラメータ $\hat{a}, \hat{b}$ を計算する。

# 事後分布のパラメータを計算:式(3.15)
a_hat <- sum(x_n) + a
b_hat <- N - sum(x_n) + b
a_hat; b_hat

[1] 31
[1] 71

　事後分布のパラメータは、次の式で計算できる。

$\displaystyle \begin{aligned} \hat{a} &= \sum_{n=1}^N x_n + a \\ \hat{b} &= N - \sum_{n=1}^N x_n + b \end{aligned} \tag{3.15}$

　事後分布の確率密度を計算する。

# 事後分布の確率密度を計算:式(2.41)
posterior_df <- tibble::tibble(
  mu   = mu_vec, # 確率変数
  a    = a_hat, # 形状パラメータ
  b    = b_hat, # 形状パラメータ
  dens = dbeta(x = mu, shape1 = a, shape2 = b) # 確率密度
)
posterior_df

# A tibble: 1,001 × 4
      mu     a     b     dens
   <dbl> <dbl> <dbl>    <dbl>
 1 0        31    71 0       
 2 0.001    31    71 2.77e-63
 3 0.002    31    71 2.77e-54
 4 0.003    31    71 4.95e-49
 5 0.004    31    71 2.58e-45
 6 0.005    31    71 1.95e-42
 7 0.006    31    71 4.30e-40
 8 0.007    31    71 4.09e-38
 9 0.008    31    71 2.09e-36
10 0.009    31    71 6.68e-35
# ℹ 991 more rows

　事前分布のときと同様にして、ベータ分布に従う確率密度 $\mathrm{Beta}(\mu \mid \hat{a}, \hat{b})$ を計算する。

　事後分布のグラフを作成する。

# 事後分布のラベルを作成
posterior_param_lbl <- paste0(
  "list(", 
  "N == ",         N, ", ", 
  "mu[truth] == ", round(mu_truth, digits = 2), ", ", 
  "hat(a) == ",    round(a_hat,    digits = 1), ", ", 
  "hat(b) == ",    round(b_hat,    digits = 1), 
  ")"
) |> 
  parse(text = _)

# 事後分布を作図
ggplot() + 
  geom_vline(
    mapping = aes(xintercept = mu_truth, linetype = "model"), 
    color = "red", linewidth = 1
  ) + # 真のパラメータ
  geom_line(
    data    = posterior_df, 
    mapping = aes(x = mu, y = dens, linetype = "posterior"), 
    color = "purple", linewidth = 1
  ) + # 事後分布
  scale_x_continuous(
    sec.axis = sec_axis(
      transform = ~ ., 
      breaks    = mu_truth, 
      labels    = expression(mu[truth])
    ) # パラメータラベル
  ) + 
  scale_linetype_manual(
    breaks = c("model", "posterior"), 
    values = c("dashed", "solid"), 
    labels = c("true parameter", "posterior distribution"), 
    name = ""
  ) + # (凡例表示用)
  guides(
    linetype = guide_legend(override.aes = list(linewidth = 0.5)), 
  ) + 
  theme(
    legend.title           = element_text(size = 0), 
    legend.position        = "inside", 
    legend.position.inside = c(1, 1), 
    legend.justification   = c(1, 1), 
    legend.background      = element_rect(fill = alpha("white", alpha = 0.8))
  ) + 
  labs(
    title = "Beta distribution", 
    subtitle = posterior_param_lbl, 
    x = expression(mu), 
    y = "density"
  )

　真のパラメータを赤色の垂直線(破線)、事後分布(ベータ分布)を紫色の曲線(実線)で示す。

　データ数 $N$ が十分に大きいと、パラメータ $\mu$ の事後分布 $\mathrm{Beta}(\mu \mid \hat{a}, \hat{b})$ のピークが真のパラメータ $\mu_{\mathrm{truth}}$ に近付く。

予測分布の計算

　観測データ $\mathbf{X}$ から、未観測データ $x_{*}$ の予測分布(ベルヌーイ分布) $p(x_{*} \mid \mathbf{X}, a, b) = \mathrm{Bern}(x_{*} \mid \hat{\mu}_{*})$ のパラメータ $\hat{\mu}_{*}$ を求める。

　予測分布のパラメータ $\hat{\mu}_{*}$ を計算する。

# 事後分布のパラメータにより予測分布のパラメータを計算:式(3.19')
mu_star_hat <- a_hat / (a_hat + b_hat)
mu_star_hat

[1] 0.3039216

# 観測データにより予測分布のパラメータを計算:式(3.19')
mu_star_hat <- (sum(x_n) + a) / (N + a + b)
mu_star_hat

[1] 0.3039216

　予測分布のパラメータは、事後分布のパラメータ $\hat{a}, \hat{b}$ または観測データ $\mathbf{X}$ を用いて、次の式で計算できる。

$\displaystyle \begin{aligned} \hat{\mu}_{*} &= \frac{\hat{a}}{\hat{a} + \hat{b}} \\ &= \frac{\sum_{n=1}^N x_n + a}{N + a + b} \end{aligned} \tag{3.19'}$

　予測分布の確率を計算する。

# 予測分布の確率を計算:式(2.16)
predict_df <- tibble::tibble(
  x    = x_vec, # 確率変数
  mu   = mu_star_hat, # 成功確率パラメータ
  prob = c(1-mu_star_hat, mu_star_hat) # 確率
)
predict_df

# A tibble: 2 × 3
      x    mu  prob
  <dbl> <dbl> <dbl>
1     0 0.304 0.696
2     1 0.304 0.304

　生成分布のときと同様にして、ベルヌーイ分布に従う確率 $\mathrm{Bern}(x \mid \hat{\mu}_{*})$ を計算する。

　予測分布のグラフを作成する。

# 予測分布のラベルを作成
predict_param_lbl <- paste0(
  "list(", 
  "N == ",            N, ", ", 
  "mu[truth] == ",    round(mu_truth,    digits = 2), ", ", 
  "hat(mu)['*'] == ", round(mu_star_hat, digits = 5), 
  ")"
) |> 
  parse(text = _)

# 予測分布を作図
ggplot() + 
  geom_bar(
    data    = model_df, 
    mapping = aes(x = x, y = prob, linetype = "model"), 
    stat = "identity", position = "identity",
    fill = NA, color = "red", linewidth = 1
  ) + # 真の分布
  geom_bar(
    data    = predict_df, 
    mapping = aes(x = x, y = prob, linetype = "predict"), 
    stat = "identity", position = "identity", 
    fill = "purple", color = NA, alpha = 0.5
  ) + # 予測分布
  scale_x_continuous(breaks = x_vec, minor_breaks = FALSE) + # x軸目盛
  scale_linetype_manual(
    breaks = c("model", "predict"), 
    values = c("dashed", "blank"), 
    labels = c("true model", "predict distribution"), 
    name   = ""
  ) + # (凡例表示用)
  guides(
    linetype = guide_legend(override.aes = list(linewidth = 0.5))
  ) + 
  theme(
    legend.title           = element_text(size = 0), 
    legend.position        = "inside", 
    legend.position.inside = c(1, 1), 
    legend.justification   = c(1, 1), 
    legend.background      = element_rect(fill = alpha("white", alpha = 0.8))
  ) + 
  labs(
    title = "Bernoulli distribution", 
    subtitle = predict_param_lbl, 
    x = expression(x), 
    y = "probability"
  )

　真の分布(ベルヌーイ分布)を赤色(白抜き)の棒グラフ(破線)、予測分布(ベルヌーイ分布)を紫色(塗りつぶし)の棒グラフ(実線)で示す。

　データ数 $N$ が十分に大きいと、未観測データ $x_{*}$ の予測分布 $\mathrm{Bern}(x_{*} \mid \hat{\mu}_{*})$ の形状が真の分布 $\mathrm{Bern}(x_{n} \mid \mu_{\mathrm{truth}})$ に近付く。

　以上で、ベルヌーイモデルのベイズ推論を実装した。

スポンサードリンク

学習の推移

　次は、ベルヌーイ分布に対するベイズ推論を図で確認する。
　データ数を増やして分布の変化をアニメーションで確認する。
　作図コードについては「Suyama-Bayes/code/bernoulli_model/bayesian_inference/plot_parameter_updates.R at master · anemptyarchive/Suyama-Bayes · GitHub」を参照のこと。

データ数と分布の関係

　データ数 $N$ を変化させたときの事後分布 $p(\mu \mid \mathbf{X}, a, b)$ の変化をアニメーションにする。

　 $n$ 個のデータから求めた( $n$ 回更新した)事後分布(ベータ分布) $\mathrm{Beta}(\mu \mid a^{(n)}, b^{(n)})$ を紫色の曲線(実線)、 $n$ 番目の観測データ $x_n$ に対応する位置 $\mu = x_n$ を桃色の点で示す。

　「ベイズ推論の実装」では、 $N$ 個(複数)のデータ $\mathbf{X}$ を用いて、事後分布のパラメータ $\hat{a}, \hat{b}$ を一括更新した。
　 $n$ 番目(1つ)のデータ $x_n$ を用いて逐次更新する場合、 $n$ 回目の事後分布のパラメータ $a^{(n)}, b^{(n)}$ は、次の式で計算できる。

$\displaystyle \begin{aligned} a^{(n)} &\leftarrow x_n + a^{(n-1)} \\ b^{(n)} &\leftarrow 1 - x_n + b^{(n-1)} \end{aligned}$

　超パラメータの初期値(1回目の更新における事前分布のパラメータ)を $a^{(0)}, b^{(0)}$ として、 $a^{(n-1)}, b^{(n-1)}$ は、 $n-1$ 回目の更新における事後分布のパラメータであり、また $n$ 回目の更新における事前分布のパラメータを表す。

　データ数 $N$ が大きくなるのに応じて、パラメータ $\mu$ の事後分布 $\mathrm{Beta}(\mu \mid \hat{a}, \hat{b})$ のピークが真のパラメータ $\mu_{\mathrm{truth}}$ に近付いていくのを確認できる。

　データ数 $N$ を変化させたときの予測分布 $p(x_{*} \mid \mathbf{X}, a, b)$ の変化をアニメーションにする。

　 $n$ 個のデータから求めた( $n$ 回更新した)予測分布(ベルヌーイ分布) $\mathrm{Bern}(x_{*} \mid \mu_{*}^{(n)})$ を紫色(塗りつぶし)の棒グラフ(実線)、 $n$ 番目の観測データ $x_n$ を桃色の点で示す。

　「ベイズ推論の実装」では、 $N$ 個(複数)のデータ $\mathbf{X}$ を用いて、予測分布のパラメータ $\hat{\mu}_{*}$ を一括更新した。
　 $n$ 番目(1つ)のデータ $x_n$ を用いて逐次更新する場合、 $n$ 回目の予測分布のパラメータ $\mu_{*}^{(n)}$ は、次の式で計算できる。

$\displaystyle \mu_{*}^{(n)} = \frac{a^{(n)}}{a^{(n)} + b^{(n)}}$

　超パラメータの初期値 $a^{(0)}, b^{(0)}$ を用いて求めたパラメータを $\mu_{*}^{(0)}$ で表す。

　データ数 $N$ が大きくなるのに応じて、未観測データ $x_{*}$ の予測分布 $\mathrm{Bern}(x_{*} \mid \hat{\mu}_{*})$ の形状が真の分布 $\mathrm{Bern}(x_{n} \mid \mu_{\mathrm{truth}})$ に近付いていくのを確認できる。

観測データと分布の関係

　データ数 $N$ の変化による観測データと事後分布、予測分布の関係をアニメーションにする。

　真の分布の期待値 $\mathbb{E}[x_n]$ と真のパラメータ $\mu_{\mathrm{truth}}$ の各図(分布)に対応する位置を赤色の垂直線(破線)、観測データの標本平均 $\bar{x}$ を桃色の垂直線(破線)、事後分布の期待値 $\mathbb{E}[\mu]$ と予測分布の期待値 $\mathbb{E}[x_{*}]$ を紫色の垂直線(破線)で示す。

　生成分布(ベルヌーイ分布)、事後分布(ベータ分布)、予測分布(ベルヌーイ分布)の期待値は、それぞれ次の式で計算できる。

$\displaystyle \begin{align} \mathbb{E}_{\mathrm{Bern}(x_n \mid \mu_{\mathrm{truth}})}[x] &= \mu_{\mathrm{truth}} \\ \mathbb{E}_{\mathrm{Beta}(\mu \mid \hat{a}, \hat{b})}[\mu] &= \frac{\hat{a}}{\hat{a} + \hat{b}} \\ &= \hat{\mu}_{*} \\ \mathbb{E}_{\mathrm{Bern}(x_{*} \mid \hat{\mu}_{*})}[x] &= \hat{\mu}_{*} \end{align}$

　真の分布 $\mathrm{Bern}(x_n \mid \mu_{\mathrm{truth}})$ や真のパラメータ $\mu_{\mathrm{truth}}$ と、観測データ $\mathbf{X}$ 、事後分布 $\mathrm{Beta}(\mu \mid \hat{a}, \hat{b})$ 、予測分布 $\mathrm{Bern}(x_{*} \mid \hat{\mu}_{*})$ 、またそれぞれの統計量の対応関係を確認できる。

　以上で、ベルヌーイモデルのベイズ推論における学習推移を確認した。

　この記事では、ベルヌーイ分布に対するベイズ推論を実装した。次の記事では、カテゴリ分布を扱う。

参考文献

作者:須山敦志
講談社

おわりに

　大幅加筆修正の際に記事を分割しました。

2022.07.09：加筆修正しました。

　またまた書き直していきます。今回はグラフに載せる情報を増やすのが目的ですが、尤度と生成分布を勘違いしていたことに気付いたのでそれを早くどうにかしたいです。グラフを分かりやすくするために、作図コードが分かりにくくなってしまいました。その対応として、確率分布の計算・作図・乱数生成の記事を充実させたのでそちらも読んでみてください。
　さすがにこのシリーズの修正はこれで最後にしたい。