はじめに

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

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

【前節の内容】

【他の節の内容】

【この節の内容】

3.2.3 ポアソン分布のベイズ推論の実装

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

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

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

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

ベイズ推論の実装

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

生成分布の設定

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

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

# 真のパラメータを指定
lambda_truth <- 4
lambda_truth

[1] 4

　ポアソン分布の期待値パラメータ(正の値) $\lambda_{\mathrm{truth}} \gt 0$ を指定する。
　 $\lambda_{\mathrm{truth}}$ がポアソンモデルにおける真のパラメータであり、この値を求めるのがここでの目的(学習)である。

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

# x軸の範囲を設定
x_min <- 0
#x_max <- 15
u <- 5
k <- 3
x_max <- lambda_truth |> # 基準値を指定
  (\(.) {. * k})() |> # 定数倍
  #(\(.) {max(., x_n)})() |> # サンプルと比較
  (\(.) {ceiling(. /u)*u})() # u単位で切り上げ

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

[1] 0
[1] 15
[1] 0 1 2 3 4 5

　この例では、指定したパラメータ(または生成したデータ)を使って、範囲を設定している。

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

# 生成分布の確率を計算:式(2.37)
model_df <- tibble::tibble(
  x      = x_vec, # 確率変数
  lambda = lambda_truth, # 期待値パラメータ
  prob   = dpois(x = x, lambda = lambda) # 確率
)
model_df

# A tibble: 16 × 3
       x lambda      prob
   <dbl>  <dbl>     <dbl>
 1     0      4 0.0183   
 2     1      4 0.0733   
 3     2      4 0.147    
 4     3      4 0.195    
 5     4      4 0.195    
 6     5      4 0.156    
 7     6      4 0.104    
 8     7      4 0.0595   
 9     8      4 0.0298   
10     9      4 0.0132   
11    10      4 0.00529  
12    11      4 0.00192  
13    12      4 0.000642 
14    13      4 0.000197 
15    14      4 0.0000564
16    15      4 0.0000150

　 $x$ の値ごとに、ポアソン分布に従う確率 $\mathrm{Poi}(x \mid \lambda_{\mathrm{truth}})$ を計算する。
　ポアソン分布の確率質量関数 dpois() の確率変数の引数 x に $x$ 、パラメータの引数 lambda に $\lambda_{\mathrm{truth}}$ を指定する。

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

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

# 生成分布を作図
ggplot() + 
  geom_bar(
    data    = model_df, 
    mapping = aes(x = x, y = prob, fill = "model"), 
    stat = "identity", position = "identity"
  ) + # 真の分布
  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 = "Poisson distribution", 
    subtitle = model_param_lbl, 
    x = expression(x), 
    y = "probability"
  )

　真の分布(ポアソン分布)を赤色の棒グラフで示す。

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

データの生成

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

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

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

# 観測データを生成
x_n <- rpois(n = N ,lambda = lambda_truth)
head(x_n)

[1] 5 4 7 2 4 6

　データ数 $N$ を指定して、ポアソン分布に従う乱数 $\mathbf{X} = \{x_1, \cdots, x_N\}$ を生成する。
　ポアソン分布の乱数生成関数 rpois() のサンプルサイズの引数 n に $N$ 、パラメータの引数 lambda に $\lambda_{\mathrm{truth}}$ を指定する。

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

# 観測データを集計
obs_df <- tidyr::tibble(
  x = x_n # 観測値
) |> 
  dplyr::count(
    x, name = "freq" # 度数
  ) |> 
  dplyr::mutate(
    rel_freq = freq / N
  ) |> 
  tidyr::complete(
    x = x_vec, 
    fill = list(freq = 0, rel_freq = 0)
  ) # 未観測値を補完
obs_df

# A tibble: 16 × 3
       x  freq rel_freq
   <dbl> <int>    <dbl>
 1     0     0     0   
 2     1     2     0.02
 3     2    19     0.19
 4     3    12     0.12
 5     4    22     0.22
 6     5    19     0.19
 7     6    17     0.17
 8     7     7     0.07
 9     8     1     0.01
10     9     0     0   
11    10     0     0   
12    11     1     0.01
13    12     0     0   
14    13     0     0   
15    14     0     0   
16    15     0     0

　 $x$ の値ごとに、観測データ x_n に含まれる要素数をカウントして、度数 $N_x$ と相対度数 $\frac{N_x}{N}$ を求める。

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

