[avehtari]

Part of #56. #3408 fixes the remaining part, poisson_lccdf.

Merge order: please merge this PR after #3408, #3409, #3489, #3490, #3433, #3434, #3437, #3441, #3445 and #3470. Seven of them change 16 of the 18 headers in this PR. #3490 changes the four neg_binomial cdf headers, and #3433 changes binomial_lccdf and neg_binomial_lccdf. #3434 changes beta_binomial_lccdf and the beta_neg_binomial and yule_simon cdf headers. #3437 and #3445 also change beta_binomial_lccdf, and #3441 also changes the three beta_neg_binomial headers. #3408, #3409 and #3489 change none of these headers. #3490 and #3433 are stacked on #3409 and #3489. #3470 removes unused #include lines from 14 of the 18 headers, so its conflicts are in the include blocks. One of these includes must stay: #3470 removes select.hpp from bernoulli_lccdf.hpp, because develop does not use select there, but this PR does. I keep this PR as a draft until those PRs are merged. Then I merge develop into this branch, resolve the conflicts and run the tests again.

This PR was Claude-assisted.

The scripts behind the numbers in this text are in https://github.com/avehtari/stan-pr-checks/tree/6abd778ab2f5ab0f55f3a68b34015af34f1da285/pr3496. Its README.md quotes the claims of this text in sections 1 to 8, with the commands and the expected output. The quoted values come from probes built with g++ 11.4.0.

Summary

For a vector argument, a cdf returns the product of the element probabilities, and an lcdf or lccdf returns the sum of their logs. The functions in the table below had one or both of two defects. One element with probability 1 made the whole result 1 (cdf) or 0 (lcdf, lccdf), and the function ignored the other elements. Example: neg_binomial_lccdf({3, -1} | 2.5, 1.5) returned 0, but the correct value is neg_binomial_lccdf(3 | 2.5, 1.5) = -2.04. An element with probability 0 gave the correct -inf. The gradient was not 0, because the function returned after it had set the partials of earlier elements. In the table, INT_MAX is std::numeric_limits<int>::max(). On develop, the discrete cdf functions already treat a count n = INT_MAX as an infinite count. For such an element, the cdf is 1, the lcdf is 0 and the lccdf is -inf. This PR keeps that convention.

function	element	before	after
bernoulli_lccdf (CPU and OpenCL)	n < 0	0	skipped
binomial_lccdf	n < 0	0	skipped
binomial_lccdf	n >= N	-inf, gradient not 0	-inf, gradient 0
beta_binomial_lccdf	n < 0	0	skipped
beta_binomial_lccdf	n >= N	-inf, gradient not 0	-inf, gradient 0
neg_binomial_cdf, neg_binomial_2_cdf	n = INT_MAX	1	skipped
neg_binomial_lcdf	n = INT_MAX	0	skipped
neg_binomial_2_lcdf	n = INT_MAX	gradient not accurate	skipped
neg_binomial_lccdf (also used by neg_binomial_2_lccdf)	n < 0	0	skipped
neg_binomial_lccdf	n = INT_MAX	-inf, gradient not 0	-inf, gradient 0
beta_neg_binomial_cdf, beta_neg_binomial_lcdf	n = INT_MAX	1 or 0	skipped
beta_neg_binomial_lccdf	n < 0	0	skipped
beta_neg_binomial_lccdf	n = INT_MAX	-inf, gradient not 0	-inf, gradient 0
yule_simon_cdf, yule_simon_lcdf	n = INT_MAX	1 or 0	skipped
yule_simon_lccdf	n < 1	0	skipped
poisson_binomial_lccdf	y = N	finite, or -inf with gradient not 0	-inf, gradient 0

poisson_binomial_cdf and poisson_binomial_lcdf now also skip an element with y = N. Before, such an element gave a factor 1 (cdf) or a term 0 (lcdf), with a rounding error of about 1e-16.

All functions in the table now use one rule:

• A first pass looks for an element with probability 0. If it finds one, the function returns 0 or -inf before it sets any partial.
• The main loop skips each element with probability 1.

Three functions needed more than the rule:

• neg_binomial_2_lcdf calls beta_lcdf with the second shape n + 1. For n = INT_MAX that shape is 2^31, and the gradient through beta_lcdf is not accurate there. For neg_binomial_2_lcdf({3, INT_MAX} | 4, 2) on develop, the phi derivative is 19.41; the correct value, the phi derivative of neg_binomial_2_lcdf(3 | 4, 2), is -0.0604. The function now removes the elements with n = INT_MAX before the call. Large finite counts still go through beta_lcdf, and this PR does not change that path. There the phi derivative is also wrong. At mu = 4, phi = 2, it is 11.8 for n = 10^6 and 18.7 for n = 10^9. The correct derivative is below 1e-176000. At the head of Evaluate the incomplete beta complement without cancellation in the lccdf and Student t cdf functions #3433, which is stacked on More stable incomplete beta gradient roots #3409, Reverse-mode grad_reg_inc_beta (part 2 of #3409) #3489 and Move the beta, neg_binomial and student_t cdfs to the new grad_reg_inc_beta (part 3 of #3409) #3490, the derivative becomes 0.
• poisson_binomial_lccdf at y = N returned log1m_exp of a rounded sum near 0. For theta = (0.5, 0.2, 0.7) and y = 3 it returned -36.74, which is log(2^-53). The existing test prim/prob/poisson_binomial_ccdf_log_test expected that value. This PR changes the expectation to -inf, as the other y = N cases in that file expect.
• neg_binomial_lcdf computes digamma(n + 1) for each element before the main loop. With n = INT_MAX, n + 1 overflows int (UBSan reports it on develop); this PR uses n + 1.0.

The PR does not change the discrete cdf functions that already followed the rule: bernoulli_cdf and _lcdf, binomial_cdf and _lcdf, beta_binomial_cdf and _lcdf, and discrete_range.

Testing

• The new test test/unit/math/mix/prob/discrete_cdfs_edge_elements_test.cpp covers the cdf, lcdf and lccdf functions of the 8 families. An edge element is an element with probability 0 or 1. For an edge element with probability 1, the test compares f({ordinary, edge}) and f({edge, ordinary}) with f({ordinary}). It compares value and gradient, in reverse and forward mode. For an edge element with probability 0, it checks the value 0 or -inf and a zero gradient. Other cases use vector parameters, a scalar edge count with vector parameters, or (for poisson_binomial) a separate theta vector per element. With vector parameters, the gradient with respect to the parameters of the edge element must be exactly 0. expect_ad checks the new neg_binomial_2_lcdf code with scalar and vector mu and phi. All 11 test cases pass with the change and fail with the develop headers (503 failed expectations). For this test, I replaced all 17 changed prim headers with their develop versions at once.
• The existing OpenCL test opencl_matches_cpu_small_n_negative of bernoulli_lccdf (counts {0, 0, -5}) compares OpenCL with the fixed CPU function. It passes on a V100. With the develop OpenCL file and the fixed CPU file, it fails. Two new OpenCL tests check broadcasting in bernoulli_lccdf: one with a scalar n = -1, one with a scalar theta and n = {0, -1, 0}. They compare OpenCL with OpenCL and also pass on develop.

The new test and 35 existing unit suites of the 8 families (prim, rev, mix) pass. The generated test/prob tests of the 7 families that have them (all except poisson_binomial) pass. cpplint, make test-math-dependencies, a compile check of each changed header, make test-headers and clang-format 10 pass. make doxygen with doxygen 1.9.7 gives no warnings. I ran the tests with the change on develop at dc66941969, with g++ 14.3.0 (CPU) and a V100 (OpenCL).

Release notes

Discrete cdf, lcdf and lccdf functions now skip a vector element with probability 1; before, such an element made the result 1 or 0. An element with probability 0 now gives 0 or -inf with a zero gradient.

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

  • new 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