はじめに

　『ベイズ推論による機械学習入門』の学習時のノートです。「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」によって、「数式」と「プログラム」から理解するのが目標です。
　省略してある内容等ありますので、本とあわせて読んでください。

　この記事では、「尤度関数をポアソン分布」、「事前分布をガンマ分布」とした場合の「パラメータの事後分布」と「未観測値の予測分布」の計算をR言語で実装します。

【数式読解編】

www.anarchive-beta.com

【他の節の内容】

www.anarchive-beta.com

【この節の内容】

はじめに
3.2.3 ポアソン分布の学習と予測
- ベイズ推論の実装
- 学習推移の可視化
参考文献
おわりに

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

　尤度関数をポアソン分布、事前分布をガンマ分布とするモデルに対するベイズ推論を実装します。人工的に生成したデータを用いて、ポアソン分布のパラメータを推定し、また未観測データに対する予測分布を求めます。ポアソン分布については「【R】ポアソン分布の作図 - からっぽのしょこ」、ガンマ分布については「【R】ガンマ分布の作図 - からっぽのしょこ」を参照してください。

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

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

　この記事では、基本的にパッケージ名::関数名()の記法を使うので、パッケージを読み込む必要はありません。ただし、作図コードがごちゃごちゃしないようにパッケージ名を省略しているため、ggplot2は読み込む必要があります。
　magrittrパッケージのパイプ演算子%>%ではなく、ベースパイプ(ネイティブパイプ)演算子|>を使っています。%>%に置き換えても処理できます。
　学習の推移をアニメーション(gif画像)で確認するのにgganimateパッケージを利用します。不要であれば省略してください。

ベイズ推論の実装

　まずは、モデルを設定して、データを生成します。生成したデータを用いて、事後分布のパラメータを計算します。さらに、事後分布のパラメータを用いて、予測分布のパラメータを計算します。

生成分布の設定

　データ生成分布(真の分布)として、ポアソン分布$\mathrm{Poi}(x | \lambda)$を設定します。

　真の分布(ポアソン分布)のパラメータ$\lambda$を設定します。

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

　成功回数の期待値$\lambda > 0$をlambda_truthとして値を指定します。lambda_truthが真のパラメータであり、この値を求めるのがここでの目的です。

　真の分布を計算して、作図用のデータフレームを作成します。

# グラフ用のxの値を作成
x_vec <- seq(0, lambda_truth*3)

# 真の分布を計算:式(2.37)
model_df <- tibble::tibble(
  x = x_vec, # 確率変数
  prob = dpois(x = x_vec, lambda = lambda_truth) # 確率
)
model_df

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

　ポアソン分布の確率変数がとり得る値(非負の整数)$x$を作成して、値ごとに確率を計算します。この例では、0からlambda_truthの3倍を範囲とします。
　ポアソン分布の確率は、dpois()で計算できます。確率変数の引数xにx_vec、パラメータの引数lambdaにlambda_truthを指定します。

　真の分布のグラフを作成します。

# 真の分布を作図
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軸目盛
  labs(title = "Poisson Distribution", 
       subtitle = parse(text = paste0("lambda==", lambda_truth)), 
       x = "x", y = "probability")

　真のパラメータを求めることは、真の分布を求めることを意味します。

データの生成

　続いて、設定した生成分布から観測データ$\mathbf{X} = \{x_1, x_2, \cdots, x_N\}$を生成します。

　ポアソン分布に従う$N$個のデータをランダムに生成します。

# (観測)データ数を指定
N <- 100

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

## [1] 4 8 7 2 4 6

　生成するデータ数$N$を指定します。
　ポアソン分布に従う乱数は、rpois()で生成できます。データ数の引数nにN、パラメータの引数lambdaにlambda_truthを指定します。生成したN個のデータをx_nとします。

　観測データ$\mathbf{X}$の度数を集計します。

# 観測データを集計
freq_df <- tidyr::tibble(x = x_n) |> # 観測データを格納
  dplyr::count(x, name = "freq") # 度数を集計
