からっぽのしょこ

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

3.3.3:1次元ガウス分布の学習と予測:平均・精度が未知の場合【緑ベイズ入門のノート】

はじめに

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

 この記事は3.3.3項の内容です。尤度関数・事前分布を1次元ガウス分布とした場合の平均と精度パラメータの事後分布を導出し、学習した事後分布を用いた予測分布を導出します。またその推論過程をR言語で実装します。

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

【前節の内容】

www.anarchive-beta.com

【他の節一覧】

www.anarchive-beta.com

【この節の内容】

3.3.3 平均・精度が未知の場合

 ガウス分布に従うと仮定するN個の観測データ$\boldsymbol{X} = \{x_1, x_2, \cdots, x_N\}$を用いて、平均パラメータ$\mu$と精度パラメータ$\lambda$の事後分布を求めていく。

 観測データ$\boldsymbol{X}$、平均パラメータ$\mu$、精度パラメータ$\lambda$の同時分布は次のように分解できる。

$$ p(\boldsymbol{X}, \mu, \lambda) = p(\mu | \lambda, \boldsymbol{X}) p(\lambda | \boldsymbol{X}) p(\boldsymbol{X}) \tag{3.84} $$

 また事前分布(パラメータの同時分布)$p(\mu, \lambda)$をガウス・ガンマ分布と仮定すると、同じ形式の事後分布となることが知られている。

$$ \begin{align} p(\mu, \lambda) &= \mathrm{NB}(\mu, \lambda | m, \beta, a, b) \\ &= \mathcal{N}(\mu | m, (\beta \lambda)^{-1}) \mathrm{Gam}(\lambda | a, b) \\ &= \frac{1}{\sqrt{2 \pi (\beta \lambda)^{-1}}} \exp\left\{ - \frac{ (\mu - m)^2 }{ 2 (\beta \lambda)^{-1} } \right\} C_G(a, b) \lambda^{a-1} \exp(- b \lambda) \tag{3.81} \end{align} $$

 対数をとると

$$ \begin{aligned} \ln p(\mu, \lambda) &= \ln \mathrm{NB}(\mu, \lambda | m, \beta, a, b) \\ &= \ln \mathcal{N}(\mu | m, (\beta \lambda)^{-1}) + \ln \mathrm{Gam}(\lambda | a, b) \\ &= - \frac{1}{2} \left\{ (\mu - m)^2 \beta \lambda + \ln (\beta \lambda)^{-1} + \ln 2 \pi \right\} \\ &\qquad + (a - 1) \ln \lambda - b \lambda + \ln C_G(a, b) \end{aligned} $$

となる。

・平均$\mu$の事後分布の計算

 まずは$\mu$の事後分布$p(\mu | \lambda, \boldsymbol{X})$を導出する。

 図3.5のグラフィカルモデルでパラメータの依存関係を確認すると、3.3.1項と同様の手順で求められることが分かる。よって$\mu$の事後分布のパラメータ(3.53)、(3.54)について、$\beta \lambda = \lambda_{\mu}$で置き換えると

$$ \begin{align} \hat{\lambda}_{\mu} &= N \lambda + \lambda_{\mu} \tag{3.53}\\ \hat{\beta} \lambda &= N \lambda + \beta \lambda \\ \hat{\beta} &= N + \beta \tag{3.83a} \end{align} $$

また

$$ \begin{align} \hat{m} &= \frac{ \lambda \sum_{n=1}^N x_n + \lambda_{\mu} m }{ \hat{\lambda}_{\mu} } \tag{3.54}\\ &= \frac{ \lambda \sum_{n=1}^N x_n + \beta \lambda m }{ N \lambda + \beta \lambda } \\ &= \frac{1}{\hat{\beta}} \left( \sum_{n=1}^N x_n + \beta m \right) \tag{3.83b} \end{align} $$

となる。従って$\mu$の事後分布$p(\mu | \lambda, \boldsymbol{X})$は

