3.2.1：ベルヌーイ分布の学習と予測【緑ベイズ入門のノート】

はじめに

　『ベイズ推論による機械学習入門』の学習時のノートです。基本的な内容は「数式の行間を読んでみた」とそれを「Rで組んでみた」になります。「数式」と「プログラム」から理解するのが目標です。

　この記事は3.2.1項の内容になります。尤度関数をベルヌーイ分布、事前分布をベータ分布とした場合のパラメータの事後分布を導出し、その学習した事後分布を用いた予測分布を導出します。またその学習過程をR言語で実装します。

　省略してる内容等ありますので、本と併せて読んでください。初学者な自分が理解できるレベルまで落として書き下していますので、分かる人にはかなりくどくなっています。同じような立場の人のお役に立てば幸いです。

【他の節一覧】

www.anarchive-beta.com

【この節の内容】

はじめに
3.2.1 ベルヌーイ分布の学習と予測
おわりに

3.2.1 ベルヌーイ分布の学習と予測

・事後分布の計算

　ベルヌーイ分布に従うと仮定するN個の観測データ$\boldsymbol{X} = \{x_1, x_2, \cdots, x_N \}$を用いて、パラメータ$\mu$の事後分布を求めていく。

　まずは観測モデル$p(\boldsymbol{X} | \mu)$について確認する。離散値データ$\boldsymbol{X} = \{x_1, x_2, \cdots, x_N \}$はそれぞれ独立に発生しているとの仮定の下で、観測モデルは

$$ \begin{aligned} p(\boldsymbol{X} | \mu) &= p(x_1, x_2, \cdots, x_N | \mu) \\ &= p(x_1 | \mu) p(x_2 | \mu) \cdots p(x_N | \mu) \\ &= \prod_{n=1}^N p(x_n | \mu) \end{aligned} $$

と分解できる。

　これを用いて、$\mu$の事後分布$p(\mu | \boldsymbol{X})$はベイズの定理より

$$ \begin{align} p(\mu | \boldsymbol{X}) &= \frac{ p(\boldsymbol{X} | \mu) p(\mu) }{ p(\boldsymbol{X}) } \\ &= \frac{ \left\{ \prod_{n=1}^N p(x_n | \mu) \right\} p(\mu) }{ p(\boldsymbol{X}) } \\ &= \frac{ \left\{ \prod_{n=1}^N \mathrm{Bern}(x_n | \mu) \right\} \mathrm{Beta}(\mu | a, b) }{ p(\boldsymbol{X}) } \\ &\propto \left\{ \prod_{n=1}^N \mathrm{Bern}(x_n | \mu) \right\} \mathrm{Beta}(\mu | a, b) \tag{3.12} \end{align} $$

となる。分母の$p(\boldsymbol{X})$は$\mu$に影響しないため省略して、比例関係にのみ注目する。省略した部分については、最後に正規化することで対応できる。

　次にこの分布の具体的な形状を明らかにしていく。対数をとって指数部分の計算を分かりやすくして進めると

$$ \begin{align} \ln p(\mu | \boldsymbol{X}) &= \ln \prod_{n=1}^N \mathrm{Bern}(x_n | \mu) + \ln \mathrm{Beta}(\mu | a, b) - \ln p(\boldsymbol{X}) \\ &= \sum_{n=1}^N \ln \mathrm{Bern}(x_n | \mu) + \ln \mathrm{Beta}(\mu | a, b) + \mathrm{const.} \\ &= \sum_{n=1}^N \ln \left\{ \mu^{x_n} (1 - \mu)^{1-x_n} \right\} + \ln \left\{ C_B(a, b) \mu^{a-1} (1 - \mu)^{b-1} \right\} + \mathrm{const.} \\ &= \sum_{n=1}^N \left\{ x_n \ln \mu + (1 - x_n) \ln(1 - \mu) \right\} \\ &\qquad + \ln C_B(a, b) + (a - 1) \ln \mu + (b - 1) \ln (1 - \mu) + \mathrm{const.} \\ &= \sum_{n=1}^N x_n \ln \mu + \left(N - \sum_{n=1}^N x_n\right) \ln(1 - \mu) \\ &\qquad + (a - 1) \ln \mu + (b - 1) \ln (1 - \mu) + \mathrm{const.} \\ &= \left( \sum_{n=1}^N x_n + a - 1 \right) \ln \mu + \left( N - \sum_{n=1}^N x_n + b - 1 \right) \ln (1 - \mu) + \mathrm{const.} \tag{3.13} \end{align} $$

