Numerical issues in neg_binomial_lccdf if beta is very large #2031

Description

When the second argument to neg_binomial_lccdf is very large, the function can return -Inf even if the answer is finite.

Example

We take advantage of the fact the neg_binomial_lccdf(0|a, b) == log(1 - exp(neg_binomial_lpmf(0| a, b))
The formula then simplifies to neg_binomial_lccdf(0|a, b) == log1m( (b / (b + 1)) ^ a). Now for large b b/(b +1) is numerically exactly 1, but we can further rewrite as neg_binomial_lccdf(0|a, b) == log1m_exp(-a * log1p(1/b)).

So consider the following Stan program:

transformed data {
  real alpha = 1;
  real beta = 1e18;
  print(neg_binomial_lccdf(0 | alpha, beta));
  print(log1m_exp(-alpha * log1p(1/beta)));
}

Which outputs

-inf
-41.4465

Expected Output

-41.4465
-41.4465

Note that the two versions start to diverge only at around beta > 1e14, so this is unlikely to be super important in practice.

Current Version:

v3.3.0