[avehtari]

Fixes #2031. Fixes #2923. Fixes #3234. Fixes #3377. Fixes other functions. Partially addresses #2790.

Currently this PR is based on stable-inc-beta branch PR #3409, so that PR should be merged first.

This improves beta_lccdf, beta_proportion_lccdf, binomial_lccdf, neg_binomial_lccdf, neg_binomial_2_lccdf, student_t_cdf, student_t_lcdf, and student_t_lccdf.

This PR was Claude assisted

EDIT 2026-10-08: Added three more commits to fix more issues in functions touched by this PR

EDIT 2026-10-22: Rebased onto the current head of #3409, which now includes the reverse-mode and one-sided grad_reg_inc_beta from @SteveBronder. The seven commits keep their content; the shape gradients now request only the derivatives that each function needs. See the new Speed section.

Summary

Seven tail functions form the complement of the regularized incomplete beta function as 1 - inc_beta(a, b, x) in linear space: beta_lccdf, beta_proportion_lccdf, neg_binomial_lccdf, binomial_lccdf, and the q >= 2 branch of student_t_cdf, student_t_lcdf and student_t_lccdf. Once inc_beta rounds to 1, the complement is exactly 0, its log is -inf, and the gradients are inf or NaN. This happens for every tail probability below eps, and the lccdf functions are used for such tails. The q < 2 branch of the three Student t functions has a second subtraction of the same kind, 1 - r.

This PR evaluates each complement without the subtraction. A new internal function, internal::inc_beta_complement, evaluates it from whichever of x and 1 - x is at most 1/2. binomial_lccdf needs no helper and uses the symmetry relation directly. The shape gradients come from grad_reg_inc_beta on the same argument as the value. The density factors, whose powers underflow separately and gave 0 / 0, are replaced by inc_beta_ddz; in neg_binomial_lccdf the density is also evaluated at the argument at most 1/2. In the Student t q < 2 branch, 1 - r is formed without subtraction.

Which part of the PR fixes which issue:

issue	report	part of this PR
#3377	beta_lccdf(0.865169 | 1, 19) is -inf	the complement in beta_lccdf
#3234	student_t_cdf p-values of exactly 0	the complement in the Student t q >= 2 branch
#2923	student_t_cdf is 0 or 1 for large |t|	1 - r in the Student t q < 2 branch
#2031	neg_binomial_lccdf is -inf for large beta	the complement from 1 / (beta + 1), and the density factor of the beta gradient
#2790, partly	the binomial (lc)cdf underflows	tails correct down to DBL_MIN; below that log_inc_beta is needed

Depends on #3409: the shape gradients go through the repaired grad_reg_inc_beta. The value changes do not depend on it. This PR is based on the stable-inc-beta branch of #3409, so its diff shows only its own seven commits. I will retarget it to develop when #3409 merges.

Fixes #3377. On develop, the reporter's case beta_lccdf(0.865169 | 1, 19) returns -inf with infinite gradients. On this branch it returns -38.0709295957162, which equals the closed form 19 * log1m(y) to all printed digits. The three gradients match a 50-digit reference: d/dy = -140.917148134 = -19 / (1 - y), d/dalpha = 3.4106443951, d/dbeta = -2.00373313662. The cause is the defect described below: (1 - y)^19 = 3e-17 is below eps, so 1 - inc_beta(...) is exactly 0.

Fixes #3234. student_t_cdf returned exactly 0 for p-values below about 1e-16 whenever q = df / t^2 >= 2, the branch that formed 1 - inc_beta(...). At the six rows of the report, this branch matches a 40-digit reference to all printed digits. For example, 2 * student_t_cdf(-9.466995 | 196.9309, 0, 1) is 9.173253e-18, where develop gives 0.

