emceeを試してみた
pythonのMCMCライブラリとしてemceeというのがあるらしいので試してみました。
Paralell tempering(レプリカ交換モンテカルロ法)が使えるの他のライブラリとの大きな違いになります。
http://dan.iel.fm/emcee/current/
ほかにもMultiprocessing,MPIに対応している点などが特徴です。
本体はpure pythonで書かれているらしいのですが、普通のpythonと同様にCythonなどで高速化が出来るようです。
インストール
他のライブラリに依存していないので
pip install emcee
だけでインストールは完了です。ただしMPIを使う場合にはmpi4pyをインストールする必要があります。
実行例
http://dan.iel.fm/emcee/current/user/line/#maximum-likelihood-estimation
にある人工データを使った例です。
http://nbviewer.ipython.org/gist/xiangze/b75a9066162396b4a2a0
に実行結果を置きました。
事前分布,尤度のlogを関数として定義します。
事前分布は無情報分布です。複数のハイパーパラメータに対して一括で確率の値を返しています。
def lnprior(theta): m, b, lnf = theta if -5.0 < m < 0.5 and 0.0 < b < 10.0 and -10.0 < lnf < 1.0: return 0.0 return -np.inf
尤度はガウス分布なので
def lnlike(theta, x, y, yerr): m, b, lnf = theta model = m * x + b inv_sigma2 = 1.0/(yerr**2 + model**2*np.exp(2*lnf)) return -0.5*(np.sum((y-model)**2*inv_sigma2 - np.log(inv_sigma2)))
と書けます。
事後分布=事前分布*尤度+定数 (lnprob=lnlike+lnprior)のベイズの定理そのままです。
def lnprob(theta, x, y, yerr): return lnprior(theta) + lnlike(theta, x, y, yerr)
分布に関してはstan,pymcと比べると手作り感があります。
階層的なモデルのときにも同様にlnlikeを足しあわせていけば良さそうです。
定義したlnprobをEnsembleSamplerへ代入して、run_mcmcでサンプリングを行います。その際に推定すべき変数の数(ndim)とwalkerの数を指定する必要があります。walkerとはchainの数のようなものですが、chainごとに完全に独立しているわけではないようです(http://dan.iel.fm/emcee/current/user/faq/#what-are-walkers)。
初期値posとして例では scipy.optimizeで推定された値を指定しています。
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(x, y, yerr))
sampler.run_mcmc(pos, 500)
sampler.chainにサンプルされたデータがあります。
可視化
http://nbviewer.ipython.org/gist/xiangze/b75a9066162396b4a2a0#plot,可視化
例ではtriangle_plotというライブラリを使って分布を可視化しています。作者が同じ人らしいです。
Paralell tempering
emceeの特長であるParalell temperingを行うにはPTSamplerを使用します。例では混合ガウス分布を推定しています。
http://nbviewer.ipython.org/gist/xiangze/b75a9066162396b4a2a0#Paralell-tempering
ntempsで温度の数を指定します。EnsampleSamplerと違い事前分布と尤度を別々に指定するようです
(http://dan.iel.fm/emcee/current/api/#the-parallel-tempered-ensemble-sampler)。
walkers数 = 100
temps数 = 20
10000 stepで300秒くらいかかりました。
低温部を見ると多峰がサンプルできているようです。
http://nbviewer.ipython.org/gist/xiangze/b75a9066162396b4a2a0#低温部の結果
この例では
walkers数 = 20
temps数 = 20
でもサンプリングがうまくいっているようです。
