[avehtari]

EDIT: This PR text was updated 2026-10-21

This PR was Claude assisted

Issues

In addition of great accuracy, stability and speed improvement

• 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 13 of them returned NaN 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	6 379	0
roots, worst relative error outside a ≈ z ≈ 1e4	inf	3e-11
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 10-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. The series branch fails the same way through exp(a * log z) when a log z > 709, for example at (a, z) = (150, 150).
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 = 0.5 and 9.0e-05 at a = 8. 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 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. g, dig, precision and max_steps are accepted and ignored.

Accuracy

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

region	develop grad_reg_inc_gamma	develop grad_reg_lower_inc_gamma	this PR
series, a < 1	1.5e-03	1.5e-10	8.5e-15
series, 10 < a <= 170	2.6e+29	5.8e-01	1.9e-13
asymptotic, a <= z <= 2a	5.1e-01	1.1e+07	1.7e-13
asymptotic, z >= 3a	3.7e-04	2.1e+165	3.8e-12
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	4.3e-12

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.

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. Three independent references agree to 17 digits there, so this is the algorithm, not the measurement. 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 %), 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 10-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 carries a commented-out case with the note "this still fails forward mode", which is the same defect.

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 cancelled. With an accurate root and the stale copy, the case fails by 2.4e-03; with the copy removed, 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 or 0 where log Q < -708, for example at (alpha, y) = (10, 1000) and (200, 1600) (NaN) and (150, 1200) (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:

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.
2. The density factor in the y and scale partials overflows. 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.
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.

Two 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 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. 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. 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, 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.
• test/unit/math/rev/prob/incomplete_gamma_cdfs_reference.hpp with incomplete_gamma_cdfs_reference_test.cpp: 324 rows, the value and gradient of the 18 functions at 18 points, 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:

• 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. 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, 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.

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:

old headers	result
roots and gamma_q of develop	2 of 4 suites fail, 10 cases
gamma_lccdf before its commits	1 of 18 var cases fails (8 rows); both mix cases fail (162 expectations)
the 18 cdfs before their commit	the reference test fails with 421 expectations, in all 17 changed functions
rev/fun/gamma_p.hpp and poisson_lccdf before their commit	gammaP_rev_far_upper_tail fails 7 expectations (all three overloads and the case of #2006); lccdf_negative_count_contributes_zero fails both vector cases of #56

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. Tested on develop at e50e79f188 with this PR applied, with g++ 11.5 and g++ 14.3.

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 the repair is left for a separate PR.

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.

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