題意に沿えばよいだけだろう

「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

当然，速いに決まっているわ～な(^_^)

最適化なんかしなくてよい その5

これも同じく，本家 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

どれでも同じ。遠いってどういうこと?
使いたい機能，使わなければならない機能は，何の遠慮もなく使うがよいゾッ!

よって，この珍説もガセと認定!!

最適化なんかしなくてよい その4

これも同じく，本家 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 にとっては痛くもかゆくもない。

よって，この珍説もガセと認定!!

最適化なんかしなくてよい その3

本家 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秒にも満たないのだから，あれこれ小細工を弄する必要は全くないと言っていいだろう。