階層ディリクレ過程を実装してみる (4) k の全条件付き分布と HDP-LDA の更新式を導出

来週、サイボウズ・ラボユース合宿なんてのがあったりする(昨年の様子ちら見)。
ちょっと缶詰っぽい感じでコードをガジガジ書く機会になるので、この前の ICML 読み会で紹介した DP-MRM でも実装してみるかー、でも HDP-LDA の実装の細かいところ忘れちゃった（てへ）、TokyoML で発表する約束もしてるし、ちょっと計算やり直してみよっかな、という感じで Teh さんの HDP の論文をまたぞろひっぱりだして読んでたりする。
しかし、初めて読んだときは全くわからなくて、計算も何度も間違えて、なんかもうひどく苦労した印象しかなかったのだが、今読みかえすと、あれれこんなに簡単だったっけ……？　ディリクレ分布のお化けみたいな計算を10回くらいこなしたおかげか、なんか迷わずさくさく計算できちゃってあの苦労はなんだったんだ状態。人間の適応力って大したもんですね。

で、計算をまとめた自分のブログを見返したら、t のサンプリングまで導出しておいて、「 k のサンプリングはまた次回」で終わっている。これはひどい。

というわけで、１年前から放置されていた記事の続き「 k のサンプリング式 (34) の導出」をようやく書こう。

$p(k_{jt}=k|\boldsymbol{t}, \boldsymbol{k}^{-jt}) \propto \begin{cases}m_{\cdot k}^{-jt}f_k^{-\boldsymbol{x}_{jt}}(\boldsymbol{x}_{jt})\hspace{10px}&\text{if } k \text{ previously used,}\\\gamma f_{k^{\text{new}}}^{-\boldsymbol{x}_{jt}}(\boldsymbol{x}_{jt})&\text{if }k=k^{\text{new}}\end{cases}$ 　　(34)

例によって省略されている x を左辺の given に補って、k_jt に絡んでいる {x_ji} たちを軸に展開していくことになる。
そこで $\boldsymbol{x}_{jt}=\{x_{ji}|t_{ji}=t\}$ , $-\boldsymbol{x}_{jt}=\{x_{ji}|t_{ji}\neq t\}$ と書くことにすると、

$p(k_{jt}=k|\boldsymbol{t}, \boldsymbol{k}^{-jt},\boldsymbol{x})=\frac{p(\boldsymbol{x}_{jt}|k_{jt}=k,\boldsymbol{t},\boldsymbol{k}^{-jt},-\boldsymbol{x}_{jt})p(k_{jt}=k|\boldsymbol{t},\boldsymbol{k}^{-jt},-\boldsymbol{x}_{jt})}{p(\boldsymbol{x}_{jt}|\boldsymbol{t},\boldsymbol{k}^{-jt},-\boldsymbol{x}_{jt})}$

分子第1項は、t のサンプリングのところで使った $f_k^{-x_{ji}}(x_{ji})=p(x_{ji}|\boldsymbol{t},\boldsymbol{k},\boldsymbol{x}^{-ji})$ の x_ji がちょうど $\boldsymbol{x}_{jt}$ に置き換わったものなので、 $f_{k}^{-\boldsymbol{x}_{jt}}(\boldsymbol{x}_{jt})$ と書くことにする。つまり t_ji のサンプリングはそれと結びつく x_ji 一つで済むが、k_jt のサンプリングはテーブル t に座っているすべての x_ji たちが関わってきてしまう、ということ。
実装にあたってはこの $f_{k}^{-\boldsymbol{x}_{jt}}(\boldsymbol{x}_{jt})$ の評価が一番めんどくさい。

分子第2項は CRF なり Polya なりの式に出ている通りである。

$p(k_{jt}=k|\boldsymbol{t},\boldsymbol{k}^{-jt},-\boldsymbol{x}_{jt}) \propto \begin{cases}m_{\cdot k}^{-jt}\hspace{10px}&\text{if } k \text{ previously used,}\\\gamma&\text{if }k=k^{\text{new}}\end{cases}$

分母は正規化項に押し込んでしまえば、(34) 式が得られる。

ここまで base measure H の密度関数を h(φ)、 emission F(φ) の密度関数を f として計算してきた。
つまり HDP は任意の分布 H, F(φ) に対してモデリング可能な一般的な枠組みだったわけだ（実際には $f_k^{-x_{ji}}(x_{ji})$ の積分が計算可能であるために H は F(φ) の共役事前分布であることが望ましい)。
この HDP で特に base measure H を単語分布 H=Dir(β), F(φ) をトピック-単語間の多項分布 F(φ)=Multi(φ) と置いた場合が HDP-LDA である。
というわけでいよいよ HDP-LDA の実装に必要な更新式を出してみる。

$f_{k}^{-x_{ji}}(x_{ji})$ にそれぞれの密度関数を放り込もう。

$f(x_{ji}=v|\phi_k)=\phi_{kv}$
$h(\phi)=\frac 1Z \prod_v \phi_v^{\beta-1}$ , where $Z=\frac{\prod_v \Gamma(\beta)}{\Gamma(\sum_v \beta)}$

そこで v=x_ji, $n_{kw}^{-ji}$ を「トピック k, 単語 w をもつ x_mn たちの個数(ただし x_ji を除く)」とおくと、

