ハイブリッドモンテカルロ試してみた。

PRML読書会でなんかうまくいかない的な話になったので、ちょっと書いてみた。
p(z) = N(0,1) としている。

# Hybrid Monte Carlo sampling

N <- 1000; # number of sampling
leapfrog_count <- 100;
leapfrog_epsilon <- 0.01;

# p(z) = N(0,1) = exp(-z^2/2)/sqrt(2pi), E(z) = z^2/2, Zp = sqrt(2pi)
# dE/dz = z

# Hamiltonian: H(z,r) = E(z) + K(r) = z^2/2 + r^2/2

r <- rnorm(1,0,1);
z <- 0;
zlist <- numeric(N);

for(i in 1:N) {
	# leapfrog
	e <- sample(c(-1,1), 1) * leapfrog_epsilon;
	r <- r - e * z / 2; # r(t+e/2)
	for(j in 1:leapfrog_count) {
		z <- z + e * r;
		r <- r - e * z;
	}
	r <- r - e * z / 2; # r(t+e/2)

	zlist[i] <- z;

	# resampling of r
	# p(r|z) = p(r) = N(0,1);
	r <- rnorm(1,0,1);
}
hist(zlist);

リープフロッグと r のリサンプリングを交互に行えば、大丈夫そう。
左右対称でないのは、棄却してないせいかな。サンプル数を増やしてもその傾向があるし。あとで棄却やってみる。

イメージとしては、リープフロッグは p(z|r) からのサンプリングに相当、r も一緒に変わっているけど、H 一定という条件下なので、実質変わってないのと同じ(z と相関して決まっている)。
続けて p(r|z)=p(r) から r をサンプリングすれば H を変えることが出来る。
これはギブス的に見なせばごく自然。というか、z をサンプリングするときは H を止めていて(つまり等高線の中で動かしていて)、r のリサンプリングはある意味 H のリサンプリングだと思えば、スライスサンプリングにそっくりやん！*1　と考えてスッキリしたのだが、勘違いだったら誰か教えて。

r のリサンプリングを行わないと、wk さんが書いたのと同じようなヒストグラムに。

「H 一定だとエルゴードでない」、というのはまさにこの状態。

*1:一緒、とは言ってない