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

根据 Dirk 的回复here：

有没有办法在不使用显式循环的情况下修复代码 在 C++ 代码中？

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

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

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

但问题是（再次引用 Dirk），

而且我建议在花费大量精力之前检查一下 你是否有可能比 R 函数 rbinom 做得更好。那 R 函数在 C 代码中被矢量化，你不太可能做事 使用 Rcpp 更快，除非你想使用随机变量 在另一个 C++ 函数中。

而且它实际上更慢（在我的机器上是 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，也许可以节省一些时间，避免一些开销。