はじめに

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

　この記事では、平均と精度が未知の1次元ガウス分布に対するベイズ推論をR言語でスクラッチ実装します。

【前節の内容】

www.anarchive-beta.com

【他の節の内容】

www.anarchive-beta.com

【この節の内容】

はじめに
3.3.3 1次元ガウス分布のベイズ推論の実装：平均と精度が未知の場合
- ベイズ推論の実装
- 学習の推移
  - データ数と分布の関係
  - 観測データと分布の関係
参考文献
おわりに

3.3.3 1次元ガウス分布のベイズ推論の実装：平均と精度が未知の場合

　1次元ガウスモデル(Gaussian model)に対するベイズ推論(Bayesian inference)を実装して、人工的に生成したデータを用いて、パラメータの学習と未知変数の予測を行う。この記事では、生成分布の平均パラメータ(mean parameter)と精度パラメータ(precision parameter)が未知の場合を扱う。平均と精度が未知の1次元ガウスモデルでは、尤度関数を1次元ガウス分布(Gaussian distribution・一変量正規分布・Normal distribution)、事前分布をガウス-ガンマ分布(Gaussian-Gamma distribution・正規-ガンマ分布・Normal-Gamma distribution)とする。この記事では、Rを利用して実装する。
　1次元ガウスモデルについては「3.3.0：1次元ガウスモデルの生成モデルの導出【緑ベイズ入門のノート】 - からっぽのしょこ」、ベイズ推論については「3.3.3：1次元ガウス分布のベイズ推論の導出：平均・精度が未知の場合【緑ベイズ入門のノート】 - からっぽのしょこ」、Pythonを利用する場合は「【Python】3.3.3：1次元ガウス分布のベイズ推論の実装：平均・精度が未知の場合【緑ベイズ入門のノート】 - からっぽのしょこ」を参照のこと。

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

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

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

ベイズ推論の実装

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

生成分布の設定

　データの生成分布(真の分布・ガウス分布) $p(x_n \mid \mu_{\mathrm{truth}}, \lambda_{\mathrm{truth}}) = \mathcal{N}(x_n \mid \mu_{\mathrm{truth}}, \lambda_{\mathrm{truth}}^{-1})$ のパラメータ(真のパラメータ) $\mu_{\mathrm{truth}}, \lambda_{\mathrm{truth}}$ を設定する。
　ガウス分布については「【R】1次元ガウス分布の作図 - からっぽのしょこ」を参照のこと。

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

# 真のパラメータを指定
mu_truth     <- 5
lambda_truth <- 0.25
mu_truth; lambda_truth; 1/sqrt(lambda_truth)

[1] 5
[1] 0.25
[1] 2

　ガウス分布の平均パラメータ(実数) $\mu_{\mathrm{truth}}$ 、精度パラメータ(正の値) $\lambda_{\mathrm{truth}} \gt 0$ を指定する。
　 $\mu_{\mathrm{truth}}, \lambda_{\mathrm{truth}}$ がガウスモデルにおける真のパラメータであり、この値を求めるのがここでの目的(学習)である。

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

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

# x軸の値を作成
x_vec <- seq(from = x_min, to = x_max, length.out = 1001)
x_min; x_max; head(x_vec)

[1] -5
[1] 15
[1] -5.00 -4.98 -4.96 -4.94 -4.92 -4.90

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

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

# 生成分布の確率密度を計算:式(2.64)
model_df <- tibble::tibble(
  x      = x_vec, # 確率変数
  mu     = mu_truth,       # 平均パラメータ
  lambda = lambda_truth,   # 精度パラメータ
  sigma  = 1/sqrt(lambda), # 標準偏差パラメータ
  dens   = dnorm(x = x, mean = mu, sd = sigma) # 確率密度
)
model_df

