3.2.3：ポアソン分布の学習と予測【緑ベイズ入門のノート】

はじめに

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

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

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

【前節の内容】

www.anarchive-beta.com

【他の節一覧】

www.anarchive-beta.com

【この節の内容】

はじめに
3.2.3 ポアソン分布の学習と予測
おわりに

3.2.3 ポアソン分布の学習と予測

・事後分布の計算

　ポアソン分布に従うと仮定するN個の観測データ$\boldsymbol{X} = \{x_1, x_2, \cdots, x_N\}$を用いて、$\lambda$の事後分布を求めていく。

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

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

と分解できる。

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

$$ \begin{align*} p(\lambda | \boldsymbol{X}) &= \frac{ p(\boldsymbol{X} | \lambda) p(\lambda | a, b) }{ p(\boldsymbol{X}) } \\ &= \frac{ \left\{ \prod_{n=1}^N p(x_n | \lambda) \right\} p(\lambda | a, b) }{ p(\boldsymbol{X}) } \\ &= \frac{ \left\{ \prod_{n=1}^N \mathrm{Poi}(x_n | \lambda) \right\} \mathrm{Gam}(\lambda | a, b) }{ p(\boldsymbol{X}) } \\ &\propto \left\{ \prod_{n=1}^N \mathrm{Poi}(x_n | \lambda) \right\} \mathrm{Gam}(\lambda | a, b) \tag{3.35} \end{align*} $$

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

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

$$ \begin{align*} \ln p(\lambda | \boldsymbol{X}) &= \ln \prod_{n=1}^N \mathrm{Poi}(x_n | \lambda) + \ln \mathrm{Gam}(\lambda | a, b) - \ln p(\boldsymbol{X}) \\ &= \sum_{n=1}^N \ln \mathrm{Poi}(x_n | \lambda) + \ln \mathrm{Gam}(\lambda | a, b) + \mathrm{const.} \\ &= \sum_{n=1}^N \ln \left\{ \frac{\lambda^{x_n}}{x_n!} \exp(-\lambda) \right\} + \ln \Bigl\{ C_G(a, b) \lambda^{a-1} \exp(- b \lambda) \Bigr\} + \mathrm{const.} \\ &= \sum_{n=1}^N \{ x_n \ln \lambda - \ln x_n! - \lambda \} + \ln C_G(a, b) + (a - 1) \ln \lambda - b \lambda + \mathrm{const.} \\ &= \sum_{n=1}^N x_n \ln \lambda - N \lambda + (a - 1) \ln \lambda - b \lambda + \mathrm{const.} \\ &= \left( \sum_{n=1}^N x_n + a - 1 \right) \ln \lambda - (N + b) \lambda + \mathrm{const.} \tag{3.36} \end{align*} $$

【途中式の途中式】

式(3.36)の下から2行目の式を用いている。
対数をとると積が和に変わる。また、$\lambda$に影響しない項($- \ln p(\boldsymbol{X})$)を$\mathrm{const.}$とおく。
$\mathrm{Poi}(x_n | \lambda),\ \mathrm{Gam}(\lambda | a, b)$に、それぞれポアソン分布の定義式(2.37)、ガンマ分布の定義式(2.56)を用いて具体的な確率分布を代入する。
それぞれ対数をとる。$\ln (x^a y^b) = a \ln x + b \ln y$である。
$\lambda$と関係のない$- \sum_{n=1}^N \ln x_n!,\ \ln C_G(a, b)$を$\mathrm{const.}$にまとめる。
$\ln \lambda$と$\lambda$の項をまとめて式を整理する。

となる。

　式(3.36)について

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

とおき

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

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

$$ p(\lambda | \boldsymbol{X}) = C_G(\hat{a}, \hat{b}) \lambda^{\hat{a}} \exp ( - \hat{b} \lambda ) = \mathrm{Gam}(\lambda | \hat{a}, \hat{b}) \tag{3.37} $$

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

・予測分布の計算

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

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

