作って遊ぶ機械学習。

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

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

 さて、今回は混合モデルに対する効率的な学習アルゴリズムとして、崩壊型ギブスサンプリング(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:尤度で性能を比較するのはあまり望ましくないですが、今回は完全に同じモデルからアルゴリズムを導いているので、尤度でもほぼ「良いクラス割り当て」を反映しているとしてしまってもいいんじゃないかと思います。

機械学習の4つのアプローチ

おつかれさまです。今日はちょっと趣を変えて、近年のいわゆる「機械学習」という技術のアプローチをカジュアルに少しカテゴリ分けしたいと思います。

といっても、自分はアカデミックの研究者ではなく大量の論文を読み漁るということもほとんどしないので、理論的なバックグラウンドに基づいたソリッドなカテゴリ分けはできません。ここで紹介するのはあくまで、実用上の機械学習技術者から見た視点で「こんな傾向があるかなぁ」くらいの気持ちで書いたものです。 

 

<代表的な4つのアプローチ>

1、最適化(目的関数ベース)

まず始めは最適化手法をベースにした機械学習アルゴリズムです。たぶん一番例が多いんじゃないでしょうか。

ここでは、ある解きたい課題を目的関数によって定式化し、適切な最適化手法を使って解きます。伝統的な線形回帰や線形識別はもちろん、主成分分析(PCA)や非負行列因子分解(NMF)なんかもこの枠組みで解かれることが多く、ほとんどの場合で誤差関数を定義してそれを最小化するというアプローチを取ります。他にはサポートベクターマシンSVM)なんかもマージン最大化という基準にしたがって最適化を行いますね。

この種のアルゴリズムの利点としては、比較的自由に目的関数をデザインして解くことができる点が挙げられると思います。

欠点としては、ある種の複雑な課題に対してはそもそもどうやって目的関数を設定していいかわからない場合があることと、最適化に基づいているので簡単にデータに対してオーバーフィットしてしまうことです。また、予測が「どれだけ確からしいのか」といった、予測に関する不確実性・信頼性に関する情報を上手く扱えないことが多いです。

 

2、深層学習

最近流行っているあれです。一番シンプルな例は(伝統的な)ニューラルネットワークなどで、アイデアとしては、データの表現を多層構造を使って学習することにより、データの次元間にまたがった共通の特徴を高レベルの層でシェアして汎化性能を高めようというものです。画像や音声の認識タスクで性能を発揮するようです。

利点としては、概念が比較的わかりやすいこと、ある特定の課題に対して顕著な性能を発揮していることです。

欠点としては、基本的にものすごい量のデータと計算コストを使用しなければならないことと、アルゴリズムのパフォーマンス向上(モデル設計やチューニング)に系統だった手段がないことです。あと、だいたい単純な誤差関数の最小化とかを基準にするので上の方法と同様オーバーフィットしやすいです。あとは、現状あまり理論面が洗練されていないのも課題です。

 

3、機械学習ライブラリ

これもカテゴリに入れていいのか謎ですが、実際の現場ではよく使われています。Scikit-learn*1などの機械学習ライブラリを駆使して問題を解くという、あらゆる場面で比較的にスピーディーに試せる方法です。

基本的なアプローチとしては、入力データをよく眺めて、特徴量抽出や次元削減などを行い、その結果を識別器にかけるというものです。

利点としては、ややこしい数学の知識などがなくてもプログラミングさえ少しできれば誰でも始められることです。機械学習の入門として、あるいはクイックプロトタイピングとしてお勧めかもしれません。あとは、わりと洗練されたライブラリが世に出回っているので、自分でいろいろアルゴリズムを開発してバグに悩まされるという心配もあまりありません。

欠点としては、解けるタスクが非常に限られていることです。問題自体をデザインすることがほとんどできないので、最大性能もほかのアプローチと比べて限界があります。


4、ベイズ学習

すべての問題を確率モデリングと確率推論の2ステップで解くアプローチです。

基本的には、データの発生プロセスを確率分布を使って物理的・論理的に記述し(モデリング)、計算機を使ったアルゴリズムによって事後分布を計算する(推論)という流れです。上に具体的に挙げたような各アルゴリズムの多くがベイズ化できるほか、ベイジアンノンパラメトリクスと呼ばれる、ガウス過程やディリクレ混合過程などといったベイズ手法ならではのモデルもあります。

このアプローチの利点としては、基本的にはオーバーフィットをしないため*2、少量・多次元のデータでも性能が出やすいことです。 また、確率分布を組み合わせてモデルを自由にデザインできるので、幅広い課題に対してアルゴリズムを与えることができます。

欠点としては、アルゴリズムの導出に多少の数学力を要することと、厳密な計算には非常に計算コストがかかることです。

 

 

 <ベイズ学習の汎用性>

で、自分は完全なベイズ派なので、最後に今一度ベイズ学習の汎用性と、他の3手法との関係性について触れたいと思います。

 

1、ベイズ学習と最適化

ベイズ学習でも最初に挙げた最適化手法をサブルーチンとして使うことは良くあります。というのも、ある種の最適化問題は、ベイズ学習で必要となる積分の近似と等価になるからなんですね。

例えば最尤推定(最小二乗法を使った誤差関数の最小化)なんかも、実はベイズの近似推論に対して極端な仮定を置くことで導かれます*3。主成分分析、非負行列因子分解もベイズ化できるほか、サポートベクターマシンガウス過程の関係性も指摘されています*4

 

2、ベイズ学習と深層学習

ベイズ学習は汎用的なデータ解析の枠組みであるので、もちろん深層学習の多様なモデルをベイズ学習として再定式化することができます。代表的には、ガウス過程を使ったディープなモデル*5があるほか、変分ベイズを使ったオートエンコーダの実装*6もあります。あとは畳み込みニューラルネットワーク(CNN)をベイズの枠組みで焼き直して変分近似で解くといった手法とかもあります*7

