Suspicious code in neg_binomial_lccdf

Developers
1

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
2

That does look suspicious.

3

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

4

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

2 Likes
5

Good catch, back to the c++ side!

6

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

7

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
8

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
9

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

3 Likes
10

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.

11

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