Compute the beta_neg_binomial cdfs from sums of the pmf

[avehtari]

This PR is based on #3434 and must be merged after it.

This PR improves the accuracy and speed of beta_neg_binomial_cdf, beta_neg_binomial_lcdf and beta_neg_binomial_lccdf. The speed-up is from 4x to 9e5x better.

This was assisted by Claude

Summary

beta_neg_binomial_cdf, beta_neg_binomial_lcdf and beta_neg_binomial_lccdf computed the upper tail as pmf(n + 1) * 3F2(1, beta + n + 1, r + n + 1; n + 2, alpha + r + beta + n + 1; 1). This series has the parameter excess alpha, so its terms fall only like k^(-1 - alpha). On a grid of 2352 points (alpha 0.25 to 1e4, r and beta 0.5 to 1e6, n 0 to 20000), the three functions threw at 56 to 57 % of the points and returned a wrong value or gradient at another 22 to 28 %. Only 16 to 21 % of the points were correct, and a call took 15 ms (median) and up to 14 s. References: mpmath, the cdf as the finite sum of the pmf, the ccdf as 1 - cdf with extra digits or as the sum of the pmf over k > n.

The causes: the value of the 3F2 overflows for shapes from 300; grad_F32 runs into its limit of 1e8 steps (8 to 13 s, then a throw); Boost hypergeometric_pFq returns the series wrong by a factor up to e^24 at small alpha; grad_F32 stops at a term below 1e-8, while the rest of such a series is about k t_k / alpha (gradient errors of 1e-5 even at alpha = 10); and cdf and lcdf formed 1 - ccdf, so lcdf was wrong by up to 490 where the cdf is small.

This PR computes one tail as a sum and the other one as its complement, so that no result cancels:

• the lower tail as the sum of the pmf from n down to 0, in log scale, with the gradient as the sum of pmf * d log pmf; the pmf is unimodal, so below the mode the sum stops when the remaining terms are negligible;
• the upper tail, when it is small, as the sum of the pmf over k > n when its terms fall fast, or by a Thomae transformation of the 3F2, whose new series has the excess x + n + 1 (x is r or beta) instead of alpha; for an integer x the new series is a polynomial, and its derivative with respect to x continues past the zero factor;
• the complement only where the other tail is at least 1e-2, so that it loses at most 2 digits;
• beta_neg_binomial_cdf as the exponential of the sum of the log cdfs.

Other libraries that have this cdf (SciPy betanbinom, R extraDistr::pbnbinom, the generic discrete logcdf of PyMC) sum the pmf for the lower tail and use 1 - cdf for the upper tail; that cancels where the ccdf is small. This PR adds the direct upper-tail methods for that region.

on the grid of 2352 points	before	this PR
correct points, cdf / lcdf / lccdf	501 / 379 / 376	2352 / 2352 / 2352
throws	1315 / 1341 / 1315	0
worst value error	lcdf wrong by 490	3.6e-10 relative (cdf), 4.7e-10 absolute (lcdf at -1.4e6)
worst log-shape gradient error	6.7e5 (lcdf)	9.6e-10 absolute

"Correct" means that the value and the gradients with respect to log(r), log(alpha), log(beta) are within 1e-8, absolutely or relatively.

Time per call with var arguments (value and three gradients), Xeon E5-2680 v3, median over 30 points per region, two interleaved rounds:

region	before	this PR	before / this PR
ordinary (alpha 2 to 10, shapes up to 300)	32 ms	1.3 us	6700
small alpha (0.25 to 1)	286 ms	0.7 us	2.9e5
n = 20000	31 ms	1.2 us	3700
shapes from 3000 (19 of 30 throw before)	60 ms	1.0 us	2900
alpha >= 100	3.3 us	1.4 us	4.1
8 points where a call took more than 1 s before (6 throw)	2.1 s	1.3 us	8.9e5

The slowest of the 158 timed points takes 0.72 ms (a lower sum of 2e4 terms below the mode); before, that call threw after 13.9 s. For the upper tail, the cheaper of the tail sum and the Thomae series runs first, by an estimate of their numbers of terms.

With #3434 alone, the cdf gradient is NaN at 15 grid points where the cdf is below 1e-15: the more accurate log pmf makes exp(C) F round to exactly 1, so 1 - ccdf = 0. This PR does not form 1 - ccdf there.

Two smaller changes:

• n equal to the largest int stands for infinity. beta_neg_binomial_cdf and beta_neg_binomial_lcdf returned 1 and 0 for the whole vector when one element had it; now that element contributes a factor of 1 (a term of 0).
• The optional precision argument is no longer used (the sums stop when their remaining terms are below 1e-17 of the partial sums). max_steps is now the largest number of terms of a sum; a function throws if a sum needs more. No call on the grid comes near it.

Limits that remain

• At n = 20000, r = 0.5, beta = 1e6 and alpha <= 1 the log pmf at n + 1 has an absolute error of about 2e-11, because its terms (for example lbeta(20001, 1e6), of size 1e5) cancel to about -15. lccdf inherits this, and cdf as its complement multiplies it by 17: 3.6e-10 relative. beta_neg_binomial_lpmf has the same error there.
• The time per call grows with the number of terms: a lower sum below the mode of a wide distribution, or a slowly falling tail, can take up to about 1 ms.

Tests

• New test/unit/math/rev/prob/beta_neg_binomial_cdfs_test.cpp: values and log-scale gradients of the three functions at 15 points that cover every route (lower sum, both complements, short and long tail sum, the Thomae series with a non-integer and an integer x), including small alpha, shapes of 1e6, a cdf of 2.6e-35 and n = 20000; tolerances 1e-11 (value) and 1e-10 (gradient), relative to max(1, |reference|). A second test checks vectors and the largest int. With the previous headers both tests fail.
• test/unit/math/rev/prob/large_shapes_test.cpp: the beta_neg_binomial_cdf and _lccdf rows now use the default tolerances (1e-12, 1e-11) instead of (1e-6, 1e-4); with the previous headers they fail.
• The generated tests of test/prob/beta_neg_binomial pass (69 498 cases).

Side Effects

None.

Release notes

Improved the accuracy and speed of beta_neg_binomial_cdf, beta_neg_binomial_lcdf and beta_neg_binomial_lccdf

Checklist

• Copyright holder: Aalto University

  The copyright holder is typically you or your assignee, such as a university or company. By submitting this pull request, the copyright holder is agreeing to the license the submitted work under the following licenses:
  - Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
  - Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)

• the basic tests are passing

  • unit tests pass (to run, use: ./runTests.py test/unit)
  • header checks pass, (make test-headers)
  • dependencies checks pass, (make test-math-dependencies)
  • docs build, (make doxygen)
  • code passes the built in C++ standards checks (make cpplint)

• the code is written in idiomatic C++ and changes are documented in the doxygen

• the new changes are tested