からっぽのしょこ

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

トップ > ノート > 3.3.1：1次元ガウス分布の学習と予測：平均が未知の場合【緑ベイズ入門のノート】

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

ノート ノート-緑ベイズ入門 R 機械学習

はじめに

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

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

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

【前節の内容】

www.anarchive-beta.com

【他の節一覧】

www.anarchive-beta.com

【この節の内容】

3.3.1 平均が未知の場合

・事後分布の計算

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

　観測モデル$p(\boldsymbol{X} | \mu)$を用いて、観測データ$\boldsymbol{X}$によって学習した$\mu$の事後分布$p(\mu | \boldsymbol{X})$はベイズの定理より

$$ \begin{align} p(\mu | \boldsymbol{X}) &= \frac{ p(\boldsymbol{X} | \mu) p(\mu) }{ p(\boldsymbol{X}) } \\ &= \frac{ \left \{\prod_{n=1}^N p(x_n | \mu) \right\} p(\mu) }{ p(\boldsymbol{X}) } \\ &= \frac{ \left \{\prod_{n=1}^N \mathcal{N}(x_n | \mu, \lambda^{-1}) \right\} \mathcal{N}(\mu | m, \lambda_{\mu}^{-1}) }{ p(\boldsymbol{X}) } \\ &\propto \left \{\prod_{n=1}^N \mathcal{N}(x_n | \mu, \lambda^{-1}) \right\} \mathcal{N}(\mu | m, \lambda_{\mu}^{-1}) \tag{3.49} \end{align} $$

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

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

$$ \begin{align} \ln p(\mu | \boldsymbol{X}) &= \ln \prod_{n=1}^N \mathcal{N}(x_n | \mu, \lambda^{-1}) + \ln \mathcal{N}(\mu | m, \lambda_{\mu}^{-1}) - \ln p(\boldsymbol{X}) \\ &= \sum_{n=1}^N \ln \mathcal{N}(x_n | \mu, \lambda^{-1}) + \ln \mathcal{N}(\mu | m, \lambda_{\mu}^{-1}) + \mathrm{const.} \\ &= \sum_{n=1}^N \ln \left\{ \frac{1}{\sqrt{2 \pi \lambda^{-1}}} \exp\left( - \frac{(x_n - \mu)^2}{2 \lambda^{-1}} \right) \right\} \\ &\qquad + \ln \left\{ \frac{1}{\sqrt{2 \pi \lambda_{\mu}^{-1}}} \exp\left( - \frac{(\mu - m)^2}{2 \lambda_{\mu}^{-1}} \right) \right\} + \mathrm{const.} \\ &= \sum_{n=1}^N \left\{ - \frac{1}{2} \ln 2 \pi \lambda^{-1} - \frac{(x_n - \mu)^2}{2 \lambda^{-1}} \right\} \\ &\qquad - \frac{1}{2} \ln 2 \pi \lambda_{\mu}^{-1} - \frac{(\mu - m)^2}{2 \lambda_{\mu}^{-1}} + \mathrm{const.} \\ &= - \frac{1}{2} \left\{ \sum_{n=1}^N (x_n - \mu)^2 \lambda + (\mu - m)^2 \lambda_{\mu} \right\} + \mathrm{const.} \\ &= - \frac{1}{2} \left\{ \sum_{n=1}^N ( x_n^2 \lambda - 2 x_n \lambda \mu + \lambda \mu^2 ) + \lambda_{\mu} \mu^2 - 2 m \lambda_{\mu} \mu + m^2 \lambda_{\mu} \right\} + \mathrm{const.} \\ &= - \frac{1}{2} \left( - 2 \sum_{n=1}^N x_n \lambda \mu + N \lambda \mu^2 + \lambda_{\mu} \mu^2 - 2 m \lambda_{\mu} \mu \right) + \mathrm{const.} \\ &= - \frac{1}{2} \left\{ (N \lambda + \lambda_{\mu}) \mu^2 - 2 \left( \sum_{n=1}^N x_n \lambda + m \lambda_{\mu} \right) \mu \right\} + \mathrm{const.} \tag{3.50} \end{align} $$

