for ループが遅いなんて，誰が言った？ その3

わたしは，R には何の義理もないので，好き勝手言います...

> ベンチマークで利用しているデータは10,000行1,000列のデータフレームですが、for 文が有利なように恣意的に設定されている気がします。
> colMeans のボトルネックは as.matrix であり、for 文のボトルネックは R 側で行われる繰り返し処理そのものです。
> よって、同じデータサイズであっても1,000行10,000列のデータフレームであれば for 文が極端に遅くなります。

恣意的に選んだつもりはないし，普通，データ行列というのは，行数は列数より多いもの（統計データはそんな風です）。

最初に，システムの記述。
=====================================================================
Mac OS X
バージョン 10.7.3
プロセッサ 2.4 GHz Intel Core 2 Duo
メモリ 4 GB　667 MHz DDR2 SDRAM  要するに，非力なノートパソコンです（泣）
=====================================================================
> sessionInfo()
R version 2.14.1 RC (2011-12-21 r57955)
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

locale:
[1] ja_JP.UTF-8/ja_JP.UTF-8/ja_JP.UTF-8/C/ja_JP.UTF-8/ja_JP.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods 
[7] base    
=====================================================================

> よって、同じデータサイズであっても1,000行10,000列のデータフレームであれば for 文が極端に遅くなります。

確かに，1,000 行 10,000 列のデータフレームだと，次のようになりますね。

> gc(); gc()
         used (Mb) gc trigger (Mb) max used (Mb)
Ncells 340340 18.2     407500 21.8   350000 18.7
Vcells 455195  3.5     905753  7.0   879046  6.8
         used (Mb) gc trigger (Mb) max used (Mb)
Ncells 340297 18.2     467875   25   350000 18.7
Vcells 455104  3.5     905753    7   879046  6.8
> set.seed(123454321)
> nr <- 1000
> nc <- 10000
> x <- matrix(rnorm(nr*nc), nr, nc)
> d <- as.data.frame(x)
> gc(); gc(); system.time({
           used  (Mb) gc trigger  (Mb) max used  (Mb)
Ncells   360991  19.3     667722  35.7   420021  22.5
Vcells 20486276 156.3   30561555 233.2 30560771 233.2
           used  (Mb) gc trigger  (Mb) max used  (Mb)
Ncells   360998  19.3     667722  35.7   420021  22.5
Vcells 20486303 156.3   30561555 233.2 30560771 233.2
+     m1 <- numeric(nc)
+     for (i in 1:nc) {
+         m1[i] <- mean(d[, i])
+     }
+ })
   ユーザ   システム       経過 
     1.120      0.123      1.253
> gc(); gc(); system.time({
           used  (Mb) gc trigger  (Mb) max used  (Mb)
Ncells   362956  19.4     667722  35.7   667722  35.7
Vcells 20502981 156.5   30561555 233.2 30560771 233.2
           used  (Mb) gc trigger  (Mb) max used  (Mb)
Ncells   362954  19.4     667722  35.7   667722  35.7
Vcells 20502987 156.5   30561555 233.2 30560771 233.2
+     m3 <- colMeans(d)
+ })
   ユーザ   システム       経過 
     0.612      0.095      0.700             # ★★ 仰るとおり，colMeans が倍くらい速い（パチパチ）
                                                       # ★★ しかし「極端に遅くなる」わけではない

=====================================================================

でも，これは普通の統計データの状況とはちょっと違うので，せめて，行数は列数と同じにしましょうよ。
ということで，10,000 行 10,000 列のデータフレームで，やってみましょう。

> gc(); gc()
           used  (Mb) gc trigger  (Mb) max used  (Mb)
Ncells   363975  19.5     667722  35.7   667722  35.7
Vcells 20524784 156.6   45054625 343.8 40571168 309.6
           used  (Mb) gc trigger  (Mb) max used  (Mb)
Ncells   363973  19.5     667722  35.7   667722  35.7
Vcells 20524790 156.6   45054625 343.8 40571168 309.6
> set.seed(123454321)
> nr <- 10000 ★★ 行数を 10 倍にしました(^_^;)
> nc <- 10000
> x <- matrix(rnorm(nr*nc), nr, nc)
> d <- as.data.frame(x)
> gc(); gc(); system.time({
            used   (Mb) gc trigger   (Mb)  max used   (Mb)
Ncells    363966   19.5     667722   35.7    667722   35.7
Vcells 200524763 1529.9  312110644 2381.3 312104760 2381.2
            used   (Mb) gc trigger   (Mb)  max used   (Mb)
Ncells    363973   19.5     667722   35.7    667722   35.7
Vcells 200524790 1529.9  312110644 2381.3 312104760 2381.2
+     m1 <- numeric(nc)
+     for (i in 1:nc) {
+         m1[i] <- mean(d[, i])
+     }
+ })
   ユーザ   システム       経過 
     1.402      0.040      1.451 # ★★ for は，原理的に言ってそう増えるわけではない。