# A tibble: 1,001 × 5
       x    mu lambda sigma        dens
   <dbl> <dbl>  <dbl> <dbl>       <dbl>
 1 -5        5   0.25     2 0.000000743
 2 -4.98     5   0.25     2 0.000000781
 3 -4.96     5   0.25     2 0.000000821
 4 -4.94     5   0.25     2 0.000000863
 5 -4.92     5   0.25     2 0.000000907
 6 -4.9      5   0.25     2 0.000000953
 7 -4.88     5   0.25     2 0.00000100 
 8 -4.86     5   0.25     2 0.00000105 
 9 -4.84     5   0.25     2 0.00000111 
10 -4.82     5   0.25     2 0.00000116 
# ℹ 991 more rows

　 $x$ の値ごとに、ガウス分布に従う確率密度 $\mathcal{N}(x \mid \mu_{\mathrm{truth}}, \lambda_{\mathrm{truth}}^{-1})$ を計算する。
　ガウス分布の確率密度関数 dnorm() の確率変数の引数 x に $x$ 、平均の引数 mean に $\mu_{\mathrm{truth}}$ 、標準偏差の引数 sd に $\sigma_{\mathrm{truth}} = \frac{1}{\sqrt{\lambda_{\mathrm{truth}}}}$ を指定する。

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

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

# 生成分布を作図
ggplot() + 
  geom_line(
    data    = model_df, 
    mapping = aes(x = x, y = dens, color = "model"), 
    linewidth = 1
  ) + # 生成分布
  scale_color_manual(
    breaks = "model", 
    values = "red", 
    labels = "true model (generator)", 
    name   = ""
  ) + # (凡例表示用)
  guides(
    color = 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 = "Gaussian distribution", 
    subtitle = model_param_lbl, 
    x = expression(x), 
    y = "density"
  )

　真の分布(ガウス分布)を赤色の曲線で示す。

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

データの生成

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

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

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

# 観測データを生成
x_n <- rnorm(n = N , mean = mu_truth, sd = 1/sqrt(lambda_truth))
head(x_n)

[1] 6.417591 7.590121 4.821934 4.140894 6.268684 3.460541

　データ数 $N$ を指定して、ガウス分布に従う乱数 $\mathbf{X} = \{x_1, \cdots, x_N\}$ を生成する。
　ガウス分布の乱数生成関数 rnorm() のサンプルサイズの引数 n に $N$ 、平均の引数 mean に $\mu_{\mathrm{truth}}$ 、標準偏差の引数 sd に $\sigma_{\mathrm{truth}} = \frac{1}{\sqrt{\lambda_{\mathrm{truth}}}}$ を指定する。

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

# 階級数を指定
bin_num <- 40

# 階級幅を計算
bin_size <- (x_max - x_min) / bin_num

# 境界値の範囲を設定
bin_min <- x_min - 0.5*bin_size
bin_max <- x_max + 0.5*bin_size

# 観測データを集計
obs_df <- tibble::tibble(
  x = x_n # サンプル値
) |> 
  dplyr::mutate(
    bin_i  = (x - bin_min) %/% bin_size,        # 階級番号
    center = bin_min + (bin_i + 0.5) * bin_size # 階級値
  ) |> 
  dplyr::count(
    center, name = "freq" # 度数
  ) |> 
  dplyr::mutate(
    dens = freq / (bin_size * N) # 密度
  ) |> 
  tidyr::complete(
    center = seq(from = x_min, to = x_max, by = bin_size), 
    fill = list(freq = 0, dens = 0)
  ) # 未観測値を補完
obs_df

# A tibble: 41 × 3
   center  freq  dens
    <dbl> <int> <dbl>
 1   -5       0     0
 2   -4.5     0     0
 3   -4       0     0
 4   -3.5     0     0
 5   -3       0     0
 6   -2.5     0     0
 7   -2       0     0
 8   -1.5     0     0
 9   -1       0     0
10   -0.5     0     0
# ℹ 31 more rows

