在 R 中模拟 - 我怎样才能让它更快？

Question

我正在模拟类似Jim Berger's applet 的东西。

模拟工作如下：我将从零分布 N(0,1) 或替代分布 N( θ，1）。我将假设空值的先验概率是某个比例prop（因此替代方案的先验是1-prop）并且theta 在替代方案中的分布是N(0,2) （我可以更改所有这些参数，但这只是开始）。

我想从上述模拟场景中获得一定范围内的大量 pvalue（例如 0.049 和 0.05 之间的 2000 个 pvalue，在模拟中这将相当于 z stats arround 1.96 和 1.97），并查看有多少来自 null，有多少来自替代。

到目前为止，我想出了一个这样的解决方案：

berger <- function(prop, n){
  z=0
  while(z<=1.96|z>=1.97){
    u <- runif(1)
    if(u<prop){
      H0 <- TRUE
      x<-rnorm(n, 0, 1)
    }else{
      H0 <- FALSE
      theta <- rnorm(1, 0, 2)
      x <- rnorm(n, theta, 1)
    }
    z <- sqrt(n)*abs(mean(x))
  }
  return(H0)
}

results<-replicate(2000, berger(0.1, 100))
sum(results)/length(results) ## approximately 25%

大约需要 3.5 分钟。有可能加快这个速度吗？如何？欢迎所有答案，包括与 C 的集成。

更新：并行化可以稍微加快速度。但是，我在 Julia 中尝试过相同的代码，并且在没有任何并行化的情况下只需要 14 秒（代码如下）。

更新 2：使用 Rcpp 和并行化可以将模拟时间缩短到 8 秒。查看新答案。

function berger(prop, n)
       z = 0 
       h0 = 0
       while z<1.96 || z > 1.97

              u = rand()

              if u < prop
                     h0 = true;
                     x = randn(n)             
              else
                     h0 = false
                     theta = randn()*2
                     x = randn(n) + theta
              end

              z = sqrt(n)*abs(mean(x))
       end

       h0
end

results = [0]

for i in 1:2000
       push!(results, berger(0.1, 100))
end

sum(results)/length(results)

Answer 1

可能有一些方法可以让这个函数更快一点（例如通过并行化），但你不会得到数量级的差异（edit：in R em>）。关键问题是您从正态分布中抽取了大约 4 亿次。

这是一个函数，它返回您的函数通过 while 的平均运行次数：

f<-function(prop,n){
  i<-0
  z<-0
  while(z<=1.96|z>=1.97){
    i<-i+1
    u <- runif(1)
    if(u<prop){
      H0 <- TRUE
      x<-rnorm(n, 0, 1)
    }else{
      H0 <- FALSE
      theta <- rnorm(1, 0, 2)
      x <- rnorm(n, theta, 1)
    }
    z <- sqrt(n)*abs(mean(x))
  }
  return(i)
}

现在我们可以计算你的函数运行了多少次：

set.seed(1)
runs<-replicate(200,f(prop=0.1, n=100))
mean(runs) # 2034
sd(runs) # 2121

所以，从正态分布中计算抽奖次数：

# number of replicates
# times normal distributions per replicate
# draws from each distribution
2000*mean(runs)*100
# 406,853,000 normal distribution draws

rnorm 函数调用已编译的 C 函数，并且可能接近最佳速度。您可以测试在您自己的机器上进行这么多抽奖的“下限”：

system.time(rnorm(406853000))
# My machine:
#   user  system elapsed 
#  53.78    2.39   56.62 

相比之下，您的函数运行速度大约慢了四倍：

system.time(replicate(2000,berger(prop=0.1,n=100)))
#    user  system elapsed 
#  210.40    0.03  211.12 

因此，您的函数实际上并没有那么慢，尤其是当您考虑到每次调用 rnorm 都会产生开销时。如果你提高这个函数的速度非常关键，并且你有几个内核，你可以很容易地在 R 中并行化它：

library(parallel)
mclapply(1:2000,function(x) berger(prop=0.1,n=100))

除此之外，您可以用 C 编写一个超级优化的函数并节省几分钟，但这可能不值得。

Answer 2

使用 Rcpp 来加快速度实际上很简单。将 Rcpp 与并行化相结合，我能够将时间缩短到 8 秒。

.cpp 文件是这样的（使用 Rcpp "sugars" 使这项任务变得非常容易 - 因为这是我第一次使用 Rcpp，也许这段代码不是最佳的，但它完成了工作！）：

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]


int RcppBerger(double prop, int n) {

  double z=0,theta=0, u=0;
  int h = 0;
  NumericVector x;
    while (z<1.96 || z > 1.97){
      u = R::runif(0, 1);
      if(u < prop){
        h = 1;
        x = rnorm(n);
        }else{
          h = 0;
          theta = R::rnorm(0, 2);
          x = rnorm(n, theta, 1);
          }
          z = sqrt(n)*mean(x);
          if(z<0){z = -1*z;};
    }
  return h;
}

然后，在没有并行化的情况下，您可以只使用sourceCpp 函数，RcppBerger 将在工作区中可用：

library(Rcpp)
sourceCpp("RcppBerger.cpp")
results<-replicate(2000, RcppBerger(0.1, 100))
sum(results)/length(results) ## approximately 25%

这已经将时间从 3.5 分钟缩短到 40 秒左右。之后我们就可以并行化了。

在Windows中，这有点棘手，看来您必须先创建一个包。但是 Rcpp 提供了一个很好的功能来做到这一点Rcpp.package.skeleton。只需将源文件放入其中，它将创建所有必要的文档和文件夹：

Rcpp.package.skeleton("RcppBerger", cpp_files = "RcppBerger.cpp")

然后，安装包后，就可以与foreach和doParallel并行了：

library(foreach)
library(doParallel)
library(RcppBerger)
registerDoParallel(cores=8)
results<- foreach(1:2000, .packages="RcppBerger") %dopar% RcppBerger(0.1, 100)

现在模拟只需要 8 秒。