# 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>
digammaN_vec(size(alpha));
VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
digammaAlpha_vec(size(alpha));
VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
digammaSum_vec(size(alpha));

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.

5 Likes

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.

2 Likes

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 https://github.com/stan-dev/math/issues/1531

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 https://github.com/stan-dev/math/issues/1662 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
9 Likes

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

3 Likes

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.