トップ > 機械学習 > Multinomial DPMを実装してみた

Multinomial DPMを実装してみた

機械学習 ベイズ統計 C++

ちょっと前に実装してたんだけど、メモを書くがてら公開してみる。やりたいこととしてはnested Chinse Restaurant Processまで行きたいんだけど、ノンパラベイズ初心者なので一番取りかかりやすいであろうDirichlet Process Mixture(DPM)を文書モデルでやってみたという感じです。HDPではなくDPMとしてコーパスをモデル化するので、いくつ文書があろうがそれらは全部一つの文書として取り扱います(というか、そういう形でしか取り扱えません。扱いたかったらHDPの世界へ行こう)。

目的

やってみようと思った理由はいくつかあって

  • 実際に自分で把握できるミニマムな*1ノンパラベイズのプログラムをgetする
  • 実際に書いてみることでノンパラベイズのプログラムではまりやすいところを知る
  • コード書く段階までやってみないと分かったつもりになっていることが多い

などなどです。実際に動かしてみるとサンプリングとかよりクラスタ削除した際のラベルの付け換えに一番時間くったりしているっぽいけど、とりあえず動いているのでまあいいことにしておこう(ぇ。

DPMの簡単な復習

DPMのグラフィカルモデルは以下のように書き表わすことができる。
http://www.cs.brown.edu/~sudderth/papers/sudderthPhD.pdf
Hが基底分布で、今回は文書を取り扱うのでHは具体的ディリクレ分布となる。\alphaは集中度パラメータで、こいつが大きいほど新しいクラスタができやすくなる。今回コードを書いて実験したのはこのパラメータの影響度が実際にどんなもんかを身を持って確かめたかったというのもある。(連続分布である)Hから生成した分布Gは離散分布となる。Gが離散分布となることからサンプルx_i毎にパラメータ\theta_iを用意しておいても実際には同じ値が使われやすいような効果が起きる。

モデルとしてはこれで何も不足していないのだが、実際にGibbs Samplingをするときには補助変数を導入したり、基底分布を具体的に書き換えたほうがやりやすいことが多いので書き換える。
http://www.cs.brown.edu/~sudderth/papers/sudderthPhD.pdf
文書の区切りがないようなLDAでトピック数は無限大といった感じのイメージと言えば伝わりやすいだろうか。ベイズなので、パラメータ\piと\theta_kについて積分消去したものを考える。LDA Gibbsのときと同様に考えて、サンプリングの式は以下の2つのfactorに分解して考えることができる。
p(z_i | \mathbf{z}_{-i}, \mathbf{x}, \alpha, \lambda) \propto p(z_i | \mathbf{z}_{-i}, \alpha) p(x_i | \mathbf{z}, \mathbf{x}_{-i}, \lambda)
2つのfactorを、文書の場合について書きくだしてみたのがこれ。新しいクラスタができるところも実際に書きくだしてみるとそんなに怖くない、ということが分かる。

対数尤度

Gibbs Samplingの最中は対数尤度がある程度単調に増加していくことを確かめたいので、その式の導出。まず、事後分布は
p(\mathbf{z} | \mathbf{w}, \alpha, \gamma) \propto p(\mathbf{w} | \mathbf{z}, \gamma) p(\mathbf{z}, \alpha)
というように2つのfactorに分解することができる。最初のfactorについてはLDAの要領で、以下のポリア分布が出てくる形で書ける。
p(\mathbf{w} | \mathbf{z}) = \int \prod_{i=1}^N p(w_i | z_i, \phi) p(\phi | \gamma) d \phi = \left(\frac{\Gamma(V \gamma)}{\Gamma(\gamma)^V}\right)^K \prod_{k=1}^K \frac{\prod_{v=1}^V \Gamma(N_{k, w} + \gamma)}{\Gamma(N_k + V \gamma)}
ここで、N_{k, w}はクラスタkに属している単語wの数。次に第二項だが、ここはLDAのように行かない。CRPの定義を生かすために隠れ変数のN個の同時分布をchain ruleを使ってばらすと以下のように書き表わすことができる。
p(\mathbf{z}) = p(z_1) p(z_2 | z_1) \cdots p(z_n | z_1, \cdots, z_{N-1}) = \frac{\alpha^K \prod_{k=1}^K (N_k - 1)!}{\alpha (\alpha + 1) \cdots (\alpha + N - 1)}
ここで、N_kはクラスタkに属している単語の数。プログラムを書くときにはこの対数のものが分かっていると便利なのでメモ。階乗とかが出てくるので3秒くらい焦るが普通にやればできる。
\log{p(\mathbf{z})} = K \log{\alpha} + \sum_{k=1}^K \log{(N_k - 1)!} - \sum_{i=1}^N \log{\alpha + i - 1} = K \log{\alpha} + \sum_{k=1}^K \sum_{i=1}^{N_k - 1} \log{i} - \sum_{i=1}^N \log{(\alpha + i - 1)}

テストセットパープレキシティ

新しいデータが新しいクラスタに所属するかそうでないかのどっちにすればいいかについて結構悩むのだが、Bleiの"Variational methods for the Dirichlet Process"とかを見ていると"note that the next component can take on K + 1 possible values"と書いてあるので新しいデータがK+1番目のクラスタを作ることを考慮してテストセットパープレキシティを導出してみる。

テストセットパープレキシティの定義式は
\mbox{PP}(D_{\mbox{test}} | D_{\mbox{train}}}) = \exp{\left(- \frac{1}{N} \sum_{n=1}^N \log{p(w_n)}\right)}
であるが、これを隠れ変数を周辺化したものだと考えると
\mbox{PP}(D_{\mbox{test}} | D_{\mbox{train}}}) = \exp{\left(- \frac{1}{N} \sum_{n=1}^N \sum_{z=1}^{K+1} \log{p(w_n, z)}\right)} = \exp{\left(- \frac{1}{N} \sum_{n=1}^N \sum_{z=1}^{K+1} \log{p(w_n | z) p(z)}\right)} = \exp{\left(- \frac{1}{N} \sum_{n=1}^N \left(\left(\sum_{z=1}^K \log{\frac{N_{z, w} + \lambda}{N_z + V \alpha} \frac{N_z}{N + \alpha}}\right) + \log{\frac{1}{V} \frac{\alpha}{N + \alpha}}\right)\right)}
と書ける(p(z)のところがちょっと自身ない...)。

実験

コードはgithubに上げてあります(DPMなのにdirichlet_processというプロジェクトの名前のままだった...)。

ややこしいですが、このエントリでのハイパーパラメータとコードでのハイパーパラメータがちょっと違ってて

  • エントリ => コード
  • \alpha => \gamma
  • \lambda => \beta

という風になっているので注意。プロの方々はハイパーパラメータも最適化するらしいですが、今回はそれは置いとく。

実験の結果では、iteration毎の

  • クラスタの数
  • 対数尤度
  • テストセットパープレキシティ

を表示しています。クラスタ数は大体のところで頭打ちになっていて、対数尤度も途中から安定してるので多分当ってて、テストセットパープレキシティもほどほどな値になっている。ということできっと動いているんでしょう、たぶん。
Quartz 2 [*]
ただ、バグっているのかそもそもこういうものなのかよく分かっていないところが一点。p(w|z)で各クラスタで特有な単語を列挙してみた際にそれぞれのクラスタ間であまち違いが出ていないという結果になってしまう。なんでだろー。

*1:DPを使ったコードでいくつか公開されているものがあるけど、DPM以上に複雑なものが多いので、把握しずらいことが多かった。あとmatlabで公開されているものも多かったが、C++での例が欲しかった