Fixes #2923. In the q < 2 branch, r = 1 / (1 + q) rounds to 1 once q is below eps, so 1 - r is 0 and the tail probability is 0. The reported case student_t_cdf(-1e8 | 0.1, 0, 1) is 0 on develop for a true 0.0661503218, and its mirror at 1e8 is 1 for a true 0.9338496782. Before that point, 1 - r loses eps / q relative: 1e-13 at |t| = 10 with nu = 0.1. Every nu fails once |t| exceeds about 6.7e7 sqrt(nu); the Cauchy student_t_lcdf(-1e10 | 1, 0, 1) is -inf for a true -24.17. Over the reported table (nu = 0.1, |t| from 1 to 1e29, both signs), this branch is within 4e-16 of a 60-digit reference in the value and within 5e-15 in the y, nu and sigma gradients.

Fixes #2031. The reported case neg_binomial_lccdf(0 | 1, 1e18) is -inf on develop for a true -41.4465316739, which is log1m_exp(-alpha * log1p(1 / beta)). On this branch the value and both gradients are correct to 1e-15. For n >= 1, the beta gradient also needs the density factor at 1 / (beta + 1). At p = beta / (beta + 1) that factor is exactly 0 once p rounds to 1 (from beta near 1e16), where the true gradient is about -(n + 1) / beta; at beta = 1e8 it loses 8 digits.

call	develop	true value	this PR
beta_lccdf(0.9 | 50, 50)	-inf, gradients ±inf	-54.0888847	correct to 1e-10
beta_lccdf(0.758 | 352, 590)	-inf, gradients ±inf	-316.1955742	correct to 1e-10
beta_proportion_lccdf(0.758 | 0.4, 900)	-inf, gradients ±inf	-264.2052063	correct to 1e-7
neg_binomial_lccdf(50 | 2, 99)	-inf, gradients NaN	-230.9222919	correct to 1e-9
neg_binomial_lccdf(1 | 2, 1e18)	-inf, gradients inf	-81.7944510591	correct to 2e-15
binomial_lccdf(500 | 1000, 0.2)	-inf, gradient NaN	-227.9265046	correct to 1e-8
student_t_lccdf(10 | 1000, 0, 1)	-inf, gradients ∓inf	-50.8389515	correct to 1e-9
student_t_lcdf(-10 | 1000, 0, 1)	-inf, gradients ∓inf	-50.8389515	correct to 1e-9
student_t_cdf(-10 | 1000, 0, 1)	value correct, gradients NaN	8.34e-23	gradients finite and correct
student_t_cdf(-1e8 | 0.1, 0, 1)	0, gradients NaN	0.0661503218	correct to 1e-15
student_t_lcdf(-1e10 | 1, 0, 1)	-inf, gradients inf	-24.1705808158	correct to 1e-15

The defect on develop

This section describes develop. The reference is mpmath at 50 or 60 digits, with each value computed by two independent routes. For the beta family, the reference computes the complement on the small side of the symmetry relation. For Student t, the tail is 0.5 I_{nu / (nu + t^2)}(nu / 2, 1 / 2), differentiated numerically. I confirmed the numbers by compiling and running the real headers.

1. The complement cancels. Pn = 1.0 - inc_beta(alpha, beta, y) is 0 whenever inc_beta returns exactly 1. That happens whenever the true complement is below eps = 2.2e-16, that is, whenever the true lccdf is below log(eps) = -36.7. beta_lccdf(0.9 | 50, 50), with both shapes 50 and y at 0.9, returns -inf for a true -54.1. A Student t with 1000 degrees of freedom at 10 standard deviations returns -inf for a true -50.8. None of the parameter sets in the table above is unusual; the only condition is a tail probability below 1e-16.

2. The gradients follow. inv_Pn is inf, so the shape and y partials are inf where the density factor is finite, and NaN where it is not. student_t_cdf(-10 | 1000, 0, 1) is the worst case: its value 8.3e-23 is representable and correct, but both gradients are NaN.

