Создание вложенных циклов в R более эффективно

Question

Я работаю над исследовательским проектом, где я хочу определить эквивалентность двух распределений. В настоящее время я использую тест Mann-Whitney для эквивалентности, а код, который я запускаю (ниже), был предоставлен в книге «Тестирование статистических гипотез эквивалентности и неинфекции» Стефана Веллека (Stefan Wellek, 2010). Перед запуском моих данных я тестирую этот код со случайными нормальными распределениями, которые имеют одинаковое среднее и стандартное отклонение. Моя проблема заключается в том, что существует три вложенных цикла и при работе с большими размерами распределений (как в примере ниже) код выполняется навсегда для запуска. Если бы мне только пришлось запускать его один раз, это не было бы такой проблемой, но я выполняю симуляционный тест и создаю кривые мощности, поэтому мне нужно запустить много итераций этого кода (около 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 и запустил ее на моем ноутбуке. Просто байт, компилирующий ваш текущий код, он работает в два раза быстрее.

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 

Вы получаете примерно 50x:

ускорение

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 

Остальные петли сложнее. для еще лучшего предложения