Compute the beta_binomial cdfs from sums of the pmf

[avehtari]

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

This PR improves the speed of beta_binomial_cdf, beta_binomial_lcdf and beta_binomial_lccdf. The speed-up is from 1.3x to 176x.

This was assisted by Claude

Summary

Motivation to speed-up comes from comments outside of gtihub issues. For example, in the forum thread Beta binomial truncated takes 1000x more time than the non truncated version (2021), a truncated model with beta_binomial_lcdf took 193 s, and the same model without truncation took 1.5 s.

With #3437, beta_binomial_lcdf and beta_binomial_lccdf evaluate a terminating 3F2 series with N - n terms (for the lcdf where the cdf is below 1/2, also the mirrored series with n + 1 terms), and its two shape derivatives with grad_F32; beta_binomial_cdf is exp(beta_binomial_lcdf(...)). These loops work in log scale with exp and log calls in every term, and most calls sum nearly the whole polynomial. A call with var shapes takes about 0.1 ms at N = 1000 and 0.9 ms at N = 10000. A series that needs more than 1e5 terms throws (the step limit of hypergeometric_3F2).

Most likely users use these cdfs this via brms. For the beta_binomial() family, brms generates beta_binomial_lccdf and beta_binomial_lcdf for cens() (R/stan-likelihood.R) and for trunc() (same file); the brms tests check these calls (tests.stancode.R). The zero_inflated_beta_binomial family builds its lcdf and lccdf on beta_binomial_lccdf (fun_zero_inflated_beta_binomial.stan).

This PR computes the cdfs as sums of the pmf, pmf(k) = choose(N, k) B(k + alpha, N - k + beta) / B(alpha, beta):

• the lower tail from k = n down to 0 and the upper tail from k = n + 1 up to N, in log scale with rescaling, by the term ratio pmf(k + 1) / pmf(k) = (N - k)(k + alpha) / ((k + 1)(N - k - 1 + beta)); the first term is the log pmf by Avoid the cancellation of lgamma, lbeta and digamma differences at large shape parameters #3434's internal::log_beta_ratio, which does not cancel at large shapes;
• the gradients as sums of pmf(k) d log pmf(k), with digamma_diff at the first term and the recurrences 1 / (k + alpha) and 1 / (N - k - 1 + beta) after it;
• a sum stops early only where all its remaining terms are smaller than the current one: pmf(k + 1) / pmf(k) - 1 has the sign of A + B k with A = N alpha - N + 1 - beta and B = 2 - alpha - beta, so the signs at the two ends of the remaining range show this. The sum then stops when a bound of the remaining terms is below 1e-17 of the partial sum. The pmf is unimodal for alpha + beta > 2 and U-shaped for alpha + beta < 2; a U-shaped pmf can need all terms of a tail;
• the tail with fewer terms (n + 1 below, N - n above) is summed first. If it is at most 1/2, the other tail is its complement (log1m_exp), which then loses no digits. Else the other tail is summed too, and if that one is at most 1/2, the first tail is replaced by its complement: a sum near 1 that starts at a tiny term keeps the rounding error of the log of that term (about 1e-13 for a log of -500).

No 3F2 function, grad_F32 or step limit is used. Most of the code changes are in

file	+	−	added code / comment
prim/prob/beta_binomial_log_cdfs.hpp (new)	272	0	203 / 57
prim/prob/beta_binomial_lcdf.hpp	13	90	9 / 4
prim/prob/beta_binomial_lccdf.hpp	19	63	11 / 8

In addition prim/prob.hpp gets one include line, and beta_binomial_cdf is unchanged. Tests add 157 lines.

Accuracy on the grid of #3437 (3801 points per function: N from 1 to 117, all n up to N = 10 and 10 values of n above, shapes from 0.1 to 1e4, and the point of #2508), against mpmath (sums of the pmf at 50 digits); before is #3437. Each cell: worst absolute error of the value / d/dalpha / d/dbeta:

function	before	this PR
beta_binomial_cdf	1.4e-14 / 1.1e-13 / 1.4e-13	9.4e-15 / 2.8e-14 / 2.8e-14
beta_binomial_lcdf	2.4e-13 / 2.9e-13 / 4.3e-13	2.4e-13 / 4.9e-14 / 5.6e-14
beta_binomial_lccdf	7.7e-13 / 8.9e-14 / 2.8e-12	2.4e-13 / 1.1e-13 / 1.1e-13

Neither version throws or returns NaN on the grid. The largest value errors of the lcdf and lccdf (2.4e-13) are at log values near -620, a relative error of 4e-16, in both versions. Per point and function (the largest error over the value and both gradients, counted where it is above 1e-13), this PR is more than 2 times more accurate at 276 of the 11 403 cases and more than 2 times less accurate at 6.

