からっぽのしょこ

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

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

はじめに

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

 この記事では、「尤度関数を平均が未知の1次元ガウス分布(一変量正規分布)」、「事前分布を1次元ガウス分布」とした場合の「平均パラメータの事後分布」と「未観測値の予測分布」を導出します。

【実装編】

www.anarchive-beta.com

www.anarchive-beta.com

【前の節の内容】

www.anarchive-beta.com

【他の節一覧】

www.anarchive-beta.com

【この節の内容】

3.3.1 平均が未知の場合

 尤度関数を1次元ガウス分布(Gaussian Distribution)・一変量正規分布(Normal Distribution)、事前分布を1次元ガウス分布とするモデルに対するベイズ推論を導出する。
 ガウス分布については「1次元ガウス分布の定義式 - からっぽのしょこ」を参照のこと。

尤度関数の確認

 ガウス分布に従うと仮定する$N$個の観測データ$\mathbf{X} = \{x_1, x_2, \cdots, x_N\}$を用いて、平均パラメータ$\mu$の事後分布と未観測データ$x_{*}$の予測分布を求めていく。
 まずは、尤度関数を確認する。

 $\mathbf{X}$の各データが独立に発生しているとの仮定の下で、尤度$p(\mathbf{X} | \mu)$は

$$ \begin{aligned} p(\mathbf{X} | \mu) &= p(x_1, x_2, \cdots, x_N | \mu) \\ &= p(x_1 | \mu) p(x_2 | \mu) \cdots p(x_N | \mu) \\ &= \prod_{n=1}^N p(x_n | \mu) \end{aligned} $$

と分解できる。

事後分布の計算

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

 $\mathbf{X}$が与えられた下での$\mu$の事後分布は、ベイズの定理より

$$ \begin{align} p(\mu | \mathbf{X}) &= \frac{ p(\mathbf{X} | \mu) p(\mu) }{ p(\mathbf{X}) } \\ &= \frac{ \left\{ \prod_{n=1}^N p(x_n | \mu) \right\} p(\mu) }{ p(\mathbf{X}) } \\ &= \frac{ \left \{\prod_{n=1}^N \mathcal{N}(x_n | \mu, \lambda^{-1}) \right\} \mathcal{N}(\mu | m, \lambda_{\mu}^{-1}) }{ p(\mathbf{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(\mathbf{X})$は$\mu$に影響しないため省略して、比例関係のみに注目する。省略した部分については、最後に正規化することで対応できる。

 この分布の具体的な形状を明らかにしていく。対数をとって指数部分の計算を分かりやすくして、$\mu$に関して整理すると

$$ \begin{align} \ln p(\mu | \mathbf{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\} + \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{1}{2} (x_n - \mu)^2 \lambda \right\} - \frac{1}{2} \ln 2 \pi \lambda_{\mu}^{-1} - \frac{1}{2} (\mu - m)^2 \lambda_{\mu} + \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行目の式を用いている。ただし、$\mu$と関係のない$- \ln p(\mathbf{X})$を$\mathrm{const.}$とおく。
  2. 尤度と事前分布に、それぞれ具体的な確率分布の式を代入する。
  3. それぞれ対数をとる。自然対数の性質より、$\ln \sqrt{\frac{1}{x}} = \ln x^{-\frac{1}{2}} = -\frac{1}{2} \ln x$、$\ln \exp(x) = x$である。
  4. $-\frac{1}{2}$を括り出す。また、$\mu$と関係のない項を$\mathrm{const.}$に含める。
  5. 2乗に関する括弧を展開する。
  6. $\sum_{n=1}^N$に関する括弧を展開し、$\mu$と関係のない項を$\mathrm{const.}$に含める。
  7. $\mu^2$と$\mu$の項をそれぞれまとめて式を整理する。

となる。

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

$$ \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 | \mathbf{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 + \Bigl( \frac{\tilde{m}}{\hat{\lambda}_{\mu}} \Bigr)^2 - \Bigl( \frac{\tilde{m}}{\hat{\lambda}_{\mu}} \Bigr)^2 \right\} + \mathrm{const.} \\ &= - \frac{1}{2} \hat{\lambda}_{\mu} \left\{ \Bigl( \mu - \frac{\tilde{m}}{\hat{\lambda}_{\mu}} \Bigr)^2 - \Bigl( \frac{\tilde{m}}{\hat{\lambda}_{\mu}} \Bigr)^2 \right\} + \mathrm{const.} \\ &= - \frac{1}{2} \Bigl( \mu - \frac{\tilde{m}}{\hat{\lambda}_{\mu}} \Bigr)^2 \hat{\lambda}_{\mu} + \frac{\tilde{m}^2}{2 \hat{\lambda}_{\mu}} + \mathrm{const.} \end{aligned} $$

【途中式の途中式】

  1. $\hat{\lambda}_{\mu}$を括り出す。
  2. 括弧内に$(\frac{\tilde{m}}{\hat{\lambda}_{\mu}})^2 - (\frac{\tilde{m}}{\hat{\lambda}_{\mu}})^2 = 0$を加える。
  3. 2乗の括弧でまとめる。
  4. 波括弧を展開する。

 後の項は$\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} $$

とおき

$$ \ln p(\mu | \mathbf{X}) = - \frac{1}{2} (\mu - \hat{m})^2 \hat{\lambda}_{\mu} + \mathrm{const.} $$

さらに$\ln$を外し、$\mathrm{const.}$を正規化項に置き換える(正規化する)と

$$ p(\mu | \mathbf{X}) = \frac{1}{\sqrt{2 \pi \hat{\lambda}_{\mu}^{-1}}} \exp \left( - \frac{1}{2} (\mu - \hat{m})^2 \hat{\lambda}_{\mu} \right) = \mathcal{N}(\mu | \hat{m}, \hat{\lambda}_{\mu}^{-1}) \tag{3.51} $$

事後分布は式の形状から、平均$\hat{m}$、精度$\hat{\lambda}_{\mu}$の1次元ガウス分布であることが分かる。
 また、式(3.53)、(3.54)が超パラメータ$\hat{m}, \hat{\lambda}_{\mu}$の計算式(更新式)である。

予測分布の計算

 続いて、1次元ガウス分布に従う未観測のデータ$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} $$

が成り立つ。この式の両辺の対数をとり

$$ \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} $$

