ダメ出し：ifelse はなるべく使わない

2012-02-10　Rでトービット・モデルによる打ち切りデータの推定　の「2. optim()関数で最尤推定」で ifelse を使っているのは，ダメダメ。

掛け算の方がまだ速い。これで，元のプログラムより 25% 速くなった（元のプログラムが 0.104 sec. で，このプログラムが 0.079 sec. まあ，全然どうってことはないのです(^_^)）

> system.time({
+     df <- read.table("kanazawa_accumulation_of_snow.txt",
+         header = TRUE, sep = "\t")
+     frml <- AS ~ TEMP
+     r_ols <- lm(frml, data = df)
+     p <- c(sqrt(deviance(r_ols) / df.residual(r_ols)), coefficients(r_ols))
+     X <- model.matrix(frml, data = df)
+     y <- df$AS
+     print(optim(p, function(p) {
+         sigma <- p[1]
+         xb <- X %*% p[-1]
+         -sum((y <= 0) * log(1 - pnorm(xb / sigma)) +
+              (y > 0)  * log(dnorm((y - xb) / sigma) / sigma))
+     }, method = "BFGS", hessian = TRUE))
+ })
$par
            (Intercept)        TEMP
   24.99719   135.28077   -17.73532

$value
[1] 170.5886

$counts
function gradient
     105       99

$convergence
[1] 0

$message
NULL

$hessian
                        (Intercept)       TEMP
            0.118920177 0.007176688 0.05931874
(Intercept) 0.007176688 0.062764165 0.34602183
TEMP        0.059318737 0.346021832 2.11862555

   ユーザ   システム       経過 
     0.079      0.004      0.094

ダメ出し：「もっといい方法があるはず」と常に考える

配列の拡張とパフォーマンス

差を強調したいがために，基本を忘れている。

> num = 50000
> ### 配列を拡張していく ###
> gc(); gc()
         used (Mb) gc trigger (Mb) max used (Mb)
Ncells 340948 18.3     597831 32.0   597831 32.0
Vcells 602119  4.6    1162592  8.9  1162586  8.9
         used (Mb) gc trigger (Mb) max used (Mb)
Ncells 340955 18.3     597831 32.0   597831 32.0
Vcells 602146  4.6    1162592  8.9  1162586  8.9
> system.time({
+     x <- numeric(0)
+     i <- 1
+     while (i <= num) {
+         x[i] <- i
+         i <- i + 1
+     }
+ })
   ユーザ   システム       経過 
     5.788      0.317      6.078
> ### 配列を一度に初期化 ###
> gc(); gc()
         used (Mb) gc trigger (Mb) max used (Mb)
Ncells 340957 18.3     597831 32.0   597831 32.0
Vcells 602140  4.6    1162592  8.9  1162592  8.9
         used (Mb) gc trigger (Mb) max used (Mb)
Ncells 340955 18.3     597831 32.0   597831 32.0
Vcells 602146  4.6    1162592  8.9  1162592  8.9
> system.time({
+     x <- numeric(num)
+     i <- 1
+     while (i <= num) {
+         x[i] <- i
+         i <- i + 1
+     }
+ })
   ユーザ   システム       経過 
     0.181      0.001      0.193

while は，初期化，判定，加算が必要なので，無駄が多い。
for の替わりに使うべきものではない。

> ### 配列を一度に初期化 ### for があるんだからfor を使う
> gc(); gc()
         used (Mb) gc trigger (Mb) max used (Mb)
Ncells 340957 18.3     597831 32.0   597831 32.0
Vcells 621526  4.8    1162592  8.9  1162592  8.9
         used (Mb) gc trigger (Mb) max used (Mb)
Ncells 340955 18.3     597831 32.0   597831 32.0
Vcells 621532  4.8    1162592  8.9  1162592  8.9
> system.time({
+     x <- numeric(num)
+    for (i in 1:num) {
+         x[i] <- i
+     }
+ })
   ユーザ   システム       経過 
     0.128      0.000      0.165

> ### 配列を一度に初期化 ### for,  seq.int を使う
> gc(); gc()
         used (Mb) gc trigger (Mb) max used (Mb)
Ncells 340957 18.3     597831 32.0   597831 32.0
Vcells 621526  4.8    1162592  8.9  1162592  8.9
         used (Mb) gc trigger (Mb) max used (Mb)
Ncells 340955 18.3     597831 32.0   597831 32.0
Vcells 621532  4.8    1162592  8.9  1162592  8.9
> system.time({
+     x <- numeric(num)
+    for (i in seq.int(num)) {
+         x[i] <- i
+     }
+ })
   ユーザ   システム       経過 
     0.105      0.001      0.105

> # 瞬殺
> gc(); gc()
         used (Mb) gc trigger (Mb) max used (Mb)
Ncells 340957 18.3     597831 32.0   597831 32.0
Vcells 644789  5.0    1162592  8.9  1162592  8.9
         used (Mb) gc trigger (Mb) max used (Mb)
Ncells 340955 18.3     597831 32.0   597831 32.0
Vcells 644795  5.0    1162592  8.9  1162592  8.9
> system.time({x <- seq(num)}) # これを忘れちゃいけないでしょう
   ユーザ   システム       経過 
         0          0          0
> system.time({x <- seq.int(num)})
   ユーザ   システム       経過 
         0          0          0
> system.time({x <- seq_len(num)})
   ユーザ   システム       経過 
     0.001      0.000      0.000

ダメ出し：マルコフ連鎖モンテカルロ法

2012-01-28  Rでマルコフ連鎖モンテカルロ法を試す の後半「3. MCMC（M-Hアルゴリズム）で推定してみる」で，

 

for (n in 1:10000) を 100 倍の for (n in 1:1000000) にしたら，10 分経っても終わらない。Stop させてみるとやっと予定の半分くらいしかできていなかった。for (n in 1:100000) では 15 秒ほどで終わる。逐次的に  s <- c(s, s0) をやっているのがその原因。

 

以下のプログラムで，burn-in の 500 回を除いて，実質 1000000 個の s を求めると，36 秒で終わる。ちなみに，元のプログラムでは n=1000000 としても，有効な s はその 8 割くらいしか記録されない。

 

system.time({
set.seed(2658817)
s0 <- 0.1
LL0 <- exp(-mlpoi(s0)) # -1 を掛ける必要はない
burn.in <- 500
k <- burn.in + 1000000
s <- numeric(k) # 前もって必要個数を確保
runif1 <- runif(k) # 乱数も前もって必要個数確保
for (n in seq_len(k)) {
    while ((s1 <- s0 + rnorm(1)) <= 0) { } # ループ本体は「空」
    LL1 <- exp(-mlpoi(s1))
    if (min(1, LL1/LL0) > runif1[n]) { # 不必要な一時変数は使わない
        s0 <- s1
        LL0 <- LL1
    }
    s[n] <- s0
}
s <- s[-(1:burn.in)]
cat(sprintf("lambda:%.5f variance:%.5f\n", mean(s), var(s))) # %1.5f というのは変な指定
})
length(s)

# lambda:0.85257 variance:0.00850
#    ユーザ   システム       経過 
#     36.474      0.787     36.969
# > length(s)
# [1] 1000000