Rコードの最適化（初心者用メモ）の編集

//mase 2004.5.20

SIZE(18){COLOR(blue){Rコードの最適化（初心者用メモ）}}

以下は r-help の投稿記事です。参考になると思うので、記録しておきます。

新人 R プログラマーへ:

この長いノートは R におけるコードの最適化の例を与えます。初心者向けですから、そうでない人は無視して下さい。私の希望は、R におけるコードの改良、ループ処理、ベクトル化といったいくつかの基本的な注意点を明らかにすることを手助けできたらいいなということです。しかしながら、私は、このノートを拡張し、修正し、そして改良するような他の人からのコメントを歓迎します。私は依然として R プログラミングの初心者です。

参考までに述べれば、私の使用機械／ＯＳは 1.5 GHz Windows 2000 マシンで R 1.9.0 を使っています。

作業は以下の通りです：mydat という非負の値からなるデータセットと、それらの値よりはるかに大きな定数 K があります。抽出された値の和がちょうど K を越えるのに必要な mydat からのブートストラップ抽出（無作為復元抽出）の数の分布をシミュレーションしたいとします。 
//The task is this: I have a data set, mydat, consisting of values >=0, and a //fixed constant, K, that is "much" larger than any of the values. I want to //simulate the distribution of the number of bootstrap draws
//(random sampling with replacement) from mydat that is needed  to make
//the sum of the values drawn just exceed K. 
つまり、結果の和が K を越えるまで mydat からブートストラップで標本抽出を繰り返します。
これが７回だったとしましょう。すると、７が私の関心のある分布からのサイズ１の無作為標本ということになります。分布を推定するために、この手順を多数回, nsim 回、繰り返します。  
//That is, repeatedly bootstrap
//sample from mydat until the sum of the results exceeds K. Suppose this
//requires 7 draws. Then 7 is a random sample of size 1 from the
//distribution I'm interested in. I want to repeat this procedure a
//"large" number, nsim, of times to estimate the distribution.

先に進む前に、R の初心者にどうしたら、もしくはどうすればうまく、これを実行できるコードを書けるか考えてみることをお薦めします。
//Before reading on, I invite R novices to consider how -- or better yet
//write code -- to do this.

最初のアプローチ（素朴／直接的）：入れ子になったループを用い、上に解説した通りにシミュレーションを実行する。外側のループは 1 から nsim を動き、必要だった抽出回数を記録する。
内側のループは繰返し mydat からサイズ 1 の標本を抽出し、総和を取り、それが K を越えるまでチェックし、必要だった回数を外側のループに返す。
//Approach 1 (Naive/ straightforward): Perform the simulation exactly as
//described by using a nested loop. The outer loop goes from 1 to nsim to
//record the number of draws needed. The inner loop repeatedly draws a
//sample of size 1 from mydat and accumulates and checks the sum until it
//exceeds K, returning the number of times it took to do this to the outer
//loop.

これは内側のループに "while" ループを使えば簡単にプログラムでき、私の環境では、適当な nsim (5000-10000) で短い（数秒）時間で実行できます。しかし、もっとうまくやれるはず...。
//This is easy to program up using a "while" loop for the inner loop and
//takes about a short look's time out my window (a few seconds) to run for
//modest nsim (5-10K). But surely one can do better...