さらに,ユニット数が無限のニューラルネットワークガウス過程と等価になります.

深層学習はガウス過程 - 作って遊ぶ機械学習。

 

3、ベイズ学習とライブラリ

ベイズ学習の1つの欠点は、近似推論手法の開発に非常に数学力と手間がかかることです。この手間を省いてくれるものとして、種々の統計確率言語があります。代表的なものにはStan*8があり、モデルの記述さえすればあとは勝手に近似推論アルゴリズムMCMCや変分推論)を実行してくれます。残念ながら、複雑なモデルに対する推論性能は現在のところあまりよろしくないようです。

他にはベイズの枠組みで作られたアルゴリズムがそのままライブラリ化されている場合ももちろんあります。最近人気なのはガウス過程を使ったベイズ的最適化(Bayesian Optimization)で、深層学習のチューニングなんかに使われているようです。

 

- - -

さて、ちょっとした落書きのような記事になってしまいましたが、少しでもたくさんある機械学習技術のアプローチとその関係性がクリアになっていただければ幸いです。

*1:Scikit-learn: http://scikit-learn.org/stable/

*2:PRMLとかではオーバーフィットしないと言い切ってますが、実際の多くの場面では近似アルゴリズムを使うので、近似性能が悪ければオーバーフィットもアンダーフィットも起こり得ます。

*3:最尤推定、MAP推定、ベイズ推論 - 作って遊ぶ機械学習。

*4:Gaussian Process for Machine Learning: http://www.gaussianprocess.org/gpml/

*5:Deep Gaussian Process: http://www.jmlr.org/proceedings/papers/v31/damianou13a.pdf

*6:Auto enoding variational Bayes: https://arxiv.org/abs/1312.6114

*7:Bayesian CNN: https://arxiv.org/abs/1506.02158

*8:Stan: http://mc-stan.org/

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

さて、今回はMCMCの代表的な手法であるギブスサンプリングを使って、混合モデルによるクラスタリングを行いたいと思います。今回も前回の変分近似の記事と同様、ポアソン混合分布を具体的なモデル例として使っていきます。

 

[必要な知識]

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

 

・ギブスサンプリング

さて、ギブスサンプリングでは、サンプルを取りたい同時分布{p(z_1,z_2,z_3)}に対して、下記のような繰り返し手続きを用いて変数{z_1,z_2,z_3}を順にサンプリングします。

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

他の変数を既にサンプリングされた値で固定し、残りを確率分布からサンプルする、というのを繰り替えしていきます。それぞれの確率分布が、注目している変数をサンプリングするのに十分なほど簡単になっているのがギブスサンプリングが適用できる条件です。証明は今回省きますが、このような手続きで得られたサンプルは、数が十分に多ければ真の分布から得られたものとみなすことができます。

ギブスサンプリングは変分近似と比べ、細かい積分計算(期待値計算)を必要としない分だけいくらか変分近似より楽に導出できるかと思います。また、サンプリング手法が元々積分計算を近似するために開発されたという経緯もあり、変分近似と比べて若干適用できるモデルの範囲が広いです。変分近似とギブスサンプリングの簡単な比較はこちらの記事をご参照ください。

 

ポアソン混合モデルに対するギブスサンプリング

さて、前回と同様、お題としてポアソン混合モデルを使用します。利便上もう一度確率モデルを書いておくことにします。 

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

{x}{z}はそれぞれ観測データと隠れ変数の集合で、{\lambda}{\pi}はパラメータです。

それぞれのパーツは具体的に次の様に設定します。

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

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

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

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

早速このモデルの事後分布({x}が観測された後の残りの確率変数の分布)の推定にギブスサンプリングを適用してみましょう。今回はパラメータ{\{\lambda, \pi \}}と隠れ変数{z}を別々に分けることによって、十分容易にサンプリングが出来る確率分布が得られることを示します。

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

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

まずパラメータである{\lambda}{\pi}のサンプルから見てましょう。すべての隠れ変数{z}がすでにサンプルされたある値で固定された状態になっているので、パラメータの分布はちょこちょこっと計算すると次のようなガンマ分布とディリクレ分布の独立した積になることがわかります(表記を簡単にするためインデックス{i}は取り除きます)。

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

ガンマ分布に関しては、さらにそれぞれのクラスタ{k})と次元({d})ごとに独立にサンプルすればOKだということがわかりますね。実装上では、お使いのプログラミング言語のライブラリでガンマ分布とディリクレ分布からそれぞれ値をサンプルすればOKです。

さて、パラメータがサンプルされたので、今度は隠れ変数{z}をサンプリングしてみることにします。{\lambda}{\pi}が一旦サンプルされた値で固定されると、次の様に{z}はそれぞれデータ点ごとに独立にサンプリングできることがわかります。

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

このようにしてすべての{z_n}をサンプルし終えたら、また{\lambda}{\pi}のサンプリングに戻ります。

 

・簡単な実験

さて、導出したポアソン混合モデルのギブスサンプリングを使って、簡単なトイデータのクラスタリングをしてみましょう。可視化のためデータの次元{D}は2に固定し、データ数は300とします。サンプリングのイタレーションは50回とします。

次の結果はクラスタ{K=2}とし、多少クラスタ間のオーバーラップがあるデータから隠れ変数をサンプリングしてプロットしてみたものです。

f:id:sammy-suyama:20160730013101p:plainf:id:sammy-suyama:20160730012717g:plain

左の図が正解のクラスタリングで、右図はサンプリングのプロセスをアニメにしています。右図を見ると、サンプリング数が増えていっても2つのクラスタのちょうど中間くらいにあるデータに対する割り当てはずっと安定しません。サンプリング法では、このような事後確率の不確実性がそのままサンプルされる値のバラツキ加減によって表現されるわけですね。 

 

次はクラスタ{K=6}のデータを分離してみた結果です。

