嵌套 for 循环的效率

Question

我创建了以下代码，将 for 循环嵌套在 R 中的 for 循环内。它是计算 Power 的模拟。我读到 R 不太适合执行 for 循环，但我想知道是否可以应用任何效率来使其运行得更快一些。我对 R 以及任何类型的编程都很陌生。现在我看到的运行时间是：

m=10 我得到 0.17 秒

m=100 我得到 3.95 秒

m=1000 我得到 246.26 秒

m=2000 我得到 1003.55 秒

我希望将采样次数设置为 100K 以上，但恐怕甚至无法将其设置为 10K

这是代码：

m = 1000                        # number of times we are going to  take samples
popmean=120                     # set population mean at 120
popvar=225                      # set known/established population 
variance at 225
newvar=144                      # variance of new methodology 
alpha=.01                       # set alpha
teststatvect = matrix(nrow=m,ncol=1)    # empty vector to populate with test statistics
power = matrix(nrow=200,ncol=1)     # empty vector to populate with power

system.time(                    # not needed - using to gauge how long this takes
    for (n in 1:length(power))          # begin for loop for different sample sizes
      for(i in 1:m){                # begin for loop to take "m" samples
      y=rnorm(n,popmean,sqrt(newvar))   # sample of size n with mean 120 and var=144
      ts=sum((y-popmean)^2/popvar)      # calculate test statistic for each sample
      teststatvect[i]=ts            # loop and populate the vector to hold test statistics
      vecpvals=pchisq(teststatvect,n)   # calculate the pval of each statistic
      power[n]=length(which(vecpvals<=alpha))/length(vecpvals) # loop to populate      power vector. Power is the proportion lessthan ot equal to alpha
        }
   }
 )

Answer 1

我稍微重新组织了你的代码并摆脱了内部循环。

• 对一个长随机数向量进行采样（然后将其折叠成一个矩阵）比重复对短向量进行采样要快得多（replicate正如另一个答案中所建议的那样，这对于可读性很好，但在这种情况下，您可以通过对以下中的随机数进行采样来做得更好一个块）
• colSums比for循环内求和或使用apply.
• 它只是糖（即它实际上并没有更有效），但你可以mean(pvals<=alpha)使用sum(pvals<=alpha)/length(alpha)
• 我定义了一个函数来返回一组指定参数（包括样本大小）的幂，然后用于sapply范围化大小向量（不比for循环快，但更干净并且可能更容易泛化）。

代码：

powfun <- function(ssize=100,
                   m=1000,      ## samples per trial
                   popmean=120, ## pop mean
                   popvar=225,  ## known/established pop variance
                   newvar=144,  ## variance of new methodology
                   alpha=0.01,
                   sampchisq=FALSE)  ## sample directly from chi-squared distrib?
{
    if (!sampchisq) {
      ymat <- matrix(rnorm(ssize*m,popmean,sd=sqrt(newvar)),ncol=m)
      ts <- colSums((ymat-popmean)^2/popvar)          ## test statistic
    } else {
      ts <- rchisq(m,df=ssize)*newvar/popvar
    }
    pvals <- pchisq(ts,df=ssize)                    ## pval
    mean(pvals<=alpha)                              ## power
}

您是否真的需要样本大小的每个整数值的功效，或者间隔更宽的样本就可以（如果您需要精确的值，插值可能会非常准确）

ssizevec <- seq(10,250,by=5)
set.seed(101)
system.time(powvec <- sapply(ssizevec,powfun,m=5000))  ## 13 secs elapsed

这是相当快的，如果你需要的话，可能会让你达到目标m=1e5，但我不太清楚为什么你需要那么精确的结果——功率曲线相当平滑m=5000......

sapply(ssizevec,powfun,m=5000)如果您不耐烦地等待长时间的模拟，您还可以通过替换为来打印进度条library(plyr); aaply(ssizevec,.margins=1,powfun,.progress="text",m=5000)

最后，我认为您可以通过直接采样卡方值或通过进行分析功效计算（！）来大大加快整体速度。我认为这rchisq(m,df=ssize)*newvar/popvar相当于循环的前两行，您甚至可以直接对卡方密度进行数值计算......

system.time(powvec2 <- sapply(ssizevec,powfun,m=5000,sampchisq=TRUE))
## 0.24 seconds elapsed

（我刚刚尝试过，m=1e5对从 1 到 200 的样本大小的每个值进行采样......需要 24 秒......但我仍然认为这可能是不必要的。）

照片：

par(bty="l",las=1)
plot(ssizevec,powvec,type="l",xlab="sample size",ylab="power",
     xlim=c(0,250),ylim=c(0,1))
lines(ssizevec,powvec2,col="red")