はじめに

pythonでトピックモデルの変分ベイズ推定を実装しました．
教科書として『トピックモデル』を使いました．
※ 以前書いた pythonでトピックモデルの最尤推定実装も合わせて読んでみてください．

本記事の構成

はじめに
トピックモデル
準備
- ベータ関数
- ディリクレ分布
- ディリクレ分布に従う確率変数
変分ベイズ推定
- 変分下限
- トピックの変分事後分布
- トピック分布と単語分布の分解
- トピック分布の変分事後分布
- 単語分布の変分事後分布
パラメータ推定
実装
結果
おわりに

トピックモデル

トピックモデルについての説明は，pythonでトピックモデルの最尤推定実装に記載したので割愛します．

準備

必要な式を先に導出しておきます．

ベータ関数
ディリクレ分布
ディリクレ分布に従う確率変数

ベータ関数

ベータ関数は以下のように表されます．

\begin{align}
\int_0^1 \phi^{\alpha - 1} (1 - \phi)^{\beta - 1} d\phi &= \left[ \cfrac{\phi^{\alpha}}{\alpha} (1 - \phi)^{\beta - 1} \right]_0^1 + \int_0^1 \cfrac{\phi^{\alpha}}{\alpha} (\beta - 1) (1 - \phi)^{\beta - 2} d\phi \\
&= \cfrac{\beta - 1}{\alpha} \int_0^1 \phi^{\alpha} (1 - \phi)^{\beta - 2} d\phi \\
&= \cfrac{\beta - 1}{\alpha} \cdots \cfrac{1}{\alpha + \beta - 2} \int_0^1 \phi^{\alpha + \beta - 2} d\phi \\
&= \cfrac{(\beta - 1) \cdots 1}{\alpha \cdots (\alpha + \beta - 1)} \\
&= \cfrac{(\beta - 1)!}{\cfrac{(\alpha + \beta - 1)!}{(\alpha - 1)!}} \\
&= \cfrac{(\alpha - 1)!(\beta - 1)!}{(\alpha + \beta - 1)!} \\
&= \cfrac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha + \beta)} \equiv B(\alpha, \beta)
\end{align}

これより以下の式が成り立ちます．

\begin{align}
\int_0^{1 - q} x^{\alpha - 1} (1 - q - x)^{\beta - 1} dx
&= \int_0^1 \bigl( (1 - q) y \bigr)^{\alpha - 1} \bigl( (1 - q) (1 - y) \bigr)^{\beta - 1} (1 - q) dy \ \ \ \ \because x = (1 - q)y \\
&= (1 - q)^{\alpha + \beta + 1} \int_0^1 y^{\alpha - 1} (1 - y)^{\beta - 1} dy \\
&= (1 - q)^{\alpha + \beta + 1} B(\alpha, \beta) \\
\end{align}

ディリクレ分布

ディリクレ分布は以下のように表されます．

\begin{align}
{\rm Diriclet}(\boldsymbol \phi \mid \boldsymbol \beta)
= \cfrac{\displaystyle \prod_{v = 1}^V \phi_v^{\beta_v - 1}}{\displaystyle \int \prod \limits_{v = 1}^V \phi_v^{\beta_v - 1} d \boldsymbol \phi} \\
\end{align}

正規化項を展開します．