$$ \begin{align*} p(x_{*}) &= \int p(x_{*} | \lambda) p(\lambda) d\lambda \\ &= \int \mathrm{Poi}(x_{*} | \lambda) \mathrm{Gam}(\lambda | a, b) d\lambda \tag{3.39}\\ &= \int \frac{\lambda^{x_{*}}}{x_{*}!} \exp(- \lambda) C_G(a, b) \lambda^{a-1} \exp(- b \lambda) d\lambda \\ &= \frac{C_G(a, b)}{x_{*}!} \int \lambda^{x_{*}+a-1} \exp\Bigl( - (1 + b) \lambda \Bigr) d\lambda \tag{3.40} \end{align*} $$

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

$$ \int \lambda^{x_{*}+a-1} \exp\Bigl( - (1 + b) \lambda \Bigr) d\lambda = \frac{1}{C_G(x_{*} + a, 1 + b)} $$

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

$$ p(x_{*}) = \frac{ C_G(a, b) }{ x_{*}! C_G(x_{*} + a, 1 + b) } \tag{3.42} $$

となる。更にガンマ分布の正規化項(2.57)を用いて、式を整理すると

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

となる。
　この式の形状の確率分布は負の二項分布と呼ばれる。

　この式について

$$ \begin{align*} r &= a \\ p &= \frac{1}{1 + b} \tag{3.44} \end{align*} $$

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

$$ p(x_{*}) = \frac{ \Gamma(x_{*} + r) }{ x_{*}! \Gamma(r) } (1 - p)^a p^{x_{*}} = \mathrm{NB}(x_{*} | r, p) \tag{3.43} $$

パラメータ$r,\ p$を持つ負の二項分布となることが分かる。

　予測分布のパラメータ$r,\ p$を構成する$a,\ b$について、事後分布のパラメータ(3.38)を用いると

$$ \begin{aligned} \hat{r} &= \hat{a} \\ &= \sum_{n=1}^N x_n + a \end{aligned} $$

また

$$ \begin{aligned} \hat{p} &= \frac{1}{1 + \hat{b}} \\ &= \frac{1}{N + 1 + b} \end{aligned} $$

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

$$ p(x_{*} | \boldsymbol{X}) = \frac{ \Gamma(x_{*} + \hat{r}) }{ x_{*}! \Gamma(\hat{r}) } (1 - \hat{p})^a \hat{p}^{x_{*}} = \mathrm{NB}(\hat{x}_{*} | \hat{r}, \hat{p}) $$

が得られる。

・Rでやってみよう

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

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

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

・事後分布

## パラメータの初期値を設定
# 観測モデルのパラメータ
lambda_truth <- 4

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

# 試行回数
N <- 50

　パラメータ$\lambda$をlambda_truthとします。この値を推定するのが目的です。

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

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

# ポアソン分布に従うデータを生成
x_n <- rpois(n = N ,lambda = lambda_truth)

　ポアソン分布に従う乱数を発生させる関数rpois()を使って、ランダムにデータを生成します。
　試行回数の引数nにはN、パラメータ引数lambdaにはlambda_truthを指定します。

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

# 観測データを確認
x_n_df <- x_n %>% 
  table() %>% 
  as_tibble()
colnames(x_n_df) <- c("n", "x")
x_n_df$n <- as.numeric(x_n_df$n)
head(x_n_df)
## # A tibble: 6 x 2
##       n     x
##   <dbl> <int>
## 1     0     1
## 2     1     2
## 3     2     9
## 4     3     7
## 5     4    11
## 6     5     7

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

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

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

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

　$\lambda$が取り得る値をseq()を使って値を用意します。ただしlambda_truthの値によって必要な範囲が変わります。ここでは0からlambda_truthの2倍までとします。そこで第1引数(最小値)には0、第2引数(最大値)には2 * lambda_truthを指定します。
　lambdaの各値に対してガンマ分布の定義式(2.56)の計算を行い、確率密度を求めます。ただし値が大きくなるとgamma()で計算できなくなるため、対数をとって計算することにします。なので最後にexp()で値を戻します。

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

head(posterior_df)
## # A tibble: 6 x 3
##   lambda   C_G density
##    <dbl> <dbl>   <dbl>
## 1  0     -95.5       0
## 2  0.001 -95.5       0
## 3  0.002 -95.5       0
## 4  0.003 -95.5       0
## 5  0.004 -95.5       0
## 6  0.005 -95.5       0

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

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

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

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