$$ p(\mu | \lambda, \boldsymbol{X}) = \mathcal{N}( \mu | \hat{m}, (\hat{\beta} \lambda)^{-1} ) \tag{3.82} $$

平均パラメータ$\mu$、精度パラメータ$(\hat{\beta} \lambda)^{-1}$を持つ1次元ガウス分布となる。

・精度$\lambda$の事後分布の計算

 次に$\lambda$の事後分布$p(\lambda | \boldsymbol{X})$を導出する。

 $\lambda$の事後分布$p(\lambda | \boldsymbol{X})$は式(3.84)より

$$ \begin{align} p(\lambda | \boldsymbol{X}) &= \frac{ p(\boldsymbol{X}, \mu, \lambda) }{ p(\mu | \lambda, \boldsymbol{X}) p(\boldsymbol{X}) } \\ &= \frac{ p(\boldsymbol{X} | \mu, \lambda) p(\mu, \lambda) }{ p(\mu | \lambda, \boldsymbol{X}) p(\boldsymbol{X}) } \\ &= \frac{ \prod_{n=1}^N p(x_n | \mu, \lambda) p(\mu, \lambda) }{ p(\mu | \lambda, \boldsymbol{X}) p(\boldsymbol{X}) } \\ &\propto \frac{ \prod_{n=1}^N p(x_n | \mu, \lambda) p(\mu, \lambda) }{ p(\mu | \lambda, \boldsymbol{X}) } \tag{3.85} \end{align} $$

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

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

$$ \begin{align} \ln p(\lambda | \boldsymbol{X}) &= \sum_{n=1}^N \ln p(x_n | \mu, \lambda) + \ln p(\mu, \lambda) - \ln p(\mu | \lambda, \boldsymbol{X}) + \mathrm{const.} \\ &= \sum_{n=1}^N \ln \mathcal{N}(x_n | \mu, \lambda^{-1}) + \ln \mathcal{N}(\mu | m, (\beta \lambda)^{-1}) + \ln \mathrm{Gam}(\lambda | a, b) \\ &\qquad - \ln \mathcal{N}(\mu | \hat{m}, (\hat{\beta} \lambda)^{-1}) + \mathrm{const.} \\ &= - \frac{1}{2} \sum_{n=1}^N \left\{ (x_n - \mu)^2 \lambda + \ln \lambda^{-1} + \ln 2 \pi \right\} \\ &\qquad - \frac{1}{2} \left\{ (\mu - m)^2 \beta \lambda + \ln (\beta \lambda)^{-1} + \ln 2 \pi \right\} \\ &\qquad + (a - 1) \ln \lambda - b \lambda + \ln C_G(a, b) \\ &\qquad + \frac{1}{2} \left\{ (\mu - \hat{m})^2 \hat{\beta} \lambda + \ln (\hat{\beta} \lambda)^{-1} + \ln 2 \pi \right\} + \mathrm{const.} \\ &= - \frac{1}{2} \sum_{n=1}^N x_n^2 \lambda + \sum_{n=1}^N x_n \mu \lambda - \frac{N}{2} \mu^2 \lambda \\ &\qquad - \frac{1}{2} \mu^2 \beta \lambda + \mu m \beta \lambda - \frac{1}{2} m^2 \beta \lambda \\ &\qquad + (a - 1) \ln \lambda - b \lambda + a \ln b - \ln \Gamma(a) \\ &\qquad + \frac{1}{2} \mu^2 (N + \beta) \lambda - \mu \frac{1}{N + \beta} \left( \sum_{n=1}^N x_n + \beta m \right) (N + \beta) \lambda + \frac{1}{2} \hat{m}^2 \hat{\beta} \lambda + \mathrm{const.} \\ &= - \frac{1}{2} \sum_{n=1}^N x_n^2 \lambda - \frac{1}{2} m^2\beta \lambda + (a - 1) \ln \lambda - b \lambda + \frac{1}{2} \hat{m}^2 \hat{\beta} \lambda + \mathrm{const.} \\ &= \left( \frac{N}{2} + a - 1 \right) \ln \lambda - \left\{ \frac{1}{2} \left( \sum_{n=1}^N x_n^2 + \beta m^2 - \hat{\beta} \hat{m}^2 \right) + b \right\} \lambda + \mathrm{const.} \tag{3.86} \end{align} $$

