はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。基本的な内容は「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」になります。「数式」と「プログラム」から理解するのが目標です。
この記事は、4.3.4項の内容です。「観測モデルをポアソン混合モデル」、「事前分布をガンマ分布」とする混合モデルを崩壊型ギブスサンプリング(周辺化ギブスサンプリング)を用いて推論します。
省略してある内容等ありますので、本とあわせて読んでください。初学者な自分が理解できるレベルまで落として書き下していますので、分かる人にはかなりくどくなっています。同じような立場の人のお役に立てれば幸いです。
【実装編】
www.anarchive-beta.com
www.anarchive-beta.com
【他の節一覧】
www.anarchive-beta.com
【この節の内容】
4.3.4 崩壊型ギブスサンプリング
崩壊型ギブスサンプリングを用いて、ポアソン混合モデルの事後分布$p(\mathbf{S} | \mathbf{X})$から潜在変数$\mathbf{S}$をサンプルするアルゴリズムを導出する。
観測データ$\mathbf{X}$、潜在変数$\mathbf{S}$、観測モデルのパラメータ$\boldsymbol{\lambda}$、混合比率パラメータ$\boldsymbol{\pi}$の同時分布$p(\mathbf{X}, \mathbf{S}, \boldsymbol{\lambda}, \boldsymbol{\pi})$からパラメータを周辺化(積分消去)した
$$
p(\mathbf{X}, \mathbf{S})
= \iint
p(\mathbf{X}, \mathbf{S}, \boldsymbol{\lambda}, \boldsymbol{\pi})
d\boldsymbol{\lambda} d\boldsymbol{\pi}
\tag{4.63}
$$
を用いる。
・潜在変数の条件付き分布の導出
$\mathbf{s}_n$をサンプルするための条件付き分布($\mathbf{s}_n$の事後分布)$p(\mathbf{s}_n | \mathbf{X}, \mathbf{S}_{\backslash n})$を求めていく。
観測データ$\mathbf{X}$に加えて、$\mathbf{s}_n$以外の潜在変数$\mathbf{S}_{\backslash n} = \{\mathbf{s}_1, \cdots, \mathbf{s}_{n-1}, \mathbf{s}_{n+1}, \cdots, \mathbf{s}_N\}$が観測された(サンプルされた)としたとき、$\mathbf{s}_n$の条件付き分布は
$$
\begin{align}
p(\mathbf{s}_n | \mathbf{X}, \mathbf{S}_{\backslash n})
&= \frac{
p(x_n, \mathbf{X}_{\backslash n}, \mathbf{s}_n, \mathbf{S}_{\backslash n})
}{
p(x_n, \mathbf{X}_{\backslash n}, \mathbf{S}_{\backslash n})
}
\\
&\propto
p(x_n, \mathbf{X}_{\backslash n}, \mathbf{s}_n, \mathbf{S}_{\backslash n})
\tag{4.64}\\
&= p(x_n | \mathbf{X}_{\backslash n}, \mathbf{s}_n, \mathbf{S}_{\backslash n})
p(\mathbf{X}_{\backslash n} | \mathbf{s}_n, \mathbf{S}_{\backslash n})
p(\mathbf{s}_n | \mathbf{S}_{\backslash n})
p(\mathbf{S}_{\backslash n})
\tag{4.65}\\
&= p(x_n | \mathbf{X}_{\backslash n}, \mathbf{s}_n, \mathbf{S}_{\backslash n})
p(\mathbf{X}_{\backslash n} | \mathbf{S}_{\backslash n})
p(\mathbf{s}_n | \mathbf{S}_{\backslash n})
p(\mathbf{S}_{\backslash n})
\\
&\propto
p(x_n | \mathbf{X}_{\backslash n}, \mathbf{s}_n, \mathbf{S}_{\backslash n})
p(\mathbf{s}_n | \mathbf{S}_{\backslash n})
\tag{4.66}
\end{align}
$$
で求められる。適宜$\mathbf{s}_n$に影響しない項を省いて比例関係に注目する。省略した部分については、最後に正規化することで対応できる。また、$\mathbf{X} = \{x_n, \mathbf{X}_{\backslash n}\}$、$\mathbf{S} = \{\mathbf{s}_n, \mathbf{S}_{\backslash n}\}$である。
$\mathbf{s}_n$の条件付き分布の具体的な形状を明らかにしていく。後の項は、$\mathbf{S}_{\backslash n}$が与えられた下での、パラメータ$\boldsymbol{\pi}$を周辺化した$\mathbf{s}_n$の予測分布
$$
\begin{align}
p(\mathbf{s}_n | \mathbf{S}_{\backslash n})
&= \int
p(\mathbf{s}_n, \boldsymbol{\pi} | \mathbf{S}_{\backslash n})
d\boldsymbol{\pi}
\\
&= \int
p(\mathbf{s}_n | \boldsymbol{\pi})
p(\boldsymbol{\pi} | \mathbf{S}_{\backslash n})
d\boldsymbol{\pi}
\tag{4.70}
\end{align}
$$
と解釈できる。
さらに、この式の前の項は、混合モデルの定義(4.1.2項)より
$$
p(\mathbf{s}_n | \boldsymbol{\pi})
= \mathrm{Cat}(\mathbf{s}_n | \boldsymbol{\pi})
\tag{4.2}
$$
である。
また、後ろの項は、$\mathbf{S}_{\backslash n}$が与えられた下での、$\boldsymbol{\pi}$の条件付き分布(事後分布)
$$
\begin{align}
p(\boldsymbol{\pi} | \mathbf{S}_{\backslash n})
&= \frac{
p(\mathbf{S}_{\backslash n}, \boldsymbol{\pi})
}{
p(\mathbf{S}_{\backslash n})
}
\\
&\propto
p(\mathbf{S}_{\backslash n}, \boldsymbol{\pi})
\\
&= p(\mathbf{S}_{\backslash n} | \boldsymbol{\pi})
p(\boldsymbol{\pi})
\tag{4.71}\\
&= \left\{
\prod_{n' \neq n}
p(\mathbf{s}_{n'} | \boldsymbol{\pi})
\right\}
p(\boldsymbol{\pi})
\\
&= \left\{
\prod_{n' \neq n}
\mathrm{Cat}(\mathbf{s}_{n'} | \boldsymbol{\pi})
\right\}
\mathrm{Dir}(\boldsymbol{\pi} | \boldsymbol{\alpha})
\end{align}
$$
と解釈できる。
この式は、「3.2.2：カテゴリ分布の学習と予測」の事後分布(3.25)と同じ形なので、同様の手順で
$$
p(\boldsymbol{\pi} | \mathbf{S}_{\backslash n})
= C_D(\hat{\boldsymbol{\alpha}}_{\backslash n})
\prod_{k=1}^K \pi_k^{\alpha_{k \backslash n}-1}
= \mathrm{Dir}(\boldsymbol{\pi} | \hat{\boldsymbol{\alpha}}_{\backslash n})
\tag{4.72}
$$
パラメータ$\hat{\boldsymbol{\alpha}}_{\backslash n} = (\hat{\alpha}_{1 \backslash n}, \hat{\alpha}_{2 \backslash n}, \cdots, \hat{\alpha}_{K \backslash n})$を持つディリクレ分布と求められる。ここで
$$
\hat{\alpha}_{k \backslash n}
= \sum_{n'\neq n} s_{n',k} + \alpha_k
\tag{4.73}
$$
である。
つまり、式(4.70)は
$$
\begin{align}
p(\mathbf{s}_n | \mathbf{S}_{\backslash n})
&= \int
p(\mathbf{s}_n | \boldsymbol{\pi})
p(\boldsymbol{\pi} | \mathbf{S}_{\backslash n})
d\boldsymbol{\pi}
\tag{4.70}\\
&= \int
\mathrm{Cat}(\mathbf{s}_n | \boldsymbol{\pi})
\mathrm{Dir}(\boldsymbol{\pi} | \hat{\boldsymbol{\alpha}}_{\backslash n})
d\boldsymbol{\pi}
\end{align}
$$
である。この式は、「3.2.2：カテゴリ分布の学習と予測」の予測分布(3.32')と同じ形なので、同様の手順で
$$
p(\mathbf{s}_n | \mathbf{S}_{\backslash n})
= \prod_{k=1}^K
\eta_{k \backslash n}^{s_{n,k}}
= \mathrm{Cat}(\mathbf{s}_n | \boldsymbol{\eta}_{\backslash n})
\tag{4.74}
$$
パラメータ$\boldsymbol{\eta}_{\backslash n} = (\eta_{1 \backslash n}, \eta_{2 \backslash n}, \cdots, \eta_{K \backslash n})$を持つカテゴリ分布と求められる。ここで
$$
\eta_{k \backslash n}
= \frac{
\hat{\alpha}_{k \backslash n}
}{
\sum_{k'=1}^K \hat{\alpha}_{k' \backslash n}
}
\propto
\hat{\alpha}_{k \backslash n}
\tag{4.75}
$$
である。
話戻って、式(4.66)の前の項は、$\mathbf{X}_{\backslash n},\ \mathbf{S}$が与えられた下での、パラメータ$\boldsymbol{\lambda}$を周辺化した$x_n$の予測分布
$$
\begin{align}
p(x_n | \mathbf{X}_{\backslash n}, \mathbf{s}_n, \mathbf{S}_{\backslash n})
&= \int
p(x_n, \boldsymbol{\lambda} | \mathbf{X}_{\backslash n}, \mathbf{s}_n, \mathbf{S}_{\backslash n})
d\boldsymbol{\lambda}
\\
&= \int
p(x_n | \mathbf{s}_n, \boldsymbol{\lambda})
p(\boldsymbol{\lambda} | \mathbf{X}_{\backslash n}, \mathbf{S}_{\backslash n})
d\boldsymbol{\lambda}
\tag{4.76}
\end{align}
$$
と解釈できる。
さらに、この式の前の項は、ポアソン混合モデルの定義(4.3.1項)より
$$
p(x_n | \mathbf{s}_n, \boldsymbol{\lambda})
= \prod_{k=1}^K
\mathrm{Poi}(x_n | \lambda_k)^{s_{n,k}}
\tag{4.28}
$$
である。
後の項は
$$
\begin{align}
p(\boldsymbol{\lambda} | \mathbf{X}_{\backslash n}, \mathbf{S}_{\backslash n})
&= \frac{
p(\mathbf{X}_{\backslash n}, \mathbf{S}_{\backslash n}, \boldsymbol{\lambda})
}{
p(\mathbf{X}_{\backslash n}, \mathbf{S}_{\backslash n})
}
\\
&\propto
p(\mathbf{X}_{\backslash n}, \mathbf{S}_{\backslash n}, \boldsymbol{\lambda})
\\
&= p(\mathbf{X}_{\backslash n} | \mathbf{S}_{\backslash n}, \boldsymbol{\lambda})
p(\mathbf{S}_{\backslash n})
p(\boldsymbol{\lambda})
\\
&\propto
p(\mathbf{X}_{\backslash n} | \mathbf{S}_{\backslash n}, \boldsymbol{\lambda})
p(\boldsymbol{\lambda})
\tag{4.77}\\
&= \left\{
\prod_{n'\neq n}
p(x_{n'} | \mathbf{s}_{n'}, \boldsymbol{\lambda})
\right\}
\prod_{k=1}^K
p(\boldsymbol{\lambda}_k)
\\
&= \left\{
\prod_{n'\neq n} \prod_{k=1}^K
\mathrm{Poi}(x_{n'} | \boldsymbol{\lambda})^{s_{n',k}}
\right\}
\prod_{k=1}^K
\mathrm{Gam}(\boldsymbol{\lambda}_k | a, b)
\end{align}
$$
で求められる。適宜$\boldsymbol{\lambda}$に影響しない項を省いている。
対数をとって指数部分の計算を分かりやすくして、$\boldsymbol{\lambda}$に関して整理すると
$$
\begin{align}
\ln p(\boldsymbol{\lambda} | \mathbf{X}_{\backslash n}, \mathbf{S}_{\backslash n})
&= \ln p(\mathbf{X}_{\backslash n} | \mathbf{S}_{\backslash n}, \boldsymbol{\lambda})
+ \ln p(\boldsymbol{\lambda})
+ \mathrm{const.}
\\
&= \sum_{n'\neq 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'\neq 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'\neq 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'\neq n} s_{n',k} x_{n'} + a - 1
\right)
\ln \lambda_k
- \left(
\sum_{n'\neq n} s_{n',k} + b
\right)
\lambda_k
\right\}
+ \mathrm{const.}
\tag{4.78}
\end{align}
$$
となる。適宜$\boldsymbol{\lambda}$に影響しない項を$\mathrm{const.}$にまとめている。
式(4.78)について
$$
\begin{aligned}
\hat{a}_{k \backslash n}
&= \sum_{n'\neq n} s_{n',k} x_{n'} + a
\\
\hat{b}_{k \backslash n}
&= \sum_{n'\neq n} s_{n',k} + b
\end{aligned}
\tag{4.80}
$$
とおき
$$
\ln p(\boldsymbol{\lambda} | \mathbf{X}_{\backslash n}, \mathbf{S}_{\backslash n})
= \sum_{k=1}^K \left\{
(\hat{a}_{k \backslash n} - 1)
\ln \lambda_k
- \hat{b}_{k \backslash n}
\lambda_k
\right\}
+ \mathrm{const.}
$$
さらに$\ln$を外し、$\mathrm{const.}$を正規化項に置き換える(正規化する)と
$$
p(\boldsymbol{\lambda} | \mathbf{X}_{\backslash n}, \mathbf{S}_{\backslash n})
= \prod_{k=1}^K
C_G(\hat{a}_{k \backslash n}, \hat{b}_{k \backslash n})
\lambda_k^{\hat{a}_{k \backslash n}-1}
\exp(
- \hat{b}_{k \backslash n} \lambda_k
)
= \prod_{k=1}^K
\mathrm{Gam}(\lambda_k | \hat{a}_{k \backslash n}, \hat{b}_{k \backslash n})
\tag{4.79}
$$
$\lambda_k$の条件付き分布は、パラメータ$\hat{a}_{k \backslash n},\ \hat{b}_{k \backslash n}$を持つガンマ分布になることが分かる。
つまり、式(4.28)と式(4.79)を式(4.76)に代入して、計算を簡単にするため$s_{n,k} = 1$(それ以外の$s_{n,1}, \cdots, s_{n,k-1}, s_{n,k+1}, \cdots, s_{n,K}$は0)の場合を考えると
$$
\begin{align}
p(x_n | \mathbf{X}_{\backslash n}, s_{n,k} = 1, \mathbf{S}_{\backslash n})
&= \int
p(x_n | s_{n,k} = 1, \lambda_k)
p(\lambda_k | \mathbf{X}_{\backslash n}, \mathbf{S}_{\backslash n})
d\lambda_k
\tag{4.76'}\\
&= \int
\mathrm{Poi}(x_n | \lambda_k)
\mathrm{Gam}(\lambda_k | \hat{a}_{k \backslash n}, \hat{b}_{k \backslash n})
d
\end{align}
$$
となる。
この式は、「3.2.3：ポアソン分布の学習と予測」の予測分布(3.39)と同じ形なので、同様の手順で
$$
p(x_n | \mathbf{X}_{\backslash n}, s_{n,k} = 1, \mathbf{S}_{\backslash n})
= \frac{\Gamma(x_n + \hat{r}_{k \backslash n})}{x_n! \Gamma(\hat{r}_{k \backslash n})}
(1 - \hat{p}_{k \backslash n})^{\hat{r}_{k \backslash n}}
\hat{p}_{k \backslash n}^{x_n}
= \mathrm{NB}(x_n | \hat{r}_{k \backslash n}, \hat{p}_{k \backslash n})
\tag{4.81}
$$
パラメータ$\hat{r}_{k \backslash n},\ \hat{p}_{k \backslash n}$を持つ負の二項分布と求められる。ここで
$$
\begin{aligned}
\hat{r}_{k \backslash n}
&= \hat{a}_{k \backslash n}
\\
\hat{p}_{k \backslash n}
&= \frac{1}{\hat{b}_{k \backslash n} + 1}
\end{aligned}
$$
である。
よって、$x^0 = 1$であることを利用して、$k = 1, \cdots, K$の項をまとめると
$$
p(x_n | \mathbf{X}_{\backslash n}, \mathbf{s}_n, \mathbf{S}_{\backslash n})
= \prod_{k=1}^K
\mathrm{NB}(x_n | \hat{r}_{k \backslash n}, \hat{p}_{k \backslash n})^{s_{n,k}}
\tag{4.81'}
$$
となる。
したがって、式(4.81')と式(4.74)を式(4.66)に代入すると、$\mathbf{s}_n$の条件付き分布
$$
\begin{align}
p(\mathbf{s}_n | \mathbf{X}, \mathbf{S}_{\backslash n})
&\propto
p(x_n | \mathbf{X}_{\backslash n}, \mathbf{s}_n, \mathbf{S}_{\backslash n})
p(\mathbf{s}_n | \mathbf{S}_{\backslash n})
\tag{4.66}\\
&= \left\{
\prod_{k=1}^K
\mathrm{NB}(x_n | \hat{r}_{k \backslash n}, \hat{p}_{k \backslash n})^{s_{n,k}}
\right\}
\mathrm{Cat}(\mathbf{s}_n | \boldsymbol{\eta}_{\backslash n})
\end{align}
$$
が得られる。
実際に$\mathbf{s}_n$をサンプリングするには、負の二項分布(4.81)とカテゴリ分布(4.74)の積が$s_{n,k} = 1$となる確率の比を表すので、$s_{n,1}$から$s_{n,K}$までそれぞれが1となる確率の比を求めて、$K$個の値の和が1となるように正規化したものが$\mathbf{s}_n$のサンプリング確率となる。
別の表現をすると、式(4.66)は$\prod_{k=1}^K \{ \mathrm{NB}(x_n | \hat{r}_{k \backslash n}, \hat{p}_{k \backslash n}) \eta_{k \backslash n} \}^{s_{n,k}}$となるので、波括弧を新たなパラメータとおき正規化することで$\mathbf{s}_n$が従うカテゴリ分布となる。
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
…ん？式(4.66)で$\boldsymbol{s}_n$をサンプリングするんだよな？？式(4.75)と式(4.81)を組み合わせるとどうなるんだ…？という訳で実装編はまたいずれ。
2020.04.26：元こぶしファクトリーはまちゃん二十歳の誕生日！
またかっこよく歌ってる姿が観たい。
【次節の内容】
www.anarchive-beta.com
上の件、これで理解できたと思う。混合分布を$s_{n,k}$乗で設計する妙味を1年越しで掴めた気がする。要はカテゴリ分布との相性が良いんだね！
あと、一年越しで実装できた！