freq_df

## # A tibble: 12 × 2
##        x  freq
##    <int> <int>
##  1     0     2
##  2     1     7
##  3     2    15
##  4     3    21
##  5     4    19
##  6     5    12
##  7     6    11
##  8     7     5
##  9     8     3
## 10     9     2
## 11    10     2
## 12    11     1

　観測データx_nをデータフレームに格納して、count()で重複をカウントします。
　観測されなかった値は、データフレームに含まれません。

　そこで、観測データにない値も含めた集計結果を作成します。

# ラベル用に度数を整形
freq_vec <- freq_df |> 
  dplyr::right_join(tidyr::tibble(x = x_vec), by = "x") |> # グラフ用の値に結合
  dplyr::mutate(freq = tidyr::replace_na(freq, 0)) |> # 観測にない場合の欠損値を0に置換
  dplyr::arrange(x) |> # 昇順に並べ替え
  dplyr::pull(freq) # ベクトルとして抽出
head(freq_vec)

## [1]  2  7 15 21 19 12

　x_vecの値に対応した度数を求めます。x_vecの範囲外の値は含まれません。
　freq_dfの度数列を、x_vecの値を持つデータフレームにright_join()で結合します。x_nにない値が欠損値NAになるので、replace_na()で0に置換します。
　この処理は、作図自体には影響しません。サブタイトルに度数を表示するのに必要な処理です。

　観測データのヒストグラムを真の分布と重ねて確認します。

# 観測データのヒストグラムを作成
ggplot() + 
  geom_bar(data = freq_df, mapping = aes(x = x, y = freq/N, fill = "data", color = "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 = "") + # 線の色:(凡例表示用)
  guides(color = guide_legend(override.aes = list(size = c(0.5, 0.5), linetype = c(2, 1)))) + # 凡例の体裁:(凡例表示用)
  scale_x_continuous(breaks = 0:max(c(x_vec, x_n)), minor_breaks = FALSE) + # x軸目盛
  labs(title = "Poisson Distribution", 
       subtitle = parse(text = paste0("list(lambda==", lambda_truth, ", N==", N, "(", paste0(freq_vec, collapse = ", "), "))")), 
       x = "x", y = "relative frequency")

　各値の度数(freq列)をデータ数Nで割って、相対度数をバーの高さとします。

　データ数が十分に大きいと、ヒストグラムの形状が真の分布に近付きます。

事前分布の設定

　次は、尤度に対する共役事前分布を設定します。ポアソン分布のパラメータ$\lambda$の事前分布としてガンマ分布$\mathrm{Gam}(\lambda | a, b)$を設定します。

　事前分布(ガンマ分布)のパラメータ(超パラメータ)$a, b$を設定します。

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

　$a > 0, b > 0$の値を指定します。

　事前分布を計算します。

# グラフ用のlambdaの値を作成
lambda_vec <- seq(0, lambda_truth*2, length.out = 501)

# 事前分布を計算:式(2.56)
prior_df <- tibble::tibble(
  lambda = lambda_vec, # 確率変数
  dens = dgamma(x = lambda_vec, shape = a, rate = b) # 確率密度
)
prior_df

## # A tibble: 501 × 2
##    lambda  dens
##     <dbl> <dbl>
##  1  0     1    
##  2  0.024 0.976
##  3  0.048 0.953
##  4  0.072 0.931
##  5  0.096 0.908
##  6  0.12  0.887
##  7  0.144 0.866
##  8  0.168 0.845
##  9  0.192 0.825
## 10  0.216 0.806
## # … with 491 more rows

　ガンマ分布の確率変数がとり得る値$\lambda > 0$を作成して、値ごとに確率密度を計算します。この例では、0からlambda_truthの2倍を範囲とします。
　ガンマ分布の確率密度は、dgamma()で計算できます。確率変数の引数xにlambda_vec、形状の引数shapeにa、尺度の引数rateにbを指定します。

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

# 事前分布を作図
ggplot() + 
  geom_line(data = prior_df, mapping = aes(x = lambda, y = dens, color = "prior"), 
            size = 1) + # 事前分布
  geom_vline(aes(xintercept = lambda_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.5, 0.5), linetype = c(2, 1)))) + # 凡例の体裁:(凡例表示用)
  labs(title = "Gamma Distribution", 
       subtitle = parse(text = paste0("list(a==", a, ", b==", b, ")"),), 
       x = expression(lambda), y = "density")