　階級数を指定して、階級幅 $w$ を設定する。
　 $x$ の値ごとに、観測データ x_n に含まれる要素数をカウントして、度数 $N_x$ と密度 $\frac{N_x}{w N}$ を求める。

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

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

# 観測データを作図
ggplot() + 
  geom_line(
    data    = model_df, 
    mapping = aes(x = x, y = dens, linetype = "model"), 
    color = "red", linewidth = 1
  ) + # 生成分布
  geom_bar(
    data    = obs_df, 
    mapping = aes(x = center, y = dens, linetype = "data"), 
    stat = "identity", position = "identity", width = bin_size, 
    fill = "hotpink", color = NA, alpha = 0.5
  ) + # 観測データ
  scale_y_continuous(
    sec.axis = sec_axis(transform = ~ .*bin_size*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 = "Gaussian distribution", 
    subtitle = obs_param_lbl, 
    x = expression(x), 
    y = "density"
  )

　 $x$ の値ごとの確率密度(生成分布)を赤色の曲線(破線)、観測データ $\mathbf{X}$ の度数 $N_x$ (密度 $\frac{N_x}{w N}$ )を桃色の棒グラフ(実線)で示す。

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

事前分布の設定

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

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

# μの事前分布のパラメータを指定
m    <- 0
beta <- 1

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

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

　ガウス分布の平均パラメータ(実数) $m$ 、精度パラメータの係数(正の値) $\beta \gt 0$ と、ガンマ分布の形状パラメータ(正の値) $a \gt 0$ 、尺度パラメータ(正の値) $b \gt 0$ を指定する。

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

# μ軸の範囲を設定
mu_min <- x_min
mu_max <- x_max

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

[1] -5
[1] 15
[1] -5.0 -4.8 -4.6 -4.4 -4.2 -4.0

　この例では、生成分布の確率変数 $x$ の範囲を設定している。

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

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

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

[1] 0
[1] 1
[1] 0.00 0.01 0.02 0.03 0.04 0.05

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

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

# 事前分布の確率密度を計算:式(3.81)
prior_df <- tidyr::expand_grid(
  mu     = mu_vec,    # μの確率変数
  lambda = lambda_vec # λの確率変数
) |> # 格子点を作成
  dplyr::mutate(
    m         = m,                 # 平均パラメータ
    beta      = beta,              # 係数パラメータ
    lambda_mu = beta * lambda,     # 精度パラメータ
    sigma_mu  = 1/sqrt(lambda_mu), # 標準偏差パラメータ
    a         = a,                 # 形状パラメータ
    b         = b,                 # 尺度パラメータ
    N_dens    = dnorm(x = mu, mean = m, sd = sigma_mu),  # μの確率密度
    Gam_dens  = dgamma(x = lambda, shape = a, rate = b), # λの確率密度
    dens      = N_dens * Gam_dens # μ, λの確率密度
  )
prior_df

# A tibble: 10,201 × 11
      mu lambda     m  beta lambda_mu sigma_mu     a     b N_dens Gam_dens
   <dbl>  <dbl> <dbl> <dbl>     <dbl>    <dbl> <dbl> <dbl>  <dbl>    <dbl>
 1    -5   0        0     1      0      Inf        1     1 0         1    
 2    -5   0.01     0     1      0.01    10        1     1 0.0352    0.990
 3    -5   0.02     0     1      0.02     7.07     1     1 0.0439    0.980
 4    -5   0.03     0     1      0.03     5.77     1     1 0.0475    0.970
 5    -5   0.04     0     1      0.04     5        1     1 0.0484    0.961
 6    -5   0.05     0     1      0.05     4.47     1     1 0.0477    0.951
 7    -5   0.06     0     1      0.06     4.08     1     1 0.0462    0.942
 8    -5   0.07     0     1      0.07     3.78     1     1 0.0440    0.932
 9    -5   0.08     0     1      0.08     3.54     1     1 0.0415    0.923
10    -5   0.09     0     1      0.09     3.33     1     1 0.0389    0.914
# ℹ 10,191 more rows
# ℹ 1 more variable: dens <dbl>

　 $\mu, \lambda$ の格子状の点( mu_vec, lambda_vec の全ての組み合わせ)を expand_grid() で作成する。
　 $\mu, \lambda$ の値の組み合わせごとに、ガウス-ガンマ分布に従う確率密度 $\mathrm{NG}(\mu, \lambda \mid m, \beta, a, b) = \mathcal{N}(\mu \mid m, (\beta \lambda)^{-1}) \mathrm{Gam}(\lambda \mid a, b)$ を計算する。
　ガウス分布の確率密度関数 dnorm() の確率変数の引数 x に $\mu$ 、平均の引数 mean に $m$ 、標準偏差の引数 sd に $\sigma_{\mu} = \frac{1}{\sqrt{\beta \lambda}}$ を指定する。
　ガンマ分布の確率密度関数 dgamma() の確率変数の引数 x に $\lambda$ 、形状の引数 a に $a$ 、尺度の引数 scale に $b$ を指定する。
　 $\mu$ の値ごとのガウス分布の確率密度 $\mathcal{N}(\mu \mid m, (\beta \lambda)^{-1})$ と、 $\lambda$ の値ごとのガンマ分布の確率密度 $\mathrm{Gam}(\lambda \mid a, b)$ の積を計算する。

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

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

# 事前分布を作図
ggplot() + 
  geom_contour_filled(
    data    = prior_df, 
    mapping = aes(x = mu, y = lambda, z = dens, fill = after_stat(level), linetype = "prior"), 
    alpha = 0.5
  ) + # 事前分布
  geom_vline(
    mapping = aes(xintercept = mu_truth, linetype = "model"), 
    color = "red", linewidth = 1
  ) + # 真のパラメータ
  geom_hline(
    mapping = aes(yintercept = lambda_truth, linetype = "model"), 
    color = "red", linewidth = 1
  ) + # 真のパラメータ
  scale_x_continuous(
    sec.axis = sec_axis(
      transform = ~ ., 
      breaks    = mu_truth, 
      labels    = expression(mu[truth])
    ) # パラメータラベル
  ) + 
  scale_y_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), order = 1), 
    fill     = guide_legend(order = 2)
  ) + 
  labs(
    title = "Gaussian-Gamma distribution", 
    subtitle = prior_param_lbl, 
    fill = "density", 
    x = expression(mu), 
    y = expression(lambda)
  )

