作って遊ぶ機械学習。

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

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

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

混合モデルの代表的なアプリケーションはクラスタリングですが、今回ご紹介するモデルの構築方法はクラスタリングにとどまらず、もっと凝ったモデルを構築するための基本パーツになるものです。例えば時系列情報を取り扱うための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分布になります。

アニメでわかるベイズ推論によるパラメータ学習

さて、今日はガウス分布を使った簡単な実験を行って、ベイズ推論における機械学習の本質の一端を説明したいと思います。せっかくなので前回取り扱った多峰性事前分布も実験に取り入れてみたいと思います。

 

改めてベイズ学習を数式で書くと次のようになります。

{p(\theta | x) \propto p(x|\theta)p(\theta) }

パラメータに関する事前の知識{p(\theta)}が、尤度関数{p(x|\theta)}を通して、事後の知識{p(\theta|x)}に変換されるんでしたね。今回はこのプロセスをアニメーションを見ながら確認してみようというお話です。

 

で、今回は次のような平均値パラメータ{\mu=0}を持った真のガウス分布のパラメータを推論する問題を考えてみたいと思います。

{N(x|0, \Sigma)}

分散{\Sigma}は簡単化のため、既知で固定ということにしておきます。先ほどのベイズ学習の表記を使うと{\theta = \{\mu \}}ということになりますね。

 

さて、データ{x}に対する観測モデル{p(x|\mu)}ガウス分布を使うとして、事前分布{p(\mu)}は違ったものを3種類用意してみたいと思います。そして実際に真のガウス分布からサンプルされたデータを与えてあげることによって、それぞれの事前分布たちがどのような事後分布に成長していくのかを観察してみたいと思います。

 

<共役事前分布A>なだらかな分布を持つ、あまり事前情報のない分布

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

 

<共役事前分布B>ちょっと偏った思想を持つガンコな分布 

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

 

<多峰事前分布C>多峰性を持つちょっと優柔不断な分布

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

 

さて、これらの個性的な登場人物に対して、真の平均値0を持つ未知のガウス分布からサンプルされたデータ点を1つずつ与えてあげるとどうなるのでしょうか?アニメを作ったのでご覧ください。

 

・共役事前分布A

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

・共役事前分布B

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

・多峰性事前分布C

 

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

 

左の図には再びそれぞれの事前分布の形状を表示しています。右の図は観測されたデータを元に計算されたパラメータの事後分布で、ピコピコ現れる縦線は観測されたデータを表しています。得られた分布はデータ点に対する密度分布{p(x)}ではなく、平均値パラメータに関する不確かさ{p(\theta|x)}を表しているので注意してください。細くて0に集中している分布ほど、平均値パラメータが0であることの確信の高さを表しているんですね。

 

さて、改めてデータ数がある程度十分に観測されたとき(40点)の事後分布を比べてみましょう。

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

興味深いのは、それぞれが違う事前分布から学習をスタートしたにもかかわらず、ある程度学習データが集まるとみんな同じような事後分布に収束していることです。ただ、ガンコ者の事前分布Bは未だに最初の知識をほんの少し引きずっているように見えます(平均値が左にずれてますね)。優柔不断だった事前分布Cは、多峰性の形状がほとんど失われ一つのガウス分布に収まっているように見えますが、それでも事前分布Aと比べると若干あいまいさを残した広めの結果になっていますね。

この結果はベイズ学習において、もっとも重要な特性を示しています。データが少ないうちはそれぞれの推論結果は事前分布の形状に大きく左右されます。しかし証拠となるデータがたくさん集まるにつれて、最初の知識による影響は次第に薄れ、すべての事前分布たちが似たような形状の事後分布に収束していきます。

ベイズ推論が「異なった意見を持つ人たちでも証拠さえ十分に集まれば最終的に合意を得る」という、非常に理にかなった原理に基づいていることが分かります。*1

 

