使 R 中的嵌套 for 循环更有效

Question

我正在从事一个研究项目，我想确定两个分布的等价性。我目前正在使用 Mann-Whitney Test for Equivalence，我正在运行的代码（如下）由 Stefan Wellek（2010 年）的《Testing Statistical Hypotheses of Equivalence and Noninferiority》一书提供。在运行我的数据之前，我正在使用具有相同均值和标准差的随机正态分布测试此代码。我的问题是有三个嵌套的 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 行开始的函数中，并在我的笔记本电脑上运行它。只需对当前代码进行字节编译，它的运行速度就会提高一倍。

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 

编辑

我应该补充一下，这是一种不同的语言会比 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

其他循环比较棘手。