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:一緒、とは言ってない