2

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

m=10 我得到 0.17 秒

m=100 我得到 3.95 秒

m=1000 我得到 246.26 秒

m=2000 我得到 1003.55 秒

我希望将采样次数 m 设置为 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
        }
   }
 )
4

2 回答 2

3

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

  • 采样一个长的随机数向量（然后将其折叠成一个矩阵）比重复采样短向量要快得多（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")

在此处输入图像描述

于 2012-10-22T23:20:35.880 回答
0

一般来说，您希望尽可能利用矢量化，而不是为了速度，而是为了提高可读性/理解力。

power[n]为什么要在内循环内部写入（我猜也是计算vecpals）？内循环执行后不应该在外循环中吗？您可能希望将平方根的计算移到两个循环之外。

为什么teststatvect和被power初始化为矩阵（明确地是二维数组）而不是向量（或者更确切地说，作为一维数组，使用array）？只是variance at 225上一行注释的结尾吗？您可能需要检查格式。（这是作业吗？）

对于您在这里尝试做的事情，您可能希望利用非常方便的函数replicate，也许通过编写一个特定的函数来调用它。

于 2012-10-22T22:53:36.147 回答