将 R 转换为 Rcpp：迭代器

Question

我正在将一个函数从 R 翻译成 Rcpp，我已经挣扎了一段时间。我看过这样的一些材料一个，但我在RCPP一个初学者，我一直没能找出我目前做错了。

我在 Rcpp 中的最佳结果如下，

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
Rcpp::NumericVector aux_pred_C(IntegerVector m, 
                               NumericVector a, 
                               NumericVector b, 
                               int n_group, 
                               IntegerVector group_index,
                               NumericMatrix theta_sample, 
                               IntegerMatrix prior_parm, 
                               int n_mcmc){
  
  NumericVector p(n_mcmc);
  group_index = group_index - 1;

  for (int j = 0; j < n_mcmc; ++j) {
    NumericMatrix parm(2, n_group);
    NumericMatrix sample(n_mcmc, n_group);
    
    for(IntegerVector::iterator i = group_index.begin(); i != group_index.end(); ++i) {
      NumericVector temp = Rcpp::rbinom(m(i) , 1, theta_sample(j, i));
      
      parm(0,i) = prior_parm(0,i) + a(i) + sum(temp);
      parm(1,i) = prior_parm(1,i) + b(i) + m(i) - sum(temp);
      sample(_,i) = Rcpp::rbeta(n_mcmc, parm(0,i), parm(1,i));
    }

    int i1 = group_index[0];
    int i2 = group_index[1];
    LogicalVector V = sample(_,i2) > sample(_,i1);
    IntegerVector res = ifelse(V, 1, 0);
    p(j) = mean(res);
  }
  
  return(p);
}

我的错误消息在第 22、24、25 和 26 行重复出现：

sampler_predictive_distribution.cpp:22:44: error: no match for call to '(Rcpp::IntegerVector {aka Rcpp::Vector<13>}) (Rcpp::traits::storage_type<13>::type*&)'

invalid conversion from 'Rcpp::Vector<13>::iterator' {aka 'int*'} to 'size_t' {aka 'long long unsigned int'} [-fpermissive] NumericVector temp = Rcpp::rbinom(m(i) , 1, theta_sample(j, i));

我的猜测是我没有正确使用迭代器。我跟着这个例子。

最后，我在 R 中的可重现示例，

n.group <- 4
group.index <- c(1, 3)
m <- rep(NA, n.group)
m[group.index[1]] <- 50
m[group.index[2]] <- 50
a <- rep(NA, n.group)
a[group.index[1]] <- 20
a[group.index[2]] <- 20
b <- rep(NA, n.group)
b[group.index[1]] <- 30
b[group.index[2]] <- 30

n.mcmc <- 100
theta.sample.e1 <- matrix(NA, nrow = n.mcmc, ncol = n.group)
theta.sample.e1[, group.index[1]] <- rbeta(n.mcmc, a[group.index[2]], b[group.index[1]])
theta.sample.e1[, group.index[2]] <- rbeta(n.mcmc, a[group.index[2]], b[group.index[2]])
prior.parm.e1 <- matrix(1, ncol = 4, nrow = 2)

aux_pred <- function(m, a, b, n.group, group.index, theta.sample, prior.parm, n.mcmc){
  
  p <- rep(NA, n.mcmc)
  
  for (j in 1:n.mcmc){
    parm <- matrix(NA, ncol = n.group, nrow = 2)
    sample <- matrix(NA, ncol = n.group, nrow = n.mcmc)
    
    for (i in group.index){
      temp <- rbinom(m[i], size = 1, prob =  theta.sample[j, i])
      
      parm[1, i] <- prior.parm[1, i] + a[i] + sum(temp)
      parm[2, i] <- prior.parm[2, i] + b[i] + length(temp) - sum(temp)
      sample[, i] <-  rbeta(n.mcmc, parm[1, i], parm[2, i])
    }
    p[j] <- mean(sample[, group.index[2]] - sample[, group.index[1]] > 0)
  }
  
  return(p)
}

aux_pred(m, a, b, n.group, group.index, theta.sample.e1, prior.parm.e1, n.mcmc)

我错过了什么？

Answer 1

您的问题源于未使用星号取消引用迭代器。您将在此处的示例中注意到第 28 章迭代器 | Rcpp for Everyone，当用迭代器遍历向量的元素时，为了得到特定迭代器指向的值，我们必须用星号取消引用迭代器。所以，你的代码看起来像：

for(IntegerVector::iterator i = group_index.begin(); i != group_index.end(); ++i) {
    NumericVector temp = Rcpp::rbinom(m(*i) , 1, theta_sample(j, *i));
    
    parm(0,*i) = prior_parm(0,*i) + a(*i) + sum(temp);
    parm(1,*i) = prior_parm(1,*i) + b(*i) + m(*i) - sum(temp);
    sample(_,*i) = Rcpp::rbeta(n_mcmc, parm(0,*i), parm(1,*i));
}

这会很快变得混乱。对我来说，最好使用范围基础循环：

for(auto i: group_index) {
    NumericVector temp = Rcpp::rbinom(m(i) , 1, theta_sample(j, i));
    
    parm(0,i) = prior_parm(0,i) + a(i) + sum(temp);
    parm(1,i) = prior_parm(1,i) + b(i) + m(i) - sum(temp);
    sample(_,i) = Rcpp::rbeta(n_mcmc, parm(0,i), parm(1,i));
}

有一个基于范围的 for 循环的警告。也就是说，它只是因为可用的C++11，但是如果你使用的最新版本的不应该是一个令人关切的是大R这使得C++11默认。或者，您始终可以放入// [[Rcpp::plugins(cpp11)]]源文件（请参阅在 Rcpp 中使用 C++11 的第一步）。

此外，您正在使用(检查向量边界的访问器，因此效率不如[. 后者更危险，并且让开发人员有责任确保您不会访问不属于您的内存，但在实践中，保证这一承诺并不难。见第 8 章第 2 节 | 适合所有人的 Rcpp。