はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」によって、「数式」と「プログラム」から理解するのが目標です。
省略してある内容等ありますので、本とあわせて読んでください。
この記事では、「尤度関数をベルヌーイ分布」、「事前分布をベータ分布」とした場合の「パラメータの事後分布」と「未観測値の予測分布」を導出します。
【実装編】
www.anarchive-beta.com
www.anarchive-beta.com
【他の節の内容】
www.anarchive-beta.com
【この節の内容】
3.2.1 ベルヌーイ分布の学習と予測
尤度関数をベルヌーイ分布、事前分布をベータ分布とするモデルに対するベイズ推論を導出する。ベルヌーイ分布については「ベルヌーイ分布の定義式 - からっぽのしょこ」、ベータ分布については「1.2.3：ベータ分布【『トピックモデル』の勉強ノート】 - からっぽのしょこ」を参照のこと。
尤度関数の確認
ベルヌーイ分布に従うと仮定する$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}
$$
と分解できる。
各データ$x_n$は、0か1の2値をとる。また$\mu$は、$x_n = 1$となる確率を表すパラメータで、$0 \leq \mu \leq 1$である。
事後分布の計算
次に、パラメータ$\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 \mathrm{Bern}(x_n | \mu)
\right\}
\mathrm{Beta}(\mu | a, b)
}{
p(\mathbf{X})
}
\\
&\propto
\left\{
\prod_{n=1}^N \mathrm{Bern}(x_n | \mu)
\right\}
\mathrm{Beta}(\mu | a, b)
\tag{3.12}
\end{align}
$$
となる。分母の$p(\mathbf{X})$は$\mu$に影響しないため省略して、比例関係のみに注目する。省略した部分については、最後に正規化することで対応できる。
この分布の具体的な形状を明らかにしていく。対数をとって指数部分の計算を分かりやすくして、$\mu$に関して整理すると
$$
\begin{align}
\ln p(\mu | \mathbf{X})
&= \ln \prod_{n=1}^N
\mathrm{Bern}(x_n | \mu)
+ \ln \mathrm{Beta}(\mu | a, b)
- \ln p(\mathbf{X})
\\
&= \sum_{n=1}^N
\ln \mathrm{Bern}(x_n | \mu)
+ \ln \mathrm{Beta}(\mu | a, b)
+ \mathrm{const.}
\\
&= \sum_{n=1}^N \ln \left\{
\mu^{x_n}
(1 - \mu)^{1-x_n}
\right\}
+ \ln \left\{
C_B(a, b)
\mu^{a-1}
(1 - \mu)^{b-1}
\right\}
+ \mathrm{const.}
\\
&= \sum_{n=1}^N \left\{
x_n \ln \mu
+ (1 - x_n) \ln(1 - \mu)
\right\}
+ \ln C_B(a, b)
+ (a - 1) \ln \mu
+ (b - 1) \ln (1 - \mu)
+ \mathrm{const.}
\\
&= \sum_{n=1}^N
x_n \ln \mu
+ \left(N - \sum_{n=1}^N x_n\right) \ln(1 - \mu)
+ (a - 1) \ln \mu
+ (b - 1) \ln (1 - \mu)
+ \mathrm{const.}
\\
&= \left(
\sum_{n=1}^N x_n + a - 1
\right)
\ln \mu
+ \left(
N - \sum_{n=1}^N x_n + b - 1
\right)
\ln (1 - \mu)
+ \mathrm{const.}
\tag{3.13}
\end{align}
$$
【途中式の途中式】
- 式(3.12)の下から2行目の関係を用いている。
- $\mu$に影響しない$- \ln p(\mathbf{X})$を$\mathrm{const.}$とおく。
総乗部分の対数をとると、自然対数の性質$\ln x y = \ln x + \ln y$より
$$
\begin{aligned}
\ln p(\mathbf{X} | \mu)
&= \ln \prod_{n=1}^N p(x_n | \mu)
\\
&= \ln \Bigl\{
p(x_1 | \mu) * p(x_2 | \mu) * \cdots * p(x_N | \mu)
\Bigr\}
\\
&= \ln p(x_1 | \mu) + \ln p(x_2 | \mu) + \cdots + \ln p(x_N | \mu)
\\
&= \sum_{n=1}^N \ln p(x_n | \mu)
\end{aligned}
$$
総和になる。
- 尤度と事前分布に、それぞれ具体的な確率分布の式を代入する。
- それぞれ対数をとる。$\ln x^a = a \ln x$である。
- $\sum_{n=1}^N$に関する括弧を展開し、また$\mu$に影響しない$\ln C_B(a, b)$を$\mathrm{const.}$にまとめる。
- $\mu$と$1 - \mu$の項をそれぞれまとめて式を整理する。
となる。
式(3.13)について
$$
\begin{aligned}
\hat{a}
&= \sum_{n=1}^N x_n + a
\\
\hat{b}
&= N - \sum_{n=1}^N x_n + b
\end{aligned}
\tag{3.15}
$$
とおき
$$
\ln p(\mu | \mathbf{X})
= (\hat{a} - 1)
\ln \mu
+ (\hat{b} - 1)
\ln (1 - \mu)
+ \mathrm{const.}
$$
さらに$\ln$を外して、$\mathrm{const.}$を正規化項に置き換えると(正規化する)と
$$
p(\mu | \mathbf{X})
= C_B(\hat{a}, \hat{b})
\mu^{\hat{a} - 1}
(1 - \mu)^{\hat{b} - 1}
= \mathrm{Beta}(\mu | \hat{a}, \hat{b})
\tag{3.14}
$$
事後分布は式の形状から、パラメータ$\hat{a},\ \hat{b}$を持つベータ分布であることが分かる。
また、式(3.15)が超パラメータ$\hat{a},\ \hat{b}$の計算式(更新式)である。
予測分布の計算
続いて、ベルヌーイ分布に従う未観測のデータ$x_{*}$に対する予測分布を導出する。
事前分布による予測分布
事前分布(観測データによる学習を行っていない分布)$p(\mu)$を用いて、パラメータ$\mu$を周辺化することで予測分布$p(x_{*})$となる。
$$
\begin{align}
p(x_{*})
&= \int
p(x_{*} | \mu)
p(\mu)
d\mu
\\
&= \int
\mathrm{Bern}(x_{*} | \mu)
\mathrm{Beta}(\mu | a, b)
d\mu
\\
&= \int
\mu^{x_{*}}
(1 - \mu)^{1-x_{*}}
C_B(a, b)
\mu^{a-1}
(1 - \mu)^{b-1}
d\mu
\\
&= C_B(a, b)
\int
\mu^{x_{*}+a-1}
(1 - \mu)^{1-x_{*}+b-1}
d\mu
\tag{3.16}
\end{align}
$$
$\mu$と$1 - \mu$の項をまとめて積分部分に注目すると、パラメータ$x_{*} + a,\ 1 - x_{*} + b$を持つ正規化項のないベータ分布の形をしている。
これはベータ分布の正規化項の逆数に変形できる(そもそもこの部分の逆数が正規化項である。詳しくは「1.2.3：ベータ分布【『トピックモデル』の勉強ノート】 - からっぽのしょこ」を参照のこと)。
$$
\int
\mu^{x_{*}+a-1}
(1 - \mu)^{1-x_{*}+b-1}
d\mu
= \frac{1}{C_B(x_{*} + a, 1 - x_{*} + b)}
\tag{3.17}
$$
これを式(3.16)に代入すると、予測分布は
$$
p(x_{*})
= \frac{C_B(a, b)}{C_B(x_{*} + a, x_{*} + b)}
$$
となる。
さらにベータ分布の正規化項(2.42)を用いると
$$
\begin{align}
p(x_{*})
&= C_B(a, b)
\frac{1}{C_B(x_{*} + a, x_{*} + b)}
\\
&= \frac{
\Gamma(a + b)
}{
\Gamma(a)
\Gamma(b)
}
\frac{
\Gamma(x_{*} + a)
\Gamma(1 - x_{*} + b)
}{
\Gamma(x_{*} + a + 1 - x_{*} + b)
}
\\
&= \frac{
\Gamma(a + b)
}{
\Gamma(a)
\Gamma(b)
}
\frac{
\Gamma(x_{*} + a)
\Gamma(1 - x_{*} + b)
}{
\Gamma(a + b + 1)
}
\\
&= \frac{
\Gamma(a + b)
}{
\Gamma(a)
\Gamma(b)
}
\frac{
\Gamma(x_{*} + a)
\Gamma(1 - x_{*} + b)
}{
(a + b) \Gamma(a + b)
}
\\
&= \frac{
\Gamma(x_{*} + a)
\Gamma(1 - x_{*} + b)
}{
\Gamma(a)
\Gamma(b)
(a + b)
}
\tag{3.18}
\end{align}
$$
となる。ガンマ関数の性質$\Gamma(x + 1) = x \Gamma(x)$を用いている。
$x_{*}$は0か1しかとらないため、場合分けしてこの式を整理する。$x_{*} = 1$のとき
$$
\begin{align}
p(x_{*} = 1)
&= \frac{
\Gamma(1 + a)
\Gamma(1 - 1 + b)
}{
\Gamma(a)
\Gamma(b)
(a + b)
}
\tag{3.18'}\\
&= \frac{
a \Gamma(a)
\Gamma(b)
}{
\Gamma(a)
\Gamma(b)
(a + b)
}
\\
&= \frac{a}{a + b}
\tag{3.19}
\end{align}
$$
となる。また$x_{*} = 0$のとき
$$
\begin{align}
p(x_{*} = 0)
&= \frac{
\Gamma(0 + a)
\Gamma(1 - 0 + b)
}{
\Gamma(a)
\Gamma(b)
(a + b)
}
\tag{3.18'}\\
&= \frac{
\Gamma(a)
b \Gamma(b)
}{
\Gamma(a)
\Gamma(b)
(a + b)
}
\\
&= \frac{b}{a + b}
\tag{3.20}
\end{align}
$$
となる。これは
$$
\frac{b}{a + b}
= 1 - \frac{a}{a + b}
$$
と置き換えられる。
従って、$x^0 = 1$であることを利用して式(3.19)と式(3.20)をまとめると、式(3.18)は
$$
\begin{aligned}
p(x_{*})
&= \Bigl( \frac{a}{a + b} \Bigr)^{x_{*}}
\Bigl( \frac{b}{a + b} \Bigr)^{1-x_{*}}
\\
&= \Bigl( \frac{a}{a + b} \Bigr)^{x_{*}}
\Bigl( 1 - \frac{a}{a + b} \Bigr)^{1-x_{*}}
\end{aligned}
$$
と書き換えられる。($x_{*} = 0$のとき前の項が1となり、$x_{*} = 1$のとき後ろの項が1となるので、それぞれ式(3.20)と式(3.19)が成り立つ。)
この式について
$$
\mu_{*}
= \frac{a}{a + b}
$$
とおくと
$$
p(x_{*})
= \mu_{*}^{x_{*}}
(1 - \mu_{*})^{1-x_{*}}
= \mathrm{Bern}(x_{*} | \mu_{*})
\tag{3.21}
$$
予測分布は式の形状から、パラメータ$\mu_{*}$を持つベルヌーイ分布であることが分かる。ちなみに、$\frac{a}{a + b}$はベータ分布$\mathrm{Beta}(\mu | a, b)$の期待値$\mathbb{E}[\mu]$である。
事後分布による予測分布
予測分布の計算に事前分布$p(\mu)$を用いることで、観測データによる学習を行っていない予測分布$p(x_{*})$(のパラメータ$\mu_{*}$)を求めた。事後分布$p(\mu | \mathbf{X})$を用いると、同様の手順で観測データ$\mathbf{X}$によって学習した予測分布$p(x_{*} | \mathbf{X})$(のパラメータ)を求められる。
そこで、$\mu_{*}$を構成する事前分布のパラメータ$a,\ b$について、事後分布のパラメータ(3.15)に置き換えたものを$\hat{\mu}_{*}$とおくと
$$
\begin{aligned}
\hat{\mu}_{*}
&= \frac{\hat{a}}{\hat{a} + \hat{b}}
\\
&= \frac{
\sum_{n=1}^N x_n + a
}{
\sum_{n=1}^N x_n + a
+ N - \sum_{n=1}^N x_n + b
}
\\
&= \frac{\sum_{n=1}^N x_n + a}{N + a + b}
\end{aligned}
$$
となり、予測分布
$$
p(x_{*} | \mathbf{X})
= \hat{\mu}_{*}^{x_{*}}
(1 - \hat{\mu}_{*})^{1-x_{*}}
= \mathrm{Bern}(x_{*} | \hat{\mu}_{*})
\tag{3.22}
$$
が得られる。
また上の式が、予測分布のパラメータ$\hat{\mu}_{*}$の計算式(更新式)である。
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
ずっと辞書的に使っていた須山ベイズ本の精読を始めました。
実装編とで口調が変わることに深い意味はありません。正直ですます調の方が書きやすいです。
暫く続きます。宜しくお願いします。
2020年2月19日は、モーニング娘。'20森戸知沙希ちゃん二十歳の誕生日！！！
Enjoy!
- 2020/03/04：加筆修正しました。
- 2020/02/19：加筆修正して、Rで実装編とは記事を分割しました。ブログ投稿から丁度一年後に更新することになるとは。
【次節の内容】
www.anarchive-beta.com
