[avehtari]

EDIT: I updated this PR text on 2026-10-21 and on 2026-10-28, after the fact check by @stevieb-bot. A comment lists the changes.

This PR was Claude assisted

Issues

• Closes gamma_p and gamma_q gradient problems #2006: gamma_q with fvar and gamma_p with var return wrong derivatives (Part 2).
• Partly addresses Incorrect behaviour of CCDFs when there is a negative random variable #56: poisson_lccdf no longer drops the other elements of a vector when one count is negative. The same defect remains in other discrete cdfs; see "Not included".
• Related to Add Boost's tgamma_lower and tgamma(a,z) functions #2769: log P and log Q are now computed in log space inside the library, but no public function is added.
• Repairs the gamma_lccdf shape-gradient failures that @mjcarter95 reported in this thread (Part 3).

Summary

This PR repairs the derivatives of the regularized incomplete gamma functions and the 18 distribution functions that are built on them: the cdf, lcdf and lccdf of chi_square, gamma, inv_chi_square, inv_gamma, scaled_inv_chi_square and poisson. It has three parts.

1. The gradient roots grad_reg_inc_gamma and grad_reg_lower_inc_gamma (the derivatives of Q(a, z) and P(a, z) with respect to a) are replaced by Eigen::numext::igamma_der_a, which is already vendored in lib/eigen_5.0.1. The old roots were wrong by tens of percent, or NaN, in regions that ordinary models reach.
2. The forward-mode gamma_q had its own copy of the old series, which did not terminate for z above about 700. The copy is removed. The reverse-mode gamma_p returned zero derivatives for z / a > 10; it now returns the derivatives. Closes gamma_p and gamma_q gradient problems #2006.
3. The 18 cdf functions now compute log P or log Q and their derivatives in log space, through one internal function. Before, they returned -inf values or NaN gradients where the probability underflows, and 14 of them returned NaN or infinite gradients in the bulk once the shape is above about 150.

The values of the distributions change only where they were -inf or lost digits; elsewhere only the gradients change.

	develop	this PR
roots, non-finite results at 49 990 points	7 379	0
roots, worst relative error outside the band near z = a	inf	3e-11 (reference-limited)
roots, autodiff-order spread	inf in 7 of 10 regions	2.2e-14
forward-mode gamma_q, z ≳ 700	does not terminate	finite
18 cdfs, test points with a non-finite value or gradient, of 324	199	0
18 cdfs, worst relative gradient error where finite	1.6e+93	2.7e-10
speed, double roots	1x	2.2x to 4.4x
speed, cdf value plus gradients (var)	1x	1.9x to 12x

Part 1: the gradient roots

The defects on develop

All numbers in this part were measured against an mpmath reference at 80 digits, built through two independent routes that agree to at least 13 digits, with the real develop headers at 5252d51d47 compiled and run.

1. grad_reg_inc_gamma takes an invalid branch when z is near a. The test z >= a && z >= 8 selects a fixed 9-term asymptotic expansion (DLMF 8.11.2), which needs z >> a. In the whole band a <= z < 2a the terms do not decay. The relative error at z = a is 4.4e-02 at a = 20, 5.8e-01 at a = 100 and 7.2e-01 at a = 171.
2. grad_reg_inc_gamma returns NaN for a above about 171.6. Every call site passes g = tgamma(a), which is inf there, so the correction term exp(-z + (a-1) log z) * S / g is inf / inf. In the series branch, exp(a * log z) overflows when a log z > 709. The result is then -inf or +inf. At (a, z) = (170, 100) it is -inf; the true value is 6.6e-11.
3. grad_reg_lower_inc_gamma loses all accuracy for z well above a. Its delegation test sends every point with z > 42 and a >= 12 to the Gautschi branch, where emz * (log_z * sum_a - sum_b) is a difference of two large, nearly equal sums. It returns about 1e-14 of rounding noise with arbitrary sign: at (50, 150) the reference is -8.29e-22 and develop returns -4.70e-15; at (50, 500) the reference is -5.34e-148 and develop returns +1.47e-14.
4. The default precision of 1e-6 sets the error of the plain series branch: 9.2e-08 at (a, z) = (0.5, 0.25) and 9.0e-05 at (8, 4), with an exact g. It is the smallest of the four problems.

The defects reach the distribution functions. Measured with var against the same reference (the cdf values are correct; only the gradients are wrong):

function	point	reference	develop	rel. error
chi_square_lccdf	nu = 200	0.0410184911365	0.0173344546475	5.8e-01
chi_square_lccdf	nu = 342	0.0311564510080	0.0086722587902	7.2e-01
chi_square_lcdf	nu = 200	-0.0388930365624	-0.0164362354566	5.8e-01
gamma_cdf	alpha = 100	-0.0399274978578	-0.0168733998162	5.8e-01
inv_gamma_lccdf	alpha = 100	-0.0777860731247	-0.0328724709132	5.8e-01

The gradients are NaN for chi_square_lccdf and chi_square_lcdf at nu = 360, gamma_cdf and inv_gamma_lccdf at alpha = 180, and scaled_inv_chi_square_lcdf at nu = 360. With this PR, every row of the table agrees with the reference to all 12 printed digits. The rows of the table and the NaN points are in the reference test (rows marked root_table).

The change

Eigen::internal::igamma_generic_impl differentiates the Cephes power series and continued fraction term by term, stops at machine epsilon, and branches on x > 1 && x > a, which is the well-conditioned split. It is generic in the scalar type except for three helpers, lgamma_impl, digamma_impl and cephes_helper. This PR supplies those three for Stan's autodiff scalars in two new headers, stan/math/fwd/fun/Eigen_SpecialFunctions.hpp and stan/math/rev/fun/Eigen_SpecialFunctions.hpp. Both roots then become thin wrappers, and the hand-written series, the asymptotic branch and the Gautschi branch are deleted.

One algorithm now serves every scalar type and every autodiff order, so reverse mode, forward mode, fvar<fvar<double>>, fvar<var> and hessian() agree. I tried the double path alone first and rejected it: hessian() stayed on the old code, and gradient() and hessian() then returned first derivatives that differed by 5.8e-01 at alpha = 100 and by NaN at alpha = 180.

The signatures do not change. The roots ignore precision and max_steps, and check g and dig only for NaN. If one of them is NaN, the result is NaN, as before.

Accuracy

A sweep of 1 800 points in nine regions chosen by the branch structure; worst relative error per region (200 points each), against the reference ref_mpmath.dQda:

region	develop grad_reg_inc_gamma	develop grad_reg_lower_inc_gamma	this PR
series, a < 1, z < 8	1.5e-03	1.5e-10	8.5e-15
series, 1 <= a <= 10	7.7e-03	9.6e-03	2.1e-15
series, 10 < a <= 170	2.6e+29	5.8e-01	1.2e-13
asymptotic, a <= z <= 2a	5.1e-01	1.1e+07	1.7e-13
asymptotic, z >= 3a	3.7e-04	2.1e+165	5.4e-12 (4.5e-14)
seam z = 8	1.2e-02	9.9e-12	4.4e-15
integer a	1.9e+15	2.6e+01	3.3e-14
a > 171.6	176 of 200 NaN	1.1e+01	3.8e-13
z >= 3a, a >= 12	9.6e-07	1.9e+127	2.5e-11 (9.4e-14)

In the two rows with a value in parentheses, the error comes from the reference. The value in parentheses is the worst error of this PR at the 200 points of the row against mpmath's numerical derivative mp.diff.

A second sweep compared 49 990 points with a and z both in [1e-3, 1e4], dense on z ≈ a, the z = 8 seam, a = 171.62 and a log z = 709. develop returns NaN or inf at 7 379 of them for grad_reg_inc_gamma; this PR is finite at every one, with median relative errors between 1e-16 and 1e-14 in every region. Outside the band near z = a, the worst error against the reference is 3e-11. At the worst point of each region, this PR is within 1e-13 of mp.diff, so there the 3e-11 is the error of the reference.

One known limit: near z = a the error grows with a, to 2.9e-7 for a < 10 000 and 5.2e-6 at a = z = 1e4, where Eigen's Cephes recurrence loses digits near the mode. Three independent references (numerical differentiation, a quadrature and a series) agree to 16 digits at the worst point, so the error comes from the algorithm. For a < 1 000 the worst error against the reference is 3e-11; at that point this PR is within 6e-14 of mp.diff. In the far upper tail at a ≳ 1e4 the reference is certified to 11 digits only.

Speed

Function level, net ns per call, Xeon E5-2680 v3, performance governor, two rounds (spread under 3 %), Triton job 20356756, on a in [0.5, 60] and z in [0.1, 40], where develop terminates and is finite:

kernel	develop	this PR	speed-up
grad_reg_inc_gamma, double, with tgamma and digamma	805	370	2.2x
grad_reg_inc_gamma, fvar	1 749	592	3.0x
grad_reg_inc_gamma, fvar<fvar>	3 230	1 320	2.4x
grad_reg_lower_inc_gamma, double	1 135	256	4.4x
grad_reg_lower_inc_gamma, fvar	1 766	468	3.8x
grad_reg_lower_inc_gamma, fvar<fvar>	2 536	1 133	2.2x
grad_reg_inc_gamma, var, value + grad	2 690	3 727	0.72x
grad_reg_lower_inc_gamma, var, value + grad	3 921	3 618	1.1x

The var instantiation of the upper root is slower: Eigen's series and continued fraction take more operations than the 9-term asymptotic branch where that branch runs, and with var every operation is a tape node. Reverse-mode code in the library calls the roots with double values; only the fvar<var> path of the cdfs passes var.

Part 2: the derivatives of gamma_q and gamma_p

#2006 shows both problems at (a, z) = (3, 36): gamma_q with fvar<double> returns -6.7e+12 for a derivative of -11.5, and gamma_p with var returns zero derivatives.

The reverse-mode gamma_p in rev/fun/gamma_p.hpp returned without derivatives when z / a > 10, with the comment that gamma_p is flat to machine precision there. Its derivatives there are small but not zero (about 1e-13 at (3, 36)), and log1m(gamma_p(...)), the case of #2006, divides them by 1 - P, so its whole gradient was lost. With the new roots the derivatives are accurate in that region, and this PR removes the short-cut from all three overloads.

stan/math/fwd/fun/gamma_q.hpp on develop held its own copy of the series in grad_reg_inc_gamma, with while (fabs(delta) > 1e-6) and no iteration cap. For z above about 700 the alternating terms (-z)^k / k! overflow before they decay, delta stays inf, and the loop never exits; a call at a = 0.168, z = 4 748 did not return in 15 minutes. The asymptotic branch of the old root calls gamma_q with fvar arguments, so grad_reg_inc_gamma reached the loop too. The existing mix test for gamma_q carried a commented-out case with the note "this still fails forward mode", which is the same defect. This PR enables it again.

This PR replaces the copy with a call to grad_reg_inc_gamma, as fwd/fun/gamma_p.hpp already calls grad_reg_lower_inc_gamma. It also evaluates the derivative with respect to z, -z^(a-1) e^(-z) / Γ(a), in log space; it was NaN at (a, z) = (200, 150) because pow and tgamma overflow.

The copy had been hiding the defect. The generated third-order case in test/prob/chi_square compares the analytic partials against autodiff through log(gamma_q(nu * 0.5, y * 0.5)). Both routes used the same inaccurate series, so the errors canceled. With the accurate root and the stale copy (commit b172d61), the case fails by 2.4e-03; with the copy removed (commit af3205c), it passes.

Part 3: the cdf functions in log space

What failed

gamma_lccdf computed log Q in log space, but the shape derivative as dQ/dalpha / exp(log Q). That is NaN where log Q < -708, for example at (alpha, y) = (10, 1000) and (200, 1600). At (150, 1200), develop returns 2; the true value is 2.084. With only the roots of Part 1 (commit b172d61), it returns 0. @mjcarter95 shows this region in the plot in this thread. The same plot shows that gamma_lccdf inherits defects 1 and 2 of Part 1 for alpha + 1 < y <= 2 alpha. With fvar<fvar<double>> and fvar<fvar<var>>, the second derivative with respect to alpha was exactly 0 at every point: the code seeded a copy of alpha and kept only a plain number as the partial. With fvar<var>, as hessian() uses it, the second derivative was NaN at the tested points alpha = 700 and 1e4.

The other 17 functions had three defects, measured at 18 points (a, z) of the incomplete gamma function against mpmath, with the roots of Part 1 already in place (commit 8d376d5):