さて、この結果から機械学習におけるベイズ学習の利点をいくつか見ることができます。

第1に、ベイズ学習ではデータの数が少ない(最悪の場合0個!)の場合でも何かしらの予測結果をひねり出すことができます。データが少ないなら少ないなりに、事前知識を駆使して何かしらの結果を出してくれるのがベイズ学習の強みです。

第2に、ベイズ学習は複雑なモデルに対する学習においても力を発揮することができます。実際の機械学習の問題では、上記のような単純なガウス分布の平均値を推定するだけの簡単な問題はありません。混合ガウス分布隠れマルコフモデルに代表されるような、グラフィカルモデルを駆使したはるかに複雑なモデルを使って問題を解くのが通常です。このようなモデルでは学習に必要なデータ数が多くなるのが通常で、「ビッグデータ」といえどもまるで学習データが足りない、というような状況はよく起こります。*2

第3に、ベイズの定理を使ったこのような事前分布から事後分布への学習プロセスは、非常に自然で効率的な逐次学習の手段になります。これは、例えばリアルタイムで学習・予測するようなオンラインシステムを開発する場合や、他のモデルで学習した結果を別のモデルに組み合わせたりする場合に非常に便利な特性です。

 

以上、今回は簡単な実験でベイズ学習の本質に迫ってみました。細かい技術面に注目するとベイズ学習の利点は実はこれだけではありません。また別の機会でベイズ学習におけるモデル選択や予測分布に関するお話をしたいと思います。

*1:ちなみに事前分布で確率0を割り当ててしまった領域に対しては、いくら学習しても合意に至りません。今回のようなガウス分布の例では実数すべての値に対して0ではない確率を割り当てているのでそういったことは起こりませんが、実際のモデルでは、事前のモデルで表現しえないような予測結果は決して得ることができないという事実は、当たり前ですが知っておいた方がいいでしょう。

*2:ビッグデータの本質は、ただ単にデータのレコード数が増えるということではなく、データの次元数が増加することにあると個人的には考えています。予測精度を高めるためには様々なリソースから得られたデータを統合して学習に使う必要がありますが、このような多次元データに対しては、最尤推定や最適化などの手法は簡単にオーバーフィットしてしまうため、ベイズ学習が有効な解決手段になります。

第一線のAI研究者が注目する最新機械学習技術6選(NIPS2015招待講演より)

ちょっと前になりますが、昨年12月に行われた機械学習のトップカンファレンスであるNIPS2015の講演ビデオが上がっているようなのでチェックしてみました。今回ご紹介するのはケンブリッジ大学のZoubin Ghahramani教授の研究です。

datasciencereport.com

同教授は今後excitingな機械学習の基礎・応用に関する取り組みとして次のような6テーマを紹介しています。

 

・Bayesian Nonparametrics

ベイジアンノンパラメトリクスでは、無限次元を持つモデルを仮定することにより、データ量に応じて適切なモデルを学習することができます。こういったモデルは関数上の確率分布(確率過程)を考慮することによって実現できます。代表的な例はガウス過程、中華料理店過程、インド料理過程などです(ご飯ばっかりですね)。

ガウス過程は、回帰や識別、ランキングや次元削減に使われています。またある種のニューラルネットワークガウス過程の一例であることが1994年に示されているそうです。

・Probabilistic Programming

確率モデルを構築し、推論アルゴリズムを導き、実装をする・・・といったプロセスはものすごく時間がかかって大変です。Probabilistic programmingは確率モデル(あるいは確率的にデータを生成するプロセス)を記述するための言語で、推論はサンプリングや変分近似で勝手にやってくれます。代表的なのはBugsやStanです。

・Bayesian Optimization

日本語だと「ベイズ的最適化」の名前でコンセンサスが取れていたと思います。ガウス過程を使い未知の関数の最適値を探す手法です。未知の関数の評価にとても時間orコストがかかるような状況を仮定しており、いかに少ない評価回数で最適値を探せるかがポイントになります。

