ダメ出し：マルコフ連鎖モンテカルロ法

2012-01-28  Rでマルコフ連鎖モンテカルロ法を試す の後半「3. MCMC（M-Hアルゴリズム）で推定してみる」で，

 

for (n in 1:10000) を 100 倍の for (n in 1:1000000) にしたら，10 分経っても終わらない。Stop させてみるとやっと予定の半分くらいしかできていなかった。for (n in 1:100000) では 15 秒ほどで終わる。逐次的に  s <- c(s, s0) をやっているのがその原因。

 

以下のプログラムで，burn-in の 500 回を除いて，実質 1000000 個の s を求めると，36 秒で終わる。ちなみに，元のプログラムでは n=1000000 としても，有効な s はその 8 割くらいしか記録されない。

 

system.time({
set.seed(2658817)
s0 <- 0.1
LL0 <- exp(-mlpoi(s0)) # -1 を掛ける必要はない
burn.in <- 500
k <- burn.in + 1000000
s <- numeric(k) # 前もって必要個数を確保
runif1 <- runif(k) # 乱数も前もって必要個数確保
for (n in seq_len(k)) {
    while ((s1 <- s0 + rnorm(1)) <= 0) { } # ループ本体は「空」
    LL1 <- exp(-mlpoi(s1))
    if (min(1, LL1/LL0) > runif1[n]) { # 不必要な一時変数は使わない
        s0 <- s1
        LL0 <- LL1
    }
    s[n] <- s0
}
s <- s[-(1:burn.in)]
cat(sprintf("lambda:%.5f variance:%.5f\n", mean(s), var(s))) # %1.5f というのは変な指定
})
length(s)

# lambda:0.85257 variance:0.00850
#    ユーザ   システム       経過 
#     36.474      0.787     36.969
# > length(s)
# [1] 1000000