となる。$x_{*}$に影響しない$\ln p(\mu)$を$\mathrm{const.}$とおいている。
 この式から、具体的な予測分布の式を計算する。

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

$$ \ln p(\mu | x_{*}) = - \frac{1}{2} \Bigl\{ (\mu - m_{*})^2 (\lambda + \lambda_{\mu}) + \ln (\lambda + \lambda_{\mu})^{-1} + \ln 2 \pi \Bigr\} = \ln \mathcal{N}( \mu | m_{*}, (\lambda + \lambda_{\mu})^{-1} ) \tag{3.58} $$

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

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

とおく。

 尤度(3.47)と式(3.58)を式(3.57)に代入して、$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_{*}, (\lambda + \lambda_{\mu})^{-1} \right) + \mathrm{const.} \\ &= - \frac{1}{2} \Bigl\{ (x_{*} - \mu)^2 \lambda + \ln 2 \pi \lambda^{-1} \Bigr\} + \frac{1}{2} \Bigl\{ (\mu - m_{*})^2 (\lambda + \lambda_{\mu}) + \ln 2 \pi (\lambda + \lambda_{\mu})^{-1} \Bigr\} + \mathrm{const.} \\ &= - \frac{1}{2} \left\{ (x_{*} - \mu)^2 \lambda - \Bigl( \mu - \frac{ \lambda x_{*} + \lambda_{\mu} m }{ \lambda + \lambda_{\mu} } \Bigr)^2 (\lambda + \lambda_{\mu}) \right\} + \mathrm{const.} \\ &= - \frac{1}{2} \left\{ \lambda x_{*}^2 - 2 \mu \lambda x_{*} + \mu^2 \lambda - \mu^2 (\lambda + \lambda_{\mu}) + 2 \mu (\lambda x_{*} + \lambda_{\mu} m) - \frac{ (\lambda x_{*} + \lambda_{\mu} m)^2 }{ \lambda + \lambda_{\mu} } \right\} + \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.}$に含めている。

 式(3.60)について

$$ \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_{*} + \Bigl( \frac{ \tilde{\mu}_{*} }{ 2 \lambda_{*} } \Bigr)^2 - \Bigl( \frac{ \tilde{\mu}_{*} }{ 2 \lambda_{*} } \Bigr)^2 \right\} + \mathrm{const.} \\ &= - \frac{1}{2} \lambda_{*} \left\{ \Bigl( x_{*}^2 - \frac{ \tilde{\mu}_{*} }{ 2 \lambda_{*} } \Bigr)^2 - \Bigl( \frac{ \tilde{\mu}_{*} }{ 2 \lambda_{*} } \Bigr)^2 \right\} + \mathrm{const.} \\ &= - \frac{1}{2} \Bigl( x_{*}^2 - \frac{ \tilde{\mu}_{*} }{ 2 \lambda_{*} } \Bigr)^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} $$

とおき

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

さらに$\ln$を外し、$\mathrm{const.}$を正規化項に置き換える(正規化する)と

$$ p(x_{*}) = \frac{1}{\sqrt{2 \pi \lambda_{*}^{-1}}} \exp \left( - \frac{1}{2} (x_{*}^2 - \mu_{*})^2 \lambda_{*} \right) = \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} $$

である。

事後分布による予測分布

 予測分布の計算に事前分布$p(\mu)$を用いることで、観測データによる学習を行っていない予測分布$p(x_{*})$(のパラメータ$\mu_{*}, \lambda_{*}$)を求めた。事後分布$p(\mu | \mathbf{X})$を用いると、同様の手順で観測データ$\mathbf{X}$によって学習した予測分布$p(x_{*} | \mathbf{X})$(のパラメータ)を求められる。

 そこで、$\mu_{*}, \lambda_{*}$を構成する事前分布のパラメータ$m, \lambda_{\mu}$について、事後分布のパラメータ(3.53)、(3.54)に置き換えたものをそれぞれ$\hat{\mu}_{*}, \hat{\lambda}_{*}$とおくと

$$ \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} $$

となり、予測分布

$$ p(x_{*} | \mathbf{X}) = \frac{1}{\sqrt{2 \pi \hat{\lambda}_{*}^{-1}}} \exp\left( - \frac{1}{2} (x_{*} - \hat{\mu}_{*})^2 \hat{\lambda}_{*} \right) = \mathcal{N}(x_{*} | \hat{\mu}_{*}, \hat{\lambda}_{*}) $$

が得られる。
 また、上の式が予測分布のパラメータ$\hat{\mu}_{*}, \hat{\lambda}_{*}$の計算式(更新式)である。

参考文献

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

おわりに

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

2020/03/05:加筆修正しました。
2021/04/04:加筆修正しました。その際にRで実装編と記事を分割しました。これでも多少は読みやすくなったはず。

【次の節の内容】

www.anarchive-beta.com