【途中式の途中式】

式(3.12)の下から2行目の式を用いている。
$\mu$に影響しない項(ここでは$- \ln p(\boldsymbol{X})$)を$\mathrm{const.}$とおく。

　総乗部分の対数をとると

$$ \begin{aligned} \ln p(\boldsymbol{X} | \mu) &= \ln \prod_{n=1}^N p(x_n | \mu) \\ &= \ln \Bigl\{ p(x_1 | \mu) * p(x_2 | \mu) * \cdots * p(x_N | \mu) \Bigr\} \\ &= \ln p(x_1 | \mu) + \ln p(x_2 | \mu) + \cdots + \ln p(x_N | \mu) \\ &= \sum_{n=1}^N \ln p(x_n | \mu) \end{aligned} $$

のように総和になる。

$\mathrm{Bern}(x_n | \mu),\ \mathrm{Beta}(\mu | a, b)$に、それぞれベルヌーイ分布の定義式(2.16)、ベータ分布の定義式(2.41)を用いて具体的な確率分布を代入する。
それぞれ対数をとる。$\ln (x^a y^b) = a \ln x + b \ln y$である。
$\sum_{n=1}^N$に関する括弧を展開し、また$\mu$に影響しない項($\ln C_B(a, b)$)を$\mathrm{const.}$にまとめる。
$\mu$と$(1 - \mu)$の項をまとめて式を整理する。

となる。

　式(3.13)について

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

とおき

$$ \ln p(\mu | \boldsymbol{X}) = (\hat{a} - 1) \ln \mu + (\hat{b} - 1) \ln (1 - \mu) + \mathrm{const.} $$

更に$\ln$を外して、$\mathrm{const.}$を正規化項に置き換える(正規化する)と

$$ p(\mu | \boldsymbol{X}) = C_B(\hat{a}, \hat{b}) \mu^{\hat{a} - 1} (1 - \mu)^{\hat{b} - 1} = \mathrm{Beta}(\mu | \hat{a}, \hat{b}) \tag{3.14} $$

となる。式の形状から、事後分布$p(\mu | \boldsymbol{X})$がパラメータ$\hat{a},\ \hat{b}$を持つベータ分布であることを確認できる。
　また事後分布のパラメータの計算式(3.15)によって、ハイパーパラメータ$a,\ b$が更新される。

・予測分布の計算

　続いてベルヌーイ分布に従う未観測のデータ$x_{*}$に対する予測分布を求めていく。

　(観測データによる学習を行っていない)事前分布$p(\mu)$を用いて、パラメータ$\mu$を周辺化することで予測分布$p(x_{*})$となる。

$$ \begin{align} p(x_{*}) &= \int p(x_{*} | \mu) p(\mu) d\mu \\ &= \int \mathrm{Bern}(x_{*} | \mu) \mathrm{Beta}(\mu | a, b) d\mu \\ &= \int \mu^{x_{*}} (1 - \mu)^{1-x_{*}} C_B(a, b) \mu^{a-1} (1 - \mu)^{b-1} d\mu \\ &= C_B(a, b) \int \mu^{x_{*}+a-1} (1 - \mu)^{1-x_{*}+b-1} d\mu \tag{3.16} \end{align} $$

　積分部分に注目すると、パラメータ$x_{*} + a,\ 1 - x_{*} + b$を持つ正規化項のないベータ分布の形をしていることが分かる。よって、ベータ分布の定義式(2.41)を用いて

$$ \int \mu^{x_{*}+a-1} (1 - \mu)^{1-x_{*}+b-1} d\mu = \frac{1}{C_B(x_{*} + a, 1 - x_{*} + b)} \tag{3.17} $$

