作って遊ぶ機械学習。

~基礎的な確率モデルから最新の機械学習技術まで~

ベイズ混合モデルにおける近似推論③ ~崩壊型ギブスサンプリング~

 さて、今回は混合モデルに対する効率的な学習アルゴリズムとして、崩壊型ギブスサンプリング(Collapsed Gibbs Sampling)*1を紹介します。ベースとなるアイデア自体はギブスサンプリングそのものなのですが、標準的な方法と比べて種々の計算テクニックが必要になってくるのと、実際に推論精度が向上することが多いのでここで紹介したいと思います。 

また、今回の実験で使用したソースコードGithubにおいてあります。

GitHub - sammy-suyama/MLBlog: source codes used in MLBlog

 

[必要な知識]

下記をさらっとだけ確認しておくといいです。

 

1、崩壊型ギブスサンプリングのアイデア

基本的なアイデアとしては、対象の確率モデル(同時分布)から一部の確率変数をあらかじめ周辺化除去(marginalise out)しておき、より数の少ない確率変数をサンプルすることによってサンプリングの効率性を高めようとするものです*2

ここでは、崩壊型ギブスサンプリングの位置付けを明確するために、いくつかの異なるギブスサンプリングの手法を比較紹介します。

ある確率分布{p(z_1, z_2, z_3)}があり、この分布から{z_1},{z_2},{z_3}をすべて同時にサンプルするのが不可能であるとしましょう。この様な問題に対して、一般的なギブスサンプリングは次のように確率変数を1つ1つ別々にサンプリングするという作戦を取ります。

f:id:sammy-suyama:20160729175801p:plain

これは「直接同時にサンプルするのが難しいのであれば、他の値を固定したより簡易な条件付き分布から順次サンプルすればいい」というアイデアなのでしたね。

次に、ブロッキングギブスサンプリング(blocking Gibbs sampling)というのもあり、この場合はいくつかの変数をまとめてサンプリングします。

f:id:sammy-suyama:20160813140837p:plain

「すべて同時にサンプルするのは難しいが、{z_1},{z_2}の同時分布なら簡単に計算できる」というような場合ではこちらを使った方が断然性能が良いです。

さて、本題の崩壊型ギブスサンプリングでは上記2つの手法とは異なる作戦を取ります。この手法では、モデル中のいくつかの変数を周辺化除去することによって、より少ない数の変数だけをサンプルします。(Collapsed spaceでサンプルする、と呼ぶようです。)

f:id:sammy-suyama:20160813142357p:plain

f:id:sammy-suyama:20160813142409p:plain

最初の{z_3}に対する積分除去が解析的に出来ることと、それによって得られた同時分布{p(z_1,z_2)}からの交互サンプリングが簡単に行えることが適用できる条件です。単純なギブスサンプリングと比べて、そもそもサンプルする変数の種類が少ないので、より速く目的の分布からのサンプリングが得られることが期待できます。

一般的にはブロッキングギブスサンプリングや崩壊型ギブスサンプリングの方が、単純なギブスサンプリングより性能が良い傾向にあります。しかし単純ギブスサンプリグの方が1回あたりの計算が軽いので、問題によってはより早く所望の結果が得られる場合があります。実用上では、まず最初の2つの手法を検討してみて、それでも精度がイマイチだった場合には崩壊型ギブスサンプリングを試してみるのがオススメです。*3

 

2、混合モデルからパラメータを周辺化除去

まずは標準的な混合モデルを考えてみましょう。ここでは{x}{N}個の観測データの集合、{z}が各データに対する隠れ変数の集合、{\theta}{x}の分布を制御するパラメータ、{\pi}が隠れ変数の分布に対するパラメータです。f:id:sammy-suyama:20160813144049p:plain

{p(\theta)p(\pi)}がパラメータに対する事前分布です。一般的に混合モデルでは、{p(z)}はカテゴリカル分布を使います。{p(x|\theta, z)}{x}の発生をモデル化するための確率分布であり、今回は例としてポアソン分布を使いますが、用途次第ではガウス分布でも多項分布でも構いません。

さて、この式に対するグラフィカルモデルを次の様に表現してみることにします。

f:id:sammy-suyama:20160813151130p:plain

混合モデルに対する崩壊型ギブスサンプリングでは、パラメータをモデルからあらかじめ周辺化除去するのが一般的なやり方です。

f:id:sammy-suyama:20160813144102p:plain

こんな感じに積分除去すると、対応するグラフィカルモデルは次のようになります。一応これをここでは周辺化モデルって呼ぶことにします。

