混合ガウス分布の変分下界の計算式

やっぱりまじめに実装して確認しないと、ということで PRML 10.2.2 の変分下界の膨大な式、つまり (10.71)+(10.72)+(10.73)+(10.74)-(10.75)-(10.76)-(10.77) を変形＆整理していったら、打ち消しあって消えて消えて、残ったのはたったこれだけ。

$\begin{eqnarray}\mathcal{L} &=& \ln C(\bf{\alpha}_0)-\ln C(\bf{\alpha}) +\frac{D}{2}\sum_{k=1}^K(\ln\beta_0-\ln\beta_k)+\sum_{k=1}^K \{\ln B(\bf{W}_0,\nu_0)-\ln B(\bf{W}_k,\nu_k)\} \\&& - \sum_{n=1}^N\sum_{k=1}^K r_{nk}\ln r_{nk} -\frac{DN}{2}\ln(2\pi) \end{eqnarray}$

\tilde{π} や \tilde{Λ} が消えただけでもすごいすごい、と思っていたが、二次形式を計算していったら W_k^{-1} が出てきて W_k と積を取って消えて、さらにそのトレースから出てきた D も H[Λ_k] の中の ν_k D/2 と打ち消しあって、消えた。
さすがにここまでうまくいくと(しかも式の形がやたらきれい)、何か間違っているんじゃあないかと不安になってきたので、晒してみる。

【追記】実装してみた。

http://github.com/shuyo/iir/blob/master/unsupervised/vb.r

ちゃんと増え続けていくので、大丈夫？
実装で困ったのは、ディリクレ分布の正規化項を計算するのに、ガンマ関数が発散してしまうこと。
悩みながら調べていたら、R はガンマ関数の対数 lgamma() を持つことを知り、これを使って解決。R は行き届いてるね。

さっそくいろいろ感想＆疑問。

変分下界やパラメータの値を表示させながら回すと、変分下界が収束した後もパラメータはまだもう少し変動する。対数尤度は変分下界より収束が遅い？
変分下界は徐々にある値に近づいていくのではなく、急停止的に値の増加が止まる。なんとなく不思議。エントロピー最大化法とかまだ試したこと無いけど、その手の奴もこういう挙動なんだろうか？
このスクリプトでの変分下界の収束先は実行するたびごとにばらばらで、変動もそこそこある。でも PRML の図10.7 では、ばらつきは量子的な感じ。やっぱりまだどこか間違ってるんだろうか……→★解決。α_0, β_0, ν_0を固定、m_k の初回のみランダムにしたら、変分下界の収束先がほぼ一定値で、縮退数が異なる場合のみ違うという非常に直感に適合した結果が出るようになった。PRML もα_0, β_0, ν_0を固定しているのだろう。
こんなにきれいな式になるなら、PRML 10.2.2 はどうして L を整理した式が書かれていないんだろう。

次回は「縮退するぞ！」の話の予定。