f:id:anemptyarchive:20200305150802p:plain — $\lambda$の事後分布：ガンマ分布

・予測分布

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

　予測分布のパラメータの計算式の計算を行い、$\hat{r},\ \hat{p}$をそれぞれr_hat, p_hatとします。

r_hat <- sum(x_n) + a
p_hat <- 1 / (N + 1 + b)

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

# 予測分布を計算
predict_df <- tibble(
  x = seq(0, 4 * lambda_truth),  # 作図用の値
  C_NB = lgamma(x + r_hat) - lgamma(x + 1) - lgamma(r_hat), # 正規化項(対数)
  prob = exp(C_NB + r_hat * log(1 - p_hat) + x * log(p_hat)) # 確率
)

　事後分布と同様に、$x_{*}$が取り得る値をseq()を使って用意します。
　xの各値となる確率は、負の二項分布の定義式(3.43)で計算します。

　推定結果を確認してみましょう。

head(predict_df)
## # A tibble: 6 x 3
##       x  C_NB   prob
##   <int> <dbl>  <dbl>
## 1     0  0    0.0148
## 2     1  5.38 0.0617
## 3     2 10.1  0.129 
## 4     3 14.4  0.182 
## 5     4 18.4  0.192 
## 6     5 22.2  0.163

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

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

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

f:id:anemptyarchive:20200305150906p:plain — $x_{*}$の予測分布：負の二項分布

・おまけ

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

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

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

## パラメータの初期値を設定
# 観測モデルのパラメータ
lambda_truth <- 4

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

# 試行回数
N <- 50

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

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

# 初期値による予測分布を計算
predict_df <- tibble(
  x = seq(0, 4 * lambda_truth),  # 作図用の値
  C_NB = lgamma(x + r) - lgamma(x + 1) - lgamma(r), # 正規化項(対数)
  prob = exp(C_NB + r * log(1 - p) + x * log(p)), # 確率
  N = 0 # 試行回数
)

# パラメータを推定
x_n <- rep(0, N) # 受け皿
for(n in 1:N){
  
  # ポアソン分布に従うデータを生成
  x_n[n] <- rpois(n = 1 ,lambda = lambda_truth)
  
  # ハイパーパラメータを更新
  a <- sum(x_n[n] * 1) + a
  b <- 1 + b
  
  # 事後分布を推定
  tmp_posterior_df <- tibble(
    lambda = seq(0, 2 * lambda_truth, by = 0.001),  # 作図用のlamndaの値
    C_G = a * log(b) - lgamma(a),   # 正規化項(対数)
    density = exp(C_G + (a - 1) * log(lambda) - b * lambda), # 確率密度
    N = n # 試行回数
  )
  
  # 予測分布のパラメータを計算
  r <- a
  p <- 1 / (b + 1)
  
  # 予測分布を計算
  tmp_predict_df <- tibble(
    x = seq(0, 4 * lambda_truth),  # 作図用の値
    C_NB = lgamma(x + r) - lgamma(x + 1) - lgamma(r), # 正規化項(対数)
    prob = exp(C_NB + r * log(1 - p) + x * log(p)), # 確率
    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(lambda, density)) + 
  geom_line(color = "#56256E") + # 折れ線グラフ
  geom_vline(aes(xintercept = lambda_truth), 
             color = "red", linetype = "dashed") + # 垂直線
  transition_manual(N) + # フレーム
  labs(title = "Gamma Distribution", 
       subtitle = "N= {current_frame}", 
       x = expression(lambda)) # ラベル

# 描画
animate(posterior_graph)

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

# 描画
animate(predict_graph)

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

　各データによってどのように学習する(推定値が変化する)のかを確認するため、こちらはforループで1データずつ処理します。
　よって生成するデータ数として設定したNがイタレーション数になります。

　パラメータの推定値について、$\hat{a},\ \hat{b}$に対応するa_hat, b_hatを新たに作るのではなく、a, bをイタレーションごとに更新(上書き)していきます。
　それに伴い、事後分布のパラメータの計算式(3.38)の$\sum_{n=1}^N$の計算は、forループによってN回繰り返しx_n[n]を加えることで行います。n回目のループ処理のときには、n-1回分のx_n[n](bの場合は1)が既にaとbに加えられているわけです。