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

Interfaces
,
9

It’s probably a Stan math library issue math arose between Stan 2.32 and Stan 2.35 (CmdStanR is more up to date as you can see from versions). @sbronder or @WardBrian may know more as they’ve been looking into how our complex numbers are handled due to changing ground in the standard libraries.

If you can initialize from a sensible place, that can be huge help given what you’ve reported. Sometimes that only requires reducing the scale of random initialization from (-2, 2) to (-0.1, 0.1) or something like that (it samples uniformly on the unconstrained scale after constrained parameters are transformed).

If things are failing after warmup, then there are bigger problems.

exp() is dangerous due to the potential for overflow, but then with the density you’ve defined, I’m not sure how to deal with that. So I’d recommend just coding things this way for simplicity:

parameters {
  vector<lower=0> theta;  // formerly called theta_exp
}
model {
  theta ~ lognormal(0, 1);
}

You then want to vectorize your Skellam lpdf so that you can write:

diff_y ~ skellam(theta[1], theta[2]);

That way, you can compute the normalization term -(lambda1 + lambda2) + (r/2) * log(lambda1/lambda2) just once.

real skellam_lpmf(array[] int k, real lambda1, real lambda2) {
  int N = size(k);
  vector[size] r = to_vector(k);
  real norm = -(lambda1 + lambda2) + (r/2) * log(lambda1/lambda2);
  real lp = N * norm;
  for (n in 1:N) {
      lp += log(modified_bessel_first_kind(k[n], 2 * sqrt(lambda1 * lambda2)));
  }
  return lp;
}

You can then look at where the log modified Bessel is failing and whether it can be stabilized somehow or if there’s a reasonable value to set if it’s about to overflow/underflow (though that usually loses derivatives, so can be a real problem if it persists).

Maybe @bgoodri knows what’s up with the log Bessel function implementations: