Suspicious code in neg_binomial_lccdf

I’m looking at stan/math/prim/prob/neg_binomial_lccdf, and I found this bit of code:

  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>

  if (!is_constant_all<T_shape>::value) {
    for (size_t i = 0; i < size(alpha); i++) {
      const T_partials_return n_dbl = value_of(n_vec[i]);
      const T_partials_return alpha_dbl = value_of(alpha_vec[i]);

      digammaN_vec[i] = digamma(n_dbl + 1);
      digammaAlpha_vec[i] = digamma(alpha_dbl);
      digammaSum_vec[i] = digamma(n_dbl + alpha_dbl + 1);

I find it strange that we loop over the size of alpha but inside the loop we read n_vec: if alpha_vec is a scalar and n_vec is a vector, wouldn’t we be computing the wrong quantity (and ignoring a number of values from n_vec)?

Also the other *cdf functions for neg_binomial use the same pattern.


That does look suspicious.

I’ve been using expose_stan_functions in R to play about with this, and results look correct… but I don’t know why! Can anybody think of an example that shows it’s doing the wrong thing, or explain why it’s actually fine? @syclik, @Bob_Carpenter, @nhuurre, @martinmodrak

Does expose_stan_functions let you see the derivatives? The problematic code only affects the gradient.


Good catch, back to the c++ side!

It looks like we have a lot of issues with sizing, the code looks problematic, similar issue was reported in

Thanks @martinmodrak, I had forgotten about that (I thought our problems in these functions were limited to the numerics). Also neg_binomial_2_cdf has the same bug, while neg_binomial_2_lccdf is also affected since it delegates to neg_binomial_lccdf. On the other hand, neg_binomial_2_lcdf seems to be fine.

I’ve opened to track this.

1 Like

Now that the issue has been resolved, let me summarise what we found out. If functions are called with inputs of different lengths:

  • neg_binomial_cdf , neg_binomial_lcdf , neg_binomial_lccdf , neg_binomial_2_cdf and neg_binomial_2_lccdf had wrong derivatives (not anymore!)
  • neg_binomial_lpmf , neg_binomial_2_lpmf , neg_binomial_2_lcdf and neg_binomial_2_log_lpmf were fine

Fantastic! Thanks Marco for handling all these rough edges in the Math library. We are really lucky you joined the project!


Nicely done! I did some of this bug-chasing a long time ago and I still think it’s really key to making Stan the awesome tool that it is.

I agree! I’m hoping improved automatic tests will help with the combinatorics here and find this kind of bug in the future.