$\displaystyle\begin{eqnarray} f_k^{-x_{ji}}(x_{ji}) &=& \frac{\int f(x_{ji}|\phi_k)\prod_{j'i'\neq ji,z_{j'i'}=k}f(x_{j'i'}|\phi_k)h(\phi_k)d\phi_k}{\int \prod_{j'i'\neq ji,z_{j'i'}=k}f(x_{j'i'}|\phi_k)h(\phi_k)d\phi_k} \\ &=&\frac{\int\phi_{kv}\prod_{w}\phi_{kw}^{n_{kw}^{-ji}}\prod_{w}\phi_{kw}^{\beta-1}d\phi_{k}}{\int\prod_{w}\phi_{kw}^{n_{kw}^{-ji}}\prod_{w}\phi_{kw}^{\beta-1}d\phi_{k}} \\ &=&\frac{\int\prod_{w}\phi_{kw}^{I(w=v)+n_{kw}^{-ji}+\beta-1}d\phi_{k}}{\int\prod_{w}\phi_{kw}^{n_{kw}^{-ji}+\beta-1}d\phi_{k}} \end{eqnarray}$

ここで多項ベータ積分 $\int\prod_i t_i^{\alpha_i-1}d\boldsymbol{t}=\frac{\prod_i\Gamma(\alpha_i)}{\Gamma(\sum_i\alpha_i)}$ (ただし t は t_i>0, Σt_i<1 なる単体上を動く) を思い出すと、

$\displaystyle=\frac{\Gamma(\beta+n_{kv}^{-ji}+1)\prod_{w\neq v}\Gamma(\beta+n_{kw}^{-ji})}{\Gamma(V\beta+n_{k\cdot}^{-ji}+1)}\cdot\frac{\Gamma(V\beta+n_{k\cdot}^{-ji})}{\prod_{w}\Gamma(\beta+n_{kw}^{-ji})}$

w=v 以外はキャンセルして、

$\displaystyle=\frac{\Gamma(\beta+n_{kv}^{-ji}+1)\Gamma(V\beta+n_{k\cdot}^{-ji})}{\Gamma(\beta+n_{kv}^{-ji})\Gamma(V\beta+n_{k\cdot}^{-ji}+1)}$

よってガンマ関数の特徴 $\Gamma(n+1)=n\Gamma(n)$ より

$\displaystyle f_k^{-x_{ji}}(x_{ji}) = \frac{\beta+n_{kv}^{-ji}}{V\beta+n_{k\cdot}^{-ji}}$ (ただし v=x_ji)

となる。あら簡単。
同様に $f_{k}^{-\boldsymbol{x}_{jt}}(\boldsymbol{x}_{jt})$ についても計算することができる。

$n_{kw}=\#\{x_{mn}|k_{jt_{mn}}=k, x_{mn}=w \}$ 　(トピック k 単語 w である x_mn の個数)
$n_{kw}^{-jt}=\#\{x_{mn}|k_{jt_{mn}}=k, x_{mn}=w, (m,t_{mn})\neq(j,t) \}$ 　(トピック k 単語 w である x_mn で、テーブル (j,t) に座っていない個数)

とおくと、

$\displaystyle f_{k}^{-\boldsymbol{x}_{jt}}(\boldsymbol{x}_{jt}) = \frac{\prod_w\Gamma(\beta+n_{kw}^{-jt}+n_{jtw})}{\Gamma(V\beta+n_{k\cdot}^{-jt}+n_{jt})} \cdot \frac{\Gamma(V\beta+n_{k\cdot}^{-jt})}{\prod_w\Gamma(\beta+n_{kw}^{-jt})}$
$\displaystyle = \frac{\prod_w(\beta+n_{kw}^{-jt}+n_{jtw}-1)\cdots(\beta+n_{kw}^{-jt})}{(V\beta+n_{k\cdot}^{-jt}+n_{jt}-1)\cdots(V\beta+n_{k\cdot}^{-jt})}$
【2012/08/22 修正】

この２つを (32), (33), (34) 式に放り込めば、HDP-LDA の実装に必要な数式が揃う。

最後頑張ってガンマ関数を消しているが、実装で手抜きをしたいならガンマ関数を展開せず log-gamma 関数を使うのが楽かもしれない(対数を取らないと、浮動小数点の範囲をさくっとオーバーフローする)。
この HDP-LDA の実装では、ガンマ関数を展開した方の式で結果をキャッシュして使いまわしている(cur_log_base 〜 log_f_k_new_x_jt メソッドのあたり)が、Python 実装なら log-gamma 関数をばしっと呼んだほうが速かったかもしれない……。

【追記】
あ、f_{k^new} 忘れてた。

$f_{k^{\text{new}}}^{-x_{ji}}(x_{ji})=\int p(x_{ji}|\phi)p(\phi)d\boldsymbol{\phi}=\int \phi_v\frac{\Gamma(V\beta)}{\prod_w\Gamma(\beta)}\prod_w \phi_w^{\beta-1}d\boldsymbol{\phi}$
$=\frac{\Gamma(V\beta)}{\prod_w\Gamma(\beta)}\frac{\Gamma(\beta+1)\prod_{w\neq v}\Gamma(\beta)}{\Gamma(V\beta+1)}=\frac{1}{V}$

同様に

$f_{k^{\text{new}}}^{-\boldsymbol{x}_{jt}}(\boldsymbol{x}_{jt})=\frac{\Gamma(V\beta)\prod_w\Gamma(\beta+n_{jtw})}{\Gamma(V\beta+n_{jt})\Gamma(\beta)^V}$ 【2012/9/26 修正】
【／追記】

階層ディリクレ過程を実装してみる (5) HDP-LDA の更新式を導出 (パープレキシティ) - 木曜不足