r - 矢量化 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?

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

这是我对“矢量化”代码的实现:

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