![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
Rコードの最適化（初心者用メモ）
以下は r-help の投稿記事です。参考になると思うので、記録しておきます。
新人 R プログラマーへ:
この長いノートは R におけるコードの最適化の例を与えます。初心者向けですから、そうでない人は無視して下さい。私の希望は、R におけるコードの改良、ループ処理、ベクトル化といったいくつかの基本的な注意点を明らかにすることを手助けできたらいいなということです。しかしながら、私は、このノートを拡張し、修正し、そして改良するような他の人からのコメントを歓迎します。私は依然として R プログラミングの初心者です。
参考までに述べれば、私の使用機械／ＯＳは 1.5 GHz Windows 2000 マシンで R 1.9.0 を使っています。
作業は以下の通りです：mydat という非負の値からなるデータセットと、それらの値よりはるかに大きな定数 K があります。抽出された値の和がちょうど K を越えるのに必要な mydat からのブートストラップ抽出（無作為復元抽出）の数の分布をシミュレーションしたいとします。 つまり、結果の和が K を越えるまで mydat からブートストラップで標本抽出を繰り返します。 これが７回だったとしましょう。すると、７が私の関心のある分布からのサイズ１の無作為標本ということになります。分布を推定するために、この手順を多数回, nsim 回、繰り返します。
先に進む前に、R の初心者にどうしたら、もしくはどうすればうまく、これを実行できるコードを書けるか考えてみることをお薦めします。
最初のアプローチ（素朴／直接的）：入れ子になったループを用い、上に解説した通りにシミュレーションを実行する。外側のループは 1 から nsim を動き、必要だった抽出回数を記録する。 内側のループは繰返し mydat からサイズ 1 の標本を抽出し、総和を取り、それが K を越えるまでチェックし、必要だった回数を外側のループに返す。
これは内側のループに "while" ループを使えば簡単にプログラムでき、私の環境では、適当な nsim (5000-10000) で短い（数秒）時間で実行できます。しかし、もっとうまくやれるはず...。
アプローチ２：すべての抽出を「一度に」やる以外は同じことをする。つまり、（１）の内側のループの各実行は、サイズ 1 の標本を得るために R の "sample()" 関数を呼び出して実行します。 これは各繰返し毎に関数呼び出しをするという大きなオーバーヘッドがあるために、非常に非効率的です。 sample() をたった一度だけ、もしくは最大でもわずかな回数、呼び出し、大きな標本を得て、それに対して添字操作をすれば、これは避けられます。 これは、「大きな標本」とはどのくらいか見積り、添字操作を少々慎重に行なわなければならないので、アプローチ１よりも少しだけ複雑になります。 しかし、（特にポインタ操作を行なえる C プログラマにとって）これは特に困難というわけではなく、無作為抽出は非常に高速なため、安全のため単に必要なサイズより（5 から 10 倍も）大きな標本 WAYYY を生成するのに何の問題もありません。もしくは、そんなに大きくはない標本を生成し、外側のループ毎にもう少し標本を生成する必要があるかどうかチェックします。
このように、内側のループで関数呼び出しでなく添字操作を行なうことにすれば、かなりの改良 （5 から 10 倍）が得られました。分布を生成するのに必要な時間は、一瞬（1 秒以下）でした。 これはいかなる状況下でも、私の必要に対して十分過ぎるものでした。 しかしそれでもなお、自尊心のある R プログラマはループを避けるものとされ、ループの内部にはまだこの醜い（？）ループが残っています。もう少し励んでみることにしましょう。
アプローチ３：したがって、問題は、如何にして内側のループを取り除くかということです。続ける前に、再び R 初心者にこれについて考えてみることを薦めます。
ここでの核心は、内側のループが単なる累積和という非常に簡単な計算をしていることです。 なるほど、回答はすぐそこです。R は既にこれを行なう（背後にある C コードでループを行なう） cumsum() 関数をもっていますから、それをどう使うか首をひねるだけです。
これを行なう本質的なアイデアは、大きな mydat 標本から、総和が K を越えることが保証されているだけ十分大きな、小さな標本を取り出すことです。再び浪費に耽ることができます。例えば 総和が K を越えるには平均で約 10 個分の mydat データが必要とすれば、 30 ないし 40 回そうした抽出を仮に行なえば十分でしょう（そして常に検査を行なうか、もし心配でしょうがなければ try() 関数を使います）。smallsamp を私の大きなブートストラップ標本からの次の 30 ないし 40 個の標本とすれば、以下のコードは内部ループをベクトル化します：
count <- sum(cumsum(smallsamp) < K) + 1
ここで使われているトリックに注意しましょう。「< K」は K 未満のすべての累積和に対し論理値 TRUE のベクトルを返します。これは sum() 関数で加えられる時、内緒で数値 1 とみなされます。少し考えれば、何故最後に 1 を加えているか分かるでしょう。
実際、このベクトル化は、私の適当な大きさの nsim に対し、計算時間を約半分、まばたきする程度、に減らしました。
アプローチ４：最後の段階は、もちろん、毎回分布ベクトルに一標本を加える(再び R の禁じ手です)外側のループを取り除くことです。これを行なう一つの方法（私が思いつく唯一の方法ですが、より賢い R プログラマは他の可能性に挑戦して下さい）は apply() 関数を使うことで、これにより原則として常にループを取り除くことができます。残念ながら、apply() 関数がやっていることは賢い R コードからループを隠すことで、ループ自体を除くことではありませんから、これは見かけだおしです。 ここで理解できるように、そのメリットは、添字の管理を自分自身で明示的に行なわず、R 関数に任せることにより、コードの可読性を高めることです。
apply() を使うための簡単なアイデアは、行数が nsim で、十分な長さの列をもち、したがって各列が総和が K 以上の「十分大きな」mydat からの無作為抽出を含むことができるような行列を用意することです。M 行で間に合うとしましょう。そうするとすべてのシミュレーションは以下のようになります：
drawmat <- matrix( sample(mydat, M*nsim, rep = TRUE), ncol = nsim ) # (注）mydat から M*nsim 回復元抽出し、それを列数 nsim （したがって行数 M になる）の # 行列 drawmat に付値 drawmat <- apply(drawmat, 2, cumsum) < K # （注）apply 関数で行毎（apply 関数の引数 2）に累積和をとり、 # さらに K 未満かどうかで TRUE, FALSE に置き換える。drawmat は論理値行列になる result <- colSums(drawmat) + 1 # （注）colSum 関数で行毎の総和を計算し、1 を加えて、ベクトル result に付値する
最初の行はブートストラップ標本を生成します。二・三番めは前と同じトリックを使い、各列に対し、総和が K を越えるのに必要な行数を計算します。列中の 1 を加算するのに、手軽で（高速な）内部関数 colSums() を使っていることを注意しましょう。このコードはいっそう短く、先の明示的なループ使用戦略を使ったコードよりも、可読性が高いと思います。
しかしながら、驚くには当たりませんが、これはアプローチ３のコードよりも（遅くもない代わりに）早くはありません。apply() 関数に隠されたループが依然としてあるからです。
この長ったらしい例が、 とっつきにくいかのように思える言語への R 初心者の足掛かりになればと思います。他の人がそうしているように、私は皆さんに、この help ML に投稿する前に、R の豊富なドキュメントとヘルプ機能、そして Venables と Ripley の S 言語とデータ解析に関するに関する二つの貴重な本を、参考にされるように心から薦めたいと思います。
Cheers,
Bert Gunter
Non-Clinical Biostatistics
Genentech
"統計家の仕事は、科学における学習過程への触媒となることである。" -- George E.P. Box
訳注：上の三種類のアプローチによるコード例と実行時間（投稿者の述べるところとは異なり、アプローチ4は3より若干時間がかかる。恐らく巨大な配列を操作するオーバヘッド？）。
nsim <- 1000 # シミュレーション回数 M <- 500 # 十分大きな数（のはず） mydat <- runif(10000) # 大きなデータセット（のつもり） result <- numeric(nsim) # シミュレーション結果用のベクトル K <- 100 # しきい値
test1 <- function(){
for (i in 1:nsim) {
Sum <- 0
count <- 0
while (Sum < K) {
Sum <- Sum + sample(mydat, 1, replace=TRUE)
count <- count + 1
}
result[i] <<- count
}
}
test2 <- function(){
for (i in 1:nsim) {
x <- sample(mydat, M, replace=TRUE)
# 小心なユーザ向きのチェック(M 個のサンプリングで足りない時警告)
if (sum(x) < K) cat("warning! sum(x) < K at ", i,"-th simulation.\n")
result[i] <<- sum((cumsum(x) < K)) + 1
}
}
test3 <- function() {
drawmat <- matrix( sample(mydat, M*nsim, rep = TRUE), ncol = nsim )
drawmat <- ( apply(drawmat, 2, cumsum) < K)
result <<- colSums(drawmat) + 1
}
> system.time(test1())
[1] 6.37 0.00 6.38 0.00 0.00 # 約６秒
> system.time(test2())
[1] 0.20 0.00 0.21 0.00 0.00 # 約0.2秒
> system.time(test3())
[1] 0.26 0.12 0.38 0.00 0.00 # 約0.4秒
＝＝＝＝＝
原文にはコードが示されていないので，こんなもんだろうかという例を
どんなものでも，関数として書いておくと便利かもしれない
approach1 <- function(mydat, nsim, K)
{
result <- numeric(nsim) # 結果を保存する領域を確保
for (i in 1:nsim) { # nsim 回繰り返す
n <- x <- 0
while (x < K) { # 和が K を超えるまで
x <- x+sample(mydat, 1, rep=TRUE) # 和を計算
n <- n+1 # 何個発生させたか
}
result[i] <- n # 結果を保存
}
return(result) # シミュレーション結果を返す
}
# テストデータを作成
mydat <- runif(100000)
# 実行してみる
result <- approach1(mydat, 1000, 100)
hist(result)
ーーーーー
# アプローチ2に対するものがどんなものなのか，イメージしにくかったのだ
approach3 <- function(mydat, nsim, K, M=100)
{
result <- numeric(nsim)
for (i in 1:nsim) {
result[i] <- sum(cumsum(sample(mydat, M, rep=TRUE)) < K)+1
}
return(result)
}
ーーーーー
# ついでに関数に仕立てた例も
approach4 <- function(mydat, nsim, K, M=100)
{
drawmat<-matrix(sample(mydat,M*nsim,rep=TRUE),ncol=nsim)
drawmat <-apply(drawmat,2,cumsum)<K
colSums(drawmat)+1
}
ーーーーー
コメント
コーディング例はだぶっちゃいました。
result を大域変数にして <<- で付値するのは初心者には勧めない方がいいのかもしれません。
関数内で使用する変数の値は，引数で渡すようにするのがいいと思います。
私は，アプローチ1のプログラムを推奨したい。
