[avehtari]

EDIT: 2026-10-29 updated after revision

Fixes #3153. Fixes #1578. Addresses #2508 partially.

This PR improves beta_binomial_lpmf, beta_binomial_cdf, beta_binomial_lcdf, beta_binomial_lccdf, beta_neg_binomial_lpmf, beta_neg_binomial_cdf, beta_neg_binomial_lcdf, beta_neg_binomial_lccdf, neg_binomial_lpmf, neg_binomial_2_lpmf, neg_binomial_2_log_lpmf, neg_binomial_2_log_glm_lpmf, dirichlet_multinomial_lpmf, lkj_corr_lpdf, lkj_corr_cholesky_lpdf, lkj_cov_lpdf, student_t_lpdf, multi_student_t_lpdf, multi_student_t_cholesky_lpdf, yule_simon_lpmf, yule_simon_cdf, yule_simon_lcdf, yule_simon_lccdf, and the derivatives of lbeta.

Claude assisted with this PR.

The scripts behind the numbers in this text are in stan-pr-checks/pr3434. Its README has one section per claim, with the commands, the expected output and the commits. "Develop" in this text is the merge base 510d263078.

Summary

Several distributions form a difference f(x + k) - f(x). Here f is lgamma, lbeta or digamma, x is a shape, concentration, dispersion or degrees-of-freedom parameter, and k is a count or an offset. When x is large, the two terms are large and almost equal, and the difference keeps no correct digits. In the values the terms have the size of x, so the error is about eps * x. In the gradients the terms are about log(x), and a model with the parameter on the log scale multiplies the error by x.

In the model of #3153, beta_binomial_lpmf(57 | 117, s, s) with s = exp(lc) / 2, develop returns a log probability 81 too high from lc = 42. Its gradient in lc is 0 or noise from lc = 24, so warmup stops in a region where the log density is flat. With this PR, for lc from 0 to 44 in steps of 0.25, the value is within 4.3e-14 of the exact value. The gradient is within 3.7e-14.

lc	develop: diff to binomial limit	this PR	exact	develop d/dlc	this PR	exact
20	-4.4e-08	-1.1130e-07	-1.1130e-07	1.7e-06	1.1130e-07	1.1130e-07
28	-4.6e-05	-3.7309e-11	-3.7338e-11	0	3.7357e-11	3.7338e-11
32	+4.5e-03	-6.5937e-13	-6.8386e-13	-0.14	6.8570e-13	6.8386e-13
42	+81.098	+2.3e-14	-3.1e-17	0	0	3.1e-17

From lc = 34, the exact values are below the rounding error of any form (about 2e-14), because 57 log p + 60 log q is about -81 there.

#2508 reports the same defect: beta_binomial_lpmf(500 | 1000, s, s) for s = 0.5e10, 0.5e13, 0.5e19. Develop is wrong by 5.3e-8, 2.8e-4 and 693 there, and returns +689.47 at the largest s. This PR fixes all three points, and they are rows of prim/prob/beta_binomial_test.cpp. This PR does not close #2508. The comment there about beta_binomial_lcdf(3 | 10, 298.3, 146.75) (NaN, exact -3.93335) has a different cause. hypergeometric_3F2 does not stop at the zero term of the terminating series 3F2({1, 302.3, -6}; {5, -151.75}; 1) (exact 17.63), and a later 0 * inf gives NaN. #3437 repairs that and closes #2508; its test has this point.

#1578 reports wrong derivatives of neg_binomial_lpmf for large alpha. Its example (n = 0, alpha = 1e10 + 1, beta = 0.8) is already correct on develop, because develop no longer uses the Poisson approximation that the issue describes. The example is now a row of large_shapes_test.cpp, and it passes on develop too. This PR repairs the cancellation that remains in d/dalpha = digamma(alpha + n) - digamma(alpha) + log(beta / (1 + beta)). The cancellation matters when the mean alpha / beta stays fixed. There, the error of alpha * d/dalpha on develop reaches 3.0 at alpha = 1e15. With this PR it is at most 1.0e-14 over 42 points with the mean fixed at 3 or 50.

The same root cause is in more functions. The table gives the worst absolute error over a sweep of each function. Where the value or the gradient is large, it also gives the relative error at that point. For shape-type arguments, "gradient" is the log-shape gradient x * d/dx, the gradient that a sampler sees when the model parameter is log(x).

function	develop	this PR
beta_binomial_lpmf (shapes up to 1e18)	value 220, gradient 7.7e3	value 7.1e-9 (rel 2.8e-16), gradient 3.5e-10 (rel 2.3e-14)
beta_neg_binomial_lpmf	value 10, gradient 13	value 2.6e-2 (rel 4.9e-17), gradient 8.2e-2 (rel 2.8e-16)
neg_binomial_lpmf, neg_binomial_2_lpmf, neg_binomial_2_log_lpmf	gradient 3.0, 27, 26	1.0e-14, 1.7e-14, 2.0e-13
neg_binomial_2_log_glm_lpmf	value 23, gradient 7.1	4.5e-13, 2.8e-14
dirichlet_multinomial_lpmf	gradient 2.8	8.1e-15
lkj_corr_lpdf, lkj_corr_cholesky_lpdf	value 3.9, gradient 8.4 (rel 1.5e-14); at the identity matrix gradient 1.5 (rel 1.0)	value 0.15 (rel 2.7e-16), gradient 0.14 (rel 2.6e-16)
yule_simon_lccdf	value 7.6, gradient 2.1	1.7e-13, 8.2e-15
yule_simon_cdf, yule_simon_lcdf	gradient rel 1.0 (abs 5e-15)	gradient rel 1.3e-13
beta_binomial_{cdf,lcdf,lccdf}	gradient 19, 4.1e4, 22; lcdf NaN	value 1.8e-7, gradient 2.3e-6 (3F2 series)
beta_neg_binomial_{cdf,lccdf}	gradient 15, 13	value 2.0e-9, gradient 3.4e-6 (3F2 series)
beta_neg_binomial_lcdf	gradient 22 at shape 1e15, 43 at shape 100	2.0e-6 at shape 1e15; 69 at shape 100, where both are wrong because of the complement (see "Limits")
student_t_lpdf	value 0.16, nu gradient 0.14	1.3e-11 (rel 1.2e-16), 4.5e-12
multi_student_t_lpdf, multi_student_t_cholesky_lpdf	value 1.8e-3, nu gradient 8.9e-5	8.1e-15, 6.3e-15
yule_simon_lpmf	value 3.9e-3, gradient 3.6e-3	value 6.9e-10 (rel 2.1e-16), gradient 1.7e-10 (rel 1.9e-16)
lbeta (var, fvar), derivative	0.63	1.0e-14

The parameters go up to 1e15, with these exceptions. The beta_binomial_lpmf sweep goes up to 1e18, because the model of #3153 reaches 8.7e17 at lc = 42. Up to 1e15, the develop errors of this row are 0.19 and 7.6. The student_t_lpdf and lbeta rows go up to 1e14. The multi_student_t and yule_simon_lpmf rows go up to 1e12, with counts up to 1e6. For the lpmf and lpdf rows, where the PR column has an absolute error above 1e-12, the value or the gradient is of order 1e5 to 1e15. There the relative error is at the rounding level. yule_simon_cdf and _lcdf are 1 - c and log1p(-c) with a tiny c = n B(n, alpha + 1), so their errors show only in relative terms. That row comes from 24 points per function with alpha up to 1e15 and n up to 20 (review_rows/). The cdfs of beta_binomial and beta_neg_binomial keep the tolerance of their 3F2 series.

All references come from two independent mpmath routes that agree. The values use loggamma, the gradients closed-form digamma differences and mp.diff, and the cdfs sums of the pmf. The references are at 80 to 140 digits; those of student_t_lpdf, multi_student_t, yule_simon_lpmf and lbeta are at 50 digits, checked at 110. At the sweep points of the table, no PR error is above 4 times the develop error. The exceptions are at the rounding level of the point or within the series tolerance of a cdf. At these points the PR has no more NaN values or throws than develop (beta_binomial_lcdf: 16 NaN on develop, 0 with this PR). For the beta_neg_binomial_cdf grid, see "Limits that remain".

neg_binomial_lpmf, neg_binomial_2_lpmf and neg_binomial_2_log_lpmf also lost the digits of phi at small phi, and the GLM partials overflowed at large theta. This PR repairs both (items 5 and 6 of "Related defects fixed in this PR").

propto. In some of the new forms, a constant that propto drops cancels against terms that depend on a parameter. The code cannot skip such a constant. For example, lbeta(n, alpha + 1) contains lgamma(n). To keep the propto behavior, the code subtracts such a constant from the result when propto drops it. In yule_simon_lpmf this is lgamma(n). In neg_binomial_2_log_glm_lpmf it is lgamma(y + 1), phi log(phi) (data phi) and y theta (data x, alpha and beta), also in OpenCL. For the GLM with data phi, this adds a rounding error of about eps * phi log(phi) to the propto value. In a probe with five observations, the error at phi = 1e12 is 5.2e-3, against eps |propto value| = 3.1e-2; develop's error has the same size. A probe of all 14 changed densities uses every combination of data and parameter arguments. It matches the propto and full values of develop in all 456 values. A new test checks full - propto against the dropped terms.

Most of the added source lines are in these files:

file	added lines
prim/fun/log_rising_factorial_ratio.hpp (new helpers)	276
prim/fun/digamma_diff.hpp (new)	199
opencl/kernels/device_functions/log_beta_ratio.hpp (new)	168
prim/prob/neg_binomial_2_log_glm_lpmf.hpp	155
opencl/kernels/device_functions/digamma_diff.hpp (new)	115
opencl/kernel_generator/elt_function_cl.hpp (registrations)	104
student_t_lpdf (series), CPU and OpenCL	89 + 41

The other 25 distribution files (CPU, OpenCL and the GLM kernel) add 3 to 56 lines each. The rev and fwd lbeta files add 28 and 27.

What this PR adds

digamma_diff(x, d) = digamma(x + d) - digamma(x), a new function in prim/fun. It is public, but not exposed in the Stan language, because 19 distribution files and the derivatives of lbeta use it. log_rising_factorial can use it in a follow-up. It has three paths:

• a non-autodiff integer d from 0 to 8 (a count): sum_{j < d} 1 / (x + j). All terms are positive, and d = 1 gives exactly 1 / x;
• x < 10 and d >= 10: the plain difference, which loses only a few ulp there (the result is at least 0.72 and digamma(x) <= 2.25);
• otherwise: the recurrence shifts both arguments to at least 10, then log1p(d / y) plus the asymptotic series, with the differences y^-2i - (y + d)^-2i formed as (u - v) h_i so that nothing cancels.

Over 7 281 points, x from 1e-3 to 1e18 and d from 1e-6 to 1e7, the error is at most 4.2 ulp; the plain difference reaches 6.5e16 ulp. The function is vectorized and works with all autodiff types through generic arithmetic, with expect_ad tests. The PR also adds an OpenCL device function with the same paths and a kernel-generator registration.

internal::log_rising_factorial_ratio, internal::log_beta_ratio and their helpers in prim/fun/log_rising_factorial_ratio.hpp. They write each lgamma through lgamma_stirling_diff, so that the terms linear in the arguments cancel exactly. They form the rest without cancellation: they split off the x log1p(k / x) - k parts for k < x and add the integer parts exactly. log_beta_ratio(alpha, beta, n, m) = lbeta(alpha + n, beta + m) - lbeta(alpha, beta) selects the Stirling form or the plain lbeta form by the size of their terms. The Stirling form has rounding errors of the size of n + m, and the lbeta form of the size of the smaller shape. With the Stirling form at all points, the error rose to 6.1e-11 relative (2.8e-10 absolute) at 131 of the 13 177 sweep points where develop was correct. All 131 are at N = 1000 and N = 1e6, and the selection removes this loss.

The rev and fwd derivatives of lbeta. The partial in an argument of at least 10 is now -digamma_diff(a, b). Below 10 it stays digamma(a) - digamma(a + b), and the two partials share digamma(a + b) as before. lbeta itself is accurate at large arguments, but its derivative had the same cancellation. So every model that calls lbeta with a large parameter had inaccurate gradients, and so did multi_student_t_lpdf, which gets its nu gradient by autodiff through lbeta. Below 10 the plain difference has an absolute error of a few eps, which stays at the rounding level on the log scale too.

