# Initialization failure in negative-binomial model with exact and right-censored likelihood contributions

**URL:** <https://discourse.mc-stan.org/t/initialization-failure-in-negative-binomial-model-with-exact-and-right-censored-likelihood-contributions/41644>\
**Category:** Modeling\
**Tags:** fitting-issues, cmdstanr\
**Created:** [October 5, 2026, 2:57pm UTC](https://discourse.mc-stan.org/t/initialization-failure-in-negative-binomial-model-with-exact-and-right-censored-likelihood-contributions/41644 "2026-10-05T14:57:45Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![arumadri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/arumadri/32/20838_2.png) [@arumadri](https://discourse.mc-stan.org/u/arumadri)\
**Post date:** [October 5, 2026, 2:57pm UTC](https://discourse.mc-stan.org/t/initialization-failure-in-negative-binomial-model-with-exact-and-right-censored-likelihood-contributions/41644/1 "2026-10-05T14:57:45Z")

</div>

I’m trying to fit a large spatial-temporal Bayesian model with close to 1,700 areas across 10 years in Stan using a negative-binomial model with exact and right-censored likelihood contributions. Predictors include meterological variables built using DLNM statistical framework.

Here is the Stan code of the model ```

```stan

data {
  // ---------------------------------------------------------------------------
  // Exact observed data
  // ---------------------------------------------------------------------------

  int<lower=0> N;
  // Number of exactly observed observations

  array[N] int<lower=0> Y;
  // Exactly observed case counts

  int<lower=2> K;
  // Intercept and at least one basis column

  matrix[N, K] X;
  // Design matrix including the intercept column

  vector[N] log_offset;
  // Log population offset

  int<lower=1> N_month_lfo;
  // Number of calendar months

  int<lower=1> S;
  // Number of subcounties

  array[N] int<lower=1, upper=S> subcounty_id;
  // Subcounty index for exact observations

  // ---------------------------------------------------------------------------
  // Right-censored observed data
  // ---------------------------------------------------------------------------

  int<lower=0> N_cens;
  // Number of right-censored observations

  array[N_cens] int<lower=1> Y_cens;
  // Censored lower-bound counts

  matrix[N_cens, K] X_cens;
  // Design matrix for censored observations

  vector[N_cens] log_offset_cens;
  // Log population offset for censored observations

  array[N_cens] int<lower=1, upper=S> subcounty_id_cens;
  // Subcounty index for censored observations

  // ---------------------------------------------------------------------------
  // Held-out exact data for M-step-ahead LFO-CV
  // ---------------------------------------------------------------------------

  int<lower=0> N_test;
  // Number of held-out exact observations

  array[N_test] int<lower=0> Y_test;
  // Held-out exact case counts

  matrix[N_test, K] X_test;
  // Held-out exact design matrix

  vector[N_test] log_offset_test;
  // Held-out exact log population offset

  array[N_test] int<lower=1, upper=N_month_lfo> month_id_test;
  // Held-out exact LFO-CV month index

  array[N_test] int<lower=1, upper=S> subcounty_id_test;
  // Held-out exact subcounty index

  // ---------------------------------------------------------------------------
  // Held-out right-censored data for M-step-ahead LFO-CV
  // ---------------------------------------------------------------------------

  int<lower=0> N_test_cens;
  // Number of held-out censored observations

  array[N_test_cens] int<lower=1> Y_test_cens;
  // Held-out censored lower-bound counts

  matrix[N_test_cens, K] X_test_cens;
  // Held-out censored design matrix

  vector[N_test_cens] log_offset_test_cens;
  // Held-out censored log population offset

  array[N_test_cens] int<lower=1, upper=N_month_lfo>
    month_id_test_cens;
  // Held-out censored LFO-CV month index

  array[N_test_cens] int<lower=1, upper=S>
    subcounty_id_test_cens;
  // Held-out censored subcounty index

  // ---------------------------------------------------------------------------
  // LFO-CV indices
  // ---------------------------------------------------------------------------

  array[N] int<lower=1, upper=N_month_lfo> month_id;
  // LFO-CV month index for exact observations

  array[N_cens] int<lower=1, upper=N_month_lfo> month_id_cens;
  // LFO-CV month index for censored observations

  int<lower=0, upper=1> prior_only;
  // 1 = prior predictive sampling; 0 = fit to observations

  // ---------------------------------------------------------------------------
  // Priors
  // ---------------------------------------------------------------------------

  vector[K - 1] prior_b_mean;
  // Prior means for fixed-effect coefficients

  vector[K - 1] prior_b_sd;
  // Prior standard deviations for fixed-effect coefficients

  real prior_intercept_mean;
  // Prior mean for the intercept

  real prior_intercept_sd;
  // Prior standard deviation for the intercept

  real prior_shape;
  // Rate parameter for the exponential prior on NB shape

  real<lower=0> prior_sigma_subcounty;
  // Prior standard deviation for between-subcounty variation
}

parameters {
  // ---------------------------------------------------------------------------
  // Fixed effects
  // ---------------------------------------------------------------------------

  vector[K - 1] b;
  // Fixed-effect coefficients

  real Intercept;
  // Overall model intercept

  real<lower=0> shape;
  // Negative-binomial overdispersion/precision parameter

  // ---------------------------------------------------------------------------
  // Subcounty effects
  // ---------------------------------------------------------------------------

  vector[S] z_subcounty;
  // Standardized subcounty effects

  real<lower=0> sigma_subcounty;
  // Between-subcounty standard deviation
}

transformed parameters {
  // ---------------------------------------------------------------------------
  // Non-centred subcounty random intercepts
  // ---------------------------------------------------------------------------

  vector[S] subcounty_effect =
    sigma_subcounty * z_subcounty;
}

model {
  // ---------------------------------------------------------------------------
  // Priors
  // ---------------------------------------------------------------------------

  b ~ normal(prior_b_mean, prior_b_sd);
  Intercept ~ normal(prior_intercept_mean, prior_intercept_sd);
  shape ~ exponential(prior_shape);
  z_subcounty ~ std_normal();
  sigma_subcounty ~ normal(0, prior_sigma_subcounty);

  // ---------------------------------------------------------------------------
  // Likelihood
  // ---------------------------------------------------------------------------

  if (!prior_only) {
    vector[N] eta =
      Intercept +
      subcounty_effect[subcounty_id] +
      log_offset +
      X[, 2:K] * b;

    if (N > 0) {
      Y ~ neg_binomial_2_log(eta, shape);
    }

    if (N_cens > 0) {
      vector[N_cens] eta_cens =
        Intercept +
        subcounty_effect[subcounty_id_cens] +
        log_offset_cens +
        X_cens[, 2:K] * b;

      for (n in 1:N_cens) {
        target += neg_binomial_2_lccdf(
          Y_cens[n] - 1 |
          exp(eta_cens[n]),
          shape
        );
      }
    }
  }
}

generated quantities {
  // ---------------------------------------------------------------------------
  // Training-data pointwise log likelihoods
  // ---------------------------------------------------------------------------

  vector[N] log_lik_obs;
  // Exact-observation log likelihoods

  vector[N_cens] log_lik_cens;
  // Censored-observation log likelihoods

  vector[N_month_lfo] log_lik_month;
  // Training log likelihood summed by calendar month

  // ---------------------------------------------------------------------------
  // Held-out pointwise log likelihoods
  // ---------------------------------------------------------------------------

  vector[N_test] log_lik_test;
  // Held-out exact log likelihoods

  vector[N_test_cens] log_lik_test_cens;
  // Held-out censored log likelihoods

  vector[N_month_lfo] log_lik_test_month;
  // Held-out log likelihood summed by calendar month

  // ---------------------------------------------------------------------------
  // Posterior predictive counts
  // ---------------------------------------------------------------------------

  array[N] int Y_rep;
  // Posterior predictive exact training counts

  array[N_cens] int Y_rep_cens;
  // Posterior predictive latent counts for censored training rows

  array[N_test] int Y_rep_test;
  // Posterior predictive exact held-out counts

  array[N_test_cens] int Y_rep_test_cens;
  // Posterior predictive latent counts for censored held-out rows

  // ---------------------------------------------------------------------------
  // Posterior expected counts
  // ---------------------------------------------------------------------------

  vector[N] mu_hat;
  // Posterior expected counts for exact training observations

  vector[N_cens] mu_hat_cens;
  // Posterior expected latent counts for censored training observations

  vector[N_test] mu_hat_test;
  // Posterior expected counts for exact held-out observations

  vector[N_test_cens] mu_hat_test_cens;
  // Posterior expected latent counts for censored held-out observations

  // Initialize monthly log-likelihood totals
  log_lik_month = rep_vector(0, N_month_lfo);
  log_lik_test_month = rep_vector(0, N_month_lfo);

  {
    // -------------------------------------------------------------------------
    // Exact training observations
    // -------------------------------------------------------------------------

    vector[N] eta =
      Intercept +
      subcounty_effect[subcounty_id] +
      log_offset +
      X[, 2:K] * b;

    if (N > 0) {
      for (n in 1:N) {
        mu_hat[n] = exp(eta[n]);

        Y_rep[n] =
          neg_binomial_2_log_rng(eta[n], shape);

        log_lik_obs[n] =
          neg_binomial_2_log_lpmf(
            Y[n] | eta[n], shape
          );

        log_lik_month[month_id[n]] +=
          log_lik_obs[n];
      }
    }

    // -------------------------------------------------------------------------
    // Right-censored training observations
    // -------------------------------------------------------------------------

    if (N_cens > 0) {
      vector[N_cens] eta_cens =
        Intercept +
        subcounty_effect[subcounty_id_cens] +
        log_offset_cens +
        X_cens[, 2:K] * b;

      for (n in 1:N_cens) {
        mu_hat_cens[n] = exp(eta_cens[n]);

        Y_rep_cens[n] =
          neg_binomial_2_log_rng(
            eta_cens[n], shape
          );

        log_lik_cens[n] =
          neg_binomial_2_lccdf(
            Y_cens[n] - 1 |
            exp(eta_cens[n]),
            shape
          );

        log_lik_month[month_id_cens[n]] +=
          log_lik_cens[n];
      }
    }

    // -------------------------------------------------------------------------
    // Exact held-out observations
    // -------------------------------------------------------------------------

    if (N_test > 0) {
      vector[N_test] eta_test =
        Intercept +
        subcounty_effect[subcounty_id_test] +
        log_offset_test +
        X_test[, 2:K] * b;

      for (n in 1:N_test) {
        mu_hat_test[n] = exp(eta_test[n]);

        Y_rep_test[n] =
          neg_binomial_2_log_rng(
            eta_test[n], shape
          );

        log_lik_test[n] =
          neg_binomial_2_log_lpmf(
            Y_test[n] | eta_test[n], shape
          );

        log_lik_test_month[month_id_test[n]] +=
          log_lik_test[n];
      }
    }

    // -------------------------------------------------------------------------
    // Right-censored held-out observations
    // -------------------------------------------------------------------------

    if (N_test_cens > 0) {
      vector[N_test_cens] eta_test_cens =
        Intercept +
        subcounty_effect[subcounty_id_test_cens] +
        log_offset_test_cens +
        X_test_cens[, 2:K] * b;

      for (n in 1:N_test_cens) {
        mu_hat_test_cens[n] =
          exp(eta_test_cens[n]);

        Y_rep_test_cens[n] =
          neg_binomial_2_log_rng(
            eta_test_cens[n], shape
          );

        log_lik_test_cens[n] =
          neg_binomial_2_lccdf(
            Y_test_cens[n] - 1 |
            exp(eta_test_cens[n]),
            shape
          );

        log_lik_test_month[
          month_id_test_cens[n]
        ] += log_lik_test_cens[n];
      }
    }
  }
}

```

I’m using Pathfinder for initialization as

`init <- if (prior_only) NULL else model$pathfinder(data = sd, seed = 124, num_paths = 4, output_dir = "cmdstan_csv")`

The error arises from the right-censored likelihood. The call to `neg_binomial_2_lccdf(Y_cens[n] - 1 | exp(eta_cens[n]), shape)` either evaluates the probability as numerically zero during initialization or causes an overflow in the internal `hypergeometric_pFq` calculation during sampling. Pathfinder finds a finite initial point, but Chain 1 subsequently fails and CmdStanR cannot retrieve the fit metadata.

`Rejecting initial value:`  
`Log probability evaluates to log(0), i.e. negative infinity.`  
`Stan can't start sampling from this initial value.`  
`Path [3] :Initial log joint density = -571168.607012`  
`Chain 1 Exception: Error in function boost::math::hypergeometric_pFq(double,double,double): numeric overflow (in '/var/folders/qf/bzkp2r5d0zz_5lxfrbnc1wdh0000gn/T/Rtmplay7A8/model-7de54f298a8e.stan', line 124, column 8 to line 126, column 10)`

Is there a way to resolve this? Thanks in advance!

---

<div class="post-metadata">

**Author:** ![LucC](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/lucc/32/4475_2.png) [@LucC](https://discourse.mc-stan.org/u/LucC)\
**Post date:** [October 6, 2026, 5:34am UTC](https://discourse.mc-stan.org/t/initialization-failure-in-negative-binomial-model-with-exact-and-right-censored-likelihood-contributions/41644/2 "2026-10-06T05:34:48Z")

</div>

A similar underflow seems to have been reported and solved (by @sakrejda) for left-censored data and the **[`neg_binomial_2_lcdf`](https://mc-stan.org/docs/functions-reference/functions_index.html#neg_binomial_2_lcdf "Jump to index entry")** function:

> [@\`neg\_binomial\_2\_lcdf()\` "grad\_2F1: k (internal counter) exceeded 100000 iterations"](https://discourse.mc-stan.org/t/neg-binomial-2-lcdf-grad-2f1-k-internal-counter-exceeded-100000-iterations/8723):
>
> 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 grea…

Not sure if it matters or whether the fix was also applied to the lccdf function at the time as well?

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [October 6, 2026, 7:20pm UTC](https://discourse.mc-stan.org/t/initialization-failure-in-negative-binomial-model-with-exact-and-right-censored-likelihood-contributions/41644/3 "2026-10-06T19:20:07Z")

</div>

Sorry about that. Our cdfs and ccdfs are generally not very stable as most of them work on the linear scale and take a log because nobody’s worked out robust and efficient expansions of them. So this is mostly doing what @LuuC suggests. If there’s something that should be patched in Stan, you can open an issue in `stan-dev/math` and we’ll try to fix it.

How did you conclude that Pathfinder found a finite initial point? I ask because it looks like the sampler is rejecting the initial value.

This isn’t relevant for fixing your initialization problems, but you can make a lot of this more efficient by vectorizing operations. For example, this

```stan
    if (N_cens > 0) {
      vector[N_cens] eta_cens =
        Intercept +
        subcounty_effect[subcounty_id_cens] +
        log_offset_cens +
        X_cens[, 2:K] * b;

      for (n in 1:N_cens) {
        target += neg_binomial_2_lccdf(
          Y_cens[n] - 1 |
          exp(eta_cens[n]),
          shape
        );
      }
    }

```

can just be

```stan
transformed data {
  matrix[N_cens, K - 1] X_cens_2 = X_cens[, 2:K];
  array[N_cens] Y_cens_minus_1 = Y_cens;
  for (n in 1:N_cens) {
    y_cens_minus_1[n] -= 1;
  }
}

vector[N_cens] eta_cens
   = Intercept
     + subcounty_effect[subcounty_id_cens]
     + log_offset_cens
     + X_cens_2 * b;
target += neg_binomial_2_lccdf(
    Y_cens_minus_1 |
    Intercept + subcounty_effect[subcounty_id_cens] + log_offset_cens + ba * X_cens_2K,
    shape);

```

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [October 10, 2026, 5:49pm UTC](https://discourse.mc-stan.org/t/initialization-failure-in-negative-binomial-model-with-exact-and-right-censored-likelihood-contributions/41644/4 "2026-10-10T17:49:57Z")

</div>

This PR [Evaluate the incomplete beta complement without cancellation in the lccdf and Student t cdf functions by avehtari · Pull Request #3433 · stan-dev/math · GitHub](https://github.com/stan-dev/math/pull/3433) might fix this, at least it is improving `neg_binomial_2_lccdf()` a lot. Hopefully merged before the next release.

> [@Bob\_Carpenter](#):
>
> Our cdfs and ccdfs are generally not very stable as most of them work on the linear scale and take a log because nobody’s worked out robust and efficient expansions of them

There are now [10 `stan-dev/math` PRs greatly improving our cdfs, lcdfs and lccdfs](https://github.com/stan-dev/math/pulls?q=is%3Apr+state%3Aopen+label%3Anumerics)

---

<div class="post-metadata">

**Author:** ![LucC](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/lucc/32/4475_2.png) [@LucC](https://discourse.mc-stan.org/u/LucC)\
**Post date:** [October 21, 2026, 1:57pm UTC](https://discourse.mc-stan.org/t/initialization-failure-in-negative-binomial-model-with-exact-and-right-censored-likelihood-contributions/41644/5 "2026-10-21T13:57:01Z")

</div>

While waiting for these updates to be merged and released, is there anything that could be manually coded as a custom function. I was looking at a [PR](https://github.com/stan-dev/math/issues/2031) about the NB lccdf function, but that is only for the special case when the data (count) equals zero.

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [October 22, 2026, 4:52pm UTC](https://discourse.mc-stan.org/t/initialization-failure-in-negative-binomial-model-with-exact-and-right-censored-likelihood-contributions/41644/6 "2026-10-22T16:52:52Z")

</div>

> [@LucC](#):
>
> While waiting for these updates to be merged and released, is there anything that could be manually coded as a custom function.

Probably yes, but it is possible that even if the value would be finite, the autodiffed gradient might be non-finite. Many of the numerical stability fixes have been specifically for gradients and those can be fixed only in C++. If you need these before next release, it is possible to build CmdStan using specific github branches and get these updates (getting updates from different PRs before they are merged is a bit more complicated)