餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。

M-Hアルゴリズムによるポアソン分布推定コードのチューニング

Rでマルコフ連鎖モンテカルロ法を試す」で書いたMCMCのコードが、NR法による最尤推定とのつながりで混乱していたせいか、それをチューニングしようと言うエントリーが出てしまったので、ちょっと補足をします。
以前のエントリーでは対数尤度関数を使っているので引数が負になるとエラーになる事から、無駄な分岐が増えている上に、サンプル数を切り捨て、配列サイズが動的に決まる事になっています。問題点をフォローしてもらっているわけですが、そもそもMCMC対数尤度関数を使う必要がありません。尤度比でマルコフ過程を導出するので対数化したものを指数化しているとか奇妙です。
必要なのは事後確率1/事後確率0=尤度1/尤度0=目的関数の値1/目的関数の値0を満たす目的関数なので、チューニング・ポイントとしては目的関数を簡素化してしまう方が効果的です。

# セットアップ
set.seed(1631697)
lambda <- 1
x <- rpois(100, lambda)

# 目的関数を定義
sum_x <- sum(x)
lpoi <- function(p){
# (p>0)*exp(-length(x)*p)*p^(sum(x))/prod(factorial(x))の分母と固定係数を簡素化
(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))
  for(n in 2:length(s)){
    s1 <- s0 + r[n]
    LL1 <- lpoi(s1)
    if(min(1, LL1/LL0) > runif(1)){
      s0 <- s1
      LL0 <- LL1
    }
    s[n] = s0
  }
  s <- tail(s, length(s) - 500)
  print(sprintf("lambda:%1.5f variance:%1.5f", mean(s), var(s)))
})

上記のように簡素化すると手元の環境では1.67秒になり、チューニングしてくれたエントリーの5.35秒よりも高速になります。
グダグダなコードを改良してもらって恐縮なのですが、そもそもの目的関数がちょっと問題が多かったと言う事のようです。