からっぽのしょこ

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

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

はじめに

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

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

【実装編】

www.anarchive-beta.com

www.anarchive-beta.com

【前の節の内容】

www.anarchive-beta.com

【他の節一覧】

www.anarchive-beta.com

【この節の内容】

3.3.2 精度が未知の場合

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

尤度関数の確認

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

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

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

と分解できる。

事後分布の計算

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

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

$$ \begin{align} p(\lambda | \mathbf{X}) &= \frac{ p(\mathbf{X} | \lambda) p(\lambda) }{ p(\mathbf{X}) } \\ &= \frac{ \left\{ \prod_{n=1}^N p(x_n | \lambda) \right\} p(\lambda) }{ p(\mathbf{X}) } \\ &= \frac{ \left\{ \prod_{n=1}^N \mathcal{N}(x_n | \mu, \lambda^{-1}) \right\} \mathrm{Gam}(\lambda | a, b) }{ p(\mathbf{X}) } \\ &\propto \left\{ \prod_{n=1}^N \mathcal{N}(x_n | \mu, \lambda^{-1}) \right\} \mathrm{Gam}(\lambda | a, b) \tag{3.66} \end{align} $$

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

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

$$ \begin{align} \ln p(\lambda | \mathbf{X}) &= \sum_{n=1}^N \ln \mathcal{N}(x_n | \mu, \lambda^{-1}) + \ln \mathrm{Gam}(\lambda | a, b) + \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 \Bigl\{ C_G(a, b) \lambda^{a-1} \exp(- b \lambda) \Bigr\} + \mathrm{const.} \\ &= \sum_{n=1}^N \left\{ - \frac{1}{2} \ln 2 \pi - \frac{1}{2} \ln \lambda^{-1} - \frac{1}{2} (x_n - \mu)^2 \lambda \right\} \\ &\qquad + \ln C_G(a, b) + (a - 1) \ln \lambda - b \lambda + \mathrm{const.} \\ &= \frac{N}{2} \ln \lambda - \frac{1}{2} \sum_{n=1}^N (x_n - \mu)^2 \lambda + (a - 1) \ln \lambda - b \lambda + \mathrm{const.} \\ &= \left( \frac{N}{2} + a - 1 \right) \ln \lambda - \left\{ \frac{1}{2} \sum_{n=1}^N (x_n - \mu)^2 + b \right\} \lambda + \mathrm{const.} \tag{3.67} \end{align} $$

【途中式の途中式】

  1. 式(3.66)の下から2行目の式を用いている。ただし、$\lambda$と関係のない$- \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. $\sum_{n=1}^N$に関する括弧を展開し、$\lambda$と関係のない項を$\mathrm{const.}$に含める。
  5. $\ln \lambda$と$\lambda$の項をそれぞれまとめて式を整理する。

となる。

 式(3.67)について

$$ \begin{aligned} \hat{a} &= \frac{N}{2} + a \\ \hat{b} &= \frac{1}{2} \sum_{n=1}^N (x_n - \mu)^2 + b \end{aligned} \tag{3.69} $$

とおき

$$ \ln p(\lambda | \mathbf{X}) = (\hat{a} - 1) \ln \lambda - \hat{b} \lambda + \mathrm{const.} $$

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

$$ p(\lambda | \mathbf{X}) = C_G(\hat{a}, \hat{b}) \lambda^{\hat{a} - 1} \exp(- \hat{b} \lambda) = \mathrm{Gam}(\lambda | \hat{a}, \hat{b}) \tag{3.68} $$

事後分布は式の形状から、パラメータ$\hat{a}, \hat{b}$を持つガンマ分布であることが分かる。
 また、式(3.69)が超パラメータ$\hat{a}, \hat{b}$の計算式(更新式)である。

予測分布の計算

 続いて、1次元ガウス分布に従う未観測のデータ$x_{*}$に対する予測分布を導出する。

事前分布による予測分布

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

$$ \begin{align} p(x_{*}) &= \int p(x_{*} | \lambda) p(\lambda) d\lambda \\ &= \int \mathcal{N}(x_{*} | \mu, \lambda^{-1}) \mathrm{Gam}(\lambda | a, b) d\lambda \tag{3.70} \end{align} $$

 3.3.1項のときと同様に、この計算を直接するのではなく予測分布$p(x_{*})$と事前分布$p(\lambda)$の関係を考える。ベイズの定理より

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

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

$$ \ln p(\lambda | x_{*}) = \ln p(x_{*} | \lambda) + \ln p(\lambda) - \ln p(x_{*}) $$

予測分布に関して整理すると

$$ \ln p(x_{*}) = \ln p(x_{*} | \lambda) - \ln p(\lambda | x_{*}) + \mathrm{const.} \tag{3.72} $$

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

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

$$ \ln p(\lambda | x_{*}) = \left( \frac{1}{2} + a - 1 \right) \ln \lambda - b_{*} \lambda + \ln C_G \left( \frac{1}{2} + a, b_{*} \right) = \ln \mathrm{Gam} \left( \lambda \middle| \frac{1}{2} + a, b_{*} \right) \tag{3.73} $$

パラメータ$\frac{1}{2} + a, b_{*}$を持つガンマ分布となる。ただし

$$ b_{*} = \frac{1}{2} (x_{*} - \mu)^2 + b \tag{3.74} $$

とおく。

 尤度(3.64)と式(3.73)を式(3.72)に代入する。