　真のパラメータを赤色の直線(破線)、事前分布(ガウス-ガンマ分布)を塗りつぶしの等高線で示す。

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

事後分布の計算

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

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

# μの事後分布のパラメータを計算:式(3.83)
beta_hat <- N + beta
m_hat    <- (sum(x_n) + beta * m) / beta_hat
m_hat; beta_hat

[1] 5.222297
[1] 101

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

$\displaystyle \begin{aligned} \hat{m} &= \frac{1}{\hat{\beta}} \left( \sum_{n=1}^N x_n + \beta m \right) \\ \hat{\beta} &= N + \beta \end{aligned} \tag{3.83}$

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

# λの事後分布のパラメータを計算:式(3.88)
a_hat <- 0.5 * N + a
b_hat <- 0.5 * (sum(x_n^2) + beta * m^2 - beta_hat * m_hat^2) + b
a_hat; b_hat

[1] 51
[1] 181.0194

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

$\displaystyle \begin{aligned} \hat{a} &= \frac{N}{2} + a \\ \hat{b} &= \frac{1}{2} \left( \sum_{n=1}^N x_n^2 + \beta m^2 - \hat{\beta} \hat{m}^2 \right) + b \end{aligned} \tag{3.88}$

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

# 事前分布の確率密度を計算:式(3.81)
posterior_df <- tidyr::expand_grid(
  mu     = mu_vec,    # μの確率変数
  lambda = lambda_vec # λの確率変数
) |> # 格子点を作成
  dplyr::mutate(
    m         = m_hat,             # 平均パラメータ
    beta      = beta_hat,          # 係数パラメータ
    lambda_mu = beta * lambda,     # 精度パラメータ
    sigma_mu  = 1/sqrt(lambda_mu), # 標準偏差パラメータ
    a         = a_hat,             # 形状パラメータ
    b         = b_hat,             # 尺度パラメータ
    N_dens    = dnorm(x = mu, mean = m, sd = sigma_mu),  # μの確率密度
    Gam_dens  = dgamma(x = lambda, shape = a, rate = b), # λの確率密度
    dens      = N_dens * Gam_dens # μ, λの確率密度
  )