f:id:sammy-suyama:20160730012556p:plain f:id:sammy-suyama:20160730012834g:plain

サンプリングの途中で誤ったクラスタ分けに一度収束しかけていますが、30回目くらいでなんとか運良く脱出し、最終的には正しい結果が得られたようですね。ギブスサンプリングはじっくり待てば理論的には真の事後分布からのサンプルを得ることが出来るのですが*1、実際には局所的に尤度の高い位置にひっかかってしまい抜け出せなくなってしまう場合も多いです。

 

・まとめ

ということで、今回はポアソン混合分布に対するギブスサンプリングを導出し、簡単な実験をやってみました。ガウス混合分布と比べてパラメータ数も1つ少なく行列計算も必要ないので、実際に手で計算してみてアルゴリズムを導出してみることをオススメします。変分近似と比べてもいくぶんか導出は楽ですね*2

次回は1つ発展として崩壊型ギブスサンプリング(Collapsed Gibbs sampling)を紹介する予定です。

*1:真の事後分布が得られることと、想定している正解のクラスタリング結果が得られることとは同値ではないので注意してください。

*2:前回の変分近似の記事と結果の式を比べてみると、期待値計算を除いてほとんど同じになっていると思います。実を言うと、前回の記事のtexをそのままコピーして今回の記事の式を作っています。

ベイズ混合モデルにおける近似推論① ~変分近似~

今回から数回に分けてベイズ混合モデルの構築法と種々の近似推論(変分近似、ギブスサンプリング、崩壊型ギブスサンプリング)に関してお話ししたいと思います。

混合モデルの代表的なアプリケーションはクラスタリングですが、今回ご紹介するモデルの構築方法はクラスタリングにとどまらず、もっと凝ったモデルを構築するための基本パーツになるものです。例えば時系列情報を取り扱うためのHMM(Hidden Markov Model)自然言語処理でよく使われているLDA(Latent Dirichlet Allocation)などはこれらの混合モデルを拡張することにより得ることができます。レゴブロックのように、アイデア次第ではそれ以上に複雑で豊かな表現力を持ったモデルを構築することも可能です*1