【途中式の途中式】

  1. 式(3.49)の下から2行目の式を用いている。
  2. 対数をとると積が和に変わる。また$\mu$と関係のない項($- \ln p(\boldsymbol{X})$)を$\mathrm{const.}$とおく。
  3. $\mathcal{N}(x_n | \mu, \lambda^{-1}),\ \mathcal{N}(\mu | m, \lambda_{\mu}^{-1})$に、1次元ガウス分布の定義式(2.64)を用いて具体的な確率分布を代入する。
  4. それぞれ対数をとる。$\ln \sqrt{\frac{1}{x}} = \ln x^{-\frac{1}{2}} = -\frac{1}{2} \ln x$、また$\ln \exp(x) = x$である。
  5. $-\frac{1}{2}$を括り出し、$\mu$と関係のない項を$\mathrm{const.}$に含める。
  6. 2乗部分を展開する。
  7. $\sum_{n=1}^N$に関する括弧を展開し、$\mu$と関係のない項を$\mathrm{const.}$に含める。
  8. $\mu^2,\ \mu$の項をそれぞれまとめて式を整理する。


となる。

　これまでと同様に、この式が(対数をとった)1次元ガウス分布の形(2.65)をしていることを確認したい。そこで

$$ \begin{align} \hat{\lambda}_{\mu} &= N \lambda + \lambda_{\mu} \tag{3.53}\\ \tilde{m} &= \lambda \sum_{n=1}^N x_n + \lambda_{\mu} m \end{align} $$

とおき、平方完成を行う。

$$ \begin{aligned} \ln p(\mu | \boldsymbol{X}) &= - \frac{1}{2} \left\{ \hat{\lambda}_{\mu} \mu^2 - 2 \tilde{m} \mu \right\} + \mathrm{const.} \\ &= - \frac{1}{2} \hat{\lambda}_{\mu} \left\{ \mu^2 - 2 \frac{ \tilde{m} }{ \hat{\lambda}_{\mu} } \mu \right\} + \mathrm{const.} \\ &= - \frac{1}{2} \hat{\lambda}_{\mu} \left\{ \mu^2 - 2 \frac{ \tilde{m} }{ \hat{\lambda}_{\mu} } \mu + \left( \frac{\tilde{m}}{\hat{\lambda}_{\mu}} \right)^2 - \left( \frac{\tilde{m}}{\hat{\lambda}_{\mu}} \right)^2 \right\} + \mathrm{const.} \\ &= - \frac{1}{2} \hat{\lambda}_{\mu} \left\{ \left( \mu - \frac{\tilde{m}}{\hat{\lambda}_{\mu}} \right)^2 - \left( \frac{\tilde{m}}{\hat{\lambda}_{\mu}} \right)^2 \right\} + \mathrm{const.} \\ &= - \frac{1}{2} \left( \mu - \frac{\tilde{m}}{\hat{\lambda}_{\mu}} \right)^2 \hat{\lambda}_{\mu} + \frac{\tilde{m}^2}{2 \hat{\lambda}_{\mu}} + \mathrm{const.} \end{aligned} $$

後の項は$\mu$と無関係であるため$\mathrm{const.}$に含めて

$$ \begin{align} \hat{m} &= \frac{\tilde{m}}{\hat{\lambda}_{\mu}} \\ &= \frac{ \lambda \sum_{n=1}^N x_n + \lambda_{\mu} m }{ \hat{\lambda}_{\mu} } \tag{3.54} \end{align} $$

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

$$ \ln p(\mu | \boldsymbol{X}) = - \frac{1}{2} ( \mu - \hat{m} )^2 \hat{\lambda}_{\mu} + \mathrm{const.} = \ln \mathcal{N}(\mu | \hat{m}, \hat{\lambda}_{\mu}^{-1}) \tag{3.51} $$