　真のパラメータの値を赤色の破線で示します。

事後分布の計算

　観測データ$\mathbf{X}$からパラメータ$\lambda$の事後分布を求めます(パラメータ$\lambda$を分布推定します)。

　観測データ$\mathbf{X}$を用いて、事後分布(ガンマ分布)のパラメータ$\hat{a}, \hat{b}$を計算します。

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

## [1] 411
## [1] 101

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

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

　事前分布のパラメータ$a, b$に、それぞれ観測データの総和$\sum_{n=1}^N x_n$とデータ数$N$を加えています。

　事後分布を計算します。

# 事後分布を計算:式(2.56)
posterior_df <- tibble::tibble(
  lambda = lambda_vec, # 確率変数
  dens = dgamma(x = lambda_vec, shape = a_hat, rate = b_hat) # 確率密度
)
posterior_df

## # A tibble: 501 × 2
##    lambda  dens
##     <dbl> <dbl>
##  1  0         0
##  2  0.024     0
##  3  0.048     0
##  4  0.072     0
##  5  0.096     0
##  6  0.12      0
##  7  0.144     0
##  8  0.168     0
##  9  0.192     0
## 10  0.216     0
## # … with 491 more rows

　更新した超パラメータa_hat, b_hatを使って、事前分布のときと同様にして処理します。

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

# 事後分布を作図
ggplot() + 
  geom_line(data = posterior_df, mapping = aes(x = lambda, y = dens, color = "posterior"), 
            size = 1) + # 事後分布
  geom_vline(mapping = aes(xintercept = lambda_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.5, 0.5), linetype = c(2, 1)))) + # 凡例の体裁:(凡例表示用)
  labs(title = "Gamma Distribution", 
       subtitle = parse(text = paste0("list(N==", N, ", hat(a)==", a_hat, ", hat(b)==", b_hat, ")")), 
       x = expression(lambda), y = "density")

　$\lambda$の真の値付近をピークとする分布になっています。

予測分布の計算

　最後に、観測データ$\mathbf{X}$または(観測データから求めた)事後分布のパラメータ$\hat{a}, \hat{b}$から未観測のデータ$x_{*}$の予測分布を求めます。

　事後分布のパラメータ$\hat{a}, \hat{b}$、または観測データ$\mathbf{X}$と事前分布のパラメータ$a, b$を用いて、予測分布(負の二項分布)のパラメータを計算します。

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

## [1] 411
## [1] 0.009803922

　予測分布のパラメータは、次の式で計算できます。

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

　1つ目の式だと、事後分布のパラメータa_hat, b_hatを使って計算できます。2つ目の式だと、観測データx_nと事前分布のパラメータa, bを使って計算できます。

　予測分布を計算します。

# 予測分布を計算:式(3.43)
predict_df <- tibble::tibble(
  x = x_vec, # 確率変数
  prob = dnbinom(x = x_vec, size = r_hat, prob = 1-p_hat) # 確率
)
predict_df

## # A tibble: 13 × 2
##        x     prob
##    <int>    <dbl>
##  1     0 0.0174  
##  2     1 0.0703  
##  3     2 0.142   
##  4     3 0.191   
##  5     4 0.194   
##  6     5 0.158   
##  7     6 0.107   
##  8     7 0.0628  
##  9     8 0.0322  
## 10     9 0.0147  
## 11    10 0.00604 
## 12    11 0.00227 
## 13    12 0.000782

