Evaluate the incomplete beta complement without cancellation in the lccdf and Student t cdf functions - #3433

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

Summary

Seven tail functions form the complement of the regularized incomplete beta 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. That happens as soon as the tail probability is below eps, which is the regime the lccdf functions exist to represent. 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, from whichever of x and 1 - x is at most 1/2, through a new internal function internal::inc_beta_complement (in binomial_lccdf directly by the symmetry relation, which there needs no helper), and takes the shape gradients from grad_reg_inc_beta on the same argument. It replaces the density factors, whose powers underflow separately and gave 0 / 0, with inc_beta_ddz; in neg_binomial_lccdf that density is also evaluated at the argument at most 1/2. It forms 1 - r in the Student t q < 2 branch without subtraction.

Which part of the PR fixes which issue:

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

Fixes #3377. The reporter's case, beta_lccdf(0.865169 | 1, 19), returns -inf with infinite gradients on develop; on this branch it returns -38.0709295957162, equal to the closed form 19 * log1m(y) to all printed digits, and 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 one 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 that 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, 1 - r loses eps / q relative: 1e-13 at |t| = 10 with nu = 0.1. Any 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 in the value and 5e-15 in the y, nu and sigma gradients of a 60-digit reference.

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) it is exactly 0 once p rounds to 1 (from beta near 1e16), where the true gradient is about -(n + 1) / beta, and it loses 8 digits at beta = 1e8.

The defect on develop

Everything in this section describes develop. The reference is mpmath at 50 or 60 digits, each value by two independent routes; for the beta family the complement is computed on the small side of the symmetry relation, and for Student t the tail is 0.5 I_{nu / (nu + t^2)}(nu / 2, 1 / 2), differentiated numerically. The numbers were confirmed 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, which is 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), 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, so nothing looks wrong, while both gradients are NaN.

3. The density factor underflows on its own. pow(theta, n) * pow(1 - theta, N - n - 1) / beta(N - n, n + 1) in binomial_lccdf, and the equivalent in the other files, is a ratio of quantities that each underflow long before the ratio does. binomial_lccdf(500 | 1000, 0.2) has a correct value after the complement is repaired, and still a NaN gradient: 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.

Two regimes must be kept apart. Where inc_beta underflows to 0 the complement is 1 and the old code is right: 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 needed because Boost's ibetac for the arcsine case a = b = 1/2 (the Cauchy, nu = 1) evaluates asin(sqrt(1 - x)), which loses digits for small x. For 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 on (a, b, x), as on develop, whose outputs are the negatives of the derivatives of the complement; for x > 1/2 the call on (b, a, 1 - x), whose two outputs are the derivatives of the complement itself with respect to b and a. grad_reg_inc_beta applies the symmetry relation internally above the mean, 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 is unchanged.
• 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 is not changed itself: 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; continuity across q = 2 is tested.

The density factor. beta_lccdf, beta_proportion_lccdf, neg_binomial_lccdf and binomial_lccdf use inc_beta_ddz for the density instead of 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 returns 0 at x == 1 when its second shape is above 1; in neg_binomial_lccdf the density at p is therefore evaluated as the equal reflected density at 1 - p when p > 1/2. 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, the same explicit bound handling that neg_binomial_lccdf and binomial_lccdf already have (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 residual relative errors are 1e-15 to 1e-10, which is the accuracy of Boost's inc_beta and of the gradient root; nothing in this PR adds error of its own.

The change keeps the accuracy of develop where develop was correct. At 117 points, in the bulk (|t| down to 1e-7, y down to 1e-10, beta down to 1e-9) and in the tails, compared point by point, no value or gradient is more than 4 times less accurate than on develop above a relative error of 1e-14, at every point where develop is finite. 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) 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.

Testing

Two new files, with 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 and one control, and at 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 cdf's gradients were NaN.

Regression swap: with develop's versions of the changed headers installed under these two files, 11 of the 14 cases fail. The three that pass are 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 this branch's headers all 14 pass, before and after the swap.

Full suites, all passing on the final head: the root inc_beta unit tests in prim, rev, fwd and mix; the unit prob tests of beta, beta_proportion, student_t and neg_binomial* (which includes neg_binomial_2); and the test/prob generated tests of beta, beta_proportion, student_t, neg_binomial, neg_binomial_2 and binomial.

The generated tests found one thing the fixed-reference tests did not. inc_beta_ddz for double 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 it. This is the reason for the explicit bound at y = 1. The value at the bound is unchanged, -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 at all. neg_binomial_lccdf(200 | 2, 99) and binomial_lccdf(500 | 1000, 0.01) are still -inf, for true values -920.34 and -1622.73. That needs log_inc_beta, the incomplete beta in log space, which the tree does not have. The same function would repair the mirror defect in beta_lcdf and its siblings, where inc_beta underflows to 0 for y far below the mean. It 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, where 0 is the correct rounding.

Side Effects

No

Release notes

More stable computation in beta_lccdf, beta_proportion_lccdf, binomial_lccdf, neg_binomial_lccdf, neg_binomial_2_lccdf, student_t_cdf, student_t_lcdf, 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