1. The value underflows. 11 of the 12 log-scale functions computed log(gamma_p(...)) or log(gamma_q(...)). Where the probability is below 1e-308, the value is -inf and the gradient is NaN. Where it is subnormal, the value has relative error 1.5e-4, and the gradient is NaN or has relative error 0.1 to 1. Example: chi_square_lccdf(2000 | 20) is -inf; the true value is -950.6.
2. The density factor in the y, scale and lambda partials overflows. 14 functions wrote it as exp(-z) pow(z, a - 1) / tgamma(a) (poisson_cdf as exp(-lambda) pow(lambda, n) / tgamma(n + 1)). tgamma overflows above 171.6 and pow once (a - 1) log z > 709, so the partial is NaN or infinite also in the bulk, for example the y gradient of chi_square_cdf, chi_square_lccdf and gamma_cdf at the median for nu = 400.
3. The probability-scale cdfs divide by the probability, so their gradient is NaN where the cdf underflows to 0, although the true gradient is 0.

The change

One internal function, log_gamma_tail(a, z) in log_gamma_q_dgamma.hpp, returns log P or log Q with its derivatives with respect to a and z, all /in log space and generic in the scalar type:

• For z > a + 1, the continued fraction gives log Q, and with the derivative carried through the recurrence also d/da log Q. It stops only when the value and the derivative have converged. log P = log1m_exp(log Q).
• Otherwise gamma_p gives P, and grad_reg_lower_inc_gamma its derivative. Where P is below DBL_MIN / EPSILON, a power series in log space gives log P and its derivative; it stops on a bound of the remaining terms.
• d/dz is the density divided by the probability, exp(log density - log prob).

Each of the 18 functions keeps its argument checks and applies its own chain rule. The probability-scale cdf sums the logs and partials and multiplies by exp(sum) at the end, so it has neither the density overflow nor the 0/0. A derivative that is differentiated again is computed in the autodiff type, so nested forward mode and hessian() get correct second derivatives.

Three behavior changes at the edges of the support:

• gamma_lcdf and chi_square_lcdf at y = inf, and the lccdf of chi_square, inv_gamma, inv_chi_square and scaled_inv_chi_square at y = 0, returned 0 for the whole vector when one element had probability 1, so the other elements were lost. Such an element now contributes 0.
• poisson_lccdf returned 0 for the whole vector when one count was negative (Incorrect behaviour of CCDFs when there is a negative random variable #56). Such an element has probability 1 and now contributes 0.
• The Poisson cdfs treat n = INT_MAX as an infinite count, as neg_binomial, beta_neg_binomial and yule_simon already do. Before, poisson_lccdf returned -inf there only because gamma_p underflowed.

Accuracy

At the 324 test points (the 18 functions at 18 points), no value or gradient is non-finite. The worst relative error is 2.4e-13 for values and 8e-13 for gradients, with one exception: the nu partial of scaled_inv_chi_square near the mode has a relative error up to 2.7e-10. There a and z both grow with nu, and the two terms of the chain rule nearly cancel; the absolute error is of rounding size. For gamma_lccdf in all autodiff modes, first derivatives are within 5.6e-12 and second derivatives within 1.9e-10 (7e-9 near y = alpha + 1 at alpha = 5). On the grid of the plot by @mjcarter95, all 266 cells are finite and the worst relative error of the shape derivative is 3.4e-12.

Speed

Value plus all gradients with var arguments, ns per call, AMD EPYC 9965, g++ 11.5, -O3, medians of 4 alternating runs, Roihu job 2430495. Points with a = 5 and a = 50 and z from 0.2 a to 3 a, mapped to each function's arguments:

function	develop, a = 5	this PR, a = 5	speed-up, a = 5	develop, a = 50	this PR, a = 50	speed-up, a = 50
gamma_cdf	366	178	2.1x	1 358	262	5.2x
gamma_lcdf	934	174	5.4x	3 156	254	12.4x
gamma_lccdf	475	161	3.0x	1 511	235	6.4x
chi_square_cdf	330	172	1.9x	1 256	252	5.0x
chi_square_lccdf	602	162	3.7x	1 638	236	6.9x
inv_gamma_cdf	632	174	3.6x	1 740	252	6.9x
inv_gamma_lcdf	640	167	3.8x	1 743	241	7.2x
scaled_inv_chi_square_lcdf	598	170	3.5x	1 639	247	6.6x
poisson_lcdf	351	95	3.7x	527	147	3.6x

Most of the gain comes from no longer calling tgamma, pow and a root with tgamma and digamma per element.

gamma_lccdf alone, shape derivative only, with alpha from 5 to 1000 (ns per call):

mode	region	develop	this PR	speed-up
var	y > alpha + 1, normal Q	658	105	6.3x
var	Q underflows	774	68	11.4x
var	y <= alpha + 1	13 420	385	34.9x
fvar<double>	y > alpha + 1, normal Q	626	98	6.4x
fvar<var>	y > alpha + 1, normal Q	2 659	1 655	1.6x
fvar<var>	Q underflows	2 804	805	3.5x
fvar<var>	y <= alpha + 1	60 997	7 926	7.7x

Where Q underflows, develop returns NaN or 0, so its time there is no baseline. One of the four runs of this PR was about 35 % slower in three rows of the cdf table; the medians do not depend on it.

Testing

New tests, all with fixed references from mpmath; expect_ad compares against finite differences and cannot see errors of this size:

• test/unit/math/prim/fun/grad_reg_inc_gamma_accuracy_test.cpp: the roots at the points of Part 1, including the series-branch points of defects 2 and 4, references at 80 digits.
• ProbInternalMath.gradRegIncGamma_gradient_matches_hessian: the first derivative of gamma_cdf is the same from reverse mode and from hessian().
• gammaQ_fwd_shape_derivative_large_z and gammaQ_fwd_z_derivative_large_a in test/unit/math/mix/fun/gamma_q_test.cpp: the forward-mode gamma_q derivatives. A test that hangs is not usable in CI, so the first one uses (a, z) = (1.5, 300), where the old series terminates but returns 1.2e+112 for a derivative of 5.7e-129. gammaQ_fwd_log_upper_tail_composite is the fvar<double> example of gamma_p and gamma_q gradient problems #2006, log(gamma_q(alpha, beta * y)) at y = 12, alpha = 3, beta = 3 with all tangents 1.
• test/unit/math/rev/prob/incomplete_gamma_cdfs_reference.hpp with incomplete_gamma_cdfs_reference_test.cpp: 334 rows. They hold the value and gradient of the 18 functions at 18 points, and at the 10 points of the Part 1 table, at 60 digits by two routes (the incomplete gamma function with numerical differentiation, and a quadrature of the scaled density with the chain rule), which agree to 1e-25. Value to 1e-12, gradient to 1e-9 relative.
• test/unit/math/rev/prob/gamma_lccdf_tail_reference.hpp: 13 rows for gamma_lccdf with second derivatives, used by a var case in rev/prob/gamma_lccdf_test.cpp and by fvar<double>, fvar<var>, fvar<fvar<double>> and fvar<fvar<var>> cases in mix/prob/gamma_lccdf_test.cpp.
• mix/prob/incomplete_gamma_cdfs_test.cpp and mix/prob/inv_incomplete_gamma_cdfs_test.cpp: expect_ad for all 18 functions at two points each, one on each branch.
• gammaP_rev_far_upper_tail in test/unit/math/mix/fun/gamma_p_test.cpp: the reverse-mode derivatives of gamma_p(3, 36) in all three overloads, and the gradient of log1m(gamma_p(alpha, beta * y)) from gamma_p and gamma_q gradient problems #2006.
• lccdf_negative_count_contributes_zero in test/unit/math/prim/prob/poisson_ccdf_log_test.cpp: the example of Incorrect behaviour of CCDFs when there is a negative random variable #56.

Changed tests:

• Four constants in test/unit/math/mix/fun/grad_reg_inc_gamma_test.cpp expected 0.38984156 for d/da Q(0.5, 1.0); the true value is 0.38983726432851057. The old constants combined the rounded g = 1.77245 that the test passes in and the 1e-6 series truncation. The two second- and third-derivative cases of that file were not changed and still pass.
• lccdf_log1m_exp_lcdf_rounds_to_inf (prim and rev) asserted that gamma_lcdf(20000 | 800, 0.1) rounds to exactly 0. gamma_lcdf now returns its true value, -1.48e-205, so the case checks this value and that log1m_exp(gamma_lcdf) equals gamma_lccdf there, and is renamed to lccdf_far_tail_matches_log1m_exp_lcdf.

The new tests fail without the changes. I installed the develop versions of the 25 changed headers that exist on develop, all at once, and kept the new tests. Then I restored the headers of this PR; their checksums equal those before. 23 of the 24 new test cases fail with the old headers. The exception is mathMixScalFun_poisson_cdfs, an expect_ad case that passes on develop too. The reference test fails with 588 expectations, in all 18 functions. The four changed constants of grad_reg_inc_gamma_test.cpp fail too.

The existing unit suites of the roots, gamma_p, gamma_q and all 18 functions pass, and so do the 90 generated test/prob cdf programs of the six families. The existing grad_reg_lower_inc_gamma_test checks 9 600 values from Mathematica with a 1e-10 absolute tolerance; Eigen satisfies all of them, but that grid (a up to 19.75, z up to 29.75) contains no point of defect 3. I tested the merge of this PR with develop at 510d263078 (Triton, g++ 14.3). I also tested the earlier commits with g++ 11.5.

Reproduction

The scripts behind every number in this text are in https://github.com/avehtari/stan-pr-checks/tree/27838b5236/pr3408. They cover the mpmath generators of the test references (two routes each), the two sweeps of the roots, the survey of the 18 cdfs, the grid of the plot by @mjcarter95, the forward-mode hang, the swap test and the speed benchmarks. The README in that directory gives one command per claim, the expected output and the commit that I measured it on. On 2026-10-28 I ran the accuracy measurements again on develop (5252d51d47) and on the merged head. I did not run the speed benchmarks again.

Dependency on Eigen

The two new headers specialize lgamma_impl, digamma_impl and cephes_helper inside Eigen::internal. Stan Math already specializes inside that namespace, in stan/math/rev/core/Eigen_NumTraits.hpp. An Eigen upgrade that renames those three helpers breaks the build; that failure is immediate, and the alternative was to copy about 300 lines of the Cephes algorithm into Stan Math. The headers belong to the fun layer, not core, because they need digamma and lgamma; in core they break the include order of every mix target.

Not included

OpenCL is unaffected: no OpenCL code calls these functions, and the six families have OpenCL implementations only for their lpdf or lpmf.

The incomplete beta gradient roots have their own defects; they are in #3409.

#2769 asks for log-scale incomplete gamma functions, because gamma_lcdf lost resolution at the boundaries. This PR computes log P and log Q in log space inside the library (internal::log_gamma_tail), which meets that motivation for the cdfs, but it adds no public function such as log_gamma_p or log_gamma_q.

Other discrete cdfs have the defect of #56: one element with probability 1 decides the result for the whole vector. bernoulli_lccdf, binomial_lccdf, beta_binomial_lccdf, neg_binomial_lccdf (and so neg_binomial_2_lccdf), beta_neg_binomial_lccdf and yule_simon_lccdf return 0 when one count is below the support; neg_binomial_lcdf, beta_neg_binomial_lcdf and yule_simon_lcdf return 0, and neg_binomial_cdf, neg_binomial_2_cdf, beta_neg_binomial_cdf and yule_simon_cdf return 1, when one count is INT_MAX. Most of these files are changed by other open PRs, so I'll make the PR for this after those other PRs have been merged.

The missing poisson_log_cdf, poisson_log_lcdf and poisson_log_lccdf (#90, #2657) can be written as short wrappers of the same log-space function, log Q(n + 1, exp(eta)); #3476 adds them.

g, dig, precision and max_steps of the roots are now unused. Removing them changes the remaining call sites (gamma_q in rev and fwd, and log_gamma_q_dgamma) and the public signature, so it is left for a separate PR. log_gamma_q_dgamma itself is no longer called by the library and is kept for compatibility.

Release notes

Use Eigen's igamma_der_a for the incomplete gamma gradient roots, and compute the cdf, lcdf and lccdf of chi_square, gamma, inv_chi_square, inv_gamma, scaled_inv_chi_square and poisson in log space, so that their values and gradients stay finite in the tails and for large shapes.

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

[avehtari]

I'll be away for two weeks, but did create this PR to avoid overlap in work. Similar PR should be made for incomplete beta gradient roots (mentioned in PR text)

[stan-buildbot]

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan	0.23	0.23	1.0	0.05% faster
stat_comp_benchmarks/benchmarks/gp_regr/gen_gp_data.stan	0.06	0.06	1.0	-0.18% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	6.34	6.35	1.0	-0.15% slower
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	1.02	2.22% faster
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	8.44	8.57	0.99	-1.47% slower
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	4.49	4.49	1.0	0.07% faster
stat_comp_benchmarks/benchmarks/sir/sir.stan	165.55	170.86	0.97	-3.21% slower
stat_comp_benchmarks/benchmarks/garch/garch.stan	0.89	0.89	1.0	-0.01% slower
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.71	0.71	1.0	-0.48% slower
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	42.71	42.61	1.0	0.24% faster
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.6	0.58	1.04	3.51% faster
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	20.8	20.9	1.0	-0.49% slower
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.11	0.11	1.0	-0.19% slower
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.21	3.19	1.01	0.55% faster
performance.compilation	384.38	382.76	1.0	0.42% faster
Mean result: 1.0007941730987455

Jenkins Console Log
Jenkins Build Stages
Commit hash: c865b6f141c9d526e21e2f2e20a3edf2aa0e1fe0

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:                                 1494.753
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

[mjcarter95]

@avehtari A colleague of mine noticed instability in the gradient of Stan's gamma_lccdf function which was occurring when a right-censored datum was larger than the mean of the Gamma distribution. This PR seems to have addressed some of the instability, however, the gradient of the gamma_lccdf method still appears to be unstable when moving deeper into the tail. Are there any plans to investigate / patch this?

[avehtari]

Are there any plans to investigate / patch this?

@mjcarter95 There was no plan before your comment! Thanks for testing and reporting, I'll investigate a possible fix

[mjcarter95]

Are there any plans to investigate / patch this?

@mjcarter95 There was no plan before your comment! Thanks for testing and reporting, I'll investigate a possible fix

Thank you, Aki. Claude helped me to identify that the issue arises because Stan computes the shape gradient as dQ/dα divided by exp(log Q), and once log Q falls below about −745, both terms underflow to zero and the result is 0/0 = NaN.

Are there guidelines for how to submit an AI-assisted PR? If so, I'd be happy to supervise Claude in coming up with a fix and submitting a PR that piggybacks off this one.

[avehtari]

I'm already fixing this and what you mentioned is only part of the problem. It's just that running the accuracy, speed and unit tests takes some time and also checking which other functions can be improved, so please wait a few hours.

[avehtari]

@mjcarter95 Just an intermediate report that I'm making fixes for 17 other functions in addition to gamma_lccdf. Better to fix all those functions at the same time which were already affected by the roots part of this PR and which now gain also similar fixes as gamma_lccdf after you asked stability for wider input range, and which are not handled by other open PRs

[mjcarter95]

@mjcarter95 Just an intermediate report that I'm making fixes for 17 other functions in addition to gamma_lccdf. Better to fix all those functions at the same time which were already affected by the roots part of this PR and which now gain also similar fixes as gamma_lccdf after you asked stability for wider input range, and which are not handled by other open PRs

Brilliant, thank you for following up so quickly!

[avehtari]

@mjcarter95 The cdf, lcdf and lccdf of chi_square, gamma, inv_chi_square, inv_gamma, scaled_inv_chi_square and poisson had the same tail failures, and most of them also had NaN gradients at large shape. They now use the same log-space computation, and the updated PR text has now the details and the timings. Thanks again for pointing out the issue!

Here's my version of the figure (with probably more accurate reference as when using the same color scale as you used, it is very boring)

[mjcarter95]

Re-ran the tests my end, all looks good! I am comparing against mpmath.

Thank you for fixing this @avehtari

[avehtari]

I am comparing against mpmath.

I guess you are using lower accuracy in mpmath, as you report worst relative error 4.1e-6 and I have worst relative error 3.4e-12 (I used 40 to 80 digits accuracy)

[mjcarter95]

I am comparing against mpmath.

I guess you are using lower accuracy in mpmath, as you report worst relative error 4.1e-6 and I have worst relative error 3.4e-12 (I used 40 to 80 digits accuracy)

I used precision 60 for all cells.

[andrjohns]

@avehtari can you split this PR into smaller batches? There's a bit too much to review properly, and doesn't all need to land in the one PR

[avehtari]

an you split this PR into smaller batches? There's a bit too much to review properly, and doesn't all need to land in the one PR

OK. I have worried that I have too many PRs, but I can do more.

[avehtari]

@andrjohns how many PRs would you prefer? There are 5 logical parts and thus it would be possible to split this in 5 PRs

part	content	base	files	lines
G1 (#3408 reused)	gradient roots by Eigen igamma_der_a: grad_reg_inc_gamma, grad_reg_lower_inc_gamma, fwd/rev Eigen_SpecialFunctions, root tests	develop	8	+293 -239
G2	fwd gamma_q (no duplicated series, log-space z derivative), rev gamma_p without the z / a > 10 short cut (Closes #2006), their tests	G1	4	+107 -60
G3	internal::log_gamma_tail in log_gamma_q_dgamma.hpp, gamma_cdf, gamma_lcdf, gamma_lccdf (mjcarter95's report), tail reference test, gamma rows of the reference table	G1	~10	~+800 -350
G4	the same pattern for inv_gamma, chi_square, inv_chi_square, scaled_inv_chi_square (12 cdfs, generated by gen_cdfs.py), their reference rows and mix tests	G3	~13	~+700 -700
G5	poisson_cdf, _lcdf, _lccdf on log_gamma_tail; negative counts (part of #56); Poisson reference rows	G3	~5	~+200 -50

G2 and G3 are independent of each other; G4 and G5 are independent of each other.
In addition #3476 follows G5.

Instead of 5 PRs it would possible to combine G1 and G2 in one PR and/or G3, G4 and G5 to one PR

Let me know what you prefer, and I'll do the split

[stevieb-bot]

Fact check, part 1 of 2.

Action: Not Ready for Human Review

Commit checked: 22d781f290cb

Summary

Thank you for this PR and for the detailed write-up of the defects on develop.

The PR replaces the incomplete gamma gradient roots with Eigen's igamma_der_a, removes the copied series in the forward-mode gamma_q and the z / a > 10 short cut in the reverse-mode gamma_p, and moves the cdf, lcdf and lccdf of six families onto one log-space function.
It reports accuracy and speed tables for Parts 1 and 3.

I read the source and the tests at 22d781f and ran the PR's new and changed tests in a container without network access.
I ran each test file twice, once at the PR head and once at the merge base 5252d51 with the PR's test files copied in.
I also ran the existing unit tests of the roots and of the 18 functions, and the generated test/prob cdf programs of the six families, at the head.

The output of each run is in this gist.
Its README.md gives the command and commit behind each file.

• Claims and Evidence: 8 Confirmed, 19 Partially Confirmed, 6 Not Verified.
• Description vs Code: 15 Confirmed, 4 Partially Confirmed, 1 Not Verified.

Claims and Evidence

Each section below lists facts from the PR summary and discussion that an agent was able to confirm, partially confirm, or not confirm at all, or that a run of the code in the PR contradicted.
Please review the "Partially Confirmed", "Not Verified", and "Contradicted" facts for facts missing evidence found by this agent.

Confirmed Facts

FACT: The reverse-mode gamma_p derivatives at (3, 36) are about 1e-13, not zero

Quote:

Its derivatives there are small but not zero (about 1e-13 at (3, 36)), and log1m(gamma_p(...)), the case of #2006, divides them by 1 - P, so its whole gradient was lost.

(PR description, Part 2)

STATUS: Confirmed

gammaP_rev_far_upper_tail checks the derivatives of gamma_p(3, 36) in the vv, vd and dv overloads against -4.27e-13 in a and 1.50e-13 in z.
On develop each of the four adjoints is 0.
At the head each is within 1e-14 absolute of the stored value, because expect_near_rel with a tolerance of 1e-10 uses an absolute floor of 1e-14 at this size (relative_tolerance.hpp#L29-L30).
The log1m case in the same test uses the reference gradients from #2006, and develop returns 0 for all three of them.

Evidence:

The PR adds a test with mpmath references for the three overloads and the case of #2006.
Code: gamma_p_test.cpp.

What I ran: python3 runTests.py -j8 --run-all on mix/fun/gamma_q_test.cpp and mix/fun/gamma_p_test.cpp, at the head and at the merge base with the PR's tests copied in.
Output in the gist: gamma_pq-head.out, gamma_pq-base.out.

---

FACT: The forward-mode gamma_q derivative in z was NaN at (200, 150)

Quote:

It also evaluates the derivative with respect to z, -z^(a-1) e^(-z) / Γ(a), in log space; it was NaN at (a, z) = (200, 150) because pow and tgamma overflow.

(PR description, Part 2)

STATUS: Confirmed

gammaQ_fwd_z_derivative_large_a takes the fvar<double> derivative of gamma_q(200, 150) in z.
On develop it is not finite, and at the head it is within 1e-12 relative of -2.00507038377120400e-05, the closed form -z^(a-1) e^(-z) / Gamma(a).

Evidence:

The PR adds a test at the named point, with the closed-form reference.
Code: gamma_q_test.cpp.

What I read: Develop's fwd/fun/gamma_q.hpp#L39 computes der2 = -exp(-x2.val_) * pow(x2.val_, x1.val_ - 1.0) / g with g = tgamma(x1.val_).

What I ran: python3 runTests.py -j8 --run-all on mix/fun/gamma_q_test.cpp and mix/fun/gamma_p_test.cpp, at the head and at the merge base with the PR's tests copied in.
Output in the gist: gamma_pq-head.out, gamma_pq-base.out.

---

FACT: The second derivative of gamma_lccdf in alpha was 0 with nested fvar, and NaN with fvar

Quote:

With fvar<fvar<double>> and fvar<fvar<var>>, the second derivative with respect to alpha was exactly 0 at every point: the code seeded a copy of alpha and kept only a plain number as the partial.
With fvar<var>, as hessian() uses it, the second derivative was NaN at the tested points alpha = 700 and 1e4.

(PR description, Part 3, What failed)

STATUS: Confirmed

mathMixScalFun_gamma_lccdf_upper_tail_reference runs the 13 tail rows with fvar<double>, fvar<var>, fvar<fvar<double>> and fvar<fvar<var>>.
With develop's headers, the second derivative in alpha from fvar<fvar<double>> and fvar<fvar<var>> is 0 at every row (relative error 1), and the fvar<var> second derivative is -nan at alpha = 700 and alpha = 1e4, among others.
The expect_ad test in the same file also fails on develop, because the fvar<fvar<double>> Hessian is 0 against a finite-difference value of -1.34.
At the head both tests pass.

Evidence:

The PR adds a mix test that runs the 13 tail rows in four autodiff types.
Code: mix/prob/gamma_lccdf_test.cpp.

What I read: Develop's gamma_lccdf.hpp#L59-L64 seeds alpha_unit.d_ = 1 on a copy of alpha and keeps log_Q_fvar.d_ as the partial.

What I ran: python3 runTests.py -j8 --run-all on the prim, rev and mix gamma_lccdf_test.cpp files, at the head and at the merge base with the PR's tests copied in.
Output in the gist: gamma_lccdf-head.out, gamma_lccdf-base.out.

---

FACT: Repairs the gamma_lccdf shape-gradient failures reported by @mjcarter95

Quote:

• Repairs the gamma_lccdf shape-gradient failures that @mjcarter95 reported in this thread (Part 3).

(PR description, Issues)

STATUS: Confirmed

The reported failure is the shape gradient of gamma_lccdf deep in the upper tail.
The 13 rows of the tail reference test include 8 points with log Q < log(DBL_MIN). On develop the var shape derivative is -nan at 7 of the 13 rows and off by up to 4 % at others, and at the head all 13 rows are within 1e-11 relative of the reference in var, fvar<double>, fvar<var>, fvar<fvar<double>> and fvar<fvar<var>>.
The 266-cell grid of the plot is a separate fact check.

Evidence:

The PR adds a tail reference test of 13 rows in five autodiff types.
Code: gamma_lccdf_tail_reference.hpp, rev/prob/gamma_lccdf_test.cpp, mix/prob/gamma_lccdf_test.cpp.

What I ran: python3 runTests.py -j8 --run-all on the prim, rev and mix gamma_lccdf_test.cpp files, at the head and at the merge base with the PR's tests copied in.
Output in the gist: gamma_lccdf-head.out, gamma_lccdf-base.out.

---

FACT: The probability-scale cdfs return NaN gradients where the cdf underflows to 0

Quote:

3. The probability-scale cdfs divide by the probability, so their gradient is NaN where the cdf underflows to 0, although the true gradient is 0.

(PR description, Part 3, What failed)

STATUS: Confirmed

In the reference test on develop, chi_square_cdf has -nan in both gradients at (200, 2000), (60, 1000) and (0.002, 200), where the stored cdf is 1.03e-611, 2.97e-409 and 1.07e-458.
Develop's chi_square_cdf.hpp#L100-L110 divides both partials by Pn.
In the first two rows a is 1000 and 500, so tgamma(alpha_dbl) in the y partial also overflows.
The row (0.002, 200) has a = 100, so its NaN does not involve that overflow.
At the head these rows pass with an absolute error below 1e-300, which is how the test compares references below 1e-300.

Evidence:

The reference test has rows where the cdf is below 1e-300.
Code: incomplete_gamma_cdfs_reference_test.cpp.

What I read: Develop's chi_square_cdf.hpp L100-L110 divides the y and nu partials by Pn.

What I ran: python3 runTests.py -j8 --run-all on rev/prob/incomplete_gamma_cdfs_reference_test.cpp, prim/prob/poisson_ccdf_log_test.cpp and the two mix cdf tests, at the head and at the merge base with the PR's tests copied in.
Output in the gist: cdfs-head.out, cdfs-base.out.

---

FACT: 199 of 324 test points are non-finite on develop, and none with the PR

Quote:

At the 324 test points (the 18 functions at 18 points), no value or gradient is non-finite.

(PR description, Part 3, Accuracy)

STATUS: Confirmed

The reference test has 324 rows, 18 functions at 18 points each.
In the base run, 199 rows have at least one value or gradient reported as "not finite", which matches the Summary row "18 cdfs, test points with a non-finite value or gradient, of 324" (199 against 0).
At the head the test passes, and it asserts that each value and gradient is finite.

Evidence:

The reference test asserts that each value and gradient is finite.
Code: incomplete_gamma_cdfs_reference_test.cpp.

What I ran: python3 runTests.py -j8 --run-all on rev/prob/incomplete_gamma_cdfs_reference_test.cpp, prim/prob/poisson_ccdf_log_test.cpp and the two mix cdf tests, at the head and at the merge base with the PR's tests copied in.
I counted the distinct (function, point) pairs with a "not finite" line in cdfs-base.out.
Output in the gist: cdfs-head.out, cdfs-base.out.

---

FACT: Partly addresses #56 for poisson_lccdf

Quote:

• Partly addresses Incorrect behaviour of CCDFs when there is a negative random variable #56: poisson_lccdf no longer drops the other elements of a vector when one count is negative.

(PR description, Issues)

STATUS: Confirmed

lccdf_negative_count_contributes_zero checks that poisson_lccdf(-1, 5.0) and poisson_lccdf(-1, 0.0) are 0, and that the vector calls with counts {-1, 1} equal poisson_lccdf(1, 5.0).
On develop both vector cases return 0 against -0.0413, and at the head the test passes.

Evidence:

The PR adds a test with the example of #56.
Code: poisson_ccdf_log_test.cpp.

What I ran: python3 runTests.py -j8 --run-all on rev/prob/incomplete_gamma_cdfs_reference_test.cpp, prim/prob/poisson_ccdf_log_test.cpp and the two mix cdf tests, at the head and at the merge base with the PR's tests copied in.
Output in the gist: cdfs-head.out, cdfs-base.out.

---

FACT: The existing unit suites and the 90 generated test/prob cdf programs pass

Quote:

The existing unit suites of the roots, gamma_p, gamma_q and all 18 functions pass, and so do the 90 generated test/prob cdf programs of the six families.

(PR description, Testing)

STATUS: Confirmed

At the head, 24 existing unit test files pass, with 67 tests.
The files are the prim tests of gamma_p, gamma_q and the two roots, the prim tests of the six families including their cdf_log and ccdf_log files, rev/prob/gamma_lcdf_test.cpp, rev/prob/inv_gamma_test.cpp and mix/prob/gamma_cdf_log_test.cpp.
The PR's changed tests pass too, as the runs above show.
make generate-tests produces 90 cdf programs for the six families, and all 90 pass at the head, with 59 940 tests.

Evidence:

The suites and the generator are in the repository, and the PR names them.
Code: test/prob/chi_square, grad_reg_lower_inc_gamma_test.cpp.

What I ran: python3 runTests.py -j8 --run-all on the 24 existing files.
For test/prob, make -j8 generate-tests -s, then python3 runTests.py -j8 --run-all on the generated *cdf*_generated_*_test.cpp files of the six families.
review.py run cut the middle of the long test/prob output, so I ran the built programs again with runTests.py -t --run-all and kept only the command and summary lines.
Output in the gist: existing-head.out, prob-cdfs-head.out, prob-cdfs-summary-head.out.

---

Partially Confirmed Facts

FACT: The error near z = a grows with a, to 5.2e-6 at a = z = 1e4

Quote:

One known limit: near z = a the error grows with a, to 3e-11 for a < 1 000, 2.9e-7 for a < 10 000, and 5.2e-6 at a = z = 1e4, where Eigen's Cephes recurrence loses digits near the mode.

(PR description, Part 1, Accuracy)

STATUS: Partially Confirmed

The accuracy test checks a = z at 20, 100, 171, 180, 300 and 500 with a relative tolerance of 1e-10, and it passes at the head.
That shows the a < 1 000 part only as an error below 1e-10 at those six points.
The PR does not include code that measures the error near z = a above a = 500, so the PR does not include code for 2.9e-7, 5.2e-6, or the three references.

Evidence:

The accuracy test covers a = z up to 500 with a bound of 1e-10.
No script or test covers larger a.
Code: grad_reg_inc_gamma_accuracy_test.cpp.

What I ran: python3 runTests.py -j8 --run-all on prim/fun/grad_reg_inc_gamma_accuracy_test.cpp and mix/fun/grad_reg_inc_gamma_test.cpp, at the head and at the merge base with the PR's tests copied in.
Output in the gist: roots-head.out, roots-base.out.

Request: Could you add the script that measures the error near z = a for a up to 1e4, with the three references?

---

FACT: grad_reg_inc_gamma on develop is wrong near z = a, by 4.4e-02 to 7.2e-01

Quote:

The relative error at z = a is 4.4e-02 at a = 20, 5.8e-01 at a = 100 and 7.2e-01 at a = 171.

(PR description, Part 1, The defects on develop)

STATUS: Partially Confirmed

grad_reg_inc_gamma_diagonal_z_equals_a checks a = z at 20, 100 and 171 against stored mpmath values with a relative tolerance of 1e-10.
With develop's headers it fails with relative errors 0.0443, 0.577 and 0.722, the same as the PR's 4.4e-02, 5.8e-01 and 7.2e-01.
At the head it passes.
grad_reg_inc_gamma_agrees_with_lower, which compares the two roots with each other and needs no reference, also fails on develop (relative differences from 1.6e-4 at a = z = 5 to 2.59 at a = z = 171) and passes at the head.
I could not check the reference values, because the PR does not include the mpmath code that produced them.

Evidence:

The PR does not include the mpmath code behind the expected values of the test.
Code: grad_reg_inc_gamma_accuracy_test.cpp.

What I ran: python3 runTests.py -j8 --run-all on the accuracy test and mix/fun/grad_reg_inc_gamma_test.cpp, at the head and at the merge base with the PR's tests copied in.
The head passes all 14 tests.
The base fails 5 of 5 in the accuracy test and 5 of 9 in the mix test.
With --run-all, runTests.py exits 0 even when a test fails, so the base output shows exit status 0 above the failures.
Output in the gist: roots-head.out, roots-base.out.

Request: Could you add the mpmath script that produced the reference values in grad_reg_inc_gamma_accuracy_test.cpp, with the two routes the comment mentions?

---

FACT: grad_reg_inc_gamma on develop returns NaN for a above about 171.6

Quote:

2. grad_reg_inc_gamma returns NaN for a above about 171.6.
   Every call site passes g = tgamma(a), which is inf there, so the correction term exp(-z + (a-1) log z) * S / g is inf / inf.
   The series branch fails the same way through exp(a * log z) when a log z > 709, for example at (a, z) = (150, 150).

(PR description, Part 1, The defects on develop)

STATUS: Partially Confirmed

grad_reg_inc_gamma_above_tgamma_overflow checks a = z at 180, 300 and 500.
With develop's headers each result is -nan, and at the head each is finite and within 1e-10 relative of the stored value.
gradRegIncGamma_gradient_matches_hessian also fails on develop because the reverse-mode gradient of gamma_cdf is not finite at alpha = 180.
The third sentence names (a, z) = (150, 150) as a failure of the series branch.
On develop, though, grad_reg_inc_gamma.hpp#L65 sends every point with z >= a && z >= 8 to the asymptotic branch.
So (150, 150) does not reach the series branch, where exp(a * l) * S / g is at L106.
The new tests do not include that point.

Evidence:

The test covers the tgamma overflow above 171.6.
The example point of the series branch is in the asymptotic branch on develop.
Code: grad_reg_inc_gamma_accuracy_test.cpp, grad_reg_inc_gamma_test.cpp.

What I read: Develop's grad_reg_inc_gamma.hpp L65 sends (150, 150) to the asymptotic branch (L66-L89).
The series branch, with exp(a * l) * S / g, is at L90-L113.

What I ran: python3 runTests.py -j8 --run-all on prim/fun/grad_reg_inc_gamma_accuracy_test.cpp and mix/fun/grad_reg_inc_gamma_test.cpp, at the head and at the merge base with the PR's tests copied in.
The base output has a=z=180: result is not finite, got -nan, the same for 300 and 500, and reverse mode, alpha = 180 for the hessian test.
Output in the gist: roots-head.out, roots-base.out.

Request: At (150, 150) develop's L65 takes the asymptotic branch.
Could you correct the example point, or name the branch it fails in?

---

FACT: grad_reg_lower_inc_gamma on develop returns rounding noise for z well above a

Quote:

It returns about 1e-14 of rounding noise with arbitrary sign: at (50, 150) the reference is -8.29e-22 and develop returns -4.70e-15; at (50, 500) the reference is -5.34e-148 and develop returns +1.47e-14.

(PR description, Part 1, The defects on develop)

STATUS: Partially Confirmed

grad_reg_lower_inc_gamma_far_upper_tail checks (50, 150), (50, 250), (50, 500) and (170, 510).
With develop's headers the results are -4.697913701728197e-15, 1.3362552414382279e-14, 1.4659036750884937e-14 and 2.0769270875699226e-14, which match the PR's -4.70e-15 at (50, 150) and +1.47e-14 at (50, 500).
The stored references are -8.29e-22 and -5.34e-148 at those two points, and the head agrees with them to 1e-10.
The reference values come from mpmath code that is not in the PR.

Evidence:

The PR does not include the mpmath code behind the expected values of the test.
Code: grad_reg_inc_gamma_accuracy_test.cpp.

What I ran: python3 runTests.py -j8 --run-all on prim/fun/grad_reg_inc_gamma_accuracy_test.cpp and mix/fun/grad_reg_inc_gamma_test.cpp, at the head and at the merge base with the PR's tests copied in.
Output in the gist: roots-head.out, roots-base.out.

Request: Could you add the mpmath script behind the reference values in grad_reg_inc_gamma_accuracy_test.cpp?

---

FACT: The 1e-6 precision of develop's series branch gives errors of 9.2e-08 at a = 0.5 and 9.0e-05 at a = 8

Quote:

4. The default precision of 1e-6 sets the error of the plain series branch: 9.2e-08 at a = 0.5 and 9.0e-05 at a = 8.

(PR description, Part 1, The defects on develop)

STATUS: Partially Confirmed

On develop the series branch stops when a term falls below precision, which defaults to 1e-6, so its error depends on that setting.
The tests do not check a = 0.5 or a = 8 with an exact g.
grad_reg_inc_gamma_unchanged_points fails on develop with relative error 3.98e-08 at (1.1, 0.2), which is in the series branch, and 3.17e-06 at (9, 10), which is in develop's asymptotic branch.
The changed constant in gradRegIncGamma_typical at (0.5, 1.0) shows develop's 0.38984156 against the new 0.38983726, but that test passes the rounded g = 1.77245, so it does not separate the series truncation from the rounding.

Evidence:

The stopping rule follows from develop's code.
The PR does not include code for the two stated numbers.
Code: grad_reg_inc_gamma_accuracy_test.cpp, grad_reg_inc_gamma_test.cpp.

What I read: Develop's grad_reg_inc_gamma.hpp#L53 sets precision = 1e-6, and L105 stops the series when log_delta <= log(precision).
L65 puts (1.1, 0.2) in the series branch and (9, 10) in the asymptotic branch.

What I ran: python3 runTests.py -j8 --run-all on prim/fun/grad_reg_inc_gamma_accuracy_test.cpp and mix/fun/grad_reg_inc_gamma_test.cpp, at the head and at the merge base with the PR's tests copied in.
Output in the gist: roots-head.out, roots-base.out.

Request: Could you add the points a = 0.5 and a = 8 in the series branch, with an exact g, to the accuracy test, or the script that measured these two errors?

---

FACT: The root defects reach the distribution gradients, as in the five-row table

Quote:

The gradients are NaN for chi_square_lccdf and chi_square_lcdf at nu = 360, gamma_cdf and inv_gamma_lccdf at alpha = 180, and scaled_inv_chi_square_lcdf at nu = 360.
With this PR, every row of the table agrees with the reference to all 12 printed digits.

(PR description, Part 1, The defects on develop)

STATUS: Partially Confirmed

The new tests check the roots at the points of the table's rows.
They do not check the five distribution functions at those points.
On develop each of the five functions takes its shape partial from grad_reg_inc_gamma, so an error in the root is an error in the shape partial.
For gamma_cdf at alpha = 100 the root is grad_reg_inc_gamma(100, 100), and develop's value in the base run is 0.016873399816234952, against the stored 0.03992749785778605.
These are the table's -0.0168733998162 and -0.0399274978578 up to sign.
The NaN at alpha = 180 shows in gradRegIncGamma_gradient_matches_hessian, where develop's reverse-mode gradient of gamma_cdf is not finite.
The reference values of the five rows, and the statement that the PR agrees with all five rows to 12 digits, are not covered by a test.

Evidence:

The roots under the table are tested.
The five distribution rows are not.
Code: grad_reg_inc_gamma_accuracy_test.cpp, grad_reg_inc_gamma_test.cpp.

What I read: On develop the shape partial calls grad_reg_inc_gamma in gamma_cdf.hpp#L116, chi_square_lccdf.hpp#L110, chi_square_lcdf.hpp#L110, inv_gamma_lccdf.hpp#L103 and scaled_inv_chi_square_lcdf.hpp#L107.

What I ran: python3 runTests.py -j8 --run-all on prim/fun/grad_reg_inc_gamma_accuracy_test.cpp and mix/fun/grad_reg_inc_gamma_test.cpp, at the head and at the merge base with the PR's tests copied in.
Output in the gist: roots-head.out, roots-base.out.

Request: Could you add the five rows of this table as a test of the distribution functions with var, with the references?

---

FACT: Reverse mode, forward mode and hessian() now agree on the first derivative

Quote:

One algorithm now serves every scalar type and every autodiff order, so reverse mode, forward mode, fvar<fvar<double>>, fvar<var> and hessian() agree.

(PR description, Part 1, The change)

STATUS: Partially Confirmed

gradRegIncGamma_gradient_matches_hessian compares the reverse-mode gradient of gamma_cdf in alpha with the gradient from hessian() at alpha = y from 5 to 300, with a relative tolerance of 1e-12.
It passes at the head.
On develop it fails at alpha = 180, where the reverse-mode gradient is NaN, and the two agree up to alpha = 171.
The tail reference test of gamma_lccdf checks fvar<double>, fvar<var>, fvar<fvar<double>> and fvar<fvar<var>> against one reference and passes at the head, so those modes agree there too.
The Summary row "roots, autodiff-order spread" (2.2e-14 against inf in 7 of 10 regions) and the 5.8e-01 difference of the rejected double-only version are not backed by code in the PR.

Evidence:

Two tests cover gamma_cdf and gamma_lccdf.
The PR does not include code for the spread over regions or the rejected version.
Code: grad_reg_inc_gamma_test.cpp, mix/prob/gamma_lccdf_test.cpp.

What I ran: python3 runTests.py -j8 --run-all on prim/fun/grad_reg_inc_gamma_accuracy_test.cpp and mix/fun/grad_reg_inc_gamma_test.cpp, at the head and at the merge base with the PR's tests copied in.
Also python3 runTests.py -j8 --run-all on the prim, rev and mix gamma_lccdf_test.cpp files, at the head and at the merge base with the PR's tests copied in.
Output in the gist: roots-head.out, roots-base.out, gamma_lccdf-head.out, gamma_lccdf-base.out.

Request: Could you add the script behind the autodiff-order spread row, covering the ten regions?

---

FACT: Closes #2006

Quote:

• Closes gamma_p and gamma_q gradient problems #2006: gamma_q with fvar and gamma_p with var return wrong derivatives (Part 2).

(PR description, Issues)

STATUS: Partially Confirmed

#2006 has two examples at y = 12, alpha = 3, beta = 3.
The first is the reverse-mode gradient of log1m(gamma_p(alpha, beta * y)).
gammaP_rev_far_upper_tail covers it with the reference gradients from the issue (-2.83796, 2.68924, -11.3518).
On develop all three adjoints are 0, and at the head they match within the test's 1e-2 relative tolerance.
The second example is the fvar<double> derivative of log(gamma_q(alphaf, betaf * yf)) with all three tangents set to 1, which the issue gives as -6.70499e+12 against -11.5005.
The tests in the PR do not evaluate that expression.
The new forward-mode tests use (1.5, 300), (3, 650) and (200, 150).

Evidence:

The gamma_p example of #2006 has a test.
The gamma_q example does not.
Code: gamma_p_test.cpp, gamma_q_test.cpp.

What I ran: python3 runTests.py -j8 --run-all on mix/fun/gamma_q_test.cpp and mix/fun/gamma_p_test.cpp, at the head and at the merge base with the PR's tests copied in.
The head passes all 8 tests.
The base fails gammaQ_fwd_shape_derivative_large_z, gammaQ_fwd_z_derivative_large_a and gammaP_rev_far_upper_tail.
Output in the gist: gamma_pq-head.out, gamma_pq-base.out.

Request: Could you add the gamma_q example of #2006, the fvar<double> derivative of log(gamma_q(alpha, beta * y)) at (12, 3, 3) with all tangents 1, to gamma_q_test.cpp?

---

FACT: The forward-mode gamma_q on develop does not terminate for z above about 700

Quote:

For z above about 700 the alternating terms (-z)^k / k! overflow before they decay, delta stays inf, and the loop never exits; a call at a = 0.168, z = 4 748 did not return in 15 minutes.

(PR description, Part 2)

STATUS: Partially Confirmed

The develop loop while (fabs(delta) > 1e-6) does not have an iteration cap, so a delta that stays inf keeps it running.
The tests do not run z above 700, as the PR explains.
gammaQ_fwd_shape_derivative_large_z runs (1.5, 300) and (3, 650), and on develop the derivatives are 1.2142709859604167e+112 and 2.2482228680714907e+267, against 5.7e-129 and 6.0e-277.
That matches the PR's 1.2e+112 at (1.5, 300).
At the head both are within 1e-10 relative.
The call at a = 0.168, z = 4 748 that did not return in 15 minutes is not backed by a test or script.

Evidence:

The non-termination follows from develop's code, and the test covers z up to 650.
The PR does not include code for the 15-minute call.
Code: gamma_q_test.cpp.

What I read: In develop's fwd/fun/gamma_q.hpp#L23-L35, s *= -x2.val_ / k grows like z^k / k!.
Once s overflows it stays inf, so delta stays inf and fabs(delta) > 1e-6 never fails.
The second loop at L56 is the same.
The reviewed commit removes both loops.

What I ran: python3 runTests.py -j8 --run-all on mix/fun/gamma_q_test.cpp and mix/fun/gamma_p_test.cpp, at the head and at the merge base with the PR's tests copied in.
Output in the gist: gamma_pq-head.out, gamma_pq-base.out.

Request: Could you add the call at a = 0.168, z = 4 748 as a small program that times out?

---

FACT: The gamma_lccdf shape derivative on develop is NaN or 0 where log Q < -708

Quote:

That is NaN or 0 where log Q < -708, for example at (alpha, y) = (10, 1000) and (200, 1600) (NaN) and (150, 1200) (0).

(PR description, Part 3, What failed)

STATUS: Partially Confirmed

ProbDistributionsGamma_lccdf_upper_tail_reference checks 13 rows of gamma_lccdf with var.
With develop's headers, d/dalpha is -nan at (10, 1000) and (200, 1600), as the PR says, and also at (20, 1000), (1000, 3000), (1e6, 1.1e6), (700, 2100) and (1e4, 10100).
At (150, 1200), however, develop returns 2 against the reference 2.0837 (relative error 0.040), where the PR says 0.
Develop is also finite and not 0 at two more rows with log Q < -708: 1.6111 at (300, 1500) (relative error 5.2e-4) and 2 at (150, 600, beta = 2).
At the head all 13 rows pass.
The value 0 may come from an intermediate state of the branch, with the new roots and the old gamma_lccdf, which I did not run.

Evidence:

The tail reference test is in the PR and covers the three named points.
Code: gamma_lccdf_tail_reference.hpp, rev/prob/gamma_lccdf_test.cpp.

What I ran: python3 runTests.py -j8 --run-all on prim/prob/gamma_lccdf_test.cpp, rev/prob/gamma_lccdf_test.cpp and mix/prob/gamma_lccdf_test.cpp, at the head and at the merge base with the PR's tests copied in.
The head passes all 39 tests.
The base fails the two renamed far-tail tests, the var tail reference test and both mix tests.
Output in the gist: gamma_lccdf-head.out, gamma_lccdf-base.out.

Request: Could you check the (150, 1200) example?
On develop at 5252d51 the shape derivative there is 2, where the PR says 0.
If 0 is the value with the new roots and the old gamma_lccdf, please say so in the text.

---

FACT: The other log-scale cdfs underflow: chi_square_lccdf(2000 | 20) is -inf on develop

Quote:

The other 17 functions had three defects, measured at 18 points (a, z) of the incomplete gamma function against mpmath, with the roots of Part 1 already in place: [...]
1.
The value underflows.
11 of the 12 log-scale functions computed log(gamma_p(...)) or log(gamma_q(...)).
Where the probability is below 1e-308, the value is -inf and the gradient NaN; where it is subnormal, the value has relative error 1.5e-4 and the gradient 0.1 to 1.
Example: chi_square_lccdf(2000 | 20) is -inf; the true value is -950.6.

(PR description, Part 3, What failed)

STATUS: Partially Confirmed

The develop code of 11 of the 12 log-scale functions takes log(Pn) of gamma_p or gamma_q (for example chi_square_lccdf.hpp#L100), and gamma_lccdf is the twelfth.
In the reference test on develop, chi_square_lccdf at (2000, 20) has a non-finite value, and at the head it matches the stored -950.62299837015608.
At the subnormal rows, such as chi_square_lccdf at (2400, 300), develop's value has relative error 1.53e-4, as the PR says.
The PR measured the gradient errors of 0.1 to 1 with the roots of Part 1 already in place, which is the state of commit 8d376d5.
I did not run that commit.
With develop's roots the base run has -nan, 0.040 and 0.108 at such rows.

Evidence:

The reference test covers the named point and the subnormal rows.
The gradient range was measured at an intermediate commit.
Code: incomplete_gamma_cdfs_reference_test.cpp, incomplete_gamma_cdfs_reference.hpp.

What I read: Develop's chi_square_lcdf, chi_square_lccdf, inv_chi_square_lcdf, inv_chi_square_lccdf, inv_gamma_lcdf, inv_gamma_lccdf, scaled_inv_chi_square_lcdf and scaled_inv_chi_square_lccdf add log(Pn), gamma_lcdf takes log(Pn), and poisson_lcdf and poisson_lccdf take log(gamma_q(...)) and log(gamma_p(...)).

What I ran: python3 runTests.py -j8 --run-all on the reference test, prim/prob/poisson_ccdf_log_test.cpp and the two mix cdf tests, at the head and at the merge base with the PR's tests copied in.
The head passes all 9 tests.
The base fails the reference test, the Poisson negative-count test and the mix tests of the gamma, inv_gamma, inv_chi_square and scaled_inv_chi_square families.
Output in the gist: cdfs-head.out, cdfs-base.out.

Request: Could you add the script or the commit, 8d376d5 or another, at which the reference test shows the gradient errors of 0.1 to 1?

---

FACT: The density factor overflows, so the y gradient is NaN in the bulk at large shape

Quote:

13 functions wrote it as exp(-z) pow(z, a - 1) / tgamma(a).
tgamma overflows above 171.6 and pow once (a - 1) log z > 709, so the partial is NaN also in the bulk, for example the y gradient of chi_square_cdf, chi_square_lccdf and gamma_cdf at the median for nu = 400.

(PR description, Part 3, What failed)

STATUS: Partially Confirmed

In the reference test on develop, the y gradient is -nan for chi_square_cdf and chi_square_lccdf at (400, 400) and (800, 800), and for gamma_cdf at (200, 400, 2), which are medians at large shape.
At the head these rows pass.
On develop I count 14 files that write the density as exp(-z) pow(z, a - 1) / tgamma(a) or its Poisson form.
Thirteen are the files of chi_square, gamma_cdf, inv_chi_square, inv_gamma, and scaled_inv_chi_square.
The fourteenth is poisson_cdf.hpp#L61-L62, which computes -exp(-lambda) * pow(lambda, n) / (tgamma(n + 1.0) * Pi) * P.

Evidence:

The reference test covers the named medians.
The count of 13 files is a statement about develop's code.
Code: incomplete_gamma_cdfs_reference_test.cpp.

What I read: Develop's chi_square_cdf.hpp L101-L102 and the other files named in the finding use pow and tgamma in the density factor.

What I ran: python3 runTests.py -j8 --run-all on rev/prob/incomplete_gamma_cdfs_reference_test.cpp, prim/prob/poisson_ccdf_log_test.cpp and the two mix cdf tests, at the head and at the merge base with the PR's tests copied in.
Output in the gist: cdfs-head.out, cdfs-base.out.

Request: Could you check the count of 13?
poisson_cdf on develop has the same form for the lambda partial.

---

FACT: The worst errors at the 324 points are 2e-13 for values, 8e-13 for gradients, and 2.7e-10 for one partial

Quote:

The worst relative error is 2e-13 for values and 8e-13 for gradients, with one exception: the nu partial of scaled_inv_chi_square near the mode has a relative error up to 2.7e-10.

(PR description, Part 3, Accuracy)

STATUS: Partially Confirmed

At the head the reference test passes with a relative tolerance of 1e-12 for values and 1e-9 for gradients, so the errors are below those bounds.
The test does not print the worst errors, so the numbers 2e-13, 8e-13 and 2.7e-10 are not shown by any output.
On develop the largest relative gradient error among the finite results in the base run is 1.6353632961763142e+93, which matches the Summary row "18 cdfs, worst relative gradient error where finite" (1.6e+93 against 2.7e-10).
The stored references come from mpmath code that is not in the PR.

Evidence:

The test checks the bounds 1e-12 and 1e-9.
The worst errors and the mpmath code are not in the PR.
Code: incomplete_gamma_cdfs_reference_test.cpp.

What I ran: python3 runTests.py -j8 --run-all on rev/prob/incomplete_gamma_cdfs_reference_test.cpp, prim/prob/poisson_ccdf_log_test.cpp and the two mix cdf tests, at the head and at the merge base with the PR's tests copied in.
I took the largest "relative error" on a gradient line in cdfs-base.out.
Output in the gist: cdfs-head.out, cdfs-base.out.

Request: Could you add the mpmath script that produced incomplete_gamma_cdfs_reference.hpp by the two routes?
A way to print the worst errors would also let others reproduce the numbers in this paragraph.

---

FACT: gamma_lccdf in all autodiff modes: first derivatives within 5.6e-12, second within 1.9e-10

Quote:

For gamma_lccdf in all autodiff modes, first derivatives are within 5.6e-12 and second derivatives within 1.9e-10 (7e-9 near y = alpha + 1 at alpha = 5).

(PR description, Part 3, Accuracy)

STATUS: Partially Confirmed

The mix tail test checks first derivatives with a relative tolerance of 1e-11 and second derivatives with 1e-9, or 5e-8 when y < alpha + 3.
It passes at the head, so the errors are below those bounds.
The test does not print the worst errors, so 5.6e-12, 1.9e-10 and 7e-9 are not shown by any output.

Evidence:

The test checks bounds that are looser than the stated numbers.
Code: gamma_lccdf_tail_reference.hpp, mix/prob/gamma_lccdf_test.cpp.

What I ran: python3 runTests.py -j8 --run-all on the prim, rev and mix gamma_lccdf_test.cpp files, at the head and at the merge base with the PR's tests copied in.
Output in the gist: gamma_lccdf-head.out, gamma_lccdf-base.out.

Request: Could you add the mpmath script behind gamma_lccdf_tail_reference.hpp, and the code that measured these worst errors?

---

FACT: The cdf values change only where they were -inf or lost digits

Quote:

The values of the distributions change only where they were -inf or lost digits; elsewhere only the gradients change.

(PR description, Summary)

STATUS: Partially Confirmed

In the base run, 66 of the 324 rows fail on the value.
Each of them is non-finite on develop, or has a relative error of 1.08e-9 or more, which is the loss of digits the PR describes.
The other 258 values agree with the reference within 1e-12 on develop.
This covers the 324 test points only, and the PR does not include code that compares values over a wider range.

Evidence:

The reference test covers 324 points.
Code: incomplete_gamma_cdfs_reference_test.cpp.

What I ran: python3 runTests.py -j8 --run-all on rev/prob/incomplete_gamma_cdfs_reference_test.cpp, prim/prob/poisson_ccdf_log_test.cpp and the two mix cdf tests, at the head and at the merge base with the PR's tests copied in.
I read the "value of" lines in cdfs-base.out.
Output in the gist: cdfs-head.out, cdfs-base.out.

Request: If the statement rests on more than the 324 points, could you add the code that compared the values?

---

FACT: The new tests fail without the changes

Quote:

The new tests fail without the changes.
I installed the old headers, kept the new tests, and restored afterwards; the header checksums after the restore equal those before:

(PR description, Testing)

STATUS: Partially Confirmed

I put the PR's test files on develop's headers, all at once, and ran them.
Every new test fails there except mathMixScalFun_chi_square_cdfs and mathMixScalFun_poisson_cdfs, the expect_ad tests of two families, which pass on develop.
The failing tests of the roots are the 5 accuracy tests, gradRegIncGamma_gradient_matches_hessian, and the four tests with the changed constant (gradRegIncGamma_typical, _fd, _ffd and _fv).
The other failing tests are the two forward-mode gamma_q tests, gammaP_rev_far_upper_tail, the var and mix tail tests, and the mix expect_ad tests of gamma_lccdf, gamma, inv_gamma, inv_chi_square and scaled_inv_chi_square.
The reference test, the Poisson negative-count test, and the two renamed far-tail tests fail too.
The PR's table installs old headers one part at a time, so its counts are not comparable with mine.
For example, the reference test fails in all 18 functions in my run, with 569 failing expectations against the PR's 421 in 17.

Evidence:

The tests are in the PR.
The per-part header swaps of the table are not.
Code: incomplete_gamma_cdfs_test.cpp, inv_incomplete_gamma_cdfs_test.cpp.

What I ran: The four runTests.py -j8 --run-all runs at the merge base with the PR's 13 test files copied in, and the same four at the head.
The gist README lists the files of each run.
Output in the gist: roots-base.out, gamma_pq-base.out, gamma_lccdf-base.out, cdfs-base.out.

Request: Could you say in the text that the expect_ad tests of chi_square and poisson in incomplete_gamma_cdfs_test.cpp also pass on develop?
If the per-part counts in the table matter, a script that installs each set of old headers would let others rerun them.

---

FACT: The true d/da Q(0.5, 1.0) is 0.38983726432851057, and the old constant was off

Quote:

• Three constants in test/unit/math/mix/fun/grad_reg_inc_gamma_test.cpp expected 0.38984156 for d/da Q(0.5, 1.0); the true value is 0.38983726432851057.

(PR description, Testing)

STATUS: Partially Confirmed

In gradRegIncGamma_typical, _fd, _ffd and _fv, develop returns 0.38984156 (0.38984153 for fvar<var>), and the head matches 0.38983726432851057 under EXPECT_FLOAT_EQ.
gradRegIncGamma_fv_1stderiv and gradRegIncGamma_fv_2ndderiv, which the PR did not change, pass at the head.
The PR changes four constants, and the text says three.
The value 0.38983726432851057 is from mpmath, and the PR does not include that code.
The split of the old error between the rounded g and the series truncation is not shown by any test.

Evidence:

The test is in the PR.
The mpmath code behind the new constant is not.
Code: grad_reg_inc_gamma_test.cpp.

What I ran: python3 runTests.py -j8 --run-all on prim/fun/grad_reg_inc_gamma_accuracy_test.cpp and mix/fun/grad_reg_inc_gamma_test.cpp, at the head and at the merge base with the PR's tests copied in.
Output in the gist: roots-head.out, roots-base.out.

Request: Could you change "Three" to "Four"?
As with the other references, could you also add the mpmath line that gives 0.38983726432851057?

---

FACT: gamma_lcdf(20000 | 800, 0.1) is now about -1e-204, so log1m_exp(lcdf) equals lccdf

Quote:

gamma_lcdf now returns its true value, about -1e-204, so the case checks that log1m_exp(gamma_lcdf) equals gamma_lccdf there, and is renamed to lccdf_far_tail_matches_log1m_exp_lcdf.

(PR description, Testing)

STATUS: Partially Confirmed

lccdf_far_tail_matches_log1m_exp_lcdf in prim and rev asserts lcdf < 0 and that log1m_exp(lcdf) is within 1e-12 relative of gamma_lccdf at that point.
On develop lcdf is 0 and log1m_exp(lcdf) is -inf against an lccdf of -471.636, so both tests fail.
At the head both pass.
The tests do not print or assert the value of the lcdf, so "about -1e-204" is not shown.

Evidence:

The PR changes the prim and rev tests at this point.
Code: prim/prob/gamma_lccdf_test.cpp, rev/prob/gamma_lccdf_test.cpp.

What I ran: python3 runTests.py -j8 --run-all on the prim, rev and mix gamma_lccdf_test.cpp files, at the head and at the merge base with the PR's tests copied in.
Output in the gist: gamma_lccdf-head.out, gamma_lccdf-base.out.

Request: Could you assert or print the lcdf value of -1e-204 at this point?

---

FACT: The generated third-order chi_square case failed with the stale gamma_q copy and passes without it

Quote:

With an accurate root and the stale copy, the case fails by 2.4e-03; with the copy removed, it passes.

(PR description, Part 2)

STATUS: Partially Confirmed

test/prob/chi_square/chi_square_ccdf_log_test.hpp computes the reference through log(gamma_q(nu * 0.5, y * 0.5)).
The three programs generated from it pass at the head.
The failure by 2.4e-03 is in a state with the new roots and the old forward-mode gamma_q, which is the PR's commit b172d61.
I ran only the reviewed commit and the merge base.

Evidence:

The passing half is shown by the generated tests.
The failing half is at an intermediate commit.
Code: chi_square_ccdf_log_test.hpp.

What I ran: make -j8 generate-tests -s, then python3 runTests.py -j8 --run-all on the generated cdf programs of the six families at the head, and the summary rerun with -t.
prob-cdfs-summary-head.out lists the three chi_square_ccdf_log_00000_generated_*_test programs among the 90 it ran, and it has 90 PASSED lines and no FAILED line.
Output in the gist: prob-cdfs-summary-head.out.

Request: Is b172d61 the state that shows the 2.4e-03 failure?
If so, could you say so in the text?

---

[stevieb-bot]

Fact check, part 2 of 2.

Claims and Evidence (continued)

Not Verified Facts

FACT: The 1 800-point sweep of the roots

Quote:

A sweep of 1 800 points in nine regions chosen by the branch structure; worst relative error per region (200 points each):

(PR description, Part 1, Accuracy)

STATUS: Not Verified

The PR does not include the code of this sweep, so I could not rerun it.
The new accuracy test makes 18 checks at fixed points, which are not the 1 800 points of the table.
The sentence names nine regions of 200 points, and the table has seven rows.

Evidence:

The PR does not include the script of the table or the mpmath reference for the sweep.

Request: Could you add the script that draws the 1 800 points, computes the mpmath reference, and prints the worst relative error per region, for example under a scripts/ or test/ folder or linked from the PR text?

---

FACT: The 49 990-point sweep: 6 379 non-finite on develop and none with the PR

Quote:

A second sweep compared 49 990 points with a and z both in [1e-3, 1e4], dense on z ≈ a, the z = 8 seam, a = 171.62 and a log z = 709.
develop returns NaN or inf at 6 379 of them for grad_reg_inc_gamma; this PR is finite at every one, with median relative errors between 1e-16 and 1e-14 in every region.

(PR description, Part 1, Accuracy)

STATUS: Not Verified

The PR does not include the code of this sweep.
The first Summary row ("roots, non-finite results at 49 990 points", 6 379 against 0) matches the numbers of this sweep.
The second row ("worst relative error outside a ≈ z ≈ 1e4", 3e-11) matches the 3e-11 of the limit near z = a.

Evidence:

The PR does not include the script of the sweep.

Request: Could you add the script of the 49 990-point sweep, with the reference it compares against?

---

FACT: The roots are 2.2x to 4.4x faster with double, and the upper root is 0.72x with var

Quote:

Function level, net ns per call, Xeon E5-2680 v3, performance governor, two rounds (spread under 3 %), on a in [0.5, 60] and z in [0.1, 40], where develop terminates and is finite:

(PR description, Part 1, Speed)

STATUS: Not Verified

The PR does not include the benchmark behind this table, so I could not rerun it.
The Summary row "speed, double roots" comes from the same table.

Evidence:

The PR does not include a benchmark file or script.

Request: Could you add the benchmark, for example a google-benchmark file that times both roots with double, fvar<double>, fvar<fvar<double>> and var on the stated grid?
Please also say how the develop build was made, so it can be run against develop on the same machine.

---

FACT: The 18 cdfs are 1.9x to 12x faster with var

Quote:

Value plus all gradients with var arguments, ns per call, AMD EPYC 9965, g++ 11.5, -O3, medians of 4 alternating runs.
Points with a = 5 and a = 50 and z from 0.2 a to 3 a, mapped to each function's arguments:

(PR description, Part 3, Speed)

STATUS: Not Verified

The PR does not include the benchmark behind this table, so I could not rerun it.
The Summary row "speed, cdf value plus gradients (var)" comes from the same table.

Evidence:

The PR does not include a benchmark file or script.

Request: Could you add the benchmark code for the cdf table, with the points it maps to each function's arguments?

---

FACT: The gamma_lccdf shape derivative alone is 1.6x to 34.9x faster

Quote:

gamma_lccdf alone, shape derivative only, with alpha from 5 to 1000 (ns per call):

(PR description, Part 3, Speed)

STATUS: Not Verified

The PR does not include the benchmark behind this table, so I could not rerun it.

Evidence:

The PR does not include a benchmark file or script.

Request: Could you add the benchmark code for this table, with the points of each region and the autodiff types?

---

FACT: On the grid of the plot by @mjcarter95, the worst relative error is 3.4e-12

Quote:

On the grid of the plot by @mjcarter95, all 266 cells are finite and the worst relative error of the shape derivative is 3.4e-12.

(PR description, Part 3, Accuracy)

STATUS: Not Verified

The PR does not include the grid or the script that computes the 266 cells.
The new tail reference test checks 13 points of gamma_lccdf, which are not the grid of the plot.
In the conversation @mjcarter95 compared the fixed code against mpmath, and the author replied that @mjcarter95 reports a worst relative error of 4.1e-6 against the PR's 3.4e-12. @mjcarter95 answered that every cell used precision 60.
A shared script would resolve the difference.

Evidence:

The grid and the mpmath reference are not included.

Request: Could you add the script that computes the 266 cells and the reference, or link it?

---

Description vs Code

Each section below lists facts from the PR summary and discussion that an agent was able to confirm, partially confirm, or not confirm at all, or that a run of the code in the PR contradicted.
Please review the "Partially Confirmed", "Not Verified", and "Contradicted" facts for facts missing evidence found by this agent.

Confirmed Facts

FACT: The roots call Eigen's igamma_der_a, which is vendored

Quote:

The gradient roots grad_reg_inc_gamma and grad_reg_lower_inc_gamma (the derivatives of Q(a, z) and P(a, z) with respect to a) are replaced by Eigen::numext::igamma_der_a, which is already vendored in lib/eigen_5.0.1.

(PR description, Summary)

STATUS: Confirmed

At the reviewed commit grad_reg_inc_gamma returns -Eigen::numext::igamma_der_a(TP(a), TP(z)) and grad_reg_lower_inc_gamma returns Eigen::numext::igamma_der_a(TP(a), TP(z)).
The function is defined in lib/eigen_5.0.1/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h, which is the same on develop.
The series, asymptotic, and Gautschi branches of the old roots are deleted.

---

FACT: Two new headers supply the three Eigen helpers for autodiff types

Quote:

This PR supplies those three for Stan's autodiff scalars in two new headers, stan/math/fwd/fun/Eigen_SpecialFunctions.hpp and stan/math/rev/fun/Eigen_SpecialFunctions.hpp.

(PR description, Part 1, The change)

STATUS: Confirmed

fwd/fun/Eigen_SpecialFunctions.hpp specializes lgamma_impl, digamma_impl and cephes_helper for fvar<T>, and rev/fun/Eigen_SpecialFunctions.hpp does the same for var.
Both are inside namespace Eigen::internal and are included from stan/math/fwd/fun.hpp and stan/math/rev/fun.hpp.
stan/math/rev/core/Eigen_NumTraits.hpp already opens namespace internal inside Eigen on develop.

---

FACT: The forward-mode gamma_q calls grad_reg_inc_gamma and evaluates d/dz in log space

Quote:

This PR replaces the copy with a call to grad_reg_inc_gamma, as fwd/fun/gamma_p.hpp already calls grad_reg_lower_inc_gamma.

(PR description, Part 2)

STATUS: Confirmed

On develop, both fvar overloads of gamma_q that differentiate in a have their own series loop, while (fabs(delta) > 1e-6), without an iteration cap.
At the reviewed commit the loops are gone, and der1 comes from grad_reg_inc_gamma.
der2 is -exp(-z + (a - 1) log z - lgamma(a)) in both overloads that need it.
fwd/fun/gamma_p.hpp already calls grad_reg_lower_inc_gamma on develop.

---

FACT: The old root's asymptotic branch reached the fwd gamma_q loop

Quote:

The asymptotic branch of the old root calls gamma_q with fvar arguments, so grad_reg_inc_gamma reached the loop too.
The existing mix test for gamma_q carries a commented-out case with the note "this still fails forward mode", which is the same defect.

(PR description, Part 2)

STATUS: Confirmed

On develop the asymptotic branch of grad_reg_inc_gamma returns gamma_q(a, z) * (l - dig) + ..., with a and z of the caller's types, so an fvar call goes to the forward-mode gamma_q.
The develop gamma_q_test.cpp has the commented-out case with the note "this still fails forward mode".

---

FACT: The z / a > 10 short cut is removed from all three gamma_p overloads

Quote:

With the new roots the derivatives are accurate in that region, and this PR removes the short-cut from all three overloads.

(PR description, Part 2)

STATUS: Confirmed

On develop, gamma_p_vv_vari, gamma_p_vd_vari and gamma_p_dv_vari each return without derivatives when the ratio is above 10 (L37-L41, L65-L69, L89-L93).
At the reviewed commit the three checks are deleted, and rev/fun/gamma_p.hpp goes straight to the derivatives.

---

FACT: log_gamma_tail chooses the continued fraction or gamma_p and gets d/dz from the density

Quote:

One internal function, log_gamma_tail(a, z) in log_gamma_q_dgamma.hpp, returns log P or log Q with its derivatives with respect to a and z

(PR description, Part 3, The change)

STATUS: Confirmed

In log_gamma_q_dgamma.hpp, internal::log_gamma_tail takes the continued fraction log_q_gamma_cf_and_da when z > a + 1 and sets log P = log1m_exp(log Q).
Otherwise it calls gamma_p and grad_reg_lower_inc_gamma, and switches to the log-space series log_p_gamma_series_and_da when P < DBL_MIN / EPSILON (L282).
The continued fraction stops when both delta and d_delta are within the tolerance (L144-L148), and the series stops on the bound t q / (1 - q) of the remaining terms.
d/dz is exp(log_density - log_prob) (L295-L296).

---

FACT: Each of the 18 functions keeps its argument checks, and the cdfs multiply by exp(sum) at the end

Quote:

Each of the 18 functions keeps its argument checks and applies its own chain rule.
The probability-scale cdf sums the logs and partials and multiplies by exp(sum) at the end, so it has neither the density overflow nor the 0/0.

(PR description, Part 3, The change)

STATUS: Confirmed

All 18 headers call internal::log_gamma_tail.
In each of the 18 files the check_* calls are the same at both commits.
The only change is the label "Scale Parameter" to "Scale parameter" in check_consistent_sizes of the three inv_gamma files.
In gamma_cdf.hpp the loop adds tail.log_prob to log_P and the partials, then computes P = exp(log_P), multiplies each partial by P, and returns P.

---

FACT: An element with probability 1 now contributes 0 instead of ending the loop

Quote:

• gamma_lcdf and chi_square_lcdf at y = inf, and the lccdf of chi_square, inv_gamma, inv_chi_square and scaled_inv_chi_square at y = 0, returned 0 for the whole vector when one element had probability 1, so the other elements were lost.
  Such an element now contributes 0.

(PR description, Part 3, The change)

STATUS: Confirmed

On develop, gamma_lcdf and chi_square_lcdf return ops_partials.build(0.0) for the whole vector when one y is infinite, and the four lccdfs do the same when one y is 0.
poisson_lccdf returns 0 when any count is negative.
At the reviewed commit each of these is a continue for that element, for example gamma_lcdf.hpp#L85-L86, chi_square_lccdf.hpp#L79-L80 and poisson_lccdf.hpp#L77-L78.
poisson_cdf and poisson_lcdf skip n = INT_MAX (L84, L82), and poisson_lccdf returns -inf for it (L69).

---

FACT: OpenCL code does not call these functions

Quote:

OpenCL is unaffected: no OpenCL code calls these functions, and the six families have OpenCL implementations only for their lpdf or lpmf.

(PR description, Not included)

STATUS: Confirmed

A search of stan/math/opencl at the reviewed commit does not find a call to grad_reg_inc_gamma, grad_reg_lower_inc_gamma, gamma_p or gamma_q.
The OpenCL files for the six families are chi_square_lpdf, gamma_lpdf, inv_chi_square_lpdf, inv_gamma_lpdf, scaled_inv_chi_square_lpdf, poisson_lpmf, poisson_log_lpmf and poisson_log_glm_lpmf.

---

FACT: log P and log Q are internal, with no new public function

Quote:

• Related to Add Boost's tgamma_lower and tgamma(a,z) functions #2769: log P and log Q are now computed in log space inside the library, but no public function is added.

(PR description, Issues)

STATUS: Confirmed

log_gamma_tail is in namespace internal in log_gamma_q_dgamma.hpp.
A search of stan/ at the reviewed commit finds no log_gamma_p or log_gamma_q.

---

FACT: log_gamma_q_dgamma is no longer called by the library

Quote:

log_gamma_q_dgamma itself is no longer called by the library and is kept for compatibility.

(PR description, Not included)

STATUS: Confirmed

On develop the only caller is gamma_lccdf.hpp#L42.
At the reviewed commit a search of stan/ finds log_gamma_q_dgamma only at its definition, log_gamma_q_dgamma.hpp#L327.

---

FACT: Other discrete cdfs still let one element decide the whole result

Quote:

Other discrete cdfs have the defect of #56: one element with probability 1 decides the result for the whole vector.

(PR description, Not included)

STATUS: Confirmed

At the reviewed commit, bernoulli_lccdf, binomial_lccdf, beta_binomial_lccdf, neg_binomial_lccdf, beta_neg_binomial_lccdf and yule_simon_lccdf return 0 for the whole call when one count is below the support, and neg_binomial_2_lccdf calls neg_binomial_lccdf.
neg_binomial_lcdf, beta_neg_binomial_lcdf and yule_simon_lcdf return 0, and neg_binomial_cdf, neg_binomial_2_cdf, beta_neg_binomial_cdf and yule_simon_cdf return 1, when one count is INT_MAX.

---

FACT: The existing lower-root test grid has no point of defect 3

Quote:

The existing grad_reg_lower_inc_gamma_test checks 9 600 values from Mathematica with a 1e-10 absolute tolerance; Eigen satisfies all of them, but that grid (a up to 19.75, z up to 29.75) contains no point of defect 3.

(PR description, Testing)

STATUS: Confirmed

grad_reg_lower_inc_gamma_test.cpp has n_tests = 9600, with a up to 19.7501 and z up to 29.7501, so it does not have a point with z > 42.
The tolerance is 1e-10 except for the first point, which uses 1.5e-5.

---

FACT: On develop, every point with z > 42 and a >= 12 goes to the Gautschi branch of the lower root

Quote:

Its delegation test sends every point with z > 42 and a >= 12 to the Gautschi branch, where emz * (log_z * sum_a - sum_b) is a difference of two large, nearly equal sums.

(PR description, Part 1, The defects on develop)

STATUS: Confirmed

Develop's grad_reg_lower_inc_gamma.hpp#L131-L137 delegates to grad_reg_inc_gamma when a < 0.8 && z > 15, a < 12 && z > 30, or a < sqrt(-756 - z^2 + 60 z). For z > 42 the argument of the square root is negative, so the third test is false, and with a >= 12 the first two are false too.
Such a point reaches the Gautschi branch, which returns emz * (log_z * sum_a - sum_b).

---

FACT: The Poisson cdfs treat n = INT_MAX as an infinite count, as three other families do

Quote:

• The Poisson cdfs treat n = INT_MAX as an infinite count, as neg_binomial, beta_neg_binomial and yule_simon already do.

(PR description, Part 3, The change)

STATUS: Confirmed

At the reviewed commit poisson_cdf and poisson_lcdf skip an element with n = INT_MAX (L84, L82), and poisson_lccdf returns -inf for it (L69).
neg_binomial_lcdf, beta_neg_binomial_lcdf and yule_simon_lcdf test for INT_MAX too (L92, L85, L61).

---

Partially Confirmed Facts

FACT: The root signatures are unchanged and the extra arguments are ignored

Quote:

The signatures do not change.
g, dig, precision and max_steps are accepted and ignored.

(PR description, Part 1, The change)

STATUS: Partially Confirmed

The signatures match develop.
grad_reg_inc_gamma(T1 a, T2 z, T1 g, T1 dig, double precision = 1e-6, int max_steps = 1e5) and grad_reg_lower_inc_gamma(const T1& a, const T2& z, double precision = 1e-10, int max_steps = 1e5) are the same at both commits.
precision and max_steps are unused.
g and dig are still read, though, because grad_reg_inc_gamma.hpp#L43 returns NaN when either of them is NaN.

Request: Could you say that g and dig are still checked for NaN?

---

FACT: There are two behavior changes at the edges of the support

Quote:

Two behavior changes at the edges of the support:

(PR description, Part 3, The change)

STATUS: Partially Confirmed

The list under this sentence has three items, and the code has all three changes (see the fact check on elements with probability 1).

Request: Could you change "Two" to "Three"?

---

FACT: On develop, z >= a && z >= 8 selects a fixed asymptotic expansion

Quote:

The test z >= a && z >= 8 selects a fixed 10-term asymptotic expansion (DLMF 8.11.2), which needs z >> a.

(PR description, Part 1, The defects on develop)

STATUS: Partially Confirmed

Develop's grad_reg_inc_gamma.hpp#L65-L89 takes the asymptotic branch when z >= a && z >= 8, with the comment "asymptotic expansion http://dlmf.nist.gov/8.11#E2", and sums a fixed number of terms without a convergence test.
The loop is for (int k = 1; k < 10; ++k) with one S += delta per pass, so it adds 9 terms.

Request: The loop adds 9 terms.
Could you check "10-term"?

---

FACT: The mix tests run expect_ad for all 18 functions at two points each, one on each branch

Quote:

• mix/prob/incomplete_gamma_cdfs_test.cpp and mix/prob/inv_incomplete_gamma_cdfs_test.cpp: expect_ad for all 18 functions at two points each, one on each branch.

(PR description, Testing)

STATUS: Partially Confirmed

The two mix files call expect_ad for all 18 functions at two points each.
log_gamma_tail takes the continued fraction only when z > a + 1.
For the gamma family, z = beta * y and a = alpha, so the points (2.5, 3.0, 1.5) and (1.0, 4.0, 1.5) give z = 3.75 with a + 1 = 4, and z = 1.5 with a + 1 = 5.
For chi_square, z = y / 2 and a = nu / 2, so (5.0, 3.0) and (2.0, 6.0) give z = 2.5 with a + 1 = 2.5, and z = 1 with a + 1 = 4.
Both points of these two families are therefore on the gamma_p side.
The poisson, inv_gamma, inv_chi_square, and scaled_inv_chi_square families have one point on each side.

Request: Could you move one point of the gamma and chi_square checks to z > a + 1, or change the comment and the text?

---

Not Verified Facts

FACT: The poisson_log cdfs are now short wrappers of the log-space function

Quote:

The missing poisson_log_cdf, poisson_log_lcdf and poisson_log_lccdf (#90, #2657) are now short wrappers of the same log-space function, log Q(n + 1, exp(eta)); they need new signatures in stanc3 and are left for a separate PR.

(PR description, Not included)

STATUS: Not Verified

The reviewed commit does not have poisson_log_cdf, poisson_log_lcdf or poisson_log_lccdf.
A search of stan/ does not find them, and stan/math/prim/prob/ has only poisson_log_lpmf, poisson_log_rng and poisson_log_glm_lpmf.
The sentence also says they are left for a separate PR, and #3476 adds them.

Request: Could you reword this to say that they can be written as short wrappers, and that #3476 adds them?

---

Not Checked

• The 49 990-point sweep, the 1 800-point sweep, the autodiff-order spread, the accuracy limit near z = a, and the plot grid of 266 cells.
  The PR does not include the scripts or the mpmath code behind them.
• The speed table of the roots and the two speed tables of the cdfs.
  The PR does not include the benchmark code.
• The script that produced the mpmath reference values in the new tests.
  It is not in the PR, so I could not rerun it.
  I ran the tests that use the values.
• The red table in Testing installs old headers one part at a time.
  My base runs put all of develop's headers under the new tests at once, so the failure counts are not comparable row by row.
• The reason the new headers are in the fun layer (that in core they break the include order of every mix target).
  I did not build that alternative.
• The test suites with g++ 11.5 and g++ 14.3.
  The container has g++ 13.3.
• Intermediate commits of the branch (b172d61, the commits before 9341c4c, and 8d376d5) behind the 2.4e-03 failure, the 0 at (150, 1200), and the 0.1 to 1 gradient errors.
  I ran only the reviewed commit and the merge base.
• Whether reverse-mode code in the library calls the roots only with double values.
  I did not search every call site.

[avehtari]

I updated the PR text, added tests and published the scripts behind the numbers. The PR had a conflict with develop, so I merged develop at 510d263078. Below, "develop" means 5252d51d47, the merge base before the merge, as in your runs. The items follow the order of the fact check.

Scripts

The scripts are in https://github.com/avehtari/stan-pr-checks/tree/27838b5236/pr3408. The README gives one command per claim, the expected output and the commit I measured it on. The directory has:

• the mpmath generators of grad_reg_inc_gamma_accuracy_test.cpp, gamma_lccdf_tail_reference.hpp and incomplete_gamma_cdfs_reference.hpp (two routes each), and of the constant 0.38983726432851057;
• the 1 800-point and 49 990-point sweeps of the roots, with the autodiff-order spread;
• the survey of the 18 cdfs at the 324 points, the gamma_lccdf autodiff modes and the 266-cell grid;
• the forward-mode hang as a program under timeout, and a swap-test script;
• the three speed benchmarks with their job scripts.

I ran the accuracy measurements again on develop and on the merged head. I did not run the speed benchmarks again.

Numbers that changed in the rerun

• develop is non-finite at 7 379 of the 49 990 points. The PR text said 6 379; the per-region counts sum to 7 379, and the old total was an addition error.
• In the 1 800-point sweep, the column of this PR changed in three regions: 1.2e-13, 5.4e-12 and 2.5e-11 (before: 1.9e-13, 3.8e-12, 4.3e-12). I had changed the mpmath reference script after I made the table. The two larger values come from the reference. Against mp.diff at all 200 points of each of those regions, the worst error of this PR is 4.5e-14 and 9.4e-14. The table now has all nine regions.
• Outside the band near z = a, the worst error against the reference is 3e-11. At the worst point of each region this PR is within 1e-13 of mp.diff, so the PR text now calls the 3e-11 reference-limited.
• The "3e-11 for a < 1 000" near z = a was also the reference: at its worst point this PR is within 6e-14 of mp.diff. I removed it from the stated limit. The 2.9e-7 and the 5.2e-6 stay. At the worst point near a = z = 1e4, mp.diff, the quadrature and a series (ref_mpmath.dPda_C) agree to 16 digits, and the scripts include all three.
• The worst value error at the 324 points is 2.4e-13; the text said 2e-13.

Partially confirmed and not verified facts

• These items now have scripts (see Scripts above):
  • the references of the accuracy test and of the far upper tail of the lower root;
  • the error near z = a up to a = 1e4, and the autodiff-order spread;
  • the worst errors at the 324 points and the gamma_lccdf modes;
  • the 266-cell grid and the speed tables.
• (150, 150) is in the asymptotic branch on develop, as you found. In the series branch, when a log z > 709, develop returns -inf or +inf; at (170, 100) it returns -inf, and the true value is 6.6e-11. The text now uses that point. grad_reg_inc_gamma_series_region tests it and (165, 80).
• grad_reg_inc_gamma_series_region also checks (0.5, 0.25) and (8, 4) with an exact g. develop fails there with 9.2e-8 and 9.0e-5.
• The reference test has 10 new rows, marked root_table, at the points of the Part 1 table: the five finite rows and the five NaN points. Each of the 10 rows fails on develop.
• The 5.8e-01 of the rejected double-only version is the error of the develop root at a = z = 100, which hessian() still used in that version. The accuracy test shows that error (0.577). I did not keep the prototype.
• gammaQ_fwd_log_upper_tail_composite is the second example of gamma_p and gamma_q gradient problems #2006: the fvar<double> derivative of log(gamma_q(alpha, beta * y)) at y = 12, alpha = 3, beta = 3 with all tangents 1. develop returns -6.7e+12, and the reference is -11.5005397823109.
• The scripts run the call at a = 0.168, z = 4 748 under timeout 60. The timeout stops develop, and this PR returns. I also enabled the commented-out case expect_value(f, 8.01006, 2.47579e+215) of the gamma_q mix test; develop hangs on it.
• At (150, 1200), develop returns 2, as you found; the true value is 2.084. The 0 in the text came from commit b172d61, which has the Part 1 roots and the old gamma_lccdf. The text now says this.
• I measured the gradient errors of 0.1 to 1 at commit 8d376d5. At the subnormal rows the finite gradient errors there are 0.108 and 1.0, and the other partials are NaN. The text now says "NaN or 0.1 to 1" and names the commit.
• The density factor is in 14 functions, including poisson_cdf. On develop its lambda gradient at poisson_cdf(n | n) is -inf from n = 150 and NaN from n = 172.
• The scripts compare the develop values with the reference at the 324 points: 28 are -inf, 38 have a relative error of 1e-9 or more, and the other 258 agree to 1e-12.
• I replaced the per-part swap table with one run that installs all develop headers at once, as you did. 23 of the 24 new test cases fail there. mathMixScalFun_poisson_cdfs passes on develop, and the text says so. mathMixScalFun_chi_square_cdfs now fails on develop because it has a new point (see Description vs code).
• Four constants changed; the text said three.
• The gamma_lcdf value is -1.48e-205, and the prim and rev cases now assert it.
• The 2.4e-03 failure is commit b172d61 (new roots, old forward-mode gamma_q). The text names it, and af3205c for the fix.

Description vs code

• The text now says that the roots check g and dig for NaN.
• There are three behavior changes at the edges of the support, and the loop of the asymptotic branch adds 9 terms. I fixed both counts (the 9 in two places).
• The gamma and chi_square checks in incomplete_gamma_cdfs_test.cpp have a new point with z > a + 1, (4, 2, 1.5) and (8, 3). The old points stay.
• The text now says that the poisson_log cdfs can be written as short wrappers, and that Add poisson_log_cdf, poisson_log_lcdf and poisson_log_lccdf #3476 adds them.

Other changes

I renamed grad_reg_inc_gamma_unchanged_points to grad_reg_inc_gamma_small_a, because two of its three points have z > a, where Eigen uses the continued fraction. gamma_cdf.hpp no longer includes the unused value_of.hpp, as on develop.

On the merged head with the new tests (Triton, g++ 14.3), these pass: 35 unit suites, the 123 generated test/prob programs of the six families (90 of them cdf, lcdf or lccdf), cpplint, make test-math-dependencies, the header tests of the 29 changed headers, and doxygen. Test changes: 8e50c47.

@SteveBronder, let me know if you want the responses for @stevieb-bot to have some specific structure