　負の二項分布の確率は、dnbinom()で計算できます。失敗回数の引数xにx_vec、成功回数の引数sizeにr_hat、成功確率の引数probに1-p_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 = "model", 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 = "") + # 線の色:(凡例表示用)
  guides(color = guide_legend(override.aes = list(size = c(0.5, 0.5), linetype = c(2, 1)))) + # 凡例の体裁:(凡例表示用)
  scale_x_continuous(breaks = x_vec, minor_breaks = FALSE) + # x軸目盛
  labs(title = "Negative Binomial Distribution", 
       subtitle = parse(text = paste0("list(N==", N, ", hat(r)==", r_hat, ", hat(p)==", round(p_hat, 3), ")")), 
       x = "x", y = "probability")

　観測データが増えると、予測分布が真の分布に近付きます。

　以上で、ポアソン分布のベイズ推論を実装できました。

学習推移の可視化

　gganimateパッケージを利用して、パラメータの推定値(更新値)の推移(分布の変化)をアニメーション(gif画像)で確認します。

モデルの設定

　真の分布のパラメータ$\lambda$と、事前分布のパラメータ(の初期値)$a, b$を設定します。

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

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

# データ数(試行回数)を指定
N <- 150

# グラフ用のlambdaの値を作成
lambda_vec <- seq(0, lambda_truth*3, length.out = 501)

# グラフ用のxの値を作成
x_vec <- seq(0, lambda_truth*3)

　実装時と同じく、パラメータを指定して、確率変数がとり得る値(x軸の値)を作成します。グラフが粗い場合や処理が重い場合は、lambda_vec, x_vecの要素数(値の間隔)を調整してください。

　事後分布と予測分布のパラメータの計算(更新)処理について、2通りの方法を載せます。目的に応じて使い分けてください。

推論処理：for関数による処理

　1つ目の処理方法では、for()を使って、1データずつ生成してパラメータの更新を繰り返し処理します。こちらの方が、前ステップで求めた事後分布(のパラメータ)を次ステップの事前分布(のパラメータ)として用いる逐次学習をイメージしやすいです。

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

　超パラメータの初期値を使って、事前分布を計算します。

# 事前分布(ガンマ分布)を計算
anime_posterior_df <- tibble::tibble(
  lambda = lambda_vec, # 確率変数
  dens = dgamma(x = lambda_vec, shape = a, rate = b), # 確率密度
  param = paste0("N=", 0, "=(", paste0(rep(0, times = length(x_vec)), collapse = ", "), "), a=", a, ", b=", b) |> 
    as.factor() # フレーム切替用ラベル
)
anime_posterior_df

## # A tibble: 501 × 3
##    lambda  dens param                                                
##     <dbl> <dbl> <fct>                                                
##  1  0     1     N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), a=1, b=1
##  2  0.024 0.976 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), a=1, b=1
##  3  0.048 0.953 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), a=1, b=1
##  4  0.072 0.931 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), a=1, b=1
##  5  0.096 0.908 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), a=1, b=1
##  6  0.12  0.887 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), a=1, b=1
##  7  0.144 0.866 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), a=1, b=1
##  8  0.168 0.845 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), a=1, b=1
##  9  0.192 0.825 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), a=1, b=1
## 10  0.216 0.806 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), a=1, b=1
## # … with 491 more rows

　現在のパラメータを文字列結合し因子型に変換して、フレーム切替用のラベル列とします。

　また、事前分布のパラメータを使って、予測分布のパラメータと予測分布を計算します。

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

# 初期値による予測分布(負の二項分布)を計算
anime_predict_df <- tibble::tibble(
  x = x_vec, # 確率変数
  prob = dnbinom(x = x_vec, size = r, prob = 1-p), # 確率
  param = paste0("N=", 0, "=(", paste0(rep(0, times = length(x_vec)), collapse = ", "), "), r=", r, ", p=", round(p, 5)) |> 
    as.factor() # フレーム切替用ラベル
)
anime_predict_df

