「Rコードの最適化例:混合正規乱数の発生コード」 にあるんだが,なんだかなあ
結局 N1 に従う乱数の個数 n1 が決まれば良いと喝破(なるほど)
rmixnorm3 <- function(n, p1, N1, N2) {
n1 <- sum(runif(n) < p1)
c( N1[2]*rnorm(n1) + N1[1], N2[2]*rnorm(n-n1)+N2[1] )
}
rmixnorm4 <- function(n, p1, N1, N2) {
n1 <- sum(runif(n) < p1)
c( rnorm(n1, m=N1[1], s=N1[2]), rnorm(n-n1, m=N2[1], s=N2[2]) )
}
何が喝破か。
「確率 0.6 で混合」と書いてあるのに,何が悲しくて「その個数を決める」のに,sum(runif(n) < p1) なんかやるの!!
n*p1 でよいでしょうが(必要なら整数に丸めないといけないけど)。ということで,rmixnorm5 を定義して,ベンチマーク。
rmixnorm5 <- function(n, p1, N1, N2) {
n1 <- round(n*p1)
c( rnorm(n1, m=N1[1], s=N1[2]), rnorm(n-n1, m=N2[1], s=N2[2]) )
}
library(rbenchmark)
benchmark(rmixnorm(1e6, 0.4, c(-1,1), c(3,2)), rmixnorm2(1e6, 0.4, c(-1,1), c(3,2)), rmixnorm3(1e6, 0.4, c(-1,1), c(3,2)), rmixnorm4(1e6, 0.4, c(-1,1), c(3,2)), rmixnorm5(1e6, 0.4, c(-1,1), c(3,2)), replications=100)
test replications elapsed relative user.self sys.self
1 rmixnorm(1e+06, 0.4, c(-1, 1), c(3, 2)) 100 19.705 2.039644 17.734 2.418
2 rmixnorm2(1e+06, 0.4, c(-1, 1), c(3, 2)) 100 22.562 2.335369 20.642 2.233
3 rmixnorm3(1e+06, 0.4, c(-1, 1), c(3, 2)) 100 13.959 1.444881 13.148 0.945
4 rmixnorm4(1e+06, 0.4, c(-1, 1), c(3, 2)) 100 13.103 1.356278 12.720 0.467
5 rmixnorm5(1e+06, 0.4, c(-1, 1), c(3, 2)) 100 9.661 1.000000 9.333 0.392
当然,速いに決まっているわ~な(^_^)
これも同じく,本家 Rjpwiki の「Rコード最適化のコツと実例集」にあるのだけど。。。みんな~,騙されるんじゃないぞ~!
> 永続付値(他の環境中の変数への付値)は手間がかかる。遠くの x よりも近くの x
>
> > test1 <- function(){x <- runif(10000)} # 自前の環境中の x へ代入
> > test2 <- function(){x <<- runif(10000)} # 親環境中の x へ代入
> > rm(x); c = 0; for(i in 1:100) c <- c + system.time(test1()); c/100
> [1] 0.0015 0.0000 0.0020 0.0000 0.0000
> > x
> Error: Object "x" not found # 変数 x は関数終了時に消える
> > x <- numeric(10000); c = 0; for(i in 1:100) c <- c + system.time(test2()); c/100
> [1] 0.0018 0.0000 0.0021 0.0000 0.0000
> > x[1:5] # 当然 x は存在
> [1] 0.02132324 0.36266310 0.50504673 0.89662908 0.67079898
> [1] 0.48 0.03 0.50 0.00 0.00
はいはい。やってみましょうね~~。
> test1 <- function() {x <- runif(1000000)} # 自前の環境中の x へ代入
> test2 <- function() {x <<- runif(1000000)} # 親環境中の x へ代入
> library(rbenchmark)
> benchmark(test1(), test2(), replications=1000)
test replications elapsed relative user.self sys.self user.child sys.child
1 test1() 1000 28.783 1.000000 27.833 1.282 0 0
2 test2() 1000 29.294 1.017754 28.239 1.256 0 0
どれでも同じ。遠いってどういうこと?
使いたい機能,使わなければならない機能は,何の遠慮もなく使うがよいゾッ!
よって,この珍説もガセと認定!!
これも同じく,本家 Rjpwiki の「Rコード最適化のコツと実例集」にあるのだけど。。。みんな~,騙されるんじゃあね~ぞ~!
> 既存変数を再利用(変数の構造・サイズが変わる場合はどうかはわかりません)。環境にやさしい資源リサイクル
>
> > test1 <- function() {x <- runif(1000000); y <- x^2} # 変数 y は無駄
> > system.time(test1())
> [1] 0.62 0.04 0.67 0.00 0.00
> > test2 <- function() {x <- runif(1000000); x <- x^2} # 変数 x を再利用
> > system.time(test2())
> [1] 0.48 0.03 0.50 0.00 0.00
この場合は,そもそも代入する必要がないのだから,以下も含めてベンチマーク・テストしてみよう
test3 <- function() {x <- runif(1000000); x^2}
> library(rbenchmark)
> benchmark(test1(), test2(), test3(), replications=1000)
test replications elapsed relative user.self sys.self user.child sys.child
1 test1() 1000 37.084 1.000000 30.681 6.671 0 0
2 test2() 1000 37.149 1.001753 30.557 6.778 0 0
3 test3() 1000 37.144 1.001618 30.558 6.804 0 0
どれでも同じ。「資源はたくさんある」ので,たかが 1000000 個の要素を持つ numeric 変数を新たに使っても,R にとっては痛くもかゆくもない。
よって,この珍説もガセと認定!!
本家 Rjpwiki の「Rコード最適化のコツと実例集」にあるのだけど。。。みんな~,騙されるんじゃないぞ~!
> 変数初期化のコツ。同一のベクトル、配列変数を全要素 0 で繰返し初期化する(良くあること)場合、既存の同じサイズの変数に 0 を掛ければ済む。x-x でも同じこと(少し早くなるとのコメント)。全要素 1 で初期化したければ、0*x + 1 とする。0 は無ならず
> 確かに,y <- 0*x は速いですね。y <- x-x の方がもう少しだけ速いようです。 -- 2004-06-23 (水) 10:41:05
> 1で初期化する場合は0*x+1よりx^0が高速。さらに論理値で良ければx==xがもっと高速 -- aaa? 2009-05-01 (金) 06:31:21
ガセ 全部,ウソッパチ
> x <- runif(100000000)
> gc();gc();system.time(y <- 0*x)
ユーザ システム 経過
0.223 0.216 0.437
> gc();gc();system.time(y <- x-x) # 0*x より速くなることもない
ユーザ システム 経過
0.229 0.220 0.447
奇をてらう必要はない。素直にやるのが一番速い。 <-- 四則演算も結構高くつく
> gc();gc();system.time(y <- numeric(length(x)))
ユーザ システム 経過
0.117 0.223 0.338
===============================================
> gc();gc();system.time(y <- 0*x+1)
ユーザ システム 経過
0.477 0.433 0.906
> gc();gc();system.time(y <- x^0) # 全然,速くない
ユーザ システム 経過
0.613 0.216 0.825
> gc();gc();system.time(y <- x==x) # 全然,速くない
ユーザ システム 経過
0.629 0.112 0.739
奇をてらう必要はない。素直にやるのが一番速い。
> gc();gc();system.time(y <- numeric(length(x))+1)
ユーザ システム 経過
0.370 0.443 0.809
昔は本当だったトリビアも,今ではガセになりはてているものもある。
マシンやコンパイラの違いによることもあるけど。
それにしても,毎度いうことだが,これだけの計算量で違いは1秒にも満たないのだから,あれこれ小細工を弄する必要は全くないと言っていいだろう。