Rでネストされたforループをより効率的に作る

Question

私は2つのディストリビューションの等価性を決定したい研究プロジェクトに取り組んでいます。私は現在、Mann-Whitney Test for Equivalenceを使用しています。実行しているコードは、Stefan Wellek（2010）の「等価性と非劣性の統計的仮説検定」の本で提供されています。私のデータを実行する前に、私は同じ平均と標準偏差を持つランダムな正規分布でこのコードをテストしています。私の問題は、3つのネストされたforループがあり、より大きなディストリビューションサイズ（以下の例のように）を実行するとコードが永遠に実行されることです。一度だけ実行しなければならないのであれば、それほど問題はありませんが、私はシミュレーションテストを行い、パワーカーブを作成していますので、このコードを繰り返し実行する必要があります（約10,000）。現時点では、配布サイズをどのように変更するかによって、10,000回の繰り返しを実行するには数日かかります。

これのパフォーマンスを向上させるための助けがあれば、大歓迎です。

x <- rnorm(n=125, m=3, sd=1) 
y <- rnorm(n=500, m=3, sd=1) 

alpha <- 0.05 
m <- length(x) 
n <- length(y) 
eps1_ <- 0.2 #0.1382 default 
eps2_ <- 0.2 #0.2602 default 

eqctr <- 0.5 + (eps2_-eps1_)/2 
eqleng <- eps1_ + eps2_ 

wxy <- 0 
pihxxy <- 0 
pihxyy <- 0 

for (i in 1:m) 
for (j in 1:n) 
    wxy <- wxy + trunc(0.5*(sign(x[i] - y[j]) + 1)) 

for (i in 1:m) 
for (j1 in 1:(n-1)) 
    for (j2 in (j1+1):n) 
    pihxyy <- pihxyy + trunc(0.5*(sign(x[i] - max(y[j1],y[j2])) + 1)) 

for (i1 in 1:(m-1)) 
for (i2 in (i1+1):m) 
    for (j in 1:n) 
    pihxxy <- pihxxy + trunc(0.5*(sign(min(x[i1],x[i2]) - y[j]) + 1)) 

wxy <- wxy/(m*n) 
pihxxy <- pihxxy*2/(m*(m-1)*n) 
pihxyy <- pihxyy*2/(n*(n-1)*m) 
sigmah <- sqrt((wxy-(m+n-1)*wxy**2+(m-1)*pihxxy+(n-1)*pihxyy)/(m*n)) 

crit <- sqrt(qchisq(alpha,1,(eqleng/2/sigmah)**2)) 

if (abs((wxy-eqctr)/sigmah) >= crit) rej <- 1 
if (abs((wxy-eqctr)/sigmah) < crit) rej <- 0 

if (is.na(sigmah) || is.na(crit)) rej <- 1 

MW_Decision <- rej 

cat(" ALPHA =",alpha," M =",m," N =",n," EPS1_ =",eps1_," EPS2_ =",eps2_, 
    "\n","WXY =",wxy," SIGMAH =",sigmah," CRIT =",crit," REJ=",MW_Decision) 

Answer 1

は、スピードブーストのビットを取得するための一つの簡単な提案はあなたのコードをbyte compileすることです

さらに優れた提案については、以下の編集を参照してください。

たとえば、コードをalpha <- 0.05行から始まる関数にラップし、ラップトップで実行しました。現在のコードを単純にバイトコンパイルすると、2倍の速さで実行されます。私は追加する必要があります

set.seed(1234) 
x <- rnorm(n=125, m=3, sd=1) 
y <- rnorm(n=500, m=3, sd=1) 

# f1 <- function(x,y){ ...your code...} 

system.time(f1(x, y)) 
# user system elapsed 
# 33.249 0.008 33.278 

library(compiler) 
f2 <- cmpfun(f1) 

system.time(f2(x, y)) 

# user system elapsed 
# 17.162 0.002 17.170 

EDIT

、これは、異なる言語がRよりもはるかに良いのだろう物事の種類は、あなたがRcppとinlineパッケージを見てきましたでしょうか？

