矢量化 Rcpp 随机二项式绘制 [英] Vectorised Rcpp random binomial draws

问题描述

这是这个问题的后续问题:Generating sameRcpp和R中的随机变量

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

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

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

 cppFunction("NumericVector cpprbinom(int n, double size, double prob) { \返回(rbinom(n，大小，概率))；}")设置种子(42)；cpprbinom(10, 1, 0.5)

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

 cpprbinom(length(x), 1, x)

如何修改 Rcpp 代码来实现这一点?

谢谢

解决方案

Following Dirk's response 这里:

<块引用><块引用>

有没有一种不使用显式循环来修复代码的方法在 C++ 代码中?

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

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

库(Rcpp)cppFunction("NumericVector cpprbinom(int n, double size, NumericVector prob) {数值向量 v(n);for (int i=0; i<n; i++) {v[i] = as<double>(rbinom(1, size, prob[i]));}退货(五)；}")r <- runif(1e6)all.equal({set.seed(42); rbinom(length(r), 1, r)},{set.seed(42);cpprbinom(长度(r)，1，r)})#真的

但问题是(再次引用德克)，

<块引用>

我建议在为此花费大量精力之前先检查一下你是否有可能比 R 函数 rbinom 做得更好.那R 函数在 C 代码中被向量化，你不太可能制作东西使用 Rcpp 快得多，除非你想使用随机变量在另一个 C++ 函数中.

而且它实际上更慢(在我的机器上为 x3)，所以至少像我这样幼稚的实现无济于事:

库(微基准)微基准(rbinom(长度(r)，1，r)，cpprbinom(长度(r)，1，r))单位:毫秒expr min lq 平均中位数 uq max nevalrbinom(length(r), 1, r) 55.50856 56.09292 56.49456 56.45297 56.65897 59.42524 100cpprbinom(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); });返回(v);}")r <- runif(1e6)all.equal({set.seed(42); rbinom(length(r), 1, r)},{set.seed(42);cpprbinom(长度(r), 1, r)},{set.seed(42);cpprbinom2(长度(r), 1, r)})#真的微基准(rbinom(长度(r)，1，r)，cpprbinom(长度(r)，1，r)，cpprbinom2(长度(r)，1，r))单位:毫秒expr min lq 平均中位数 uq max nevalrbinom(length(r), 1, r) 55.26412 56.00314 56.57814 56.28616 56.59561 60.01861 100cpprbinom(length(r), 1, r) 113.72513 115.94758 122.81545 117.24708 119.95134 168.47246 100cpprbinom2(length(r), 1, r) 36.67589 37.12182 38.95318 37.37436 37.97719 84.73516 100

This is a follow-on question from this one: Generating same random variable in Rcpp and R

I'm trying to speed up a vectorised call to rbinom of this form:

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

In the live code of x is a vector of variable length (but typically numbering in the millions). I have no experience with Rcpp but I was wondering could I use Rcpp to speed this up. From the linked question this Rcpp code was suggested for non-vectorised rbinom calls by @Dirk Eddelbuettel :

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

....and is about twice as fast as the non Rcpp option, but can't handle my vectorised version

    cpprbinom(length(x), 1, x)

How can the Rcpp code be modified to implement this?

Thanks

解决方案

Following Dirk's response here:

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

I don't think so. The code currently has this hard-wired: <...> so until one of us has sufficient [time] to extend this (and test it) will have to do the loop at your end.

Here's my implementation of a "vectorised" 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

But the problem is (again citing Dirk),

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.

And it is actually slower (x3 on my machine), so at least such naive implementation as mine won't help:

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

EDIT: according to Romain's comment below, here's an advanced version, which is faster!

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

这篇关于矢量化 Rcpp 随机二项式绘制的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！