Sampling in cmdstan often (but not always) fails with error, while rstan always works

Hi everyone,

I’m trying to fit a model with a custom likelihood function (Skellam distribution).
The model works (compiles, sampling seems to be fine, estimates are ok) in rstan.
In cmdstanr the model compiles, too. But sampling fails occasionally (actually most of the time) with two different error messages, but it also occasionally samples fine (including estimates that look ok).

The error messages I get point to issues in boost:math, but that’s about all I can extract from these messages (reprex below).

Is there anything that can be done to make this model work using cmdstanr?

The reason I’m asking is because I would really like to use cmdstanr, but if need be I can switch to rstan for this project.

library(cmdstanr)
library(rstan)

# data:
set.seed(1)
sdat <- list(N = 50)
theta1 <- 0.5
theta2 <- 2
sdat$diff_y <- rpois(sdat$N, 0.5) - rpois(sdat$N, 2)

scode <- "
functions {
  // from https://github.com/LeoEgidi/footBayes/blob/master/R/stan_foot.R
  real skellam_lpmf(int k, real lambda1, real lambda2) {
    real r = k;
    return -(lambda1 + lambda2) + (r/2) * log(lambda1/lambda2) +
      log(modified_bessel_first_kind(k, 2 * sqrt(lambda1 * lambda2)));
  }
}
data {
  int N;
  array[N] int diff_y; // frequency differential
}
parameters {
  real theta1;
  real theta2;
}
model {
  for (n in 1:N){
    target+= skellam_lpmf(diff_y[n] | exp(theta1), exp(theta2));
  }
  theta1 ~ normal(0,1);
  theta2 ~ normal(0,1);
}
generated quantities {
  real theta1_exp = exp(theta1);
  real theta2_exp = exp(theta2);
}
"

# compile models
mod_cmdstan <- cmdstan_model(write_stan_file(scode))
mod_rstan <- stan_model(write_stan_file(scode))

# sampling fails
mod_cmdstan$sample(data = sdat, chains = 1, iter_warmup = 200, iter_sampling = 200, seed = 4)
# Error in function boost::math::bessel_ik<long double>(long double,long double) in CF2_ik: Series evaluation exceeded 1000000 iterations, giving up now.

# sampling fails
mod_cmdstan$sample(data = sdat, chains = 1, iter_warmup = 200, iter_sampling = 200, seed = 7)
# Error in function boost::math::asymptotic_bessel_i_large_x<long double>(long double,long double): Overflow Error

# sampling is fine (estimates look ok)
mod_cmdstan$sample(data = sdat, chains = 1, iter_warmup = 200, iter_sampling = 200, seed = 3)

# rstan samples without complaints
sampling(mod_rstan, data = sdat, chains = 1, iter = 1000)

Operating System: MacOS
Interface Version: cmdstanr (0.8.1.9000), CmdStan (2.35.0), rstan (2.32.6) with Stan version 2.32.2
Compiler/Toolkit:

These 2 errors are common when you using complex functions in boost. If for most time the sampler runs well and you can obtail reasonable posterior then don’t worry about these.

1 Like

Thank you for the response. That’s good to know.

The problem remains that sampling in cmdstanr fails unpredictably, i.e. I got it to work by trying different seeds and by using a low number of iterations. If I increase the number of iterations the probability of these errors goes up, too. And if I use a reasonable number (e.g. 4*1000) it essentially fails in 100% of the attempts.