私はこれを使う方法を知りたいので、これは良いチャンスだと思いました。

inlineパッケージとFortranを使用したコードの調整があります（これはCよりも快適です）。 FortranやCを知っていれば、それほど難しいことではありませんでした。私はちょうどcfunctionにリストされた例を辿った。

まずは、あなたのループを再書き込みし、それらをコンパイルしてみましょう：

library(inline) 

# Fortran code for first loop 
loop1code <- " 
    integer i, j1, j2 
    real*8 tmp 
    do i = 1, m 
     do j1 = 1, n-1 
     do j2 = j1+1, n 
      tmp = x(i) - max(y(j1),y(j2)) 
      if (tmp > 0.) pihxyy = pihxyy + 1 
     end do 
     end do 
    end do 
"  
# Compile the code and turn loop into a function 
loop1fun <- cfunction(sig = signature(x="numeric", y="numeric", pihxyy="integer", m="integer", n="integer"), dim=c("(m)", "(n)", "", "", ""), loop1code, language="F95") 

# Fortran code for second loop 
loop2code <- " 
    integer i1, i2, j 
    real*8 tmp 
    do i1 = 1, m-1 
     do i2 = i1+1, m 
     do j = 1, n 
      tmp = min(x(i1), x(i2)) - y(j) 
      if (tmp > 0.) pihxxy = pihxxy + 1 
     end do 
     end do 
    end do 
"  
# Compile the code and turn loop into a function 
loop2fun <- cfunction(sig = signature(x="numeric", y="numeric", pihxxy="integer", m="integer", n="integer"), dim=c("(m)", "(n)", "", "", ""), loop2code, language="F95") 

それでは、これらを使用する新しい関数を作成してみましょう。だから、それはあまりにも長くはないですが、私はちょうど私があなたのコードから変更キーパーツスケッチます：

f3 <- function(x, y){ 

    # ... code ... 

# Remove old loop 
## for (i in 1:m) 
## for (j1 in 1:(n-1)) 
## for (j2 in (j1+1):n) 
##  pihxyy <- pihxyy + trunc(0.5*(sign(x[i] - max(y[j1],y[j2])) + 1)) 

# Call new function from compiled code instead 
pihxyy <- loop1fun(x, y, pihxyy, m, n)$pihxyy 

# Remove second loop 
## for (i1 in 1:(m-1)) 
## for (i2 in (i1+1):m) 
## for (j in 1:n) 
##  pihxxy <- pihxxy + trunc(0.5*(sign(min(x[i1],x[i2]) - y[j]) + 1)) 

# Call new compiled function for second loop 
pihxxy <- loop2fun(x, y, pihxxy, m, n)$pihxxy 

# ... code ... 
} 

をそして今、我々はそれを実行して、出来上がり、私たちは巨大なスピードブーストを取得します！:)

system.time(f3(x, y)) 
# user system elapsed 
    0.12 0.00 0.12 

あなたのコードと同じ結果が得られたことを確認しましたが、場合によっては追加のテストを実行したいと思うかもしれません。

Answer 2

代わりに最初の二重ループのouterを使用することができます。

set.seed(42) 

f1 <- function(x,y) { 
wxy <- 0 
for (i in 1:m) 
    for (j in 1:n) 
    wxy <- wxy + trunc(0.5*(sign(x[i] - y[j]) + 1)) 
wxy 
} 

f2 <- function(x,y) sum(outer(x,y, function(x,y) trunc(0.5*(sign(x-y)+1)))) 

f1(x,y) 
[1] 32041 
f2(x,y) 
[1] 32041 

あなたはおよそ50倍のスピードアップを取得：

library(microbenchmark) 
microbenchmark(f1(x,y),f2(x,y)) 
Unit: milliseconds 
    expr  min   lq  median   uq  max neval 
f1(x, y) 138.223841 142.586559 143.642650 145.754241 183.0024 100 
f2(x, y) 1.846927 2.194879 2.677827 3.141236 21.1463 100 

他のループはトリッキーです。