となる。適宜$\lambda$と関係のない項を$\mathrm{const.}$に含めている。

 この式について

$$ \begin{align} \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 \tag{3.88} \end{align} $$

とおくと、$\lambda$の事後分布$p(\lambda | \boldsymbol{X})$は

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

パラメータ$\hat{a},\ \hat{b}$を持つガンマ分布となる。
 また事後分布のパラメータの計算式(3.88)によって、ハイパーパラメータが更新される。

・予測分布の計算

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

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

$$ p(x_{*}) = \iint p(x_{*} | \mu, \lambda) p(\mu, \lambda) d\mu d\lambda \tag{3.89} $$

 ここでもこの計算を直接するのではなく、予測分布$p(x_{*})$と事前分布$p(\mu, \lambda)$の関係を考えてみる。ベイズの定理より

$$ p(\mu, \lambda | x_{*}) = \frac{ p(x_{*} | \mu, \lambda) p(\mu, \lambda) }{ p(x_{*}) } $$

という関係が成り立つことが分かる。両辺で対数をとり、予測分布に関して整理すると

$$ \begin{align} \ln p(\mu, \lambda | x_{*}) &= \ln p(x_{*} | \mu, \lambda) + \ln p(\mu, \lambda) - \ln p(x_{*}) \\ \ln p(x_{*}) &= \ln p(x_{*} | \mu, \lambda) - p(\mu, \lambda | x_{*}) + \mathrm{const.} \tag{3.90} \end{align} $$

となる。$x_{*}$と関係のない$\ln p(\mu, \lambda)$を$\mathrm{const.}$とおいている。
 つまり式(3.90)の計算によっても予測分布を求められることが分かる。

 $p(\mu, \lambda | x_{*})$は、1つのデータ$x_{*}$が与えられた下での$\mu,\ \lambda$の条件付き同時分布である。つまり$p(\mu, \lambda | x_{*})$は、N個の観測データが与えられた下での事後分布$p(\mu | \boldsymbol{X}),\ p(\lambda | \boldsymbol{X})$の導出と同様の手順で求められる。
 従って、事後分布のパラメータ(3.82)、(3.87)を用いると、$p(\mu, \lambda | x_{*})$は$N = 1$より

$$ \begin{align} p(\mu, \lambda | x_{*}) &= p(\mu | \lambda, x_{*}) p(\lambda | x_{*}) \\ &= \mathcal{N}(\mu | m(x_{*}), \{(1 + \beta) \lambda\}^{-1}) \mathrm{Gam} \left( \lambda \middle| \frac{1}{2} + a, b(x_{*}) \right) \tag{3.91} \end{align} $$

となる。ただし

$$ \begin{align} m(x_{*}) &= \frac{ x_{*} + \beta m }{ 1 + \beta } \\ b(x_{*}) &= \frac{ \beta }{ 2 (1 + \beta) } (x_{*} - m)^2 + b \tag{3.92} \end{align} $$

とおく。

 $p(\mu, \lambda | x_{*})$が求まったので、予測分布$p(x_{*})$の計算に戻る。
 式(3.92)を式(3.90)に代入して$x_{*}$について整理すると