Changes per function

• beta_binomial_lpmf: value by log_beta_ratio, partials by digamma_diff.
• beta_neg_binomial_lpmf: value as -log(n) - lbeta(n, beta) + log((r)_n / (r + alpha + beta)_n) + log((alpha)_r / (alpha + beta)_r), partials by digamma_diff.
• neg_binomial_lpmf: alpha partial by digamma_diff. The beta partial term alpha / beta - alpha / (1 + beta) is now alpha / beta / (1 + beta). The value term lchoose(n + alpha - 1, n) is now -lbeta(alpha, n + 1) - log(n + alpha).
• neg_binomial_2_lpmf, neg_binomial_2_log_lpmf: phi partial by digamma_diff. The value term lchoose(n + phi - 1, n) is now -lbeta(phi, n + 1) - log(n + phi). When phi is autodiff, neg_binomial_2_lpmf now evaluates log(mu + phi) once, and the value and the phi partial both read it; the lazy expression computed the log twice per element. The value and every adjoint are bitwise identical before and after this change (2 000 instances, five call types).
• neg_binomial_2_log_glm_lpmf: value per instance as -lbeta(phi, y + 1) - log(y + phi) - y log1p_exp(log(phi) - theta) - phi log1p_exp(theta - log(phi)), the form of neg_binomial_2_log_lpmf, instead of phi log(phi) - lgamma(phi) + lgamma(y + phi) - (y + phi) log(mu + phi) + y theta. The partials do not form exp(theta). With u = theta - log(phi) and e = exp(-|u|), mu / (mu + phi) is 1 / (1 + e) for u > 0 and e / (1 + e) otherwise, and phi / (mu + phi) is the other of the two. The theta partial is y - (y + phi) mu / (mu + phi), and the phi partial is mu / (mu + phi) - (y / phi) phi / (mu + phi) - log1p_exp(u) + digamma_diff(phi, y) (formulas from @jachymb's review). The value and the partials share e and log1p(e): log1p_exp(±u) = max(±u, 0) + log1p(e). With a scalar phi and a vector y, the terms that depend on y only (lchoose, lgamma(y + 1) under propto, digamma_diff(phi, y)) are evaluated once per distinct count, through internal::per_count_values, if the largest count is at most 4 N + 64. This is in the CPU version only.
• dirichlet_multinomial_lpmf: partials by digamma_diff.
• lkj_corr_lpdf (and lkj_corr_cholesky_lpdf, lkj_cov_lpdf through do_lkj_constant): the constant is a sum of lgamma(x + k/2) - lgamma(x). Each term is now lgamma(k/2) - lbeta(k/2, x), and the code sets the eta partial explicitly with digamma_diff.
• yule_simon_{cdf,lcdf,lccdf}: lgamma(alpha + 1) + lgamma(n + 1) - lgamma(n + alpha + 1) is now log(n) + lbeta(n, alpha + 1); partial by digamma_diff.
• yule_simon_lpmf: value log(alpha) + lbeta(n, alpha + 1); partial 1 / alpha - digamma_diff(alpha + 1, n).
• beta_binomial_{cdf,lcdf,lccdf}: the leading lbeta(nu, mu) - lbeta(alpha, beta) by log_beta_ratio, the digamma differences by digamma_diff.
• beta_neg_binomial_{cdf,lcdf,lccdf}: the leading term is the log pmf at n + 1, now computed as in beta_neg_binomial_lpmf; partials by digamma_diff.
• student_t_lpdf: unchanged for nu < 40. For nu >= 40, the terms lgamma((nu + 1) / 2) - lgamma(nu / 2) - log(nu) / 2 and digamma((nu + 1) / 2) - digamma(nu / 2) come from the asymptotic series in h = nu / 2, -log(2) / 2 - 1/(8h) + 1/(192h^3) - 1/(640h^5) + 17/(14336h^7) - 31/(18432h^9), and its derivative. Against mpmath at 2 012 points from h = 20 to 1e15, the errors of the two terms are at most 0.51 and 1.76 eps. The series is also cheaper than two lgamma or two digamma calls (see "Speed").
• multi_student_t_lpdf, multi_student_t_cholesky_lpdf: the constant lgamma((nu + p) / 2) - lgamma(nu / 2) is now lgamma(p / 2) - lbeta(p / 2, nu / 2). The cholesky version sets the nu partial with digamma_diff(nu / 2, p / 2); the other version gets it through lbeta.
• OpenCL: beta_binomial_lpmf, neg_binomial_lpmf, neg_binomial_2_lpmf, neg_binomial_2_log_lpmf, neg_binomial_2_log_glm_lpmf (prim file and kernel) and student_t_lpdf compute the same formulas.

Related defects fixed in this PR

1. lkj_corr_lpdf dropped d constant / d eta at eta == 1.0. The special case for eta == 1 returned the constant as a number without a derivative. At eta = 1 with K = 3 and the matrix in the new test, develop returns d/deta = -0.165; the exact value is 1.2214. eta = exp(0) = 1 is the initial value of a positive parameter with init = 0.
2. beta_binomial_lcdf computed log(1 - C), with C the complementary series. Where the cdf is small, this keeps no correct digits, and with the accurate C it can become log of a value at or below 0. Where C > 1/2, the function now uses the mirror relation X <= n if and only if N - X > N - n - 1. N - X has the distribution beta_binomial(N, beta, alpha). That is the same series with the shapes swapped, without the complement. On develop this function returned NaN at 16 of 81 sweep points and was wrong by up to 74 at the others. With this PR it returns no NaN there, and its errors are those of the series (value 7.7e-8, gradient 2.3e-6). For alpha = 1 exactly, the mirrored series is 0 / 0 at one term, so that case keeps log(1 - C).
3. OpenCL beta_binomial_lpmf did not check 0 <= n <= N. The code built the check but never computed its result, so the function never returned LOG_ZERO. Without propto the value was -inf by accident (through binomial_coefficient_log), but with a nonzero gradient. Under propto it was finite: -8.216 for n = 6, N = 5 in the new test. The CPU version returns LOG_ZERO without a gradient, and now the OpenCL version does too.
4. OpenCL beta_binomial_lpmf with scalar autodiff shapes did not compile. On develop, the compiler stops at line 81 of opencl/prim/beta_binomial_lpmf.hpp with "no matching function for call to lbeta(const double&, const double&)". The new test file beta_binomial_lpmf_scalar_shapes_test.cpp calls the function with scalar var shapes.
5. neg_binomial_lpmf, neg_binomial_2_lpmf and neg_binomial_2_log_lpmf lost the digits of phi at small phi. They computed lchoose(n + phi - 1, n) as binomial_coefficient_log(n + phi - 1, n), which rounds phi into n + phi - 1. The value then has an error of about eps * n / phi. At n = 1000, phi = 1e-8, the error is 7.9e-7 (5.0e-9 in neg_binomial_lpmf). At n = 1, phi = 1e-17, the value is -inf; the exact value is -39.1439. These three functions and the GLM, CPU and OpenCL, now use -lbeta(phi, n + 1) - log(n + phi), and the errors at these points are at most 5.3e-13. The partials were already correct, because they use the exact phi.
6. The neg_binomial_2_log_glm_lpmf partials overflowed at large theta. On develop, the theta partial y - exp(theta) (y + phi) / (exp(theta) + phi) is NaN at theta = 720 (inf / inf), and -inf at theta = 709, where exp(theta) (y + phi) overflows. The exact value is -1.5 at both points (y = 5 and 3, phi = 1.5). The e form in "Changes per function" has no overflow, and its second derivative does not cancel. At y = 1, phi = 0.07, the fvar<var> Hessian entry in theta is now -4.7225e-17 at theta = 35 (develop: 0) and -1.40777e-13 at theta = 27 (develop: -1.40970e-13). A new test checks both.

Speed

For nu < 40, student_t_lpdf computes the same formulas as develop, plus the test for nu >= 40. This adds 2.3 % instructions for the value and 3.7 % for value and gradient, and the time is 1.01 and 1.06 times that of develop. From nu = 40 it executes fewer instructions. At nu >= 100 it executes 72 % fewer for the value and 54 % fewer for value and gradient. The time is 0.33 and 0.54 times that of develop. The instruction counts are from perf stat and valgrind cachegrind, which agree. The times are from the machine below.

The table gives the time of this PR divided by the time of develop. Shapes are below 10 unless the row says otherwise. The times are medians of 5 interleaved rounds on one whole node (AMD EPYC 9965, g++ 11.5, -O2). The spread between rounds is at most 2.8 % for develop and 6.9 % for this PR.

benchmark	ratio
neg_binomial, neg_binomial_2, neg_binomial_2_log gradient, counts 0 to 8	0.88 to 0.90
neg_binomial_2 gradient, phi from 10 to 1e4, counts 0 to 8	0.83
neg_binomial_2 gradient, counts 9 to 100	0.92
neg_binomial_2 gradient, 1000 counts from 0 to 20, one phi	1.07
beta_binomial gradient, N = 20	1.09
lbeta gradient, both arguments below 10	1.01
lbeta gradient, one argument from 10 to 1000	1.21
neg_binomial_2_log_glm_lpmf value, N = 1000, K = 5, counts with mean 3	0.28
same, gradient in alpha, beta	0.30
same, gradient in alpha, beta, scalar phi	0.30
same, with propto	0.39
same, with propto, all 1000 counts distinct	1.45
same, with propto, counts from 5000 to 5999	1.33
same, gradient in alpha, beta, vector phi	0.69
same, gradient in x	0.42

With a scalar phi, the GLM evaluates lchoose, lgamma(y + 1) and digamma_diff(phi, y) once per distinct count. Where the counts rarely repeat or exceed 4 N + 64, the forms that repair the GLM (items 5 and 6) make it slower than develop.

Limits that remain, by a different cause

• beta_binomial_{cdf,lcdf,lccdf} stop their 3F2 series at a term below 1e-6. This leaves errors up to 2.3e-6 in the gradients. The functions also return NaN or throw at some arguments, for example at beta = 1, at alpha + beta = 1 and at the point in beta_binomial_lpmf is unstable for large shape arguments #2508. Fix and speed up the terminating 3F2 series of the beta_binomial cdfs #3437 repairs these.
• beta_neg_binomial_{cdf,lcdf,lccdf} stop their series at a term below 1e-8. They throw at some arguments from shape 1e4 and can take seconds per call. beta_neg_binomial_cdf and beta_neg_binomial_lcdf still form the complement 1 - ccdf, so they are wrong where the cdf is small. With the more accurate log pmf of this PR, 1 - ccdf rounds to exactly 0 where the cdf is below about 1e-15. On a wide grid of 2 352 points, 1 546 points finished within 2 s on both trees. At 15 of them, beta_neg_binomial_cdf then returns 0 with NaN gradients. There, develop returns wrong finite values (from -7.2e-11 to 5.4e-11) with finite gradients. Compute the beta_neg_binomial cdfs from sums of the pmf #3441 replaces the series and the complement and repairs all of these.
• The log-shape gradients have a floor that no implementation of these functions can remove. The function returns d/dalpha and d/dbeta, and autodiff forms alpha d/dalpha + beta d/dbeta outside it. For tied shapes in numerical instability in beta_binomial_lpmf at large shape parameters #3153, the two products are -1.5 and +1.5, and their sum is 27 / s = 54 / exp(lc). Each partial is accurate to a few eps of its own terms, so the sum is accurate to a few eps of 1.5, not of its own size.
• The neg_binomial_2_log_lpmf Hessian has the quotient-rule cancellation of item 6 at large |eta| (NaN entries at eta = 720 and -400, n = 3, phi = 0.07). Develop has the same defect. A follow-up PR can use the e form of the GLM there too.
• log_rising_factorial is still lgamma(x + n) - lgamma(x). No Stan Math function calls it. The rev and fwd derivatives of log_rising_factorial, rising_factorial, log_falling_factorial, falling_factorial and beta still form a plain digamma difference. A follow-up can use the new helpers there.

Tests

• test/unit/math/prim/fun/digamma_diff_test.cpp, test/unit/math/mix/fun/digamma_diff_test.cpp: special cases, exact finite sums, 17 mpmath values, the count path and the plain path (digamma_diff(x, 1) == 1 / x exactly from x = 1e-300 to 1e300); expect_ad at 8 points.
• test/unit/math/prim/prob/beta_binomial_test.cpp: values at the points of numerical instability in beta_binomial_lpmf at large shape parameters #3153 and beta_binomial_lpmf is unstable for large shape arguments #2508 and at large unequal shapes; the difference to the binomial limit at lc = 42.
• test/unit/math/rev/prob/beta_binomial_lpmf_test.cpp (new): log-shape gradients at 8 points; the lc gradient of the model in numerical instability in beta_binomial_lpmf at large shape parameters #3153 for lc = 16 ... 32. From lc = 34, develop returns exactly 0, which is within 1e-13 of the exact gradient, so a test there cannot see the defect.
• test/unit/math/rev/prob/large_shapes_test.cpp (new): 97 rows over the other functions, including student_t_lpdf, both multi_student_t functions, yule_simon_lpmf and lbeta. For most functions, the rows are points where develop was wrong by more than 100 times the tolerance, plus one point at a small shape where develop was correct. Three kinds of rows differ. The neg_binomial_lpmf row at the point of The derivatives of neg_binomial_lpmf are incorrect for large alpha #1578 passes on develop. The GLM rows at phi = 1e-17 and 1e-8 pass on develop too: they check a regression of an earlier commit of this PR. The yule_simon_cdf and _lcdf rows use a tolerance relative to the reference, because the cdf is 1 - c with a tiny c, and an absolute tolerance cannot see an error in c. The rows include the four negative binomial functions at phi = 1e-8 and 1e-17 and the GLM at theta = 720 and 709. The file also tests the lkj_corr gradient at eta = 1, and full - propto against the dropped terms for student_t_lpdf, yule_simon_lpmf, beta_neg_binomial_lpmf and three cases of the GLM.
• test/unit/math/rev/prob/neg_binomial_2_log_glm_lpmf_test.cpp: a vector call against the sum of single-instance calls, value and partials, with and without propto, once with repeated counts (per-count path) and once with a count above 4 N + 64.
• test/unit/math/mix/prob/neg_binomial_2_log_glm_lpmf_test.cpp: the fvar<var> second derivative in theta at theta = 35 and 27.
• test/unit/math/mix/fun/lbeta_test.cpp: the var and fvar<double> derivatives of lbeta at arguments from 1e8 to 1e14, in both argument orders, compared as x * d/dx with mpmath references.
• OpenCL: an absolute-reference test per changed distribution in test/unit/math/opencl/rev/, because the existing differential tests compare OpenCL with the CPU and cannot see an error that both share. The four negative binomial reference tests include the small-phi and large-theta rows. The beta_binomial_lpmf rows are in beta_binomial_lpmf_reference.hpp: beta_binomial_lpmf_test.cpp uses them with vector shapes, and beta_binomial_lpmf_scalar_shapes_test.cpp with scalar shapes (item 4). Also: a device-function test for digamma_diff; a test that beta_binomial_lpmf gives LOG_ZERO without a gradient for n outside 0...N, also under propto; and a comparison of the GLM with the CPU for a scalar alpha with a vector phi and the reverse, at N = 153 (more than one work group).

Swap tests. CPU: I replaced the 24 changed CPU files that exist on develop with their develop versions, all at once; the new digamma_diff.hpp and log_rising_factorial_ratio.hpp stay. With these files, 7 of the 9 new test cases fail: the 4 new beta_binomial cases, large_shapes_value_and_log_scale_gradients, lkj_corr_eta_gradient_at_one and the Hessian test. The other 2 pass, as intended. large_shapes_propto_dropped_terms checks full - propto against the dropped terms, which develop also satisfies. The GLM vector test checks that a vector call equals the single calls. In large_shapes_test.cpp, each of the 22 function tags has failing assertions. With the branch files restored, all tests pass again. A separate run with the develop rev and fwd lbeta files shows that the new lbeta test fails there too, and the two old tests of that file pass. OpenCL: I replaced the 32 changed files that exist on develop (24 CPU and 8 OpenCL files) with their develop versions. Then every new reference test fails (beta_binomial, the four negative binomial functions, student_t), and so does opencl_n_outside_support. beta_binomial_lpmf_scalar_shapes_test.cpp does not compile on develop (item 4). The differential tests pass, but with the CPU files swapped too they compare develop with develop.

Full suites. CPU builds on Triton batch nodes, g++ 14.3.0. At the PR head, 56 unit test files of the touched functions pass (prim, rev, fwd, mix, including lbeta, digamma_diff, student_t, multi_student_t, lkj_cov, glm_forward_mode_tape_test and the Laplace tests that use neg_binomial_2_log). The generated test/prob tests of beta_binomial, beta_neg_binomial, neg_binomial, neg_binomial_2, yule_simon and student_t pass (271 306 cases at 64407093ba; neg_binomial and neg_binomial_2 again at the head, 23 490 cases). The 9 OpenCL suites of the six distributions, digamma_diff and the kernel generator pass on a V100 (115 cases). make cpplint, make test-math-dependencies, make test-headers, the header checks of the changed files and Doxygen 1.9.7 (no warnings, at 64407093ba) pass.

Side effects

None

Release notes

Improved stability of beta_binomial_lpmf, beta_binomial_cdf, beta_binomial_lcdf, beta_binomial_lccdf, beta_neg_binomial_lpmf, beta_neg_binomial_cdf, beta_neg_binomial_lcdf, beta_neg_binomial_lccdf, neg_binomial_lpmf, neg_binomial_2_lpmf, neg_binomial_2_log_lpmf, neg_binomial_2_log_glm_lpmf, dirichlet_multinomial_lpmf, lkj_corr_lpdf, lkj_corr_cholesky_lpdf, lkj_cov_lpdf, student_t_lpdf, multi_student_t_lpdf, multi_student_t_cholesky_lpdf, yule_simon_lpmf, yule_simon_cdf, yule_simon_lcdf, yule_simon_lccdf, and the derivatives of lbeta.

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

[stan-buildbot]

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan	0.23	0.22	1.01	1.35% faster
stat_comp_benchmarks/benchmarks/gp_regr/gen_gp_data.stan	0.06	0.06	0.99	-0.85% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	6.31	6.28	1.0	0.35% faster
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	1.01	0.95% faster
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	8.41	8.41	1.0	0.04% faster
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	4.46	4.45	1.0	0.32% faster
stat_comp_benchmarks/benchmarks/sir/sir.stan	169.67	168.85	1.0	0.49% faster
stat_comp_benchmarks/benchmarks/garch/garch.stan	0.9	0.89	1.01	0.94% faster
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.63	0.62	1.02	1.53% faster
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	43.03	42.74	1.01	0.68% faster
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.59	0.59	1.01	0.63% faster
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	20.88	20.82	1.0	0.29% faster
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.11	0.11	1.0	0.18% faster
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.19	3.18	1.0	0.21% faster
performance.compilation	389.62	386.82	1.01	0.72% faster
Mean result: 1.0052754019026933

Jenkins Console Log
Jenkins Build Stages
Commit hash: 44a50912ff3f2e71023ea29a6f91adf7c6b1cf87

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.237
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

[jachymb]

@avehtari Hi! I was working on something related and took a similar, but a bit different approach. Here are some notes on how this can be further improved.

I have this on a branch with Hessian tests (η from −800 to 720) and the Laplace test. I can fold it into this PR or rebase it as a follow-up once this is merged, or you can paste these to your agent. lmk what is a preferred option. :)

