矢量化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來加快速度。從鏈接的問題，這RCPP代碼有人建議對非矢量化rbinom調用由@Dirk Eddelbuettel：

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：

有固定不使用顯式循環 在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 

但問題是，（再次引用德克），

而且我認爲上花費了很大的力氣才這你檢查 你是否可能比R函數更好。 R函數是用C代碼矢量化的，並且除非要在其他C++函數中使用隨機變量 ，否則您不太可能通過使用Rcpp使事情變得更加快速 。

而且它實際上是比較慢（我的機器上三次），所以至少這種天真的實現礦山不會幫助：

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 

編輯：根據以下羅曼的評論，這裏是一個高級版本，這是更快！

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，避免一些開銷。