3. The density factor underflows on its own. In binomial_lccdf, pow(theta, n) * pow(1 - theta, N - n - 1) / beta(N - n, n + 1) is a ratio of quantities that each underflow long before the ratio does; the other files have the equivalent expression. After the repair of the complement, binomial_lccdf(500 | 1000, 0.2) has a correct value but still a NaN gradient, because pow(0.2, 500) is 0 and beta(500, 501) is 0. In neg_binomial_lccdf, the factor pow(1 - p, n) is formed from p = beta / (beta + 1), which rounds to 1 for large beta.

4. 1 - r in the Student t q < 2 branch. The branch forms 1 - r by subtraction, with r = 1 / (1 + q) and q = nu / t^2. See #2923 above.

Where inc_beta underflows to 0, the complement is 1 and the old code is right. For example, neg_binomial_lccdf(3 | 900, 0.02) returns 0 with zero gradients on both trees, and that is the correct answer.

What this PR does

The complement from the argument at most 1/2. Boost's ibeta forms 1 - x from its argument, so an argument close to 1 loses relative precision in 1 - x. The symmetry relation 1 - I_x(a, b) = I_{1-x}(b, a) removes the cancellation of a tail below eps, but it passes 1 - x, which is close to 1 where x is small. The new internal::inc_beta_complement(a, b, x, one_m_x) in prim/fun/inc_beta.hpp takes both x and 1 - x, and calls Boost with the one that is at most 1/2:

• x > 1/2: inc_beta(b, a, 1 - x), by the symmetry relation;
• x below the mean a / (a + b), while inc_beta(a, b, x) is at most 1/2: 1 - inc_beta(a, b, x), which then has full precision;
• otherwise: Boost's complement ibetac(a, b, x).

The middle case is necessary because, for the arcsine case a = b = 1/2 (the Cauchy, nu = 1), Boost's ibetac evaluates asin(sqrt(1 - x)), which loses digits for small x. For an autodiff T_partials_return, which occurs only in second-order autodiff, the function uses the symmetry relation through inc_beta.

The shape gradients come from grad_reg_inc_beta on the same argument. For x <= 1/2, the call is on (a, b, x), as on develop, and its outputs are the negatives of the derivatives of the complement. For x > 1/2, the call is on (b, a, 1 - x), and its two outputs are the derivatives of the complement itself with respect to b and a. Above the mean, grad_reg_inc_beta applies the symmetry relation internally, so on an argument at most 1/2 it keeps its relative precision in the deep tail too.

• beta_lccdf: x = y, 1 - x = 1 - y.
• beta_proportion_lccdf: the same with a = mu kappa and b = kappa - mu kappa; the chain rule to mu and kappa does not change.
• neg_binomial_lccdf: x = p = beta / (beta + 1) and 1 - x = 1 / (beta + 1), which the function already has as inv_beta_p1; both are formed without cancellation.
• binomial_lccdf: the complement is inc_beta(n + 1, N - n, theta). Boost forms 1 - theta exactly from theta, so the function needs no helper. The derivative is closed form.
• neg_binomial_2_lccdf does not change: it returns neg_binomial_lccdf(n, phi, phi / mu), so it gets every change of neg_binomial_lccdf.
• student_t_cdf, student_t_lcdf, student_t_lccdf: the q >= 2 branch takes z = 1 - I_r(1/2, nu/2) from the helper, with r <= 1/3; its gradient code is the same as on develop. The q < 2 branch forms 1 - r as 1 / (1 + t^2 / nu), which has no cancellation and is 1 at t = 0. The two branches compute the same quantity, and a test checks continuity across q = 2.

The density factor. beta_lccdf, beta_proportion_lccdf, neg_binomial_lccdf and binomial_lccdf use inc_beta_ddz for the density, not the product of powers over a beta function. It is the same quantity; for double it is Boost's ibeta_derivative, which scales internally. beta_cdf already used it, so this also removes a duplicated expression. ibeta_derivative also forms 1 - x from its argument, and it returns 0 at x == 1 when its second shape is above 1. For p > 1/2, neg_binomial_lccdf therefore evaluates the density at p as the equal reflected density at 1 - p. The three Student t files keep their d_ibeta, which does not underflow in their parameter range.

