からっぽのしょこ

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

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

はじめに

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

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

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

【実装編】

www.anarchive-beta.com

www.anarchive-beta.com

【他の節一覧】

www.anarchive-beta.com

【この節の内容】

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

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

 まずは、事前分布と事後分布について確認する。

 パラメータ$\mu,\ \lambda$の事前分布(パラメータの同時分布)$p(\mu, \lambda)$としてガウス・ガンマ分布$\mathrm{NG}(\mu, \lambda | m, \beta, a, b)$を仮定すると、同じ形式の事後分布になることが知られている。

 また、事前分布をガウス・ガンマ分布とすると、尤度の平均パラメータ$\mu$の精度は$\beta \lambda$となる。

$$ \begin{align} p(\mu, \lambda) &= \mathrm{NG}(\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{NG}(\mu, \lambda | m, \beta, a, b) \\ &= \ln \mathcal{N}(\mu | m, (\beta \lambda)^{-1}) + \ln \mathrm{Gam}(\lambda | a, b) \\ &= - \frac{1}{2} \Bigl\{ (\mu - m)^2 \beta \lambda + \ln (\beta \lambda)^{-1} + \ln 2 \pi \Bigr\} + (a - 1) \ln \lambda - b \lambda + \ln C_G(a, b) \end{aligned} $$

である。

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

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

となる。分母の$p(\mathbf{X})$は$\mu,\ \lambda$に影響しないため省略して、比例関係のみに注目する。省略した部分については、最後に正規化することで対応できる。
 また、図3.5のグラフィカルモデルでパラメータの依存関係を確認すると、左辺の事後分布は

$$ p(\mu, \lambda | \mathbf{X}) = p(\mu | \lambda, \mathbf{X}) p(\lambda | \mathbf{X}) \tag{2} $$

と分解できる。

・平均の事後分布の計算

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

 $\mathbf{X}$が与えられた下での$\mu$の事後分布は、式(1)、(2)より

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

となる。$\mathrm{Gam}(\lambda | a, b)$と$p(\lambda | \mathbf{X})$も$\mu$に影響しないため省く。
 $p(\mu | \lambda, \mathbf{X})$と3.3.1項「3.3.1:1次元ガウス分布の学習と予測:平均が未知の場合【緑ベイズ入門のノート】 - からっぽのしょこ」の式(3.49)を比べると、事前分布の精度パラメータが$\lambda_{\mu}$から$\beta \lambda$に置き換わっただけで同じ形状である。よって同様の手順で求められるので、パラメータを置き換えればよい。

 平均が未知の場合の$\mu$の事後分布の精度パラメータ

$$ \hat{\lambda}_{\mu} = N \lambda + \lambda_{\mu} \tag{3.53} $$

について、$\lambda_{\mu}$を$\beta \lambda$に置き換えると

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

精度パラメータの係数$\hat{\beta}$の計算式(更新式)が得られる。
 $\mu$の事後分布の平均パラメータ

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

についても同様に置き換えると

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

平均パラメータ$\hat{m}$の計算式(更新式)が得られる。

 従って、$\mu$の事後分布は

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

平均$\hat{m}$、精度$\hat{\beta} \lambda$の1次元ガウス分布となることが分かる。

・精度の事後分布の計算

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

 $\mathbf{X}$が与えられた下での$\lambda$の事後分布も、式(1)、(2)より

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

となる。

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

$$ \begin{align} \ln p(\lambda | \mathbf{X}) &= \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) - \ln \mathcal{N}(\mu | \hat{m}, (\hat{\beta} \lambda)^{-1}) + \mathrm{const.} \\ &= - \frac{1}{2} \sum_{n=1}^N \Bigl\{ (x_n - \mu)^2 \lambda + \ln \lambda^{-1} + \ln 2 \pi \Bigr\} \\ &\qquad - \frac{1}{2} \Bigl\{ (\mu - m)^2 \beta \lambda + \ln (\beta \lambda)^{-1} + \ln 2 \pi \Bigr\} \\ &\qquad + (a - 1) \ln \lambda - b \lambda + \ln C_G(a, b) \\ &\qquad + \frac{1}{2} \Bigl\{ (\mu - \hat{m})^2 \hat{\beta} \lambda + \ln (\hat{\beta} \lambda)^{-1} + \ln 2 \pi \Bigr\} + \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 + \frac{N}{2} \ln \lambda \\ &\qquad - \frac{1}{2} \mu^2 \beta \lambda + \mu m \beta \lambda - \frac{1}{2} m^2 \beta \lambda + \frac{1}{2} \ln \beta + \frac{1}{2} \ln \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}{\hat{\beta}} \left( \sum_{n=1}^N x_n + \beta m \right) \hat{\beta} \lambda + \frac{1}{2} \hat{m}^2 \hat{\beta} \lambda - \frac{1}{2} \ln \hat{\beta} - \frac{1}{2} \ln \lambda + \mathrm{const.} \\ &= - \frac{1}{2} \sum_{n=1}^N x_n^2 \lambda + \frac{N}{2} \ln \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.}$に含めている。

 式(3.86)について

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

とおき

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

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

・予測分布の計算

 最後に、未観測のデータ$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_{*}) } $$

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

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

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

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

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

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

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

とおく。

 尤度(3.80)と式(3.91)を式(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}) - \ln \mathcal{N}(\mu | m_{*}, \{(1 + \beta) \lambda\}^{-1}) - \ln \mathrm{Gam} \left( \lambda \middle| \frac{1}{2} + a, b_{*} \right) + \mathrm{const.} \\ &= - \frac{1}{2} \Bigl\{ (x_{*} - \mu)^2 \lambda + \ln \lambda^{-1} + \ln 2 \pi \Bigr\} \\ &\qquad + \frac{1}{2} \left[ \Bigl( \mu - \frac{ x_{*} + \beta m }{ 1 + \beta } \Bigr)^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} \mu^2 (1 + \beta) \lambda - \mu (x_{*} + \beta m) \lambda + \frac{ (x_{*} + \beta m)^2 \lambda }{ 2 (1 + \beta) } \\ &\qquad + \frac{ (x_{*} - m)^2 \beta \lambda }{ 2 (1 + \beta) } + \frac{ b \beta \lambda }{ 2 (1 + \beta) } - \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_{*} - \mu m \beta \lambda + \frac{ \lambda x_{*}^2 + 2 m \beta \lambda x_{*} + m^2 \beta^2 \lambda }{ 2 (1 + \beta) } \\ &\qquad + \frac{ \beta \lambda x_{*}^2 - 2 m \beta \lambda x_{*} + m^2 \beta \lambda }{ 2 (1 + \beta) } - \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{1 + 2a}{2} \ln \left\{ 1 + \frac{\beta}{2 (1 + \beta) b} (x_{*} - m)^2 \right\} + \mathrm{const.} \tag{3.75} \end{align} $$

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

 式(3.75)について、$\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} $$

とおき

$$ \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{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.94} $$

予測分布は式の形状から、パラメータ$\mu_s,\ \lambda_s,\ \nu_s$を持つスチューデントのt分布となることが分かる。

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

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

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

となり、予測分布

$$ p(x_{*} | \mathbf{X}) = \mathrm{St}(x_{*} | \hat{\mu}_s, \hat{\lambda}_s, \hat{\nu}_s) $$

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

参考文献

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

おわりに

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

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

【次節の内容】

www.anarchive-beta.com

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

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

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

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

 現在表記・表現の統一中です。

 途中式が非常にゴチャゴチャしていますが、プラスマイナスで打ち消されて項が綺麗に消えていくのは気持ちいいですよ。いやうんたぶん。