Comments

Checked at 468ff26 with a probe (var gradients, stan::math::hessian for fvar<var> and fvar<fvar<double>>), the same probe on develop 6110e53, and mpmath as reference. The large-φ rows reproduce. These fail:

quantity	inputs	this branch	develop	exact
GLM ∂/∂φ, (theta_exp - y) / (theta_exp + phi), CPU and kernel	y=3, x=1, α=0, β=720 (data), φ=1.5	NaN	−717.2421539395109	−717.24215393951088
GLM ∂/∂β, y - theta_exp (y+φ) / (theta_exp+φ), CPU and kernel	y=5, x=1, α=0, β=720, φ=1.5	NaN	NaN	−1.5
same	β=709	−inf	−inf	−1.5
GLM ∂²/∂θ² (fvar<var> and fvar<fvar<double>>)	y=1, θ=35, φ=0.07	0	0	−4.7225e−17
same	y=1, θ=27, φ=0.07	−1.4097e−13	same	−1.40777e−13
GLM value, binomial_coefficient_log(y + phi - 1, y)	y=1, θ=0, φ=1e−17	−inf	−39.143946580898778	−39.143946580898777
neg_binomial_2_log_lpmf Hessian	n=3, φ=0.07, η=720 / η=−400	NaN entries	same	−4.4e−314 / −8.4e−173

Rows 1 and 6 are regressions of this branch (inf/inf; y + phi is rounded to ulp(y) before the call). The rest is exp(theta) on the tape, as on develop. The second-derivative error (quotient rule, relative error ≈ eps·e^θ/φ) makes block_matrix_sqrt reject a rank-deficient W in the embedded Laplace approximation: a Laplace test of mine (user likelihood calling the GLM, e^θ/φ ≈ 6e9) throws on this branch and on develop.

All of it goes away by never forming exp(theta): with $u = \theta - \log\phi$ and $e = e^{-|u|}$,

$$s=\frac{e^{\theta}}{e^{\theta}+\phi}=\begin{cases}1/(1+e) & u>0\ e/(1+e) & u\le 0\end{cases},\qquad 1-s=\begin{cases}e/(1+e) & u>0\ 1/(1+e) & u\le 0\end{cases}$$ $$\frac{\partial\log p}{\partial\theta}=y-(y+\phi)\,s,\qquad \frac{\partial\log p}{\partial\phi}=s-\frac{y}{\phi}\,(1-s)-\log\left(1+e^{u}\right)+\psi(y+\phi)-\psi(\phi)$$

where log(1+e^u) is your log1p_exp_theta_m_log_phi and ψ(y+φ)−ψ(φ) your digamma_diff(phi, y), both from the same e (log1p_exp(±u) = max(±u, 0) + log1p(e)); same in the kernel and in neg_binomial_2_log_lpmf. For the value, -lbeta(phi, y + 1) - log(y + phi) is what binomial_coefficient_log evaluates for y + φ ≥ 10, with φ unrounded.

[avehtari]

@jachymb thanks for the review and suggestions. I'll check these

[avehtari]

@jachymb Rows 1 and 6 are regressions of this PR, and I could reproduce them. Your formulas are correct, and I used them. In the other PR checking your comments lead to additional findings, so I did ask Claude to check your comments, too, but I guess here the fix is the same as you had. Below there is one request for you

• neg_binomial_2_log_glm_lpmf (CPU and OpenCL kernel): the partials use 1 / (1 + e) and e / (1 + e) with e = exp(-|theta - log(phi)|) and do not form exp(theta). This fixes rows 1 to 3. It also fixes the GLM second derivatives of rows 4 and 5, because the derivative of e / (1 + e) does not cancel: the fvar<var> Hessian entry is now -4.7225e-17 at theta = 35 (before: 0) and -1.40777e-13 at theta = 27 (before: -1.40970e-13).
• The value term lchoose(y + phi - 1, y) is now -lbeta(phi, y + 1) - log(y + phi). The defect is wider than row 6: binomial_coefficient_log(y + phi - 1, y) rounds phi into y + phi - 1, so the value has an error of about eps * y / phi. Develop has this form in neg_binomial_2_lpmf, neg_binomial_2_log_lpmf and neg_binomial_lpmf (with alpha), and this PR had brought it into the GLM. At y = 1000, phi = 1e-8 the error was 7.9e-7 (5.0e-9 in neg_binomial_lpmf), and at phi = 1e-17 the value was -inf. The four functions, CPU and OpenCL, now use the lbeta form; at these points the errors are at most 5.3e-13.
• New test rows with mpmath references in large_shapes_test.cpp and in the four OpenCL reference tests: the GLM at theta = 720 and 709, and the four functions at phi = 1e-8 and 1e-17. With the previous head of this PR exactly these rows fail. With the new commit all unit tests of the negative binomial family, the Laplace tests that use neg_binomial_2_log, the generated test/prob tests (23 490 cases) and the OpenCL suites (V100) pass.

Row 7, the neg_binomial_2_log_lpmf Hessian, has the same quotient-rule cancellation, and develop has it too. To keep this PR reviewable, could you (@jachymb) make it a follow-up PR after this one is merged, together with your Hessian tests and the Laplace test? Your e form applies there in the same way.

