-1*f()より-f()

Rで「-1*f()が気持ち悪い」と言う指摘を受けて、-f()だとSyntax errorになりそうで逆に気持ち悪いんじゃないかと思っていたら、年寄りの戯言でした。
どういう事かと言うと、この二つがRではS式が異なるわけです。

f <- function(n){ n + 1 }
as.list(quote(-f(2)))
as.list(quote(-1*f(2)))

実行結果は以下のようになります。Rにはオペラント一つの演算子"-"があるのか！*1

	1	2	3
-1*f()	`-`	f(2)
-f()	`*`	-1	f(2)

速度的には大きく影響する部分にならないと思いますが、以下のような単純な100万回ループで2.42秒が2.05秒になったりするので、細かいチューニング・ポイントとして覚えておいても良いかも知れません。

gc(); gc()
n <- 0
system.time({
  for(i in 1:1000000){
    n <- n - f(7)
  }
})
print(n)

gc(); gc()
n <- 0
system.time({
  for(i in 1:1000000){
    n <- n - 1 * f(7)
  }
})
print(n)

*1:Lispでも(- 1)と書けばいいわけだし、驚く事は無いです

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

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

# セットアップ
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秒よりも高速になります。
グダグダなコードを改良してもらって恐縮なのですが、そもそもの目的関数がちょっと問題が多かったと言う事のようです。