はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。「数式の行間を読んでみた」とそれを「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}
$$
【途中式の途中式】
- 式(3.66)の下から2行目の式を用いている。ただし、$\lambda$と関係のない$- \ln p(\mathbf{X})$を$\mathrm{const.}$とおく。
- 尤度と事前分布に、それぞれ具体的な確率分布の式を代入する。
- それぞれ対数をとる。自然対数の性質より、$\ln \sqrt{\frac{1}{x}} = \ln x^{-\frac{1}{2}} = -\frac{1}{2} \ln x$、$\ln \exp(x) = x$である。
- $\sum_{n=1}^N$に関する括弧を展開し、$\lambda$と関係のない項を$\mathrm{const.}$に含める。
- $\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}
$$
【途中式の途中式】
- $\frac{1}{2} (x_{*} - \mu)^2 \lambda$の項が打ち消される。また、括弧から$b$を括り出す。
- $a \ln x y = a (\ln x + \ln y) = a \ln x + a \ln y$の変形を行う。
- 式を整理する。
となる。この式の形状の確率分布はスチューデントの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
