Expectation Propagation Methodの実装

はじめに

Expectation propagation method(EP法)はベイズ推論の近似解を求めるための手法で、2001にMinkaが提案したようです。PRMLで紹介されているほか、GPMLでもロジスティック回帰の良い解を得るために導入されています。

手法

データ$D = \{x_1, \ldots , x_N \}$とパラメータ$\theta$の同時分布を因子の積で表現します。

ここで、パラメータの事前分布と$N$個の尤度の積として表すと分かりやすいのではないでしょうか。

次に、事後分布

を以下のように近似因子の積で近似します。

EP法ではこの近似因子を一つずつ正確な尤度に近づくように$q(\theta)$を更新していきます。 まずサイト$i$の近似因子を最適化することを考えます。$q(\theta)$から近似因子を除いた分布は次式から得られます。

この分布(cavity distribution)と正確な尤度の積からなる分布

に近づくように$q(\theta)$を更新します。ここで、$Z_i$は正規化項であり、

を計算することで得られます。したがって、最適化は以下のKLダイバージェンスを最小化することです。

この最小化は$q(\theta)$が指数型分布族であるとき、一次および二次モーメントを一致させることで求められます。 上式の左側の分布の平均および分散を用いて近似事後分布を$q^{new}(\theta)$として更新できたら、サイト$i$の近似因子を更新します。

以上の手続きを$i=1,\ldots ,N$について繰り返し計算を行なっていきます。

実験

Minkaの論文やPRMLと同様に、clutter problemを例に実装しました。 2次元ベクトルの観測データが100個あるとき、本当に知りたいデータの平均パラメータをベイズ 推論から求めます。 問題設定の詳細は文献をあたってください。

下図に今回用意したデータを示します。青のプロットが知りたいデータの分布(平均を(1,1)とした)から生成された観測データ、赤がノイズデータです。

EP法では近似因子に対して初期値を設定する必要があります。次のように初期値を与えるのが普通でしょう。

学習の結果を下図に示します。青のプロットが各サイトの平均パラメータ、赤のxが事後分布の平均パラメータを表しています。 推論値は(0.96, 0.82)となりました。

実装コード

jupyter notebookをgithubに公開しました。計算の安定性を欠いているかもしれません。 http://nbviewer.jupyter.org/github/KDOG08/kdogBlog/blob/master/src/EPmethod.ipynb

参考文献

Minka, T. Expectation Propagation for Approximate Bayesian Inference. https://arxiv.org/abs/1301.2294

Minka, T. A family of algorithms for approximate bayesian inference.

Bishop, C. Pattern Recognition and Machine Learning.

Rasmussen, C. and Williams, C. Gaussian Processes for Machine Learning. http://www.gaussianprocess.org

Barthelmé, S. and Chopin, N. Expectation-Propagation for Likelihood-Free Inference. https://arxiv.org/abs/1107.5959