More stable incomplete beta gradient roots by avehtari · Pull Request #3409 · stan-dev/math

avehtari · 2026-09-21T07:10:10Z

EDIT: 2026-10-29 updated after revision

Part 1 of 3. This PR holds the prim grad_reg_inc_beta (continued fraction and GRAD_RIB API), inc_beta_dda, inc_beta_ddb, inc_beta_ddz, fwd inc_beta, and their tests. #3489 (part 2) adds the reverse-mode grad_reg_inc_beta, rev inc_beta and their rev and mix tests. #3490 (part 3, stacked on part 2) moves the downstream cdfs to the new API; there beta_cdf calls the root once instead of twice. The pre-split head is on the branch stable-inc-beta-presplit. This text describes part 1 only.

Replace the incomplete beta gradient roots with a continued-fraction implementation from TensorFlow Probability.

Fixes #2417.

This PR was Claude assisted.

The scripts behind the numbers in this text are in https://github.com/avehtari/stan-pr-checks/tree/8b80a636364b7751bbf3553e63cf18d010f394c2/pr3409. Its README.md has one section per claim, with the commands and the expected output.

Summary

grad_reg_inc_beta, inc_beta_dda and inc_beta_ddb compute the derivatives of the regularized incomplete beta function I_z(a, b) with respect to the shapes a and b. This text calls such a function a gradient root, or a root. develop has two separate implementations of these derivatives. beta_cdf, neg_binomial_cdf and neg_binomial_2_cdf call path A (inc_beta_dda, inc_beta_ddb). Path B (grad_reg_inc_beta, through grad_2F1) has these callers: the lcdf and lccdf of beta, beta_proportion, neg_binomial and student_t; student_t_cdf; and the autodiff of inc_beta. At large shapes, at b ≫ a and near I = 1, both paths return NaN gradients or gradients wrong by O(1). Both paths also throw at some points.

This PR replaces both paths with one implementation: the Boik and Robinson-Cox (1998) algorithm as TensorFlow Probability implements it. It uses a power series near the endpoints and a modified-Lentz continued fraction elsewhere, and differentiates both term by term. It evaluates the prefactor in log space. inc_beta_dda and inc_beta_ddb become wrappers. The callers in prim/prob stay as they are, and so do the values of the distributions: the PR changes gradients only.

I compiled the merge base 579253050c ("develop" below) and this PR, and ran both on the same 18 000 points in nine regions:

	develop	this PR
throws, all 18 000 points	path A 98, path B 24	0
non-finite (NaN or throw), path A, all regions (1)	76 of 17 472	0
non-finite (NaN or throw), path B, all regions (1)	3104 of 17 472	0
non-finite, path B, large shapes near the mode	1325 of 2000	0
gradient of `log I` wrong by > 1e-2, path A, `b ≫ a`	408 of 2000	0
gradient of `log I` wrong by > 1e-2, path B, negbin `n > 200`	1203 of 1999	0
worst error of the gradient of `log I` (2)	n/a	1e-12
worst relative error where the derivative is not negligible (2)	n/a	3.8e-12
speed, `grad_reg_inc_beta` in `double`	1x	54x
speed, `beta_lcdf` per observation, `α = β = 5`	1x	56x
speed, `beta_cdf` per observation, `α = β = 50`	1x	2.1x

(1) The table leaves out the 528 points where the double I is 0.
(2) The table leaves out the 14 points where I is subnormal.

The defects in the current code on `develop`

This section describes develop at the merge base 579253050c, compiled. The reference is mpmath, computed by two independent routes: quadrature of the density at 80 digits, and numerical differentiation at 110 digits. At the 21 points of the unit tests the two routes agree to at least 60 digits.

1. `grad_reg_inc_beta` returns NaN for every parameter when `beta(a, b)` underflows

Its callers pass beta(a, b) as betaAB, and the function divides by it. beta(a, b) underflows when a is near b and a + b is above about 1100. The shape gradients of beta_lcdf(0.5 | 600, 600) are then NaN, and so are those of inc_beta(600, 600, 0.5) with autodiff shapes. In the sweep region with a, b in [100, 12500] and z near the mode, where I is of order 1, path B is NaN at 1325 of 2000 points.

The function also forms c3 = B(a, b) I and divides by B(a, b) at the end. When c3 underflows, both partials are 0. At neg_binomial_lcdf(132 | 313, 0.068), I = 1e-262 and B(a, b) I = 2e-381, and develop returns d/dα = 0 for a true -2.40.

2. `grad_reg_inc_beta` drops the hypergeometric term when `z^a (1-z)^b` underflows

When C = z^a (1-z)^b / a is 0, the function skips grad_2F1 and returns I (log z - 1/a - ψ(a) + ψ(a+b)). Near I = 1 the dropped term is the whole answer, because 2F1(a+b, 1; a+1; z) grows like (1-z)^(-b). beta_lcdf(0.758 | 352, 590) returns d/dα = 0.705 for a true -3.4e-138.

3. `inc_beta_dda` and `inc_beta_ddb` stop on an absolute threshold that the first term can meet

Both stop the series when fabs(summand) < 1e-10. The first summand contains the factor ((a+1)/(a+b))^3, which is below 1e-10 for b ≫ a, so the loop exits at once and returns the first term. beta_cdf(1.8e-4 | 1.145, 6786) returns d/dα = -0.203 for a true -0.386, at I = 0.65. In the sweep region with b ≫ a and small z, the d/dα part of the gradient of log I has an error above 1e-6 at 1301 of 2000 points, and above 1e-2 at 408.

4. `inc_beta_dda` and `inc_beta_ddb` reflect on fixed `z` thresholds

For z > 0.75, and in other fixed bands, they evaluate 1 - I_{1-z}(b, a) whether or not I is small. In the lower tail, the result is I' ≈ 1 times a bracket that must cancel to 1e-47. The computed bracket is rounding noise of size 1e-15. beta_cdf(0.76 | 446, 5) returns d/dα = -1.8e-15 for a true -1.06e-47, which makes the gradient of log(beta_cdf(...)) -4.4e+31.

At 98 of the 18 000 sweep points the loop reaches its iteration limit and throws "did not converge", inside the range that the docstring calls tested. Path B throws at 24 other points, from hypergeometric_2F1.

5. The existing unit tests assert wrong values