I also saw #3446 for the inv_logit derivatives. The GLM uses the explicit e / (1 + e) form, so it does not depend on #3446.

[WardBrian]

@avehtari could you resolve the conflicts in this PR to let jenkins run?

[avehtari]

conflict resolved (it was only in the include block)

[stan-buildbot]

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.99	-1.21% slower
stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan	0.16	0.17	0.96	-4.08% slower
stat_comp_benchmarks/benchmarks/garch/garch.stan	0.87	0.86	1.01	1.24% faster
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.01	0.01	0.94	-6.19% slower
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.11	0.11	1.02	2.42% faster
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	7.97	7.98	1.0	-0.11% slower
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	3.69	3.75	0.98	-1.71% slower
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	32.92	32.96	1.0	-0.1% slower
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.38	0.38	1.01	0.59% faster
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	6.79	6.83	0.99	-0.61% slower
stat_comp_benchmarks/benchmarks/sir/sir.stan	136.74	136.78	1.0	-0.03% slower
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.21	3.23	0.99	-0.77% slower
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.6	0.61	0.99	-0.78% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	25.57	25.26	1.01	1.23% faster
performance.compilation	387.23	379.43	1.02	2.01% faster
Mean result: 0.9950799155442522

Jenkins Console Log
Jenkins Build Stages
Commit hash: 29e32c4215464fd5550540ce9a448016d22be56f

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.526
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

[stevieb-bot]

Action: Not Ready for Human Review

Commit checked: c962588c6d0e

Summary

Thank you for this PR and for the mpmath reference values in the tests.

The PR says that differences of lgamma, lbeta and digamma at large shape parameters lose all digits on develop in 24 functions, and that two new helpers, digamma_diff and log_beta_ratio, repair the values and gradients.
It also lists six related defects that it fixes, a speed table, and runs of the new tests against develop.

I read the source and test changes at c962588.
In a container without network access, I ran the new and changed CPU tests at the head, and at the merge base d12f488 on develop with the PR's test files copied in.
The base run leaves out the two digamma_diff tests, which call a function develop does not have.
I ran the OpenCL tests the same way on an NVIDIA GPU.
At the head I also ran the 25 unit test files of the negative binomial family.
I also ran the generated test/prob tests of neg_binomial and neg_binomial_2 there.
The script runTests.py with --run-all exits with status 0 even when a test fails, so the results below come from the gtest lines in each output.

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

• Claims and Evidence: 9 Confirmed, 9 Partially Confirmed, 10 Not Verified, 1 Contradicted.
• Description vs Code: 12 Confirmed, 1 Partially Confirmed, 1 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 three points of #2508 are wrong on develop by 5.3e-8, 2.8e-4 and 693, and the PR fixes them

Quote:

#2508 reports the same defect: beta_binomial_lpmf(500 | 1000, s, s) for s = 0.5e10, 0.5e13, 0.5e19.
Develop is wrong by 5.3e-8, 2.8e-4 and 693 there, and returns +689.47 at the largest s.
This PR fixes all three points, and they are rows of the beta_binomial_lpmf test.

(PR description, Summary)

STATUS: Confirmed

The three points are rows 4 to 6 of large_shapes in the prim beta_binomial test, which checks them to 1e-13 relative.
They pass at the head.
On develop with the test copied in they fail by 5.27e-8, 2.80e-4 and 693.147, and develop returns 689.467 at s = 0.5e19, as the description says.

Evidence:

The test rows are in the PR, with reference values that agree with the binomial limit -3.67992 given in #2508.
Code: beta_binomial_test.cpp, large_shapes.

What I ran: python3 runTests.py -j8 --run-all test/unit/math/prim/prob/beta_binomial_test.cpp, with the other tests, at the head and on develop with the PR's tests copied in.
Output in the gist: cpu_tests-head.out, cpu_tests-base.out.

---

FACT: digamma_diff(x, 1) equals 1 / x exactly from x = 1e-300 to 1e300

Quote:

the count path and the plain path (digamma_diff(x, 1) == 1 / x exactly from x = 1e-300 to 1e300)

(PR description, Tests)

STATUS: Confirmed

digamma_diff_count_and_plain_paths checks with EXPECT_EQ that digamma_diff(x, 1) and digamma_diff(x, 1.0) equal 1.0 / x at x = 1e-300, 0.1, 3.7, 9.99, 1e8 and 1e300.
It passes at the head.

Evidence:

The test is in the PR.
Code: digamma_diff_test.cpp, count and plain paths.

What I read: digamma_diff.hpp L118-L125 sends an arithmetic d = 1 through the count path, which returns 0 + inv(x) for any positive x.

What I ran: python3 runTests.py -j8 --run-all test/unit/math/prim/fun/digamma_diff_test.cpp, with the other tests, at the head.
Output in the gist: cpu_tests-head.out.

---

FACT: lkj_corr_lpdf dropped d/deta at eta == 1, where develop returns -0.165 and the exact value is 1.2214

Quote:

The special case for eta == 1 returned the constant as a number without a derivative.
At eta = 1 with K = 3 and the matrix in the new test, develop returns d/deta = -0.165; the exact value is 1.2214.

(PR description, Related defects fixed on the way, item 1)

STATUS: Confirmed

ProbDistributions_lkj_corr_eta_gradient_at_one checks d/deta at eta = 1 against 1.2214197179296566 to 1e-14.
It passes at the head.
On develop with the test copied in it fails, with eta.adj() = -0.16487.
At the head, do_lkj_constant sets the eta partial in both branches, including eta == 1.

Evidence:

The test is in the PR.
Code: large_shapes_test.cpp, lkj_corr_eta_gradient_at_one.

What I read: lkj_corr_lpdf.hpp L53-L61 computes the eta partial with digamma_diff after the eta == 1 branch, so both branches get it.

What I ran: python3 runTests.py -j8 --run-all test/unit/math/rev/prob/large_shapes_test.cpp at the head and on develop with the test copied in.
Output in the gist: cpu_tests-head.out, cpu_tests-base.out.

---

FACT: OpenCL neg_binomial_2_log_glm_lpmf with a scalar alpha and a vector phi failed on develop

Quote:

The kernel decided from alpha instead of phi whether to sum the phi derivative.
It returned work-group sums where per-instance derivatives were necessary, and the gradient assignment threw a size mismatch.

(PR description, Related defects fixed on the way, item 4)

STATUS: Confirmed

opencl_matches_cpu_mixed_alpha_phi compares OpenCL with the CPU for a scalar alpha with a vector phi and for a vector alpha with a scalar phi, at N = 3 and N = 153.
It passes at the head on the GPU.
On develop with the test copied in, it fails with the exception "binary_operator+: Rows of a (3) and rows of b (1) must match in size", as the description says.

Evidence:

The test opencl_matches_cpu_mixed_alpha_phi is in the PR.
Code: opencl/rev/neg_binomial_2_log_glm_lpmf_test.cpp, opencl_matches_cpu_mixed_alpha_phi.

What I read: On develop the prim function set need_phi_derivative_sum = !is_alpha_vector (L127 at d12f488749).
At the head, neg_binomial_2_log_glm_lpmf.hpp L128-L129 sets it from is_phi_vector and computes the phi derivative only when phi is autodiff.

What I ran: python3 runTests.py -j8 --run-all test/unit/math/opencl/rev/neg_binomial_2_log_glm_lpmf_test.cpp, with the other OpenCL tests, at the head and on develop with the tests copied in.
Output in the gist: opencl_tests-head.out, opencl_tests_no_bb-base.out.

---

FACT: At small phi the partials were already correct on develop

Quote:

The partials were already correct, because they use the exact phi.

(PR description, Related defects fixed on the way, item 5)

STATUS: Confirmed

On develop with large_shapes_test.cpp copied in, the rows at phi = 1e-8 and 1e-17 of neg_binomial_lpmf, neg_binomial_2_lpmf and neg_binomial_2_log_lpmf fail only in the value.
Their gradients, compared on the log scale to 1e-11, pass on develop and at the head.

Evidence:

The test rows are in the PR.
Code: large_shapes_test.cpp, NB rows.

What I ran: python3 runTests.py -j8 --run-all test/unit/math/rev/prob/large_shapes_test.cpp at the head and on develop with the test copied in.
Output in the gist: cpu_tests-base.out, cpu_tests-head.out.

---

FACT: The GLM theta partial is NaN at theta = 720 and -inf at theta = 709 on develop, and the exact value is -1.5

Quote:

On develop the theta partial y - exp(theta) (y + phi) / (exp(theta) + phi) is NaN at theta = 720 (inf / inf) and -inf at theta = 709; the exact value is -1.5 at both points (y = 5 and 3, phi = 1.5).

(PR description, Related defects fixed on the way, item 6)

STATUS: Confirmed

The two theta rows of large_shapes_test.cpp check the alpha and beta partials, which both equal the theta partial at x = 1, against -1.5 to 1e-11.
On develop with the test copied in, both partials are NaN at theta = 720 and -inf at theta = 709, as jachymb's review comment reports.
At the head they pass, on the CPU and in the OpenCL reference test of the GLM.
On develop the OpenCL GLM reference test fails with the same exception as the mixed alpha and phi test, so it does not show these rows on develop.

Evidence:

The test rows are in the PR.
Code: large_shapes_test.cpp, theta rows.

What I ran: python3 runTests.py -j8 --run-all test/unit/math/rev/prob/large_shapes_test.cpp at the head and on develop with the test copied in, and the OpenCL GLM test at the head.
Output in the gist: cpu_tests-head.out, cpu_tests-base.out, opencl_tests-head.out, opencl_tests_no_bb-base.out.

---

FACT: A new test checks full - propto against the dropped terms

Quote:

a new test checks full - propto against the dropped terms

(PR description, Summary, propto)

STATUS: Confirmed

ProbDistributions_large_shapes_propto_dropped_terms checks full - propto to 1e-13 for student_t_lpdf (nu on both sides of 40), yule_simon_lpmf, beta_neg_binomial_lpmf and three cases of the GLM.
It passes at the head.
It also passes on develop with the test copied in, which matches the description's note that the propto case states the develop rule.

Evidence:

The test is in the PR.
Code: large_shapes_test.cpp, propto_dropped_terms.

What I ran: python3 runTests.py -j8 --run-all test/unit/math/rev/prob/large_shapes_test.cpp at the head and on develop with the test copied in.
Output in the gist: cpu_tests-head.out, cpu_tests-base.out.

---

FACT: With develop's headers, the new beta_binomial and large-shape cases fail, and the old cases and the propto case pass

Quote:

With the 23 develop headers swapped in (21 prim/prob files and the rev/fwd lbeta), all 4 new beta_binomial cases and both large-shape cases fail.
The 3 old cases and the propto case pass (the propto case states the develop rule).
With the branch headers restored, all tests pass again.

(PR description, Tests)

STATUS: Confirmed

On develop with the PR's test files copied in, which uses all of develop's headers, the 4 new beta_binomial cases fail (large_shapes, binomial_limit, log_shape_gradients and log_concentration_gradient).
Both large-shape cases fail (value_and_log_scale_gradients and lkj_corr_eta_gradient_at_one).
The 3 old cases of beta_binomial_test.cpp and the propto case pass.
At the head all of them pass.
The new GLM test NegBinomial2LogGLM_vector_matches_instances passes on develop too, which fits its purpose of comparing a vector call with single calls.

Evidence:

The tests are in the PR.
Code: beta_binomial_test.cpp, beta_binomial_lpmf_test.cpp, large_shapes_test.cpp.

What I ran: python3 runTests.py -j8 --run-all on the four CPU test files, at the head and on the merge base d12f488 with the test files copied in.
Output in the gist: cpu_tests-head.out, cpu_tests-base.out.

---

FACT: At the current head, the negative binomial unit tests, the generated test/prob tests (23 490 cases) and the four OpenCL negative binomial suites pass

Quote:

At the current head, after the merge of develop d12f488749, these pass: the 25 unit test files of the negative binomial family (including the Laplace tests that use neg_binomial_2_log and the new glm_forward_mode_tape_test from develop), the generated test/prob/neg_binomial and test/prob/neg_binomial_2 tests (23 490 cases), and the four OpenCL negative binomial suites (V100).

(PR description, Tests, Full suites)

STATUS: Confirmed

