メトロポリス法による正規分布からのサンプリング

TokyoNLP #9 で tkng さんが MCMC で正規分布からサンプリングしてみたけど、1000件ぐらいだとなかなかきれいな釣鐘型にならない、という話をされていた。
途中の系列を捨ててないということだったので、そのせいでは？　と思ってちょっと試してみたが、結論としては、そもそも 1000件くらいではきれいなヒストグラムにならないということがわかった。

それは R で「ちゃんとした正規分布からのサンプリング」を使ってヒストグラムを書けばすぐわかる。あえてグラフは載せないので、自分で実行してみてね。

hist(rnorm(1000), breaks=30)

これで終わったら芸がないので、「相関を抑えるにはどれくらい捨てればいいか」を確認してみた。
まずはメトロポリス法による正規分布サンプリングのソース。
スライドで一瞬見ただけなので多分いろいろ違っているだろうがまあ大丈夫？

#!/usr/bin/env ruby

def gauss(x)
    Math.exp( - x * x / 2 )
end

x = 0
samples = []
tries = 0
N = 10000
while samples.length < N
    x1 = x + (rand * 2 - 1)
    if rand < gauss(x1) / gauss(x)
        samples << x1
        x = x1
    end
    tries += 1
end
puts samples
$stderr.puts "rejects / tries = #{tries - N} / #{tries}"

出力を R に取り込んで、acf() で自己相関を出してみる。

x <- read.table('リダイレクトを出力したファイルのパス')
acf(x$V1[seq(1,10000,by=1)])
acf(x$V1[seq(1,10000,by=5)])
acf(x$V1[seq(1,10000,by=10)])
acf(x$V1[seq(1,10000,by=12)])