PRML 10章の変分ベイズによる混合ガウス分布推論の検証（フォロー編)

「パターン認識と機械学習(PRML)」10.2 章に従って変分ベイズ(Variational Bayes, VB)を R で実装してみて、PRML に書いてある内容通りか確認してみたところ、なんか違う。
「『変分混合ガウス分布は、余った混合要素は勝手にゼロになるから K が大きくてもいいよ』とか書いてあるけど全然縮退しないよ。ベイズ職人でないとうまくいかないらしいよ」
「初期値について『対称性から、通常 m_0=0 とおく』と書いてあるけど、ほんとに m_0=0 にしたら、全パラメータが k に対して同じ値になっちゃうよ」
と言いふらしていたら、スクリプトのバグだった。

ので、罪滅ぼしにまじめにもうちょっといろいろ検証してみたよ、というお話。

経緯

変分ベイズ実装(PRML 10.2): https://shuyo.hatenablog.com/entry/20100306/variational
PRML 読書会 #12 ９章ＥＭアルゴリズム＆１０章変分ベイズ: https://shuyo.hatenablog.com/entry/20100309/vb

事の起こり。PRML 読書会で「縮退しないよ」「m_0=0 は孔明の罠だって」と言いふらす。

PRML 読書会 #13 「10.2 変分混合ガウス分布」資料(1): https://shuyo.hatenablog.com/entry/20100406/variational
PRML 読書会 #13 「10.2 変分混合ガウス分布」資料(2): https://shuyo.hatenablog.com/entry/20100410/mixture_gauss
PRML 読書会 #13 10章近似推論法(変分ベイズ): https://shuyo.hatenablog.com/entry/20100412/prml

次の PRML 読書会でも同じことを繰り返す。
が、現場で id:tsubosaka さんにスクリプトのバグを見つけてもらった(感謝)。バグをつぶすと、「縮退するじゃない！」
というわけで、ウソでした、ごめんなさい、と謝る。

むしろきれいに縮退しすぎて N_k が０に → N_k で割るところで NaN が頻発 → 結構あちこち直さなきゃ。
初期値についても、m_0=0 にして、初回の m_k だけランダムに散らしてみた → うまくいく*1

混合ガウス分布の変分下界の計算式: https://shuyo.hatenablog.com/entry/20100419/lowerbound

PRML 10.2.2 の変分下界 L の式が膨大なのでさぼって実装していなかったのだが、これもやっぱりまじめに検証しないとあかん。
で、実装する前に式を整理してみたら、めちゃめちゃ簡単できれいな式になった(ある程度は整理できるだろうと思っていたが、まさかここまでとは)。
最初からこの式も載せといてくれたら、絶対に実装さぼったりしなかったのに……(逆恨み)。

混合ガウス VB 検証

検証に用いた R スクリプトはこちら。

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

このスクリプトは、Old Faithful に対して以下のように初期値を変えつつ推論を行う。

K: 2 から 6 まで
α_0: コマンドラインから与える
β_0: 0.05 から 1.00 まで 0.05 刻み
ν_0: D-0.5 から D-9 まで 0.5 刻み(D はデータ点の次元)
m_0: つねに 0
W_0: つねに単位行列

推論は変分下界 L の上昇が 0.0001 を下回ったら停止。
各初期値について推論を一定回数行い、変分下界がもっとも大きかった結果を summary として出力する。
summary には PRML 10.2.4 に従い、変分下界に ln K! を加えた値もあわせて出力する。

このスクリプトを用い、α_0 に対して、PRML 10.2.1(p193) に記載のある 0.001, 1, 10 を与え、結果を検証してみた。
例えば、K=6, α_0=0.001, β_0=0.10, ν_0=10.0 のときの推論結果がこちら。

#1638-1: convergence=199, L=-428.065, N_k = 3(7.292, 0.000, 169.119, 0.000, 0.000, 95.589)
#1638-2: convergence=53, L=-428.065, N_k = 3(169.119, 95.593, 0.000, 0.000, 7.289, 0.000)
#1638-3: convergence=212, L=-428.065, N_k = 3(0.000, 95.588, 7.291, 0.000, 0.000, 169.122)
#1638-4: convergence=412, L=-428.065, N_k = 3(0.000, 169.115, 0.000, 0.000, 95.590, 7.295)
#1638-5: convergence=27, L=-422.168, N_k = 2(0.000, 96.909, 0.000, 0.000, 175.091, 0.000)
#1638: K=6, alpha=0.001000, beta=0.100000, nu=10.000000, max_L=-422.168, p(D|K)=-415.588

各行で N_k として表示されている値は、先頭が消えずに残った混合要素数、続けて各混合要素の N_k の値(=その混合要素に負担されているデータ点の個数の期待値。0.5 未満になると消えたと判断)が並ぶ。convergence は収束回数。
推論を５回行い、残った混合要素数は 2 または 3 、変分下界 L はその残った要素数によって決まった値に収束する、といったことが読み取れるだろう。また各 N_k も、順番は違うが、ほぼ同じ値に収束している。
ちなみに残った混合要素数が小さいほど L は大きくなる傾向にあるが、必ずそうなると決まっているわけではない(頻度は少ないが、逆転するパターンもあった)。

こういったデータが 6000件あるので、いろいろ遊べる。
パラメータに対してどういう傾向があるのかを知るのが、ベイズ職人への第一歩。

