ハイブリッドモンテカルロをもっと試してみた

ハイブリッドモンテカルロ試してみた。 - Mi manca qualche giovedi`? の続き。
PRML 読書会 #15の真っ最中に取り急ぎで書き散らかしたコードだったので、Metropolis-Hastings で決まる棄却もやってないし、もうちょっと別の確率分布でも試してみたい。

というわけで改良版コード。

http://github.com/shuyo/iir/blob/master/sampling/hmc.r
- 棄却を行うようにした
- リープフロッグ積分のεと繰り返し回数を毎回変化
- N(0,1) の他にガンマ分布と多峰分布を試す

前回記事でも考察してたように、リープフロッグ積分は p(z|H) からのサンプリングに相当、と考えれば、εと繰り返し回数が固定なのはむしろおかしい。
z の値やイテレーションの回数とかに関係なく固定の範囲で変動させれば、詳細釣り合い条件も壊さないから問題ないだろう。

ハイブリッドモンテカルロの実装上のポイントは、E(z) によらず r が N(0,1) に従う、というところ。
そのため、E(z) と $\frac{\partial E}{\partial z}$ さえあれば、ハイブリッドモンテカルロが適用可能となる(正規化係数すらいらない)。

ガンマ分布は $\rm{Gamma}(a,b) = \frac{1}{Zp}\exp((a-1)\ln z - bz)$ なので

$E(z)=bz-(a-1)\ln z$
$\frac{\partial E}{\partial z}=b-\frac{a-1}{z}$

となる。

正規分布 N(0,1) とガンマ分布 Gamma(3,2) のハイブリッドモンテカルロによるサンプルのヒストグラムが以下の通り。

こちらは正規分布の自己相関。ガンマ分布でも同様のグラフになる。

通常のサンプリング法では、サンプル列から何回かおきに取り出さなければ十分な独立性を確保できないものなのだが、ハイブリッドモンテカルロでは捨てなくても良さそうなくらい最初から自己相関が低い。
ハイブリッドモンテカルロが、ランダムウォーク的振る舞い(移動距離が時間の平方根に比例)という足の短さに対する一つの解決方法である、というのがとても納得できる結果。

もう一つ、やっぱり多峰な分布を試しておきたい。
E(z) を適当にでっちあげて多峰性のある分布を作ってみる(サンプリングするだけなら正規化係数が要らないからラクチン)。1/12 は峰の尖り具合/広がり具合を調整するために付いている。

$E(z) = \frac{1}{12} z(z-1)(z-4)(z-6)$

分布とサンプリング結果のヒストグラム、そして自己相関を見る。

ヒストグラムはちゃんと分布に沿っているが、自己相関がかなり大きい＆なかなか小さくならない。
さては、とステップをプロットすると、やはり峰からなかなか出られない状態に陥っていた。

考察のために p(z,r) を描いてみる。

E <- function(z)z*(z-1)*(z-4)*(z-6)/12;
exp_H <- function(z,r)exp(-E(z)-r^2/2);
z <- seq(-1, 7, length=30);
r <- seq(-3, 3, length=30);
par(mai=c(0,0.1,0,0.1));
persp(z, r, outer(z, r, exp_H), theta=-30, phi=15, zlab="exp(-H(z,r))");

峰のしばりはハイブリッドモンテカルロでもそれなりにキツイ。
が、r が大きくなれば z の値に寄らず等高線が峰の外に出ることもわかる。これがハイブリッドモンテカルロの嬉しいところなのだろう。
分布の形に応じてパラメータを調整する必要も出てくるだろうが、リープフロッグ積分をもっとちょっと長めに回すとか工夫次第でいろいろ改善できそうだ。