平均$\hat{\mu}$、精度$\hat{\lambda}_{\mu}^{-1}$の1次元ガウス分布であることが確認できる。
　また事後分布のパラメータの計算式(3.53)、(3.54)によって、ハイパーパラメータ$m,\ \lambda_{\mu}$が更新される。

・予測分布の計算

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

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

$$ \begin{align} p(x_{*}) &= \int p(x_{*} | \mu) p(\mu) d\mu \\ &= \int \mathcal{N}(x_{*} | \mu, \lambda^{-1}) \mathcal{N}(\mu | m, \lambda_{\mu}^{-1}) d\mu \tag{3.55} \end{align} $$

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

$$ p(\mu | x_{*}) = \frac{ p(x_{*} | \mu) p(\mu) }{ p(x_{*}) } \tag{3.56} $$

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

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

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

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

$$ \begin{align} \ln p(\mu | x_{*}) &= - \frac{1}{2} \left( \mu - \frac{ \lambda x_{*} + \lambda_{\mu} m }{ \lambda + \lambda_{\mu} } \right)^2 (\lambda + \lambda_{\mu}) + \mathrm{const.} \\ &= \ln \mathcal{N}( \mu | m(x_{*}), (\lambda + \lambda_{\mu})^{-1} ) \tag{3.58} \end{align} $$

平均$m(x_{*})$、精度$\lambda + \lambda_{\mu}$を持つ1次元ガウス分布となる。ただし

$$ m(x_{*}) = \frac{ \lambda x_{*} + \lambda_{\mu} m }{ \lambda + \lambda_{\mu} } \tag{3.59} $$

とおく。

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

$$ \begin{align} \ln p(x_{*}) &= \ln p(x_{*} | \mu) - \ln p(\mu | x_{*}) + \mathrm{const.} \tag{3.57}\\ &= \ln \mathcal{N}(x_{*} | \mu, \lambda^{-1}) - \ln \mathcal{N}\left( \mu | m(x_{*}), (\lambda + \lambda_{\mu})^{-1} \right) + \mathrm{const.} \\ &= - \frac{1}{2} \Bigl\{ \left( x_{*} - \mu \right)^2 \lambda + \ln 2 \pi \lambda^{-1} \Bigr\} \\ &\qquad + \frac{1}{2} \left\{ \left( \mu - \frac{ \lambda x_{*} + \lambda_{\mu} m }{ \lambda + \lambda_{\mu} } \right)^2 (\lambda + \lambda_{\mu}) + \ln 2 \pi (\lambda + \lambda_{\mu})^{-1} \right\} + \mathrm{const.} \\ &= - \frac{1}{2} \left\{ \left( x_{*} - \mu \right)^2 \lambda - \left( \mu - \frac{ \lambda x_{*} + \lambda_{\mu} m }{ \lambda + \lambda_{\mu} } \right)^2 (\lambda + \lambda_{\mu}) \right\} + \mathrm{const.} \\ &= - \frac{1}{2} \Biggl\{ \lambda x_{*}^2 - 2 \mu \lambda x_{*} + \mu^2 \lambda \Biggr.\\ &\qquad \Biggl. - \mu^2 (\lambda + \lambda_{\mu}) + 2 \mu (\lambda x_{*} + \lambda_{\mu} m) - \frac{ (\lambda x_{*} + \lambda_{\mu} m)^2 }{ \lambda + \lambda_{\mu} } \Biggr\} + \mathrm{const.} \\ &= - \frac{1}{2} \left\{ \lambda x_{*}^2 - 2 \mu \lambda x_{*} + 2 \mu \lambda x_{*} - \frac{ \lambda^2 x_{*}^2 + 2 m \lambda \lambda_{\mu} x_{*} }{ \lambda + \lambda_{\mu} } \right\} + \mathrm{const.} \\ &= - \frac{1}{2} \left\{ \frac{ \lambda (\lambda + \lambda_{\mu}) - \lambda^2 }{ \lambda + \lambda_{\mu} } x_{*}^2 - \frac{ 2 m \lambda \lambda_{\mu} }{ \lambda + \lambda_{\mu} } x_{*} \right\} + \mathrm{const.} \\ &= - \frac{1}{2} \left\{ \frac{ \lambda \lambda_{\mu} }{ \lambda + \lambda_{\mu} } x_{*}^2 - \frac{ 2 m \lambda \lambda_{\mu} }{ \lambda + \lambda_{\mu} } x_{*} \right\} + \mathrm{const.} \tag{3.60} \end{align} $$

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

