読者です 読者をやめる 読者になる 読者になる

xiangze's sparse blog

機械学習、ベイズ統計、コンピュータビジョンと関連する数学について

emceeを試してみた

pythonMCMCライブラリとして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
でもサンプリングがうまくいっているようです。