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?