\begin{align}
\int \prod_{v = 1}^V \phi_v^{\beta_v - 1} d\boldsymbol \phi
&= \int_0^1 \phi_1^{\beta_1 - 1} \int_0^{1 - \phi_1} \phi_2^{\beta_2 - 1} \cdots \int_0^{1 - \sum\limits_{v = 1}^{V - 2} \phi_v} \phi_{V - 1}^{\beta_{V - 1} - 1} \left( 1 - \sum_{v = 1}^{V - 1}\phi_v \right)^{\beta_V - 1} d\phi_{V - 1} \cdots d\phi_2 d\phi_1 \\
&= \int_0^1 \phi_1^{\beta_1 - 1} \int_0^{1 - \phi_1} \phi_2^{\beta_2 - 1} \cdots \int_0^{1 - \sum\limits_{v = 1}^{V - 3} \phi_v} \phi_{V - 2}^{\beta_{V - 2} - 1} \left( 1 - \sum_{v = 1}^{V - 2}\phi_v \right)^{\beta_{V - 1} + \beta_V - 1} B(\beta_{V - 1}, \beta_V) d\phi_{V - 2} \cdots d\phi_2 d\phi_1 \\
&= \int_0^1 \phi_1^{\beta_1 - 1} \int_0^{1 - \phi_1} \phi_2^{\beta_2 - 1} \cdots \int_0^{1 - \sum\limits_{v = 1}^{V - 4} \phi_v} \phi_{V - 3}^{\beta_{V - 3} - 1} \left( 1 - \sum_{v = 1}^{V - 3}\phi_v \right)^{\beta_{V - 2} + \beta_{V - 1} + \beta_V - 1} B(\beta_{V - 2}, \beta_{V - 1} + \beta_V) B(\beta_{V - 1}, \beta_V) d\phi_{V - 3} \cdots d\phi_2 d\phi_1 \\
&= B \left( \beta_1, \sum_{v = 2}^V \beta_v \right) \cdots B(\beta_{V - 2}, \beta_{V - 1} + \beta_V) B(\beta_{V - 1}, \beta_V) \\
&= \cfrac{\Gamma(\beta_1) \Gamma \left( \sum\limits_{v = 2}^V \beta_v \right)}{\Gamma \left( \sum\limits_{v = 1}^V \beta_v \right)} \cdots \cfrac{\Gamma(\beta_{V - 2}) \Gamma(\beta_{V - 1} + \beta_V)}{\Gamma(\beta_{V - 2} + \beta_{V - 1} + \beta_V)} \cfrac{\Gamma(\beta_{V - 1}) \Gamma(\beta_V)}{\Gamma(\beta_{V - 1} + \beta_V)} \\
&= \cfrac{\prod\limits_{v = 1}^V \Gamma(\beta_v)}{\Gamma \left( \sum\limits_{v = 1}^V \beta_v \right)} \\
\end{align}

正規化項を置き換えます．

\begin{align}
{\rm Diriclet}(\boldsymbol \phi \mid \boldsymbol \beta)
= \cfrac{\displaystyle \Gamma \left( \sum_{v = 1}^V \beta_v \right)}{\displaystyle \prod_{v = 1}^V \Gamma(\beta_v)} \prod_{v = 1}^V \phi_v^{\beta_v - 1} \\
\end{align}

ディリクレ分布に従う確率変数

ディリクレ分布に従う確率変数 $\boldsymbol \phi$ の $v$ 番目の要素 $\phi_v$ の対数の期待値は以下のように表せます．