まずは「縮退するっぽい」じゃあなくて、本当に縮退するのかどうか。
α_0 ごと、K ごとに残った混合要素それぞれが何件ずつあったかをまとめたのがこちら。

α_0=0.001 では混合要素２個または３個に、α_0=1 では混合要素３個に、α_0=10 ではほとんど縮退せず、ということで、まさに PRML p193 に書いてある通りの結果に。
ちょっと気になるのは、α_0=0.001, K=6 で３個になる場合が多いこと。これは他のパラメータの初期値によってコントロールできるのでは？
β_0 の違いはほとんど影響がなかったが、α_0=0.001, K=6 に対し、ν_0 ごとに残った混合要素数を見るとこちら。

ν_0 が小さいと２個に、大きいと３個になりやすいことがわかる。
これは ν_0 はウィシャート分布の自由度であり、それが大きいということはΛの自由度が大きくて、Λが大きくなりやすいということは混合要素の分散が小さくなりやすくて、だから小さくまとまった混合要素が多めにできやすい……のかな？
ウィシャート分布の自由度はもともとの定義では整数であることを考え合わせると、ν_0 = D または D+1 あたりを用いるのが適していそうだ。

が、実は ν_0 が大きい方が変分下界の収束先は大きくなる。
ν_0 ごと、K ごとの変分下界の平均がこちら。

ν_0=9 あたりで L が最大になる。もちろん K は 2 が最適だ。
というわけで、ベイズなら関連度自動決定はちゃんと働くことがわかったとは言え、それに頼りきるのではなく、できれば最適な K をきちんと評価推定することは(きっと)大切ってことで。

関連度自動決定に頼り切るわけにはいかないかも、というのはもう１つ。
PRML 10.2.4 では、K が異なるモデル間の尤度を比較するには、混合要素の入れ替えを考慮して、尤度(≒変分下界で近似)に ln K! を加えたものをモデルエビデンスとすれば良い、とある。

というわけで α_0=0.001, β_0=0.10, ν_0=10.0 という同じ初期値で、K=2 と K=6 で推論した結果がこちら。p(D|K) が変分下界に ln K! を加えたものである。

#38-1: convergence=12, L=-421.044, N_k = 2(175.091, 96.909)
#38-2: convergence=21, L=-421.044, N_k = 2(175.091, 96.909)
#38-3: convergence=4, L=-567.607, N_k = 1(0.000, 272.000)
#38-4: convergence=13, L=-421.044, N_k = 2(175.091, 96.909)
#38-5: convergence=8, L=-421.044, N_k = 2(175.091, 96.909)
#38: K=2, alpha=0.001000, beta=0.100000, nu=10.000000, max_L=-421.044, p(D|K)=-420.351

#1638-1: convergence=199, L=-428.065, N_k = 3(7.292, 0.000, 169.119, 0.000, 0.000, 95.589)
#1638-2: convergence=53, L=-428.065, N_k = 3(169.119, 95.593, 0.000, 0.000, 7.289, 0.000)
#1638-3: convergence=212, L=-428.065, N_k = 3(0.000, 95.588, 7.291, 0.000, 0.000, 169.122)
#1638-4: convergence=412, L=-428.065, N_k = 3(0.000, 169.115, 0.000, 0.000, 95.590, 7.295)
#1638-5: convergence=27, L=-422.168, N_k = 2(0.000, 96.909, 0.000, 0.000, 175.091, 0.000)
#1638: K=6, alpha=0.001000, beta=0.100000, nu=10.000000, max_L=-422.168, p(D|K)=-415.588

もちろん K=2 の方がよいモデルであるという結果を期待するわけだが、上の結果では K=6 の p(D|K) の方が大きくなってしまっている。

ふと気付くのだが、N_k = 0 となってしまった混合要素については、もはや区別が付かないのだから、入れ替えても「同等だが別の設定」が出来たりはしない。PRML p197 にも「(縮退している特別な場合を除いて)」とさりげなく書いている。
つまり、縮退してしまっている場合は ln K! ではなく $\ln {}_K P_r$ を足すべきでは？(ただし r は非ゼロの混合要素数)
まあ、上のケースではそのように変えてもまだ K=6 のエビデンスの方が少し大きいのだけれど……でも、他のケースでは K=2 の方がちゃんと良くなることがままあった。

他にも面白い切り口が色々あるかもしれない。β_0 まわりはまだ見てないし。
んが、きりがないので、summary データを github に commit しておいた。

https://github.com/shuyo/iir/blob/master/unsupervised/vb_result.csv

フィールドは左から K, α_0, β_0, ν_0, 変分下界, モデルエビデンス(の近似値。=L + ln K!), 非０な混合要素数, モデルエビデンスその２(= L + ln _K P _r)。
Excel で読み込んでピボットテーブルでうはうはすると楽しいかもしれない。

なお、今回のスクリプトは３次元以上のデータにも対応できるよう汎用的に作られており(Old Faithful しか食わせたことないが)、PRML と対応させやすいように最適化もしていない。
そのため、この 6000 件のデータを取るのに3晩かかるくらい、とても遅い。

*1:バグを直す前も同様のことを試したが、そのときは m_k が常に 0 に収束する、という結果だった

経緯

混合ガウスVB 検証

混合ガウス VB 検証