はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」によって、「数式」と「プログラム」から理解するのが目標です。
省略してある内容等ありますので、本とあわせて読んでください。
この記事では、「尤度関数をポアソン分布」、「事前分布をガンマ分布」とした場合の「パラメータの事後分布」と「未観測値の予測分布」を導出します。
【実装編】
www.anarchive-beta.com
www.anarchive-beta.com
【前の節の内容】
www.anarchive-beta.com
【他の節の内容】
www.anarchive-beta.com
【この節の内容】
3.2.3 ポアソン分布の学習と予測
尤度関数をポアソン分布(Poisson Distribution)、事前分布をガンマ分布(Gamma Distribution)とするモデルに対するベイズ推論を導出する。
ポアソン分布については「【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