posterior_df

# A tibble: 10,201 × 11
      mu lambda     m  beta lambda_mu sigma_mu     a     b    N_dens Gam_dens
   <dbl>  <dbl> <dbl> <dbl>     <dbl>    <dbl> <dbl> <dbl>     <dbl>    <dbl>
 1    -5   0     5.22   101      0     Inf        51  181. 0         0       
 2    -5   0.01  5.22   101      1.01    0.995    51  181. 4.84e- 24 7.49e-51
 3    -5   0.02  5.22   101      2.02    0.704    51  181. 8.28e- 47 1.38e-36
 4    -5   0.03  5.22   101      3.03    0.574    51  181. 1.23e- 69 1.44e-28
 5    -5   0.04  5.22   101      4.04    0.498    51  181. 1.71e- 92 4.16e-23
 6    -5   0.05  5.22   101      5.05    0.445    51  181. 2.31e-115 4.77e-19
 7    -5   0.06  5.22   101      6.06    0.406    51  181. 3.06e-138 7.10e-16
 8    -5   0.07  5.22   101      7.07    0.376    51  181. 3.99e-161 2.59e-13
 9    -5   0.08  5.22   101      8.08    0.352    51  181. 5.16e-184 3.36e-11
10    -5   0.09  5.22   101      9.09    0.332    51  181. 6.61e-207 1.98e- 9
# ℹ 10,191 more rows
# ℹ 1 more variable: dens <dbl>

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

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

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

# 事後分布を作図
ggplot() + 
  geom_contour_filled(
    data    = posterior_df, 
    mapping = aes(x = mu, y = lambda, z = dens, fill = after_stat(level), linetype = "posterior"), 
    alpha = 0.5
  ) + # 事後分布
  geom_vline(
    mapping = aes(xintercept = mu_truth, linetype = "model"), 
    color = "red", linewidth = 1
  ) + # 真のパラメータ
  geom_hline(
    mapping = aes(yintercept = lambda_truth, linetype = "model"), 
    color = "red", linewidth = 1
  ) + # 真のパラメータ
  scale_x_continuous(
    sec.axis = sec_axis(
      transform = ~ ., 
      breaks    = mu_truth, 
      labels    = expression(mu[truth])
    ) # パラメータラベル
  ) + 
  scale_y_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), order = 1), 
    fill     = guide_legend(order = 2)
  ) + 
  labs(
    title = "Gaussian-Gamma distribution", 
    subtitle = posterior_param_lbl, 
    fill = "density", 
    x = expression(mu), 
    y = expression(lambda)
  )

　真のパラメータを赤色の直線(破線)、事後分布(ガウス-ガンマ分布)を塗りつぶしの等高線で示す。
　下の図は、パラメータの真の値から標準偏差1つ分の範囲を拡大したものである。(そのため、作図用の点の数によっては、等高線が粗くなる。)

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

予測分布の計算