Beyond the grid:

• beta_binomial_lccdf(1000 | 10000, 100, 30): the ccdf is practically 1; the value is -1.8e-71 and the derivative with respect to alpha is 3.5e-71. Before, it returned -2.6e-12 and -4.0e-12; this PR returns both with a relative error of 6e-13 (the absolute error of its lcdf at -163).
• N - n above 1e5: before, the functions threw where the series does not stop early, for example at (n, N, alpha, beta) = (50000, 200000, 0.5, 0.5) for all three functions (the lcdf evaluates the N - n series before it chooses the mirror). This PR returns the values and gradients there (50 001 terms).

Time per call with var shapes (value and both gradients), Xeon E5-2680 v3, on the 60 points per function of the timing of #3437 (N in {10, 100, 1000, 10000}, n in {0.1, 0.5, 0.9} N, and the shape pairs (0.5, 2), (2, 3), (10, 10), (100, 30), (1000, 3000)), median of 5 repetitions per point, two interleaved rounds. Medians over the 45 calls per N (15 points, three functions); the last column is the geometric mean [min, max]:

N	before	this PR	before / this PR
10	1.56 us	0.533 us	2.5 [1.3, 4.6]
100	10.9 us	0.891 us	6.4 [2.9, 34]
1000	92.5 us	3.46 us	11 [4.7, 128]
10000	873 us	29.9 us	13 [4.1, 176]

The smallest speed-up is 1.3x (beta_binomial_lccdf(5 | 10, 100, 30): 1.56 us to 1.19 us). The time of this PR grows with the number of summed terms, about 30 ns per term; at these points that is mostly the shorter tail, 0.1 N terms.

One smaller change: beta_binomial_lccdf returned 0 for the whole vector when one element had n < 0; now that element contributes a term of 0 (log(1)), as an element with n >= N does in beta_binomial_lcdf.

Limits that remain

• The time per call grows linearly with the number of summed terms. A U-shaped pmf, or a wide pmf with n near its middle, needs most of the shorter tail: 50 000 terms take about 1 to 1.5 ms.
• Far in a tail, the lcdf and lccdf values have the absolute error of the log pmf, about 1e-16 |log pmf| (2.4e-13 at -620), as before.
• At large N the log pmf of the first term is a sum of terms of size up to 1e5 (for example binomial_coefficient_log(200000, 50000)), and every summed term inherits its rounding error. At (50000, 200000, 0.5, 0.5) that log pmf has an absolute error of 3.9e-12, and so has the lcdf value; beta_binomial_lpmf returns the same log pmf. The gradients agree with mpmath to 6e-15.

Tests

• New test/unit/math/rev/prob/beta_binomial_log_cdfs_test.cpp: values and log-scale gradients of the three functions at 15 points that cover every route (each tail first, each complement, both tails summed, and the complement of the second tail in both directions), the unimodal and the U-shaped pmf, an early stop at shapes (1000, 3000), the point of beta_binomial_lpmf is unstable for large shape arguments #2508, tails of 2e-276, N = 1, an lccdf gradient of 3.5e-71 at N = 10000, and N = 200000. References: mpmath at 50 digits, sums of the pmf. Tolerances: 1e-11 relative for the values, 1e-10 relative to max(1, |reference|) for the gradients. A second test checks vectors and the edge cases n < 0 and n >= N.
• With the previous headers both tests fail: the value check at 8 points where the lcdf or lccdf is between -7.8e-9 and -2e-276 (absolute errors from 1e-18 to 2.6e-12, above the relative tolerance), the lccdf gradient at N = 10000, the throw at N = 200000, and the n < 0 vector case.
• The existing unit tests of the beta_binomial family and large_shapes_test pass, and so do the generated tests of test/prob/beta_binomial (31 600 cases).

Side Effects

None.

Release notes

Improved the speed 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

[WardBrian]

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

Would you mind marking these as drafts when this is the case? It's difficult to tell from the PR listing what is ready to review at the moment...

[avehtari]

Would you mind marking these as drafts when this is the case? It's difficult to tell from the PR listing what is ready to review at the moment...

Sure. I had not thought some of these PRs as drafts although they are based on other PRs, but I can do that. I don't know will happen now that @SteveBronder stacked them. Claude did warn that using github stacking can cause difficulties.

[avehtari]

Actually now that @SteveBronder stacked this PRs, you can easily see in the PR listing a small stack symbol and 1/2 or 2/2, so all with 2/2 require that some other PR is merged first. @WardBrian, do you still want that I change all 2/2:s in stacks to drafts?