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?
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.