The PR does not list its 25 files, so I chose the 25 files of the negative binomial family by name.
At the head they run 150 tests in 25 binaries.
114 pass and none fail.
The other 36 are in the two Laplace marginal tests, which skip parameter combinations that the test utility marks as invalid (laplace_utility.hpp L310-L315).
The grep filter drops the SKIPPED lines, so this count is the difference between the tests run and the tests passed.
The generated tests of test/prob/neg_binomial and test/prob/neg_binomial_2 run 23 490 tests in 36 binaries and all pass.
The description gives the same count.
The four OpenCL negative binomial suites pass on the GPU of this machine (the description used a V100).

Evidence:

The test files are in the repository and the PR.
Code: test/prob/neg_binomial, test/prob/neg_binomial_2.

What I ran: python3 runTests.py -j8 --run-all on the 25 unit test files, then on test/prob/neg_binomial and test/prob/neg_binomial_2, at the head.
The first run printed more than the size limit of review.py, which cut the middle of the output, so I ran the same tests again with the gtest summary lines filtered by grep.
The OpenCL suites are in the OpenCL run at the head.
Output in the gist: nb_suites_summary-head.out, nb_suites-head.out, opencl_tests-head.out.

---

Partially Confirmed Facts

FACT: In the model of #3153, develop's value and lc gradient are wrong, and the PR is within 2.3e-14 and 3.7e-14 for lc from 0 to 44

Quote:

In the model of #3153, beta_binomial_lpmf(57 | 117, s, s) with s = exp(lc) / 2, develop returns a log probability 81 too high from lc = 42, and its gradient in lc is 0 or noise from lc = 24.
Warmup then finds a flat plateau and cannot leave it.
With this PR the value is within 2.3e-14 of the exact value and the gradient within 3.7e-14, for lc from 0 to 44.

(PR description, Summary)

STATUS: Partially Confirmed

The tests check this model at some of the points.
The prim test large_shapes checks the value at lc = 32, 36 and 42 to 1e-13 relative, and binomial_limit checks the difference to the binomial limit at lc = 42.
The rev test log_concentration_gradient checks the lc gradient at lc = 16, 18, 20, 26, 30 and 32 to 1e-12 absolute.
All of them pass at the head.

On develop with these tests copied in, they fail.
At lc = 42 the value is 78.451 against the reference -2.647, which is 81.098 too high, and at lc = 32 it is 4.47e-3 too high (the table gives +4.5e-3).
The lc gradient on develop is 1.72e-6 at lc = 20 and -0.140 at lc = 32, as in the develop column of the table (1.7e-06 and -0.14).

The bounds 2.3e-14 and 3.7e-14 over lc from 0 to 44, the row at lc = 28, and the other exact entries come from a sweep that is not in the PR.
The tests hold the exact d/dlc at lc = 20 and 32.

Evidence:

The tests have the points above, with mpmath reference values.
The sweep over lc from 0 to 44 that gives the table and the two bounds is not included.
Code: beta_binomial_test.cpp, large_shapes and binomial_limit, beta_binomial_lpmf_test.cpp, log_concentration_gradient.

What I ran: python3 runTests.py -j8 --run-all on the prim and rev beta_binomial tests, at the head (all pass) and on develop with the PR's tests copied in (the four new cases fail).
Output in the gist: cpu_tests-head.out, cpu_tests-base.out.

Request: Could you add or link the script behind the table and the bounds 2.3e-14 and 3.7e-14 for lc from 0 to 44 (the mpmath reference and the Stan Math driver)?
Then the rows that the tests do not cover can be rerun.

---

FACT: Fixes #1578: the error of alpha * d/dalpha of neg_binomial_lpmf reaches 3.0 on develop at alpha = 1e15 and is 1e-14 with this PR

Quote:

#1578 reports wrong derivatives of neg_binomial_lpmf for large alpha.
Its example (n = 0, alpha = 1e10 + 1, beta = 0.8) is already correct on develop, because develop no longer uses the Poisson approximation that the issue describes.
This PR repairs the loss that remains in d/dalpha = digamma(alpha + n) - digamma(alpha) + log(beta / (1 + beta)).
The loss matters when the mean alpha / beta stays fixed: the error of alpha * d/dalpha on develop reaches 3.0 at alpha = 1e15, and with this PR it is 1e-14 (rows of the new large-shape test).

(PR description, Summary)

STATUS: Partially Confirmed

The NB rows of large_shapes_test.cpp put alpha at 1e15 with the mean alpha / beta fixed at 3 and at 50.
The test compares alpha * d/dalpha to 1e-11 relative, or 1e-11 absolute where the value is below 1.
These rows pass at the head.
On develop they fail, and the largest error of a log-scale NB gradient is 3.0, as the description says.

The example in #1578 (n = 0, alpha = 1e10 + 1, beta = 0.8, where d/dalpha is about -0.81093 and d/dbeta about 6.94e9) is not a test row.
The Poisson branch that the issue describes is gone on develop.
neg_binomial_alpha_cutoff is still defined but no longer used.
So only the values at the issue's point (-0.81093 and 6.94e9) lack a test.
The figure 1e-14 for this PR is below the test tolerance, which shows an error below 1e-11.

Evidence:

The large-shape rows are in the PR.
The example from the issue and the measurement behind 1e-14 are not.
Code: large_shapes_test.cpp, NB rows, large_shapes_test.cpp, the comparison.

What I read: At d12f488, neg_binomial_lpmf.hpp L26 defines neg_binomial_alpha_cutoff = 1e10, and the rest of the file does not use it.

What I ran: python3 runTests.py -j8 --run-all test/unit/math/rev/prob/large_shapes_test.cpp, with the other tests, at the head and on develop with the test copied in.
Output in the gist: cpu_tests-head.out, cpu_tests-base.out.

Request: Could you add the example from #1578 as a test row?
The test would then show the issue's own case on develop and with this PR.
A script behind the 1e-14 figure would also let it be rerun.

---

FACT: The table of worst errors at large shapes, on develop and with this PR

Quote:

function	develop, worst error at large shape	this PR
beta_binomial_lpmf	value 81, log-shape gradient 6e3	rounding level
[...]
lbeta (var, fvar), derivative	log-scale derivative 0.63	1e-14

(PR description, Summary)

STATUS: Partially Confirmed

The rows of large_shapes_test.cpp and the beta_binomial tests are at these shapes.
On develop with the tests copied in, the largest errors at those rows compare with the table as follows:

• About the size in the table: beta_binomial_lpmf value 81, beta_neg_binomial_lpmf value 10.1 and gradient 12.5, neg_binomial_lpmf gradient 3.0, the GLM value 23.0 and gradient 7.1, and lkj_corr value 3.88 and gradient 1.5.
  The same holds for yule_simon_lccdf (value 7.6, gradient 2.1), the beta_binomial cdfs (gradient 22.1, lcdf NaN), the beta_neg_binomial cdfs (gradient 14.9), the student_t nu gradient (0.145) and lbeta (0.632).
• Smaller at the test rows than the worst case in the table: beta_binomial_lpmf log-shape gradient 1.76e3 (table 6e3), neg_binomial_2_lpmf gradient 3.6e-5 (table 27), neg_binomial_2_log_lpmf gradient 3.0 (table 26), dirichlet_multinomial gradient 1.16 (table 2.8), and the student_t value 0.052 (table 0.16).
• Without a row at a large shape: yule_simon_cdf and yule_simon_lcdf, which have one row each at alpha = 100.
• Larger than the table: multi_student_t and yule_simon_lpmf, in a separate fact check below.

At the head every row passes, to 1e-12 relative in the value and 1e-11 in the gradient (1e-6 and 1e-4 for the beta_binomial and beta_neg_binomial cdfs).
That bounds the column for this PR at the test tolerance.
The smaller entries of that column, such as 1e-14 and 2e-16, come from a sweep that is not in the PR.

Evidence:

The test rows, with mpmath reference values, are in the PR.
The script that found the worst errors over a range of shapes is not.
Code: large_shapes_test.cpp, beta_binomial_lpmf_test.cpp.

What I ran: python3 runTests.py -j8 --run-all on large_shapes_test.cpp and the beta_binomial tests, at the head and on develop with the tests copied in.
The develop errors above are the largest differences that gtest printed for each function.
Output in the gist: cpu_tests-head.out, cpu_tests-base.out.

Request: Could you add or link the script that measured the worst errors in this table, for develop and for this PR?
It would also help to add a row at a large shape for yule_simon_cdf and yule_simon_lcdf.

---

FACT: On develop, beta_binomial_lcdf returned NaN at 16 of 81 test points and was wrong by up to 61

Quote:

On develop this function returned NaN at 16 of 81 test points and was wrong by up to 61.
Now it has the accuracy of the series everywhere.

(PR description, Related defects fixed on the way, item 2)

STATUS: Partially Confirmed

The BBLC rows of large_shapes_test.cpp have four points.
On develop with the test copied in, the value is NaN at the three large-shape rows (alpha = beta = 1e4 and 1e6), and the log-shape gradients there are off by up to 70.
At the head the four rows pass to 1e-6 in the value and 1e-4 in the gradient.
The PR does not include the 81 test points behind the 16 NaN values and the error of 61, or a check of the accuracy everywhere.

Evidence:

Four test rows.
The 81 points are not in the PR.
Code: large_shapes_test.cpp, BBLC rows.

What I ran: python3 runTests.py -j8 --run-all test/unit/math/rev/prob/large_shapes_test.cpp at the head and on develop with the test copied in.
Output in the gist: cpu_tests-head.out, cpu_tests-base.out.

Request: Could you add the 81 points, or the script that checks them?
Then the count of NaN values on develop can be rerun.

---

FACT: OpenCL beta_binomial_lpmf did not check 0 <= n <= N, and returned -8.216 under propto for n = 6, N = 5

Quote:

The code built the check but never computed its result, so the function never returned LOG_ZERO.
Without propto the value was -inf by accident (through binomial_coefficient_log), but with a nonzero gradient.
Under propto it was finite: -8.216 for n = 6, N = 5 in the new test.

(PR description, Related defects fixed on the way, item 3)

STATUS: Partially Confirmed

The code part follows from reading.
On develop the kernel results do not include return_neg_inf (L96-L100 at d12f488749), so the check is never computed.
At the head it is computed.
opencl_n_outside_support checks LOG_ZERO and a zero gradient for n outside 0 to N, with and without propto, and it passes at the head on the GPU.

The red run cannot show the values on develop.
The test file test/unit/math/opencl/rev/beta_binomial_lpmf_test.cpp does not compile against develop, because the new reference test calls beta_binomial_lpmf with scalar var shapes (line 336).
Develop's OpenCL version then calls lbeta(const double&, const double&), which has no matching overload.
So no run reproduces the value -8.216 under propto or the nonzero gradient on develop.

Evidence:

The test opencl_n_outside_support is in the PR.
Code: opencl/rev/beta_binomial_lpmf_test.cpp, opencl_n_outside_support.

What I read: On develop the kernel results do not include return_neg_inf (beta_binomial_lpmf.hpp L96-L100 at d12f488749), so return_neg_inf_cl is never filled.
At the head, beta_binomial_lpmf.hpp L97-L106 computes it and returns LOG_ZERO.

What I ran: python3 runTests.py -j8 --run-all on the six OpenCL distribution tests, at the head (all pass) and on develop with the tests copied in, where the beta_binomial file fails to compile.
Output in the gist: opencl_tests-head.out, opencl_tests-base.out.

Request: Could you put the defect-3 test where it compiles against develop, for example in its own file or with the scalar-shape reference rows in a separate test file?
Then the red run can show the -8.216 under propto.

---

FACT: The negative binomial functions lost phi at small phi on develop, and the PR's errors at these points are at most 5.3e-13

Quote:

At n = 1000, phi = 1e-8, the error is 7.9e-7 (5.0e-9 in neg_binomial_lpmf).
At n = 1, phi = 1e-17, the value is -inf; the exact value is -39.1439.
These three functions and the GLM, CPU and OpenCL, now use -lbeta(phi, n + 1) - log(n + phi), and the errors at these points are at most 5.3e-13.

(PR description, Related defects fixed on the way, item 5)

STATUS: Partially Confirmed

The rows marked small phi and small alpha in large_shapes_test.cpp check these points for neg_binomial_lpmf, neg_binomial_2_lpmf, neg_binomial_2_log_lpmf and the GLM, to 1e-12 relative.
On develop with the test copied in, the value at phi = 1e-8 is off by 7.93e-7 for neg_binomial_2_lpmf and neg_binomial_2_log_lpmf and by 5.02e-9 for neg_binomial_lpmf, and at phi = 1e-17 it is -inf, as the description says.
The GLM rows pass on develop, which fits your reply to jachymb that the rounding in the GLM came from an earlier commit of this PR.
At the head all these rows pass.
In OpenCL, the reference tests of neg_binomial_lpmf, neg_binomial_2_lpmf and neg_binomial_2_log_lpmf fail on develop at phi = 1e-8 and 1e-17 and pass at the head.

The test tolerance is 2.5e-11 to 3.9e-11 at these values, so the tests do not show the bound of 5.3e-13.

Evidence:

The test rows are in the PR.
The measurement behind 5.3e-13 is not.
Code: large_shapes_test.cpp, NB rows.

What I read: The -lbeta(., n + 1) - log(n + .) form is in neg_binomial_lpmf.hpp L90-L96, neg_binomial_2_lpmf.hpp L74-L83, the GLM L223-L251 and opencl/prim/neg_binomial_2_lpmf.hpp L85-L100.

What I ran: python3 runTests.py -j8 --run-all test/unit/math/rev/prob/large_shapes_test.cpp at the head and on develop with the test copied in, and the same command on the OpenCL distribution tests.
Output in the gist: cpu_tests-head.out, cpu_tests-base.out, opencl_tests-head.out, opencl_tests_no_bb-base.out.

Request: Could you add the script behind the bound of 5.3e-13, or tighten the tolerance of these rows to that bound?

---

FACT: neg_binomial_2_lpmf evaluates log(mu + phi) once when phi is autodiff, and the results are bitwise identical

Quote:

When phi is autodiff, neg_binomial_2_lpmf evaluates log(mu + phi) once: the value and the phi partial both read it, and the lazy expression computed the log twice per element.
The results are bitwise identical.

(PR description, Changes per function)

STATUS: Partially Confirmed

The first part follows from the code.
At the head, log_mu_plus_phi is stored with to_ref_if<is_autodiff_v<T_precision>>, and the value and the phi partial both read it.
The PR does not include a comparison of the results before and after this change, so the bitwise identity cannot be rerun.

Evidence:

The code change is in commit 1c8afb7.
There is no comparison of results.
Code: neg_binomial_2_lpmf.hpp.

What I read: neg_binomial_2_lpmf.hpp L63-L65 stores the log, and L72 and L94 read it.

Request: Could you add or link the check that compared the results bit for bit, for example a test that runs both forms on the same inputs?

---

FACT: With develop's OpenCL files, every new absolute-reference test and the two tests of defects 3 and 4 fail, and the old differential tests pass except opencl_matches_cpu_big of the GLM

Quote:

In OpenCL, with the develop OpenCL files every new absolute-reference test fails (6 files), and so do the two new tests of related defects 3 and 4.
Every old differential test passes, including the propto comparisons of the GLM and student_t_lpdf, except opencl_matches_cpu_big of the GLM.

(PR description, Tests)

STATUS: Partially Confirmed

My base tree is all of develop at d12f488 with the PR's tests copied in, which differs from your swap of the develop OpenCL files into the branch.
In that tree the build of test/unit/math/opencl/rev/beta_binomial_lpmf_test.cpp fails, as described in the fact check of item 3.
So the reference test and the test of item 3 in that file do not run.
With the other five files:

• The absolute-reference tests of neg_binomial_lpmf, neg_binomial_2_lpmf, neg_binomial_2_log_lpmf and student_t_lpdf fail at the large-shape and small-phi rows and pass at the small-shape rows.
• opencl_matches_cpu_mixed_alpha_phi (item 4) fails with a size mismatch.
• opencl_matches_cpu_big and opencl_large_shapes_reference of the GLM then fail with the same exception, as you describe for opencl_matches_cpu_big.
  So the GLM reference test does not show its own rows on develop.
• Every other old differential test in these five files passes.

At the head all tests in the six files pass.

Evidence:

The six OpenCL test files are in the PR.
Code: test/unit/math/opencl/rev.

What I ran: python3 runTests.py -j8 --run-all on the six OpenCL distribution tests at the head, and on develop with the tests copied in, first with all six files (the build stops) and then without beta_binomial_lpmf_test.cpp.
Output in the gist: opencl_tests-head.out, opencl_tests-base.out, opencl_tests_no_bb-base.out.

Request: Could you say which files your swap replaced, or give a version of the beta_binomial OpenCL tests that compiles against develop?
Then the red run can cover the beta_binomial reference test and the test of item 3.

---

FACT: beta_binomial_lcdf(3 | 10, 298.3, 146.75) is NaN, with exact value -3.93335, for a different cause

Quote:

This PR does not close #2508, because the comment there about beta_binomial_lcdf(3 | 10, 298.3, 146.75) (NaN, exact -3.93335) has a different cause: hypergeometric_3F2 does not stop at the zero term of the terminating series 3F2({1, 302.3, -6}; {5, -151.75}; 1) (exact 17.63), and a later 0 * inf gives NaN.

(PR description, Summary)

STATUS: Partially Confirmed

The stated mechanism follows from the code of hypergeometric_3F2 at the head.
A zero parameter is replaced by 1 before the log, so the loop does not stop at the zero term, and t_sign becomes 0 after it and stays 0.
The PR does not include a test at this point, so the NaN and the exact value -3.93335 cannot be rerun here.
The comment in #2508 shows the lcdf only as a screenshot.

Evidence:

The PR does not include test at this point.
Code: hypergeometric_3F2.hpp.

What I read: hypergeometric_3F2.hpp L45-L50 replaces zero parameters by 1, so p is never -inf, and L62-L64 sets t_sign from the signs, which is 0 once a parameter is 0.

Request: Could you add this point as a test that records the current result, or link the test in #3437 that covers it?

---

Not Verified Facts

FACT: The references are mpmath at 80 to 140 digits, and two independent routes agree

Quote:

All references are mpmath at 80 to 140 digits, and two independent routes agree: loggamma for the values, closed-form digamma differences and mp.diff for the gradients, and sums of the pmf for the cdfs.

(PR description, Summary)

STATUS: Not Verified

Comments in the test files say that the reference values come from mpmath, but the PR does not include the mpmath scripts or the second route.
The reference values are in the tests, and the way they were produced cannot be rerun.

Evidence:

Comments in the tests.
The scripts are not in the PR.

Request: Could you add the mpmath scripts that produce the reference values in the tests, with both routes, or link them in a gist?

---

FACT: Over 7 281 points the error of digamma_diff is at most 4.2 ulp

Quote:

Over 7 281 points, x from 1e-3 to 1e18 and d from 1e-6 to 1e7, the error is at most 4.2 ulp; the plain difference reaches 6.5e16 ulp.

(PR description, What this PR adds)

STATUS: Not Verified

The prim test checks 17 mpmath values to 1e-15 relative, about 4.5 ulp, and passes at the head (cpu_tests-head.out).
The sweep over 7 281 points, and the script that measured the ulp errors of digamma_diff and of the plain difference, are not in the PR.

Evidence:

17 test values.
The sweep is not in the PR.

Request: Could you add or link the sweep script, with the reference values?

---

FACT: The student_t series has errors of at most 0.9 and 2.4 eps from h = 20 to 1e15

Quote:

Against mpmath, from h = 20 to 1e15, the errors are at most 0.9 and 2.4 eps.

(PR description, Changes per function)

STATUS: Not Verified

The ST rows of large_shapes_test.cpp and the OpenCL reference test check student_t_lpdf at four points to 1e-12 and 1e-11, and they pass at the head (cpu_tests-head.out and opencl_tests-head.out).
The comparison against mpmath from h = 20 to 1e15 that gives 0.9 and 2.4 eps is not in the PR.

Evidence:

Four test points.
The sweep is not in the PR.

Request: Could you add or link the script that compares the series with mpmath over this range?

---

FACT: Without the selection in log_beta_ratio, points with N = 1e6 and one small shape lost up to 4e-11

Quote:

Without the selection, points with N = 1e6 and one small shape lost up to 4e-11 where develop was correct.

(PR description, What this PR adds)

STATUS: Not Verified

The rev beta_binomial test has rows with N = 1e6 and one small shape, and they pass at the head to 1e-13 relative in the value (cpu_tests-head.out).
They test the code with the selection.
The run without the selection that lost 4e-11 is not in the PR.

Evidence:

The rows with the selection are in the PR.
The comparison without it is not.

Request: Could you add or link the comparison that shows the loss of 4e-11 without the selection?

---

FACT: The fvar Hessian entry of the GLM in theta is now -4.7225e-17 at theta = 35 and -1.40777e-13 at theta = 27

Quote:

At y = 1, phi = 0.07, the fvar<var> Hessian entry in theta is now -4.7225e-17 at theta = 35 (develop: 0) and -1.40777e-13 at theta = 27 (develop: -1.40970e-13).

(PR description, Related defects fixed on the way, item 6)

STATUS: Not Verified

The PR does not include a test or a script for the second derivatives.
The Hessian tests mentioned in the conversation are on jachymb's branch and are not part of this PR.

Evidence:

The PR does not include code for these numbers.

Request: Could you add a test of these two Hessian entries, or link the probe that computed them?

---

FACT: student_t_lpdf runs 2.4 % more instructions for nu < 40, 52 % fewer for nu >= 100, and the scalar value takes 0.35 times as long

Quote:

For nu < 40, student_t_lpdf executes the same operations as develop plus one comparison (+2.4 % instructions for value and gradient). From nu = 40 it executes fewer: -52 % instructions for value and gradient at nu >= 100, and the scalar value takes 0.35 times as long.

(PR description, Speed)

STATUS: Not Verified

The PR does not include the benchmark or the instruction counts behind these numbers.

Evidence:

The PR does not include benchmark.

Request: Could you add or link the benchmark, for example a google-benchmark file and the command that counted instructions?
Then it can be rerun at the head and on develop.

---

FACT: Time ratios to develop for the other changed functions

Quote:

Time ratio to develop for the other changed functions, at ordinary arguments (shapes below 10 unless stated).
Median of 5 interleaved rounds on one whole node (AMD EPYC 9965, g++ 11.5, -O2); the spread between rounds is at most 2.9 %, except 4.3 % for one develop row.

(PR description, Speed)

STATUS: Not Verified

The PR does not include the benchmark that produced the table.
The stan-buildbot end-to-end table in the conversation times whole models, and it does not isolate these kernels.

Evidence:

The PR does not include benchmark.

Request: Could you add or link the benchmark source and the script that interleaves the rounds?
Then the ratios can be rerun at the head and on develop, including the rows where the PR is slower, such as 1.45 and 1.33 for the GLM with propto.

---

FACT: A probe of the 14 changed densities matches the propto and full values of develop at all 456 points

Quote:

A probe of all 14 changed densities, with every combination of data and parameter arguments, matches the propto and full values of develop at all 456 points on the current head

(PR description, Summary, propto)

STATUS: Not Verified

The probe is not in the PR.
The new test of full - propto, in the Confirmed fact check "A new test checks full - propto against the dropped terms", covers four of the densities.

Evidence:

The probe is not in the PR.

Request: Could you add or link the probe?
Then the comparison with develop at the 456 points can be rerun.

---

FACT: For the GLM with data phi, the propto value has a rounding error of about eps * phi log(phi)

Quote:

For the GLM with data phi, this brings back a rounding error of about eps * phi log(phi) in the propto value.

(PR description, Summary, propto)

STATUS: Not Verified

The PR does not include a test or a script that shows the size of this error.

Evidence:

The PR does not include code for this statement.

Request: Could you add a short test or script that shows the size of this error at a large phi?

---

FACT: beta_neg_binomial_cdf now returns 0 with NaN gradients at 15 of 2352 points where develop returns finite values

Quote:

beta_neg_binomial_cdf then returns 0 with NaN gradients (15 of 2352 points of a wide grid), where develop returns a wrong finite value with finite gradients.

(PR description, Limits that remain)

STATUS: Not Verified

The PR does not include the grid of 2 352 points or a test at one of the 15 points, so this change from develop cannot be rerun.

Evidence:

The PR does not include code for this statement.

Request: Could you add one of the 15 points as a test, or link the grid script?
It would show reviewers the behavior that changes from develop.

---

Contradicted Facts

FACT: On develop, the worst error of multi_student_t and yule_simon_lpmf at large shapes is about 2e-7

Quote:

| multi_student_t_lpdf, multi_student_t_cholesky_lpdf | value 2e-7, nu gradient 2e-7 | 8e-15, 5e-15 |
| yule_simon_lpmf | value 2e-7, gradient 3e-7 | 3e-14, 2e-15 |

(PR description, Summary)

STATUS: Contradicted

The PR's own rows show larger errors on develop.
At nu = 1e12, the value of multi_student_t_cholesky_lpdf is off by 1.78e-3 and the value of multi_student_t_lpdf by 1.58e-3, and nu * d/dnu of the cholesky version is off by 8.9e-5.
For yule_simon_lpmf at alpha = 1e12, the value is off by 3.9e-3 and alpha * d/dalpha by 3.6e-3.
These shapes are within the range of at most 1e15 that the description gives for the last four rows, so a worst error of 2e-7 or 3e-7 on develop does not hold.
The rows pass at the head.

Evidence:

The rows MST, MSTF and YS of large_shapes_test.cpp.
Code: large_shapes_test.cpp, MST, MSTF and YS rows.

What I ran: python3 runTests.py -j8 --run-all test/unit/math/rev/prob/large_shapes_test.cpp on develop with the test copied in, and at the head.
The errors above are the differences that gtest printed.
Output in the gist: cpu_tests-base.out, cpu_tests-head.out.

Request: Could you check the develop column of these two rows?
With the PR's test, develop is off by about 2e-3 at nu = 1e12 and about 4e-3 at alpha = 1e12.

---

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: digamma_diff has three paths: a count from 0 to 8, the plain difference for x < 10 and d >= 10, and a shift with an asymptotic series otherwise

Quote:

otherwise: the recurrence shifts both arguments to at least 10, then log1p(d / y) plus the asymptotic series, with the differences y^-2i - (y + d)^-2i formed as (u - v) h_i so that nothing cancels.

(PR description, What this PR adds)

STATUS: Confirmed

digamma_diff.hpp L118-L126 sums 1 / (x + j) from the smallest term when d is arithmetic and an integer from 0 to 8.
L128-L130 returns the plain difference for x < 10 and d >= 10.
L145-L173 shifts y up to at least 10 and adds log1p(d / y) and the series with u - v and h_i as described.

---

FACT: digamma_diff is public, not in the Stan language, and used by 19 distribution files and the derivatives of lbeta

Quote:

It is public, but not exposed in the Stan language, because 19 distribution files and the derivatives of lbeta use it.

(PR description, What this PR adds)

STATUS: Confirmed

digamma_diff is in stan::math and is included from prim/fun.hpp.
It is called in 19 files under stan/math/prim/prob (multi_student_t_lpdf.hpp only names it in a comment and gets the derivative through lbeta), and in rev/fun/lbeta.hpp and fwd/fun/lbeta.hpp.
The PR does not change stanc3.

---

FACT: digamma_diff is vectorized, works with the autodiff types, and has an OpenCL device function and a kernel-generator registration

Quote:

The function is vectorized and works with all autodiff types through generic arithmetic, with expect_ad tests.
The PR also adds an OpenCL device function with the same paths and a kernel-generator registration.

(PR description, What this PR adds)

STATUS: Confirmed

digamma_diff.hpp L186-L194 vectorizes it with apply_scalar_binary, and the scalar version uses generic arithmetic on return_type_t.
The mix test calls expect_ad at 8 points and passes at the head (cpu_tests-head.out).
The device function in opencl/kernels/device_functions/digamma_diff.hpp has the same three paths, and elt_function_cl.hpp L417-L418 registers it.

---

FACT: log_beta_ratio selects the Stirling form or the plain lbeta form by the size of their terms

Quote:

log_beta_ratio(alpha, beta, n, m) = lbeta(alpha + n, beta + m) - lbeta(alpha, beta) selects this form or the plain lbeta form by the size of their terms.

(PR description, What this PR adds)

STATUS: Confirmed

log_rising_factorial_ratio.hpp L202-L239 uses the Stirling form when both shapes are large and size_stirling <= size_lbeta, and the plain lbeta difference otherwise.

---

FACT: The lbeta partial in an argument of at least 10 is -digamma_diff, and below 10 it is the plain difference

Quote:

The partial in an argument of at least 10 is now -digamma_diff(a, b).
Below 10 it stays digamma(a) - digamma(a + b), and the two partials share digamma(a + b) as before.

(PR description, What this PR adds)

STATUS: Confirmed

In rev/fun/lbeta.hpp L21-L45, lbeta_partial returns -digamma_diff(x, y) for x >= 10 and digamma(x) - digamma(x + y) otherwise.
The two partials share digamma(a + b) when both arguments are below 10.
fwd/fun/lbeta.hpp L22-L55 does the same.

---

FACT: The GLM evaluates the terms that depend on y only once per distinct count, if the largest count is at most 4 N + 64, in the CPU version only

Quote:

With a scalar phi and a vector y, the terms that depend on y only (lchoose, lgamma(y + 1) under propto, digamma_diff(phi, y)) are evaluated once per distinct count, through internal::per_count_values, if the largest count is at most 4 N + 64.
This is in the CPU version only.

(PR description, Changes per function)

STATUS: Confirmed

neg_binomial_2_log_glm_lpmf.hpp L211-L222 turns on the table for a vector y and a scalar phi when the largest count is at most 4 * N_instances + 64.
The table is used for the lchoose term (L226-L238), for lgamma(y + 1) under propto (L255-L268) and for digamma_diff (L339-L348).
The partials use e = exp(-|u|) and do not form exp(theta) (L186-L193, L288-L295).
The OpenCL kernel computes each instance separately and does not use a table.

---

FACT: student_t_lpdf is unchanged for nu < 40 and uses the stated series from nu = 40

Quote:

student_t_lpdf: unchanged for nu < 40. For nu >= 40, the terms lgamma((nu + 1) / 2) - lgamma(nu / 2) - log(nu) / 2 and digamma((nu + 1) / 2) - digamma(nu / 2) come from the asymptotic series in h = nu / 2, -log(2) / 2 - 1/(8h) + 1/(192h^3) - 1/(640h^5) + 17/(14336h^7) - 31/(18432h^9), and its derivative.

(PR description, Changes per function)

STATUS: Confirmed

student_t_lpdf.hpp L38-L101 uses the lgamma and digamma differences below 40 and the series with the coefficients in the description from 40 on.
Below 40 the same functions are called as on develop.
The value now sums each element's lgamma(h + 1/2) - lgamma(h) - log(nu) / 2 (L188-L193), where develop formed three separate sums.

---

FACT: Each term of the lkj constant is lgamma(k/2) - lbeta(k/2, x), and the eta partial is set with digamma_diff

Quote:

Each term is now lgamma(k/2) - lbeta(k/2, x), and the code sets the eta partial explicitly with digamma_diff.

(PR description, Changes per function)

STATUS: Confirmed

lkj_corr_lpdf.hpp L44-L62 adds lgamma(0.5 * k) - lbeta(0.5 * k, eta_val + 0.5 * (Km1 - k)) per term and sets the eta partial from a sum of digamma_diff values.

---

FACT: The changes per function for beta_binomial, beta_neg_binomial, neg_binomial, dirichlet_multinomial, yule_simon and multi_student_t

Quote:

beta_binomial_lpmf: value by log_beta_ratio, partials by digamma_diff.
[...]
multi_student_t_lpdf, multi_student_t_cholesky_lpdf: the constant lgamma((nu + p) / 2) - lgamma(nu / 2) is now lgamma(p / 2) - lbeta(p / 2, nu / 2).

(PR description, Changes per function)

STATUS: Confirmed

The code at the head matches each listed change, for example beta_binomial_lpmf.hpp L118-L141 (log_beta_ratio and digamma_diff), neg_binomial_lpmf.hpp L80-L105 (alpha / beta / (1 + beta) and -lbeta(alpha, n + 1) - log(n + alpha)), yule_simon_lpmf.hpp L66-L79 and multi_student_t_lpdf.hpp L105-L112.
The diff of the other listed files shows the same forms.

---

FACT: beta_binomial_lcdf uses the mirror relation where C > 1/2, except for alpha = 1

Quote:

Where C > 1/2, the function now uses the mirror relation: X <= n if and only if N - X > N - n - 1, and N - X has the distribution beta_binomial(N, beta, alpha).

(PR description, Related defects fixed on the way, item 2)

STATUS: Confirmed

beta_binomial_lcdf.hpp L118-L150 takes the mirrored series with the shapes swapped when C > 0.5 && alpha_dbl != 1.0, and otherwise keeps log(1 - C).

---

FACT: The OpenCL versions of six functions compute the same formulas as the CPU

Quote:

OpenCL: beta_binomial_lpmf, neg_binomial_lpmf, neg_binomial_2_lpmf, neg_binomial_2_log_lpmf, neg_binomial_2_log_glm_lpmf (prim file and kernel) and student_t_lpdf compute the same formulas.

(PR description, Changes per function)

STATUS: Confirmed

The OpenCL prim files of beta_binomial_lpmf, neg_binomial_lpmf, neg_binomial_2_lpmf and neg_binomial_2_log_lpmf use log_beta_ratio, -lbeta(., n + 1) - log(n + .) and digamma_diff as on the CPU (for example neg_binomial_2_lpmf.hpp L85-L100).
The GLM kernel uses the e form and stan_lbeta, and opencl/prim/student_t_lpdf.hpp L104-L153 switches to the series at the same threshold.

---

FACT: The code adds back the constants that propto drops in yule_simon_lpmf and the GLM, also in OpenCL

Quote:

To keep the propto behavior, the code adds such constants back with the opposite sign: lgamma(n) in yule_simon_lpmf, and lgamma(y + 1), phi log(phi) (data phi) and y theta (data x, alpha and beta) in neg_binomial_2_log_glm_lpmf, also in OpenCL.

(PR description, Summary, propto)

STATUS: Confirmed

yule_simon_lpmf.hpp L71-L74 subtracts lgamma(n) under propto.
neg_binomial_2_log_glm_lpmf.hpp L253-L282 adds lgamma(y + 1), subtracts phi log(phi) for data phi, and subtracts y theta for data x, alpha and beta.
The OpenCL prim file does the same with need_add_lgamma_y1, need_sub_phi_log_phi and need_sub_y_theta (L134-L146).

---

Partially Confirmed Facts

FACT: large_shapes_test.cpp has 92 rows, with points where develop was wrong plus one small-shape point per function

Quote:

test/unit/math/rev/prob/large_shapes_test.cpp (new): 92 rows over the other functions, including student_t_lpdf, both multi_student_t functions, yule_simon_lpmf and lbeta.
For each function, the rows are points where develop was wrong by more than 100 times the tolerance, plus one point at a small shape where develop was correct.

(PR description, Tests)

STATUS: Partially Confirmed

The file has 92 rows.
On develop with the file copied in (python3 runTests.py -j8 --run-all, output in cpu_tests-base.out), the last row of each function passes and the other rows fail for most functions.
Two groups differ.
The functions yule_simon_cdf and yule_simon_lcdf have one row each, at alpha = 100, so neither has a row where develop is wrong.
The GLM rows at phi = 1e-17 and 1e-8 pass on develop, which matches jachymb's table in the conversation, where develop gives the correct value at phi = 1e-17.

Request: Could you update the sentence to say that the GLM rows at small phi check a regression of an earlier commit?
It could also say that the yule_simon_cdf and _lcdf rows are small-shape rows only, or you could add rows at large alpha for those two functions.

---

Contradicted Facts

FACT: Each of the other 26 distribution files adds 3 to 61 lines

Quote:

Each of the other 26 distribution files (CPU, OpenCL and the GLM kernel) adds 3 to 61 lines.

(PR description, Summary)

STATUS: Contradicted

git diff --numstat from the merge base d12f488 to c962588 shows that stan/math/prim/prob/neg_binomial_2_log_glm_lpmf.hpp adds 154 lines, 107 of them in the last commit c962588.
That file has the fourth largest number of added source lines and is not in the table.
The other 25 distribution files add 3 to 56 lines.

Request: Could you add prim/prob/neg_binomial_2_log_glm_lpmf.hpp (154 added lines) to the table and update the range for the other files?

---

Not Checked