$$ \begin{align} \ln p(x_{*}) &= \ln p(x_{*} | \mu, \lambda) - p(\mu, \lambda | x_{*}) + \mathrm{const.} \tag{3.90}\\ &= \ln \mathcal{N}(x_{*} | \mu, \lambda^{-1}) \\ &\qquad - \ln \mathcal{N}(\mu | m(x_{*}), \{(1 + \beta) \lambda\}^{-1}) - \ln \mathrm{Gam} \left( \lambda \middle| \frac{1}{2} + a, b(x_{*}) \right) + \mathrm{const.} \\ &= - \frac{1}{2} \Bigl\{ (x_{*} - \mu)^2 \lambda + \ln \lambda^{-1} + \ln 2 \pi \Bigr\} \\ &\qquad + \frac{1}{2} \left[ \left( \mu - \frac{ x_{*} + \beta m }{ 1 + \beta } \right)^2 (1 + \beta) \lambda + \ln \{(1 + \beta) \lambda\}^{-1} + \ln 2 \pi \right] \\ &\qquad - \left(\frac{1}{2} + a - 1\right) \ln \lambda + \frac{ \beta }{ 2 (1 + \beta) } \{(x_{*} - m)^2 + b\} \lambda \\ &\qquad - \left(\frac{1}{2} + a\right) \ln \left\{ \frac{ \beta }{ 2 (1 + \beta) } (x_{*} - m)^2 + b \right\} + \ln \Gamma \left(\frac{1}{2} + a\right) + \mathrm{const.} \\ &= - \frac{1}{2} \lambda x_{*}^2 + \mu \lambda x_{*} - \frac{1}{2} \mu^2 \lambda \\ &\qquad + \frac{1}{2} (1 + \beta) \mu^2 \lambda - (x_{*} + \beta m) \mu \lambda + \frac{ (x_{*} + \beta m)^2 \lambda }{ 2 (1 + \beta) } \\ &\qquad + \frac{ \beta (x_{*} - m)^2 \lambda }{ 2 (1 + \beta) } + \frac{ b \beta \lambda }{ 2 (1 + \beta) } \\ &\qquad - \left(\frac{1}{2} + a\right) \ln \left\{ \frac{ \beta }{ 2 (1 + \beta) b } (x_{*} - m)^2 + 1 \right\} b + \mathrm{const.} \\ &= - \frac{ (1 + \beta) \lambda x_{*}^2 }{ 2 (1 + \beta) } + \mu \lambda x_{*} \\ &\qquad - \mu \lambda x_{*} - \beta m \mu \lambda + \frac{ \lambda x_{*}^2 + 2 \beta m \lambda x_{*} + \beta^2 m^2 \lambda }{ 2 (1 + \beta) } \\ &\qquad + \frac{ \beta \lambda x_{*}^2 - 2 \beta m \lambda x_{*} + \beta m^2 \lambda }{ 2 (1 + \beta) } \\ &\qquad - \left(\frac{1}{2} + a\right) \ln \left\{ \frac{ \beta }{ 2 (1 + \beta) b } (x_{*} - m)^2 + 1 \right\} - \left(\frac{1}{2} + a\right) \ln b + \mathrm{const.} \\ &= - \frac{ \lambda x_{*}^2 + \beta \lambda x_{*}^2 }{ 2 (1 + \beta) } + \frac{ \lambda x_{*}^2 }{ 2 (1 + \beta) } + \frac{ \beta \lambda x_{*}^2 }{ 2 (1 + \beta) } \\ &\qquad - \left(\frac{1}{2} + a\right) \ln \left\{ \frac{ \beta }{ 2 (1 + \beta) b } (x_{*} - m)^2 + 1 \right\} + \mathrm{const.} \\ &= - \frac{1 + 2a}{2} \ln \left\{ 1 + \frac{ \beta }{ 2 (1 + \beta) b } (x_{*} - m)^2 \right\} + \mathrm{const.} \end{align} $$

となる。適宜$x_{*}$と関係のない項を$\mathrm{const.}$にまとめている。

 この式について、$\frac{\beta}{2 (1 + \beta) b} = \frac{\beta a}{2 a (1 + \beta) b}$として

$$ \begin{align} \mu_s &= m \\ \lambda_s &= \frac{ \beta a }{ (1 + \beta) b } \\ \nu_s &= 2 a \tag{3.95} \end{align} $$

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

