`neg_binomial_2_lcdf()` "grad_2F1: k (internal counter) exceeded 100000 iterations"


I am running a regression with a negative binomial error. The outcome is occasionally censored from below, but I am getting the following error on the line with the neg_binomial_2_lcdf call:

Exception: grad_2F1: k (internal counter) exceeded 100000 iterations, hypergeometric function gradient did not converge. (in ‘model4923325fd331_nb_std_censored’ at line 30)

The fact that not all chains fail makes me think it might be an issue with the prior or the initialization. Any ideas would be greatly appreciated. Model and fake data below.

data {
  int N; //the number of observations
  int<lower=0,upper=1> detected[N]; // whether the outcome is detected (i.e. not censored)
  int K; //the number of columns in the design matrix
  int<lower=0> y[N]; //the response
  matrix[N,K] X; //the design matrix
  real reciprocal_phi_scale;
parameters {
  vector[K] beta; //the regression parameters
  real<lower=0> reciprocal_phi;
transformed parameters {
  real phi;
  phi = 1. / reciprocal_phi;

model {
  vector[N] Xbeta;
  Xbeta = X * beta;
  beta ~ cauchy(0,10);
  reciprocal_phi ~ cauchy(0,reciprocal_phi_scale);
  for (i in 1:N) {
    if (detected[i]) {
      y[i] ~ neg_binomial_2_log(Xbeta[i], phi);
    } else {
      target += neg_binomial_2_lcdf(y[i] | exp(Xbeta[i]), phi);
generated quantities {
  vector[N] y_sim;
  for (i in 1:N) {
    y_sim[i] = neg_binomial_2_log_rng(X[i]*beta, phi);

Without covariates rarely gives an error. Note that I’m taking observations \leq 3 and observing them as censored \leq 3.

N <- 100
y <- rnbinom(N, mu=6, size=2)
nb_data <- list(N=N,
             y=ifelse(y > 3, y, 3),
             detected= ifelse(y > 3, 1, 0),

And some data that occasionally gives the error. My real data gives the error reliably but I can’t share it.

N <- 100
x <- rnorm(N, sd=1)
beta <- c(1, .5)
y <- rnbinom(N, mu=exp(beta[1]+beta[2]*x), size=2)

nb_data <- list(N=N,
             y=ifelse(y > 3, y, 3),
             detected= ifelse(y > 3, 1, 0),
             X=data.frame(c=rep(1,N), x=x),

Thanks for the thorough example. For this particular example could you add a print statement that outputs the neg_binomial_2_lcdf parameters immediately prior to the usage of that function? I’m guessing the gradient is just underflowing. While that function is now pretty robust it may be that we could change how it is used and avoid this error. Do you get it consistently during sampling or just occasionally?


Thanks for the quick reply. Here is an example that gives the error:


I get the error about half the times I fit my data (less often with the simulated data above). As soon as I get the error, that chain crashes.

Interestingly, I don’t get it with BRMS, which can also fit this model. I guess the priors used there avoid it.


That shouldn’t be crashing a chain, this is good to know. I’ll check what is…


Actually brms is getting the same error resulting in chains crashing, too, so I’m stuck. Is there a GitHub issue where I can track this? (Or should I make one, and if so in which repo)?


Yeah, this is in the Stan math library, it’ll be all the interfaces. There’s no issue yet because in the past I could not reproduce it using just the math library. I’ll post here when I make some progress…