The bound y = 1. beta_lccdf and beta_proportion_lccdf return negative_infinity() at y >= 1 before the loop. neg_binomial_lccdf and binomial_lccdf already have this explicit bound handling (see Testing).

After the change, the -inf threshold moves from eps to DBL_MIN, that is, from a tail probability of e^-36.7 to e^-709.

Accuracy after the change

Every value and every gradient in the table above matches the reference to all printed digits. The remaining relative errors are 1e-15 to 1e-10, which is the accuracy of Boost's inc_beta and of the gradient root; this PR adds no error of its own.

Where develop was correct, the change keeps its accuracy. I compared 117 points point by point, in the bulk (|t| down to 1e-7, y down to 1e-10, beta down to 1e-9) and in the tails. At every point where develop is finite, no value or gradient is more than 4 times less accurate than on develop above a relative error of 1e-14. For example, the nu gradient of student_t_cdf(-1e-7 | 1, 0, 1) is within 1e-15, and the alpha gradient of beta_lccdf(1e-10 | 0.5, 3) is within 6e-16.

Continuity across the Student t branch seam: at y = 1.9999, 2.0, 2.0001 with nu = 8, the value and the nu gradient are linear in y to 1e-8 on both sides of q = 2.

Speed

AMD EPYC 9965, g++ 11.5, -O3, one core, two rounds within 1 %. N = 1000 observations with shared shapes, value plus all gradients, ns per observation. The fvar<var> rows add the tangent and its reverse sweep, as in hessian(). Most of the gain comes from the gradient roots of #3409, which this PR is based on; the last column isolates the effect of the reverse-mode and one-sided roots from @SteveBronder on this PR's functions.

function	shapes	develop	this PR	speed-up	before the rebase / after
beta_lccdf	5, 5	12 410	212	59x	1.17
beta_lccdf, alpha only	5, 5	12 389	157	79x	1.58
beta_lccdf	50, 50	10 205	536	19x	1.11
beta_proportion_lccdf	0.5, 10	12 415	216	57x	1.16
neg_binomial_lccdf	5, 0.3	3 734	1 221	3.1x	1.12
neg_binomial_2_lccdf	16, 5	3 832	1 210	3.2x	1.12
student_t_lccdf	nu 3	2 354	330	7.1x	1.10
student_t_lccdf	nu 30	3 189	751	4.2x	1.23
student_t_cdf	nu 30	3 189	748	4.3x	1.23
beta_lccdf, fvar<var>	50, 50	70 123	1 437	49x	9.7
neg_binomial_lccdf, fvar<var>	5, 0.3	26 113	2 889	9.0x	3.3
student_t_lccdf, fvar<var>	nu 3	15 238	3 000	5.1x	0.98

binomial_lccdf does not call the gradient root and is not in the table: on this grid develop returns -inf at many points, so its time is no baseline.

Testing

Two new files have fixed references from the high-precision computation. Finite differences cannot see these defects: at the failing points the finite difference of -inf is NaN, and expect_ad would report agreement.

