wrong gradient computations in neg_binomial_*cdf functions

[mcol]

Description

All cumulative distribution functions for neg_binomial contain code such as the following:

math/stan/math/prim/prob/neg_binomial_lccdf.hpp

Lines 62 to 71 in dc096aa

	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);
	}
	}

Note that while the loop goes up to the size(alpha), the code also uses elements from n_vec. This means that if n is an array while alpha is a scalar, only the first element of n will be used, while the rest are going to be ignored. The quantities computed there are then used in gradients, so those will be wrong.

A similar bug is present also in neg_binomial_2_cdf, 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.

This was initially discussed in https://discourse.mc-stan.org/t/suspicious-code-in-neg-binomial-lccdf/12889.

Current Version:

v3.1.0