正規化項の逆数に変形できる(そもそもこの部分の逆数が正規化項である)。これを式(3.16)に代入すると、予測分布$p(x_{*})$は

$$ p(x_{*}) = \frac{C_B(a, b)}{C_B(x_{*} + a, x_{*} + b)} $$

となる。更にベータ分布の正規化項(2.42)を用いて、式を整理すると

$$ \begin{align} p(x_{*}) &= C_B(a, b) \frac{1}{C_B(x_{*} + a, x_{*} + b)} \\ &= \frac{ \Gamma(a + b) }{ \Gamma(a) \Gamma(b) } \frac{ \Gamma(x_{*} + a) \Gamma(1 - x_{*} + b) }{ \Gamma(x_{*} + a + 1 - x_{*} + b) } \\ &= \frac{ \Gamma(a + b) }{ \Gamma(a) \Gamma(b) } \frac{ \Gamma(x_{*} + a) \Gamma(1 - x_{*} + b) }{ \Gamma(a + b + 1) } \\ &= \frac{ \Gamma(a + b) }{ \Gamma(a) \Gamma(b) } \frac{ \Gamma(x_{*} + a) \Gamma(1 - x_{*} + b) }{ (a + b) \Gamma(a + b) } \\ &= \frac{ \Gamma(x_{*} + a) \Gamma(1 - x_{*} + b) }{ \Gamma(a) \Gamma(b) (a + b) } \tag{3.18} \end{align} $$

となる。ガンマ関数の性質$\Gamma(x + 1) = x \Gamma(x)$を用いている。

　$x_{*}$は0か1しかとらないため、場合分けしてこの式を整理する。$x_{*} = 1$のとき

$$ \begin{align} p(x_{*} = 1) &= \frac{ \Gamma(1 + a) \Gamma(1 - 1 + b) }{ \Gamma(a) \Gamma(b) (a + b) } \tag{3.18}\\ &= \frac{ a \Gamma(a) \Gamma(b) }{ \Gamma(a) \Gamma(b) (a + b) } \\ &= \frac{a}{a + b} \tag{3.19} \end{align} $$

となる。また$x_{*} = 0$のとき

$$ \begin{align} p(x_{*} = 0) &= \frac{ \Gamma(0 + a) \Gamma(1 - 0 + b) }{ \Gamma(a) \Gamma(b) (a + b) } \tag{3.18}\\ &= \frac{ \Gamma(a) b \Gamma(b) }{ \Gamma(a) \Gamma(b) (a + b) } \\ &= \frac{b}{a + b} \tag{3.20} \end{align} $$

となる。これは

$$ \frac{b}{a + b} = 1 - \frac{a}{a + b} $$

と置き換えることができる。

　従って、$x^0 = 1$であることを利用して式(3.19)と式(3.20)をまとめると、式(3.18)は

$$ \begin{aligned} p(x_{*}) &= \Bigl( \frac{a}{a + b} \Bigr)^{x_{*}} \Bigl( \frac{b}{a + b} \Bigr)^{1-x_{*}} \\ &= \Bigl( \frac{a}{a + b} \Bigr)^{x_{*}} \Bigl( 1 - \frac{a}{a + b} \Bigr)^{1-x_{*}} \end{aligned} $$

と書き換えることができる。($x_{*} = 0$のとき前の項が1となり、$x_{*} = 1$のとき後ろの項が1となる。)

　この式について

$$ \mu_{*} = \frac{a}{a + b} $$

とおくと、予測分布$p(x_{*})$は

$$ p(x_{*}) = \mu_{*}^{x_{*}} (1 - \mu_{*})^{1-x_{*}} = \mathrm{Bern}(x_{*} | \mu_{*}) \tag{3.21} $$

となる。式の形状から、$p(x_{*})$はパラメータ$\mu_{*}$を持つベルヌーイ分布であることが分かる。
　ちなみに$\frac{a}{a + b}$は、ベータ分布$\mathrm{Beta}(\mu | a, b)$の期待値$\mathbb{E}[\mu]$である。

