**黒木玄 Gen Kuroki** @genkuroki · 2017-07-12T19:16:21Z

黒木玄 Gen Kuroki @genkuroki

データにフィットしないことを承知でPoisson分布でMCMCしてみた。NUTS以外のサンプラーも試しています。

ううむ、Juliaのことをよく理解していないせいだと思いますが、色々奇妙なことがあってよくわからない。エラーがどうして出るのか？エラーがどうして消えるのか？リンク先を見れば何の話かわかるようになっています。

Jul 12, 2017, 19:16 · Web · 0 · 1

**黒木玄 Gen Kuroki** @genkuroki · Jul 13, 2017, 05:34

**黒木玄 Gen Kuroki** @genkuroki · Jul 13, 2017, 05:34

Julia Mamba による階層ベイズモデルでのMCMCの例も作りました。手抜きの方法で予測分布とWAICも計算してあります。

そこで扱っている階層ベイズモデルは階層ベイズでない負に二項分布によるベイズモデルに等価なので、そのことを利用して、手抜きの易しい方法(自前で積分の計算を行わない方法)で予測分布やWAICを計算してあります。

添画像の1つ目がサンプラーが生成した事後分布の散布部で、2つ目が予測分布とサンプルのヒストグラムを重ねたものです。色がシアンの線はサンプルから直接KDEで作ったPDF。

**黒木玄 Gen Kuroki** @genkuroki · Jul 13, 2017, 05:44

**黒木玄 Gen Kuroki** @genkuroki · Jul 13, 2017, 05:44

階層ベイズでMCMCしてみた感想

・等価な非階層モデルでのMCMCは3秒ほどで終わるのですが、階層ベイズの方は約1000秒かかりました。かなり重いです。(どちらもiterationsは1000)

・少々手抜きでも軽い方がうれしいと思いました。

・もっと速く計算する方法はないのでしょうか？

・それとも私のやり方に何か問題があるのでしょうか？

・まだJuliaのMambaパッケージの挙動をよく理解していません。

・JuliaならJulia内で Pkg.add("Mamba"); using Mamba に成功すれば、誰でも気軽にMCMCによるベイズ統計で遊べます。

・どしろうとの私でもすぐにできたくらいなので、きっと誰にでもできると思います。

・Juliaはカジュアルな使い方にも向いている実用的な言語。

・注意！各種パッケージが毎日にように更新されていて、安定していません！こういうのを楽しめる精神的余裕は必要かも。

**黒木玄 Gen Kuroki** @genkuroki · Jul 13, 2017, 09:36

**黒木玄 Gen Kuroki** @genkuroki · Jul 13, 2017, 09:36

階層ベイズのケースで NUTS サンプラーを使用したら、1000 iterations で 977 秒もかかったので、別のサンプラーも試してみました。

NUTS (no-u-turn sumpler)では、初期値はものすごくテキトーでいいし、パラメーターの調節もいらず、チェーンの質もすごくよいです。しかし、遅い。

そこで質は下がったり、初期値の範囲やパラメーターの調節だが、速いサンプラーを色々試してみました。

アルゴリズムがノータイムで理解できるほど分かり易いスライスサンプラーなら、10000x3が60秒で終わりました。これならイライラせずにすみます。

**黒木玄 Gen Kuroki** @genkuroki · Jul 13, 2017, 14:41

**黒木玄 Gen Kuroki** @genkuroki · Jul 13, 2017, 14:41

階層ベイズの確率モデルの記述部分の無駄を省いたら少し速くなった。

ベクトル値の確率分布のサイズまできちんと指定しなくていいんですね。

NUTS 1000x1 で 977秒→638秒
Slice 10000x3 で 60秒→47秒

この手の道具の使い方を勉強するときには、ある程度答えを知っているものを選んで遊んでみるのが普通だと思いますが、

・スピード
・アルゴリズムのわかりやすさ
・失敗のわかりやすさ

などから、最初にスライスサンプラーを試してみるのはありだと思いました。

他のMCMCパッケージがどれだけのサンプラーに対応しているかは知りませんが。

Mamba.jl (Juliaで使えるMCMCパッケージ、StanのJulia版みたいなもの)のExamplesを見るとSliceが結構使われている。