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
※コメント投稿者のブログIDはブログ作成者のみに通知されます