Fix and speed up the terminating 3F2 series of the beta_binomial cdfs - #3437

Fixes #2508.

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

This PR improves stability of beta_binomial_cdf, beta_binomial_lcdf and beta_binomial_lccdf

This was assisted by Claude

Summary

beta_binomial_cdf, beta_binomial_lcdf and beta_binomial_lccdf return NaN or throw at ordinary arguments. The comment by @wds15 in #2508 gives beta_binomial_lcdf(3 | 10, 298.30333970, 146.75306521): develop and #3434 return NaN, and the exact value is -3.9333500266686489. beta_binomial_lcdf(8 | 10, 3, 1) throws; an existing test expected this throw, with a FIXME that the point should be defined.

The three functions evaluate F = 3F2({1, mu, 1 - (N - n)}; {n + 2, 1 - nu}; 1) with mu = alpha + n + 1 and nu = beta + N - n - 1, and the derivatives of F with grad_F32. Because 1 - (N - n) is a non-positive integer, F is a polynomial with N - n terms. #3434 adds the mirrored series for the lcdf, which is a polynomial of the same kind.

Most of the code changes are in

In addition there are 2-6 lines changed per distribution file (except in beta_binomial_cdf also removing 90+ lines). Tests add 471 lines.

1. Four defects at the end of a terminating series.

1. hypergeometric_3F2_infsum (used at z = 1 when sum(b) < sum(a)) did not stop when a numerator parameter reached zero. It set the sign of the term to 0 and continued. The term magnitudes then overflowed, and 0 * inf gave NaN; or they did not overflow, and the loop ran up to max_steps = 1e5 steps. This is the case of beta_binomial_lpmf is unstable for large shape arguments #2508: 3F2({1, 302.30333970, -6}; {5, -151.75306521}; 1), exact 17.629673941225640. The loop now returns the sum when a numerator parameter is zero.
2. At z = 1 with sum(b) == sum(a), hypergeometric_3F2 called Boost's hypergeometric_pFq, which throws there. In the cdfs this is alpha + beta = 1. These cases now also use hypergeometric_3F2_infsum.
3. check_3F2_converges treated a denominator parameter b = -m as a pole when the polynomial ends at the term k = m. (b)_k is zero only for k > m, so the polynomial is defined. In the cdfs this is beta = 1. For example, 3F2({1, 12, -1}; {10, -1}; 1) is 1 + 1.2 = 2.2; develop threw. The Mathematica gradient for this case in the comment of grad_F32_test.cpp is already finite.
4. After the repair of 3, grad_F32 computed the term ratio 0 / 0 at the step where a numerator and a denominator parameter both reach zero. It now returns before the ratio when a numerator parameter is zero.

2. Smaller defects of the same code. check_3F2_converges took the end of a polynomial from the largest non-positive integer numerator instead of the smallest, so it threw for some defined polynomials (3F2({-1, -5, 1}; {-2, 1}; z)); and for a series that does not terminate it counted only b = 0 as a pole, not b = -1, -2, .... hypergeometric_3F2_infsum returned its partial sum without an error when the steps ran out (the check k == max_steps after the loop while (k <= max_steps ...) was never true), and it omitted the sign of z from every term after the first (no caller in Stan Math uses z < 0).

3. A stop that cannot cut the series. hypergeometric_3F2_infsum and grad_F32 stopped at a term below 1e-6 in absolute value. This cut the polynomial before its end: F had relative errors up to 1e-6 and its derivatives up to 3e-5. A plain relative stop is not safe either. In 3F2({1e-23, 300, -100}; {1, -150.5}; 1) all terms are positive, the term after 1 is 2e-21, and the sum is 5.8e36; both stops return 1, and grad_F32 with only the derivatives in a2 and b2 requested stops in the same way. A terminating series now stops only when the current term and a bound of the sum of all remaining terms are below 1e-17 of the partial sum (for grad_F32: also the bounds of every requested gradient). The bound is in a new internal header: the term ratio is a product of three factors |a_i + j| / |d_l + j|; where neither parameter changes sign, such a factor is monotone in j, so its maximum over the remaining terms is at an end of the range, and geometric sums then bound the tail. It is computed only at steps where the current term is already negligible. For an infinite series grad_F32 keeps its precision argument; the documentation says so.

4. The lcdf and lccdf. They ask grad_F32 for the two gradients they use (dF/da2 and dF/db2) instead of all six, which provides significant speed-up. #3434's lcdf kept log(1 - C) for alpha = 1, because its mirrored series then has b2 = a3 and threw by defect 3. With the repair, the lcdf uses the mirror for every alpha.

5. The cdf. beta_binomial_cdf multiplied the values 1 - C. Where the cdf is small this keeps no digits, and below about 1e-16 it rounds to 0, so the gradient (which divides by it) was NaN. It now returns exp(beta_binomial_lcdf(...)), after its own argument checks; the file shrinks from 159 to 66 lines.

