向量化 Rcpp 随机二项式绘制

Question

这是这个问题的后续问题：Generating same random variable in Rcpp and R

我正在尝试加快对这种形式的 rbinom 的矢量化调用：

    x <- c(0.1,0.4,0.6,0.7,0.8)
    rbinom(length(x),1 ,x)

在 x 的实时代码中是一个可变长度的向量（但通常以数百万计）。我没有使用 Rcpp 的经验，但我想知道我可以使用 Rcpp 来加快速度。从链接的问题来看，@Dirk Eddelbuettel 建议将此 Rcpp 代码用于非矢量化 rbinom 调用：

    cppFunction("NumericVector cpprbinom(int n, double size, double prob) { \
         return(rbinom(n, size, prob)); }")
    set.seed(42); cpprbinom(10, 1, 0.5)

...并且大约是非 Rcpp 选项的两倍，但无法处理我的矢量化版本

    cpprbinom(length(x), 1, x)

如何修改 Rcpp 代码来实现这个？

谢谢

Answer 1

继德克的回应后 here：

Is there a way of fixing the code without using an explicit loop in the C++ code?

我不这么认为。该代码目前具有此硬连线：<...> 所以直到我们中的一个人有足够的 [时间] 来扩展它（并测试它）在你的最后做循环。

这是我对 "vectorised" 代码的实现：

library(Rcpp)
cppFunction("NumericVector cpprbinom(int n, double size, NumericVector prob) { 
    NumericVector v(n);            
    for (int i=0; i<n; i++) {v[i] = as<double>(rbinom(1, size, prob[i]));} 
    return(v); }")
r <- runif(1e6)
all.equal({set.seed(42); rbinom(length(r), 1, r)}, 
          {set.seed(42); cpprbinom(length(r), 1, r)})
#TRUE

但问题是（再次引用德克），

And I suggest that before expending a lot of effort on this you check whether you are likely to do better than the R function rbinom. That R function is vectorized in C code and you are unlikely to make things much faster by using Rcpp, unless you want to use the random variates in another C++ function.

而且它实际上更慢（在我的机器上是 x3），所以至少像我这样天真的实现不会有帮助：

library(microbenchmark)
microbenchmark(rbinom(length(r), 1, r), cpprbinom(length(r), 1, r))

Unit: milliseconds
                       expr       min        lq      mean    median        uq       max neval
    rbinom(length(r), 1, r)  55.50856  56.09292  56.49456  56.45297  56.65897  59.42524   100
 cpprbinom(length(r), 1, r) 117.63761 153.37599 154.94164 154.29623 155.37247 225.56535   100

编辑：根据下面 Romain 的评论，这是一个高级版本，速度更快！

cppFunction(plugins=c("cpp11"), "NumericVector cpprbinom2(int n, double size, NumericVector prob) { 
    NumericVector v = no_init(n);
    std::transform( prob.begin(), prob.end(), v.begin(), [=](double p){ return R::rbinom(size, p); }); 
    return(v);}")
r <- runif(1e6)
all.equal({set.seed(42); rbinom(length(r), 1, r)}, 
          {set.seed(42); cpprbinom(length(r), 1, r)}, 
          {set.seed(42); cpprbinom2(length(r), 1, r)})
#TRUE
microbenchmark(rbinom(length(r), 1, r), cpprbinom(length(r), 1, r), cpprbinom2(length(r), 1, r))

Unit: milliseconds
                        expr       min        lq      mean    median        uq       max neval
     rbinom(length(r), 1, r)  55.26412  56.00314  56.57814  56.28616  56.59561  60.01861   100
  cpprbinom(length(r), 1, r) 113.72513 115.94758 122.81545 117.24708 119.95134 168.47246   100
 cpprbinom2(length(r), 1, r)  36.67589  37.12182  38.95318  37.37436  37.97719  84.73516   100

Answer 2

不是通用解决方案，但我注意到您在调用 rbinom 时将 size 参数设置为 1。如果情况总是如此，您可以绘制 length(x) 个统一值，然后与 x 进行比较。例如：

 set.seed(123)
 #create the values
 x<-runif(1000000)
 system.time(res<-rbinom(length(x),1 ,x))   
 # user  system elapsed 
 #0.068   0.000   0.070
 system.time(res2<-as.integer(runif(length(x))<x))   
 # user  system elapsed 
 #0.044   0.000   0.046

收获不大，但如果从 C++ 调用 runif 可能会节省一些时间，从而避免一些开销。

向量化 Rcpp 随机二项式绘制

Vectorised Rcpp random binomial draws

r

rcpp