はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」によって、「数式」と「プログラム」から理解するのが目標です。
省略してある内容等ありますので、本とあわせて読んでください。
この記事では、「尤度関数をポアソン分布」、「事前分布をガンマ分布」とした場合の「パラメータの事後分布」と「未観測値の予測分布」を導出します。
【実装編】
www.anarchive-beta.com
www.anarchive-beta.com
【他の節の内容】
www.anarchive-beta.com
【この節の内容】
3.2.3 ポアソン分布の学習と予測
尤度関数をポアソン分布、事前分布をガンマ分布とするモデルに対するベイズ推論を導出する。ポアソン分布については「【R】ポアソン分布の作図 - からっぽのしょこ」、ガンマ分布については「【R】ガンマ分布の作図 - からっぽのしょこ」を参照のこと。
尤度関数の確認
ポアソン分布に従うと仮定する$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 | a, b)
}{
p(\mathbf{X})
}
\\
&= \frac{
\left\{ \prod_{n=1}^N
p(x_n | \lambda)
\right\}
p(\lambda | a, b)
}{
p(\mathbf{X})
}
\\
&= \frac{
\left\{ \prod_{n=1}^N
\mathrm{Poi}(x_n | \lambda)
\right\}
\mathrm{Gam}(\lambda | a, b)
}{
p(\mathbf{X})
}
\\
&\propto
\left\{ \prod_{n=1}^N
\mathrm{Poi}(x_n | \lambda)
\right\}
\mathrm{Gam}(\lambda | a, b)
\tag{3.35}
\end{align}
$$
となる。分母の$p(\mathbf{X})$は$\lambda$に影響しないため省略して、比例関係のみに注目する。省略した部分については、最後に正規化することで対応できる。
この分布の具体的な形状を明らかにしていく。対数をとって指数部分の計算を分かりやすくして、$\lambda$に関して整理すると
$$
\begin{align}
\ln p(\lambda | \mathbf{X})
&= \ln \prod_{n=1}^N
\mathrm{Poi}(x_n | \lambda)
+ \ln \mathrm{Gam}(\lambda | a, b)
- \ln p(\mathbf{X})
\\
&= \sum_{n=1}^N
\ln \mathrm{Poi}(x_n | \lambda)
+ \ln \mathrm{Gam}(\lambda | a, b)
+ \mathrm{const.}
\\
&= \sum_{n=1}^N
\ln \left\{
\frac{\lambda^{x_n}}{x_n!}
\exp(-\lambda)
\right\}
+ \ln \Bigl\{
C_G(a, b)
\lambda^{a-1}
\exp(- b \lambda)
\Bigr\}
+ \mathrm{const.}
\\
&= \sum_{n=1}^N \Bigl\{
x_n
\ln \lambda
- \ln x_n!
- \lambda
\Bigr\}
+ \ln C_G(a, b)
+ (a - 1)
\ln \lambda
- b \lambda
+ \mathrm{const.}
\\
&= \sum_{n=1}^N
x_n
\ln \lambda
- N \lambda
+ (a - 1)
\ln \lambda
- b \lambda
+ \mathrm{const.}
\\
&= \left(
\sum_{n=1}^N x_n
+ a - 1
\right)
\ln \lambda
- (N + b) \lambda
+ \mathrm{const.}
\tag{3.36}
\end{align}
$$
【途中式の途中式】
- 式(3.36)の下から2行目の式を用いている。
- 自然対数の性質より、$\ln x y = \ln x + \ln y$である。また、$\lambda$に影響しない$- \ln p(\mathbf{X})$を$\mathrm{const.}$とおく。
- 尤度と事前分布に、それぞれ具体的な確率分布の式を代入する。
- それぞれ対数をとる。$\ln x^a = a \ln x$、$\ln \frac{x}{y} = \ln x - \ln y$、$\ln \exp(x) = x$である。
- $\lambda$と関係のない$- \sum_{n=1}^N \ln x_n!$と$\ln C_G(a, b)$を$\mathrm{const.}$にまとめる。
- $\ln \lambda$と$\lambda$の項をそれぞれまとめて式を整理する。
となる。
式(3.36)について
$$
\begin{aligned}
\hat{a}
&= \sum_{n=1}^N
x_n
+ a
\\
\hat{b}
&= N + b
\end{aligned}
\tag{3.38}
$$
とおき
$$
\ln p(\lambda | \mathbf{X})
= \hat{a}
\ln \lambda
- \hat{b}
\lambda
+ \mathrm{const.}
\tag{3.36}
$$
さらに$\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.37}
$$
事後分布は式の形状から、パラメータ$\hat{a},\ \hat{b}$を持つガンマ分布であることが分かる。
また、式(3.38)が超パラメータ$\hat{a},\ \hat{b}$の計算式(更新式)である。
予測分布の計算
続いて、ポアソン分布に従う未観測のデータ$x_{*}$に対する予測分布を導出する。
事前分布による予測分布
事前分布(観測データによる学習を行っていない分布)$p(\lambda)$を用いて、パラメータ$\lambda$を周辺化することで予測分布$p(x_{*})$となる。
$$
\begin{align}
p(x_{*})
&= \int
p(x_{*} | \lambda)
p(\lambda)
d\lambda
\\
&= \int
\mathrm{Poi}(x_{*} | \lambda)
\mathrm{Gam}(\lambda | a, b)
d\lambda
\tag{3.39}\\
&= \int
\frac{\lambda^{x_{*}}}{x_{*}!}
\exp(- \lambda)
C_G(a, b)
\lambda^{a-1}
\exp(- b \lambda)
d\lambda
\\
&= \frac{C_G(a, b)}{x_{*}!}
\int
\lambda^{x_{*}+a-1}
\exp \Bigl(
- (1 + b) \lambda
\Bigr)
d\lambda
\tag{3.40}
\end{align}
$$
積分部分に注目すると、パラメータ$x_{*} + a,\ 1 + b$を持つ正規化項のないガンマ分布の形をしている。
これはガンマ分布の正規化項の逆数に変形できる(そもそもこの部分の逆数が正規化項である)。
$$
\int
\lambda^{x_{*}+a-1}
\exp \Bigl(
- (1 + b) \lambda
\Bigr)
d\lambda
= \frac{1}{C_G(x_{*} + a, 1 + b)}
$$
これを式(3.40)に代入すると、予測分布は
$$
p(x_{*})
= \frac{
C_G(a, b)
}{
x_{*}!
C_G(x_{*} + a, 1 + b)
}
\tag{3.42}
$$
となる。
さらにガンマ分布の正規化項(2.57)を用いると
$$
\begin{aligned}
p(x_{*})
&= \frac{C_G(a, b)}{x_{*}!}
\frac{
1
}{
C_G(x_{*} + a, 1 + b)
}
\\
&= \frac{
b^a
}{
x_{*}!
\Gamma(a)
}
\frac{
\Gamma(x_{*} + a)
}{
(1 + b)^{x_{*} + a}
}
\\
&= \frac{
\Gamma(x_{*} + a)
}{
x_{*}!
\Gamma(a)
}
\frac{
b^a
}{
(1 + b)^{a}
}
\frac{1}{(1 + b)^{x_{*}}}
\\
&= \frac{
\Gamma(x_{*} + a)
}{
x_{*}!
\Gamma(a)
}
\Bigl(
\frac{
b
}{
1 + b
}
\Bigr)^{a}
\Bigl(
\frac{1}{1 + b}
\Bigr)^{x_{*}}
\\
&= \frac{
\Gamma(x_{*} + a)
}{
x_{*}!
\Gamma(a)
}
\Bigl(
1
- \frac{
1
}{
1 + b
}
\Bigr)^{a}
\Bigl(
\frac{1}{1 + b}
\Bigr)^{x_{*}}
\end{aligned}
$$
となる。この式の形状の確率分布は負の二項分布である。
この式について
$$
\begin{aligned}
r &= a
\\
p &= \frac{1}{1 + b}
\end{aligned}
\tag{3.44}
$$
とおくと
$$
p(x_{*})
= \frac{
\Gamma(x_{*} + r)
}{
x_{*}!
\Gamma(r)
}
(1 - p)^r
p^{x_{*}}
= \mathrm{NB}(x_{*} | r, p)
\tag{3.43}
$$
予測分布は式の形状から、パラメータ$r,\ p$を持つ負の二項分布であることが分かる。
事後分布による予測分布
予測分布の計算に事前分布$p(\lambda)$を用いることで、観測データによる学習を行っていない予測分布$p(x_{*})$(のパラメータ$r,\ p$)を求めた。事後分布$p(\lambda | \mathbf{X})$を用いると、同様の手順で観測データ$\mathbf{X}$によって学習した予測分布$p(x_{*} | \mathbf{X})$(のパラメータ)を求められる。
そこで、$r,\ p$を構成する事前分布のパラメータ$a,\ b$について、事後分布のパラメータ(3.38)に置き換えたものをそれぞれ$\hat{r},\ \hat{p}$とおくと
$$
\begin{aligned}
\hat{r}
&= \hat{a}
\\
&= \sum_{n=1}^N x_n
+ a
\end{aligned}
$$
また
$$
\begin{aligned}
\hat{p}
&= \frac{1}{1 + \hat{b}}
\\
&= \frac{1}{N + 1 + b}
\end{aligned}
$$
となり、予測分布
$$
p(x_{*} | \mathbf{X})
= \frac{
\Gamma(x_{*} + \hat{r})
}{
x_{*}!
\Gamma(\hat{r})
}
(1 - \hat{p})^{\hat{r}}
\hat{p}^{x_{*}}
= \mathrm{NB}(x_{*} | \hat{r}, \hat{p})
$$
が得られる。
また、上の式が予測分布のパラメータ$\hat{r},\ \hat{p}$の計算式(更新式)である。
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
順調である。
2020/03/05：加筆修正しました。
2021/04/04：加筆修正しました。その際にRで実装編と記事を分割しました。と言ってもこの記事は表記を統一したくらいですかね。
【次節の内容】
www.anarchive-beta.com
