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

Posted 2023-03-15

【中文标题】在 R 中模拟 - 我怎样才能让它更快？

【英文标题】：Simulating in R- how can I make this faster?

【发布时间】：2014-07-05 08:09:40

【问题描述】：

我正在模拟类似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)

【问题讨论】：

我真的不明白为什么只考虑0.049 < p < 0.05 很重要？我看到论文中提到了它，但对我个人来说没有意义。

@PascalvKooten 您想计算 H0 为真的概率，给定 p 约为 0.05。阅读 Jim Berger 的页面了解更多详情。

【参考方案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 编写一个超级优化的函数并节省几分钟，但这可能不值得。

【讨论】：

谢谢，真的很有趣，我会尝试并行化，看看效果如何！

虽然，我刚才运行了那个小程序，它运行得非常快，这让我觉得你的函数没有准确地反映他们做了什么。我强烈怀疑他们只从正态分布中提取了 2000 个样本，而不是 400 万个。

仅 2000 次模拟无法获得 0.049-0.05 范围内的 2000 个 p 值，它们必须运行更多。您可能使用 n=1 （默认值）运行，因此速度要快得多。当 n=100 时，applet 大约需要 1.5 分钟！我刚刚尝试了并行化，而 R 花了同样的时间！

您对模拟是正确的。关于并行化的三个问题：1）您是否使用options(mc.cores=x) 设置了核心数，其中x 是您机器上的核心数？ 2）你正在运行什么操作系统？如果是 Windows，则必须使用 foreach 而不是 parallel 3) 您是否监控 CPU 并确保它们在运行时都处于活动状态？

很高兴你在 Julia 上取得了成功。为了交叉检查，我在 C 中运行了一些代码来生成 4 亿个数字，发现只用了 3 秒。我怀疑你能否在 R 中达到这些速度。

【参考方案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 秒。

【讨论】：

写得很好，从sourceCpp() 到一个包的段落做得很好。您现在可以尝试更快的 N(0,1) 生成器：或者从 R 的默认值切换到更快的生成器，或者插入更快的 RNG，例如来自 RcppZiggurat 的一个。

在 RcppZiggurat 小插图中有一个图表——“Ahrens-Dieter”和“Kinderman-Ramage”都已经在 R 中，并且比默认反转快了可能两倍。 RcppZiggurat 可以是另一个两倍。让我们看看你能不能把这个从 8 秒缩短到 5 秒 :)

@CarlosCinelli 非常好。我受到这个问题的启发，想了解更多关于 N(0,1) 生成器的种类。我打算更新我的答案，但现在我看到你的答案，我很高兴我等了。小心RZiggurat 的实现；它在 Linux 机器上对我不起作用。另外，我认为您可能应该将复选标记切换到此答案，它比我的答案更完整且更有帮助。

谢谢@nograpes，我想是时候开始研究更多的数值方法和RNG了！