さて、今回の記事ではまず始めに混合モデルの例としてポアソン混合モデルを導入します。そして、以前にご紹介した変分近似(=変分推論、変分ベイズを使うことによって簡単なデータのクラスタリング実験を行ってみたいと思います。

 

[必要な知識]

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

 

ポアソン混合モデル

ベイズ学習の枠組みでは、まず始めにデータがどのような確率分布に従って発生しているのかをモデル化する作業から入ります。さっそくではありますが、今回の例では次のようなベイズポアソン混合モデルを考えてみたいと思います*2。今回はデータ数を{N}、データの次元数を{D}クラスタの数を{K}とおくことにします。

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

まずこの式(1)ですが、これは今回のモデルのすべての登場人物たちとその関係性を表しています。{x}{D \times N}のデータ行列です。{z}はサイズが{K \times N}の隠れ変数で、{z_n^{(k)} \in \{0,1\}}かつ{\sum_{k=1}^{K}z_n^{(k)}=1}を満たします。この離散変数{z}が「データがどこから発生するのか」を表す「スイッチ」の役割になるんですね。{\lambda}ポアソン分布の平均値パラメータで、今回は{K}個それぞれのクラスタに対して{D}個の平均値パラメータが必要なので、サイズは{K \times D}になります。{\pi}は混合係数と呼ばれる{K}次元の確率変数で、それぞれのクラスタが選ばれる確率を示しています。

さて、それぞれのパーツを具体的に定義していきましょう。

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

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

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

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

式(2)は、混合係数パラメータ{\pi}に関する事前分布を表しています。ベイズ学習では推定したい変数はすべて確率分布を持つので、このようにパラメータも確率分布を使って表現するわけですね。今回{\pi}は後ほどカテゴリカル分布に対するパラメータに使われるので、ここではその共役事前分布であるディリクレ分布を指定しています。{\alpha}は超パラメータと呼ばれているもので、今回はある固定値を持つものとします。

同様に、式(3)はポアソン分布の平均値パラメータ{\lambda}の確率分布を表しており、{a}{b}はさきほどと同じように事前分布を決める超パラメータです。こちらはポアソン分布に対する共役事前分布を選ぶ必要があるので、今回はガンマ分布からサンプルされるようにモデル化されています。こういった指数分布族の共役関係は覚えておいて損はないです。

式(4)は混合モデルで最も大事な部分です。ここではカテゴリカル分布を使うことによって所属クラスタを示す変数{z}をサンプルしているわけですね。今回{z_n}は1 of K表現を用いることにします。例えば{z_n}が3番目のクラスタを示している時は{z_n = [0,0,1,0]}などと書きます。

最後に、式(5)でようやくデータ{x}が出力されます。若干、{K}に関する掛け算の部分がトリッキーなので注意してください。{z_n}はある{k}番目の要素だけ1をとり、それ以外は0になるので、{z_n}で指定された1つのポアソン分布だけが消されずに残るという「スイッチ」の仕組みになっています。{z_n}に何が格納されているかによって、{x}がサンプルされる分布が切り替わるわけですね。

 

さて、上式のモデルはデータの発生過程を表しているので、実際に人工的なデータをシミュレートすることができます。3回ほどプログラムを回して異なるデータを発生させてみると次のような感じになります。(D=2、K=4、N=300)

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

隠れ変数{z}のもつ値は色で判別できるようにしています。このように、我々の持っている事前の仮説(モデル)から具体的な例(データサンプル)をシミュレートすることができるというのがベイズ学習の大きな特徴です。

 

・変分近似による隠れ変数Zの推定

さて、我々がやりたいのはただ単にデータをシミュレートして遊ぶだけではありません。例えばクラスタリング問題を混合モデル上の推論問題として解釈すると、「(色分けされていない)データ{x}だけが与えられたときに、それぞれのデータ{x_n}に対応する隠れ変数{z_n}は何か?」ということになります。発生源{z_n}をデータごとに特定したいわけです。

さて、この問題ですが、なかなか一筋縄ではいきません。なぜかというと、一般的なクラスタリングの問題では隠れ変数{z}だけではなくそれぞれのクラスタのパラメータ(平均値など)も未知であるからです。この問題はベイズ学習では次のような事後分布を計算することに対応します。

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

この分布はパラメータと隠れ変数が複雑に絡み合った、非常にやっかいな同時確率分布になっています。なんとかして「答えの確率分布は、ぐちゃぐちゃです」とは言わずに、「点{x_n}クラスタ◯◯◯に属し、その平均値は△△△です」のように我々が解釈できる答えがほしいところです。

 

ということで、今回はこの事後分布を近似推定する手段として変分近似を導入することにします。変分近似の基本に関しては別の記事があるのでそちらを参考にしてください。

 変分近似(Variational Approximation)の基本(1)

ここでやりたいことは真の事後分布をもっと簡単な近似事後分布で表現するということです。ここでは真の(同時)事後確率分布を直接推定することは諦め、次のような分解された分布によって近似することを考えます。

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

「隠れ変数とパラメータの複雑な相関関係は無視しちゃおう!」ってアイデアになります。で、具体的にどうやって近似分布を求めるかなんですが、次のKLダイバージェンスを近似事後分布{q}に関して最小化する最適化問題を考えることにします。

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

KLダイバージェンスは2つの分布間の「距離」のようなものを測るんでしたね。分解されたより単純な確率分布を仮定して、KLダイバージェンスの基準でなるべく真の分布に近いものを求めてあげようというのが狙いです。

さて、このような分解の仮定を置くと、近似事後分布は次のように分解された分布ごとに繰り返し更新できることがわかります(過去の記事をご参照ください。)。 

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

ここで、{\langle \rangle}は期待値操作を行うことを表しています。ここまでくればあとはモデルの式(1)~(5)を代入して頑張って積分計算(期待値計算)をやれば更新式が導けます。

まず式(9)におけるパラメータの近似事後分布ですが、下記のようにポアソン分布の平均値パラメータ{\lambda}と混合係数{\pi}の分布は、運が良いことに、関連する項が自然に分解されます。

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

最初に定義した具体的な確率分布(式(2)~(5))を当てはめて計算をすると、それぞれの近似事後分布は次のようになります。

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

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

共役事前分布を導入したおかげで、結果が事前分布と同じ形状の確率分布になりましたね。

次に、肝心の隠れ変数{z}の近似事後分布を式(10)を使って計算します。

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

ここでも具体的な確率分布を当てはめて計算すれば,次のようなカテゴリカル分布が得られます。

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

実装上は{\hat{\pi}_n}の正規化を忘れないでください。

 

・簡単な実験

さて、適当に作った非負2次元データに対して今回の変分近似アルゴリズムを動かしてみましょう。(D=2、K=6,N=300)

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

左の図が正解のクラスタ割り当てを表し、右の図が変分近似による100イテレーション後の推定結果になります。今回はクラスタリングをすることが目的なので色の組み合わせまで当てる必要はありません。このくらい明確にクラスタ構造が見えるようなデータでは、変分近似はほぼ確実に正解にたどり着くことができます。

 

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

入力データを変えてもう1度やってみました。こちらのデータはちょっと難しかったかもしれません。左図の左下のクラスタにはかなりオーバーラップが見えますね。右図の推定結果は、正解とは異なる分け方になってしまっています。この場合はデータの数を増やすか、データの次元数(使えるデータの種類数) を増やすことによって精度向上が期待できます。

 

・まとめ

ということで、今回はベイズポアソン混合分布を例にとって、対応する変分近似アルゴリズムを導出し、簡単な実験を行ってみました。今回ご紹介したように、カテゴリカル分布を使った混合モデルは確率分布をスイッチする役割を果たします。ちょうどプログラミング言語におけるif分のようなものであると考えると、その重要性がおわかりいただけるのではないかと思います。

変分近似の導出のところは実際に手を動かしてやってみると良いです。行列計算や複雑な同時分布の計算がない分だけ、PRML等に載っているガウス混合モデルより簡単だと思います。

次回以降は、ギブスサンプリングや崩壊型ギブスサンプリングを今回と同じモデルに対して適用し、また簡単な実験をしてみたいと思います。

*1:例えば深層学習と呼ばれる技術(多層パーセプトロンオートエンコーダ畳み込みニューラルネットなど)もすべてベイズモデルで書き下すことができます。ちょっと手間がかかりますが。

*2:多くの教科書ではガウス混合モデルを紹介しているかと思いますが、多次元分布なのでちょっと計算が面倒なのと、同じ例では芸がないと思ったので今回はポアソン分布にしました。実際、非負のデータって結構いっぱいあるので便利なんじゃないかと思います。

ベイズ推論による機械学習の基本

今回は基本的なベイズ学習の概念と流れを説明したいと思います。まず始めに、ベイズ学習のすべての基本となる2つの計算規則(和の規則積の規則)を取り上げます。また、ベイズ学習に関わるややこしい用語たち(データ尤度関数事前分布事後分布エビデンス予測分布、などなど)に関しても念のためここで整理しておきたいと思います。そして最後に、簡単な多次元のガウス分布とウィシャート分布を使ったベイズ推論の例を取り上げ、それぞれの用語や概念との具体的な結びつきについて触れたいと思っています。

 

 

ベイズ学習の基本概念

さて、確率モデルを使ったベイズ推論を行う上で最小限必要なのは次のたった2つの計算ルールです。

<和の規則>

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

<積の規則>

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

{p(x,y)}同時分布(joint distribution){p(y|x)}条件付き分布(conditional distribution)と呼ぶんでした。極端な言い方をしてしまうと、ベイズ推論による機械学習のすべては、この2つのルールをいかに効率よく用いることによって未知の確率分布を計算するかということに尽きてしまいます。

さて、もう少し踏み込んでみましょう。まず最初に、上の2つの規則を組み合わせることによって次のような有名なベイズの定理(Bayes's theorem)を得ることができます。

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

{p(x|y)}で表される確率分布を、逆にした{p(y|x)}やその他の項の組み合わせで表現できてしまうというのがこの式の言わんとすることです。このことから、ベイズ推論は逆確率の法則なんて呼ばれたりもするみたいです。

 

ここでようやく機械学習らしい概念を導入してみることにします。といってもちょっと用語を導入してみるだけです。{x}の代わりに{\theta}を、{y}の代わりに{D}を文字として用いると、ベイズの定理は次のようになります。

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

 ここで{D}データ{\theta}パラメータと呼ぶことにします。

さて、データに関しては説明は要らないでしょう。写真データとか、人の声の録音データとか、月間の製品の売り上げデータとかなんでも想像しやすいものでOKです。データは複数そろっていることを基本的には想定しますが、実は1個でも0個でも大丈夫だったりします。

{\theta}は一般的にはパラメータと呼ばれます。 ベイズ学習においてパラメータは推論すべき未知の変数として扱われます*1

さて、改めてこの重要な式をじっくり見てみることにしましょう。まず{p(D|\theta)}尤度関数(likelihood function)と呼ばれています。パラメータがある特定の値に条件づけされたときに、データ{D}がどれだけモデルから発生しやすいかを表現しています。その後ろについている{p(\theta)}事前分布(prior distribution)と呼ばれており、{\theta}の分布に関する我々の事前の仮説(どんな{\theta}が出やすいのか?)を表現しています*2。それに対して一番左の式{p(\theta|D)}事後分布(posterior distribution)と呼ばれています。こちらも事前分布と同様パラメータの不確かさを表す確率分布ですが、データ{D}によって条件づけされているところが違いますね。この分布は尤度関数を通すことによって更新された{\theta}の分布というふうに解釈することができます。ベイズ学習の主要な課題はこの事後分布をいかに効率よく、精度高く計算するかにあります。

さて、尤度関数と事前分布の掛け算をしただけでは事後分布の計算は完了しません。分母にもう1つ項{p(D)}がついていますね。これはエビデンス(evidence)あるいは周辺尤度(marginal likelihood)呼ばれており、事後確率がちゃんと正規化される(=分布を積分すると1になる)ことを保証するために必要になってくる項です。エビデンスは和の規則を用いることによって次のように計算することができます。

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

「分子の項{p(D|\theta)p(\theta)}における{\theta}を和の規則で除去しなさい」と言っているんですね。簡単なガウス分布などを使ったモデルであればこの項は解析的に計算できるのですが、少し複雑なモデルになってくるとこの和が計算できなくなってしまいます。離散の場合は組み合わせ爆発を起こし、連続の場合には積分が解析不可能になってしまいます。ちなみにこのモデルエビデンスの値を使って複数のモデルの良し悪しを直接決めることもあるようです(モデル選択*3)。

 

さて、最後にもう1つ大事な概念を紹介してみましょう。機械学習では、ある学習データを使ってモデルを訓練させ、新しいデータに対して予測を行う、というのが代表的なアプリケーションになります。したがって、次のような新しいデータ{x}に関する予測分布(predictive distribution)を計算することも重要な課題の1つです。

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

事後分布(学習済みのパラメータの分布)を使って、まだ入ってきていない新データ{x}に関する確率分布(平均とか分散とか)を計算しようというわけです。

 

ここまでのベイズ学習における基本的な流れをまとめてみると次のようになります。*4

1、尤度関数(観測データの表現方法)を定義する。

2、事前分布(パラメータの散らばり具合)を設定する。

3、ベイズの定理を使って事後分布を更新する。(=学習)

4、事後分布から予測分布を計算し、未知の値の確率分布を求める。(=予測)

 

 

ガウス分布とウィシャート分布を使った具体例

ちょっと抽象的な話が続いてしまったので、次は具体的なモデルの例を当てはめて話を進めたいと思います。今回は、D次元、N個の観測データを{X}とおくことにします。尤度関数と事前分布は具体的には次のように定義してみたいと思います。

<尤度関数(ガウス分布)>

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

<事前分布(ガウス・ウィシャート分布)>

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

ガウス分布(Gauss distribution)はデータに対するノイズを表現する上で最も一般的に使われているモデルの1つです。ガウス分布には2つのパラメータ{\mu}(平均)と{\Lambda}(精度行列)がありますが、ベイズ学習ではこれらは確率変数として扱われます。したがって確率分布が必要になってくるのですが、今回は計算が解析的になるという理由で共役事前分布(conjugate prior distribution)ガウス・ウィシャート分布(Gauss-Wishart distribution)を用いることにします。ちなみに{\beta_0}{m_0}{\nu_0}{W_0}といった文字は超パラメータ(hyper parameter)と呼ばれており、事前分布の形状を決定するものになります。これらは確率変数として扱われず、固定された値を持ちます。

 

さて、ここまでが準備段階です。次にこのモデルを使ってデータを「学習」してみたいと思います。もっと丁寧に言うと、ベイズ学習における「学習」とは、データ{D}を観測した後の、パラメータの事後分布{p(\theta|D)}を計算することに対応します。今は非常に基本的なガウスモデルを使っているので、特に最適化やサンプリングなど凝ったテクニックは必要とせずに解析的に手計算することができ、得られる分布は事前分布と同じ形状(ガウス・ウィシャート分布)を持ったものになります。

具体的には、次のようにベイズの定理の分子だけを使って{\mu}{\Lambda}の関数を求めてあげると、いくらか計算が楽です。

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

得られた式を正規化してあげることにより、ガウス・ウィシャート分布が得られれます。

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

分布の形状を決めるそれぞれの超パラメータは事後分布では次のような計算結果になります。

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

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

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

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

計算過程は省きましたが、このような結果が自分ですらすら計算できるようになるとベイズ学習の入門コース卒業になります。

 

それでは次に予測分布も求めてみましょう。一般的にはパラメータの分布を学習しただけで満足するケースはあまりなく、学習後のモデルを使って何らかの予測を行うことが機械学習における主要なテーマになってきます。

さて、予測分布の計算は未観測のデータ{x}に対する尤度関数と事後分布を掛け合わせ、パラメータだけを積分除去すれば得られるんでした。さっそく計算してみましょう。

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

ガウス分布の平均パラメータと精度パラメータを両方とも積分除去しなければならないので若干計算がめんどくさいですが、結果は解析的に計算でき、スチューデントのt分布(Student's t-distribution)が得られます。それぞれの超パラメータは次のようになります。

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

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

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

 

・まとめ

今回はベイズ推論を使った機械学習の基本的な流れに関して説明しました。単純な線形回帰モデルなんかは今日の結果を利用すれば比較的簡単に構築することができます。一方で線形識別(ロジスティック回帰など)は前述の積分計算が解析的に実行することができず、ラプラス近似などの近似手法を使う必要があります。またクラスタリング等で使われる混合モデルや時系列データを取り扱うためのHMMなんかも、一般的には解析計算ができないので変分近似やサンプリングなどの近似手法に頼ることになります。ただこれらは単に計算上の問題なので、基本的には上述のベイズの定理を使った学習の枠組みを使って事後分布を計算するという流れは変わりません。

今回の記事がよくわからん!という方には,次のような入門書もあります.

books.rakuten.co.jp

*1: ”パラメータ”という用語はなんだか固定された値のように聞こえてしまうことがあります。最近では今回の記事のような{\theta}を(すべてのデータ点に共有される確率変数という意味で)”大域確率変数(global variable)”と呼んだりすることもあるようです。

*2:事前分布は基本的にデータを観測する前に決めてあげなければいけません。具体的にどのような事前分布を選んだらいいのかは、また別の機会でお話しします。

*3:ちなみに究極のベイジアンはモデル選択すらしません。仮説(あるいはモデル)Hが複数個あった場合、P(X) = ΣP(X|H)P(H)のようにすべての仮説を周辺化してデータを表現するからです。この場合、P(H)はそれぞれの仮説に対する信頼度のような事前知識を表します。

*4:ちなみに機械学習界隈では、訓練データをモデルに食わせることを「学習」と呼び、未観測の値に対して何らかの予測推定を行うことを「推論」と呼ぶこともあるようでうす。しかし、ベイズ学習においてはパラメータの学習にしろ未来値の予測にしろ欠損値の補間にしろ、すべて「未観測値に対する条件付き分布の推論」で言葉は片付いてしまいます。

ややこしい離散分布に関するまとめ

今回は離散分布(discrete distribution)の代表格である多項分布(multinomial distribution)や、その共役事前分布であるディリクレ分布(Dirichlet distribution)との関係性や計算方法を整理したいと思います。

離散分布というと、本来はポアソン分布(Poisson distribution)なども含めた離散値を出力するような分布全般のこと指します。しかし実際に論文などを読んでいると、くじ引きのように単純に出目の比率が与えられたような分布を離散分布と名付けてしまっている場合もよく見られます。まぁ文脈的に誤解を招くことはあまりないと思うのですが、くじ引きの分布をもっとキッチリ表現するなら、複数あるカテゴリーから1つを抽出するという意味でカテゴリカル分布(categorical distribution)と呼ぶのが適切かと思います。あるいは、観測回数が1の場合の多項分布とかって呼ぶこともできますね。

とまぁ別に名前はどうでもいいのですが、これらの離散分布の関係性をまとめてみると次の図のようになります。

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

ピンクの線が多次元変数への一般化で、オレンジの線が複数回の試行を行った場合に対する一般化です。グレーの点線矢印は、それぞれの確率分布に対する共役の関係性を示しています(矢印の元に当たるのが共役事前分布(conjugate prior)です)。 

結論としては、多項分布とその共役事前分布であるディリクレ分布の計算だけある程度馴染んでおけば、後はその特殊な例なのでいちいち気にする必要はないということになりますね。

 

 

<各確率分布の自己紹介>

さて、これからいろいろ計算確認をする前に、それぞれの確率分布の定義式を列挙してみます。

 

・ベルヌーイ分布(Bernoulli)

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

お馴染みのひしゃげたコインの分布です。{x}は0か1のどちらかの値を必ず取り、パラメータは{\mu \in [0,1]} をみたします。*1

・二項分布(Binomial)

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

コイントス{N}回行った時の、表面の回数に関する分布です。したがって{m}は0以上{N}以下の整数値を取ります。ベルヌーイ分布と同様、{\mu \in [0,1]} です。

・カテゴリカル分布(Categorical)

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

こちらもお馴染みの{K}面を持つひしゃげたサイコロの分布ですね。{K}次元ベクトルである{\bf{x}}の要素{x_k}は0か1のどちらかの値を必ず取り、いずれか1つの{x_k}しか1にならないので、{\sum_{k=1}^{K} x_k = 1} をみたします。それぞれの出目の確率は{\mu}で与えられ、{\sum_{k=1}^{K}\mu_k = 1}をみたすように設定します。

・多項分布(Multinomial)

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

{K}面を持つひしゃげたサイコロを{N}回投げた時の、出目の回数に関する分布です。したがって{K}次元ベクトルである{\bf{m}}の要素{m_k} は非負の整数値を必ず取り、{\sum_{k=1}^{K} m_k = N} をみたします。カテゴリカル分布と同様{\sum_{k=1}^{K}\mu_k = 1}です。

・ベータ分布(Beta) 

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

コイントスの共役事前分布に当たる分布です。コインの「ひしゃげ具合」に関する分布になります。この分布から出てくる値{\mu}は過ならず0と1の間に収まってくれるように作られています。分布の特徴を決める{a}{b}は、0より大きい実数値を設定する必要があります。ちなみに正規化項に当たる部分にはイカつい格好のガンマ関数(gamma function){\Gamma()}の塊がいます。ガンマ関数は階乗「!」の実数値バージョンですが、これはライブラリとか使えば簡単に計算できるのであまり怖がる必要はないです。*2

・ディリクレ分布(Dirichlet)

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

サイコロの共役事前分布に当たる分布です。サイコロの「ひしゃげ具合」に関する分布ですが、混合モデル(mixture model)に使われるほか、自然言語処理で使われるLDA(Latent Dirichlet Allocation)の由来にもなっていますね。出てくる値は必ず{\sum_{k=1}^{K}\mu_k = 1}をみたしてくれるように分布が作られています。また、すべての{k}に対して{\alpha_k}は0より大きい実数値を設定する必要があります。ここでもまたイカつい{\Gamma()}の塊が先頭にありますが、多くの場合で無視して計算できるので気にしなくても大丈夫です。

 

 

<各確率分布の関係性> 

さて、自己紹介が終わったところでそれぞれの分布の関係性を式で確認してみましょう。ここでは、多項分布やディリクレ分布がその他の分布の特別な場合であることを確認してみます。

 

・多項分布=>カテゴリカル分布

多項分布において>{N=1}と設定してあげれば、

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

となり、カテゴリカル分布に一致します。{N=1}と設定したことにより{m_k}は0か1かしか取らなくなったことと、{1!=0!=1}の計算に注意してください。

・多項分布=>二項分布

多項分布において次元を{K=2}と設定してあげれば、

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

となり、二項分布になります。{m_1+m_2=N}の制約があるので、実質{m_1}のみの確率分布になるのがポイントです。パラメータもわざわざ2個持つ必要はなく、{\mu_1+\mu_2=1}を考慮すれば{\mu_1}だけ持っておけば大丈夫ですね。

・ディリクレ分布=>ベータ分布

ディリクレ分布において次元を{K=2}と設定してあげれば、

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

となり、ベータ分布になります。ここでも、{\mu_1+\mu_2=1}の制約条件を利用することによって実質{\mu_1}のみの確率分布にすることができます。

 

同様に、ベルヌーイ分布はカテゴリカル分布から{K=2}とおくか、二項分布から{N=1}とおくと導けます。簡単なので省略します。

 

 

<共役事前分布とベイズ推論>

さて最後に、ディリクレ分布が多項分布の共役事前分布であることを確認してみましょう。共役事前分布とは、観測モデル{p(m|\mu)}とかけ合わせた後に同じ関数形状が出てくるような分布のことです。つまり、ベイズの定理

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

を適用したときに、事前分布{p(\mu)}と事後分布{p(\mu|m)}が確率変数{\mu}に関して同じ形式の関数になるように{p(\mu)}を定めるということですね。このような事前分布を選んでおくと、各種推論計算が解析的に簡単に計算できるようになるほか、データを小分けにして学習させるオンライン学習(online learning)を構築することも容易になります。オンライン学習は逐次学習(sequential learning)追加学習(incremental learning)とも呼ばれるようです。

ということで、ディリクレ分布を事前分布とし、多項分布に従うデータ{m}を観測した後の事後分布を計算してみましょう。

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

正規化項(ガンマ関数とか階乗とか含んでる項)を無視して計算するのがラクする秘訣です。さて、結果の式を見てみると、ディリクレ分布と同じ関数形状をしていることがわかります。したがって{\hat{\alpha}_k=\alpha_k+m_k}とおいてあげれば、

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

ということになり、事後分布もディリクレ分布になることがわかりました。*3

また、上のベイズ推論で{N=1}と置けば、カテゴリカル分布のパラメータに関する学習も同様にディリクレ分布を使って行うことができることがわかります。これは特に混合モデルクラスタリングなど)を扱うときに出てくる計算なので、ある程度慣れておくと便利です。

今回の記事がよくわからん!という方には,次のような入門書もあります.

books.rakuten.co.jp

*1:ちなみに英語で「べるぬうい」と言っても通じません。「バヌーリ」という方が近いです。

*2:ガンマ関数はプログラム上実装することはそれほど多くはないのですが、どうしても計算する必要がある場合は、対数を返してくれる関数(lgammaなど)を使うと良いです。通常はものすごく大きい値になってしまうので。

*3:「事後分布の計算結果がカテゴリカル分布にも見えるんだけど・・・」という方は、ちょっと落ち着いてみてください。今は{m}ではなく{\mu}に関する分布をベイズの定理を使って計算しています。したがってこの関数形状はディリクレ分布一択になります。

グラフィカルモデルを使いこなす!~転移学習を表現してみる~

さて、今日は以前ご紹介したグラフィカルモデルを使って、転移学習(Transfer Learning)の一例をモデル化してみたいと思います。この記事を読んでいただければグラフィカルモデルを使ったベイズ学習が、機械学習における様々な問題設定に対して柔軟なアプローチを与えてくれることがわかっていただけるかと思います。

 

[必要な知識]

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

 

・グラフィカルモデルによる表現

今日は転移学習の一例をグラフィカルモデルを使って表現し、さらにグラフ上での推論を考えてみたいと思います。最後に、具体的なガウス分布を使った例を用いて簡単な実験をしてみたいと思います。

自分は転移学習に関して取り立てて何かやった経験がないので正直なところあまり深淵な議論はできないのですが、ここではとりあえず

「ある確率変数を推定する課題(ドメイン)が2つあり、それぞれに共通する情報と共通でない情報が存在するような問題設定」

という程度で勘弁してください。そして目標は、2つのそれぞれのドメインを独立に扱うのではなく、片方のドメインのデータを上手に利用して残りのドメインの推定精度を向上させることになります。

 

というわけで、早速ですが次のようなモデルを考えてみました。

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

ちょっとややこしいですね。ここでは2つの異なるドメインを取り扱うことにします。一つずつ変数を説明していくと、{x_1}{x_2}はそれぞれのドメインにおける観測データです。それに対して、星のついている{x_1^*}{x_2^*}はこれから予測したい未観測データを表します。{\theta_1}{\theta_2}はそれぞれのドメイン固有のパラメータで、{\theta_s}は両方のドメインでシェアされているようなパラメータです。例えば「平均値はそれぞれのドメインで違うんだけど、分散は共有されている」みたいなモデルを想像してみてください。*1

データが観測されていない状態でのモデルを書き下してみると、

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

のようになります。ちょっと長くて申し訳ないですが、上のグラフィカルモデルとしっかり対応が取れているか(黒丸を白丸にして)確認してみてください。

 

・推論

さて、次に観測されていない変数たちの事後分布を計算してみることにしましょう。

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

長いですね。一個ずつちゃんと確認すれば難しくはありません。一番最初の右辺は、ただ単に確率のproduct ruleに従って確率変数たちをバラバラに分解してあげただけです。次の式からは、このモデルに対して以前にご紹介した有向分離を適用して変数の独立性を一つ一つ確認し、独立な変数を赤い斜線で消していったものになります。

さて、最後の式を一つ一つじっくり見ていってみましょう。まず一番後ろの{p(\theta_s|x_1,x_2)}ですが、これはドメイン間共有のパラメータである{\theta_s}の分布が両ドメインのデータ{x_1,x_2}から学習できることがわかりますね。次にドメイン{k}に対するパラメータの分布{p(\theta_k|\theta_s, x_k)}は、どちらのドメインでも共有パラメータ{\theta_s}が必要ですが、相手ドメインのデータに関しては興味がないと言っています。最終的な未知変数{x^*_k}の予測は共有パラメータとドメイン固有のパラメータに依存して決まるようですね。

というわけで、上のモデルでは共有パラメータ{\theta_s}が両ドメイン間を仲介することにより、それぞれのドメインで予測したい{x^*_k}の分布がお互いの学習データに依存しあうことがわかりました。

ちなみに、多くの場合では我々はパラメータの確率分布には関心がありません。最終的に求めたいのは、未観測データである{x_1^*}やの予測分布です。この予測分布を式で表現すると次のようになります。

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

確率のsum ruleを用いて、関心のないパラメータの分布を積分除去(integrate outまたはmarginalize out)してあげただけですね。このようにして求めた分布は周辺分布(marginal distribution)とかって呼ばれたりしています。次の実験ではこの積分操作が解析的に実行できるような確率分布を選んで使うことにします。

 

・実験

さて、上の結果に対して具体的な確率モデルを当てはめて実験してみましょう。今回は画像で見やすいように2次元のガウス分布のモデルを考えます。ドメインが2つあり、それぞれで平均値は違うんだけど、分散パラメータは共通の分布からサンプルされている、というような設定です。ドメイン1における学習データを増やすことにより、ドメイン2の学習結果(分散の値)も変化していく様子を表現してみたいと思っています。できれば具体的な式展開を示してご説明したいのですが、ものすごく長くなってしまうので今回は割愛します。*2

さて、早速ですが実験結果のアニメーションです。

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

ご覧いただけるように、左側のドメイン1ではデータ数{N_1=1}から順に{N_1=50}まで与えていってあげています。右側のドメイン2はデータ数{N_2=5}で固定したままです。データは左側にしか与えていないにもかかわらず、右のドメインでは分散の形状(楕円)が同じように学習されていっていることに注目してください。

ちなみに、非常に見づらくて申し訳ないのですが、分散の大きさは左右で同じではないです({N_1=1}のときにかろうじてわかります。)。これは、もともとの分散行列に対する事前分布は同じ設定なのですが、事後分布はデータの与え方によって当然異なってくるためです。データ数が少ないときは左側の分散の方が大きく、データ数が多くなってくると左側の分散の方が小さくなっていきます。

また右側のドメインでは平均値は変化しません。なぜなら、左側で与えられているデータたちは右側のドメインの平均値に対する情報を与えないからです(まぁ、そういう風にモデル化してるからなんですが)。

 

・まとめ

というわけで、今回は転移学習の一例をグラフィカルモデルで表現してみました。今回のようなモデル以外にも例えば「ドメイン1とドメイン2でパラメータは同じだが、データがある関数によって{x_2=f(x_1)}のように変換されている」みたいなモデルを考えることもできます。この場合では、データを生成するパラメータたちと、変換関数のパラメータをデータから学習することになります。モノクロ画像のデータをRGB画像のデータと組み合わせるとか、アングルの異なるセンサーから得たデータを対応付けるとかに使えそうですね。

いずれにしても重要な結論としては、確率モデルを使ったベイズ学習ではいちいち「転移学習」といった固有の課題に対する問題意識は持ちません。実際のデータ解析の現場において、異なるデータ間で何かしら共通する情報があるというのであれば、ただ単に「それをモデル化すればいい」というだけの話になります。

データやデータの裏にある背景を入念に調べ上げ、丁寧にモデル化して推論をする―というのがベイズ推論を使った機械学習の本質です。

*1:うーむ、どんな状況でしょう。例えば世界の都市の人口密度なんかは、平均値(緯度と経度)こそバラバラに散らばっていますが、分散(人口の広がり具合)には共通の傾向があるかもしれないですね。もちろんそこまで単純ではないですが。

*2:今回の実験では、データ観測モデルに対してガウス分布を、平均値と分散のパラメータに対してガウス・ウィシャート分布を使用しています。また、最終的にプロットしている予測分布は多次元のt分布になります。