2chのRプログラム

2chに，以下のようなプログラムが投稿されていた。ブラウン運動のシミュレーションだという。for文を使っているので時間が掛かるよし。
series_plus <- function(x) {
len <- length(x)
y <- numeric(len)
for (i in 0:(len-1)) {
y <- c(numeric(i),x[1:(len-i)]) + y
}
y
}
中澤さんが次のようなプログラムを書いた。
http://phi.med.gunma-u.ac.jp/index.html　（第940回）
xの長さが大きくなると，メモリーがきつくなるよし。
series_plus <- function(x) { as.vector( x %*% ifelse(lower.tri(diag(length(x))),0,1) )}
どのようなことをやっているプログラムなのか，よーく読んでみると，長さnのベクトルを受け取って同じ長さのベクトルを返す。戻り値の第i要素は入力値の第1～第i要素までの和。
最初の2chのプログラムは，最初は戻り値を0ベクトルにおいて，入力ベクトルを一つずつ右にシフトしてそのたびに戻り値ベクトルと合計している。それは無駄だ。ということで，以下のようなプログラムを書く。
series_plus <- function(x) {
n <- length(x)
y <- numeric(n)
for (i in 1:length(x)) {
y[i] <- sum(x[1:i])
}
y
}
for 文も決して悪くはないし，場合によってはベストの選択にもなる。
しかし，ここではちょっと遊んでみよう。
sequence(1:n) というのは，へんちょこりんな数列を返す。
> sequence(1:10)
[1] 1 1 2 1 2 3 1 2 3 4 1 2 3 4 5 1 2 3 4 5
[21] 6 1 2 3 4 5 6 7 1 2 3 4 5 6 7 8 1 2 3 4
[41] 5 6 7 8 9 1 2 3 4 5 6 7 8 9 10
これを見ると入力値の和を取るための添え字がちゃんと現れることに気づく。でも，どこからどこまでがその添え字かわからないよな。ということで，sequenceのプログラムリストを見てみると，
> sequence
function (nvec)
unlist(lapply(nvec, seq_len))
そうそう，結果を unlist するから区切れ目がなくなっちゃうんだから，unlist しなけりゃいいんだよね。
> lapply(1:5, seq_len)
[[1]]
[1] 1

[[2]]
[1] 1 2

[[3]]
[1] 1 2 3

[[4]]
[1] 1 2 3 4

[[5]]
[1] 1 2 3 4 5
これを入力値の添え字に使って sum 関数を働かせて，ベクトルで返せばよいので，
series_plus <- function(x) {
sapply(lapply(1:length(x), seq_len), function(i) sum(x[i]))
}
まあ，2重のループが隠されているけどいいか。
2chのプログラムでは，
> x <- rnorm(10000)
> system.time(plot(series_plus(x)))
ユーザ システム 経過
2.887 5.984 8.863
三番目の，無駄を省いた for 使用プログラムでは，
> system.time(plot(series_plus(x)))
ユーザ システム 経過
1.685 1.356 3.037
最後のプログラムでは，
> x <- rnorm(10000)
> system.time(plot(series_plus(x)))
ユーザ システム 経過
1.684 1.563 3.245
まあ，そんなものかな？for を使ったのと，ほぼ同じかな。
ただし，このプログラム，メモリの消費は少ないと思いきや，入力ベクトルが長くなると，計算途中でガーベージ・コレクションがはじまるのか，時間ばかりが過ぎてゆく。

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

そもそも「戻り値の第i要素は入力値の第1～第i要素までの和。」と書いたときに気づくべきだった。「これはまさに，cumsum だ！」と。実行速度は当然ながら，そうとうの速さだ。