・Data Compression

シャノンの定理によると、すべての圧縮アルゴリズムは確率モデルに基づくそうです。ベイジアンノンパラメトリクスを使ったアダプティブな圧縮アルゴリズムが近年非常によい圧縮性能を出しているようです。

・Automatic Statistician

データを食わせるだけで何かしらレポートを返してくれるような野心的なアルゴリズムの話です。アイデアは非常にシンプルで、いくつかのカーネルを組み合わせて周辺尤度を評価することによって、データをうまく説明するような適切なモデルを探索します。

・Rational Allocation of Computational Resources

大量のデータと複数のモデルがある中、リソース(時間、CPU、メモリ、ディスクサイズ)が限られた状況下でいかに最適な機械学習システムを構築するか、という課題設定です。これを不確実性がある中での逐次的な意思決定問題として扱います。

 

・個人的な雑感

Bayesian nonparametricsは従来のベイズ学習に対する正当な拡張であるため、今後も基礎的な研究は発展していくかと思います。そろそろ大きなアプリケーションがほしいところですが、ディープラーニングなどの成功しているモデルもかなりの部分はこれに置き換わり、より洗練されたアルゴリズムになるのではないかと予測しています。

Probabilistic programmingは基本的に推論性能はあまりよくないようです。僕は使ったことがないのですが、とりあえず複数のモデルを素早く試したいときなんかには便利だと思います。普通にモデル構築→アルゴリズム導出→実装ってやってるとそれだけで1日使いますからね・・・。

ベイズ的最適化は非常に概念もわかりやすく応用先も広いのでぜひ多くの人に使っていただきたいと思います。日本の京などのでっかいコンピュータで何かシミュレーションを行う場合は、Bayesian optimizationで全体を包んであげると解探索が早くなるかもしれません。あまりいないと思いますが、例えば遺伝的アルゴリズムとかを使っている人は早々にこちらへの置き換えを検討した方が良いと思います。

データ圧縮も面白い領域です。結局のところ、jpegに代表されるように、非可逆圧縮というのはデータの中の捨ててもいい情報と残さなきゃいけない情報をうまく見つけることが基本になっています。これは機械学習における特徴量抽出といっしょです。画像だったら、人間が認知に必要でない情報は捨ててもいい可能性が高いと言えるわけですね。データ圧縮もかなりドメイン依存が強いかと思いますが、汎用的なアルゴリズムでどこまでいけるのかは興味深いですね。 (追記:すみません、招待講演の内容では「非可逆圧縮」ではなく、情報損失のない「可逆圧縮」の話だったようです。データに関する分布をベイジアンノンパラメトリクスを使って推定することにより、データを表現するための符号長を統計的に短くすることができるという話です。研究のポイントとなる点は、zipなどで使われている性能の良い可逆圧縮アルゴリズムと同等のものがベイジアンノンパラメトリクスに基づいたアルゴリズムで解釈できることを示し、かつ性能も向上させたことです。すごいですね。)

Automatic statistiicianはデータサイエンスにおける究極的な課題の1つなんじゃないかと思います。ただ、質問でもありましたが、決してデータサイエンティストの職業を奪うようなものではなく、あくまで便利なツールであるという位置づけだそうです。

リソースアロケーションの話は実応用で非常に重要です。まぁそれ以前に日本だったらちゃんと複数の確率モデルを作れる人がいないと話にならないですが・・・。

 

以上です。

多峰性の事前分布ってどうやって作るの?

ベイズ学習では、複雑なモデルにおけるパラメータの学習を効率的に行うために、しばしば観測モデルに対する共役な事前分布を仮定します。例えばベルヌーイ分布のパラメータの事後分布を推定するために、事前分布をベータ分布にしたりします。ウェブを検索すると様々な確率分布のパラメータに対応する共役事前分布のリストを見つけることができますが、そもそも、データを表現する観測モデルが決まってしまうと自動的にそのパラメータの事前分布が決まってしまうのは何だか窮屈な感じがしてしまいます。

