Rコードの最適化（初心者用メモ）
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
検索
|
最終更新
|
ヘルプ
]
開始行:
//mase 2004.5.20
SIZE(18){COLOR(blue){Rコードの最適化（初心者用メモ）}}
以下は r-help の投稿記事です。参考になると思うので、記録...
新人 R プログラマーへ:
この長いノートは R におけるコードの最適化の例を与えます。...
参考までに述べれば、私の使用機械／ＯＳは 1.5 GHz Windows ...
作業は以下の通りです：mydat という非負の値からなるデータ...
//The task is this: I have a data set, mydat, consisting ...
//(random sampling with replacement) from mydat that is n...
//the sum of the values drawn just exceed K.
つまり、結果の和が K を越えるまで mydat からブートストラ...
これが７回だったとしましょう。すると、７が私の関心のある...
//That is, repeatedly bootstrap
//sample from mydat until the sum of the results exceeds ...
//requires 7 draws. Then 7 is a random sample of size 1 f...
//distribution I'm interested in. I want to repeat this p...
//"large" number, nsim, of times to estimate the distribu...
先に進む前に、R の初心者にどうしたら、もしくはどうすれば...
//Before reading on, I invite R novices to consider how -...
//write code -- to do this.
最初のアプローチ（素朴／直接的）：入れ子になったループを...
内側のループは繰返し mydat からサイズ 1 の標本を抽出し、...
//Approach 1 (Naive/ straightforward): Perform the simula...
//described by using a nested loop. The outer loop goes f...
//record the number of draws needed. The inner loop repea...
//sample of size 1 from mydat and accumulates and checks ...
//exceeds K, returning the number of times it took to do ...
//loop.
これは内側のループに "while" ループを使えば簡単にプログラ...
//This is easy to program up using a "while" loop for the...
//takes about a short look's time out my window (a few se...
//modest nsim (5-10K). But surely one can do better...
アプローチ２：すべての抽出を「一度に」やる以外は同じこと...
これは各繰返し毎に関数呼び出しをするという大きなオーバー...
//Approach 2: Was essentially the same, except all the sa...
//"at once." That is, each execution of the inner loop in...
//to call the "sample()" function to get a bootstrap samp...
//This is very inefficient, because there is a high overh...
//that function call at each iteration.
sample() をたった一度だけ、もしくは最大でもわずかな回数、...
//It can be avoided by calling
//sample() only once -- or at most a small number of time...
//large sample and then using indexing to work one's way ...
//sample.
これは、「大きな標本」とはどのくらいか見積り、添字操作を...
//This turns out to be a little more complicated than app...
//because you first have to figure how large a "large" sa...
//and then have to be a little careful about setting up y...
しかし、（特にポインタ操作を行なえる C プログラマにとって...
//But it's not very difficult (especially for any C progr...
//manipulate pointers), and random sampling is so fast th...
//problem to generate a sample WAYYY bigger than you need...
//as large, even) just to be safe. Or generate a not so b...
//just check at each outer loop iteration whether you nee...
//some more.
このように、内側のループで関数呼び出しでなく添字操作を行...
（5 から 10 倍）が得られました。分布を生成するのに必要な...
これはいかなる状況下でも、私の必要に対して十分過ぎるもの...
//Doing it this way -- now using indexing, not function c...
//loop -- made a considerable improvement (5 - 10 fold). ...
//one quick glance out the window time to generate the di...
//than a second). This was more than adequate for my need...
//conceivable situation.
しかしそれでもなお、自尊心のある R プログラマはループを避...
//But still, self-respecting R programmers are
//supposed to eschew looping, and here was this ugly(?) l...
//loop! So let us persevere.
アプローチ３：したがって、問題は、如何にして内側のループ...
//Approach 3: So the question is -- how to remove that in...
//I invite novice R programmers to think about that befor...
//on...
ここでの核心は、内側のループが単なる累積和という非常に簡...
なるほど、回答はすぐそこです。R は既にこれを行なう（背後...
//The key here is that the inner loop is doing a very sim...
//just accumulating a sum. Forsooth! -- the answer fairl...
//already has a cumsum() function that does this (looping...
//underlying C code), so one only has to figure out how t...
これを行なう本質的なアイデアは、大きな mydat 標本から、総...
総和が K を越えるには平均で約 10 個分の mydat データが必...
//The essential idea to do this is to get a small sample ...
//mydat sample that is guaranteed to be bigger than one n...
//to more than K. Again, one can be profligate: if I need...
//values on average to sum up to K, if I sample 30 or 40 ...
//thing, I'm sure to have enough (and can always add a ch...
//"try()" to make sure if I am faint of heart). If smalls...
//sample of the next 30 or 40 values from my large bootst...
//the following code vectorizes the inner loops:
count <- sum(cumsum(smallsamp) < K) + 1
ここで使われているトリックに注意しましょう。「< K」は K ...
//Note a trick here: The <K produces a vector of logical ...
//cumsum values <K. This is silently treated as numeric 1...
//sum(). A moment's thought will make clear why 1 must be...
//end..
実際、このベクトル化は、私の適当な大きさの nsim に対し、...
//Sure enough, this vectorization reduced the time about...
//-- to an eyeblink -- for my modest nsim sizes.
アプローチ４：最後の段階は、もちろん、毎回分布ベクトルに...
//Approach 4: The last stage, of course, is to get rid of...
//which fills the vector for the distribution one element...
//another R no-no. One way to do this -- the ONLY way I c...
//a challenge to smarter R programmers) -- is to make use...
//which basically always can remove loops. Unfortunately,...
//illusion, because what the apply() functions actually d...
//loops in clever R code, not remove them.
ここで理解できるように、そのメリットは、添字の管理を自分...
//Their virtue, as will be seen
//here, is making the code more transparent by transferri...
//bookkeeping of indexing to the R functions rather than ...
//explicitly handle it yourself.
apply() を使うための簡単なアイデアは、行数が nsim で、十...
//To use apply, the simple idea is just to make up a matr...
//columns with enough rows so that each column contains "...
//draws from mydat to sum to more than K. Let's say that ...
//will do it. Then the code for the whole simulation is:
drawmat <- matrix( sample(mydat, M*nsim, rep = TRUE), nc...
# (注）mydat から M*nsim 回復元抽出し、それを列数 nsim ...
# 行列 drawmat に付値
drawmat <- apply(drawmat, 2, cumsum) < K
# （注）apply 関数で行毎（apply 関数の引数 2）に累積和を...
# さらに K 未満かどうかで TRUE, FALSE に置き換える...
result <- colSums(drawmat) + 1
# （注）colSum 関数で行毎の総和を計算し、1 を加えて、ベ...
最初の行はブートストラップ標本を生成します。二・三番めは...
//The first statement generates the matrix of bootstrap s...
//second and third use the same trick as previously to ge...
//rows needed for each column to total more than K. Note ...
//handy (and fast) internal colSums function to add up th...
//columns. This code is much shorter, and I would say mor...
//then the code produced by the previous explicit looping...
しかしながら、驚くには当たりませんが、これはアプローチ３...
//However, not surprisingly, it does not run any faster (...
//the Approach 3 code.The looping is still there hidden i...
この長ったらしい例が、 とっつきにくいかのように思える言語...
//I hope that this long-winded example will help R newcom...
//entry into what may seem like a daunting programming la...
//have, I would urge all to consult R's quite remarkable ...
//Help facility (BEFORE they post to this list), as well ...
//Ripley's two indispensable books on S language programm...
//analysis.
Cheers,~
Bert Gunter~
~
Non-Clinical Biostatistics~
Genentech~
"統計家の仕事は、科学における学習過程への触媒となることで...
訳注：上の三種類のアプローチによるコード例と実行時間（投...
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,"-t...
result[i] <<- sum((cumsum(x) < K)) + 1
}
}
test3 <- function() {
drawmat <- matrix( sample(mydat, M*nsim, rep = TRUE), ...
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)) < ...
}
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コード最適化のコツと実例集]]を参照。 -- &new{2004-06...
#comment
終了行:
//mase 2004.5.20
SIZE(18){COLOR(blue){Rコードの最適化（初心者用メモ）}}
以下は r-help の投稿記事です。参考になると思うので、記録...
新人 R プログラマーへ:
この長いノートは R におけるコードの最適化の例を与えます。...
参考までに述べれば、私の使用機械／ＯＳは 1.5 GHz Windows ...
作業は以下の通りです：mydat という非負の値からなるデータ...
//The task is this: I have a data set, mydat, consisting ...
//(random sampling with replacement) from mydat that is n...
//the sum of the values drawn just exceed K.
つまり、結果の和が K を越えるまで mydat からブートストラ...
これが７回だったとしましょう。すると、７が私の関心のある...
//That is, repeatedly bootstrap
//sample from mydat until the sum of the results exceeds ...
//requires 7 draws. Then 7 is a random sample of size 1 f...
//distribution I'm interested in. I want to repeat this p...
//"large" number, nsim, of times to estimate the distribu...
先に進む前に、R の初心者にどうしたら、もしくはどうすれば...
//Before reading on, I invite R novices to consider how -...
//write code -- to do this.
最初のアプローチ（素朴／直接的）：入れ子になったループを...
内側のループは繰返し mydat からサイズ 1 の標本を抽出し、...
//Approach 1 (Naive/ straightforward): Perform the simula...
//described by using a nested loop. The outer loop goes f...
//record the number of draws needed. The inner loop repea...
//sample of size 1 from mydat and accumulates and checks ...
//exceeds K, returning the number of times it took to do ...
//loop.
これは内側のループに "while" ループを使えば簡単にプログラ...
//This is easy to program up using a "while" loop for the...
//takes about a short look's time out my window (a few se...
//modest nsim (5-10K). But surely one can do better...
アプローチ２：すべての抽出を「一度に」やる以外は同じこと...
これは各繰返し毎に関数呼び出しをするという大きなオーバー...
//Approach 2: Was essentially the same, except all the sa...
//"at once." That is, each execution of the inner loop in...
//to call the "sample()" function to get a bootstrap samp...
//This is very inefficient, because there is a high overh...
//that function call at each iteration.
sample() をたった一度だけ、もしくは最大でもわずかな回数、...
//It can be avoided by calling
//sample() only once -- or at most a small number of time...
//large sample and then using indexing to work one's way ...
//sample.
これは、「大きな標本」とはどのくらいか見積り、添字操作を...
//This turns out to be a little more complicated than app...
//because you first have to figure how large a "large" sa...
//and then have to be a little careful about setting up y...
しかし、（特にポインタ操作を行なえる C プログラマにとって...
//But it's not very difficult (especially for any C progr...
//manipulate pointers), and random sampling is so fast th...
//problem to generate a sample WAYYY bigger than you need...
//as large, even) just to be safe. Or generate a not so b...
//just check at each outer loop iteration whether you nee...
//some more.
このように、内側のループで関数呼び出しでなく添字操作を行...
（5 から 10 倍）が得られました。分布を生成するのに必要な...
これはいかなる状況下でも、私の必要に対して十分過ぎるもの...
//Doing it this way -- now using indexing, not function c...
//loop -- made a considerable improvement (5 - 10 fold). ...
//one quick glance out the window time to generate the di...
//than a second). This was more than adequate for my need...
//conceivable situation.
しかしそれでもなお、自尊心のある R プログラマはループを避...
//But still, self-respecting R programmers are
//supposed to eschew looping, and here was this ugly(?) l...
//loop! So let us persevere.
アプローチ３：したがって、問題は、如何にして内側のループ...
//Approach 3: So the question is -- how to remove that in...
//I invite novice R programmers to think about that befor...
//on...
ここでの核心は、内側のループが単なる累積和という非常に簡...
なるほど、回答はすぐそこです。R は既にこれを行なう（背後...
//The key here is that the inner loop is doing a very sim...
//just accumulating a sum. Forsooth! -- the answer fairl...
//already has a cumsum() function that does this (looping...
//underlying C code), so one only has to figure out how t...
これを行なう本質的なアイデアは、大きな mydat 標本から、総...
総和が K を越えるには平均で約 10 個分の mydat データが必...
//The essential idea to do this is to get a small sample ...
//mydat sample that is guaranteed to be bigger than one n...
//to more than K. Again, one can be profligate: if I need...
//values on average to sum up to K, if I sample 30 or 40 ...
//thing, I'm sure to have enough (and can always add a ch...
//"try()" to make sure if I am faint of heart). If smalls...
//sample of the next 30 or 40 values from my large bootst...
//the following code vectorizes the inner loops:
count <- sum(cumsum(smallsamp) < K) + 1
ここで使われているトリックに注意しましょう。「< K」は K ...
//Note a trick here: The <K produces a vector of logical ...
//cumsum values <K. This is silently treated as numeric 1...
//sum(). A moment's thought will make clear why 1 must be...
//end..
実際、このベクトル化は、私の適当な大きさの nsim に対し、...
//Sure enough, this vectorization reduced the time about...
//-- to an eyeblink -- for my modest nsim sizes.
アプローチ４：最後の段階は、もちろん、毎回分布ベクトルに...
//Approach 4: The last stage, of course, is to get rid of...
//which fills the vector for the distribution one element...
//another R no-no. One way to do this -- the ONLY way I c...
//a challenge to smarter R programmers) -- is to make use...
//which basically always can remove loops. Unfortunately,...
//illusion, because what the apply() functions actually d...
//loops in clever R code, not remove them.
ここで理解できるように、そのメリットは、添字の管理を自分...
//Their virtue, as will be seen
//here, is making the code more transparent by transferri...
//bookkeeping of indexing to the R functions rather than ...
//explicitly handle it yourself.
apply() を使うための簡単なアイデアは、行数が nsim で、十...
//To use apply, the simple idea is just to make up a matr...
//columns with enough rows so that each column contains "...
//draws from mydat to sum to more than K. Let's say that ...
//will do it. Then the code for the whole simulation is:
drawmat <- matrix( sample(mydat, M*nsim, rep = TRUE), nc...
# (注）mydat から M*nsim 回復元抽出し、それを列数 nsim ...
# 行列 drawmat に付値
drawmat <- apply(drawmat, 2, cumsum) < K
# （注）apply 関数で行毎（apply 関数の引数 2）に累積和を...
# さらに K 未満かどうかで TRUE, FALSE に置き換える...
result <- colSums(drawmat) + 1
# （注）colSum 関数で行毎の総和を計算し、1 を加えて、ベ...
最初の行はブートストラップ標本を生成します。二・三番めは...
//The first statement generates the matrix of bootstrap s...
//second and third use the same trick as previously to ge...
//rows needed for each column to total more than K. Note ...
//handy (and fast) internal colSums function to add up th...
//columns. This code is much shorter, and I would say mor...
//then the code produced by the previous explicit looping...
しかしながら、驚くには当たりませんが、これはアプローチ３...
//However, not surprisingly, it does not run any faster (...
//the Approach 3 code.The looping is still there hidden i...
この長ったらしい例が、 とっつきにくいかのように思える言語...
//I hope that this long-winded example will help R newcom...
//entry into what may seem like a daunting programming la...
//have, I would urge all to consult R's quite remarkable ...
//Help facility (BEFORE they post to this list), as well ...
//Ripley's two indispensable books on S language programm...
//analysis.
Cheers,~
Bert Gunter~
~
Non-Clinical Biostatistics~
Genentech~
"統計家の仕事は、科学における学習過程への触媒となることで...
訳注：上の三種類のアプローチによるコード例と実行時間（投...
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,"-t...
result[i] <<- sum((cumsum(x) < K)) + 1
}
}
test3 <- function() {
drawmat <- matrix( sample(mydat, M*nsim, rep = TRUE), ...
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)) < ...
}
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コード最適化のコツと実例集]]を参照。 -- &new{2004-06...
#comment
ページ名:
![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
