2012-02-15 M-Hアルゴリズムによるポアソン分布推定コードのチューニング において
「目的関数を簡素化してしまう方が効果的です。」 ごもっともです。
その他， 何回も同じ引数で関数を呼ぶとか，関数呼び出しのオーバーヘッドも馬鹿にならないこともあり，以下のようにすると約 40% のスピードアップ
set.seed(1631697)
lambda <- 1
length_x <- 100
x <- rpois(length_x, lambda)
sum_x <- sum(x)
lpoi <- function(p) {
(p > 0) * exp(-length_x * p) * p^sum_x
}
gc()
gc()
system.time({
set.seed(2658817)
s0 <- 0.1
LL0 <- lpoi(s0) ## ここはこのままにしておこう
s <- numeric(100000 + 500)
s[1] <- s0
r <- rnorm(length(s))
ru <- runif(length(s)) ## 追加
for (n in 2:length(s)) {
s1 <- s0 + r[n]
LL1 <- (s1 > 0) * exp(-length_x * s1) * s1^sum_x ## lpoi(s1)
if (min(1, LL1 / LL0) > ru[n]) {
s0 <- s1
LL0 <- LL1
}
s[n] <- s0
}
s <- s[-(1:500)]
cat(sprintf("lambda:%.5f, variance:%.5f\n", mean(s), var(s))) ## %1.5f はおかしい
})