今回はグラフィカルモデルを工夫することによって、もっと自由な、例えば多峰性の事前分布などを簡単に作れることをお話ししたいと思います。

 

[必要な知識]

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

 

さて、改めてベイズ学習におけるパラメータの学習を書いてみると次のようになります。

{ p(\theta| x) \propto p(x|\theta)p(\theta) }

ここで{x}は観測データ、{\theta}はモデルのパラメータです。さらに{p(\theta|x)}はパラメータの事後分布{p(x|\theta)}尤度関数{p(\theta)}はパラメータの事前分布と呼ばれています。ベイズ学習を言葉に翻訳すると、

「パラメータに関する我々の事前知識{p(\theta)}は、観測データに対するモデルの尤もらしさ{p(x|\theta)}を通して、事後の知識{p(\theta|x)}に変換される」

というプロセスになります。

で、今回は尤度関数を計算するモデルとして次のようにガウス分布を仮定します。

{ p(x| \theta) = \mathcal{N}(x| \mu, \Sigma_x ) }

ここで今回推定したいパラメータは平均値{\theta=\mu}であり、分散{\Sigma_x}に関しては既知で固定であるものとします。

 

・単峰性事前分布によるガウス分布の推定

まず始めに、ガウス分布の平均値パラメータに対する一番シンプルな事前分布の選び方は、共役事前分布であるガウス分布を使うことです。*1

{ p(\mu) = \mathcal{N}(\mu | m, \Sigma_{\mu}) }

ここで{m}{\Sigma_{\mu}}はそれぞれ平均値と分散パラメータであり、{\Sigma_x}と同様固定値であるものとします。

このモデルに対して{N}個のデータを観測した後の事後分布は、グラフィカルモデルで表現すると次のようになります。

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

さて、さっそく事後分布を計算してみましょう。事後分布はベイズの定理から次のように計算できます。

{ p(\mu|x) \propto p(x|\mu)p(\mu) }

両辺対数を取って計算を進めることにします。{\mu}に関する分布を求めたいので、無関係な項はすべて定数項{c}に吸収させてしまいます。

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

結果はやはりガウス分布になり、新たに文字{\hat{m}}{\hat{\Sigma}_{\mu}}を導入すれば次のように書けます。

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

・多峰性事前分布によるガウス分布の推定

さて、今度は先ほどの単純な共役事前分布と違い、{\mu}に対して{K}個の異なる平均値の分布からなる、多峰性の事前分布を構築してみましょう。感覚的な例をいうと、推理小説で「犯人が一人なのはわかっているけど、容疑者が複数いる」というような状況設定です。

これは次のグラフィカルモデルで示すように、新たな潜在変数{z}を導入することによって実現できます。

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

まず始めに新しく登場した{z}に対する確率分布を定義してみましょう。ここではもっとも単純な離散変数の分布であるカテゴリカル分布を使います。

{\begin{eqnarray} p(z) = Cat (z|\pi) = \prod_{k=1}^{K} \{ \pi^{(k)} \}^{z^{(k)}} \end{eqnarray} }

ただし{\sum_{k=1}^{K} \pi^{(k)}=1 }です。さらに、{\mu}の分布は次のように{z}の条件付き分布になります。

{ \begin{eqnarray} p(\mu|z) = \prod_{k=1}^{K} {\mathcal{N}(\mu| m^{(k)}, \Sigma_{\mu}^{(k)})} ^{z^{(k)}} \end{eqnarray} }

このような{z}に対する記法は1 of K表記というように呼ばれたりします。例えば{z=(0,0,1,0)}のように、{z}{K}個ある要素のうち一つだけが1を取ることによって、どの分布が使われるのかを指示(indicate)することができます。

