(自分用メモ)
観測データ
\begin{eqnarray}
W & = & (w_1, w_2, ..., w_D) \\
w_d & = & (w_{d1}, w_{d2}, ..., w_{d N_d}) \\
N_{dv} & = & \left| \{ i \mid i \in \{ 1, 2, ..., N_d \}, w_{di} = v \} \right| \quad (v \in \{ 1, 2, ..., V \})
\end{eqnarray}
$W$ が文書集合、$w_d$ が文書、$w_{dn}$ が単語 (のインデックス)、$D$ が文書数、$N_d$ が $w_d$ 内の単語数、$V$ が全単語数。
モデル
まず、文書 $w_d$ にはそれぞれ隠れた「トピック」 $z_d$ があると考える。そして、トピックはサイコロを振って決められて、文書の単語はそれの属するトピック専用のサイコロを振って決められるという世界観。つまり、
\begin{eqnarray}
z_d & \sim & {\rm Categorical}(\cdot | \theta) \\
w_{dn} & \sim & {\rm Categorical}(\cdot | \phi_{z_d})
\end{eqnarray}
ただし、$K$ をトピック数として、
\begin{eqnarray}
\theta & = & (\theta_1, \theta_2, ..., \theta_K) \\
\Phi & = & (\phi_1, \phi_2, ..., \phi_K) \\
\phi_k & = & (\phi_{k1}, \phi_{k2}, ..., \phi_{kV})
\end{eqnarray}
である。というわけで、
\begin{eqnarray}
p(W | \theta, \Phi) & = & \prod_{d=1}^D p(w_d | \theta, \phi_d) \\
p(w_d | \theta, \phi_d) & = & \sum_{k=1}^K p(w_d | \phi_k) p(z_d = k | \theta)
= \sum_k^K \theta_k \prod_{v=1}^V \phi_{kv}^{N_{dv}}
\end{eqnarray}
である。
最尤推定
$\arg \max_{\theta, \Phi} \log p(W | \theta, \Phi)$ を求めたいという想い。でもまともにやろうとすると、
\begin{eqnarray}
\log p(W | \theta, \Phi) & = & \sum_d \log \left( \sum_k \theta_k \prod_v \phi_{kv}^{N_{dv}} \right)
\end{eqnarray}
となって、せっかく $\log$ かけたのに中に和があって辛い感じになる。まあ最急降下法使えるけど、EM アルゴリズムってのがたぶんきっともっと収束が速くて汎用的な手法だからよく使われている (あと積が浮動小数だと辛いみたいな話もあるのかも)。
EM アルゴリズム
突然だけど、潜在変数であるトピック $Z = (z_1, z_2, ..., z_D)$ の確率分布 $q(Z)$ を考えて、こんなふうに変形してみる。
\begin{eqnarray}
\log p(W | \theta, \Phi) & = & \sum_Z q(Z) \log \frac{p(W, Z | \theta, \Phi)}{q(Z)} - \sum_Z q(Z) \log \frac{p(Z | W, \theta, \Phi)}{q(Z)}
\end{eqnarray}
右辺第 2 項は「パラメータを固定した」$Z$ の事後分布と $q$ のカルバック・ライブラー情報量で非負になるので、これを $0$ にするいい感じの $q$ を見つける (E ステップ)。第 2 項が $0$ になって、もともとの対数尤度の最大化は第 1 項の最大化で実現できるような雰囲気になるので、これを最大化するパラメータを見つける (M ステップ)。M ステップで第 2 項が復活するので、また E ステップに戻る。このステップを繰り返して収束したパラメータは、多くの実用的なケースではもともとの対数尤度を最大化するものになる (らしい; [Wu, 1983])。
事後分布は
\begin{eqnarray}
p(Z | W, \theta, \Phi) & = & \prod_{d=1}^D \sum_{k=1}^K \delta_{z_d k} p(z_d = k | w_d, \theta, \phi_d)
\end{eqnarray}
と展開できるので、$q(Z)$ は $D K$ 個の変数 $q_{dk} \equiv p(z_d = k | w_d, \theta, \Phi)$ で表せる。というわけで、E ステップでは以下を求める。
\begin{eqnarray}
q_{dk} & = & \frac{\theta_k \prod_v \phi_{kv}^{N_{dv}}}{\sum_{k'} \theta_{k'} \prod_v \phi_{k'v}^{N_{dv}}}
\end{eqnarray}
M ステップでは以下を最大化する。
\begin{eqnarray}
\sum_Z q(Z) \log \frac{p(W, Z | \theta, \Phi)}{q(Z)}
& = & \sum_d \sum_{z_d} q(z_d) ( \log p(w_d | \phi_{z_d}) + \log p(z_d | \theta) ) + {\rm const.} \\
& = & \sum_d \sum_k q_{dk} \left( \sum_v N_{dv} \log \phi_{kv} + \log \theta_k \right) + {\rm const.}
\end{eqnarray}
ただし、$\sum_k \theta_k = 1$, $\sum_v \phi_{dv} = 1$ という条件があるので、ラグランジュ未定乗数法を使う:
\begin{eqnarray}
\sum_d \sum_k q_{dk} \left( \sum_v N_{dv} \log \phi_{kv} + \log \theta_k \right) + \lambda_0 \left( \sum_k \theta_k - 1 \right) + \sum_d \lambda_d \left( \sum_v \phi_{dv} - 1 \right)
\end{eqnarray}
あとは普通に計算すると、
\begin{eqnarray}
\theta_k & = & \frac{\sum_d q_{dk}}{D} \\
\phi_{kv} & = & \frac{\sum_d q_{dk} N_{dv}}{\sum_d q_{dk} N_d}
\end{eqnarray}
となる。