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 はおかしい
})
最新の画像[もっと見る]
-
算額(その2059) 44分前
-
算額(その2057) 4時間前
-
算額(その2056) 19時間前
-
算額(その2055) 20時間前
-
算額(その2054) 1日前
-
算額(その2053) 1日前
-
算額(その2052) 1日前
-
算額(その2051) 1日前
-
算額(その2051) 1日前
-
算額(その2050) 2日前
※コメント投稿者のブログIDはブログ作成者のみに通知されます