\begin{align}
\int p(\boldsymbol \phi \mid \boldsymbol \beta) \log \phi_v d \boldsymbol \phi
&= \int \cfrac{\Gamma(\hat \beta)}{\prod\limits_{v' = 1}^V \Gamma(\beta_{v'})} \prod_{v' = 1}^V \phi_{v'}^{\beta_{v'} - 1} \log \phi_v d \boldsymbol \phi \\
&= \cfrac{\Gamma(\hat \beta)}{\prod\limits_{v' = 1}^V \Gamma(\beta_{v'})} \int \cfrac{\partial \prod\limits_{v' = 1}^V \phi_{v'}^{\beta_{v'} - 1}}{\partial \beta_v} d \boldsymbol \phi \\
&= \cfrac{\Gamma(\hat \beta)}{\prod\limits_{v' = 1}^V \Gamma(\beta_{v'})} \cfrac{\partial}{\partial \beta_v} \int \prod_{v' = 1}^V \phi_{v'}^{\beta_{v'} - 1} d \boldsymbol \phi \\
&= \cfrac{\Gamma(\hat \beta)}{\prod\limits_{v' = 1}^V \Gamma(\beta_{v'})} \cfrac{\partial}{\partial \beta_v} \cfrac{\prod\limits_{v' = 1}^V \Gamma(\beta_{v'})}{\Gamma(\hat \beta)} \\
&= \cfrac{\partial}{\partial \beta_v} \log \cfrac{\prod\limits_{v' = 1}^V \Gamma(\beta_{v'})}{\Gamma(\hat \beta)} \\
&= \cfrac{\partial}{\partial \beta_v} \left( \log \prod_{v' = 1}^V \Gamma(\beta_{v'}) - \log \Gamma(\hat \beta) \right) \\
&= \cfrac{\partial \log \Gamma(\beta_v)}{\partial \beta_v} - \cfrac{\partial \log \Gamma(\hat \beta)}{\partial \hat \beta} = \Psi(\beta_v) - \Psi(\hat \beta) \\
\end{align}

ここで，$\hat \beta$ はパラメータの総和，$\Psi(x)$ はディガンマ関数です．

\begin{align}
& \hat \beta = \sum_{v = 1}^V \beta_v \\
& \Psi(x) = \cfrac{d}{dx} \log \Gamma(x) \\
\end{align}

変分ベイズ推定

変分ベイズ推定では，未知変数の事後分布を推定します．
トピックモデルにおける未知変数は，以下の3種類になります．

未知変数	表記
トピック集合	$\boldsymbol Z = (z_{11}, \cdots, z_{1N_1}, z_{21}, \cdots, z_{DN_D}), \ \ z_{dn} \in \{1, \cdots, K\}$
トピック分布集合	$\boldsymbol \Theta = (\boldsymbol \theta_1, \cdots, \boldsymbol \theta_D), \ \ \boldsymbol \theta_d = (\theta_{d1}, \cdots, \theta_{dK})$
単語分布集合	$\boldsymbol \Phi = (\boldsymbol \phi_1, \cdots, \boldsymbol \phi_K), \ \ \boldsymbol \phi_k = (\phi_{k1}, \cdots, \phi_{kV})$

$D$ は文書数，$K$ はトピック数，$V$ は語彙数，$N_d$ は文書 $d$ の文書長を表します．
$z_{dn}$ は文書 $d$ の $n$ 番目の単語のトピック，
$\theta_{dk}$ は文書 $d$ にトピック $k$ が割り当てられる確率，
$\phi_{kv}$ はトピック $k$ で語彙 $v$ が生成される確率を表します
また，トピック分布集合 $\Theta$ はパラメータ $\alpha$ を持つディリクレ分布から生成され，
単語分布集合 $\Phi$ はパラメータ $\beta$ を持つディリクレ分布から生成されます．

変分下限

イェンセンの不等式を用いて $\log p(\boldsymbol W \mid \alpha, \beta)$ の変分下限 $F$ を求めます．

\begin{align}
\log p(\boldsymbol W \mid \alpha, \beta)
&= \log \int \int \sum_{\boldsymbol Z} p(\boldsymbol W, \boldsymbol Z, \boldsymbol \Theta, \boldsymbol \Phi \mid \alpha, \beta) d\boldsymbol \Theta d \boldsymbol \Phi \\
&= \log \int \int \sum_{\boldsymbol Z} q(\boldsymbol Z, \boldsymbol \Theta, \boldsymbol \Phi) \cfrac{p(\boldsymbol W, \boldsymbol Z, \boldsymbol \Theta, \boldsymbol \Phi \mid \alpha, \beta)}{q(\boldsymbol Z, \boldsymbol \Theta, \boldsymbol \Phi)} d\boldsymbol \Theta d \boldsymbol \Phi \\
&\geq \int \int \sum_{\boldsymbol Z} q(\boldsymbol Z, \boldsymbol \Theta, \boldsymbol \Phi) \log \cfrac{p(\boldsymbol W, \boldsymbol Z, \boldsymbol \Theta, \boldsymbol \Phi \mid \alpha, \beta)}{q(\boldsymbol Z, \boldsymbol \Theta, \boldsymbol \Phi)} d\boldsymbol \Theta d \boldsymbol \Phi \\
&= \int \int \sum_{\boldsymbol Z} q(\boldsymbol Z) q(\boldsymbol \Theta, \boldsymbol \Phi) \log \cfrac{p(\boldsymbol W, \boldsymbol Z, \boldsymbol \Theta, \boldsymbol \Phi \mid \alpha, \beta)}{q(\boldsymbol Z) q(\boldsymbol \Theta, \boldsymbol \Phi)} d\boldsymbol \Theta d \boldsymbol \Phi \equiv F \\
\end{align}

トピックの変分事後分布

$F(q(\boldsymbol Z))$ を以下のように定義し，$q(\boldsymbol Z)$ で変分します．

\begin{align}
F(q(\boldsymbol Z))
&\equiv \int \int \sum_{\boldsymbol Z} q(\boldsymbol Z) q(\boldsymbol \Theta, \boldsymbol \Phi) \log \cfrac{p(\boldsymbol W, \boldsymbol Z, \boldsymbol \Theta, \boldsymbol \Phi \mid \alpha, \beta)}{q(\boldsymbol Z) q(\boldsymbol \Theta, \boldsymbol \Phi)} d \boldsymbol \Theta d \boldsymbol \Phi + \lambda \left( \sum_{\boldsymbol Z} q(\boldsymbol Z) - 1 \right) \\
\cfrac{\partial F(q(\boldsymbol Z))}{\partial q(\boldsymbol Z)}
&= \int \int q(\boldsymbol \Theta, \boldsymbol \Phi) \bigl( \log p(\boldsymbol W, \boldsymbol Z, \boldsymbol \Theta, \boldsymbol \Phi \mid \alpha, \beta) - \log q(\boldsymbol Z) - \log q(\boldsymbol \Theta, \boldsymbol \Phi) - 1 \bigr) d\boldsymbol \Theta d \boldsymbol \Phi + \lambda \\
&= \mathbb{E}_{q(\boldsymbol \Theta, \boldsymbol \Phi)} \bigl[ \log p(\boldsymbol W, \boldsymbol Z, \boldsymbol \Theta, \boldsymbol \Phi \mid \alpha, \beta) \bigr] - \log q(\boldsymbol Z) + \lambda - C \\
\end{align}

$\displaystyle \cfrac{\partial F(q(\boldsymbol Z))}{\partial q(\boldsymbol Z)} = 0$ を解きます．

\begin{align}
q(\boldsymbol Z)
&\propto \exp \Bigl( \mathbb{E}_{q(\boldsymbol \Theta, \boldsymbol \Phi)} \bigl[ \log p(\boldsymbol W, \boldsymbol Z, \boldsymbol \Theta, \boldsymbol \Phi \mid \alpha, \beta) \bigr] \Bigr) \\
&= \exp \Bigl( \mathbb{E}_{q(\boldsymbol \Theta, \boldsymbol \Phi)} \bigl[ \log p(\boldsymbol \Theta \mid \alpha) + \log p(\boldsymbol \Phi \mid \beta) + \log p(\boldsymbol Z \mid \boldsymbol \Theta) + \log p(\boldsymbol W \mid \boldsymbol Z, \boldsymbol \Phi) \bigr] \Bigr) \\
&\propto \exp \Bigl( \mathbb{E}_{q(\boldsymbol \Theta, \boldsymbol \Phi)} \bigl[ \log p(\boldsymbol Z \mid \boldsymbol \Theta) + \log p(\boldsymbol W \mid \boldsymbol Z, \boldsymbol \Phi) \bigr] \Bigr) \\
&= \exp \Bigl( \mathbb{E}_{q(\boldsymbol \Theta)} \bigl[ \log p(\boldsymbol Z \mid \boldsymbol \Theta) \bigr] + \mathbb{E}_{q(\boldsymbol \Phi)} \bigl[ \log p(\boldsymbol W \mid \boldsymbol Z, \boldsymbol \Phi) \bigr] \Bigr) \\
&= \exp \Bigl( \mathbb{E}_{q(\boldsymbol \Theta)} \bigl[ \log \prod_{d = 1}^D \prod_{n = 1}^{N_d} \theta_{dz_{dn}} \bigr] + \mathbb{E}_{q(\boldsymbol \Phi)} \bigl[ \log \prod_{d = 1}^{D} \prod_{n = 1}^{N_d} \phi_{z_{dn}w_{dn}} \bigr] \Bigr) \\
&= \exp \Bigl( \mathbb{E}_{q(\boldsymbol \Theta)} \bigl[ \sum_{d = 1}^D \sum_{n = 1}^{N_d} \log \theta_{dz_{dn}} \bigr] + \mathbb{E}_{q(\boldsymbol \Phi)} \bigl[ \sum_{d = 1}^D \sum_{n = 1}^{N_d} \log \phi_{z_{dn}w_{dn}} \bigr] \Bigr) \\
&= \prod_{d = 1}^D \prod_{n = 1}^{N_d} \exp \Bigl( \mathbb{E}_{q(\boldsymbol \Theta)} \bigl[ \log \theta_{dz_{dn}} \bigr] + \mathbb{E}_{q(\boldsymbol \Phi)} \bigl[ \log \phi_{z_{dn}w_{dn}} \bigr] \Bigr) \\
&= \prod_{d = 1}^D \prod_{n = 1}^{N_d} \exp \Bigl( \Psi \bigl( \alpha_{dz_{dn}} \bigr) - \Psi \bigl( \sum_{k = 1}^K \alpha_{dk} \bigr) + \Psi \bigl( \beta_{z_{dn}w_{dn}} \bigr) - \Psi \bigl( \sum_{v = 1}^V \beta_{z_{dn}v} \bigr) \Bigr) \\
\end{align}

$\displaystyle q(\boldsymbol Z) = \prod_{d = 1}^D \prod_{n = 1}^{N_d} q_{dnz_{dn}}$ より，$q_{dnk}$ は以下のようになります．

\begin{align}
q_{dnk}
&\propto \exp \Bigl( \Psi \bigl( \alpha_{dk} \bigr) - \Psi \bigl( \sum_{k' = 1}^K \alpha_{dk'} \bigr) + \Psi \bigl( \beta_{kw_{dn}} \bigr) - \Psi \bigl( \sum_{v = 1}^V \beta_{kv} \bigr) \Bigr) \\
\end{align}

トピック分布と単語分布の分解

$F(q(\boldsymbol \Theta, \boldsymbol \Phi))$ を以下のように定義し，$q(\boldsymbol \Theta, \boldsymbol \Phi)$ で変分します．

\begin{align}
F(q(\boldsymbol \Theta, \boldsymbol \Phi))
&\equiv \int \int \sum_{\boldsymbol Z} q(\boldsymbol Z) q(\boldsymbol \Theta, \boldsymbol \Phi) \log \cfrac{p(\boldsymbol W, \boldsymbol Z, \boldsymbol \Theta, \boldsymbol \Phi \mid \alpha, \beta)}{q(\boldsymbol Z) q(\boldsymbol \Theta, \boldsymbol \Phi)} d \boldsymbol \Theta d \boldsymbol \Phi + \lambda \left( \int \int q(\boldsymbol \Theta, \boldsymbol \Phi) d \boldsymbol \Theta d \boldsymbol \Phi - 1 \right) \\
\cfrac{\partial F(q(\boldsymbol \Theta, \boldsymbol \Phi))}{\partial q(\boldsymbol \Theta, \boldsymbol \Phi)}
&= \sum_{\boldsymbol Z} q(\boldsymbol Z) \bigl( \log p(\boldsymbol W, \boldsymbol Z, \boldsymbol \Theta, \boldsymbol \Phi \mid \alpha, \beta) - \log q(\boldsymbol Z) - \log q(\boldsymbol \Theta, \boldsymbol \Phi) - 1 \bigr) + \lambda \\
&= \mathbb{E}_{q(\boldsymbol Z)} \bigl[ \log p(\boldsymbol W, \boldsymbol Z, \boldsymbol \Theta, \boldsymbol \Phi \mid \alpha, \beta) \bigr] - \log q(\boldsymbol \Theta, \boldsymbol \Phi) + \lambda - C \\
\end{align}

$\displaystyle \cfrac{\partial F(q(\boldsymbol \Theta, \boldsymbol \Phi))}{\partial q(\boldsymbol \Theta, \boldsymbol \Phi)} = 0$ を解きます．

\begin{align}
q(\boldsymbol \Theta, \boldsymbol \Phi)
&\propto \exp \Bigl( \mathbb{E}_{q(\boldsymbol Z)} \bigl[ \log p(\boldsymbol W, \boldsymbol Z, \boldsymbol \Theta, \boldsymbol \Phi \mid \alpha, \beta) \bigr] \Bigr) \\
&= \exp \Bigl( \mathbb{E}_{q(\boldsymbol Z)} \bigl[ \log p(\boldsymbol \Theta \mid \alpha) + \log p(\boldsymbol \Phi \mid \beta) + \log p(\boldsymbol Z \mid \boldsymbol \Theta) + \log p(\boldsymbol W \mid \boldsymbol Z, \boldsymbol \Phi) \bigr] \Bigr) \\
&= \exp \Bigl( \mathbb{E}_{q(\boldsymbol Z)} \bigl[ \log p(\boldsymbol Z \mid \boldsymbol \Theta) \bigr] + \log p(\boldsymbol \Theta \mid \alpha) \Bigr) \times \exp \Bigl( \mathbb{E}_{q(\boldsymbol Z)} \bigl[ \log p(\boldsymbol W \mid \boldsymbol Z, \boldsymbol \Phi) \bigr] + \log p(\boldsymbol \Phi \mid \beta) \Bigr) \\
\end{align}

上式より，$q(\boldsymbol \Theta, \boldsymbol \Phi) = q(\boldsymbol \Theta) q(\boldsymbol \Phi)$ と分解でき，それぞれ以下のように表せます．

\begin{align}
& q(\boldsymbol \Theta) \propto \exp \Bigl( \mathbb{E}_{q(\boldsymbol Z)} \bigl[ \log p(\boldsymbol Z \mid \boldsymbol \Theta) \bigr] + \log p(\boldsymbol \Theta \mid \alpha) \Bigr) \\
& q(\boldsymbol \Phi) \propto \exp \Bigl( \mathbb{E}_{q(\boldsymbol Z)} \bigl[ \log p(\boldsymbol W \mid \boldsymbol Z, \boldsymbol \Phi) \bigr] + \log p(\boldsymbol \Phi \mid \beta) \Bigr) \\
\end{align}

トピック分布の変分事後分布

$q(\boldsymbol \Theta)$ は以下のように計算できます．

\begin{align}
q(\boldsymbol \Theta) &\propto \exp \Bigl( \mathbb{E}_{q(\boldsymbol Z)} \bigl[ \log p(\boldsymbol Z \mid \boldsymbol \Theta) \bigr] + \log p(\boldsymbol \Theta \mid \alpha) \Bigr) \\
&= \exp \Bigl( \mathbb{E}_{q(\boldsymbol Z)} \bigl[ \sum_{d = 1}^D \sum_{n = 1}^{N_d} \log \theta_{dz_{dn}} \bigr] + \log \prod_{d = 1}^D p(\boldsymbol \theta_d \mid \alpha) \Bigr) \\
&= \exp \Bigl( \sum_{d = 1}^D \sum_{n = 1}^{N_d} \sum_{k = 1}^{K} q_{dnk} \log \theta_{dk} + \sum_{d = 1}^D \log \cfrac{\Gamma(\alpha K)}{\Gamma(\alpha)^K} \prod_{k = 1}^K \theta_{dk}^{\alpha - 1} \Bigr) \\
&\propto \exp \Bigl( \sum_{d = 1}^D \sum_{n = 1}^{N_d} \sum_{k = 1}^K q_{dnk} \log \theta_{dk} + \sum_{d = 1}^D \sum_{k = 1}^K \log \theta_{dk}^{\alpha - 1} \Bigr) \\
&= \exp \Bigl( \sum_{d = 1}^D \sum_{k = 1}^K \log \theta_{dk}^{\sum \limits_{n = 1}^{N_d} q_{dnk}} + \sum_{d = 1}^D \sum_{k = 1}^K \log \theta_{dk}^{\alpha - 1} \Bigr) \\
&= \prod_{d = 1}^D \prod_{k = 1}^K \theta_{dk}^{\alpha + \sum \limits_{n = 1}^{N_d} q_{dnk} - 1} \\
q(\boldsymbol \Theta) &= \prod_{d = 1}^D {\rm Diriclet}(\boldsymbol \theta_d \mid \alpha_{d1}, \cdots, \alpha_{dK}) \\
\end{align}

ここで，$\alpha_{dk}$ は以下のように定義しています．

\begin{align}
\alpha_{dk} &= \alpha + \sum_{n = 1}^{N_d} q_{dnk} \\
\end{align}

単語分布の変分事後分布

$q(\boldsymbol \Phi)$ は以下のように計算できます．

\begin{align}
q(\boldsymbol \Phi) &\propto \exp \Bigl( \mathbb{E}_{q(\boldsymbol Z)} \bigl[ \log p(\boldsymbol W \mid \boldsymbol Z, \boldsymbol \Phi) \bigr] + \log p(\boldsymbol \Phi \mid \beta) \Bigr) \\
&= \exp \Bigl( \mathbb{E}_{q(\boldsymbol Z)} \bigl[ \sum_{d = 1}^D \sum_{n = 1}^{N_d} \log \phi_{z_{dn}w_{dn}} \bigr] + \log \prod_{k = 1}^K p(\boldsymbol \phi_k \mid \beta) \Bigr) \\
&= \exp \Bigl( \sum_{d = 1}^D \sum_{n = 1}^{N_d} \sum_{k = 1}^K q_{dnk} \log \phi_{kw_{dn}} + \sum_{k = 1}^K \log \cfrac{\Gamma(\beta V)}{\Gamma(\beta)^V} \prod_{v = 1}^V \phi_{kv}^{\beta - 1} \Bigr) \\
&\propto \exp \Bigl( \sum_{d = 1}^D \sum_{n = 1}^{N_d} \sum_{k = 1}^K q_{dnk} \log \phi_{kw_{dn}} + \sum_{k = 1}^K \sum_{v = 1}^V \phi_{kv}^{\beta - 1} \Bigr) \\
&= \exp \Bigl( \sum_{k = 1}^K \sum_{v = 1}^V \log \phi_{kv}^{\sum \limits_{d = 1}^D \sum \limits_{n: w_{dn}= v} q_{dnk}} + \sum_{k = 1}^K \sum_{v = 1}^V \log \phi_{kv}^{\beta - 1} \Bigr) \\
&= \prod_{k = 1}^K \prod_{v = 1}^V \phi_{kv}^{\beta + \sum \limits_{d = 1}^D \sum \limits_{n: w_{dn}= v} q_{dnk} - 1} \\
q(\boldsymbol \Phi) &= \prod_{k = 1}^K {\rm Diriclet}(\boldsymbol \phi_k \mid \beta_{k1}, \cdots, \beta_{kV}) \\
\end{align}

ここで，$\beta_{kv}$ は以下のように定義しています．

\begin{align}
\beta_{kv} &= \beta + \sum_{d = 1}^D \sum_{n: w_{dn}= v} q_{dnk} \\
\end{align}

パラメータ推定

導出した結果を用いてパラメータを推定します．

$\alpha_{dk}$ および $\beta_{kv}$ を乱数で初期化
$q_{dnk} \propto \exp \Bigl( \Psi \bigl( \alpha_{dk} \bigr) - \Psi \bigl( \sum_{k' = 1}^K \alpha_{dk'} \bigr) + \Psi \bigl( \beta_{kw_{dn}} \bigr) - \Psi \bigl( \sum_{v = 1}^V \beta_{kv} \bigr) \Bigr)$ よりトピックの変分事後分布を計算
$q_{dnk}$ を正規化
パラメータを更新 $\alpha_{dk} \leftarrow \alpha_{dk} + q_{dnk}$, $\beta_{kw_{dn}} \leftarrow \beta_{kw_{dn}} + q_{dnk}$

実装

pythonでトピックモデルの変分ベイズ推定のトイプログラムを実装しました．
理解と式展開は苦労しましたが，得られる結果が綺麗なのでコードも短く済んでいます．

import numpy as np
from scipy.special import digamma

def normalize(ndarray, axis):
    return ndarray / ndarray.sum(axis = axis, keepdims = True)

def normalized_random_array(d0, d1):
    ndarray = np.random.rand(d0, d1)
    return normalize(ndarray, axis = 1)

if __name__ == "__main__":

    # initialize parameters
    D, K, V = 1000, 2, 6
    alpha0, beta0 = 1.0, 1.0
    alpha = alpha0 + np.random.rand(D, K)
    beta = beta0 + np.random.rand(K, V)
    theta = normalized_random_array(D, K)
    phi = normalized_random_array(K, V)

    # for generate documents
    _theta = np.array([theta[:, :k+1].sum(axis = 1) for k in range(K)]).T
    _phi = np.array([phi[:, :v+1].sum(axis = 1) for v in range(V)]).T

    # generate documents
    W, Z = [], []
    N = np.random.randint(100, 300, D)
    for (d, N_d) in enumerate(N):
        Z.append((np.random.rand(N_d, 1) < _theta[d, :]).argmax(axis = 1))
        W.append((np.random.rand(N_d, 1) < _phi[Z[-1], :]).argmax(axis = 1))

    # estimate parameters
    T = 30
    for t in range(T):
        dig_alpha = digamma(alpha) - digamma(alpha.sum(axis = 1, keepdims = True))
        dig_beta = digamma(beta) - digamma(beta.sum(axis = 1, keepdims = True))

        alpha_new = np.ones((D, K)) * alpha0
        beta_new = np.ones((K, V)) * beta0
        for (d, N_d) in enumerate(N):
            # q
            q = np.zeros((V, K))
            v, count = np.unique(W[d], return_counts = True)
            q[v, :] = (np.exp(dig_alpha[d, :].reshape(-1, 1) + dig_beta[:, v]) * count).T
            q[v, :] /= q[v, :].sum(axis = 1, keepdims = True)

            # alpha, beta
            alpha_new[d, :] += count.dot(q[v])
            beta_new[:, v] += count * q[v].T

        alpha = alpha_new.copy()
        beta = beta_new.copy()

    theta_est = np.array([np.random.dirichlet(a) for a in alpha])
    phi_est = np.array([np.random.dirichlet(b) for b in beta])

結果

上記トイプログラムで求めた $\boldsymbol \beta$ から生成した $\boldsymbol \Phi$ を掲載します．
本来はパープレキシティや対数尤度を評価尺度とするべきですが，面倒なのでやってません．

phi
[[ 0.26631554  0.04657097  0.29425041  0.1746378   0.03077238  0.1874529 ]
 [ 0.2109456   0.01832505  0.30360253  0.09073456  0.14039401  0.23599826]]
phi_est
[[ 0.27591967  0.0424522   0.26088712  0.18220604  0.02874477  0.2097902 ]
 [ 0.19959096  0.02327517  0.34107528  0.08462041  0.14291539  0.20852279]]

おわりに

トピックモデルの変分ベイズ推定の導出と実装ができました．
本記事は結果があまり面白くありませんでしたが，
LDA for Pokemon analysisを見て金銀版のポケモンの分類をやってみたので，
そちらも後日Qiitaにまとめます．

pythonでトピックモデルの変分ベイズ推定実装