使 R 中的嵌套 for 循环更高效
Making nested for loops in R more efficient
我正在进行一个研究项目，我想确定两个分布的等价性。我目前正在使用 Mann-Whitney 等价性检验，Stefan Wellek (2010) 的《等价性和非劣效性检验统计假设》一书提供了我 运行ning 的代码（如下）。在 运行 收集我的数据之前，我使用具有相同均值和标准差的随机正态分布测试此代码。我的问题是存在三个嵌套的 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)
您可以使用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
其他循环比较棘手。
请参阅下面的编辑以获得更好的建议
一个提高速度的简单建议是 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 and the inline 包了吗？
我一直很想知道如何使用它们，所以我认为这是一个很好的机会。
这是使用 inline 包和 Fort运行 对您的代码进行的调整（因为我对它比 C 更熟悉）。这一点都不难（前提是你知道 Fort运行 或 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
我确实检查过它得到了与您的代码相同的结果，但您可能想要 运行 一些额外的测试以防万一。
我正在进行一个研究项目，我想确定两个分布的等价性。我目前正在使用 Mann-Whitney 等价性检验，Stefan Wellek (2010) 的《等价性和非劣效性检验统计假设》一书提供了我 运行ning 的代码（如下）。在 运行 收集我的数据之前，我使用具有相同均值和标准差的随机正态分布测试此代码。我的问题是存在三个嵌套的 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)
您可以使用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
其他循环比较棘手。
请参阅下面的编辑以获得更好的建议
一个提高速度的简单建议是 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 and the inline 包了吗？
我一直很想知道如何使用它们，所以我认为这是一个很好的机会。
这是使用 inline 包和 Fort运行 对您的代码进行的调整（因为我对它比 C 更熟悉）。这一点都不难（前提是你知道 Fort运行 或 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
我确实检查过它得到了与您的代码相同的结果，但您可能想要 运行 一些额外的测试以防万一。