　この式について

$$ \begin{align} \lambda_{*} &= \frac{ \lambda \lambda_{\mu} }{ \lambda + \lambda_{\mu} } \tag{3.62a}\\ \tilde{\mu}_{*} &= \frac{ 2 m \lambda \lambda_{\mu} }{ \lambda + \lambda_{\mu} } \end{align} $$

とおき、平方完成を行うと

$$ \begin{aligned} \ln p(x_{*}) &= - \frac{1}{2} \{ \lambda_{*} x_{*}^2 - \tilde{\mu}_{*} x_{*} \} + \mathrm{const.} \\ &= - \frac{1}{2} \lambda_{*} \left\{ x_{*}^2 - \frac{ \tilde{\mu}_{*} }{ \lambda_{*} } x_{*} \right\} + \mathrm{const.} \\ &= - \frac{1}{2} \lambda_{*} \left\{ x_{*}^2 - \frac{ \tilde{\mu}_{*} }{ \lambda_{*} } x_{*} + \left( \frac{ \tilde{\mu}_{*} }{ 2 \lambda_{*} } \right)^2 - \left( \frac{ \tilde{\mu}_{*} }{ 2 \lambda_{*} } \right)^2 \right\} + \mathrm{const.} \\ &= - \frac{1}{2} \lambda_{*} \left\{ \left( x_{*}^2 - \frac{ \tilde{\mu}_{*} }{ 2 \lambda_{*} } \right)^2 - \left( \frac{ \tilde{\mu}_{*} }{ 2 \lambda_{*} } \right)^2 \right\} + \mathrm{const.} \\ &= - \frac{1}{2} \left( x_{*}^2 - \frac{ \tilde{\mu}_{*} }{ 2 \lambda_{*} } \right)^2 \lambda_{*} + \frac{\tilde{\mu}_{*}^2}{4 \lambda_{*}} + \mathrm{const.} \end{aligned} $$

となる。後の項は$x_{*}$と無関係であるため$\mathrm{const.}$に含めて

$$ \begin{align} \mu_{*} &= \frac{ \tilde{\mu}_{*} }{ 2 \lambda_{*} } \\ &= \frac{1}{2} \frac{ \lambda + \lambda_{\mu} }{ \lambda \lambda_{\mu} } \frac{ 2 m \lambda \lambda_{\mu} }{ \lambda + \lambda_{\mu} } \\ &= m \tag{3.62b} \end{align} $$

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

$$ \ln p(x_{*}) = - \frac{1}{2} ( x_{*}^2 - \mu_{*} )^2 \lambda_{*} + \mathrm{const.} = \mathcal{N}(x_{*} | \mu_{*}, \lambda_{*}) \tag{3.61} $$

平均$\mu_{*}$、精度$\lambda_{*}$の1次元ガウス分布となることが分かる。

　ちなみに精度パラメータ$\lambda_{*}$の逆数である分散$\sigma^2 = \lambda_{*}^{-1}$に注目すると

$$ \begin{align} \lambda_{*}^{-1} &= \frac{ \lambda + \lambda_{\mu} }{ \lambda \lambda_{\mu} } \\ &= \frac{ \lambda_{\mu} }{ \lambda \lambda_{\mu} } + \frac{ \lambda }{ \lambda \lambda_{\mu} } \\ &= \frac{1}{\lambda} + \frac{1}{\lambda_{\mu}} \\ &= \lambda^{-1} + \lambda_{\mu}^{-1} \tag{3.63} \end{align} $$