$$ \ln p(x_{*}) = - \frac{1 + \nu_s}{2} \ln \left\{ 1 + \frac{ \lambda_s }{ \nu_s } (x_{*} - \mu_s)^2 \right\} + \mathrm{const.} = \ln \mathrm{St}(x_{*} | \mu_s, \lambda_s, \nu_s) \tag{3.94} $$

パラメータ$\mu_s,\ \lambda_s,\ \nu_s$を持つスチューデントのt分布となることが分かる。

 予測分布のパラメータ$\mu_s,\ \lambda_s,\ \nu_s$を構成する事前分布のパラメータ$\hat{\beta},\ \hat{m},\ a,\ b$について、事後分布のパラメータ(3.83)、(3.88)を用いると

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

また

$$ \begin{aligned} \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 - \hat{\beta} \hat{m}^2 \right) + b \right\} } \end{aligned} $$

また

$$ \begin{aligned} \hat{\nu}_s &= 2 \hat{a} \\ &= N + 2 a \end{aligned} $$

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

$$ \begin{align} \ln p(x_{*} | \boldsymbol{X}) &= - \frac{1 + \hat{\nu}_s}{2} \ln \left\{ 1 + \frac{ \hat{\lambda}_s }{ \hat{\nu}_s } (x_{*} - \hat{\mu}_s)^2 \right\} + \mathrm{const.} \\ &= \ln \mathrm{St}(x_{*} | \hat{\mu}_s, \hat{\lambda}_s, \hat{\nu}_s) \tag{3.94} \end{align} $$

が得られる。

・Rでやってみよう

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

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

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

・事後分布
## パラメータの初期値を設定
# 観測モデルのパラメータ
mu_truth <- 25
lambda_truth <- 0.01

# 平均muの事前分布のパラメータ
m <- 20
beta <- 1

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

# 試行回数
N <- 50

 観測モデルの精度パラメータ$\lambda$をlambda_truth、平均パラメータ$\mu$をmuとします。この値を推定するのが目的です。

 続いて事前分布のパラメータ(ハイパーパラメータの初期値)を指定します。$m$をm、$\beta$をbeta、$a$をa、$b$をbとします。

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

 $\lambda$は分散の逆数であるため、値が大きくなるほど散らばり具合が小さくなります。

# ガウス分布に従うデータを生成
x_n <- rnorm(n = N, mean = mu_truth, sd = sqrt(lambda_truth^(-1)))

 ガウス分布に従う乱数を発生させる関数rnorm()を使って、ランダムにデータを生成します。
 試行回数の引数nにはN、平均の引数meanにはmu_truthを指定します。
 しかし標準偏差の引数sdにはlambda_truthのままでは指定できません。そこで標準偏差$\sigma = \sqrt{\lambda^{-1}}$に置き換えるため、sqrt(lambda_truth^(-1))とします。

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

# 観測データを確認
summary(x_n)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.3328 17.2163 25.7692 25.4834 33.3597 42.5555

このデータを用いてハイパーパラメータを更新します。

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

 $\mu$の事後分布のパラメータ(3.83)の計算を行い、$\hat{\beta},\ \hat{m}$をそれぞれbeta_hat, m_hatとします。
 $\lambda$の事後分布のパラメータ(3.88)の計算を行い、$\hat{a},\ \hat{b}$をそれぞれa_hat, b_hatとします。
 また$\lambda$については分布推定しているため、ガンマ分布の期待値$\mathbb{E}[\lambda] = \frac{\hat{a}}{\hat{b}}$を用いることにします。計算結果はlambda_barとします。

