ダメ出し：小さなことからコツコツと

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 はおかしい
})