• rev/prob/beta_ccdf_log_test.cpp, 8 cases: beta_lccdf at three points below eps, one control, and four small-y points (two in the bulk, two with the complement below eps); beta_proportion_lccdf at one deep-tail point and two small-y points; neg_binomial_lccdf at two deep-tail points, six large-beta points (Numerical issues in neg_binomial_lccdf if beta is very large #2031: four with beta from 1e8 to 1e18, one control on each side of beta = 1) and four points with beta < 1; binomial_lccdf with its gradient.
• rev/prob/student_t_ccdf_log_test.cpp, 6 cases: student_t_lccdf at five points in both branches; four student_t_lccdf points at large t (Student T CDF precision for large values #2923); the reported Student T CDF precision for large values #2923 case, its mirror and the Cauchy lcdf at -1e10, with the value and all three gradients; three points near t = 0 with the value and all three gradients (nu = 1, nu = 30); the branch-seam continuity check; student_t_lcdf and student_t_cdf at the point where the gradients of the cdf were NaN.

Regression swap: with the develop versions of the changed headers installed under these two files, 11 of the 14 cases fail. Three pass: the binomial_lccdf case (its header was not swapped in this run), the seam-continuity check, and the near-t = 0 Student t case, where develop is accurate. With the headers of this branch, all 14 pass, before and after the swap.

Full suites, all passing on the rebased head: the root inc_beta unit tests in prim, rev, fwd and mix; the unit prob tests of beta, beta_proportion, binomial, student_t and neg_binomial* (which includes neg_binomial_2); and the test/prob generated cdf, lcdf and lccdf tests of beta, beta_proportion, student_t, neg_binomial, neg_binomial_2 and binomial (185 programs). Tested with g++ 11.5 on the #3409 head 3fd3fe2f6a with this PR's commits applied.

The generated tests found one defect that the fixed-reference tests did not. For double, inc_beta_ddz is Boost's ibeta_derivative, which throws on overflow where the old pow product returned inf. At y = 1 with a second shape below 1 the density is infinite, and the generated UpperBound case with an autodiff y reaches that point. This is the reason for the explicit bound at y = 1. The value at the bound does not change, -inf; the y partial there is 0 instead of inf or NaN.

Known limits, not addressed here

• Below DBL_MIN the complement is not representable. neg_binomial_lccdf(200 | 2, 99) and binomial_lccdf(500 | 1000, 0.01) are still -inf, for true values -920.34 and -1622.73. A repair needs log_inc_beta, the incomplete beta function in log space, which the tree does not have. The same function would repair the mirror defect in beta_lcdf and the related functions, where inc_beta underflows to 0 for y far below the mean. That is a separate PR: a new function with its own overloads and tests.
• In the Student t functions, t^2 overflows above |t| = 1.3e154; that range would need a log form.
• An lccdf value that is the log of a probability that rounds to 1 is returned as 0, as on develop: beta_proportion_lccdf(1e-8 | 0.3, 10) is 0 for a true -8.4e-23. The gradients there are correct.
• beta_cdf, neg_binomial_cdf and neg_binomial_2_cdf are not affected: they return the probability itself, and there 0 is the correct rounding.

Side Effects

No

Release notes

More stable and faster beta_lccdf, beta_proportion_lccdf, binomial_lccdf, neg_binomial_lccdf, neg_binomial_2_lccdf, student_t_cdf, student_t_lcdf, and student_t_lccdf.

Checklist

• Copyright holder: Aalto University

  The copyright holder is typically you or your assignee, such as a university or company. By submitting this pull request, the copyright holder is agreeing to the license the submitted work under the following licenses:
  - Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
  - Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)

• the basic tests are passing

  • unit tests pass (to run, use: ./runTests.py test/unit)
  • header checks pass, (make test-headers)
  • dependencies checks pass, (make test-math-dependencies)
  • docs build, (make doxygen)
  • code passes the built in C++ standards checks (make cpplint)

• the code is written in idiomatic C++ and changes are documented in the doxygen

• the new changes are tested

[jaburgoyne]

It’s hard to know the right place to do this, but thank you, @avehtari . I was on vacation and came back to see that almost every issue I’ve submitted to Stan is suddenly fixed with one or another of your PRs. I thought I’d try some of these easier ones to learn the system better before diving into my modified-Bessel wish…but now I guess modified Bessel is all that’s left!

[avehtari]

@jaburgoyne thanks for creating those issues, they did help to prioritize the work. I could do that modified Bessel, too, but happy to leave it also for you so you can learn more about Stan