# muの事後分布を計算
posterior_mu_df <- tibble(
  mu = seq(
    mu_truth - 3 * sqrt((beta_hat * lambda_bar)^(-1)), 
    mu_truth + 3 * sqrt((beta_hat * lambda_bar)^(-1)), 
    by = 0.01
  ),  # 作図用の値
  density = dnorm(mu, mean = m_hat, sd = sqrt((beta_hat * lambda_bar)^(-1))) # 確率密度
)

 $\mu$が取り得る値をseq()を使って値を用意します。
 ただしmu_truthの値によって作図に必要な範囲が変わります。ここでは$\mu$を中心に標準偏差$\sigma = \sqrt{(\beta \lambda)^{-1}}$の3倍をプロットすることにします。そこでseq()の第1引数(最小値)、第2引数(最大値)に、mu_truthから3 * sqrt((beta_hat * lambda_bar)^(-1))を引いた値と足した値をそれぞれ指定します。

 muの各値の確率密度は、ガウス分布の定義式(2.64)で計算できます。しかし計算式が複雑なので、ここでは手を抜いてみましょう。
 ガウス分布に従い確率密度を返す関数dnorm()を使います。引数についてはrnorm()と同様です。ただしここでは事前分布のパラメータm_hatbeta_hat * lambda_barを用います。
 次のガンマ分布についてもこのような関数が用意されていますが、そちらは計算式に従って計算しましょう。

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

head(posterior_mu_df)
## # A tibble: 6 x 2
##      mu density
##   <dbl>   <dbl>
## 1  20.8 0.00137
## 2  20.8 0.00140
## 3  20.8 0.00143
## 4  20.8 0.00146
## 5  20.9 0.00150
## 6  20.9 0.00153

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

# 作図
ggplot(posterior_mu_df, aes(mu, density)) + 
  geom_line(color = "#56256E") + # 折れ線グラフ
  geom_vline(aes(xintercept = mu_truth), 
             color = "red", linetype = "dashed") + # 垂直線
  labs(title = "Gaussian Distribution", 
       subtitle = paste0("N=", N, ", m_hat=", round(m_hat, 2), ", beta_hat=", beta_hat), 
       x = expression(mu)) # ラベル

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

 推定値は折れ線グラフgeom_line()を用いて描きます。
 パラメータ$\mu$の真の値を垂直線geom_vline()で示します。
 事前分布のパラメータm_hatr m_hatのように細かい値となるので、round()で丸めておきます。

f:id:anemptyarchive:20200305192009p:plain
$\mu$の事後分布:ガウス分布


