はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。基本的な内容は「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」になります。「数式」と「プログラム」から理解するのが目標です。
この記事は、4.3.2項の内容です。「観測モデルをポアソン混合分布」、「事前分布をガンマ分布とディリクレ分布」とするポアソン混合モデルの事後分布をギブスサンプリングを用いて推論します。
省略してある内容等ありますので、本とあわせて読んでください。初学者な自分が理解できるレベルまで落として書き下していますので、分かる人にはかなりくどくなっています。同じような立場の人のお役に立てれば幸いです。
【実装編】
www.anarchive-beta.com
www.anarchive-beta.com
【他の節一覧】
www.anarchive-beta.com
【この節の内容】
4.3.2 ギブスサンプリング
ギブスサンプリングを用いて、ポアソン混合モデルの事後分布$p(\mathbf{S}, \boldsymbol{\lambda}, \boldsymbol{\pi} | \mathbf{X})$から潜在変数$\mathbf{S}$とパラメータ$\boldsymbol{\lambda},\ \boldsymbol{\pi}$をサンプルするアルゴリズムを導出する。
・潜在変数の条件付き分布の導出
始めに、潜在変数$\mathbf{S}$をサンプルするための条件付き分布($\mathbf{S}$の事後分布)$p(\mathbf{S} | \mathbf{X}, \boldsymbol{\lambda}, \boldsymbol{\pi})$を求めていく。
観測データ$\mathbf{X}$に加えて、パラメータ$\boldsymbol{\lambda},\ \boldsymbol{\pi}$が観測された(サンプルされた)としたとき、$\mathbf{S}$の条件付き分布は
$$
\begin{align}
p(\mathbf{S} | \mathbf{X}, \boldsymbol{\lambda}, \boldsymbol{\pi})
&= \frac{
p(\mathbf{X}, \mathbf{S}, \boldsymbol{\lambda}, \boldsymbol{\pi})
}{
p(\mathbf{X}, \boldsymbol{\lambda}, \boldsymbol{\pi})
}
\\
&\propto
p(\mathbf{X}, \mathbf{S}, \boldsymbol{\lambda}, \boldsymbol{\pi})
\\
&= p(\mathbf{X} | \mathbf{S}, \boldsymbol{\lambda})
p(\mathbf{S} | \boldsymbol{\pi})
p(\boldsymbol{\lambda})
p(\boldsymbol{\pi})
\\
&\propto
p(\mathbf{X} | \mathbf{S}, \boldsymbol{\lambda})
p(\mathbf{S} | \boldsymbol{\pi})
\\
&= \prod_{n=1}^N
p(x_n | \mathbf{s}_n, \boldsymbol{\lambda})
p(\mathbf{s}_n | \boldsymbol{\pi})
\tag{4.33}\\
&= \prod_{n=1}^N
\left\{
\prod_{k=1}^K \mathrm{Poi}(x_n | \lambda_k)^{s_{n,k}}
\right\}
\mathrm{Cat}(\mathbf{s}_n | \boldsymbol{\pi})
\end{align}
$$
で求められる。4.3.1項「ポアソン混合モデル」で確認した各変数の生成過程(依存関係)に従い項を分解している。また、適宜$\mathbf{S}$に影響しない項を省略して比例関係に注目する。省略した部分については、最後に正規化することで対応できる。
$n$番目の潜在変数(ある1つのデータのクラスタ)$\mathbf{s}_n$の条件付き分布の具体的な形状を明らかにしていく。対数をとって指数部分の計算を分かりやすくすると、前の項は
$$
\begin{align}
\ln p(x_n | \mathbf{s}_n, \boldsymbol{\lambda})
&= \sum_{k=1}^K
s_{n,k}
\ln \mathrm{Poi}(x_n | \lambda_k)
\\
&= \sum_{k=1}^K
s_{n,k}
\ln \frac{\lambda_k^{x_n}}{x_n!}
\exp(- \lambda_k)
\\
&= \sum_{k=1}^K
s_{n,k} (
x_n \ln \lambda_k
- \ln x_n!
- \lambda_k
)
\\
&= \sum_{k=1}^K
s_{n,k} (
x_n \ln \lambda_k
- \lambda_k
)
+ \mathrm{const.}
\tag{4.34}
\end{align}
$$
となる。$\sum_{k=1}^K s_{n,k} = 1$なので、$\sum_{k=1}^K - s_{n,k} \ln x_n! = - \ln x_n!$となり$\mathbf{s}_n$の影響を受けなくなるので$\mathrm{const.}$に含める。
後の項は
$$
\begin{align}
\ln p(\mathbf{s}_n | \boldsymbol{\pi})
&= \ln \mathrm{Cat}(\mathbf{s}_n | \boldsymbol{\pi})
\\
&= \ln \prod_{k=1}^K
\pi_k^{s_{n,k}}
\\
&= \sum_{k=1}^K
s_{n,k}
\ln \pi_k
\tag{4.35}
\end{align}
$$
となる。
よって、式(4.34)と式(4.35)を$n$番目のデータに関する項を取り出した式(4.33)に代入すると
$$
\begin{align}
\ln p(\mathbf{s}_n | \mathbf{X}, \boldsymbol{\lambda}, \boldsymbol{\pi})
&= \ln p(x_n | \mathbf{s}_n, \boldsymbol{\lambda})
+ \ln p(\mathbf{s}_n | \boldsymbol{\pi})
+ \mathrm{const.}
\tag{4.33'}\\
&= \sum_{k=1}^K
s_{n,k} (
x_n \ln \lambda_k - \lambda_k
)
+ \sum_{k=1}^K s_{n,k} \ln \pi_k
+ \mathrm{const.}
\\
&= \sum_{k=1}^K
s_{n,k} (
x_n \ln \lambda_k - \lambda_k + \ln \pi_k
)
+ \mathrm{const.}
\tag{4.36}
\end{align}
$$
となる。適宜$\mathbf{s}_n$に影響しない項を$\mathrm{const.}$にまとめている。
式(4.36)について
$$
\eta_{n,k}
\propto
\exp(
x_n \ln \lambda_k
- \lambda_k
+ \ln \pi_k
)
\tag{4.38}
$$
とおき
$$
\ln p(\mathbf{s}_n | \mathbf{X}, \boldsymbol{\lambda}, \boldsymbol{\pi})
= \sum_{k=1}^K
s_{n,k} \ln \eta_{n,k}
+ \mathrm{const.}
$$
さらに$\ln$を外し、$\sum_{k=1}^K \eta_{n,k} = 1$となるように正規化する($\mathrm{const.}$を正規化項に置き換える)と
$$
p(\mathbf{s}_n | \mathbf{X}, \boldsymbol{\lambda}, \boldsymbol{\pi})
= \prod_{k=1}^K
\eta_{n,k}^{s_{n,k}}
= \mathrm{Cat}(\mathbf{s}_n | \boldsymbol{\eta}_n)
$$
$\mathbf{s}_n$の条件付き分布は、パラメータ$\boldsymbol{\eta}_n = (\eta_{n,1},\eta_{n,2} \cdots, \eta_{n,K})$を持つカテゴリ分布になることが分かる。
・パラメータの条件付き分布の導出
次に、パラメータ$\boldsymbol{\lambda},\ \boldsymbol{\pi}$の(同時)条件付き分布($\boldsymbol{\lambda},\ \boldsymbol{\pi}$の事後分布)$p(\boldsymbol{\lambda}, \boldsymbol{\pi} | \mathbf{X}, \mathbf{S})$から、各パラメータをサンプルするための条件付き分布$p(\boldsymbol{\lambda} | \mathbf{X}, \mathbf{S}),\ p( \boldsymbol{\pi} | \mathbf{X}, \mathbf{S})$を求めていく。
観測データ$\mathbf{X}$に加えて、潜在変数$\mathbf{S}$が観測された(サンプルされた)としたとき、$\boldsymbol{\lambda},\ \boldsymbol{\pi}$の条件付き分布は
$$
\begin{align}
p(\boldsymbol{\lambda}, \boldsymbol{\pi} | \mathbf{X}, \mathbf{S})
&= \frac{
p(\mathbf{X}, \mathbf{S}, \boldsymbol{\lambda}, \boldsymbol{\pi})
}{
p(\mathbf{X}, \mathbf{S})
}
\\
&\propto
p(\mathbf{X}, \mathbf{S}, \boldsymbol{\lambda}, \boldsymbol{\pi})
\\
&= p(\mathbf{X} | \mathbf{S}, \boldsymbol{\lambda})
p(\boldsymbol{\lambda})
p(\mathbf{S} | \boldsymbol{\pi})
p(\boldsymbol{\pi})
\tag{4.39}\\
&= \left\{
\prod_{n=1}^N
p(x_n | \mathbf{s}_n, \boldsymbol{\lambda})
\right\}
\left\{
\prod_{k=1}^K
p(\lambda_k)
\right\}
\left\{
\prod_{n=1}^N
p(\mathbf{s}_n | \boldsymbol{\pi})
\right\}
p(\boldsymbol{\pi})
\\
&= \left\{
\prod_{n=1}^N \prod_{k=1}^K
\mathrm{Poi}(x_n | \lambda_k)^{s_{n,k}}
\right\}
\left\{
\prod_{k=1}^K
\mathrm{Gam}(\lambda_k | a, b)
\right\}
\left\{
\prod_{n=1}^N
\mathrm{Cat}(\mathbf{s}_n | \boldsymbol{\pi})
\right\}
\mathrm{Dir}(\boldsymbol{\pi} | \boldsymbol{\alpha})
\end{align}
$$
で求められる。こちらも生成過程に従い項を分解して、$\boldsymbol{\lambda},\ \boldsymbol{\pi}$に影響しない項を省く。
また、左辺の(同時)条件付き分布は
$$
\begin{aligned}
p(\boldsymbol{\lambda}, \boldsymbol{\pi} | \mathbf{X}, \mathbf{S})
&= p(\boldsymbol{\lambda} | \mathbf{X}, \mathbf{S})
p(\boldsymbol{\pi} | \mathbf{X}, \mathbf{S})
\\
&= \left\{
\prod_{k=1}^K
p(\boldsymbol{\lambda}_k | \mathbf{X}, \mathbf{S})
\right\}
p(\boldsymbol{\pi} | \mathbf{X}, \mathbf{S})
\end{aligned}
$$
と分解できる。
この式を用いて、$\boldsymbol{\lambda},\ \boldsymbol{\pi}$それぞれの条件付き分布の具体的な形状を明らかにしていく。
・観測モデルのパラメータの条件付き分布
式(4.39)の対数をとって、$\boldsymbol{\lambda}$に関して整理する($\boldsymbol{\lambda}$に影響しない項を$\mathrm{const.}$にまとめる)と
$$
\begin{align}
p(\boldsymbol{\lambda} | \mathbf{X}, \mathbf{S})
&= \ln p(\mathbf{X} | \mathbf{S}, \boldsymbol{\lambda})
+ \ln p(\boldsymbol{\lambda})
- \ln p(\boldsymbol{\pi} | \mathbf{X}, \mathbf{S})
+ \mathrm{const.}
\\
&= \sum_{n=1}^N \sum_{k=1}^K
s_{n,k} \ln \mathrm{Poi}(x_n | \lambda_k)
+ \sum_{k=1}^K
\ln \mathrm{Gam}(\lambda_k | a, b)
+ \mathrm{const.}
\\
&= \sum_{k=1}^K \left\{
\sum_{n=1}^N s_{n,k}
\ln \frac{\lambda_k^{x_n}}{x_n!}
\exp(- \lambda_k)
+ \ln C_G(a, b)
\lambda_k^{a-1}
\exp(- b \lambda_k)
\right\}
+ \mathrm{const.}
\\
&= \sum_{k=1}^K \left\{
\sum_{n=1}^N s_{n,k} (
x_n \ln \lambda_k
- \ln x_n!
- \lambda_k
)
+ \ln C_G(a, b)
+ (a - 1) \ln \lambda_k
- b \lambda_k
\right\}
+ \mathrm{const.}
\\
&= \sum_{k=1}^K \left\{
\left(
\sum_{n=1}^N s_{n,k} x_n + a - 1
\right)
\ln \lambda_k
- \left(
\sum_{n=1}^N s_{n,k} + b
\right)
\lambda_k
\right\}
+ \mathrm{const.}
\tag{4.40}
\end{align}
$$
となる。$\ln p(\boldsymbol{\pi} | \mathbf{X}, \mathbf{S})$は左辺から移項したものである。
式(4.40)について
$$
\begin{aligned}
\hat{a}_k
&= \sum_{n=1}^N s_{n,k} x_n + a
\\
\hat{b}_k
&= \sum_{n=1}^N s_{n,k} + b
\end{aligned}
\tag{4.42}
$$
とおき
$$
\ln p(\boldsymbol{\lambda} | \mathbf{X}, \mathbf{S})
= \sum_{k=1}^K \Bigl\{
(\hat{a}_k - 1)
\ln \lambda_k
- \hat{b} \lambda_k
\Bigr\}
+ \mathrm{const.}
$$
さらに$\ln$を外し、$\mathrm{const.}$を正規化項に置き換える(正規化する)と
$$
p(\boldsymbol{\lambda} | \mathbf{X}, \mathbf{S})
= \prod_{k=1}^K
C_G(\hat{a}_k, \hat{b}_k)
\lambda_k^{\hat{a}_k-1}
\exp(- \hat{b}_k \lambda_k)
= \prod_{k=1}^K
\mathrm{Gam}(\lambda_k | \hat{a}_k, \hat{b}_k)
$$
$\boldsymbol{\lambda}_k$の条件付き分布は、パラメータ$\hat{a}_k,\ \hat{b}_k$を持つガンマ分布になることが分かる。
・混合比率の条件付き分布
同様に、対数をとった式(4.39)を$\boldsymbol{\pi}$に関して整理する($\boldsymbol{\pi}$に影響しない項を$\mathrm{const.}$にまとめる)と
$$
\begin{align}
\ln p(\boldsymbol{\pi} | \mathbf{X}, \mathbf{S})
&= \ln p(\mathbf{S} | \boldsymbol{\pi})
+ \ln p(\boldsymbol{\pi})
- \ln p(\boldsymbol{\lambda} | \mathbf{X}, \mathbf{S})
+ \mathrm{const.}
\\
&= \sum_{n=1}^N
\ln \mathrm{Cat}(\mathbf{s}_n | \boldsymbol{\pi})
+ \ln \mathrm{Dir}(\boldsymbol{\pi} | \boldsymbol{\alpha})
+ \mathrm{const.}
\\
&= \sum_{n=1}^N
\ln \prod_{k=1}^K \pi_k^{s_{n,k}}
+ \ln C_D(\boldsymbol{\alpha})
\prod_{k=1}^K \pi_k^{\alpha_k-1}
+ \mathrm{const.}
\\
&= \sum_{n=1}^N \sum_{k=1}^K
s_{n,k} \ln \pi_k
+ \ln C_D(\boldsymbol{\alpha})
+ \sum_{k=1}^K
(\alpha_k - 1) \ln \pi_k
+ \mathrm{const.}
\\
&= \sum_{k=1}^K
\left(\sum_{n=1}^N s_{n,k} + \alpha_k - 1\right)
\ln \pi_k
+ \mathrm{const.}
\tag{4.43}
\end{align}
$$
となる。$p(\boldsymbol{\lambda} | \mathbf{X}, \mathbf{S})$は左辺から移項したものである。
式(4.43)について
$$
\hat{\alpha}_k
= \sum_{n=1}^N s_{n,k} + \alpha_k
\tag{4.45}
$$
とおき
$$
\ln p(\boldsymbol{\pi} | \mathbf{X}, \mathbf{S})
= \sum_{k=1}^K
(\hat{\alpha}_k - 1)
\ln \pi_k
+ \mathrm{const.}
$$
さらに$\ln$を外し、$\mathrm{const.}$を正規化項に置き換える(正規化する)と
$$
p(\boldsymbol{\pi} | \mathbf{X}, \mathbf{S})
= C_D(\boldsymbol{\alpha})
\prod_{k=1}^K
\pi_k^{\hat{\alpha}_k-1}
= \mathrm{Dir}(\boldsymbol{\pi} | \hat{\boldsymbol{\alpha}})
$$
$\boldsymbol{\pi}$の条件付き分布は、パラメータ$\hat{\boldsymbol{\alpha}} = (\hat{\alpha}_1, \hat{\alpha}_2, \cdots, \hat{\alpha}_K)$を持つディリクレ分布になることが分かる。
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
白トピ本3章の実装を進めていました(ほぼほぼ完成しました(でも3月中にやりきるはずだった))が、ギブスサンプリングの理解不足を感じて色々飛ばしてこれをやることにしました。
世界中で制限された日々が続いていますが、そんな中で立ち上がったOsaka.Rのオンライン朝モク会に参加しています。この記事の半分はそのもくもく会中に書きました。
普段ならほげっとしてるないし寝ている時間に勉強すると、1日を長く感じてとても良いです。それだけでなく、朝に頭を使うと他の時間の頭の働きも良くなったのも嬉しいです。「やる気を出すにはまずやること」の最初の1歩のきっかけをもらえるのは、(意思の弱い)自分にとってとてもありがたいことです。折角なのでこのまま習慣にしたいなぁ。
というわけで意思の弱い人もそうでない人も、大阪関係なく参加OKとのことなので是非ご一緒しましょう！
osaka-r.connpass.com
【次節の内容】
www.anarchive-beta.com
- 2021/04/27：加筆修正しました。また、加筆修正の際にRで実装編と記事を分割しました。
朝モクも1年間続いてます！