アプローチ２：すべての抽出を「一度に」やる以外は同じことをする。つまり、（１）の内側のループの各実行は、サイズ 1 の標本を得るために R の "sample()" 関数を呼び出して実行します。
これは各繰返し毎に関数呼び出しをするという大きなオーバーヘッドがあるために、非常に非効率的です。
//Approach 2: Was essentially the same, except all the sampling is done
//"at once." That is, each execution of the inner loop in (1) required R
//to call the "sample()" function to get a bootstrap sample of size 1.
//This is very inefficient, because there is a high overhead of making
//that function call at each iteration. 
sample() をたった一度だけ、もしくは最大でもわずかな回数、呼び出し、大きな標本を得て、それに対して添字操作をすれば、これは避けられます。
//It can be avoided by calling
//sample() only once -- or at most a small number of times -- to get a
//large sample and then using indexing to work one's way through the large
//sample. 
これは、「大きな標本」とはどのくらいか見積り、添字操作を少々慎重に行なわなければならないので、アプローチ１よりも少しだけ複雑になります。
//This turns out to be a little more complicated than approach 1
//because you first have to figure how large a "large" sample should be
//and then have to be a little careful about setting up your indexing. 
しかし、（特にポインタ操作を行なえる C プログラマにとって）これは特に困難というわけではなく、無作為抽出は非常に高速なため、安全のため単に必要なサイズより（5 から 10 倍も）大きな標本 WAYYY を生成するのに何の問題もありません。もしくは、そんなに大きくはない標本を生成し、外側のループ毎にもう少し標本を生成する必要があるかどうかチェックします。
//But it's not very difficult (especially for any C programmer who can
//manipulate pointers), and random sampling is so fast that it's no
//problem to generate a sample WAYYY bigger than you need (5 or 10 times
//as large, even) just to be safe. Or generate a not so big version and
//just check at each outer loop iteration whether you need to generate
//some more.

このように、内側のループで関数呼び出しでなく添字操作を行なうことにすれば、かなりの改良
（5 から 10 倍）が得られました。分布を生成するのに必要な時間は、一瞬（1 秒以下）でした。
これはいかなる状況下でも、私の必要に対して十分過ぎるものでした。
//Doing it this way -- now using indexing, not function calls in the inner
//loop -- made a considerable improvement (5 - 10 fold). It now took about
//one quick glance out the window time to generate the distribution (less
//than a second). This was more than adequate for my needs under any
//conceivable situation. 
しかしそれでもなお、自尊心のある R プログラマはループを避けるものとされ、ループの内部にはまだこの醜い（？）ループが残っています。もう少し励んでみることにしましょう。
//But still, self-respecting R programmers are
//supposed to eschew looping, and here was this ugly(?) loop within a
//loop!  So let us persevere.

アプローチ３：したがって、問題は、如何にして内側のループを取り除くかということです。続ける前に、再び R 初心者にこれについて考えてみることを薦めます。
//Approach 3: So the question is -- how to remove that inner loop? Again,
//I invite novice R programmers to think about that before continuing
//on...

ここでの核心は、内側のループが単なる累積和という非常に簡単な計算をしていることです。
なるほど、回答はすぐそこです。R は既にこれを行なう（背後にある C コードでループを行なう） cumsum() 関数をもっていますから、それをどう使うか首をひねるだけです。
//The key here is that the inner loop is doing a very simple calculation:
//just accumulating a sum.  Forsooth! -- the answer fairly leaps out: R
//already has a cumsum() function that does this (looping in the
//underlying C code), so one only has to figure out how to make use of it.

これを行なう本質的なアイデアは、大きな mydat 標本から、総和が K を越えることが保証されているだけ十分大きな、小さな標本を取り出すことです。再び浪費に耽ることができます。例えば
総和が K を越えるには平均で約 10 個分の mydat データが必要とすれば、 30 ないし 40 回そうした抽出を仮に行なえば十分でしょう（そして常に検査を行なうか、もし心配でしょうがなければ  try() 関数を使います）。smallsamp を私の大きなブートストラップ標本からの次の 30 ないし 40 個の標本とすれば、以下のコードは内部ループをベクトル化します：
//The essential idea to do this is to get a small sample from the large
//mydat sample that is guaranteed to be bigger than one needs to total up
//to more than K. Again, one can be profligate: if I need about 10 mydat
//values on average to sum up to K, if I sample 30 or 40 or some such
//thing, I'm sure to have enough (and can always add a check or use
//"try()" to make sure if I am faint of heart). If smallsamp is this
//sample of the next 30 or 40 values from my large bootstrap sample, then
//the following code vectorizes the inner loops:

count <- sum(cumsum(smallsamp) < K) + 1

ここで使われているトリックに注意しましょう。「< K」は K 未満のすべての累積和に対し論理値 TRUE のベクトルを返します。これは sum() 関数で加えられる時、内緒で数値 1 とみなされます。少し考えれば、何故最後に 1 を加えているか分かるでしょう。