# lambdaの事後分布を計算
posterior_lambda_df <- tibble(
  lambda = seq(0, 3 * lambda_truth, by = 0.00001),  # 作図用の値
  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の3倍までとします。そこでseq()の第1引数(最小値)に0、第2引数(最大値)に3 * mu_truthを指定します。

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

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

head(posterior_lambda_df)
## # A tibble: 6 x 3
##    lambda   C_G  density
##     <dbl> <dbl>    <dbl>
## 1 0        146. 0.      
## 2 0.00001  146. 3.13e-62
## 3 0.00002  146. 1.02e-54
## 4 0.00003  146. 2.52e-50
## 5 0.00004  146. 3.26e-47
## 6 0.00005  146. 8.42e-45

 これも同様にggplot2パッケージを利用してグラフを描きます。

# 作図
ggplot(posterior_lambda_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=", round(b_hat, 2)), 
       x = expression(lambda)) # ラベル

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

・予測分布
# 予測分布のパラメータを計算
mu_s_hat <- m_hat
lambda_s_hat <- beta_hat * a_hat / (1 + beta_hat) / b_hat
nu_s_hat <- 2 * a_hat

 予測分布のパラメータの計算式の計算を行い、$\hat{\mu}_s,\ \hat{\lambda}_s,\ \hat{\mu}_s$をそれぞれmu_s_hat, lambda_s_hat, nu_s_hatとします。

mu_s_hat <- (sum(x_n) + beta * m) / (N + beta)
tmp_lambda_numer <- (N + beta) * (N / 2 + a)
tmp_lambda_denom <- (N + 1 + beta) * ((sum(x_n^2) + beta * m^2 - beta_hat * m_hat^2) / 2 + b)
lambda_s_hat <- tmp_lambda_numer / tmp_lambda_denom
nu_s_hat <- N + 2 * a

 このように、観測データx_nと事前分布のパラメータを使って計算することもできます(本当はbeta_hat, m_hatbeta, mに置き換える必要があります。)。

# 予測分布を計算
predict_df <- tibble(
  x = seq(-mu_s_hat, 3 * mu_s_hat, by = 0.001),  # 作図用の値
  C_St = lgamma((nu_s_hat + 1) / 2) - lgamma(nu_s_hat / 2),  # 正規化項(対数)
  term1 = log(lambda_s_hat / pi / nu_s_hat) / 2, 
  term2 = - (nu_s_hat + 1) / 2 * log(1 + lambda_s_hat / nu_s_hat * (x - mu_s_hat)^2), 
  density = exp(C_St + term1 + term2) # 確率密度
)

 $x_{*}$が取り得る値をseq()を使って用意します。
 ここでは、$-\mu_s$から$3 \mu_s$までとしました。必要に応じて調整してください。特にこのままですと、mu_sが負の値をとるとエラーとなります。

 xの各値に対してスチューデントのt分布の定義式(3.76)の計算を行い、確率密度を求めます。
 ちなみにこのpiは、前節で出てきたパラメータではなく円周率のことです。ご注意ください。

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

head(predict_df)
## # A tibble: 6 x 5
##       x  C_St term1 term2    density
##   <dbl> <dbl> <dbl> <dbl>      <dbl>
## 1 -25.4  1.62 -4.86 -10.6 0.00000102
## 2 -25.4  1.62 -4.86 -10.6 0.00000102
## 3 -25.4  1.62 -4.86 -10.6 0.00000102
## 4 -25.4  1.62 -4.86 -10.6 0.00000102
## 5 -25.4  1.62 -4.86 -10.6 0.00000102
## 6 -25.4  1.62 -4.86 -10.6 0.00000102

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

# 作図
ggplot(predict_df, aes(x, density)) + 
  geom_line(color = "#56256E") + # 折れ線グラフ
  labs(title = "Student's t Distribution", 
       subtitle = paste0("N=", N, ", mu=", round(mu_s_hat, 2), 
                         ", lambda=", round(lambda_s_hat, 2), ", nu=", nu_s_hat)) # ラベル

f:id:anemptyarchive:20200305192121p:plain
$x_{*}$の予測分布:スチューデントのt分布

・おまけ

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

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

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

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

# 平均muの事前分布のパラメータ
m <- 20
beta <- 1

# 精度lambdaの事前分布のパラメータ
a <- 1
b <- 1
lambda_bar <- a / b

# 試行回数
N <- 50

# muの事前分布を計算
posterior_mu_df <- tibble(
  mu = seq(
    mu_truth - 2 * sqrt(lambda_truth^(-1)), 
    mu_truth + 2 * sqrt(lambda_truth^(-1)), 
    by = 0.01
  ),  # 作図用の値
  density = dnorm(mu, mean = m, sd = sqrt(beta * lambda_bar)^(-1)), # 確率密度
  N = 0 # 試行回数
)

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

# 初期値による予測分布のパラメータを計算
mu_s <- m
lambda_s <- beta * a / (1 + beta) / b
nu_s <- 2 * a

# 初期値による予測分布を計算
predict_df <- tibble(
  x = seq(-mu_truth, 3 * mu_truth, by = 0.001),  # 作図用の値
  C_St = lgamma((nu_s + 1) / 2) - lgamma(nu_s / 2),  # 正規化項(対数)
  term1 = log(lambda_s / pi / nu_s) / 2, 
  term2 = - (nu_s + 1) / 2 * log(1 + lambda_s / nu_s * (x - mu_s)^2), 
  density = exp(C_St + term1 + term2), # 確率密度
  N = 0 # 試行回数
)

# パラメータを推定
x_n <- rep(0, N)
for(n in 1:N){
  
  # ガウス分布に従うデータを生成
  x_n[n] <- rnorm(n = 1, mean = mu_truth, sd = sqrt(lambda_truth^(-1)))
  
  # パラメータを更新
  beta_old <- beta
  beta <- 1 + beta
  m_old <- m
  m <- (x_n[n] + beta_old * m) / beta
  a <- 1 / 2 + a
  b <- (x_n[n]^2 + beta_old * m_old^2 - beta * m^2) / 2 + b
  lambda_bar <- a / b
  
  # muの事後分布を計算
  tmp_posterior_mu_df <- tibble(
    mu = seq(
      mu_truth - 2 * sqrt(lambda_truth^(-1)), 
      mu_truth + 2 * sqrt(lambda_truth^(-1)), 
      by = 0.01
    ),  # 作図用のmu
    density = dnorm(mu, mean = m, sd = sqrt(beta * lambda_bar)^(-1)), # 確率密度
    N = n # 試行回数
  )
  
  # lambdaの事後分布を計算
  tmp_posterior_lambda_df <- tibble(
    lambda = seq(0, 4 * lambda_truth, by = 0.00001),  # 作図用の値
    C_G = a * log(b) - lgamma(a),  # 正規化項(対数)
    density = exp(C_G + (a - 1) * log(lambda) - b * lambda), # 確率密度
    N = n # 試行回数
  )
  
  # 予測分布のパラメータを更新
  mu_s <- m
  lambda_s <- beta * a / (1 + beta) / b
  nu_s <- 2 * a
  
  # 予測分布を計算
  tmp_predict_df <- tibble(
    x = seq(-mu_truth, 3 * mu_truth, by = 0.001),  # 作図用の値
    C_St = lgamma((nu_s + 1) / 2) - lgamma(nu_s / 2),  # 正規化項(対数)
    term1 = log(lambda_s / pi / nu_s) / 2, 
    term2 = - (nu_s + 1) / 2 * log(1 + lambda_s / nu_s * (x - mu_s)^2), 
    density = exp(C_St + term1 + term2), # 確率密度
    N = n # 試行回数
  )
  
  # 推定結果を結合
  posterior_mu_df <- rbind(posterior_mu_df, tmp_posterior_mu_df)
  posterior_lambda_df <- rbind(posterior_lambda_df, tmp_posterior_lambda_df)
  predict_df <- rbind(predict_df, tmp_predict_df)
}

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

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

# 描画
animate(posterior_mu_graph)

## lambdaの事後分布
# 作図
posterior_lambda_graph <- ggplot(posterior_lambda_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_lambda_graph)

## 予測分布
# 作図
predict_graph <- ggplot(predict_df, aes(x, density)) + 
  geom_line(color = "#56256E") + # 折れ線グラフ
  transition_manual(N) + # フレーム
  labs(title = "Student's t Distribution", 
       subtitle = "N= {current_frame}") # ラベル

# 描画
animate(predict_graph)

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

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

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

 作図範囲を固定するため(更新されない)観測モデルのパラメータを使います。

f:id:anemptyarchive:20200305192420g:plain
$\mu$の事後分布の推移

f:id:anemptyarchive:20200305192447g:plain
$\lambda$の事後分布の推移

f:id:anemptyarchive:20200305192521g:plain
$x_{*}$の予測分布の推移


参考文献

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

おわりに

 ここまではめちゃくちゃ難しい訳ではないのですがとにかく大変でした。
 3月はトピックモデルの実装の予定でして、あと次の節は行列計算が出てくるっぽいので暫くお休みします。あ、でも修正作業と予測分布を組んでみたを追加するのは3月中にやるつもりです。

 4月からはPythonを触ってみる予定ですので、その時に勉強がてらPython版を追加するとともに再開できればと思います。

2020/03/05:加筆修正しました。

 ここまでの感想:データ数が増えると分布の見分けがつかない。

 各節で同じことをしていると分かりやすくするためにも表記・表現の統一に拘っていたのですが、思ったよりも大変で時間をとられてしまいました。なので以降はもう少し緩くまとめていくことにします。全体を読み終えてから、俯瞰的な知識の整理も含めてそういった作業を改めてできればと思います。

【次節の内容】

www.anarchive-beta.com