それでは事後分布を計算してみることにしましょう。

{p(\mu, z| x) \propto p(x,\mu,z) = p(x|\mu)p(\mu|z)p(z) }

今回は{\mu}{z}の2つの事後分布を求めなければなりませんが、これはそれぞれ別々に解析的に計算できます。まず上の数式から{\mu}に関わる項のみを取り出してみましょう。例によって対数計算を行うと、

{ \ln p(x,\mu,z) = \ln p(x|\mu) + \ln p(\mu|z) + c }

となります。先ほどの単峰性の例と同じように{\mu}のみに注目して計算を行うと、結果として{\mu}の事後分布は{z}に依存する形となり、次のような(条件付き)ガウス分布になります。

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

さて、今度は{z}の事後分布も求めてみましょう。これは{\ln p(z|x) = \ln p(\mu|z) - \ln p(\mu|z,x)}を考慮すれば簡単に計算できます。

{ \begin{eqnarray} \ln p(z|x) = \sum_{k=1}^{K} z^{(k)} (-\frac{1}{2} {\mu^{(k)}}' \Sigma_x^{-1} \mu^{(k)} -\frac{1}{2} \ln |\Sigma_{\mu}^{(k)}| \\ +\frac{1}{2} {\hat{\mu}^{(k)}}' \hat{\Sigma}_x^{-1} \hat{\mu}^{(k)} +\frac{1}{2} \ln |\hat{\Sigma}_{\mu}^{(k)}| +\ln \pi^{(k)})+c \end{eqnarray} }

{z}に掛け算されている部分が事後分布のパラメータになります。書き直すと次のようになります。

{ \begin{eqnarray} p(z|x) = Cat (z|\hat{\pi}) = \prod_{k=1}^{K} \{ \hat{\pi}^{(k)} \}^{z^{(k)}} \end{eqnarray} }

{\hat{\pi}}は次のように計算されます。

{ \begin{eqnarray} \hat{\pi}^{(k)} \propto exp\{ -\frac{1}{2} {\mu^{(k)}}' \Sigma_x^{-1} \mu^{(k)} -\frac{1}{2} \ln |\Sigma_{\mu}^{(k)}| \\ + \frac{1}{2} {\hat{\mu}^{(k)}}' \hat{\Sigma}_x^{-1} \hat{\mu}^{(k)} +\frac{1}{2} \ln |\hat{\Sigma}_{\mu}^{(k)}| +\ln \pi^{(k)} \} \end{eqnarray} }

実装上は{\hat{\pi}}が足して1になるように正規化してください。

以上、事後分布{p(\mu,z|x)}を見てみると、事前分布{p(\mu,z)}と同じ形状を持っていることになります。ガウス分布の平均値だからといって単純にガウス事前分布を使う必要はなく、今回のような多峰性のガウス分布を事前分布にしても、最適化のような計算をせずに解析的に事前分布と同じ関数形の事後分布が求まるんですね。*2

 

さて、ちょっと思った以上に計算式ばかりの記事になってしまったので今日は一旦ここで切ります。今回はガウス分布の平均値パラメータの事後分布を推定する例を使って、多峰性の事前分布をどのように構築したらよいのかを説明しました。次回は今回紹介した事前分布を使った簡単な実験を行い、ベイズ学習における興味深い特性に関してお話したいと思います。

 

=>続きの記事はこちらです。

アニメでわかるベイズ推論によるパラメータ学習 - 作って遊ぶ機械学習。

 

*1:ガウス分布の平均値の共役事前分布は「たまたま」ガウス分布になります。ガウス分布を観測、平均値ともに使っているので気を付けてください。

*2:ちなみにこの多峰性事前分布によるモデルは、混合ガウス分布とは異なります。混合ガウス分布の場合は潜在変数がデータ{x_n}ごとに割り当てられます。今回の例と違い、解析的に求めることは不可能で、変分近似のような近似推論法を使う必要があります。