//Note a trick here: The <K produces a vector of logical TRUE's for all
//cumsum values <K. This is silently treated as numeric 1's when added by
//sum(). A moment's thought will make clear why 1 must be added at the
//end..

実際、このベクトル化は、私の適当な大きさの nsim に対し、計算時間を約半分、まばたきする程度、に減らしました。
//Sure enough, this vectorization  reduced the time about another twofold
//-- to an eyeblink --  for my modest nsim sizes.

アプローチ４：最後の段階は、もちろん、毎回分布ベクトルに一標本を加える(再び R の禁じ手です)外側のループを取り除くことです。これを行なう一つの方法（私が思いつく唯一の方法ですが、より賢い R プログラマは他の可能性に挑戦して下さい）は apply() 関数を使うことで、これにより原則として常にループを取り除くことができます。残念ながら、apply() 関数がやっていることは賢い R コードからループを隠すことで、ループ自体を除くことではありませんから、これは見かけだおしです。
//Approach 4: The last stage, of course, is to get rid of the outer loop,
//which fills the vector for the distribution one element at a time, again
//another R no-no. One way to do this -- the ONLY way I could think of (so
//a challenge to smarter R programmers) -- is to make use of apply(),
//which basically always can remove loops. Unfortunately, this is an
//illusion, because what the apply() functions actually do is hide the
//loops in clever R code, not remove them. 
ここで理解できるように、そのメリットは、添字の管理を自分自身で明示的に行なわず、R 関数に任せることにより、コードの可読性を高めることです。 
//Their virtue, as will be seen
//here, is making the code more transparent by transferring the
//bookkeeping of indexing to the R functions rather than having to
//explicitly handle it yourself.

apply() を使うための簡単なアイデアは、行数が nsim で、十分な長さの列をもち、したがって各列が総和が K 以上の「十分大きな」mydat からの無作為抽出を含むことができるような行列を用意することです。M 行で間に合うとしましょう。そうするとすべてのシミュレーションは以下のようになります：
//To use apply, the simple idea is just to make up a matrix with nsim
//columns with enough rows so that each column contains "enough" random
//draws from mydat to sum to more than K. Let's say that we decide M rows
//will do it. Then the code for the whole simulation is:

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() を使っていることを注意しましょう。このコードはいっそう短く、先の明示的なループ使用戦略を使ったコードよりも、可読性が高いと思います。
//The first statement generates the matrix of bootstrap samples, while the
//second and third use the same trick as previously to get the number of
//rows needed for each column to total more than K. Note that I used the
//handy (and fast) internal colSums function to add up the 1's in the
//columns. This code is much shorter, and I would say more transparent,
//then the code produced by the previous explicit looping strategies.

しかしながら、驚くには当たりませんが、これはアプローチ３のコードよりも（遅くもない代わりに）早くはありません。apply() 関数に隠されたループが依然としてあるからです。
//However, not surprisingly, it does not run any faster (nor slower) than
//the Approach 3 code.The looping is still there hidden in the apply().

この長ったらしい例が、 とっつきにくいかのように思える言語への R 初心者の足掛かりになればと思います。他の人がそうしているように、私は皆さんに、この help ML に投稿する前に、R の豊富なドキュメントとヘルプ機能、そして Venables と Ripley の S 言語とデータ解析に関するに関する二つの貴重な本を、参考にされるように心から薦めたいと思います。
//I hope that this long-winded example will help R newcomers in their
//entry into what may seem like a daunting programming language. As others
//have, I would urge all to consult R's quite remarkable documentation and
//Help facility (BEFORE they post to this list), as well as Venables's and
//Ripley's two indispensable books on S language programming and data
//analysis.

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のプログラムを推奨したい。
-この記事の趣旨は「最適化」であることをお忘れなく。最適化しないでも待てる方はどうなりと。しかし、初心者の内から R プログラミングのコツに留意する癖をつけておかないと、結局 R は遅過ぎて使いものにならないと文句を言うことになりそうです。 --  &new{2004-05-23 (日) 10:56:07};
-[[Rコード最適化のコツと実例集]]を参照。 --  &new{2004-06-06 (日) 06:17:45};

#comment

タイムスタンプを変更しない