• The checklist items (the whole unit test suite, header checks, dependency checks, Doxygen and cpplint), which Jenkins runs.
• Statements about other PRs, such as Fix and speed up the terminating 3F2 series of the beta_binomial cdfs #3437 closing beta_binomial_lpmf is unstable for large shape arguments #2508 and Compute the beta_neg_binomial cdfs from sums of the pmf #3441 repairing the beta_neg_binomial cdfs.
• The items of the Limits section that do not describe a change made by this PR.
  These are the 3F2 tolerances of the beta_binomial and beta_neg_binomial cdfs, the throws from shape 1e4 and the time per call of the beta_neg_binomial cdfs, the floor of the log-shape gradients, the neg_binomial_2_log_lpmf Hessian, and log_rising_factorial.
• Two explanations of the method and the tests, without code to rerun.
  One says that develop returns exactly 0 from lc = 34, within 1e-13 of the exact gradient.
  The other says that the Stirling form has rounding errors of the size of n + m and the lbeta form of the size of the smaller shape.
• The sentence that the student_t series is cheaper than two lgamma or two digamma calls, which is the same claim as the student_t speed fact check.
• The stan-buildbot end-to-end table, which the PR does not cite as evidence.

[avehtari]

Reply to @stevieb-bot. I merged develop into the branch and added and moved tests. I put the scripts behind the numbers in a public repository, measured every number of the description again, and replaced the description. In this reply, "develop" is the new merge base 510d263078; your run used d12f488749. The branch head is now 35435f53dc, after the merge 17b0b961ec, the test commit 64407093ba and the fix commit 5005dcaa82.

Scripts

The scripts are in stan-pr-checks/pr3434. The README has one section per claim, with the commands, the expected output and the commits. The directory has:

• the mpmath reference generators, each with two routes and a print of their agreement (refs/);
• a check of every constant in the test files against these references (refs/check_test_refs.py);
• the sweeps behind the table of worst errors, the model of numerical instability in beta_binomial_lpmf at large shape parameters #3153 and digamma_diff (sweep/, dd/);
• probes for The derivatives of neg_binomial_lpmf are incorrect for large alpha #1578, small phi, the bitwise check of neg_binomial_2_lpmf, the student_t_lpdf series, log_beta_ratio without its selection, the propto checks and the beta_neg_binomial_cdf grid;
• the probes that I used to select the new test rows (review_rows/);
• the speed benchmarks with the job scripts (speed/), and the results of the new swap runs (results/swap_20261029.txt, README §"Swap tests").

I ran all accuracy scripts again on develop and on the PR (64407093ba; 5005dcaa82 changes only includes and comments in stan/). I also ran the speed benchmarks again on both trees (CSC Roihu job 2625261, AMD EPYC 9965, g++ 11.5, -O2).

Numbers that changed in the rerun

• The bound of the value in the model of numerical instability in beta_binomial_lpmf at large shape parameters #3153 is 4.3e-14 (at lc = 25.5), not 2.3e-14. The old text took 2.3e-14 from the single row lc = 42. The gradient bound 3.7e-14 holds.
• The table of worst errors: I scored every row again over all regions of its sweep. The old table used only some regions for several rows, and it gave absolute errors without the size of the value. The table now gives the relative error where the absolute error is large. As you found, the develop column of multi_student_t and yule_simon_lpmf was wrong. It did not include the region with shapes up to 1e12 and count 1, where develop is wrong by 1.8e-3 and 3.9e-3. Other rows changed for the same causes:
  • beta_binomial_lpmf, develop: value 220 and gradient 7.7e3 (before 81 and 6e3), at shapes up to 1e18. The old text said "shape 1e15", but the value 81 of numerical instability in beta_binomial_lpmf at large shape parameters #3153 is at shape 8.7e17. The text now gives the range of this row.
  • yule_simon_lpmf, PR: 6.9e-10 and 1.7e-10 at counts up to 1e6, where the value is up to 3.3e6 (relative 2e-16). Before: 3e-14 and 2e-15.
  • lkj_corr, PR: 0.15 and 0.14 absolute at a general matrix, relative 3e-16. The old numbers 4e-15 and 2e-16 were of the identity matrix only.
  • beta_neg_binomial_lpmf, PR: 2.6e-2 and 8.2e-2 absolute, relative 5e-17 and 3e-16; the old text said "rounding level".
  • GLM gradient, PR: 2.8e-14 (before 6e-14).
  • beta_binomial cdfs, PR: the gradients are at 2.3e-6, not "1e-7"; that number fits the values. Develop's lcdf gradient error is 4.1e4.
  • beta_neg_binomial_lccdf, PR: gradient 3.4e-6 (before 2e-6).
  • beta_neg_binomial_lcdf was not in the table. It has its own row now: the PR repairs it at large shapes, but not where the cdf is small, because of the complement 1 - ccdf ("Limits that remain").
• beta_binomial_lcdf on develop: NaN at 16 of the 81 points (as before), and wrong by up to 74 at the others, not 61. 61 is the second-largest error; I could not find why the old text used it.
• The student_t_lpdf series: at most 0.51 and 1.76 eps, not 0.9 and 2.4. The old numbers were of an earlier version of the series code.
• log_beta_ratio without the selection: up to 6.1e-11 relative and 2.8e-10 absolute, at 131 of 13 177 points with N = 1000 and 1e6. The old text gave 4e-11 and "one small shape".
• "Limits that remain" said that the lc gradient of numerical instability in beta_binomial_lpmf at large shape parameters #3153 is 27 / exp(lc) at tied shapes. It is 27 / s = 54 / exp(lc), as the table in the description shows; the test comment had the same error.
• The references of student_t_lpdf, multi_student_t, yule_simon_lpmf and lbeta are at 50 digits, checked at 110, not at "80 to 140".
• Speed: below nu = 40, value and gradient of student_t_lpdf execute 3.7 % more instructions (before 2.4 %). At nu >= 100, value and gradient execute 54 % fewer instructions (before 52 %). The value takes 0.33 times as long as on develop (before 0.35). The other ratios moved by at most 0.03. The spread between rounds is at most 2.8 % for develop and 6.9 % for the PR (before "2.9 %, except 4.3 %").

Partially confirmed and not verified facts

• The model of numerical instability in beta_binomial_lpmf at large shape parameters #3153: sweep/score_bb.py, region REP (lc from 0 to 44 in steps of 0.25), gives the table and the bounds. The value bound changed (above).
• The derivatives of neg_binomial_lpmf are incorrect for large alpha #1578: the example of the issue is now a row of large_shapes_test.cpp. It passes on develop too, as you found, and the test file says so. The 1.0e-14 is the worst of 42 points with the mean fixed at 3 or 50 (issue1578/, the NB rows of the sweep).
• The table of worst errors: see "Numbers that changed". I added yule_simon_cdf and yule_simon_lcdf rows at alpha = 1e8 and 1e12. These two functions are 1 - c and log1p(-c) with a tiny c = n B(n, alpha + 1). With the test tolerance tol * max(1, |reference|), no point at a large shape could fail on develop. So the test now compares these two tags relative to the reference. Develop's relative gradient error is 4e-3 at alpha = 1e12 and 1.0 at 1e15 (n = 1), and the new rows fail on develop.
• beta_binomial_lcdf, 16 of 81: the 81 points are the BBLC lines of sweep/points_sib.txt. The count holds, and the largest finite error is 74 (above).
• beta_binomial_lpmf in OpenCL (item 3): the scalar-shape reference rows are now in beta_binomial_lpmf_scalar_shapes_test.cpp, and the shared rows in beta_binomial_lpmf_reference.hpp. beta_binomial_lpmf_test.cpp now compiles against develop. There, opencl_n_outside_support fails with lp_propto.val() = -8.21632 and alpha adjoints of 1.86, 1.49 and -0.95, and opencl_large_shapes_reference fails too. The scalar-shape file does not compile against develop. That is a develop defect, which the description now names as item 4.
• Small phi, 5.3e-13: small_phi/ gives 5.3e-13 as the worst over the four functions, at the two points of item 5 and at n = 100, phi = 1e-6. I did not tighten the test tolerance.
• log(mu + phi) once: nb2_bitwise/nb2_bitwise.cpp prints the value and every adjoint for 2 000 instances and five call types. The output of the commits before and after the change is bitwise identical.
• The OpenCL swap: my new run used the develop versions of all 32 changed files that exist on develop, 24 CPU and 8 OpenCL files, as your run did. The four new files of the PR stay, because develop does not have them; the develop files do not include them. Every new reference test fails, and so does opencl_n_outside_support. The README section "Swap tests" lists the files and the results.
• beta_binomial_lcdf(3 | 10, 298.3, 146.75): the test of #3437 has this point, with the value -3.93335 and the shape gradients. A test in this PR would only record a NaN, so I linked that test.
• The GLM fvar<var> Hessian entries: a new test in mix/prob/neg_binomial_2_log_glm_lpmf_test.cpp checks both entries against mpmath (closed form and mp.diff). On develop the entry is 0 at theta = 35 and has a relative error of 1.4e-3 at 27, so the test fails.
• log_beta_ratio: lbr_select/ forces the Stirling form in the PR tree and scores it on the 13 177 points (numbers above).
• The GLM propto value with data phi: propto/glm_propto.cpp. At phi = 1e12 the error is 5.2e-3, against eps |propto value| = 3.1e-2, and develop's error has the same size. The description now gives these numbers.
• beta_neg_binomial_cdf, 15 of 2 352: bnb_grid/ finds the same 15 points as before. In this run, 806 points did not finish within 2 s per call on either tree, so the count covers the other 1 546 points. The description now says so.
• The facts whose only request was a script:
  • the references: refs/;
  • digamma_diff over 7 281 points: dd/, worst 4.24 ulp, plain difference 6.5e16 ulp;
  • the student_t_lpdf series: student_t_series/;
  • the propto probe: propto/, 0 mismatches in 456 values against develop;
  • the speed: speed/, rerun above.

Contradicted fact

• multi_student_t and yule_simon_lpmf, develop column: you are right. See "Numbers that changed".

Description vs code

• "92 rows": the file now has 97 rows. The description says which rows pass on develop by design. These are the GLM rows at small phi, which check a regression of an earlier commit of this PR, and the neg_binomial_lpmf row at the point of The derivatives of neg_binomial_lpmf are incorrect for large alpha #1578. It also describes the new yule_simon_cdf and _lcdf rows.
• "3 to 61 lines": the table now has prim/prob/neg_binomial_2_log_glm_lpmf.hpp, which adds 155 lines after the merge. The other 25 distribution files add 3 to 56 lines, and the rev and fwd lbeta files 28 and 27.

Other changes

• I merged develop 510d263078 into the branch. The conflicts came from Fix OpenCL neg_binomial_2_log_glm_lpmf with vector phi and scalar alpha #3428 and cleanup refactor: Use lmultiply throughout and make multiply_log an alias of it #3485. Fix OpenCL neg_binomial_2_log_glm_lpmf with vector phi and scalar alpha #3428 (@jachymb) repaired the OpenCL GLM defect that was item 4 of the old description, with a test of both mixed cases at N = 3. (The new item 4 is the scalar-shape compile error.) I removed the item from the description. The test opencl_matches_cpu_mixed_alpha_phi of this PR now runs only N = 153, which needs more than one work group, with and without propto. cleanup refactor: Use lmultiply throughout and make multiply_log an alias of it #3485 replaced multiply_log with lmultiply; the PR code now uses lmultiply.
• 35435f53dc adds a test of the var and fvar derivatives of lbeta at arguments from 1e8 to 1e14 (mix/fun/lbeta_test.cpp). The PR changes these derivatives, and no test checked them at large arguments before. With the develop rev and fwd lbeta files the test fails, and with the PR files it passes (Triton job 21278515).
• 5005dcaa82 removes includes that develop had removed as unused and that the merge kept. It also rewrites test and code comments that compared the code with develop, so that they say what the code computes.
• I tested the tree of 5005dcaa82 with CPU builds on Triton batch nodes (g++ 14.3.0). 56 unit test files of the touched functions pass. The generated test/prob tests of neg_binomial and neg_binomial_2 pass (23 490 cases). At 64407093ba, the generated tests of all six families passed (271 306 cases). With the 24 develop CPU files, 7 of the 9 new test cases fail; the other 2 (large_shapes_propto_dropped_terms and the GLM vector test) are meant to pass on develop. The 9 OpenCL suites pass on a V100 (115 cases). make cpplint, make test-math-dependencies and make test-headers pass, and Doxygen 1.9.7 passed at 64407093ba.