The poisson_rng, which is used by neg_binomial_2_log_rng is not very robust. I have worked around this problem with
real mean_PPD;
mean_PPD <- 0;
for (n in 1:N) {
real gamma_temp;
if (is_inf(theta)) gamma_temp <- nu[n];
else gamma_temp <- gamma_rng(theta, theta / nu[n]);
if (gamma_temp < poisson_max)
mean_PPD <- mean_PPD + poisson_rng(gamma_temp);
else mean_PPD <- mean_PPD + normal_rng(gamma_temp, sqrt(gamma_temp));
}
mean_PPD <- mean_PPD / N;
where poisson_max is pow(2.0, 30.0), theta is the overdispersion, and nu is a vector of (positive) parameters. It would be slightly different if you are utilizing the _log variants.
Ben