f:id:sammy-suyama:20160813151148p:plain

当然、積分除去された{\theta}{\pi}はもはや周辺化モデルには現れません。しかし、そのおかげでモデルが簡単になったかというとそうでもないようですね。隠れ変数{z}は元のモデルでは({\pi}が与えられたとき)独立性を保っていたもの、{\pi}積分除去された途端にすべての隠れ変数が依存関係を持つようになってしまいました。同様に、{x}{\theta}積分除去されたために、すべてのデータ点がお互いに依存するようになってしまいました*4 。このようなモデルから{z}をサンプリングする計算式を導くのが目標になります。

 

3、周辺化モデルからのサンプリング

さて、周辺化モデルは非常にタイトな変数依存関係を持ってしまってしまっているので、{z}を一度に同時にサンプルすることはできません。しかし{z_1}{z_2},・・・、{z_N}のように、1つずつサンプルすることが可能である場合があります。ある{z_n}を、{z_n}以外のサンプル{ z_{\backslash n} = \{ z_1, ..., z_{n-1}, z_{n+1}, ..., z_N \} }がすでに与えられた状態でサンプルできないか考えてみることにします。この条件付き分布は次のように計算することが出来ます。

f:id:sammy-suyama:20161103130033p:plain

1行目では、ただ単にベイズの定理の分母を無視したもので事後分布を書き直しただけです。2行目ではそれをバラバラに展開しました。真ん中の項は{x_n}を除いたグラフィカルモデルを考えれば、独立性から、

f:id:sammy-suyama:20161103164208p:plain

であることが示せるため、結果今サンプルしたい{z_n}とは無関係になるので3行目では除去されています。

さて、出来上がった最後の式をよく見てみると、目的の確率分布が{p(x_n |x_{\backslash n}, z_n, z_{\backslash n})}{p(z_n | z_{\backslash n})}の2つの分布に分解されることがわかりますね。これらの確率分布はどういう意味を持っているのでしょうか。

まず始めに{p(z_n | z_{\backslash n})}ですが、次のような{z_n}に関するある種の予測分布であると解釈することができます。

f:id:sammy-suyama:20161103130333p:plain

まず、すでにサンプルされていると仮定している{z_{\backslash n}}を観測データかのように扱うことで、{\pi}のある種の事後分布のようなもの{p(\pi|z_{\backslash n})}が計算できます。この事後分布は、次の様にベイズの定理を使って簡単に計算することができます。

f:id:sammy-suyama:20161103130508p:plain

今回は{p(z_{\backslash n}|\pi)}はカテゴリカル分布で、{p(\pi)}はディリクレ分布です。したがって、次のようなディリクレ分布事後分布が求まります。

f:id:sammy-suyama:20161103130630p:plain

さて、この{\pi}に関する事後分布を使って予測分布を計算してみることにします。結局、有限シンボルの離散分布であることがわかるので、結果はカテゴリカル分布として表現することができます。

f:id:sammy-suyama:20161103212604p:plain

同様にして、残る{p(x_n |x_{\backslash n}, z_n, z_{\backslash n})}の項も次の様に{x_n}に関する({z_n}の条件付き)予測分布になることがわかります。

f:id:sammy-suyama:20161103131603p:plain

次のセクションで解説しますが、ポアソン・ガンマモデルではこの積分は解析的に実行できます。*5

ところで{x_n}は実際にはデータとして与えられているので、より正確に言うとこの項は予測尤度とでも言った方が正しいかもしれません。肝心な点は、具体的な{z_n}を与えることによってこの尤度が評価できるという点です。{z_n}に対する{K}種類の実現値1つ1つに対して式{p(x_n |x_{\backslash n}, z_n, z_{\backslash n})p(z_n | z_{\backslash n})}を評価し、正規化してあげれば、{z_n}をサンプルするための離散分布が計算できることになります。

 

4、ポアソン・ガンマモデルにおけるベイズ推論

さて、ここまでで崩壊型ギブスサンプリングの説明はほぼ終わりなのですが、残された課題は次の条件付き予測尤度の計算ですね。

f:id:sammy-suyama:20161103131603p:plain

先ほどと同様、手順としてはまず右側の事後分布{p(\theta|x_{\backslash n},z_{\backslash n})}を計算し、その後で尤度関数{p(x_n | \theta, z_n)}と掛け合わせてから{\theta}積分除去します。

まず事後分布の計算からやってみましょう。今回は{D}次元のポアソン・ガンマモデルを仮定しているので{\theta}ポアソン分布のパラメータ{\lambda}になります。こちらもベイズの定理を使うと、