## # A tibble: 13 × 3
##        x     prob param                                                  
##    <int>    <dbl> <fct>                                                  
##  1     0 0.5      N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
##  2     1 0.25     N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
##  3     2 0.125    N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
##  4     3 0.0625   N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
##  5     4 0.0312   N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
##  6     5 0.0156   N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
##  7     6 0.00781  N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
##  8     7 0.00391  N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
##  9     8 0.00195  N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
## 10     9 0.000977 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
## 11    10 0.000488 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
## 12    11 0.000244 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
## 13    12 0.000122 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5

　同様に、フレーム切替用のラベル列を作成します。

　パラメータの更新処理をN回繰り返します。

# 観測データの受け皿を作成
x_n <- rep(NA, times = N)

# ベイズ推論
for(n in 1:N){
  
  # ポアソン分布に従うデータを生成
  x_n[n] <- rpois(n = 1 ,lambda = lambda_truth)
  
  # 観測データを集計
  freq_vec <- tidyr::tibble(x = x_n[1:n]) |> # 観測データを格納
    dplyr::count(x, name = "freq") |> # 度数を集計
    dplyr::right_join(tidyr::tibble(x = x_vec), by = "x") |> # グラフ用の値に結合
    dplyr::mutate(freq = tidyr::replace_na(freq, 0)) |> # 観測にない場合の欠損値を0に置換
    dplyr::arrange(x) |> # 昇順に並べ替え
    dplyr::pull(freq) # ベクトルとして抽出
  
  # 事後分布のパラメータを更新:式(3.38)
  a <- sum(x_n[n] * 1) + a
  b <- 1 + b
  
  # 事後分布(ガンマ分布)を計算:式(2.56)
  tmp_posterior_df <- tibble::tibble(
    lambda = lambda_vec, # 確率変数
    dens = dgamma(x = lambda_vec, shape = a, rate = b), # 確率密度
    param = paste0("N=", n, "=(", paste0(freq_vec, collapse = ", "), "), a=", a, ", b=", b) |> 
      as.factor() # フレーム切替用ラベル
  )
  
  # 予測分布のパラメータを更新:式(3.44)
  r <- a
  p <- 1 / (b + 1)
  
  # 予測分布(負の二項分布)を計算:式(3.43)
  tmp_predict_df <- tibble::tibble(
    x = x_vec, # 確率変数
    prob = dnbinom(x = x_vec, size = r, prob = 1-p), # 確率
    param = paste0("N=", n, "=(", paste0(freq_vec, collapse = ", "), "), r=", r, ", p=", round(p, 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.38)の$\sum_{n=1}^N x_n$と$N$の計算は、ループ処理によってN回繰り返しx_n[n]と1を加えることで行います。n回目のループ処理のときには、n-1回分のx_n[n]と1が既にaとbに加えられています。
　更新したパラメータを使って事後分布と予測分布を計算して、それぞれ計算結果をanime_***_dfに結合していきます。

　結果を確認します。

# 確認
head(x_n); anime_posterior_df; anime_predict_df

## [1] 4 8 7 2 4 6
## # A tibble: 50,601 × 3
##    lambda  dens param                                                
##     <dbl> <dbl> <fct>                                                
##  1  0     1     N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), a=1, b=1
##  2  0.024 0.976 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), a=1, b=1
##  3  0.048 0.953 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), a=1, b=1
##  4  0.072 0.931 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), a=1, b=1
##  5  0.096 0.908 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), a=1, b=1
##  6  0.12  0.887 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), a=1, b=1
##  7  0.144 0.866 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), a=1, b=1
##  8  0.168 0.845 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), a=1, b=1
##  9  0.192 0.825 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), a=1, b=1
## 10  0.216 0.806 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), a=1, b=1
## # … with 50,591 more rows
## # A tibble: 1,313 × 3
##        x     prob param                                                  
##    <int>    <dbl> <fct>                                                  
##  1     0 0.5      N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
##  2     1 0.25     N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
##  3     2 0.125    N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
##  4     3 0.0625   N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
##  5     4 0.0312   N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
##  6     5 0.0156   N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
##  7     6 0.00781  N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
##  8     7 0.00391  N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
##  9     8 0.00195  N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
## 10     9 0.000977 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
## # … with 1,303 more rows

　それぞれ「lambda_vec, x_vecの要素数」掛ける「N+1」行のデータフレームになります。行数が増えすぎないように注意してくださいしてください。

推論処理：tidyverseパッケージによる処理

　2つ目の処理方法では、tidyverseパッケージの関数を使って、一度の処理でN+1回分の計算をします。こちらの方が、処理時間が短いです。

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

　ポアソン分布に従う$N$個のデータを生成します。

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

## [1] 4 8 7 2 4 6

　ラベル用のテキストとして、観測データを集計します。

# 試行ごとに度数を集計:(ラベル用)
freq_vec <- tibble::tibble(
  x = c(NA, x_n), # 観測データ
  n = 0:N, # 試行回数
  freq = 1 # 集計用の値
) |> # 観測データを格納
  dplyr::right_join(tidyr::expand_grid(x = x_vec, n = 0:N), by = c("x", "n")) |> # グラフ用の値に結合
  dplyr::mutate(freq = tidyr::replace_na(freq, replace = 0)) |> # 観測にない場合の欠損値を0に置換
  dplyr::arrange(n, x) |> # 集計用に昇順に並べ替え
  dplyr::group_by(x) |> # 集計用にグループ化
  dplyr::mutate(freq = cumsum(freq)) |> # 各試行までの度数を計算
  dplyr::ungroup() |> # グループ化を解除
  tidyr::pivot_wider(
    id_cols = n, 
    names_from = x, 
    names_prefix = "x", 
    values_from = freq
  ) |> # 度数列を展開
  tidyr::unite(col = "freq", dplyr::starts_with("x"), sep = ", ") |> # 度数情報をまとめて文字列化
  dplyr::pull(freq) # ベクトルとして抽出
head(freq_vec)

## [1] "0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0"
## [2] "0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0"
## [3] "0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0"
## [4] "0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0"
## [5] "0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0"
## [6] "0, 0, 1, 0, 2, 0, 0, 1, 1, 0, 0, 0, 0"

　試行ごとに、それまでの試行における各値の度数を求めて、文字列結合したベクトルを作成します。

　(0回目に対応する欠損値を含めた)観測データx_n、試行回数(データ番号)0:N、度数計算用の値1をデータフレームに格納します。
　$x$の値x_vecとデータ番号の全ての組み合わせをexpand_grid()で作成して、そこに観測データをright_join()で結合します。x_nにない値はfreq列が欠損値NAになるので、replace_na()で0に置換します。
　x列($x$の値)でグループ化して、freq列(各試行における観測値が1でそれ以外は0)の累積和をcumsum()で計算して、各試行回数までの度数を求めます。

　pivot_longer()で、試行回数列と、x_vecに対応した度数列に展開します。
　度数列の値をunite()で文字列結合して、pull()でベクトルとして取り出します。

　事前分布を含めたN+1回分の事後分布を計算します。

# 試行ごとに事後分布(ガンマ分布)を計算
anime_posterior_df <- tidyr::expand_grid(
  n = 0:N, # 試行回数
  lambda = lambda_vec # 確率変数
) |> # 全ての組み合わせを作成
  dplyr::mutate(
    a = c(a, cumsum(x_n) + a)[n+1], 
    b = n + b
  ) |> # 事後分布のパラメータを計算:式(3.38)
  dplyr::mutate(
    dens = dgamma(x = lambda, shape = a, rate = b), # 確率密度
    param = paste0("N=", n, "=(", freq_vec[n+1], "), a=", a, ", b=", b) |> 
      (\(.){factor(., levels = unique(.))})() # フレーム切替用ラベル
  ) # 事後分布を計算:式(2.56)

## # A tibble: 1,313 × 6
##        n     x     r     p     prob param                                       
##    <int> <int> <dbl> <dbl>    <dbl> <fct>                                       
##  1     0     0     1   0.5 0.5      N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)…
##  2     0     1     1   0.5 0.25     N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)…
##  3     0     2     1   0.5 0.125    N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)…
##  4     0     3     1   0.5 0.0625   N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)…
##  5     0     4     1   0.5 0.0312   N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)…
##  6     0     5     1   0.5 0.0156   N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)…
##  7     0     6     1   0.5 0.00781  N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)…
##  8     0     7     1   0.5 0.00391  N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)…
##  9     0     8     1   0.5 0.00195  N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)…
## 10     0     9     1   0.5 0.000977 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)…
## # … with 1,303 more rows

　試行回数(観測データ数)を表す0からNまでのN+1個の整数と、x軸の値lambda_vecの全ての組み合わせをexpand_grid()で作成します。これにより、lambda_vecの各要素をN+1フレーム分に複製できます。
　事前分布とN回分の事後分布のパラメータa(のベクトル)を計算(作成)して、n列の値をインデックスとして使って各試行に対応する値を抽出します。bについてはn列を使って計算できます。
　確率変数とパラメータの組み合わせごとに確率密度を計算して、パラメータごとにフレーム切替用のラベルを作成します。
　因子型への変換処理では、無名関数function()の省略記法\()を使って、(\(x){引数xを使った具体的な処理})()としています。直前のパイプ演算子を%>%にすると、行全体(\(引数){処理})()を{}の中の処理に置き換えられます(置き換えられるように引数名を.にしています)。

　同様に、初期値を含めたN+1回分の予測分布を計算します。

# 試行ごとに予測分布(負の二項分布)を計算
anime_predict_df <- tidyr::expand_grid(
  n = 0:N, # 試行回数
  x = x_vec # 確率変数
) |> 
  dplyr::mutate(
    r = c(a, cumsum(x_n) + a)[n+1], 
    p = 1 / (n + b + 1)
  ) |> # 予測分布のパラメータを計算算:式(3.44)
  dplyr::mutate(
    prob = dnbinom(x = x, size = r, prob = 1-p), # 確率
    param = paste0("N=", n, "=(", freq_vec[n+1], "), r=", r, ", p=", round(p, 3)) |> 
      (\(.){factor(., levels = unique(.))})() # フレーム切替用ラベル
  ) # 予測分布を計算:式(3.43)
anime_predict_df

## # A tibble: 101 × 2
##        x param                                                  
##    <int> <fct>                                                  
##  1    NA N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), a=1, b=1  
##  2     4 N=1=(0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0), a=5, b=2  
##  3     8 N=2=(0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0), a=13, b=3 
##  4     7 N=3=(0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0), a=20, b=4 
##  5     2 N=4=(0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0), a=22, b=5 
##  6     4 N=5=(0, 0, 1, 0, 2, 0, 0, 1, 1, 0, 0, 0, 0), a=26, b=6 
##  7     6 N=6=(0, 0, 1, 0, 2, 0, 1, 1, 1, 0, 0, 0, 0), a=32, b=7 
##  8     2 N=7=(0, 0, 2, 0, 2, 0, 1, 1, 1, 0, 0, 0, 0), a=34, b=8 
##  9     3 N=8=(0, 0, 2, 1, 2, 0, 1, 1, 1, 0, 0, 0, 0), a=37, b=9 
## 10     5 N=9=(0, 0, 2, 1, 2, 1, 1, 1, 1, 0, 0, 0, 0), a=42, b=10
## # … with 91 more rows

　0からNまでの整数と、x軸の値x_vecの全ての組み合わせを作成します。
　N+1回分の予測分布パラメータr(のベクトル)を計算(作成)して、各試行に対応する値を抽出します。pについてはn列を使って計算できます。
　確率変数とパラメータの組み合わせごとに確率を計算して、パラメータごとにフレーム切替用のラベルを作成します。

作図処理

　事後分布と予測分布の推移をそれぞれアニメーションで可視化します。

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

　観測データと対応するラベルをデータフレームに格納します。

# 観測データを格納
anime_data_df <- tibble::tibble(
  x = c(NA, x_n), 
  param = unique(anime_posterior_df[["param"]]) # フレーム切替用ラベル
)
anime_data_df

## # A tibble: 1,313 × 3
##        x   prob param                                                  
##    <int>  <dbl> <fct>                                                  
##  1     0 0.0183 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
##  2     1 0.0733 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
##  3     2 0.147  N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
##  4     3 0.195  N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
##  5     4 0.195  N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
##  6     5 0.156  N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
##  7     6 0.104  N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
##  8     7 0.0595 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
##  9     8 0.0298 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
## 10     9 0.0132 N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5
## # … with 1,303 more rows

　事前分布には観測データが影響しないので欠損値NAとして、anime_posterior_dfのラベル列と合わせて格納します。
　参考として、各試行における観測データを表示するのに使います。

　事後分布の推移のアニメーションを作成します。

# 事後分布のアニメーションを作図
posterior_graph <- ggplot() + 
  geom_vline(mapping = aes(xintercept = lambda_truth, color = "param"), 
             size = 1, linetype = "dashed", show.legend = FALSE) + # 真のパラメータ
  geom_line(data = anime_posterior_df, mapping = aes(x = lambda, 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 = "Gamma Distribution", 
       subtitle = "{current_frame}", 
       x = expression(lambda), 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), # 確率変数
  prob = rep(dpois(x = x_vec, lambda = lambda_truth), times = N+1), # 確率
  param = anime_predict_df[["param"]] # フレーム切替用ラベル
)
anime_model_df

　棒グラフの場合は、フレームごとにデータを用意する必要があります。
　真の分布は全てのフレームで変化しないので、N+1個に複製します。

　観測データと対応するラベルをデータフレームに格納します。

# 観測データを格納
anime_data_df <- tibble::tibble(
  x = c(NA, x_n), 
  param = unique(anime_predict_df[["param"]]) # フレーム切替用ラベル
)
anime_data_df

## # A tibble: 101 × 2
##        x param                                                     
##    <int> <fct>                                                     
##  1    NA N=0=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), r=1, p=0.5   
##  2     4 N=1=(0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0), r=5, p=0.333 
##  3     8 N=2=(0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0), r=13, p=0.25 
##  4     7 N=3=(0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0), r=20, p=0.2  
##  5     2 N=4=(0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0), r=22, p=0.167
##  6     4 N=5=(0, 0, 1, 0, 2, 0, 0, 1, 1, 0, 0, 0, 0), r=26, p=0.143
##  7     6 N=6=(0, 0, 1, 0, 2, 0, 1, 1, 1, 0, 0, 0, 0), r=32, p=0.125
##  8     2 N=7=(0, 0, 2, 0, 2, 0, 1, 1, 1, 0, 0, 0, 0), r=34, p=0.111
##  9     3 N=8=(0, 0, 2, 1, 2, 0, 1, 1, 1, 0, 0, 0, 0), r=37, p=0.1  
## 10     5 N=9=(0, 0, 2, 1, 2, 1, 1, 1, 1, 0, 0, 0, 0), r=42, p=0.091
## # … with 91 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軸目盛
  coord_cartesian(ylim = c(0, 0.5)) + # 軸の表示範囲
  labs(title = "Negative Binomial Distribution", 
       subtitle = "{current_frame}", 
       x = "x", y = "probability")

# gif画像を作成
gganimate::animate(predict_graph, nframes = N+1, fps = 10, width = 800, height = 600)

　データが増えるにつれて、真のパラメータ付近の確率密度が大きくなっていくのを確認できます。

　予測分布が真の分布に近付いていくのを確認できます。

参考文献

須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.

おわりに

　加筆修正の際に記事を分割しました。もう少し修正のペースを上げたい。

2022.07.21：加筆修正しました。

【次節の内容】

からっぽのしょこ

読んだら書く！書いたら読む！同じ事は二度調べ(たく)ない

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

はじめに

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

ベイズ推論の実装

生成分布の設定

データの生成

事前分布の設定

事後分布の計算

予測分布の計算

学習推移の可視化

モデルの設定

推論処理：for関数による処理

推論処理：tidyverseパッケージによる処理

作図処理

参考文献

おわりに