Verschachtelte For-Schleifen in R effizienter machen

Question

Ich arbeite an einem Forschungsprojekt, bei dem ich die Äquivalenz zweier Verteilungen bestimmen möchte. Ich verwende derzeit den Mann-Whitney-Test für Äquivalenz, und der Code, den ich führe (unten), wurde mit dem Buch Testing Statistical Hypotheses of Equivalence and Noninferiority von Stefan Wellek (2010) bereitgestellt. Vor dem Ausführen meiner Daten teste ich diesen Code mit zufälligen Normalverteilungen, die den gleichen Mittelwert und die Standardabweichung haben. Mein Problem ist, dass es drei verschachtelte For-Schleifen gibt und wenn größere Distributionsgrößen ausgeführt werden (wie im Beispiel unten), dauert es ewig, bis der Code ausgeführt wird. Wenn ich es nur einmal ausführen müsste, wäre das kein Problem, aber ich mache einen Simulationstest und erstelle Leistungskurven, so dass ich viele Iterationen dieses Codes (etwa 10.000) ausführen muss. Im Moment dauert es, je nachdem, wie ich die Verteilungsgrößen verändere, bis zu 10.000 Iterationen.

Jede Hilfe in einer Weise, die Leistung von diesem zu erhöhen, würde sehr geschätzt werden.

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

Siehe bearbeiten unten für einen noch besseren Vorschlag

Ein einfacher Vorschlag ein bisschen ein Geschwindigkeitsschub zu bekommen, ist der Code byte compile.

Zum Beispiel habe ich Ihren Code in eine Funktion von der alpha <- 0.05 Linie gewickelt und lief es auf meinem Laptop. Einfaches Byte Kompilieren Ihres aktuellen Codes, es läuft doppelt so schnell.

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

ich hinzufügen sollte, das ist die Art von Dingen, die eine andere Sprache als R. viel besser tun würde, Haben Sie an der Rcpp und die inline Pakete geschaut?

Ich war neugierig zu lernen, wie man sie verwendet, also dachte ich, das war eine gute Chance.

Hier ist ein Tweak Ihres Codes mit dem inline Paket und Fortran (da ich mich damit wohler als C bin). Es war überhaupt nicht schwer (vorausgesetzt, Sie kennen Fortran oder C); Ich folgte nur den Beispielen, die in cfunction aufgeführt sind.

Lassen Sie uns zunächst Ihre Loops neu schreiben und sie kompilieren:

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") 

Lassen Sie uns jetzt eine neue Funktion erstellen, die diese verwendet. Es ist also nicht zu lange, werde ich skizziere nur die wichtigsten Teile, die ich aus dem Code geändert:

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 ... 
} 

Und jetzt haben wir es laufen und voila, wir bekommen einen riesiger Geschwindigkeitsschub!:)

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

Ich habe überprüfen, ob sie die gleichen Ergebnisse wie Sie Ihren Code erhalten, aber Sie wollen wahrscheinlich einige zusätzliche Tests nur für den Fall laufen.

Answer 2

können Sie verwenden outer anstelle des ersten Doppel-Loop:

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 

Sie erhalten etwa 50x Speedup:

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 

Die anderen Schleifen sind kniffliger.