$$ \begin{align} \ln p(x_{*}) &= \ln p(x_{*} | \lambda) - \ln p(\lambda | x_{*}) + \mathrm{const.} \tag{3.72}\\ &= \ln \mathcal{N}(x_{*} | \mu, \lambda^{-1}) - \ln \mathrm{Gam} \left( \lambda \middle| \frac{1}{2} + a, b_{*} \right) + \mathrm{const.} \\ &= \ln \left\{ \frac{1}{\sqrt{2 \pi \lambda^{-1}}} \exp \left( - \frac{(x_{*} - \mu)^2}{2 \lambda^{-1}} \right) \right\} \\ &\qquad - \ln \left\{ C_G \left( \frac{1}{2} + a, b_{*} \right) \lambda^{\frac{1}{2}+a-1} \exp(- b_{*} \lambda) \right\} + \mathrm{const.} \\ &= - \frac{1}{2} \ln 2 \pi \lambda^{-1} - \frac{1}{2} (x_{*} - \mu)^2 \lambda \\ &\qquad - \ln C_G \left( \frac{1}{2} + a, b_{*} \right) - \left( \frac{1}{2} + a - 1 \right) \ln \lambda + b_{*} \lambda + \mathrm{const.} \end{align} $$

 正規化項$C_G(\cdot)$と$b_{*}$を展開して、$x_{*}$と関係のない項を$\mathrm{const.}$にまとめて式を整理すると

$$ \begin{align} \ln p(x_{*}) &= - \frac{1}{2} (x_{*} - \mu)^2 \lambda \\ &\qquad - \left( \frac{1}{2} + a \right) \ln \left\{ \frac{1}{2} (x_{*} - \mu)^2 + b \right\} + \ln \Gamma \left( \frac{1}{2} + a \right) + \left\{ \frac{1}{2} (x_{*} - \mu)^2 + b \right\} \lambda + \mathrm{const.} \\ &= - \left( \frac{1}{2} + a \right) \ln \left\{ \frac{1}{2 b} (x_{*} - \mu)^2 + 1 \right\} b + \mathrm{const.} \\ &= - \left( \frac{1}{2} + a \right) \ln \left\{ \frac{1}{2 b} (x_{*} - \mu)^2 + 1 \right\} - \left( \frac{1}{2} + a \right) \ln b + \mathrm{const.} \\ &= - \frac{2 a + 1}{2} \ln \left\{ 1 + \frac{1}{2 b} (x - \mu)^2 \right\} + \mathrm{const.} \tag{3.75} \end{align} $$

【途中式の途中式】

  1. $\frac{1}{2} (x_{*} - \mu)^2 \lambda$の項が打ち消される。また、括弧から$b$を括り出す。
  2. $a \ln x y = a (\ln x + \ln y) = a \ln x + a \ln y$の変形を行う。
  3. 式を整理する。

となる。この式の形状の確率分布はスチューデントのt分布である。

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

$$ \begin{aligned} \mu_s &= \mu \\ \lambda_s &= \frac{a}{b} \\ \nu_s &= 2 a \end{aligned} \tag{3.79} $$

とおき

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

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

$$ p(x_{*}) = \frac{ \Gamma \Bigl(\frac{\nu_s + 1}{2} \Bigr) }{ \Gamma \Bigl(\frac{\nu_s}{2} \Bigr) } \Bigl( \frac{\lambda_s}{\pi \nu_s} \Bigr)^{\frac{1}{2}} \Bigl\{ 1 + \frac{\lambda_s}{\nu_s} (x_{*} - \mu_s)^2 \Bigr\}^{-\frac{\nu_s+1}{2}} = \mathrm{St}(x_{*} | \mu_s, \lambda_s, \nu_s) \tag{3.78} $$

予測分布は式の形状から、パラメータ$\mu_s, \lambda_s, \nu_s$を持つ1次元スチューデントのt分布となることが分かる。
 t分布については「1次元ガウス分布の定義式 - からっぽのしょこ」を参照のこと。

事後分布による予測分布

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

 そこで、$\lambda_s, \nu_s$を構成する事前分布のパラメータ$a, b$について、事後分布のパラメータ(3.69)に置き換えたものをそれぞれ$\hat{\lambda}_s, \hat{\nu}_s$とおくと

$$ \begin{aligned} \hat{\lambda}_s &= \frac{\hat{a}}{\hat{b}} \\ &= \frac{ \frac{N}{2} + a }{ \frac{1}{2} \sum_{n=1}^N (x_n - \mu)^2 + b } \\ &= \frac{ N + 2 a }{ \sum_{n=1}^N (x_n - \mu)^2 + 2 b } \end{aligned} $$

また

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

となり、予測分布

$$ p(x_{*} | \mathbf{X}) = \frac{ \Gamma \Bigl(\frac{\hat{\nu}_s + 1}{2} \Bigr) }{ \Gamma \Bigl(\frac{\hat{\nu}_s}{2} \Bigr) } \Bigl( \frac{\hat{\lambda}_s}{\pi \hat{\nu}_s} \Bigr)^{\frac{1}{2}} \Bigl\{ 1 + \frac{\hat{\lambda}_s}{\hat{\nu}_s} (x_{*} - \mu_s)^2 \Bigr\}^{-\frac{\hat{\nu}_s+1}{2}} = \mathrm{St}(x_{*} | \mu_s, \hat{\lambda}_s, \hat{\nu}_s) $$

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

参考文献

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

おわりに

 式(3.75)を導出できない。難しい訳ではないのですができませんでした…。どこがおかしいのか分かる方はぜひ教えてください。(アドバイスいただいて解けました!説明文を勘違いしてました。ありがとうございます!)
 ちなみに、次でも同様に式を整理できない部分があり現在止まっております。(こっちは丁寧に書き直したらできました。)

 それとは別に予測分布の方もRでやってみたいと思うので、次の項を読み(書き)終えたら一旦全編修正作業に移る予定です。(追加しました!)

2020/03/05:加筆修正しました。
2021/04/04:加筆修正しました。その際にRで実装編と記事を分割しました。

【次の節の内容】

www.anarchive-beta.com