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

Question

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

模拟的工作方式如下:我将从空分布 N(0,1) 中生成大小为 n 的样本 x或来自替代分布 N(theta, 1)。我将假设 null 的先验概率是某个比例 prop(因此备选方案的先验是 1-prop)并且 theta 的分布 替代方案是 N(0,2)(我可以更改所有这些参数，但这只是开始)。

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

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

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

可能有一些方法可以让这个函数更快一点(例如通过并行化)，但你不会得到数量级的差异(编辑:在 R 中).关键问题是你从正态分布中提取了大约 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 编写一个 super 优化的函数并节省几分钟，但这可能不值得。