2

I am translating a function from R to Rcpp, and I have been struggling for some time. I have read a few materials like this one, but I am a beginner in Rcpp and I have not been able to figure out what I am currently doing wrong.

My best result in Rcpp is below,

#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);
}

My errors messages are repeated for lines 22, 24, 25, and 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));

My guess is that I am not using iterators correctly. I followed this example.

Finally, my reproducible example in 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)

What am I missing?

0

2 Answers 2

5

Your problem is stemming from not dereferencing the iterator with an asterisk. You will note in the example here Chapter 28 Iterator | Rcpp for everyone, when one loops over the elements of a vector with an iterator, in order to get the value pointed to by the specific iterator, we must dereference the iterator with an asterisk. So, your code would look something like:

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));
}

This can get messy quickly. To me, it is preferable to use 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));
}

There is a caveat with range based for loops. That is, it is only available since C++11, however this shouldn't be that big of a concern if you are using a recent version of R which made C++11 the default. Alternatively, you can always put // [[Rcpp::plugins(cpp11)]] in your source file (See First steps in using C++11 with Rcpp).

Additionally, you are using ( accessors which check vector bounds and are thus not as efficient as [. The latter is more dangerous and puts the onus on the developer to ensure that you are not accessing memory that you do not own, however in practice it isn't that difficult to guarantee this promise. See Chapter 8 Section 2 | Rcpp for everyone.

3

I'm a simple man who has no idea what an IntegerVector::iterator is. I'd simply do this:

#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;
}

Note that R supports vectors larger than max integer. I prefer using doubles instead of ints in the for loops because of that. I don't know if that should also be a consideration for your group indices.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.