[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

[stan-buildbot]

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan	0.25	0.24	1.01	0.87% faster
stat_comp_benchmarks/benchmarks/gp_regr/gen_gp_data.stan	0.06	0.06	1.01	0.92% faster
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	6.56	6.58	1.0	-0.35% slower
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	1.03	2.5% faster
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	9.26	9.23	1.0	0.34% faster
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	4.7	4.66	1.01	0.7% faster
stat_comp_benchmarks/benchmarks/sir/sir.stan	229.97	229.76	1.0	0.09% faster
stat_comp_benchmarks/benchmarks/garch/garch.stan	1.01	1.01	1.0	0.3% faster
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.69	0.7	0.99	-0.73% slower
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	42.71	43.9	0.97	-2.79% slower
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.63	0.59	1.07	6.39% faster
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	22.22	22.15	1.0	0.31% faster
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.11	0.11	1.0	-0.26% slower
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.4	3.41	1.0	-0.23% slower
performance.compilation	374.73	367.82	1.02	1.85% faster
Mean result: 1.0070409640303366

Jenkins Console Log
Jenkins Build Stages
Commit hash: e7cab1646d71f360d212c6aecb8537bca67ebe6f

Machine information

Distributor ID:	Ubuntu
Description:	Ubuntu 20.04.3 LTS
Release:	20.04
Codename:	focal

CPU:

Architecture:                            x86_64
CPU op-mode(s):                          32-bit, 64-bit
Byte Order:                              Little Endian
Address sizes:                           43 bits physical, 48 bits virtual
CPU(s):                                  256
On-line CPU(s) list:                     0-255
Thread(s) per core:                      2
Core(s) per socket:                      64
Socket(s):                               2
NUMA node(s):                            2
Vendor ID:                               AuthenticAMD
CPU family:                              23
Model:                                   49
Model name:                              AMD EPYC 7742 64-Core Processor
Stepping:                                0
Frequency boost:                         enabled
CPU MHz:                                 1495.808
CPU max MHz:                             3416.0681
CPU min MHz:                             1500.0000
BogoMIPS:                                4491.85
Virtualization:                          AMD-V
L1d cache:                               4 MiB
L1i cache:                               4 MiB
L2 cache:                                64 MiB
L3 cache:                                512 MiB
NUMA node0 CPU(s):                       0-63,128-191
NUMA node1 CPU(s):                       64-127,192-255
Vulnerability Gather data sampling:      Not affected
Vulnerability Indirect target selection: Not affected
Vulnerability Itlb multihit:             Not affected
Vulnerability L1tf:                      Not affected
Vulnerability Mds:                       Not affected
Vulnerability Meltdown:                  Not affected
Vulnerability Mmio stale data:           Not affected
Vulnerability Old microcode:             Not affected
Vulnerability Reg file data sampling:    Not affected
Vulnerability Retbleed:                  Mitigation; untrained return thunk; SMT enabled with STIBP protection
Vulnerability Spec rstack overflow:      Mitigation; Safe RET
Vulnerability Spec store bypass:         Mitigation; Speculative Store Bypass disabled via prctl
Vulnerability Spectre v1:                Mitigation; usercopy/swapgs barriers and __user pointer sanitization
Vulnerability Spectre v2:                Mitigation; Retpolines; IBPB conditional; STIBP always-on; RSB filling; PBRSB-eIBRS Not affected; BHI Not affected
Vulnerability Srbds:                     Not affected
Vulnerability Tsa:                       Not affected
Vulnerability Tsx async abort:           Not affected
Vulnerability Vmscape:                   Mitigation; IBPB before exit to userspace
Flags:                                   fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush mmx fxsr sse sse2 ht syscall nx mmxext fxsr_opt pdpe1gb rdtscp lm constant_tsc rep_good nopl xtopology nonstop_tsc cpuid extd_apicid aperfmperf rapl pni pclmulqdq monitor ssse3 fma cx16 sse4_1 sse4_2 x2apic movbe popcnt aes xsave avx f16c rdrand lahf_lm cmp_legacy svm extapic cr8_legacy abm sse4a misalignsse 3dnowprefetch osvw ibs skinit wdt tce topoext perfctr_core perfctr_nb bpext perfctr_llc mwaitx cpb cat_l3 cdp_l3 hw_pstate ssbd mba ibrs ibpb stibp vmmcall fsgsbase bmi1 avx2 smep bmi2 cqm rdt_a rdseed adx smap clflushopt clwb sha_ni xsaveopt xsavec xgetbv1 xsaves cqm_llc cqm_occup_llc cqm_mbm_total cqm_mbm_local clzero irperf xsaveerptr rdpru wbnoinvd amd_ppin arat npt lbrv svm_lock nrip_save tsc_scale vmcb_clean flushbyasid decodeassists pausefilter pfthreshold avic v_vmsave_vmload vgif v_spec_ctrl umip rdpid overflow_recov succor smca sev sev_es

G++:

g++ (Ubuntu 9.4.0-1ubuntu1~20.04) 9.4.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

Clang:

clang version 10.0.0-4ubuntu1 
Target: x86_64-pc-linux-gnu
Thread model: posix
InstalledDir: /usr/bin