inc_beta_ddb_test.cpp passes digamma(a) as the fourth argument, digamma_b, and its expected constants come from that call (reported in #2417). inc_beta_dda_test.cpp and inc_beta_ddb_test.cpp also assert wrong output of develop: 1.8226241 for a true -3.94e-5, and 9.3959293 for a true 1.48e-8. grad_reg_inc_beta_test.cpp expects d/db = -inf at z = 1, where I ≡ 1 and the derivative is 0.

What this PR does

The new body of stan/math/prim/fun/grad_reg_inc_beta.hpp is a C++ port of _betainc_partials from TensorFlow Probability (tensorflow_probability/python/math/special.py, Apache-2.0; see the license note below). The algorithm is from Boik and Robinson-Cox (1998), Derivatives of the Incomplete Beta Function, JSS 3(1). It evaluates the DLMF 8.17.22 continued fraction with the modified Lentz method and carries the derivatives of the partial numerators through the recurrence. Near the endpoints it uses a power series, also differentiated term by term.

The repairs of the defects:

defects 1 and 2: the port evaluates the prefactor z^a (1-z)^b / (a B(a, b)) in log space;
defect 3: the series stops at eps / a and the continued fraction at 3 eps, both relative;
defect 4: the port applies the symmetry relation I_z(a, b) = 1 - I_{1-z}(b, a) for z ≥ (a + 1) / (a + b + 2) in the continued fraction, and for z at or above the mean a / (a + b) in the series.

The port differs from the TFP source in three places, and comments in the code mark each one.

Both stop tests also wait for the derivatives; TFP stops when the value has converged. For a positive integer b, the series factor (n - b) and the continued-fraction numerator d_{2m} are 0 at n = m = b. The value has then converged while the derivatives still need the rest of the series. In TFP 0.25.0, d/db is off by 4.2e-4 at (2, 3, 0.25) and by 6.7 % at (10, 1, 0.5), and d/da by 5.6 % at (1, 100, 0.01), where the symmetry relation makes a = 1 the integer shape. Integer shapes are common as initial values and as fixed shapes. I have not filed the TFP report yet.
The port uses the power series up to z = 0.5 (TFP: 0.95). Near z = 0.95 the series needs most of its 600-term limit, where the continued fraction needs tens.
The continued fraction switches to the symmetry relation at (a + 1) / (a + b + 2). TFP switches at the mode (a - 1) / (a + b - 2), which is above 1 for b < 1.

grad_reg_inc_beta<GRAD_RIB> computes only the requested partials (G1 for a, G2 for b, BOTH). inc_beta_dda uses G1 and inc_beta_ddb uses G2. Fwd inc_beta requests the partials of its autodiff shapes, and takes the z partial from inc_beta_ddz only when z is an autodiff type.

inc_beta_ddz(double, double, double) calls Boost's ibeta_derivative with a policy that returns infinity on overflow, and keeps all other defaults. Where the density is infinite (z = 1 with b < 1, z = 0 with a < 1) it returns +inf; develop throws there. I compared it with the default call on a grid of 183 600 points. At the 180 360 points where the default does not throw, the results are bit-identical. At the other 3 240 points it returns +inf.

Existing calls keep their form. The nine-argument grad_reg_inc_beta(g1, g2, ...) and the five-argument inc_beta_dda and inc_beta_ddb take the same arguments in the same order, and grad_reg_inc_beta keeps its betaAB argument without using it. The template signatures change, and the PR adds overloads:

new: grad_reg_inc_beta<GRAD_RIB>(a, b, z, digamma_a, digamma_b, digamma_ab, betaAB), which returns the requested partials;
new: container overloads of inc_beta_dda and inc_beta_ddz;
inc_beta_dda and inc_beta_ddb take five template types instead of one, and return auto;
inc_beta_ddz for double is an overload instead of an explicit specialization;
fwd inc_beta is one constrained template instead of seven overloads.

prim/fun.hpp now includes inc_beta_dda.hpp, inc_beta_ddb.hpp and inc_beta_ddz.hpp. On develop, prim.hpp reached them only through prob/beta_cdf.hpp, prob/neg_binomial_cdf.hpp and prob/neg_binomial_2_cdf.hpp.

Accuracy after the change

This section uses the 18 000 points and the mpmath reference of the summary. I compiled both trees with the same probe (g++ 11.4, -O2), which calls each root on its own, so a throw marks only that root. The metric is the absolute error of the gradient of log I. beta_lcdf returns this gradient, and autodiff of log(beta_cdf(...)) gives it too. A relative error says little for a derivative that is 1e-300 of I. A point counts once if either partial fails, and n is the number of points where the double I is not 0.

region	n	develop: non-finite / error > 1e-2	this PR: worst error
`a, b < 1`	2000	path A 0 / 0, path B 0 / 0	5e-13
`1 ≤ a, b ≤ 20`, `z < 0.75`	2000	path A 0 / 0, path B 0 / 0	3e-14
`z > 0.75`, `a < 500`	2000	path A 0 / 82, path B 3 / 125	1e-13
`b > a`, `b > 500`, `z ≤ 0.75`	1985	path A 1 / 87, path B 550 / 933	1e-12 (1)
`b ≫ a`, `z` small	2000	path A 0 / 408, path B 2 / 12	2e-13
large `a, b`, `z` uniform	1523	path A 46 / 36, path B 932 / 235	7e-13 (1)
large `a, b`, `z` near the mode	2000	path A 0 / 0, path B 1325 / 10	1e-13
negbin pattern, `n ≤ 200`	1965	path A 0 / 4, path B 0 / 26	5e-13 (1)
negbin pattern, `200 < n ≤ 5000`	1999	path A 29 / 22, path B 292 / 1203	4e-13

(1) The column leaves out the points where I is subnormal in double: 1, 11 and 2 points. There the error rises to 7e-3, because subnormal doubles have few significant bits.

Where the derivative is not negligible (|dI| / I ≥ 1e-3), the relative error of this PR is at most 3.8e-12, at large shapes. At the worst point of each region the two reference routes agree to 60 digits.

Scalar types, 2 007 points: with both shapes as autodiff arguments, the shape gradients of inc_beta in var, fvar<double> and fvar<fvar<double>> are bit-identical to the double root, on both trees. With one fvar shape, fwd inc_beta calls the one-sided root, which differs from the two-sided root by at most 2.5e-13 relative. With one var shape the gradient stays bit-identical, because rev inc_beta is unchanged in this PR.

Speed

Xeon Gold 6230 at 2.10 GHz, exclusive node, performance governor, g++ 14.3, -O3, Triton job 21271376; the two rounds agree within 1 %. Grid: a, b in [0.5, 60], z in [0.05, 0.95]. The benchmark loop adds a term that depends on the accumulator to every input, so the compiler cannot move the kernel out of the loop. I subtract the time of a baseline loop. The cycle counts of the baseline (5) and of 3 digamma + beta (330) show that the compiler removed no call.

Function level, net ns per call:

function	develop	this PR	ratio
`grad_reg_inc_beta` double, incl. 3 digamma + beta	36 860	684	54x
`grad_reg_inc_beta` var, value + grad	183 850	9 220	20x
`grad_reg_inc_beta` fvar	118 660	1 215	98x
`grad_reg_inc_beta` fvar<fvar>	335 390	3 035	110x
`inc_beta_dda` double, incl. 2 digamma	1 395	492	2.8x

Distribution level, N = 1000 observations, shared var shapes, value plus reverse sweep, ns per observation. The distribution code is the same on both trees; the gain comes from the new root.

function	shapes	develop	this PR	ratio
`beta_cdf`	5, 5	1 045	693	1.5x
`beta_lcdf`	5, 5	34 050	607	56x
`beta_lccdf`	5, 5	34 045	608	56x
`beta_cdf`	50, 50	3 566	1 659	2.1x
`beta_lcdf`	50, 50	27 720	1 279	22x
`beta_lccdf`	50, 50	27 723	1 282	22x
`neg_binomial_cdf`	5, 0.3	1 996	1 907	1.05x
`neg_binomial_lcdf`	5, 0.3	10 254	931	11x
`student_t_lcdf`	ν = 3	6 207	791	7.8x

Here beta_cdf still calls inc_beta_dda and inc_beta_ddb separately; #3490 changes it to one call. The cdf rows measure time only: the product of 1000 cdfs underflows, and the gradient is 0 on both trees.

Testing

The new tests compare with fixed references from the mpmath computation above. None of them uses finite differences, which do not detect these defects. The tests of the new overloads are in their own files, so that the defect tests compile against develop.

prim/fun/grad_reg_inc_beta_test.cpp: the two existing cases, with the z = 1 expectation corrected to 0 and the tolerance tightened to 1e-12. New cases for each defect class: beta(a, b) underflow, prefactor underflow, integer b, the tails, and a point where develop throws.
prim/fun/grad_reg_inc_beta_one_sided_test.cpp (new): the GRAD_RIB requests at the same points.
prim/fun/inc_beta_dda_test.cpp, prim/fun/inc_beta_ddb_test.cpp: rewritten from the reference at the twelve points of the old tests, with the digamma_b argument corrected, plus regression points and a point where develop throws.
prim/fun/inc_beta_dda_vectorized_test.cpp, prim/fun/inc_beta_ddz_vectorized_test.cpp (new): the container overloads.
prim/fun/inc_beta_ddz_test.cpp: +inf where the density is infinite.
rev/prob/beta_cdf_log_test.cpp (new): beta_lcdf shape gradients at (0.5 | 600, 600), at (0.758 | 352, 590) and at integer shapes.
rev/prob/neg_binomial_cdf_log_test.cpp: d/dα of neg_binomial_lcdf(132 | 313, 0.068).

Regression check (g++ 11.4): I replaced the six headers of this PR with their versions at 579253050c. Then 15 of the 19 tests in the six files that compile against develop fail. The four that pass, three older tests and the new grad_reg_inc_beta.moderate, test points where develop is right. The three files of the new overloads do not compile against develop. With the headers of this PR, all 25 tests pass.

On Triton (g++ 14.3) these also pass:

the unit suites of inc_beta and grad_2F1;
the mix tests of inc_beta and grad_reg_inc_beta;
the prob unit tests of beta, beta_proportion, neg_binomial, neg_binomial_2 and student_t;
the generated test/prob tests of beta, beta_proportion, binomial, neg_binomial, neg_binomial_2 and student_t (256 programs);
cpplint, test-math-dependencies, test-headers and doxygen.

License note

The new body of grad_reg_inc_beta.hpp derives from TensorFlow Probability, which is Apache-2.0. The file carries the TFP copyright notice and refers to licenses/tensorflow-probability-license.txt, which this PR adds, as the Boost-derived code in prim/fun/log_modified_bessel_first_kind.hpp does. The tree already contains Apache-2.0 code in vendored libraries (lib/tbb_2020.3, lib/benchmark_1.5.1). This file would be the first Apache-2.0-derived source file under stan/.

Known limits, not addressed here

beta_lccdf computes 1 - Pn in linear space after inc_beta. At (352, 590, 0.758) the double I is exactly 1, so the value is -inf and the gradients are infinite, for a true log(1 - I) ≈ -316.2. The defect is in the value code on both trees; Evaluate the incomplete beta complement without cancellation in the lccdf and Student t cdf functions #3433 repairs it.
The y partial of beta_lcdf(0.5 | 600, 600) and the β partial of neg_binomial_lcdf(132 | 313, 0.068) are NaN on both trees. Both come from the density terms of the two lcdfs. I will fix them in a separate PR.
Where I is subnormal, the error of the gradient of log I rises to 7e-3 (see the accuracy table).

Not included

grad_2F1, which other functions use, is unchanged.
The callers in prim/prob change in Move the beta, neg_binomial and student_t cdfs to the new grad_reg_inc_beta (part 3 of #3409) #3490.
The reverse-mode specialization is in Reverse-mode grad_reg_inc_beta (part 2 of #3409) #3489.

Release notes

Replaced the incomplete beta gradient roots with a continued-fraction implementation that does not return NaN or throw at 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
- 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 · 2026-09-21T07:11:31Z

I'll be away for the next two weeks, but made this PR to avoid duplicated work

avehtari · 2026-10-08T16:16:11Z

The improved grad_reg_inc_beta from this PR is used to fix also other distribution functions in #3433

avehtari · 2026-10-22T06:57:11Z

@SteveBronder, excellent additional speedup. My measurements

case	original PR/steve
first-order reverse mode, both shapes	1.02–1.10
one shape only (`beta_lcdf` alpha only, `neg_binomial_lcdf`)	1.24–1.38
`neg_binomial_2_cdf`	5.1
`fvar<var>` cdfs (one Hessian-vector product)	3.0–9.7
`fvar<var>` `inc_beta`	1.6

avehtari · 2026-10-22T08:48:26Z

@SteveBronder

The generated test/prob tests found one regression from your changes. beta_proportion_lcdf and beta_proportion_lccdf fail with fvar<fvar<var>> arguments: ibeta_derivative<long double>(...): Overflow Error. The cause is in inc_beta, not in these two files. The new rev and fwd inc_beta compute the z partial with inc_beta_ddz, which for double is Boost's ibeta_derivative. At z = 1 with b < 1 (and at z = 0 with a < 1) the density is infinite, and Boost throws. The previous code used the pow product there, which gave inf without a throw. Three direct reproductions on the head:

inc_beta(var 1.05, var 0.35, var 1.0), then grad(): throws.
beta_proportion_lcdf(1.0, fvar<var> mu = 0.75, 1.4), then d_.grad(): throws.
inc_beta(fvar<double>) at z = 1 with a zero tangent on z: throws.

I pushed a34e6f7. It calls ibeta_derivative with a policy that returns infinity on overflow, and keeps every other default. The result is bit-identical to the default call at 180 360 grid points, and is +inf at the 3 240 points where the default throws. It is the same call, so the speed does not change. At z = 1 the shape gradients of inc_beta are now 0 and the z gradient is +inf. New tests are in prim/fun/inc_beta_ddz_test, rev/fun/inc_beta_test and mix/fun/inc_beta_test; all three fail with the previous inc_beta_ddz.hpp.

In the generated tests, the UpperBound case leaves nodes at y = 1 on the stack, and a later grad() in the same process runs them. That is why the same binary passes when a filter skips the earlier cases.

The second commit, 3fd3fe2, starts the new mix test from an empty autodiff stack, because the expect_ad NaN cases earlier in that file leave nodes on the stack whose reverse pass throws.

With these two commits, the generated cdf, lcdf and lccdf tests of beta, beta_proportion, binomial, neg_binomial, neg_binomial_2 and student_t pass (185 programs, --run-all), and so do the unit suites of the incomplete beta functions. Tested on Roihu, AMD EPYC 9965, g++ 11.5.

One behavior stays as before: with a zero tangent on z, the fvar tangent of inc_beta at z = 1 is 0 * inf = NaN.

andrjohns · 2026-10-22T08:49:43Z

@avehtari can you split this one up as well? Given the amount of code changed, might be better to land this in parts. Thegrad_reg_inc_beta changes are the ones which need the most review/testing, while the distributional and other downstream function changes can land in separate PRs - much easier to review that way

avehtari · 2026-10-22T14:57:25Z

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

part	content	base	files	lines
B1 (#3409 reused)	prim `grad_reg_inc_beta` (continued fraction, `GRAD_RIB` API), `inc_beta_dda`, `inc_beta_ddb`, `inc_beta_ddz` (with `a34e6f7499`), fwd `inc_beta`, prim/fun.hpp, their prim tests	develop	10	+870 -311
B2	rev `grad_reg_inc_beta` (Hessian-vector product pass, SteveBronder's code), rev `inc_beta`, rev/fun.hpp, rev and mix tests (with `3fd3fe2f6a`)	B1	6	+975 -161
B3	the downstream cdfs on the new API: `beta_cdf`/`lcdf`/`lccdf`, `beta_proportion_lcdf`/`lccdf`, `neg_binomial_cdf`/`lcdf`/`lccdf`, `neg_binomial_2_cdf`, `student_t_cdf`/`lcdf`/`lccdf`, beta cdf tests	B1	14	+369 -528

B2 and B3 are independent of each other. B3 is is big in lines, but it's the same edit in each file.
In addition #3433 follows B3

Instead of 3 PRs, there could be 2 by having B1 and B2 in one PR. Or there could be 5 PRs by splitting B3 by family.

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

stan-buildbot · 2026-10-22T17:23:47Z

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
stat_comp_benchmarks/benchmarks/gp_regr/gen_gp_data.stan	0.04	0.04	0.98	-1.71% slower
stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan	0.17	0.16	1.04	3.52% faster
stat_comp_benchmarks/benchmarks/garch/garch.stan	0.88	0.86	1.02	2.24% faster
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.01	0.01	1.0	-0.21% slower
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.1	0.11	1.0	-0.25% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	8.23	8.0	1.03	2.77% faster
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	3.67	3.67	1.0	0.15% faster
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	32.97	32.58	1.01	1.16% faster
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.38	0.38	1.02	1.86% faster
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	6.78	6.86	0.99	-1.17% slower
stat_comp_benchmarks/benchmarks/sir/sir.stan	165.14	166.69	0.99	-0.94% slower
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.23	3.25	0.99	-0.71% slower
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.62	0.62	1.0	0.27% faster
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	25.76	25.72	1.0	0.16% faster
performance.compilation	391.31	386.77	1.01	1.16% faster
Mean result: 1.0057822753525654

Jenkins Console Log
Jenkins Build Stages
Commit hash: 8774185548440f572c1900f4a344e9ba089df121

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:                           52 bits physical, 57 bits virtual
CPU(s):                                  192
On-line CPU(s) list:                     0-191
Thread(s) per core:                      2
Core(s) per socket:                      48
Socket(s):                               2
NUMA node(s):                            2
Vendor ID:                               AuthenticAMD
CPU family:                              25
Model:                                   17
Model name:                              AMD EPYC 9474F 48-Core Processor
Stepping:                                1
Frequency boost:                         enabled
CPU MHz:                                 1497.134
CPU max MHz:                             4114.4229
CPU min MHz:                             1500.0000
BogoMIPS:                                7189.39
Virtualization:                          AMD-V
L1d cache:                               3 MiB
L1i cache:                               3 MiB
L2 cache:                                96 MiB
L3 cache:                                512 MiB
NUMA node0 CPU(s):                       0-47,96-143
NUMA node1 CPU(s):                       48-95,144-191
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:                  Not affected
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; Enhanced / Automatic IBRS; IBPB conditional; STIBP always-on; PBRSB-eIBRS Not affected; BHI Not affected
Vulnerability Srbds:                     Not affected
Vulnerability Tsa:                       Mitigation; Clear CPU buffers
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 amd_lbr_v2 nopl xtopology nonstop_tsc cpuid extd_apicid aperfmperf rapl pni pclmulqdq monitor ssse3 fma cx16 pcid 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 perfmon_v2 ibrs ibpb stibp ibrs_enhanced vmmcall fsgsbase bmi1 avx2 smep bmi2 erms invpcid cqm rdt_a avx512f avx512dq rdseed adx smap avx512ifma clflushopt clwb avx512cd sha_ni avx512bw avx512vl xsaveopt xsavec xgetbv1 xsaves cqm_llc cqm_occup_llc cqm_mbm_total cqm_mbm_local user_shstk avx512_bf16 clzero irperf xsaveerptr rdpru wbnoinvd amd_ppin cppc arat npt lbrv svm_lock nrip_save tsc_scale vmcb_clean flushbyasid decodeassists pausefilter pfthreshold avic v_vmsave_vmload vgif x2avic v_spec_ctrl vnmi avx512vbmi umip pku ospke avx512_vbmi2 gfni vaes vpclmulqdq avx512_vnni avx512_bitalg avx512_vpopcntdq la57 rdpid overflow_recov succor smca fsrm flush_l1d debug_swap

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

SteveBronder · 2026-10-22T19:08:16Z

@andrjohns sorry I made this PR a lot longer as I started fixing things adjacent / related to this PR. Some of this code has not been updated in years so I figured while we are here I would fix up nearby places

With my fixes as well a lot of cdfs are now 10x faster than they were before.

I've read over all of Aki's code and can handle the rest of the review if you are okay with that

stan-buildbot · 2026-10-23T04:46:32Z

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	8.31	7.71	1.08	7.23% faster
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	2.13	2.02	1.05	4.85% faster
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.06	0.05	1.06	5.51% faster
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	4.17	3.95	1.06	5.4% faster
stat_comp_benchmarks/benchmarks/garch/garch.stan	0.42	0.4	1.04	4.19% faster
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	2.65	2.51	1.05	5.14% faster
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.24	0.23	1.06	5.44% faster
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	19.0	17.77	1.07	6.48% faster
stat_comp_benchmarks/benchmarks/gp_regr/gen_gp_data.stan	0.02	0.02	1.04	3.4% faster
stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan	0.1	0.09	1.08	6.99% faster
stat_comp_benchmarks/benchmarks/sir/sir.stan	88.09	82.25	1.07	6.62% faster
stat_comp_benchmarks/benchmarks/arK/arK.stan	1.54	1.45	1.06	5.88% faster
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.01	0.01	1.02	2.29% faster
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.28	0.26	1.07	6.39% faster
performance.compilation	226.78	226.84	1.0	-0.03% slower
Mean result: 1.053612264726441

Jenkins Console Log
Jenkins Build Stages
Commit hash: f650e7d0b68a35569be5874169d8951cbf1b6ea6

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:                           46 bits physical, 48 bits virtual
CPU(s):                                  80
On-line CPU(s) list:                     0-79
Thread(s) per core:                      2
Core(s) per socket:                      20
Socket(s):                               2
NUMA node(s):                            2
Vendor ID:                               GenuineIntel
CPU family:                              6
Model:                                   85
Model name:                              Intel(R) Xeon(R) Gold 6248 CPU @ 2.50GHz
Stepping:                                7
CPU MHz:                                 1001.187
CPU max MHz:                             3900.0000
CPU min MHz:                             1000.0000
BogoMIPS:                                5000.00
Virtualization:                          VT-x
L1d cache:                               1.3 MiB
L1i cache:                               1.3 MiB
L2 cache:                                40 MiB
L3 cache:                                55 MiB
NUMA node0 CPU(s):                       0,2,4,6,8,10,12,14,16,18,20,22,24,26,28,30,32,34,36,38,40,42,44,46,48,50,52,54,56,58,60,62,64,66,68,70,72,74,76,78
NUMA node1 CPU(s):                       1,3,5,7,9,11,13,15,17,19,21,23,25,27,29,31,33,35,37,39,41,43,45,47,49,51,53,55,57,59,61,63,65,67,69,71,73,75,77,79
Vulnerability Gather data sampling:      Mitigation; Microcode
Vulnerability Indirect target selection: Mitigation; Aligned branch/return thunks
Vulnerability Itlb multihit:             KVM: Mitigation: Split huge pages
Vulnerability L1tf:                      Not affected
Vulnerability Mds:                       Not affected
Vulnerability Meltdown:                  Not affected
Vulnerability Mmio stale data:           Mitigation; Clear CPU buffers; SMT vulnerable
Vulnerability Old microcode:             Not affected
Vulnerability Reg file data sampling:    Not affected
Vulnerability Retbleed:                  Mitigation; Enhanced IBRS
Vulnerability Spec rstack overflow:      Not affected
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; Enhanced / Automatic IBRS; IBPB conditional; PBRSB-eIBRS SW sequence; BHI SW loop, KVM SW loop
Vulnerability Srbds:                     Not affected
Vulnerability Tsa:                       Not affected
Vulnerability Tsx async abort:           Mitigation; TSX disabled
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 dts acpi mmx fxsr sse sse2 ss ht tm pbe syscall nx pdpe1gb rdtscp lm constant_tsc art arch_perfmon pebs bts rep_good nopl xtopology nonstop_tsc cpuid aperfmperf pni pclmulqdq dtes64 monitor ds_cpl vmx smx est tm2 ssse3 sdbg fma cx16 xtpr pdcm pcid dca sse4_1 sse4_2 x2apic movbe popcnt tsc_deadline_timer aes xsave avx f16c rdrand lahf_lm abm 3dnowprefetch cpuid_fault epb cat_l3 cdp_l3 intel_ppin ssbd mba ibrs ibpb stibp ibrs_enhanced tpr_shadow flexpriority ept vpid ept_ad fsgsbase tsc_adjust bmi1 avx2 smep bmi2 erms invpcid cqm mpx rdt_a avx512f avx512dq rdseed adx smap clflushopt clwb intel_pt avx512cd avx512bw avx512vl xsaveopt xsavec xgetbv1 xsaves cqm_llc cqm_occup_llc cqm_mbm_total cqm_mbm_local dtherm ida arat pln pts hwp hwp_act_window hwp_epp hwp_pkg_req vnmi pku ospke avx512_vnni md_clear flush_l1d arch_capabilities

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

andrjohns · 2026-10-23T05:21:14Z

No worries both! I think it would be good to split at least the B1 row from Aki's table to a separate PR, since those are the big algorithmic changes

Port the Boik and Robinson-Cox continued fraction (as in TensorFlow Probability) to prim grad_reg_inc_beta, route inc_beta_dda, inc_beta_ddb and fwd inc_beta through it, and test against fixed mpmath references.

Moves this branch onto the split of #3409 (#3409, #3489, #3490). The conflicts are resolved to the tree of these commits replayed onto stable-inc-beta-cdfs.

stan-buildbot · 2026-10-27T08:46:49Z

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.01	0.6% faster
stat_comp_benchmarks/benchmarks/gp_regr/gen_gp_data.stan	0.06	0.06	0.98	-1.67% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	6.31	6.36	0.99	-0.79% slower
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	0.93	-7.35% slower
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	8.26	8.29	1.0	-0.41% slower
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	4.45	4.48	0.99	-0.65% slower
stat_comp_benchmarks/benchmarks/sir/sir.stan	217.86	218.38	1.0	-0.24% slower
stat_comp_benchmarks/benchmarks/garch/garch.stan	0.89	0.89	1.0	-0.09% slower
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.61	0.61	1.01	0.62% faster
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	42.87	42.39	1.01	1.11% faster
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.58	0.59	1.0	-0.19% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	20.97	20.76	1.01	1.02% faster
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.11	0.11	1.0	0.36% faster
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.18	3.17	1.0	0.29% faster
performance.compilation	363.83	364.57	1.0	-0.2% slower
Mean result: 0.9953205726430391

Jenkins Console Log
Jenkins Build Stages
Commit hash: 445e6af8e5d47439fea8d5d06ed11fb07f0bc4aa

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

stevieb-bot · 2026-10-28T00:08:43Z

Action: Not Ready for Human Review

Commit checked: 2764cf584d1c

Summary

Thank you for this work on the incomplete beta gradients and for the detailed write-up of the defects on develop.

The PR replaces the two implementations of the shape derivatives of the regularized incomplete beta function with one continued-fraction implementation based on TensorFlow Probability.
The description lists five defects on develop and says the PR fixes #2417.
It also reports large gains in accuracy and speed.
The description covers all three parts of the split, while this PR holds part 1 (the prim functions, fwd inc_beta, and the prim tests of these functions).

I read the changed code and tests at 2764cf5 and the replaced code on develop at the merge base 5792530.
In the container I ran the four prim test files at the head, and each of them against develop with the PR's test files copied in.
I also ran the generated test/prob tests of the six families named in a comment 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: 3 Confirmed, 9 Partially Confirmed, 9 Not Verified.
Description vs Code: 13 Confirmed, 5 Partially Confirmed, 1 Not Verified, 3 Contradicted.

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 PR fixes #2417

Quote:

Fixes #2417.

(PR description)

STATUS: Confirmed

Issue #2417 reports two things.
The first is that inc_beta_ddb_test.cpp passes digamma(a) where the function expects digamma(b).
The second is that inc_beta_ddb(15000, 1.25, 0.999) returns -2.1355524 while Wolfram and R finite differences give 2.0085e-6.
The new test passes digamma(b) in its helper, and it checks the point (15000, 1.25, 0.999) against 2.008497927784773e-6 with a relative tolerance of 1e-10.
That value agrees with the Wolfram and R values in the issue.
On develop with the PR's tests copied in, inc_beta_ddb returns -2.1355523032324224 at that point and the test fails.
At the head the test passes.
The comments on the issue mention other points.
(1.5, 1.25, 0.5) came up in a question about which incomplete beta function R computes, and the test checks it at L38.
(2, 1, 0.5) came from a check of the digamma argument, and (150, 1.25, 0.99) from a timing comparison.
At the grid points z = 0.25, a = 15, b = 12500 and z = 0.25, a = 150, b = 1250, the commenter judged the Wolfram values unstable and Stan's values close to correct.
None of these is reported as a defect of develop.

Evidence:

The PR's test inc_beta_ddb_test.cpp, which compiles on develop as well as at the head.
Code: inc_beta_ddb_test.cpp L15 (digamma(b)), inc_beta_ddb_test.cpp L46 (15000, 1.25, 0.999).

What I ran: python3 runTests.py -j8 test/unit/math/prim/fun/inc_beta_ddb_test.cpp on develop with the PR's tests copied in fails both tests, and at (15000, 1.25, 0.999) it prints got evaluates to -2.1355523032324224.
The same file passes at the head.
Output in the gist: inc_beta_ddb_test-base.out, prim_tests-head.out.

FACT: At z = 1 the shape gradients of inc_beta are 0 and the z gradient is +inf

Quote:

At z = 1 the shape gradients of inc_beta are now 0 and the z gradient is +inf.

(comment by @avehtari)

STATUS: Confirmed

At the head, grad_reg_inc_beta returns 0 for both partials when z is 0 or 1 (L435-L442).
Fwd inc_beta takes the shape partials from it and the z partial from inc_beta_ddz (L40-L57), which returns +inf at z = 1 with b < 1.
The tests grad_reg_inc_beta.1 (both partials 0 at z = 1) and inc_beta_ddz_infinite_density (+inf at (1.05, 0.35, 1.0)) pass at the head.

Evidence:

The tests grad_reg_inc_beta.1 and inc_beta_ddz_infinite_density.
Code: grad_reg_inc_beta_test.cpp L47-L63, inc_beta_ddz_test.cpp L46-L55.

What I read: The lines of grad_reg_inc_beta.hpp and fwd/fun/inc_beta.hpp linked in the finding.

What I ran: Both tests pass in the head run of the four prim test files.
Output in the gist: prim_tests-head.out.

FACT: With a zero tangent on z, the fvar tangent of inc_beta at z = 1 is NaN

Quote:

One behavior stays as before: with a zero tangent on z, the fvar tangent of inc_beta at z = 1 is 0 * inf = NaN.

(comment by @avehtari)

STATUS: Confirmed

At the head, fwd inc_beta adds x.d_ * inc_beta_ddz(a_val, b_val, x_val) whenever x is an fvar (L55-L57).
With x.d_ = 0 and an infinite density this is 0 * inf.
On develop the z term was x.d_ * d_x, with d_x computed from pow and also infinite there (L25-L27).

Evidence:

The fwd inc_beta code at the head and on develop.
Code: fwd/fun/inc_beta.hpp L55-L57.

What I read: The fwd inc_beta lines linked in the finding, at the head and on develop.

Partially Confirmed Facts

FACT: Defect 1: grad_reg_inc_beta returns NaN when beta(a, b) underflows

Quote:

1.
grad_reg_inc_beta returns NaN for every parameter when beta(a, b) underflows.
Its callers pass beta(a, b) as betaAB, and the function divides by it.

(PR description, The defects in the current code on develop)

STATUS: Partially Confirmed

The division follows from the code.
On develop, grad_reg_inc_beta divides both partials by betaAB (L56-L57).
The test grad_reg_inc_beta.beta_ab_underflow checks (600, 600, 0.5) at the head and passes.
The same file does not compile on develop, because its helper also calls the new grad_reg_inc_beta<GRAD_RIB::G1> overload.
The red run therefore shows that an API was added and does not show the NaN.
On develop, the tests of inc_beta_dda and inc_beta_ddb at (600, 600, 0.5) go through path A, which does not call grad_reg_inc_beta.
This PR does not test the examples beta_lcdf(y | 600, 600) and inc_beta with autodiff shapes.
The code behind the count of 1325 of 2000 points is missing.

Evidence:

The test grad_reg_inc_beta.beta_ab_underflow.
A test that fails on develop and the sweep that counted the NaN points are missing.
Code: grad_reg_inc_beta_test.cpp L77-L80.

What I read: On develop, grad_reg_inc_beta.hpp L40-L57 computes C and c3 = betaAB * inc_beta(a, b, z) and divides dBda and dBdb by betaAB.

What I ran: python3 runTests.py -j8 test/unit/math/prim/fun/grad_reg_inc_beta_test.cpp, run on develop with the PR's test files copied in, fails to compile with 'GRAD_RIB' has not been declared in 'stan::math'.
At the head all 7 tests in the file pass.
Output in the gist: grad_reg_inc_beta_test-base.out, prim_tests-head.out.

Request: Could you add a test of the defect that compiles on develop, for example the checks of the 9-argument grad_reg_inc_beta(g1, g2, ...) overload at (600, 600, 0.5) in a file without the GRAD_RIB calls?
A test of beta_lcdf at (600, 600) would also cover the example in the description, if it belongs in this part.
For the count of 1325 of 2000 points, please add the sweep script.

FACT: Defect 2: grad_reg_inc_beta drops the hypergeometric term when the prefactor underflows

Quote:

When C = z^a (1-z)^b / a is 0, the function skips grad_2F1 and returns I (log z - 1/a - ψ(a) + ψ(a+b)).
[...]
beta_lcdf(0.758 | 352, 590) returns d/dα = 0.705 for a true -3.4e-138.
neg_binomial_lcdf(132 | 313, 0.068) returns d/dα = 0 for a true -2.3.

(PR description, The defects in the current code on develop)

STATUS: Partially Confirmed

The skipped call follows from the code.
On develop, grad_reg_inc_beta calls grad_2F1 only when C is nonzero (L48-L51).
With the dF terms at 0, the partials reduce to the expression in the description.
The test grad_reg_inc_beta.prefactor_underflow checks (352, 590, 0.758) against -3.3772717797123242e-138 and 4.5390666074606649e-138 and passes at the head.
As in defect 1, the file does not compile on develop, so the run does not show the value 0.705 on develop.
This PR does not test the neg_binomial_lcdf(132 | 313, 0.068) example.

Evidence:

The test grad_reg_inc_beta.prefactor_underflow.
A test that fails on develop and a test of the neg_binomial_lcdf example are missing.
Code: grad_reg_inc_beta_test.cpp L82-L86.

What I read: On develop, grad_reg_inc_beta.hpp L43-L57 sets C = exp(a * c1 + b * c2) / a and skips grad_2F1 when it is 0.

What I ran: The same runs as defect 1.
The PR's test file fails to compile on develop and passes at the head.
Output in the gist: grad_reg_inc_beta_test-base.out, prim_tests-head.out.

Request: Could you add a test at (352, 590, 0.758) that compiles on develop, and a test of the neg_binomial_lcdf(132 | 313, 0.068) example, or say which of the later parts holds them?

FACT: Defect 3: inc_beta_dda and inc_beta_ddb stop on an absolute threshold that can be met before the first term

Quote:

Both stop the series when fabs(summand) < 1e-10.
The first summand contains the factor ((a+1)/(a+b))^3, which is below 1e-10 for b ≫ a.
The loop then exits at once and returns the first term.
beta_cdf(1.8e-4 | 1.145, 6786) returns d/dα = -0.203 for a true -0.386, at I = 0.65.

(PR description, The defects in the current code on develop)

STATUS: Partially Confirmed

The threshold follows from the code.
On develop, both functions loop while fabs(summand) > 1e-10, and the first summand includes the factor pow((a + 1) / (a + b), 3) (inc_beta_dda.hpp L64-L83).
At (1.145, 6786, 1.8e-4) the inc_beta_ddb test fails on develop, with -7.3310671133858816e-06 against the expected 5.8439109664263782e-05, and passes at the head.
The quoted example is the a derivative.
The matching inc_beta_dda test does not compile on develop, because the file also tests the new vectorized overload.
The run therefore does not show -0.203 on develop.
The code behind the counts of 1301 and 408 of 2000 points is missing.

Evidence:

The regression tests in inc_beta_dda_test.cpp and inc_beta_ddb_test.cpp.
An inc_beta_dda test that compiles on develop and the sweep script are missing.
Code: inc_beta_dda_test.cpp L62, inc_beta_ddb_test.cpp L62.

What I read: On develop, inc_beta_dda.hpp L64-L83 and inc_beta_ddb.hpp L59-L75 set threshold = 1e-10 and test it before the first iteration.

What I ran: The PR's inc_beta_ddb_test.cpp, run on develop with python3 runTests.py -j8 test/unit/math/prim/fun/inc_beta_ddb_test.cpp, fails at a=1.145 b=6786 z=0.00018.
The PR's inc_beta_dda_test.cpp fails to compile on develop with no matching function for call to 'inc_beta_dda(Eigen::VectorXd&, ...)'.
Both files pass at the head.
Output in the gist: inc_beta_ddb_test-base.out, inc_beta_dda_test-base.out, prim_tests-head.out.

Request: Could you move the vectorized inc_beta_dda test to a separate file?
The scalar regression tests would then compile on develop and could show the a derivative of -0.203.
For the counts of 1301 and 408 points, please add the sweep script.

FACT: Defect 4: inc_beta_dda and inc_beta_ddb reflect on fixed z thresholds

Quote:

For z > 0.75, and in other fixed bands, they evaluate 1 - I_{1-z}(b, a) wherever the mass is.
[...]
beta_cdf(0.76 | 446, 5) returns d/dα = 8.9e-16 for a true -1.06e-47, so the gradient of log(beta_cdf(...)) is 2.5e+31.
At some large-shape points the loop also reaches the k > 1e5 iteration guard and throws, inside the range that the docstring calls tested.

(PR description, The defects in the current code on develop)

STATUS: Partially Confirmed

The fixed thresholds follow from the code.
On develop, inc_beta_dda reflects for z > 0.75 && a < 500 and in other fixed bands (L44-L62), and inc_beta_ddb does the same (L47-L57). At (446, 5, 0.76) the inc_beta_ddb test fails on develop with 1.1102230246251565e-15 against the expected 1.2763067943072555e-46, which matches the rounding noise of size 1e-15 that the description gives. It passes at the head. The quoted example is the a derivative, and the inc_beta_dda test at that point does not compile on develop. The inc_beta_ddb tests do not throw on develop, and the other test files do not compile there. The runs therefore do not reach the k > 1e5 guard.

Evidence:

The regression tests at (446, 5, 0.76).
An inc_beta_dda test that compiles on develop and a test that reaches the iteration guard are missing.
Code: inc_beta_dda_test.cpp L64, inc_beta_ddb_test.cpp L63.

What I read: The reflection conditions on develop are at inc_beta_dda.hpp L44-L62, and the guard that throws for k > 1e5 is at L93-L96.

What I ran: The same runs as defect 3.
On develop the PR's inc_beta_ddb test fails at a=446 b=5 z=0.76, and the PR's inc_beta_dda file does not compile.
Output in the gist: inc_beta_ddb_test-base.out, inc_beta_dda_test-base.out, prim_tests-head.out.

Request: Could you make the inc_beta_dda test at (446, 5, 0.76) compile on develop, as in defect 3, and add a point where develop throws from the k > 1e5 guard?

FACT: Defect 5: the existing unit tests assert wrong values

Quote:

5.
The existing unit tests assert wrong values.
inc_beta_ddb_test.cpp passes digamma(a) as the fourth argument, which is digamma(b), and its expected constants come from that call (reported in #2417).
inc_beta_dda_test.cpp and inc_beta_ddb_test.cpp also assert the output of develop at points where it is wrong by orders of magnitude: 1.8226241 for a true -3.94e-5, and 9.3959293 for a true 1.48e-8.
grad_reg_inc_beta_test.cpp expects d/db = -inf at z = 1, where I ≡ 1 and the derivative is 0.

(PR description, The defects in the current code on develop)

STATUS: Partially Confirmed

The develop tests show each point.
inc_beta_ddb_test.cpp passes digamma(small_a) as the fourth argument (L18-L21).
It asserts 9.3959293 at (1.5, 12500, 0.001) (L42-L44), and inc_beta_dda_test.cpp asserts 1.8226241 at the same point (L42).
grad_reg_inc_beta_test.cpp expects -infinity for g2 at z = 1 (L19).
The new tests assert -3.937563719959171e-5 and 1.4782043391986094e-8 at that point and 0 at z = 1.
With digamma(b) passed, develop returns -0.00088048922293390845 at (1.5, 12500, 0.001) against the expected 1.4782043391986094e-8, and the head passes.
Two parts are not shown.
The true values -3.94e-5 and 1.48e-8 come from the mpmath constants, whose script is missing (see the fact check of the references below).
The PR's inc_beta_dda test does not compile on develop, so the runs do not show the develop value at the inc_beta_dda point.

Evidence:

The develop tests and the PR's rewritten tests.
Code: inc_beta_ddb_test.cpp L48, inc_beta_dda_test.cpp L47, grad_reg_inc_beta_test.cpp L47-L63.

What I read: I read the three test files on develop at 5792530 and the rewritten files at the head, at the lines linked in the finding.

What I ran: The PR's inc_beta_ddb test, run on develop, fails at a=1.5 b=12500 z=0.001.
The four prim test files pass at the head.
Output in the gist: inc_beta_ddb_test-base.out, prim_tests-head.out.

Request: The mpmath script requested below would back the true values.
Moving the vectorized inc_beta_dda test to a separate file, as requested for defect 3, would let the inc_beta_dda point run on develop.

FACT: The four scalar types give bit-identical values on both trees

Quote:

The four scalar types double, var, fvar<double> and fvar<fvar<double>> give bit-identical values on both trees.

(PR description, Accuracy after the change)

STATUS: Partially Confirmed

For var and fvar<double> with both shapes as autodiff arguments, the code calls the same double function as the double case, on both trees.
At the head, rev inc_beta calls the 9-argument grad_reg_inc_beta with double values, which forwards to grad_reg_inc_beta<GRAD_RIB::BOTH>, and fwd inc_beta calls grad_reg_inc_beta<GRAD_RIB::BOTH> on the double values.
On develop, rev and fwd inc_beta both call the 9-argument overload with double values.
For fvar<fvar<double>> the root runs in fvar<double> arithmetic, and at the head fwd inc_beta with one autodiff shape calls the G1 or G2 overload, whose stop test differs from BOTH.
These cases need a measurement, and this PR's tests do not compare the scalar types.

Evidence:

The rev and fwd inc_beta code.
Code that compares fvar<fvar<double>> and the one-shape cases with double is missing.
Code: rev/fun/inc_beta.hpp L26-L28, fwd/fun/inc_beta.hpp L36-L54.

What I read: At the head, rev/fun/inc_beta.hpp L26-L28, grad_reg_inc_beta.hpp L498-L499, and fwd/fun/inc_beta.hpp L36-L54.
On develop, fwd/fun/inc_beta.hpp L22-L24 and rev/fun/inc_beta.hpp L25-L28.

Request: Could you add or link the code that compares fvar<fvar<double>> and the one-shape cases with double, or say which test in #3489 covers it?

FACT: Regression check: the tests fail on develop and pass on the branch

Quote:

Regression check: I installed the develop headers and ran these six files.
All 6 suites fail, with 13 of 15 cases failing.
With the headers of this branch, all 15 cases pass, both before and after the swap.

(PR description, Testing)

STATUS: Partially Confirmed

Three of the six files named in the description are in this PR (grad_reg_inc_beta_test.cpp, inc_beta_dda_test.cpp, and inc_beta_ddb_test.cpp).
The three rev files are in #3489 and #3490.
At the head these three files have 12 tests (7, 3, and 2), and all of them pass.
The 3 tests of inc_beta_ddz_test.cpp, a file not among the six, also pass.
On develop, grad_reg_inc_beta_test.cpp and inc_beta_dda_test.cpp do not compile, because they call the new GRAD_RIB overload and the new vectorized inc_beta_dda.
Those two runs show that an API was added.
inc_beta_ddb_test.cpp compiles on develop, and both of its tests fail.
The count of 13 of 15 failing cases cannot be matched against this PR's files.

Evidence:

The PR's test files.
Three of the six files are in other PRs, and two of the three here do not compile on develop.
Code: grad_reg_inc_beta_test.cpp, inc_beta_dda_test.cpp, inc_beta_ddb_test.cpp.

What I ran: I ran each test file with python3 runTests.py -j8 <file> on develop with the PR's tests copied in, and all four files together at the head.
On develop three files, including inc_beta_ddz_test.cpp, fail to compile, and inc_beta_ddb_test fails 2 of 2 tests.
At the head all 15 tests of the four files pass.
Output in the gist: prim_tests-head.out, grad_reg_inc_beta_test-base.out, inc_beta_dda_test-base.out, inc_beta_ddb_test-base.out, inc_beta_ddz_test-base.out.

Request: Could you update the regression check for this part, giving the files and counts for this PR alone?
Splitting the tests of the new overloads into their own files would let the defect tests run on develop.

FACT: inc_beta_ddz returns infinity instead of throwing where the density is infinite

Quote:

I pushed a34e6f7.
It calls ibeta_derivative with a policy that returns infinity on overflow, and keeps every other default.
[...]
New tests are in prim/fun/inc_beta_ddz_test, rev/fun/inc_beta_test and mix/fun/inc_beta_test; all three fail with the previous inc_beta_ddz.hpp.

(comment by @avehtari)

STATUS: Partially Confirmed

The policy change follows from the code.
On develop, the double overload calls ibeta_derivative(a, b, z) with the default policy (L37-L41).
At the head it passes a policy with overflow_error<ignore_error> (L50-L56).
The test inc_beta_ddz_infinite_density passes at the head.
inc_beta_ddz_test.cpp does not compile on develop, because it also tests the new container overload.
The red run therefore does not show the throw.
The rev and mix tests that the comment names are in #3489.

Evidence:

The test inc_beta_ddz_infinite_density.
A version that compiles on develop is missing.
Code: inc_beta_ddz_test.cpp L46-L55.

What I read: inc_beta_ddz(double, double, double) at the head, inc_beta_ddz.hpp L50-L56.

What I ran: python3 runTests.py -j8 test/unit/math/prim/fun/inc_beta_ddz_test.cpp fails to compile on develop with no matching function for call to 'inc_beta_ddz(double, Eigen::VectorXd&, Eigen::VectorXd&)'.
At the head all 3 tests pass.
Output in the gist: inc_beta_ddz_test-base.out, prim_tests-head.out.

Request: Could you put the container tests of inc_beta_ddz in a separate file?
inc_beta_ddz_infinite_density would then run on develop and could show the throw.

FACT: The generated cdf, lcdf, and lccdf tests of six families pass

Quote:

With these two commits, the generated cdf, lcdf and lccdf tests of beta, beta_proportion, binomial, neg_binomial, neg_binomial_2 and student_t pass (185 programs, --run-all), and so do the unit suites of the incomplete beta functions.

(comment by @avehtari)

STATUS: Partially Confirmed

The comment reports a run on the branch before the split.
That branch also changed rev inc_beta and the cdf functions, which are now in #3489 and #3490, and the mix test that 3fd3fe2 changed is in #3489.
I ran the generated test/prob tests of the six families at this PR's head.
All 256 generated programs build and pass.
Of these, 185 are cdf, lcdf, or lccdf programs, which matches the count in the comment.
The four prim unit suites of the incomplete beta functions also pass at the head.
The rev and mix unit suites that the comment includes are in #3489, so they are not part of this run.

Evidence:

The generated test/prob tests in the repository and the PR's prim unit tests.
The rev and mix suites are in #3489.
Code: test/prob/beta, test/prob/student_t.

What I ran: python3 runTests.py -j8 --run-all test/prob/beta test/prob/beta_proportion test/prob/binomial test/prob/neg_binomial test/prob/neg_binomial_2 test/prob/student_t at the head.
With --run-all the exit status does not reflect test failures.
review.py run cut the middle of the output for size.
A second run read the gtest XML files that the first run wrote and reports failures="0" and errors="0" for all 256 programs.
Output in the gist: prob_tests-head.out, prob_tests-head-results.out, prim_tests-head.out.

Request: Could you note that the rev and mix suites in this sentence are part of #3489?

Not Verified Facts

FACT: The test references come from an 80-digit mpmath computation by two independent routes

Quote:

The reference is mpmath at 80 digits, computed by two independent routes: quadrature of the density, and numerical differentiation at 110 digits.
The two routes agree to 60 digits at 54 check points.
[...]
Six test files use fixed references from the 80-digit mpmath computation.

(PR description, The defects in the current code on develop and Testing)

STATUS: Not Verified

The expected values in grad_reg_inc_beta_test.cpp, inc_beta_dda_test.cpp, and inc_beta_ddb_test.cpp are fixed constants, and the comments in these three files say they come from mpmath at 80 digits.
The PR does not include or link the mpmath code that produced them, or the check that the two routes agree to 60 digits at 54 points.
As the description says, the tests do not use finite differences.

Evidence:

The constants are in the tests.
The code that produced them is missing.

Request: Could you add the mpmath script that produces the expected constants (both routes), for example under test/unit/math/prim/fun/ or in a linked gist?
A reviewer could then regenerate them.

FACT: Accuracy over 18 000 points in nine regions

Quote:

I compiled both trees and ran them on the same 18 000 points in nine regions (details below):

(PR description, Summary and Accuracy after the change)

STATUS: Not Verified

The summary table and the Accuracy section report counts of non-finite and inaccurate gradients on develop, and worst errors on the branch (1e-11 absolute, 3.7e-9 relative, 6e-13 relative outside b >> a), over a sweep of 18 000 points.
The PR does not include or link the sweep or the probe that compiled both trees.
The same holds for the statement that the relative error rises to 1e-2 where I is subnormal.

Evidence:

The sweep, the probe, and the reference values at the 18 000 points are missing.

Request: Could you add the sweep script and the probe program, with the commands that build them against develop and against this branch?

FACT: Function-level speedups of 2.2x to 68x

Quote:

Function level, net ns per call:

(PR description, Speed)

STATUS: Not Verified

The description gives net ns per call on develop and on the branch for grad_reg_inc_beta in four scalar types and for inc_beta_dda, and describes the benchmark loop.
The PR does not include or link the benchmark.

Evidence:

The benchmark code is missing.

Request: Could you add the benchmark, for example as a google-benchmark file, with the command and compiler flags used?

FACT: Distribution-level speedups of 1.1x to 40x

Quote:

Distribution level, N = 1000 observations, shared var shapes, value plus reverse sweep, ns per observation:

(PR description, Speed)

STATUS: Not Verified

The description gives ns per observation for the cdf, lcdf, and lccdf of beta, neg_binomial, and student_t.
The PR does not include or link the benchmark.
At this commit beta_cdf still calls inc_beta_dda and inc_beta_ddb separately, so its row describes the code in #3490.

Evidence:

The benchmark code is missing.

Request: Could you add the benchmark and say which rows apply to this PR alone?

FACT: Path B needs 607 607 terms at z = 0.998, and the continued fraction costs about 600 cycles

Quote:

At z = 0.998 it needs 607 607 terms.
The continued fraction of the branch costs about 600 cycles in addition to the digamma calls.

(PR description, Speed)

STATUS: Not Verified

The PR does not include code that counts the terms of the 2F1 series on develop or the cycles of the continued fraction.

Evidence:

Code that counts terms and cycles is missing.

Request: Could you add the code that produced the term count and the cycle count, or drop the two numbers?

FACT: TFP's own derivatives are off at integer shapes

Quote:

TFP's own d/db is off by 4.2e-4 at (2, 3, 0.25) and by 6.7 % at (10, 1, 0.5).
Its d/da is off by 5.6 % at (1, 100, 0.01), where the symmetry relation makes b' = 1.

(PR description, What this PR does)

STATUS: Not Verified

The PR does not include code that runs TensorFlow Probability.
The tests check this PR's port at (2, 3, 0.25), (10, 1, 0.5), and (1, 100, 0.01) and pass at the head, but they do not show the errors of TFP at those points.

Evidence:

Code that evaluates TFP at the three points is missing.

Request: Could you link the script that evaluated TFP at these points, or the upstream report once it is filed?

FACT: The new ibeta_derivative call is bit-identical to the default at 180 360 grid points

Quote:

The result is bit-identical to the default call at 180 360 grid points, and is +inf at the 3 240 points where the default throws.

(comment by @avehtari)

STATUS: Not Verified

The PR does not include or link the grid comparison of the two calls.

Evidence:

The comparison script is missing.

Request: Could you link the script that compared the two calls on the grid?

FACT: Speed ratios of the original PR over the version with the rev changes

Quote:

@SteveBronder, excellent additional speedup.
My measurements

(comment by @avehtari)

STATUS: Not Verified

The comment gives ratios between two versions of the branch for first-order reverse mode, one-shape cases, neg_binomial_2_cdf, and fvar<var>.
The PR does not include or link the benchmark.
The rev changes that the ratios compare are in #3489.

Evidence:

The benchmark code is missing.

Request: Could you link the benchmark behind this table, or move the table to #3489 where the rev code is?

FACT: Many cdfs are now 10x faster

Quote:

With my fixes as well a lot of cdfs are now 10x faster than they were before.

(comment by @SteveBronder)

STATUS: Not Verified

The comment includes a heatmap image.
The PR does not include or link the code that produced it, and the cdf changes are in #3490.

Evidence:

The benchmark behind the heatmap is missing.

Request: Could you link the benchmark that produced the heatmap?

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: This PR is part 1 of the split

Quote:

This PR is now part 1: the prim grad_reg_inc_beta (continued fraction, GRAD_RIB API), inc_beta_dda, inc_beta_ddb, inc_beta_ddz, fwd inc_beta, and their prim tests.

(PR description, opening note)

STATUS: Confirmed

The PR changes prim/fun/grad_reg_inc_beta.hpp (with the GRAD_RIB enum at L28), inc_beta_dda.hpp, inc_beta_ddb.hpp, inc_beta_ddz.hpp, fwd/fun/inc_beta.hpp, prim/fun.hpp, and the four prim test files of these functions.

FACT: The callers of path A and path B on develop

Quote:

Path A (inc_beta_dda, inc_beta_ddb) serves beta_cdf, neg_binomial_cdf and neg_binomial_2_cdf.
Path B (grad_reg_inc_beta, through grad_2F1) serves the lcdf and lccdf of beta, beta_proportion, neg_binomial and student_t, and also student_t_cdf and the autodiff of inc_beta.

(PR description, Summary)

STATUS: Confirmed

On develop, inc_beta_dda is called from beta_cdf (L124), neg_binomial_cdf (L98), and neg_binomial_2_cdf (L107), and inc_beta_ddb from beta_cdf (L130).
grad_reg_inc_beta calls grad_2F1 (L50).
It is called from the lcdf and lccdf of beta, beta_proportion, neg_binomial, and student_t, from student_t_cdf (L116), and from rev inc_beta (L26) and fwd inc_beta (L23).

FACT: Both paths go through one implementation, with inc_beta_dda and inc_beta_ddb as wrappers

Quote:

This PR replaces both paths with one implementation [...]
inc_beta_dda and inc_beta_ddb become wrappers over grad_reg_inc_beta that return one of the two derivatives.

(PR description, Summary and What this PR does)

STATUS: Confirmed

inc_beta_dda returns grad_reg_inc_beta<GRAD_RIB::G1>(...) (L37-L43), and inc_beta_ddb returns grad_reg_inc_beta<GRAD_RIB::G2>(...) (L33-L39).
The 9-argument overload that the path B callers use now calls the same templated function (L494-L502).

FACT: The new grad_reg_inc_beta evaluates a continued fraction with the modified Lentz method and a power series near the endpoints, both differentiated term by term

Quote:

It evaluates the DLMF 8.17.22 continued fraction with the modified Lentz method and carries the derivatives of the partial numerators through the recurrence.
Near the endpoints it uses the power series from Cephes, also differentiated term by term.

(PR description, What this PR does)

STATUS: Confirmed

The doc comment names Boik and Robinson-Cox (1998) as implemented in TensorFlow Probability (L362-L364).
The continued-fraction kernel evaluates DLMF 8.17.22 with the modified Lentz method and carries the logarithmic derivatives through each step (L147-L353).
The power-series kernel differentiates each term (L32-L145).
The code does not name Cephes or _betainc_partials.

FACT: The prefactor is evaluated in log space

Quote:

The port evaluates the prefactor z^a (1-z)^b / (a B(a, b)) in log space, which removes defects 1 and 2.

(PR description, What this PR does)

STATUS: Confirmed

The power series computes exp(a * log_x - lbeta(a, b)) (L128), and the continued fraction computes h * exp(a * log_x + b * log1m_x - log(a) - lbeta(a, b)) (L336-L337).
Neither divides by beta(a, b).

FACT: The power series switches to the symmetry relation at the mean

Quote:

and at the mean for the series, which removes defect 4.

(PR description, What this PR does)

STATUS: Confirmed

With mean = a / (a + b), the power series of I_{1-z}(b, a) is used when z >= mean, and the series of I_z(a, b) when z < mean (L446-L457). The series is used at all only when b z <= 1 and z <= 0.5, or a (1 - z) <= 1 and z >= 0.5.

FACT: The series stops at eps / a and the continued fraction at 3 eps

Quote:

It has no absolute threshold: the series stops at eps / a and the continued fraction at 3 eps, which removes defect 3.

(PR description, What this PR does)

STATUS: Confirmed

The power series sets tolerance = epsilon / value_of_rec(a) (L91-L92) and stops when a term is at or below it (L120).
The continued fraction sets tolerance = 3.0 * epsilon (L234) and tests C D against 1 and the derivatives relative to their size (L319-L325).

FACT: grad_reg_inc_beta keeps betaAB but no longer uses it

Quote:

grad_reg_inc_beta keeps its betaAB argument for interface compatibility but no longer uses it.

(PR description, What this PR does)

STATUS: Confirmed

The templated overload leaves the parameter unnamed (L426), and the doc comment says it is accepted for interface compatibility and not used (L417-L418).
The 9-argument overload passes it through (L498-L499).

FACT: The values of the distributions do not change

Quote:

The values of the distributions do not change.
Only the gradients change.

(PR description, Summary)

STATUS: Confirmed

The PR does not change any file under prim/prob/, and it does not change prim inc_beta.
Fwd inc_beta still takes its value from inc_beta(a_val, b_val, x_val) (L58).
The changed functions compute only derivatives.

FACT: The grad_reg_inc_beta test changes

Quote:

prim/fun/grad_reg_inc_beta_test.cpp: the two existing cases, with the z = 1 expectation corrected to 0 and the tolerance tightened to 1e-12.
New cases for each defect class: beta(a, b) underflow, prefactor underflow, integer b, and the tails.

(PR description, Testing)

STATUS: Confirmed

grad_reg_inc_beta.1 now expects 0 for both partials at z = 1 (L47-L63).
grad_reg_inc_beta.2 uses a relative tolerance of 1e-12 where develop used an absolute 1e-6 (L65-L67).
The new tests beta_ab_underflow, prefactor_underflow, integer_b, and tails are at L77-L109.

FACT: The inc_beta_dda and inc_beta_ddb tests are rewritten at the same twelve points

Quote:

prim/fun/inc_beta_dda_test.cpp, prim/fun/inc_beta_ddb_test.cpp: rewritten from the reference at the same twelve points, with the digamma_b argument corrected, plus regression points.

(PR description, Testing)

STATUS: Confirmed

Both files keep the shapes 1.5, 15000, 1.25, and 12500 and the z values of develop (0.6 for the middle z of inc_beta_dda and 0.5 for inc_beta_ddb), which give twelve points each.
The inc_beta_ddb helper passes digamma(b) (L15).
Each file adds a regressions test with nine points (inc_beta_dda_test.cpp L60-L75, inc_beta_ddb_test.cpp L61-L72).

FACT: No change to grad_2F1 or to how the neg_binomial cdfs call the root

Quote:

No change to grad_2F1, which other functions use.

No change to how the neg_binomial*_cdf functions call the root.
They need one derivative and already call the root once.

(PR description, Not included)

STATUS: Confirmed

The PR does not change prim/fun/grad_2F1.hpp, prim/prob/neg_binomial_cdf.hpp, or prim/prob/neg_binomial_2_cdf.hpp.
Each of the two cdfs calls inc_beta_dda once per observation (neg_binomial_cdf.hpp L98, neg_binomial_2_cdf.hpp L107).

FACT: The commit routes inc_beta_dda, inc_beta_ddb, and fwd inc_beta through grad_reg_inc_beta

Quote:

Port the Boik and Robinson-Cox continued fraction (as in TensorFlow Probability) to prim grad_reg_inc_beta, route inc_beta_dda, inc_beta_ddb and fwd inc_beta through it, and test against fixed mpmath references.

(commit 2764cf58)

STATUS: Confirmed

inc_beta_dda and inc_beta_ddb call grad_reg_inc_beta with G1 and G2 (inc_beta_dda.hpp L41, inc_beta_ddb.hpp L37).
Fwd inc_beta calls it with BOTH, G1, or G2, depending on which shapes are fvar (L40-L54).
The tests compare with fixed constants.

Partially Confirmed Facts

FACT: No signature changes

Quote:

No signature changes.

(PR description, Summary)

STATUS: Partially Confirmed

The 9-argument grad_reg_inc_beta(g1, g2, ...) overload keeps its parameters (L494-L497), and the five-argument calls of inc_beta_dda and inc_beta_ddb keep their order.
The templates do change.
inc_beta_dda and inc_beta_ddb move from one type T to five template types and return auto, and inc_beta_ddz changes from an explicit specialization to an overload.
Fwd inc_beta replaces its seven overloads with one constrained template.
The PR also adds the GRAD_RIB overload of grad_reg_inc_beta and container overloads of inc_beta_dda and inc_beta_ddz.

Request: Could you reword this to say that existing calls keep their form, and list the new overloads?

FACT: The code marks both changes to the TFP stop tests

Quote:

The port changes the two stop tests of the TFP source, one in the series and one in the continued fraction.
The code marks both changes.

(PR description, What this PR does)

STATUS: Partially Confirmed

Both stop tests also wait for the b derivative to converge.
The power series does this at L106-L121 and the continued fraction at L319-L325.
The doc comments explain why for integer b (L53-L59 and L195-L202).
The code does not say that these tests differ from the TFP source, and TFP is named only at L364.

Request: Could you add a note at the two stop tests saying that they differ from TFP, or reword the description?

FACT: prim/fun.hpp now includes the three derivative headers

Quote:

prim/fun.hpp now includes inc_beta_dda.hpp, inc_beta_ddb.hpp and inc_beta_ddz.hpp; before, it reached them only through prob/beta_cdf.hpp.

(PR description, What this PR does)

STATUS: Partially Confirmed

prim/fun.hpp now includes inc_beta_dda.hpp, inc_beta_ddb.hpp, and inc_beta_ddz.hpp (L137-L139).
On develop the headers were included from prob/beta_cdf.hpp and also from prob/neg_binomial_cdf.hpp and prob/neg_binomial_2_cdf.hpp, so beta_cdf.hpp was not the only route.

Request: Could you change "only through prob/beta_cdf.hpp" to name the three cdf headers?

FACT: About 180 lines of algorithm replace 123

Quote:

Net change: about 180 lines of algorithm replace 123.

(PR description, What this PR does)

STATUS: Partially Confirmed

The new code is longer than the code it replaces.
The diff of grad_reg_inc_beta.hpp is +472/-28, and inc_beta_dda.hpp and inc_beta_ddb.hpp shrink to wrappers.
Much of the new file is doc comments, and the description does not say how the lines of algorithm were counted.

Request: Could you say how the 180 and 123 lines were counted, or drop the sentence?

FACT: beta_lccdf computes 1 - Pn in linear space

Quote:

beta_lccdf computes 1 - Pn in linear space after inc_beta.
At (352, 590, 0.758) the double I is exactly 1, so the value is -inf and the gradients are infinite.

(PR description, Known limits)

STATUS: Partially Confirmed

beta_lccdf sets Pn = 1.0 - inc_beta(alpha_dbl, beta_dbl, y_dbl) and takes log(Pn) (L102-L106).
The PR does not include code that shows I rounding to 1 at (352, 590, 0.758).

Request: Could you point to a test of the value at (352, 590, 0.758), for example in #3433?

Not Verified Facts

FACT: Three new rev test files

Quote:

rev/fun/inc_beta_test.cpp (new): shape gradients of inc_beta with var arguments at the defect points.
[...]
rev/prob/beta_cdf_test.cpp (new) [...]
rev/prob/beta_cdf_log_test.cpp (new)

(PR description, Testing)

STATUS: Not Verified

test/unit/math/rev/fun/inc_beta_test.cpp, test/unit/math/rev/prob/beta_cdf_test.cpp, and test/unit/math/rev/prob/beta_cdf_log_test.cpp are not in this PR.
They belong to #3489 and #3490.

Request: Could you mark these three files as part of #3489 and #3490 in the Testing section?

Contradicted Facts

FACT: beta_cdf calls grad_reg_inc_beta once

Quote:

beta_cdf calls grad_reg_inc_beta once, not inc_beta_dda and inc_beta_ddb separately.

(PR description, What this PR does)

STATUS: Contradicted

At this commit beta_cdf still calls inc_beta_dda and inc_beta_ddb separately (L124-L131).
The opening note says the description covers all three parts, and this change is in #3490.

Request: Could you mark the statements that describe #3489 and #3490, such as this one?
A reader could then tell which code this PR holds.

FACT: The continued fraction switches to the symmetry relation at the mode

Quote:

It applies the symmetry relation I_z(a, b) = 1 - I_{1-z}(b, a) at the mode for the continued fraction

(PR description, What this PR does)

STATUS: Contradicted

The continued fraction switches when z >= (a + 1) / (a + b + 2) (L459).
The doc comment says that the mode (a - 1) / (a + b - 2) is not used as the threshold, because it exceeds 1 for b < 1 (L390-L396).

Request: Could you change "at the mode" to the threshold (a + 1) / (a + b + 2) that the code uses?

FACT: The file carries the TFP copyright notice and refers to the licence file

Quote:

The file carries the TFP copyright notice and refers to licenses/tensorflow-probability-license.txt.
It uses the same form as the Boost-derived code in prim/fun/log_modified_bessel_first_kind.hpp.

(PR description, Licence note)

STATUS: Contradicted

At this commit grad_reg_inc_beta.hpp has no TFP copyright notice and does not refer to a licence file (L1-L21).
TensorFlow Probability appears only in the doc comment at L364.
licenses/tensorflow-probability-license.txt does not exist in the tree.

Request: Could you add the TFP copyright notice to grad_reg_inc_beta.hpp and add licenses/tensorflow-probability-license.txt, as the licence note describes?

Not Checked

The checklist items in the description (header checks, dependency checks, doxygen, cpplint).
CI reports 34 successful checks, and I did not run them again.
Whether the code matches _betainc_partials in TensorFlow Probability line by line.
The PR does not include the TFP source, and this fact check reads only the code at the reviewed commit.
The stan-buildbot end-to-end benchmarks.
They do not call the changed functions, so they are not evidence for the speed claims.

…FP notice

avehtari · 2026-10-29T07:44:09Z

Reply to @stevieb-bot. I changed the PR text and the tests, and added the scripts behind the numbers. In the code, I changed only comments and the license notice of grad_reg_inc_beta.hpp, and restored the license file. The new head is 7308fa5. "Develop" in this reply is the merge base 579253050c that you used. The PR text now describes part 1 only. @SteveBronder created #3489 and #3490 for parts 1 and 2, but it is not yet clear who should update them.

Scripts

The scripts are in https://github.com/avehtari/stan-pr-checks/tree/8b80a636364b7751bbf3553e63cf18d010f394c2/pr3409, and its README has one section per claim with the commands and the expected output. The directory has:

refs/: the mpmath references (both routes) and the driver that computes the constants of the unit tests;
sweep/: the 18 000 points, their references, the probe, and the script that prints the accuracy table;
points/: the defect examples, the throws per function, and the comparison of the scalar types;
tfp/: TensorFlow Probability at the three integer-shape points;
ddz/: inc_beta_ddz against Boost's default ibeta_derivative;
speed/: the speed benchmark and its SLURM job;
swap/: the regression check;
results/: the outputs that the PR text quotes.

I ran all of them again on 579253050c and on the new code.

Numbers that changed in the rerun

claim	old text	rerun	cause
path A non-finite (throws)	100	76 of 17 472; throws: path A 98, path B 24	the old probe called the four roots of a point in one try block, so a throw in `grad_reg_inc_beta` also marked path A
develop `d/dα` of `beta_cdf(0.76 \| 446, 5)`	8.9e-16	-1.8e-15	the old value was from develop at `5252d51d47`
gradient of `log(beta_cdf(...))` there	2.5e+31	-4.4e+31	the same
true `d/dα` of `neg_binomial_lcdf(132 \| 313, 0.068)`	-2.3	-2.40	wrong in the old text
`log(1 - I)` at `beta_lccdf(0.758 \| 352, 590)`	-0.34	-316.2	wrong in the old text; computed now as `log I_{1-z}(b, a)` at 80 digits
worst error of the gradient of `log I`	1e-11	1e-12	rerun on the current code
worst relative error, derivative not negligible	3.7e-9 at `b ≫ a`, 6e-13 elsewhere	3.8e-12, at large shapes	the 3.7e-9 is from the code before `d2eeacd676`, which changed the prim continued fraction to carry logarithmic derivatives; the September scores already had 1.8e-12 to 3.8e-12 at large shapes, so the 6e-13 was a wrong summary
worst error in the first region	1e-13	5e-13	the old table gave `d/da` only; the new one covers both partials
speed, `grad_reg_inc_beta` in `double`	30x	54x	new run: this PR alone, Xeon Gold 6230, g++ 14.3, `-O3`
speed, `beta_lcdf` with shapes 5, 5	40x	56x	the same

At neg_binomial_lcdf(132 | 313, 0.068) the cause on develop also differs from the old text: I is 1e-262, and B(a, b) I underflows. The PR text now gives this example its own paragraph.

At the worst relative errors the two reference routes agree to 60 digits, so the errors are in the code. I fixed one error in the reference script. Where I > 1/2, the script differentiates 1 - I_{1-z}(b, a) and negates the result. It did the negation at 15 digits, which rounded the reference to a double. The tests and the sweep use doubles, so nothing else changed. With the fix, the two routes agree to at least 60 digits at all 21 test points.

Partially confirmed and not verified facts

Defects 1 to 5, and the regression check: as you found, the files with the defect tests did not compile on develop, because they also tested the new overloads. Those tests are now in grad_reg_inc_beta_one_sided_test.cpp, inc_beta_dda_vectorized_test.cpp and inc_beta_ddz_vectorized_test.cpp. The six other test files compile against develop, and 15 of their 19 tests fail there. Among the develop values in the failure messages are -0.2027 at (1.145, 6786, 1.8e-4), -1.8e-15 at (446, 5, 0.76) and 1.8226 at (1.5, 12500, 0.001). The Testing section gives the counts for this PR alone.
Defects 1 and 2, the distribution examples: rev/prob/beta_cdf_log_test.cpp (the file of Move the beta, neg_binomial and student_t cdfs to the new grad_reg_inc_beta (part 3 of #3409) #3490) tests the shape gradients of beta_lcdf at (0.5 | 600, 600) and (0.758 | 352, 590), and rev/prob/neg_binomial_cdf_log_test.cpp tests neg_binomial_lcdf(132 | 313, 0.068). Both fail on develop and pass with this PR alone. The y partial of beta_lcdf(0.5 | 600, 600) and the β partial of neg_binomial_lcdf(132 | 313, 0.068) stay NaN, also with Move the beta, neg_binomial and student_t cdfs to the new grad_reg_inc_beta (part 3 of #3409) #3490, because of the density terms of the lcdfs. The PR text lists them as a known limit, and I will fix them in a separate PR.
Defect 4, the iteration limit: inc_beta_dda_upper_tail_no_throw and inc_beta_ddb_upper_tail_no_throw use a point where develop throws "did not converge". The derivatives at all 98 such sweep points are below 1e-300 in magnitude; the test point has the largest, -3.9e-312 and 6.3e-312. These are subnormal, so the tests use a relative tolerance of 1e-9. grad_reg_inc_beta.upper_tail_no_throw uses a point where path B throws.
inc_beta_ddz returns infinity: inc_beta_ddz_infinite_density now compiles on develop and fails there with Boost's "Overflow Error".
The four scalar types: the claim was too broad for this PR. With both shapes as autodiff arguments, the gradients of inc_beta in the three autodiff types are bit-identical to the double root. With one fvar shape, fwd inc_beta calls the one-sided root, which differs by up to 2.5e-13 relative. The PR text now says so, and points/probe_types.cpp measures it.
The 185 generated programs: my comment of 2026-10-22 described the branch before the split, and its rev and mix suites belong to Reverse-mode grad_reg_inc_beta (part 2 of #3409) #3489. The generated tests of the six families pass with the new code (256 programs, Triton).
The references, the 18 000 points, the TFP errors and the ibeta_derivative grid: see refs/, sweep/, tfp/ and ddz/. I have not filed the TFP report yet.
Function and distribution speed: speed/probe_cmp_beta_speed.cpp with speed/run_cmp_beta_speed.sh, Triton job 21271376. Every row measures this PR alone, and the PR text notes that beta_cdf still calls inc_beta_dda and inc_beta_ddb separately here.
The 607 607 terms and the 600 cycles: I removed both numbers.
My speed comment of 2026-10-22 compares the original port with the reverse-mode code of Reverse-mode grad_reg_inc_beta (part 2 of #3409) #3489, so its table applies to Reverse-mode grad_reg_inc_beta (part 2 of #3409) #3489. The benchmark is speed/probe_speed3.cpp.

Description vs code

"No signature changes": the text now says that existing calls keep their form, and lists the new overloads and the changed templates.
The changes to the TFP stop tests: a comment at each stop test marks the change, and the doc comments of both kernels give the reason. The doc comment of the dispatch names the two other differences, the series limit z <= 0.5 and the switch point of the continued fraction.
prim/fun.hpp: the text now names the three cdf headers that included the derivative headers on develop.
"About 180 lines replace 123": removed.
beta_lccdf at (352, 590, 0.758): Evaluate the incomplete beta complement without cancellation in the lccdf and Student t cdf functions #3433 repairs this value.
The three rev test files: rev/fun/inc_beta_test.cpp and rev/prob/beta_cdf_test.cpp are removed from the list, since Reverse-mode grad_reg_inc_beta (part 2 of #3409) #3489 and Move the beta, neg_binomial and student_t cdfs to the new grad_reg_inc_beta (part 3 of #3409) #3490 hold them. The third, rev/prob/beta_cdf_log_test.cpp, is now also in this PR, unchanged from Move the beta, neg_binomial and student_t cdfs to the new grad_reg_inc_beta (part 3 of #3409) #3490.
beta_cdf calls the root once: that is part 3, and the opening note says so.
Symmetry "at the mode": the text now gives the switch point of the code, (a + 1) / (a + b + 2).
The TFP notice: the file had none. The September port had the notice and licenses/tensorflow-probability-license.txt, and @SteveBronder removed both in d2eeacd676. I restored them, so the license note matches the code again. The year is now 2020, as in the TFP source, and the license file holds the Apache-2.0 text without additions.

Tests

On Triton with g++ 14.3, these pass with the new code:

the 25 tests of the nine changed and new test files;
the unit suites of inc_beta and grad_2F1, and the mix tests of inc_beta and grad_reg_inc_beta;
the prob unit tests of the five families of the PR text;
the generated test/prob tests of the six families (256 programs);
cpplint, test-math-dependencies, test-headers and doxygen.

After that run I changed only comments, the year of the license notice and the name of one test. Then I ran the regression check on a laptop (g++ 11.4). All 25 tests pass with the new code, and 15 of 19 fail with the develop headers.

avehtari requested review from SteveBronder and andrjohns September 21, 2026 07:10

This was referenced Oct 6, 2026

lccdf of beta overflows at reasonably normal values #3377

Open

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

Draft

avehtari added the numerics Numerical issues label Oct 8, 2026

SteveBronder added this pull request to stack #3462 October 14, 2026 19:43

SteveBronder assigned avehtari Oct 14, 2026

avehtari mentioned this pull request Oct 21, 2026

Incomplete gamma gradient roots and log-space cdfs for six distributions #3408

Open

9 tasks done

Compute the incomplete beta shape gradients by continued fraction

2764cf5

Port the Boik and Robinson-Cox continued fraction (as in TensorFlow Probability) to prim grad_reg_inc_beta, route inc_beta_dda, inc_beta_ddb and fwd inc_beta through it, and test against fixed mpmath references.

SteveBronder force-pushed the stable-inc-beta branch from 6e1938c to 2764cf5 Compare October 26, 2026 19:08

This was referenced Oct 26, 2026

Reverse-mode grad_reg_inc_beta (part 2 of #3409) #3489

Draft

Move the beta, neg_binomial and student_t cdfs to the new grad_reg_inc_beta (part 3 of #3409) #3490

Draft

SteveBronder removed this pull request from stack #3462 October 26, 2026 19:14

SteveBronder added this pull request to stack #3491 October 26, 2026 19:14

SteveBronder added a commit that referenced this pull request Oct 26, 2026

Merge stable-inc-beta-cdfs into stable-beta-lccdf

709a1a4

Moves this branch onto the split of #3409 (#3409, #3489, #3490). The conflicts are resolved to the tree of these commits replayed onto stable-inc-beta-cdfs.

Split the overload tests, add throw and lcdf tests, and restore the T…

7308fa5

…FP notice

avehtari added a commit to avehtari/stan-pr-checks that referenced this pull request Oct 29, 2026

Add the scripts and results for stan-dev/math#3409

8b80a63

Conversation

avehtari commented Sep 21, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Summary

The defects in the current code on develop

1. grad_reg_inc_beta returns NaN for every parameter when beta(a, b) underflows

2. grad_reg_inc_beta drops the hypergeometric term when z^a (1-z)^b underflows

3. inc_beta_dda and inc_beta_ddb stop on an absolute threshold that the first term can meet

4. inc_beta_dda and inc_beta_ddb reflect on fixed z thresholds

5. The existing unit tests assert wrong values

What this PR does

Accuracy after the change

Speed

Testing

License note

Known limits, not addressed here

Not included

Release notes

Checklist

Uh oh!

avehtari commented Sep 21, 2026

Uh oh!

avehtari commented Oct 8, 2026

Uh oh!

avehtari commented Oct 22, 2026

Uh oh!

avehtari commented Oct 22, 2026

Uh oh!

andrjohns commented Oct 22, 2026

Uh oh!

avehtari commented Oct 22, 2026

Uh oh!

stan-buildbot commented Oct 22, 2026

Uh oh!

SteveBronder commented Oct 22, 2026

Uh oh!

stan-buildbot commented Oct 23, 2026

Uh oh!

andrjohns commented Oct 23, 2026

Uh oh!

stan-buildbot commented Oct 27, 2026

Uh oh!

stevieb-bot commented Oct 28, 2026

Summary

Claims and Evidence

FACT: The PR fixes #2417

STATUS: Confirmed

Evidence:

FACT: At z = 1 the shape gradients of inc_beta are 0 and the z gradient is +inf

STATUS: Confirmed

Evidence:

FACT: With a zero tangent on z, the fvar tangent of inc_beta at z = 1 is NaN

STATUS: Confirmed

Evidence:

FACT: Defect 1: grad_reg_inc_beta returns NaN when beta(a, b) underflows

STATUS: Partially Confirmed

Evidence:

FACT: Defect 2: grad_reg_inc_beta drops the hypergeometric term when the prefactor underflows

STATUS: Partially Confirmed

Evidence:

FACT: Defect 3: inc_beta_dda and inc_beta_ddb stop on an absolute threshold that can be met before the first term

STATUS: Partially Confirmed

Evidence:

FACT: Defect 4: inc_beta_dda and inc_beta_ddb reflect on fixed z thresholds

STATUS: Partially Confirmed

Evidence:

FACT: Defect 5: the existing unit tests assert wrong values

STATUS: Partially Confirmed

Evidence:

FACT: The four scalar types give bit-identical values on both trees

STATUS: Partially Confirmed

Evidence:

FACT: Regression check: the tests fail on develop and pass on the branch

STATUS: Partially Confirmed

Evidence:

FACT: inc_beta_ddz returns infinity instead of throwing where the density is infinite

STATUS: Partially Confirmed

Evidence:

FACT: The generated cdf, lcdf, and lccdf tests of six families pass

STATUS: Partially Confirmed

avehtari commented Sep 21, 2026 •

edited

Loading

The defects in the current code on `develop`

1. `grad_reg_inc_beta` returns NaN for every parameter when `beta(a, b)` underflows

2. `grad_reg_inc_beta` drops the hypergeometric term when `z^a (1-z)^b` underflows

3. `inc_beta_dda` and `inc_beta_ddb` stop on an absolute threshold that the first term can meet

4. `inc_beta_dda` and `inc_beta_ddb` reflect on fixed `z` thresholds