r - 嵌套 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

(我刚刚尝试了这个，在从 1 到 200 的每个样本量值上采样 m=1e5 ......这需要 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")