将 R 翻译成 Rcpp：迭代器

Question

我正在将一个函数从 R 翻译成 Rcpp，这让我苦恼了一段时间。我已经阅读了一些这样的资料one，但我是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));

我的猜测是我没有正确使用迭代器。我关注了这个 example.

最后，我在 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

我是一个简单的人，不知道 IntegerVector::iterator 是什么。我只是这样做：

#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;
  double n = group_index.length();
  int g;
  
  for (double j = 0; j < n_mcmc; ++j) {
    NumericMatrix parm(2, n_group);
    NumericMatrix sample(n_mcmc, n_group);
    
    for(double i = 0; i < n; ++i) {
      g = group_index(i);
      NumericVector temp = Rcpp::rbinom(m(g) , 1, theta_sample(j, g));
      
      parm(0,g) = prior_parm(0,g) + a(g) + sum(temp);
      parm(1,g) = prior_parm(1,g) + b(g) + m(g) - sum(temp);
      sample(_,g) = Rcpp::rbeta(n_mcmc, parm(0,g), parm(1,g));
    }
    
    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;
}

请注意，R 支持大于最大整数的向量。因此，我更喜欢在 for 循环中使用双精度而不是整数。我不知道这是否也应该成为您的组指数的考虑因素。

Answer 2

您的问题是由于没有取消对带星号的迭代器的引用。您会在此处的示例中注意到 Chapter 28 Iterator | 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));
}

这很快就会变得混乱。对我来说，最好使用 range base loops:

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)]] 放入源文件中（参见 First steps in using C++11 with Rcpp）。

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

将 R 翻译成 Rcpp：迭代器

Translating R to Rcpp: iterators

r

rcpp