Over a grid of 3801 points per function (N from 1 to 117, all n, shapes from 0.1 to 1e4), values and both shape gradients, against mpmath (sums of the pmf). Each cell: throws / NaN or inf / worst absolute error of the value or a shape gradient over the finite points.

No point that was finite on #3434 is less accurate with this PR.

Speed. Time per call with var shapes (value and both gradients), at N in {10, 100, 1000, 10000}, n in {0.1, 0.5, 0.9} N, and 5 shape pairs from (0.5, 2) to (1000, 3000) (180 points; Xeon E5-2680 v4; two interleaved rounds). At 95 of the 180 points #3434 runs the series past its end (defect 1), up to 1e5 steps: there a call takes 7.1 ms to 15.5 ms, and with this PR 0.5 µs to 1.6 ms. At the other 85 points, geometric mean of the ratio to #3434 [min, max]:

Per change (geometric means per N): the termination repair costs 2 % to 7 % from N = 100; the stop with the tail bound costs 3 % at N = 10 and 16 % to 25 % from N = 100 (most calls already summed nearly the whole polynomial; the largest increases, up to 3.4 times, are short series at shapes (1000, 3000)); asking for two gradients saves 14 % to 31 % in the lcdf and 24 % to 43 % in the lccdf; computing the cdf from the lcdf saves 21 % to 56 % in the cdf. Binaries that differ only in an unchanged function differ by up to about 6 % (code layout). The slowest point relative to #3434 is 1.7 times slower (14.9 µs to 25.8 µs, lccdf(5000 | 10000, 1000, 3000)).

Limits that remain

• hypergeometric_3F2_infsum now throws when a polynomial needs more than max_steps = 1e5 terms, as grad_F32 already did; before, it returned a truncated sum. In the cdfs this needs N - n above 1e5 and terms that do not become negligible.
• For an infinite series, grad_F32 still stops at a term below its precision argument (the beta_neg_binomial cdfs pass 1e-8). This PR does not change the beta_neg_binomial cdfs.

Tests

• test/unit/math/prim/fun/hypergeometric_3F2_test.cpp: {1, 12, -1}; {10, -1}; 1 is now the value 2.2 instead of a throw; the pole test uses {1, 12, -2}; {10, -1}; 1. New tests against exact finite sums: the series of beta_binomial_lpmf is unstable for large shape arguments #2508, a zero numerator parameter, b2 = a3 = -6, equal parameter sums at z = 1, the smallest of two non-positive integer numerators, z = -0.5 and the step limit of hypergeometric_3F2_infsum, a 117-term series whose terms fall below 1e-6 before its end, and the series with a tiny term before large terms.
• test/unit/math/prim/fun/grad_F32_test.cpp: the same 2.2 case as gradients (the Mathematica values from the file comment); the pole test uses a3 = -2. New tests: b2 = a3 = -6, the series of beta_binomial_lpmf is unstable for large shape arguments #2508, the smallest numerator, the 117-term series, and the tiny-term series with only two gradients requested.
• test/unit/math/prim/err/check_3F2_converges_test.cpp: the smallest numerator ends a polynomial; non-positive integer denominators of an infinite series are poles.
• test/unit/math/prim/fun/hypergeometric_3F2_tail_bound_test.cpp (new): the bound is not below the largest term ratio for 6 series with sign changes and a negative z, equals it for a cdf series, and the tail sums bound the exact sums.
• test/unit/math/prim/prob/beta_binomial_cdf_log_test.cpp: lcdf_matches_mathematica now checks the value log(15 / 26) instead of the throw.
• test/unit/math/rev/prob/beta_binomial_cdfs_test.cpp (new): values and shape gradients of the three cdfs at 7 points (the point of beta_binomial_lpmf is unstable for large shape arguments #2508, beta = 1 twice, alpha + beta = 1, N - n = 1, a 117-term series and its mirror); the lcdf at two points with alpha = 1 and a small cdf; the cdf at three points where it is below 1e-20. Relative tolerance 1e-12, against sums of the pmf in mpmath.

The tests find the defects. With the seven headers of #3434 swapped in, the 23 new or changed test cases fail and all other cases pass. Each commit's new tests fail with the headers of the previous commit.

Full suites. All unit tests of hypergeometric_3F2, the new tail bound, grad_F32, grad_pFq, check_3F2_converges, hypergeometric_pFq, inv_inc_beta, beta_binomial and beta_neg_binomial pass, and so does #3434's large_shapes_test (20 files: prim, rev, fwd, mix), and the generated test/prob tests of beta_binomial and beta_neg_binomial (31 600 and 69 498 cases).

Release notes

Improved stability of beta_binomial_cdf, beta_binomial_lcdf and beta_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