　観測データ $\mathbf{X}$ から、未観測データ $x_{*}$ の予測分布(スチューデントのt分布) $p(x_{*} \mid \mathbf{X}, m, \beta, a, b) = \mathrm{St}(x_{*} \mid \hat{\mu}_s, \hat{\lambda}_{s}, \hat{\nu}_s)$ のパラメータ $\hat{\mu}_s, \hat{\lambda}_{s}, \hat{\nu}_s$ を求める。
　スチューデントのt分布については「【R】1次元スチューデントのt分布の作図 - からっぽのしょこ」を参照のこと。

　予測分布のパラメータ $\hat{\mu}_s, \hat{\lambda}_{s}, \hat{\nu}_s$ を計算する。

# 事後分布のパラメータにより予測分布のパラメータを計算:式(3.95')
mu_s_hat     <- m_hat
lambda_s_hat <- beta_hat / (1 + beta_hat) * a_hat / b_hat
nu_s_hat     <- 2 * a_hat
mu_s_hat; lambda_s_hat; nu_s_hat

[1] 5.222297
[1] 0.2789756
[1] 102

# 観測データにより予測分布のパラメータを計算:式(3.95')
mu_s_hat     <- (sum(x_n) + beta * m) / (N + beta)
term_1       <- (N + beta) / (N + 1 + beta)
term_2       <- (0.5 * N + a)
term_3       <- 0.5 * (sum(x_n^2) + beta * m^2 - 1/(N + beta) * (sum(x_n) + beta * m)^2) + b
lambda_s_hat <- term_1 * term_2 / term_3
nu_s_hat     <- N + 2 * a
mu_s_hat; lambda_s_hat; nu_s_hat

[1] 5.222297
[1] 0.2789756
[1] 102

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