> gc(); gc(); system.time({
            used   (Mb) gc trigger   (Mb)  max used   (Mb)
Ncells    363975   19.5     667722   35.7    667722   35.7
Vcells 200527435 1530.0  312110644 2381.3 312104760 2381.2
            used   (Mb) gc trigger   (Mb)  max used   (Mb)
Ncells    363973   19.5     667722   35.7    667722   35.7
Vcells 200527441 1530.0  312110644 2381.3 312104760 2381.2
+     m3 <- colMeans(d)
+ })
   ユーザ   システム       経過 
     3.975     15.272    191.106 # ★★ オオッ。これはなんだ。時間は倍以上掛かるし，システムは 15 秒だって？
                                 # これ以上行数を増やして実験するなんて，もうイヤだ(^_^)

＞  データフレームを行列にしてから計算するのだから遅くてもしょうがないじゃあないか（怒）

それはそうだけど，前にも言ったけど，データフレームについて計算しようというのなら行列に変換するのに必要な時間も勘定されてもしょうがないでしょう。それがいやなら，前にも言ったけど，前もって行列にしておくこと。行列のときに colMeans は速いっていうのは否定していないのだから。怒りの矛先を私に向けられても，戸惑うばかり...

for ループが遅いなんて，誰が言った？ その2

> ベンチマークで利用しているデータは10,000行1,000列のデータフレームですが、for 文が有利なように恣意的に設定されている気がします。
> colMeans のボトルネックは as.matrix であり、for 文のボトルネックは R 側で行われる繰り返し処理そのものです。

恣意的に設定したわけではないけど。状況によってはこんなこともあるということ。シミュレーションなんて，限られた条件の下で実験するわけで，実験条件が不適切だといわれても，ネエ。

for ループを使うやり方は列数に比例して計算時間がかかる。行数が 8000 程度以下なら colMeans の方が速い。10000 程度以上になると逆転する。

as.matrix がボトルネックだということなんだから，そもそも matrix で用が足りるときには，read.table などでデータフレームに読み込まない。あるいは読んでも，すぐさま行列にして，必要なら save でもしておいて，いつも行列を使う（二度と read.table しない）と良いと言うこと。list (data frame) は便利だけど高くつく。
しかし，ああだこうだといっても，あれだけの計算をして計算時間は １ 秒未満だということ。めくじらたてるこたあない。

for ループが遅いなんて，誰が言った？

色々な条件により細かい部分は異なるが，以下のようなプログラムによって分かったことは衝撃的？

nr <- 10000
nc <- 1000
x <- matrix(rnorm(nr*nc), nr, nc)
d <- as.data.frame(x)

# データフレームに対して
## for ループ
gc();gc();system.time({
    m1 <- numeric(nc)
    for (i in 1:nc) {
        m1[i] <- mean(d[,i])
    }
})

## apply
system.time({
    m2 <- apply(d, 2, mean)
})

## colMeans
system.time({
    m3 <- colMeans(d)
})

## sapply
system.time({
    m4 <- sapply(d, mean)
})

# 行列に対して
## for ループ
system.time({
    M1 <- numeric(nc)
    for (i in 1:nc) {
        M1[i] <- mean(x[,i])
    }
})

## apply
system.time({
    M2 <- apply(x, 2, mean)
})

## colMeans
system.time({
    M3 <- colMeans(x)
})

1. 行列に対しての colMeans が一番速い

2. 行列では apply より for の方が速い

3. データフレームでは for が一番速い

4. データフレームでは colMeans はたいして速くない

5. apply はいずれに対しても一番遅い

その他，結果の通り！！

全角数字列を半角数字列に変換する

既存の関数を使おう！！

 

chartr という関数が便利。

> chartr("０１２３４５６７８９", "0123456789", "３４５０２４６８１３５７９")
[1] "3450246813579"

============= 昔の記事 ==============

どこかにあったプログラムは Windows 用で，encoding="cp932" を仮定したものだったので，汎用関数を書いてみた。 遅いかも知れないけど，こんな関数，多少時間がかかろうがなにだろうが，問題ない。 動けばよいのだから，凝る必要もない。アルファベットや記号も変換できるようにしてもよいなあ。