である。

　精度パラメータ$\lambda_{*}$を構成するパラメータ$\lambda_{\mu},\ m$について、事後分布パラメータ(3.53)、(3.54)を用いると

$$ \begin{align} \hat{\lambda}_{*} &= \frac{ \lambda \hat{\lambda}_{\mu} }{ \lambda + \hat{\lambda}_{\mu} } \tag{3.62a}\\ &= \frac{ \lambda (N \lambda + \lambda_{\mu}) }{ \lambda + N \lambda + \lambda_{\mu} } \\ &= \frac{ (N \lambda + \lambda_{\mu}) \lambda }{ (N + 1)\lambda + \lambda_{\mu} } \end{align} $$

また

$$ \begin{align} \hat{\mu}_{*} &= \hat{m} \tag{3.62b}\\ &= \frac{ \lambda \sum_{n=1}^N x_n + \lambda_{\mu} m }{ \hat{\lambda}_{\mu} } \\ &= \frac{ \lambda \sum_{n=1}^N x_n + \lambda_{\mu} m }{ N \lambda + \lambda_{\mu} } \end{align} $$

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

$$ p(x_{*} | \boldsymbol{X}) = \mathcal{N}(x_{*} | \hat{\mu}_{*}, \hat{\lambda}_{*}) $$

が得られる。

・Rでやってみよう

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

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

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

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

# 事前分布のパラメータ
m <- 20
lambda_mu <- 0.001

# 試行回数
N <- 50

　観測モデルの平均パラメータ$\mu$をmu_truthとします。この値を推定するのが目的です。
　また観測モデルの精度パラメータ$\lambda$をlambdaとします。この値は推定の中で更新されません。

　続いて事前分布のパラメータ(ハイパーパラメータの初期値)を指定します。平均$m$をm、精度$\lambda_{\mu}$をlambda_muとします。

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

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

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

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

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

# 観測データを確認
summary(x_n)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   3.958  19.825  25.161  25.193  30.541  44.753

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

# 事後分布のパラメータを計算
lambda_mu_hat <- N * lambda + lambda_mu
m_hat <- (lambda * sum(x_n) + lambda_mu * m) / lambda_mu_hat

　事後分布のパラメータの計算式(3.53)、(3.54)の計算を行い、$\hat{\lambda}_{\mu},\ \hat{m}$をそれぞれlambda_mu_hat, m_hatとします。
　m_hatの計算には、事前分布のパラメータlambda_muと事後分布のパラメータlambda_mu_hatの両方を使うので注意しましょう。

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

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

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

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

head(posterior_df)
## # A tibble: 6 x 2
##      mu density
##   <dbl>   <dbl>
## 1  20.9 0.00314
## 2  21.0 0.00320
## 3  21.0 0.00327
## 4  21.0 0.00334
## 5  21.0 0.00341
## 6  21.0 0.00349

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

# 作図
ggplot(posterior_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), 
                         ", lambda_mu_hat=", lambda_mu_hat), 
       x = expression(mu)) # ラベル

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

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

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


・予測分布
# 予測分布のパラメータを計算
lambda_hat <- lambda * lambda_mu_hat / (lambda + lambda_mu_hat)
mu_hat <- m_hat

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

lambda_hat <- (N * lambda + lambda_mu) * lambda / ((N + 1) * lambda + lambda_mu)
mu_hat <- (lambda * sum(x_n) + lambda_mu * m) / (N * lambda + lambda_mu)

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

# 予測分布の計算
predict_df <- tibble(
  x = seq(
    mu_hat - 3 * sqrt(lambda_hat^(-1)), 
    mu_hat + 3 * sqrt(lambda_hat^(-1)), 
    by = 0.01
  ),  # 作図用の値
  density = dnorm(x, mean = mu_hat, sd = sqrt(lambda_hat^(-1))) # 確率密度
)

