Hi,
I'm running a MCMC estimation of a Bayesian model, and I'm running into a
particular error during the MCMC in my likelihood function call. I've
implemented my likelihood function in C++ using RCpp, and I get this
"unimplemented type" error:
Error in loglik.func(): unimplemented type 'double' in 'coerceToReal.'
The nasty thing about this is that it is very hard to replicate, and only
occurs after a particular iteration in my MCMC. While my function stalls when
this errors occurs, all the other times it returns the correct values. My
specific question is: are there any particular functions in the Rcpp library
that calls coerceToReal?" I'm trying to understand what would cause Rcpp to
throw an error like this (eg. memory leaks, mismatching of types...etc). A
Google search returns very few relevant results on this type of error.
I think I've isolated the problem area in my code to this function:
//////////////////////////////////////////
Rcpp::NumericVector pbivpois(arma::colvec x, arma::colvec y, arma::rowvec
lambda, bool returnLog) {
// Unlike the pbivpois function, this assumes that x is a vector, and does
not error check.
// Declare Internal Variables
int N = x.n_rows;
Rcpp::NumericVector logbp(N);
Rcpp::NumericVector result(N);
int x0 = 0;
int y0 = 0;
int xymin = 0;
Rcpp::NumericVector i(2);
arma::rowvec sums = arma::zeros<arma::rowvec>(3);
double lambdaratio = 0;
double maxsums = 0;
double logsummation =0;
for(int k = 0; k < N; k++) {
x0 = x.row(k);
y0 = y.row(k);
//////////////////////////////////////////////////
// I think these lines are causing the problem.
xymin = std::min(x0, y0);
lambdaratio = lambda.col(2)/(lambda.col(0)*lambda.col(1));
i = Rcpp::seq(0,xymin);
sums = -Rcpp::lgamma(x0 - i + 1)-Rcpp::lgamma(i+1) -
Rcpp::lgamma(y0-i+1) + i*std::log(lambdaratio);
maxsums = sums.max();
sums = sums-maxsums;
logsummation = std::log(arma::accu(arma::exp(sums))) + maxsums;
logbp[k] = -arma::accu(lambda) + x0*std::log(lambda.col(0)) +
y0*std::log(lambda.col(1)) + logsummation;
////////////////////////////////////
}
if( returnLog == TRUE ) {
result = logbp;
} else {
result = Rcpp::exp(logbp);
}
return(result);
}
////////////////////////////////////////////////
Thanks for your help!
Best,
Clarence
_______________________________________________
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