　予測分布のパラメータ$\mu_{*}$を構成する$a,\ b$について、それぞれ事後分布のパラメータ(3.15)を用いると

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

となり、(観測データ$\boldsymbol{X}$によって学習した)事後分布$p(\mu | \boldsymbol{X})$を用いた予測分布$p(x_{*} | \boldsymbol{X})$

$$ p(x_{*} | \boldsymbol{X}) = \hat{\mu}_{*}^{x_{*}} (1 - \hat{\mu}_{*})^{1-x_{*}} = \mathrm{Bern}(x_{*} | \hat{\mu}_{*}) \tag{3.22} $$

が得られる。

・Rでやってみよう

　各確率分布に従いランダムに生成したデータを用いて、パラメータを推定してみましょう。

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

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

・事後分布

## パラメータの初期値を指定
# 観測モデルのパラメータ
mu_truth <- 0.25

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

# 試行回数
N <- 50

　$x_n = 1$となる確率$\mu$をmu_truthとします。この値を推定するのが目的です。

　事前分布のパラメータ(ハイパーパラメータの初期値)$a,\ b$をそれぞれa, bとします。

　生成するデータ数$N$をNとします。

# ベルヌーイ分布に従うデータを生成
x_n <- rbinom(n = N, size = 1, prob = mu_truth)

　二項分布に従う乱数を発生させる関数rbinom()を使って、ランダムにデータを生成します。これをベルヌーイ分布とするには、size引数に1を指定します。
　試行回数の引数nにはN、確率の引数probにはmu_truthを指定します。

　サンプルを確認してみましょう。

# 観測データを確認
table(x_n)
## x_n
##  0  1 
## 37 13

このデータを用いて事後分布のパラメータを計算します。

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

　事後分布のパラメータを計算式(3.15)の計算を行い、$\hat{a},\ \hat{b}$をそれぞれa_hat, b_hatとします。

# 事後分布を計算
posterior_df <- tibble(
  mu = seq(0, 1, by = 0.001),  # 作図用の値
  C_B = lgamma(a_hat + b_hat) - lgamma(a_hat) - lgamma(b_hat),  # 正規化項(対数)
  density = exp(C_B + (a_hat - 1) * log(mu) + (b_hat - 1) * log(1 - mu)) # 確率密度
)

　$\mu$が取り得る0から1までの値をseq(0, 1)を使って用意します。by引数で刻み幅を指定できるので、グラフが粗かったり処理が重くなったりする場合はこの値を調整してください。

　muの各値に対してベータ分布の定義式(2.41)の計算を行い、確率密度を求めます。ただし値が大きくなるとgamma()で計算できなくなるため、対数をとって計算することにします。なので最後にexp()で値を戻します。

　計算結果は作図用にデータフレームにまとめておきます。推定結果を確認してみましょう。

head(posterior_df)
## # A tibble: 6 x 3
##      mu   C_B  density
##   <dbl> <dbl>    <dbl>
## 1 0      30.5 0.      
## 2 0.001  30.5 1.74e-26
## 3 0.002  30.5 1.38e-22
## 4 0.003  30.5 2.58e-20
## 5 0.004  30.5 1.05e-18
## 6 0.005  30.5 1.84e-17

このデータフレームを用いてグラフを描きます。

# 作図
ggplot(posterior_df, aes(mu, density)) + 
  geom_line(color = "#56256E") + # 折れ線グラフ
  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)) # ラベル

　ggplot2パッケージを利用してグラフを作成します。

　推定値は折れ線グラフgeom_line()を用いて描きます。
　パラメータ$\mu$の真の値を垂直線geom_vline()で示します。

f:id:anemptyarchive:20200304221005p:plain — $\mu$の事後分布：ベータ分布

・予測分布

# 予測分布のパラメーターを計算
mu_hat <- a_hat / (a_hat + b_hat)

　予測分布のパラメータを計算式の計算を行い、$\hat{\mu}$をmu_hatとします。

mu_hat <- (sum(x_n) + a) / (N + a + b)