　$x_{*}$の分布も事後分布のときと同様にして計算します。

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

head(predict_df)
## # A tibble: 6 x 2
##       x  density
##   <dbl>    <dbl>
## 1 -5.12 0.000439
## 2 -5.11 0.000440
## 3 -5.10 0.000441
## 4 -5.09 0.000443
## 5 -5.08 0.000444
## 6 -5.07 0.000445

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

# 作図
ggplot(predict_df, aes(x, 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(mu_hat, 2), ", lambda_mu_hat=", round(lambda_hat, 2)), 
       x = expression(mu)) # ラベル

f:id:anemptyarchive:20200305170156p:plain
$x_{*}$の予測分布：ガウス分布


・おまけ

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

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

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

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

# 事前分布のパラメータ
m <- 20
lambda_mu <- 0.001

# 試行回数
N <- 50

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

# 初期値による予測分布のパラメータを計算
lambda_star <- lambda * lambda_mu / (lambda + lambda_mu)
mu_star <- m

# 初期値による予測分布を計算
predict_df <- tibble(
  x = seq(
    mu_truth - 3 * sqrt(lambda^(-1)), 
    mu_truth + 3 * sqrt(lambda^(-1)), 
    by = 0.01
  ),  # 作図用の値
  density = dnorm(x, mean = mu_star, sd = sqrt(lambda_star^(-1))),  # 確率密度
  N = 0 # 試行回数
)

# パラメータを推定
x_n <- rep(0, N) # 受け皿
for(n in 1:N){
  
  # ガンマ分布に従うデータを生成
  x_n[n] <- rnorm(n = 1, mean = mu_truth, sd = sqrt(lambda^(-1)))
  
  # パラメータを更新
  lambda_mu_old <- lambda_mu
  lambda_mu <- 1 * lambda + lambda_mu
  m <- (lambda * x_n[n] + lambda_mu_old * m) / lambda_mu
  
  # 事後分布を計算
  tmp_posterior_df <- tibble(
    mu = seq(
      mu_truth - 2 * sqrt(lambda^(-1)), 
      mu_truth + 2 * sqrt(lambda^(-1)), 
      by = 0.01
    ),  # 作図用の値
    density = dnorm(mu, mean = m, sd = sqrt(lambda_mu^(-1))),  # 確率密度
    N = n  # 試行回数
  )
  
  # 予測分布のパラメータを計算
  lambda_star <- lambda * lambda_mu / (lambda + lambda_mu)
  mu_star <- m
  
  # 予測分布を計算
  tmp_predict_df <- tibble(
    x = seq(
      mu_truth - 3 * sqrt(lambda^(-1)), 
      mu_truth + 3 * sqrt(lambda^(-1)), 
      by = 0.01
    ),  # 作図用の値
    density = dnorm(x, mean = mu_star, sd = sqrt(lambda_star^(-1))),  # 確率密度
    N = n # 試行回数
  )
  
  # 推定結果を結合
  posterior_df <- rbind(posterior_df, tmp_posterior_df)
  predict_df <- rbind(predict_df, tmp_predict_df)
}

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

## 事後分布
# 作図
posterior_graph <- ggplot(posterior_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_graph)

## 予測分布
predict_graph <- ggplot(predict_df, aes(x, 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(x)) # ラベル

# 描画
animate(predict_graph)

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

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

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

　$\hat{m}$の計算には、更新前の$\lambda_{\mu}$と更新後の$\hat{\lambda}_{\mu}$の両方が必要です。そこで更新前のlambda_muの値をlambda_mu_oldとして残しておきます。そしてlambda_muを更新してから、lambda_mu_oldとlambda_muを使ってmを更新します。

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

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

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


参考文献

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

おわりに

　ここまでは1日1項ペースで進んでいましたが、次とその次でかなり詰まっております…

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

【次節の内容】

www.anarchive-beta.com