f:id:sammy-suyama:20161103134221p:plain

として計算をすることができ、事後分布は当然ガンマ分布になります。

f:id:sammy-suyama:20161103134239p:plain

さて、この事後分布を使って予測分布を計算してみましょう。幸い事後分布がガンマ分布のままなので、このまま解析的にパラメータ{\lambda}積分除去することができ、得られる分布は次のような負の二項分布(Negative Binomial Distribution)になります。

f:id:sammy-suyama:20161103213242p:plain

ここでは、負の二項分布の定義として次のような表現を使いました。

f:id:sammy-suyama:20161103135740p:plain

混合モデルにおいて崩壊型ギブスサンプリングが適用できる条件としては、観測モデルの事後分布と予測分布がともに解析的に実行できることです。今回のようなポアソン・ガンマモデル以外でも、共役の関係にある分布を選べばこれは必ず実行できます。  

 

5、動作確認実験

さて、今回導いたポアソン混合モデルにおける崩壊型ギブスサンプリングアルゴリズムをJuliaで実装し、動作確認をしてみます。ソースコードGithubにあります。

観測データ数N=1000、次元数D=2、クラスタ数K=8として混合モデルを学習させてみました。次の図は100回までのサンプル結果をアニメーションにしています。 

f:id:sammy-suyama:20161028181950g:plain

左下の方のクラスタはかなり密集していますが、なんとかうまく分離できているように見えますね。右図でサンプル回数が50回を超えてもまだパタパタしている箇所がありますが、それらのデータ点はクラスの所属が「あいまい」であることを示しています。

 

さて、せっかくですので前回導出した標準的なギブスサンプリングと簡単な性能比較をしてみます。下の図は横軸をサンプル回数、縦軸を最後のサンプルに対する尤度にしてグラフを描いたものです*6。GS(ギブスサンプリング)、Collapsed(崩壊型ギブスサンプリング)ともに5回ずつ実験しています。

f:id:sammy-suyama:20161028121122p:plain

傾向として、崩壊型ギブスサンプリングの方が少ない繰り返し回数で、より高い尤度に達しているようですね。最終的な到達点も崩壊型ギブスサンプリングの方が高い傾向があります。

ちなみに計算時間で比較すると、100回のサンプリングでギブスサンプリングが平均7.03秒、崩壊型ギブスサンプリングが平均7.15秒と、あまり差がないようです。 一般的には、崩壊型ギブスサンプリングにおける予測尤度の計算が遅くなる可能性がありますが、今回の場合は特にボトルネックにはなっていないようです。

 

6、まとめ

さて、今回はベイズポアソン混合モデルを題材にして、崩壊型ギブスサンプリングのアルゴリズムを導いてみました。今回の簡易実験が示す通り、標準的なギブスサンプリングアルゴリズムよりも性能が高くなる傾向があるようです。オススメとしては、まず標準のギブスサンプリングを導いてみて、その性能で満足いかないようであれば崩壊型ギブスサンプリングを導出してみるのが良いかと思います。ただし、解析的にパラメータを積分除去するためには、今回のポアソン分布とガンマ分布の例のように、観測モデルが共役の関係性を持っている必要があります。

同様のアイデアは時系列モデルである隠れマルコフモデル(Hidden Markov Model, HMM)、トピックモデルにおけるLatent Dirichlet Allocation(LDA)に適用されているほか、中華料理店過程(Chinese Restaurant Process, CRPインド料理店過程(Indian Buffet Process, IBP)などのベイジアンノンパラメトリクスの手法にも応用されています。

 

*1:参考: Machine Learning: A Probabilistic Perspective by Kevin Murphy | The MIT Press

*2:このことから、崩壊型ギブスサンプリングは周辺化ギブスサンプリングとも呼ばれることもあるようです。

*3:ただし、モデルによっては解析的にうまくアルゴリズムを導けない場合や、解析的な周辺化のために積分近似を行う場合もあります。

*4:これは一般的には{\int \prod_{n} p(x_n|\theta)p(\theta)d \theta \neq \prod_{n} \int p(x_n|\theta)p(\theta)d \theta}であるためです。

*5:他にもガウス・ウィシャートモデルのように、共役関係が成り立っていれば解析計算できます。

*6:尤度で性能を比較するのはあまり望ましくないですが、今回は完全に同じモデルからアルゴリズムを導いているので、尤度でもほぼ「良いクラス割り当て」を反映しているとしてしまってもいいんじゃないかと思います。