　このように、観測データx_nと事前分布のパラメータa, bを使って計算することもできます。

# 予測分布を計算
predict_df <- tibble(
  x = c(0, 1),  # 作図用の値
  prob = mu_hat^x * (1 - mu_hat)^(1 - x) # 確率
)

　$x_{*}$が取り得る0と1の値を用意します。
　xの各値となる確率をベルヌーイ分布の定義式(2.16)で計算します。

　推定結果は次のようになります。

head(predict_df)
## # A tibble: 2 x 2
##       x  prob
##   <dbl> <dbl>
## 1     0 0.731
## 2     1 0.269

　こちらもggplot2パッケージを利用して作図します。

# 作図
ggplot(predict_df, aes(x, prob)) + 
  geom_bar(stat = "identity", position = "dodge", fill = "#56256E") + # 棒グラフ
  labs(title = "Bernoulli Distribution", 
       subtitle = paste0("N=", N, ", mu_hat=", round(mu_hat, 2))) # ラベル

　棒グラフはgeom_bar()を使います。

f:id:anemptyarchive:20200304221255p:plain — $x_{*}$の予測分布：ベルヌーイ分布

・おまけ

　gganimateパッケージを利用して、パラメータの推定値の推移のgif画像を作成するためのコードです。

・コード(クリックで展開)

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

## パラメータの初期値を指定
# 観測モデルのパラメータ
mu_truth <- 0.25

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

# 試行回数
N <- 50

# 事前分布を計算
posterior_df <- tibble(
  mu = seq(0, 1, by = 0.001),  # 作図用の値
  C_B = lgamma(a + b) - lgamma(a) - lgamma(b),  # 正規化項(対数)
  density = exp(C_B + (a - 1) * log(mu) + (b - 1) * log(1 - mu)), # 確率密度
  N = 0  # 試行回数
)

# 初期値による予測分布のパラメーターを計算
mu_star <- a / (a + b)

# 初期値による予測分布を計算
predict_df <- tibble(
  x = c(0, 1),  # 作図用の値
  prob = mu_star^x * (1 - mu_star)^(1 - x),  # 確率
  N = 0 # 試行回数
)

# パラメータを推定
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
    C_B = lgamma(a + b) - lgamma(a) - lgamma(b),  # 正規化項(対数)
    density = exp(C_B + (a - 1) * log(mu) + (b - 1) * log(1 - mu)), # 確率密度
    N = n  # 試行回数
  )
  
  # 予測分布のパラメーターを更新
  mu_star <- a / (a + b)
  
  # 予測分布を計算
  tmp_predict_df <- tibble(
    x = c(0, 1),  # 作図用の値
    prob = mu_star^x * (1 - mu_star)^(1 - x),  # 確率
    N = n # 試行回数
  )
  
  # 推定結果を結合
  posterior_df <- rbind(posterior_df, tmp_posterior_df)
  predict_df <- rbind(predict_df, tmp_predict_df)
}

# 観測データを確認
table(x_n)

## 事後分布
# 作図
posterior_graph <- ggplot(posterior_df, aes(mu, density)) + 
  geom_line(color = "#56256E") +  # 折れ線グラフ
  geom_vline(aes(xintercept = mu_truth), 
             color = "red", linetype = "dashed") + # 垂直線
  transition_manual(N) + # フレーム
  labs(title = "Beta Distribution", 
       subtitle = "N= {current_frame}", 
       x = expression(mu)) # ラベル

# 描画
animate(posterior_graph)

## 予測分布
# 作図
predict_graph <- ggplot(predict_df, aes(x, prob)) + 
  geom_bar(stat = "identity", position = "dodge", fill = "#56256E") + # 棒グラフ
  transition_manual(N) + # フレーム
  labs(title = "Bernoulli Distribution", 
       subtitle = "N= {current_frame}") # ラベル

# 描画
animate(predict_graph)

　異なる点のみを簡単に解説します。

　観測された各データによってどのように学習する(推定値(の分布)が変化する)のかを確認するため、こちらは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に加えられているわけです。