Zen2Han <- function(str) {
    Zen <- c("０", "１", "２", "３", "４", "５", "６", "７", "８", "９")
    Han <- c("0", "1", "2", "3", "4", "5", "6", "7", "8", "9")
    str <- unlist(strsplit(str, ""))
    result <- character(length(str))
    for (i in seq_along(str)) {
        for (j in 1:10) {
            if (str[i] == Zen[j]) {
                result[i] <- Han[j]
                break
            }
        }
    }
    return(paste(result, collapse = ""))
}
Zen2Han("１２３４５６７８９０")
Zen2Han("２１８４")

実行例

> Zen2Han("１２３４５６７８９０")
[1] "1234567890"
> Zen2Han("２１８４")
[1] "2184"

ダメ出し：ベクトルと行列の演算は，場合によって注意が必要

2011-05-21　■Ｒでクロス集計表の残差分析 にて

# クロス集計表の入力
X <- matrix(c(58, 11, 10, 35, 25, 23), nrow=3,
            dimnames=list(c("賛成", "中立", "反対"), c("男性", "女性")))
X
#      男性 女性
# 賛成   58   35
# 中立   11   25
# 反対   10   23

# 比率（縦％）を確認
round(X / apply(X,2,sum) *100, 1)
#      男性 女性
# 賛成 73.4 42.2
# 中立 13.3 31.6
# 反対 12.7 27.7

これに対して，「縦方向%が間違っているんじゃないか」というコメントがある

これは，コメントした人が正解。

round(X / apply(X,2,sum) *100, 1) は，まちがい

X / apply(X,2,sum) は，正しいベクトル演算をやるためには t(t(X) / apply(X,2,sum)) でなければならない。
t(t(foo) ...) は，冗長に見えるかも知れないが，必須。

また，apply(X,2,sum) は colSums(X) とするのが，吉。

で，正解は以下の通り。

> round(t(t(X) / colSums(X)) * 100, 1)
     男性 女性
賛成 73.4 42.2
中立 13.9 30.1
反対 12.7 27.7

蛇足ながら，あたりまえだけど，横%は t(t(foo) ...) は不要。

> round(X / rowSums(X) * 100, 1)
     男性 女性
賛成 62.4 37.6
中立 30.6 69.4
反対 30.3 69.7

ダメ出し対象募集中

なんかちがうな～

もっとすっきりかけないの？

どうしたらいいの？

 

そんな，プログラムが載っているWebページのアドレスを曝してください。お答えします。(^_^;)

ダメ出し：規則性を吸収した汎用性のあるプログラムを書こう

2012年2月12日 (日)　Multiple comparison among 4 groups using Fisher exact test
2012年2月10日 (金)　Multiple comparison among 3 groups using Fisher exact test

5グループ，6グループ，... になったらまたプログラムを書くのだろうか?

同じような内容を繰り返す部分は，combn 関数を使うと簡単に書くことができる。
このやりかただと，p.x に名前を付けるためにもう一度 combn を使わなくてはならないが，それはデメリットには勘定されない。
結果を関数内で書くのはあまりお勧めではない。
必要ならば，print メソッドを書いてやるとよい。

p.fisher <- function(tab) {
    nc <- ncol(tab)
    p.x <- combn(nc, 2, function(x) fisher.test(tab[, x])$p.value)
    names(p.x) <- combn(nc, 2, function(x) sprintf("p.fisher.%d%d", x[1], x[2]))
    n.pairs <- length(p.x)
    hochberg <- cbind(Crude = rev(sort(p.x)), "Hoch-Adj." = rev(sort(p.x)) * seq_along(p.x))
    holm <- cbind(Crude = sort(p.x), "Holm-Adj." = sort(p.x) * rev(seq_along(p.x)))
    ans <- list(hochberg = hochberg, holm = holm)
    class(ans) <- "p.fisher"
    return(ans)
}
print.p.fisher <- function(obj, digits = 4) {
    print(round(obj$hochberg, 4))
    print(round(obj$holm, 4))
}

もっと変態チックな FizzBuzz

速いけど，この先にあるのは，答えの文字列を書き出すというプログラムになりそう(^_^;)

> system.time({
+ limit <- 1000000
+ ans3 <- rep(c("d", "d", "Fizz", "d", "Buzz", "Fizz", "d", "d", "Fizz", "Buzz", "d", "Fizz", "d", "d", "FizzBuzz"), ceiling(limit/15))[1:limit]
+ temp <- ans3=="d"
+ ans3[temp]<- which(temp)
+ })
   ユーザ   システム       経過  
     0.324      0.003      0.347