$\displaystyle \begin{aligned} \hat{\mu}_s &= \hat{m} \\ &= \frac{ \sum_{n=1}^N x_n + \beta m }{ N + \beta } \\ \hat{\lambda}_s &= \frac{ \hat{\beta} \hat{a} }{ (1 + \hat{\beta}) \hat{b} } \\ &= \frac{ (N + \beta) \left( \frac{N}{2} + a \right) }{ (N + 1 + \beta) \left[ \frac{1}{2} \left\{ \sum_{n=1}^N x_n^2 + \beta m^2 - \frac{1}{N + \beta} \Bigl( \sum_{n=1}^N x_n + \beta m \Bigr)^2 \right\} + b \right] } \\ \hat{\nu}_s &= 2 \hat{a} \\ &= N + 2 a \end{aligned} \tag{3.95'}$

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

# 予測分布の確率密度を計算:式(3.76)
predict_df <- tibble::tibble(
  x        = x_vec, # 確率変数
  mu_s     = mu_s_hat,         # 位置パラメータ
  lambda_s = lambda_s_hat,     # 逆尺度パラメータ
  sigma_s  = 1/sqrt(lambda_s), # 尺度パラメータ
  nu_s     = nu_s_hat,         # 自由度
  dens     = LaplacesDemon::dst(x = x, mu = mu_s, sigma = sigma_s, nu = nu_s) # 確率密度
)
predict_df

# A tibble: 1,001 × 6
       x  mu_s lambda_s sigma_s  nu_s        dens
   <dbl> <dbl>    <dbl>   <dbl> <dbl>       <dbl>
 1 -5     5.22    0.279    1.89   102 0.000000501
 2 -4.98  5.22    0.279    1.89   102 0.000000524
 3 -4.96  5.22    0.279    1.89   102 0.000000548
 4 -4.94  5.22    0.279    1.89   102 0.000000573
 5 -4.92  5.22    0.279    1.89   102 0.000000600
 6 -4.9   5.22    0.279    1.89   102 0.000000627
 7 -4.88  5.22    0.279    1.89   102 0.000000655
 8 -4.86  5.22    0.279    1.89   102 0.000000685
 9 -4.84  5.22    0.279    1.89   102 0.000000716
10 -4.82  5.22    0.279    1.89   102 0.000000749
# ℹ 991 more rows

　 $x$ の値ごとに、スチューデントのt分布に従う確率密度 $\mathrm{St}(x \mid \hat{\mu}_s, \hat{\lambda}_s, \hat{\nu}_s)$ を計算する。
　スチューデントのt分布の確率密度関数 LaplacesDemon パッケージの dnorm() の確率変数の引数 x に $x$ 、中心の引数 mu に $\hat{\mu}_s$ 、尺度の引数 sigma に $\hat{\sigma}_s = \frac{1}{\sqrt{\hat{\lambda}_s}}$ 、自由度の引数 nu に $\hat{\nu}_s$ を指定する。

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

# 予測分布のラベルを作成
predict_param_lbl <- paste0(
  "list(", 
  "N == ", N, ", ", 
  "mu[truth] == ",      round(mu_truth,     digits = 2), ", ", 
  "lambda[truth] == ",  round(lambda_truth, digits = 5), ", ", 
  "hat(mu)[s] == ",     round(mu_s_hat,     digits = 2), ", ", 
  "hat(lambda)[s] == ", round(lambda_s_hat, digits = 5), ", ", 
  "hat(nu)[s] == ",     round(nu_s_hat,     digits = 1), 
  ")"
) |> 
  parse(text = _)

# 予測分布を作図
ggplot() + 
  geom_line(
    data    = model_df, 
    mapping = aes(x = x, y = dens, linetype = "model"), 
    color = "red", linewidth = 1
  ) + # 真の分布
  geom_line(
    data    = predict_df, 
    mapping = aes(x = x, y = dens, linetype = "predict"), 
    color = "purple", linewidth = 1
  ) + # 予測分布
  scale_linetype_manual(
    breaks = c("model", "predict"), 
    values = c("dashed", "solid"), 
    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 = "Student's t distribution", 
    subtitle = predict_param_lbl, 
    x = expression(x), 
    y = "density"
  )

　真の分布(ガウス分布)を赤色の曲線(破線)、予測分布(スチューデントのt分布)を紫色の曲線(実線)で示す。

　データ数 $N$ が十分に大きいと、未観測データ $x_{*}$ の予測分布 $\mathrm{St}(x_{*} \mid \hat{\mu}_s, \hat{\lambda}_s, \hat{\nu}_s)$ の形状が真の分布 $\mathcal{N}(x_n \mid \mu_{\mathrm{truth}}, \lambda_{\mathrm{truth}}^{-1})$ に近付く。

　以上で、平均と精度が未知の1次元ガウスモデルのベイズ推論を実装した。

スポンサードリンク

学習の推移

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

データ数と分布の関係

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

　 $n$ 個のデータから求めた( $n$ 回更新した)事後分布(ガウス-ガンマ分布) $\mathcal{N}(\mu, \lambda \mid \hat{m}, \hat{\beta}, \hat{a}, \hat{b})$ を等高線や曲面、 $n$ 番目の観測データ $x_n$ に対応する位置 $(\mu, \lambda) = (x_n, (\frac{\sigma}{x_n - \mu})^2)$ を桃色の点で示す。

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

$\displaystyle \begin{aligned} m^{(n)} &\leftarrow \frac{1}{\hat{\beta}^{(n)}} \left( x_n + \beta^{(n-1)} m^{(n-1)} \right) \\ \beta^{(n)} &\leftarrow 1 + \beta^{(n-1)} \\ a^{(n)} &\leftarrow \frac{1}{2} + a^{(n-1)} \\ b^{(n)} &\leftarrow \frac{1}{2} \left( x_n^2 + \beta^{(n-1)} (m^{(n-1)})^2 - \beta^{(n)} (m^{(n)})^2 \right) + b^{(n-1)} \end{aligned}$

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

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

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

　 $n$ 個のデータから求めた( $n$ 回更新した)予測分布(スチューデントのt分布) $\mathrm{St}(x_{*} \mid \hat{\mu}_s, \hat{\lambda}_s^{-1}, \hat{\nu}_s)$ を紫色の曲線(実線)、 $n$ 番目の観測データ $x_n$ を桃色の点で示す。

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

$\displaystyle \begin{aligned} \mu_s^{(n)} &= m^{(n-1)} \\ \lambda_s^{(n)} &= \frac{\beta^{(n)} a^{(n)}}{(1 + \beta^{(n)}) b^{(n)}} \\ \nu_s^{(n)} &\leftarrow 1 + \nu_s^{(n-1)} \end{aligned}$

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

　データ数 $N$ が大きくなるのに応じて、未観測データ $x_{*}$ の予測分布 $\mathrm{St}(x_{*} \mid \hat{\mu}_s, \hat{\lambda}_s^{-1}, \hat{\nu}_s)$ の形状が真の分布 $\mathcal{N}(x_n \mid \mu_{\mathrm{truth}}, \lambda_{\mathrm{truth}}^{-1})$ に近付いていくのを確認できる。

観測データと分布の関係

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

　真の分布の期待値 $\mathbb{E}[x_n]$ と真のパラメータ $\mu_{\mathrm{truth}}, \lambda_{\mathrm{truth}}$ の各図(分布)に対応する位置を赤色の垂直線など(破線)、観測データの標本平均 $\bar{x}$ を桃色の垂直線(破線)、事後分布の期待値 $\mathbb{E}[\mu], \mathbb{E}[\lambda]$ と予測分布の期待値 $\mathbb{E}[x_{*}]$ を紫色の垂直線など(破線)で示す。
　また、 $\sigma$ 軸と $\lambda$ 軸、 $\lambda$ 軸と $p(x)$ 軸の対応関係を恒等関数や変換曲線で示す。

　生成分布(ガウス分布)、事後分布(ガウス-ガンマ分布)、予測分布(スチューデントのt分布)の期待値は、それぞれ次の式で計算できる。

$\displaystyle \begin{aligned} \mathbb{E}_{\mathcal{N}(x_n \mid \mu_{\mathrm{truth}}, \lambda_{\mathrm{truth}}^{-1})}[x] &= \mu_{\mathrm{truth}} \\ \mathbb{E}_{\mathrm{NG}(\mu, \lambda \mid \hat{m}, \hat{\beta}, \hat{a}, \hat{b})}[\mu] &= \hat{m} \\ \mathbb{E}_{\mathrm{NG}(\mu, \lambda \mid \hat{m}, \hat{\beta}, \hat{a}, \hat{b})}[\lambda] &= \frac{\hat{a}}{\hat{b}} \\ &= \frac{1 + \hat{\beta}}{\hat{\beta}} \hat{\lambda}_s \\ \mathbb{E}_{\mathrm{St}(x_{*} \mid \hat{\mu}_s, \hat{\lambda}_s^{-1}, \hat{\nu}_s)}[x] &= \hat{\mu}_s \\ &= \hat{m} \end{aligned}$

　真の分布 $\mathcal{N}(x_n \mid \mu_{\mathrm{truth}}, \lambda_{\mathrm{truth}}^{-1})$ や真のパラメータ $\mu_{\mathrm{truth}}, \lambda_{\mathrm{truth}}$ と、観測データ $\mathbf{X}$ 、事後分布 $\mathrm{NG}(\mu, \lambda \mid \hat{m}, \hat{\beta}, \hat{a}, \hat{b})$ 、予測分布 $\mathrm{St}(x_{*} \mid \hat{\mu}_s, \hat{\lambda}_s^{-1}, \hat{\nu}_s)$ 、またそれぞれの統計量の対応関係を確認できる。

　以上で、平均と精度が未知の1次元ガウスモデルのベイズ推論における学習推移を確認した。

　この記事では、1次元ガウス分布に対するベイズ推論を扱った。次の記事では、多次元ガウス分布を扱う。