# 観測データのラベルを作成
obs_param_lbl <- paste0(
  "list(", 
  "N == ",             N, ", ", 
  "lambda[truth] == ", round(lambda_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"), #na.value = NA, 
    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 = "Poisson 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{Poi}(x_n \mid \lambda_{\mathrm{truth}})$ に近付く。

事前分布の設定

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

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

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

[1] 1
[1] 1

　ガンマ分布の形状パラメータ(正の値) $a \gt 0$ 、尺度パラメータ(正の値) $b \gt 0$ を指定する。

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

# λ軸の範囲を設定
lambda_min <- 0
#lambda_max <- 15
u <- 5
k <- 3
lambda_max <- lambda_truth |> # 基準値を指定
  (\(.) {. * k})() |> # 定数倍
  (\(.) {ceiling(. /u)*u})() # u単位で切り上げ

# λ軸の値を作成
lambda_vec <- seq(from = lambda_min, to = lambda_max, length.out = 1001)
lambda_min; lambda_max; head(lambda_vec)

[1] 0
[1] 15
[1] 0.000 0.015 0.030 0.045 0.060 0.075

　この例では、指定したパラメータを使って、範囲を設定している。

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

# 事前分布の確率密度を計算:式(2.56)
prior_df <- tibble::tibble(
  lambda = lambda_vec, # 確率変数
  a      = a, # 形状パラメータ
  b      = b, # 尺度パラメータ
  dens   = dgamma(x = lambda, shape = a, rate = b) # 確率密度
)
prior_df

# A tibble: 1,001 × 4
   lambda     a     b  dens
    <dbl> <dbl> <dbl> <dbl>
 1  0         1     1 1    
 2  0.015     1     1 0.985
 3  0.03      1     1 0.970
 4  0.045     1     1 0.956
 5  0.06      1     1 0.942
 6  0.075     1     1 0.928
 7  0.09      1     1 0.914
 8  0.105     1     1 0.900
 9  0.12      1     1 0.887
10  0.135     1     1 0.874
# ℹ 991 more rows

　 $\lambda$ の値ごとに、ガンマ分布に従う確率密度 $\mathrm{Gam}(\lambda \mid a, b)$ を計算する。
　ガンマ分布の確率密度関数 dgamma() の確率変数の引数 x に $\lambda$ 、形状パラメータの引数 a に $a$ 、尺度パラメータの引数 scale に $b$ を指定する。

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

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

# 事前分布を作図
ggplot() + 
  geom_vline(
    mapping = aes(xintercept = lambda_truth, linetype = "model"), 
    color = "red", linewidth = 1
  ) + # 真のパラメータ
  geom_line(
    data    = prior_df, 
    mapping = aes(x = lambda, y = dens, linetype = "prior"), 
    color = "purple", linewidth = 1
  ) + # 事前分布
  scale_x_continuous(
    sec.axis = sec_axis(
      transform = ~ ., 
      breaks    = lambda_truth, 
      labels    = expression(lambda[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 = "Gamma distribution", 
    subtitle = prior_param_lbl, 
    x = expression(lambda), 
    y = "density"
  )

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

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

事後分布の計算

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

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

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

[1] 430
[1] 101

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

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

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

# 事後分布の確率密度を計算:式(2.56)
posterior_df <- tibble::tibble(
  lambda = lambda_vec, # 確率変数
  a      = a_hat, # 形状パラメータ
  b      = b_hat, # 尺度パラメータ
  dens   = dgamma(x = lambda, shape = a, rate = b) # 確率密度
)
posterior_df

# A tibble: 1,001 × 4
   lambda     a     b  dens
    <dbl> <dbl> <dbl> <dbl>
 1  0       430   101     0
 2  0.015   430   101     0
 3  0.03    430   101     0
 4  0.045   430   101     0
 5  0.06    430   101     0
 6  0.075   430   101     0
 7  0.09    430   101     0
 8  0.105   430   101     0
 9  0.12    430   101     0
10  0.135   430   101     0
# ℹ 991 more rows

　事前分布のときと同様にして、ガンマ分布に従う確率密度 $\mathrm{Gam}(\lambda \mid \hat{a}, \hat{b})$ を計算する。

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

# 事後分布のラベルを作成
posterior_param_lbl <- paste0(
  "list(", 
  "N == ",             N, ", ", 
  "lambda[truth] == ", round(lambda_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 = lambda_truth, linetype = "model"), 
    color = "red", linewidth = 1
  ) + # 真のパラメータ
  geom_line(
    data    = posterior_df, 
    mapping = aes(x = lambda, y = dens, linetype = "posterior"), 
    color = "purple", linewidth = 1
  ) + # 事後分布
  scale_x_continuous(
    sec.axis = sec_axis(
      transform = ~ ., 
      breaks    = lambda_truth, 
      labels    = expression(lambda[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 = "Gamma distribution", 
    subtitle = posterior_param_lbl, 
    x = expression(lambda), 
    y = "density"
  )

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

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

予測分布の計算

　観測データ $\mathbf{X}$ から、未観測データ $x_{*}$ の予測分布(負の二項分布) $p(x_{*} \mid \mathbf{X}, a, b) = \mathrm{NB}(x_{*} \mid \hat{r}, \hat{p})$ のパラメータ $\hat{r}, \hat{p}$ を求める。
　負の二項分布については「【R】負の二項分布の作図 - からっぽのしょこ」を参照のこと。

　予測分布のパラメータ $\hat{r}, \hat{p}$ を計算する。

# 事後分布のパラメータにより予測分布のパラメータを計算:式(3.44')
r_hat <- a_hat
q_hat <- 1 / (b_hat + 1)
p_hat <- b_hat / (b_hat + 1)
r_hat; q_hat; p_hat

[1] 430
[1] 0.009803922
[1] 0.9901961

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

[1] 430
[1] 0.009803922
[1] 0.9901961

# 失敗確率パラメータから成功確率パラメータを計算
p_hat <- 1 - q_hat
p_hat

[1] 0.9901961

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

$\displaystyle \begin{aligned} \hat{r} &= \hat{a} \\ &= \sum_{n=1}^N x_n + a \\ \hat{q} &= \frac{1}{1 + \hat{b}} \\ &= \frac{1}{N + 1 + b} \\ \hat{p} &= 1 - \hat{q} \\ &= \frac{\hat{b}}{1 + \hat{b}} \\ &= \frac{N + b}{N + 1 + b} \end{aligned} \tag{3.44'}$

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

# 予測分布の確率を計算:式(3.43)
predict_df <- tibble::tibble(
  x    = x_vec, # 確率変数
  r    = r_hat, # 成功回数パラメータ
  p    = p_hat, # 成功確率パラメータ
  prob = dnbinom(x = x, size = r, prob = p) # 確率
)
predict_df

# A tibble: 16 × 4
       x     r     p      prob
   <dbl> <dbl> <dbl>     <dbl>
 1     0   430 0.990 0.0145   
 2     1   430 0.990 0.0610   
 3     2   430 0.990 0.129    
 4     3   430 0.990 0.182    
 5     4   430 0.990 0.193    
 6     5   430 0.990 0.164    
 7     6   430 0.990 0.117    
 8     7   430 0.990 0.0713   
 9     8   430 0.990 0.0382   
10     9   430 0.990 0.0182   
11    10   430 0.990 0.00784  
12    11   430 0.990 0.00307  
13    12   430 0.990 0.00111  
14    13   430 0.990 0.000369 
15    14   430 0.990 0.000115 
16    15   430 0.990 0.0000332

　 $x$ の値ごとに、負の二項分布に従う確率 $\mathrm{NB}(x \mid \hat{r}, \hat{p})$ を計算する。
　負の二項分布の確率質量関数 dnbinom() の確率変数の引数 k に $x$ 、成功回数の引数 n に $\hat{r}$ 、成功確率の引数 p に $\hat{p}$ を指定する。

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

# 予測分布のラベルを作成
predict_param_lbl <- paste0(
  "list(", 
  "N == ",             N, ", ", 
  "lambda[truth] == ", round(lambda_truth, digits = 2), ", ", 
  "hat(r) == ",        round(r_hat,        digits = 1), ", ", 
  "hat(p) == ",        round(p_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 = "Negative Binomial distribution", 
    subtitle = predict_param_lbl, 
    x = expression(x), 
    y = "probability"
  )

　真の分布(ポアソン分布)を赤色(白抜き)の棒グラフ(破線)、予測分布(負の二項分布)を紫色(塗りつぶし)の棒グラフ(実線)で示す。

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

　以上で、ポアソンモデルのベイズ推論を実装した。

スポンサードリンク

学習の推移

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

データ数と分布の関係

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

　 $n$ 個のデータから求めた( $n$ 回更新した)事後分布(ガンマ分布) $\mathrm{Gam}(\lambda \mid a^{(n)}, b^{(n)})$ を紫色の曲線(実線)、 $n$ 番目の観測データ $x_n$ に対応する位置 $\lambda = 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 + b^{(n-1)} \end{aligned}$

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

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

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

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

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

$\displaystyle \begin{aligned} r^{(n)} &\leftarrow x_n + r^{(n-1)} \\ q^{(n)} &\leftarrow \frac{1}{1 + \frac{1}{q^{(n-1)}}} \\ p^{(n)} &= 1 - q^{(n)} \end{aligned}$

　超パラメータの初期値 $a^{(0)}, b^{(0)}$ を用いて求めたパラメータを $r^{(0)}, p^{(0)}$ として、 $r^{(n-1)}, p^{(n-1)}$ は $n-1$ 回目の更新値を表す。

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

観測データと分布の関係

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

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

　生成分布(ポアソン分布)、事後分布(ガンマ分布)、予測分布(負の二項分布)の期待値は、それぞれ次の式で計算できる。

$\displaystyle \begin{align} \mathbb{E}_{\mathrm{Poi}(x_n \mid \lambda_{\mathrm{truth}})}[x] &= \lambda_{\mathrm{truth}} \\ \mathbb{E}_{\mathrm{Gam}(\lambda \mid \hat{a}, \hat{b})}[\lambda] &= \frac{\hat{a}}{\hat{b}} \\ \mathbb{E}_{\mathrm{NB}(x_{*} \mid \hat{r}, \hat{p})}[x] &= \frac{\hat{r} (1 - \hat{p})}{\hat{p}} \\ &= \frac{\hat{a}}{\hat{b}} \end{align}$

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

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

　この記事では、ポアソン分布に対するベイズ推論を実装した。次の記事では、1次元ガウス分布を扱う。