> all(ans1==ans3)
[1] TRUE

R らしいかもしれないがエレガントではない FizzBuzz

速度を求めるなら以下のように

> i <- 1:1000000
> system.time({
+ ans1 <- ifelse(i %% 15 == 0, "FizzBuzz", ifelse(i %% 3 == 0, "Fizz", ifelse(i %% 5 == 0, "Buzz", i)))
+ })
   ユーザ   システム       経過 
     3.214      0.141      3.334
> system.time({
+ ans2 <- i
+ l <- length(ans2)
+ ans2[1:(l%/%3)*3] <- "Fizz"
+ ans2[1:(l%/%5)*5] <- "Buzz"
+ ans2[1:(l%/%15)*15] <- "FizzBuzz"
+ })
   ユーザ   システム       経過 
     0.473      0.008      0.479
> all(ans1 == ans2)
[1] TRUE

R らしく FizzBuzz

「よしいずの雑記帳 FizzBuzz問題の解答例 R編」

R らしく? ベクトル演算で

i <- 1:100
ifelse(i %% 15 == 0, "FizzBuzz", ifelse(i %% 3 == 0, "Fizz", ifelse(i %% 5 == 0, "Buzz", i)))

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

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

ダメ出し：簡明直截なプログラムを書こう

このブログは，プログラムの速度を追求するのではなく，「えれがんと」なプログラミングを目指しているので誤解ないように。「えれがんと」が何を意味するかは範囲は広い。

で，「2012-01-10  Rと手作業で覚える最尤法」の「2. 尤度関数、対数尤度関数、一階条件から最尤法を試みる」について

このプログラムでは，導関数の定義を expression で行っておいてから，使うときに eval で評価している。なるほど，このようにするとある意味「えれがんと」だなあ。

でも，直接指定しても「えれがんと」でなくなるとも見えない。結果としてわずかだが速いし。

なお，連立方程式を解くとき，solve で逆用列を求めて定数項と行列掛け算するのではなく，定数項を solve の引数に渡す方がよいというのは，数値演算の定石。

もとのプログラム（整形済み）

> system.time({for (i in 1:10000) {
+ x <- c(11, 13, 23)
+ n <- length(x)
+ f1 <- expression(-sum(x - mu) / s2)
+ f2 <- expression(-n / (2 * s2) + sum((x - mu) ^ 2) / (2 * (s2 ^ 2)))
+ g11 <- expression(n / s2)
+ g12 <- expression(sum(x - mu) / (s2 ^ 2))
+ g21 <- expression(sum(mu - x) / (s2 ^ 2))
+ g22 <- expression(n / (2 * (s2 ^ 2)) - sum((x - mu) ^ 2) / (s2 ^ 3))
+ mu <- 10
+ s2 <- 10
+ for (i in 1:10) {
+     m <- matrix(c(mu, s2), 2, 1)
+     f <- matrix(c(eval(f1), eval(f2)), 2, 1)
+     j <- matrix(c(eval(g11), eval(g21), eval(g12), eval(g22)), 2, 2)
+     m <- m - solve(j) %*% f
+     mu <- m[1]
+     s2 <- m[2]
+ #    print(sprintf("[%d] (mu,s2)=(%f,%f)", i, mu, s2))
+ }
+ #print(sprintf("平均%2.3f、分散%2.3f（標準偏差%2.3f）の正規分布",
+ #    mu, s2, sqrt(s2)))
+ }})
   ユーザ   システム       経過  
    14.484      0.062     14.270

書き換えたプログラム

> system.time({for (i in 1:10000) {
+ x <- c(11, 13, 23)
+ n <- length(x)
+ mu <- 10
+ s2 <- 10
+ for (i in 1:10) {
+     m <- c(mu, s2)
+     f <- c(-sum(x - mu)/s2,
+            -n/(2 * s2) + sum((x - mu)^2)/(2 * (s2^2)))
+     j <- matrix(c(n/s2,
+                   sum(mu - x)/(s2^2),
+                   sum(x - mu)/(s2^2),
+                   n/(2 * (s2^2)) - sum((x - mu)^2)/(s2^3)), 2, 2)
+     m <- m - solve(j, f)
+     mu <- m[1]
+     s2 <- m[2]
+ #    cat(sprintf("[%d] (mu, s2) = (%f, %f)\n", i, mu, s2))
+ }
+ #cat(sprintf("平均%.3f、分散%.3f（標準偏差%.3f）の正規分布\n",
+ #    mu, s2, sqrt(s2)))
+ }})
   ユーザ   システム       経過  
     8.878      0.022      8.879

ダメ出し：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