From ad3e0d72131bda0f4c17d7ab74a2303fe7aacca5 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Fri, 9 Oct 2026 12:07:51 +0300 Subject: [PATCH 1/5] Stop hypergeometric_3F2 and grad_F32 at the end of a terminating series --- stan/math/prim/err/check_3F2_converges.hpp | 10 +- stan/math/prim/fun/grad_F32.hpp | 6 + stan/math/prim/fun/hypergeometric_3F2.hpp | 12 +- test/unit/math/prim/fun/grad_F32_test.cpp | 44 ++++++- .../math/prim/fun/hypergeometric_3F2_test.cpp | 42 +++++- .../prim/prob/beta_binomial_cdf_log_test.cpp | 12 +- .../math/rev/prob/beta_binomial_cdfs_test.cpp | 121 ++++++++++++++++++ 7 files changed, 233 insertions(+), 14 deletions(-) create mode 100644 test/unit/math/rev/prob/beta_binomial_cdfs_test.cpp diff --git a/stan/math/prim/err/check_3F2_converges.hpp b/stan/math/prim/err/check_3F2_converges.hpp index 25d9b89afe8..98260c9c133 100644 --- a/stan/math/prim/err/check_3F2_converges.hpp +++ b/stan/math/prim/err/check_3F2_converges.hpp @@ -66,8 +66,14 @@ inline void check_3F2_converges(const char* function, const T_a1& a1, num_terms = floor(fabs(value_of_rec(a3))); } - bool is_undefined = (is_nonpositive_integer(b1) && fabs(b1) <= num_terms) - || (is_nonpositive_integer(b2) && fabs(b2) <= num_terms); + // A denominator parameter b = -m (m a nonnegative integer) is a pole from + // term m + 1 on, because (b)_k is zero for k > m. A polynomial ends at + // term num_terms, so there b is a pole only if m < num_terms. + auto is_pole = [&](const auto& b) { + return is_nonpositive_integer(b) + && (is_polynomial ? fabs(b) < num_terms : fabs(b) <= num_terms); + }; + bool is_undefined = is_pole(b1) || is_pole(b2); if (is_polynomial && !is_undefined) { return; diff --git a/stan/math/prim/fun/grad_F32.hpp b/stan/math/prim/fun/grad_F32.hpp index c8f2d56d5f7..e9d7536eb8e 100644 --- a/stan/math/prim/fun/grad_F32.hpp +++ b/stan/math/prim/fun/grad_F32.hpp @@ -79,6 +79,12 @@ inline void grad_F32(T1* g, const T2& a1, const T3& a2, const T4& a3, } std::array term{0}; for (int k = 0; k <= max_steps; ++k) { + // A numerator parameter that has reached zero ends the series (a + // polynomial). Stop before the ratio, which is 0 / 0 when a denominator + // parameter reaches zero at the same k. + if (a1 + k == 0 || a2 + k == 0 || a3 + k == 0) { + return; + } T1 p = (a1 + k) * (a2 + k) * (a3 + k) / ((b1 + k) * (b2 + k) * (1 + k)); if (p == 0) { return; diff --git a/stan/math/prim/fun/hypergeometric_3F2.hpp b/stan/math/prim/fun/hypergeometric_3F2.hpp index 0cff50e8daa..a2b2e5f61ed 100644 --- a/stan/math/prim/fun/hypergeometric_3F2.hpp +++ b/stan/math/prim/fun/hypergeometric_3F2.hpp @@ -40,6 +40,13 @@ inline return_type_t hypergeometric_3F2_infsum( int k = 0; const double log_precision = log(precision); while (k <= max_steps && log_t >= log_precision) { + // A numerator parameter that has reached zero makes this term and every + // later term zero: the series is a polynomial and has ended. Without + // this stop the sign below is 0 while the magnitude keeps growing, and + // 0 * inf gives NaN. + if ((value_of_rec(a_array) == 0.0).any()) { + return t_acc; + } // Replace zero values with 1 prior to taking the log so that we accumulate // 0.0 rather than -inf const auto& abs_apk = math::fabs((a_array == 0).select(1.0, a_array)); @@ -123,8 +130,9 @@ inline auto hypergeometric_3F2(const Ta& a, const Tb& b, const Tz& z) { check_3F2_converges("hypergeometric_3F2", a_ref[0], a_ref[1], a_ref[2], b_ref[0], b_ref[1], z); // Boost's pFq throws convergence errors in some cases, fallback to naive - // infinite-sum approach (tests pass for these) - if (z == 1.0 && (sum(b_ref) - sum(a_ref)) < 0.0) { + // infinite-sum approach (tests pass for these). At z = 1 Boost also throws + // when sum(b) == sum(a), also for a terminating series. + if (z == 1.0 && (sum(b_ref) - sum(a_ref)) <= 0.0) { return internal::hypergeometric_3F2_infsum(a_ref, b_ref, z); } return hypergeometric_pFq(a_ref, b_ref, z); diff --git a/test/unit/math/prim/fun/grad_F32_test.cpp b/test/unit/math/prim/fun/grad_F32_test.cpp index 40ef862a983..f91aa0cd725 100644 --- a/test/unit/math/prim/fun/grad_F32_test.cpp +++ b/test/unit/math/prim/fun/grad_F32_test.cpp @@ -44,10 +44,22 @@ TEST(MathPrimScalFun, grad_F32_converges_by_z) { // // } +// terminate by zero numerator at k = 1; the denominator (b2)_k is zero only +// from k = 2 on, so the polynomial is defined. The derivative with respect +// to a3 is not checked: a3 is the integer that ends the series. +TEST(MathPrimScalFun, grad_F32_short_polynomial_denominator_equal_numerator) { + std::vector g = {1.2, 0.1, 0.0, -0.12, 1.2, 1.2}; + double g_calc[6]; + stan::math::grad_F32(g_calc, 1.0, 12.0, -1.0, 10.0, -1.0, 1.0); + for (int i : {0, 1, 3, 4, 5}) { + EXPECT_NEAR(g[i], g_calc[i], 1e-14); + } +} + // at pole, should throw TEST(MathPrimScalFun, grad_F32_short_polynomial_undef) { double g_calc[6]; - EXPECT_THROW(stan::math::grad_F32(g_calc, 1.0, 12.0, -1.0, 10.0, -1.0, 1.0), + EXPECT_THROW(stan::math::grad_F32(g_calc, 1.0, 12.0, -2.0, 10.0, -1.0, 1.0), std::domain_error); } @@ -106,6 +118,36 @@ TEST(MathPrimScalFun, grad_F32_double_sign_flip_2) { EXPECT_NEAR(g[i], g_calc[i], 1e-8); } +// The tests below use z = 1 and terminate by zero numerator at k = 6. The +// reference values are the derivatives of the exact finite sums. The +// derivative with respect to a3 is not checked: a3 is the integer that ends +// the series. + +// b2 = a3, so the denominator (b2)_k is zero only from k = 7 on +TEST(MathPrimScalFun, grad_F32_polynomial_denominator_equal_numerator) { + std::vector g + = {26.180979740203373, 6.4919047619047619, 0.0, + -8.0448964404555225, 15.921514834449405, 52.373725818452381}; + double g_calc[6]; + stan::math::grad_F32(g_calc, 1.0, 6.5, -6.0, 5.0, -6.0, 1.0); + for (int i : {0, 1, 3, 4, 5}) { + EXPECT_NEAR(g[i], g_calc[i], 1e-12 * std::fabs(g[i])); + } +} + +// a2 and b2 are large +TEST(MathPrimScalFun, grad_F32_polynomial_large_parameters) { + std::vector g + = {29.010777844568475, 0.16148712923118207, 0.0, + -8.114218400390939, 0.32571439023403243, 49.018778179540667}; + double g_calc[6]; + stan::math::grad_F32(g_calc, 1.0, 302.30333970, -6.0, 5.0, -151.75306521, + 1.0); + for (int i : {0, 1, 3, 4, 5}) { + EXPECT_NEAR(g[i], g_calc[i], 1e-12 * std::fabs(g[i])); + } +} + // // m = { // {1.0, 1.0, 1.0, 1.0, 1.0, 0.6, 1e-10}, diff --git a/test/unit/math/prim/fun/hypergeometric_3F2_test.cpp b/test/unit/math/prim/fun/hypergeometric_3F2_test.cpp index 1ca327223d5..087027e6a01 100644 --- a/test/unit/math/prim/fun/hypergeometric_3F2_test.cpp +++ b/test/unit/math/prim/fun/hypergeometric_3F2_test.cpp @@ -28,10 +28,17 @@ TEST(MathPrimScalFun, F32_short_polynomial) { -0.08000000000000007, stan::math::hypergeometric_3F2({1.0, 12.0, -1.0}, {10.0, 1.0}, .9), 1e-8); } +// terminate by zero numerator at k = 1; the denominator (b2)_k is zero only +// from k = 2 on, so the polynomial is defined +TEST(MathPrimScalFun, F32_short_polynomial_denominator_equal_numerator) { + EXPECT_NEAR( + 2.2, stan::math::hypergeometric_3F2({1.0, 12.0, -1.0}, {10.0, -1.0}, 1.0), + 1e-14); +} // at pole, should throw TEST(MathPrimScalFun, F32_short_polynomial_undef) { EXPECT_THROW( - stan::math::hypergeometric_3F2({1.0, 12.0, -1.0}, {10.0, -1.0}, 1.0), + stan::math::hypergeometric_3F2({1.0, 12.0, -2.0}, {10.0, -1.0}, 1.0), std::domain_error); } // converge, single sign flip via numerator @@ -56,3 +63,36 @@ TEST(MathPrimScalFun, F32_double_sign_flip) { 1.06593846110441323674, stan::math::hypergeometric_3F2({1.0, -.5, -4.5}, {10.0, 1.0}, 0.3), 1e-8); } + +// The tests below use z = 1 and sum(b) <= sum(a). The reference values are +// the exact finite sums. + +// terminate by zero numerator at k = 6, with large a2 and b2 +TEST(MathPrimScalFun, F32_polynomial_large_parameters) { + EXPECT_NEAR(17.62967394122564, + stan::math::hypergeometric_3F2({1.0, 302.30333970, -6.0}, + {5.0, -151.75306521}, 1.0), + 1e-12); +} +// a numerator parameter equal to zero: only the first term is not zero +TEST(MathPrimScalFun, F32_zero_numerator) { + EXPECT_EQ(1.0, stan::math::hypergeometric_3F2({1.0, 302.30333970, 0.0}, + {5.0, -151.75306521}, 1.0)); + EXPECT_EQ(1.0, + stan::math::hypergeometric_3F2({1.0, 101.0, 0.0}, {2.0, 0.5}, 1.0)); +} +// terminate by zero numerator at k = 6; b2 = a3, so the denominator (b2)_k +// is zero only from k = 7 on +TEST(MathPrimScalFun, F32_polynomial_denominator_equal_numerator) { + EXPECT_NEAR( + 14.118912760416667, + stan::math::hypergeometric_3F2({1.0, 6.5, -6.0}, {5.0, -6.0}, 1.0), + 1e-12); +} +// terminate by zero numerator at k = 6, with sum(b) == sum(a) +TEST(MathPrimScalFun, F32_polynomial_equal_parameter_sums) { + EXPECT_NEAR( + 9.5967841682127396, + stan::math::hypergeometric_3F2({1.0, 4.5, -6.0}, {5.0, -5.5}, 1.0), + 1e-12); +} diff --git a/test/unit/math/prim/prob/beta_binomial_cdf_log_test.cpp b/test/unit/math/prim/prob/beta_binomial_cdf_log_test.cpp index e5cf707dbf7..b1883deaf53 100644 --- a/test/unit/math/prim/prob/beta_binomial_cdf_log_test.cpp +++ b/test/unit/math/prim/prob/beta_binomial_cdf_log_test.cpp @@ -35,12 +35,8 @@ TEST(ProbBetaBinomial, lcdf_matches_mathematica) { double alpha = 3.0; double beta = 1.0; - // EXPECT_NEAR( - // -0.5500463, - // (stan::math::beta_binomial_lcdf(n, N, alpha, beta)), 1e-8); - // FIXME: this point _should_ be defined for the beta_binomial_lcdf - // to be defined over its full parameter range but the power-series - // is not defined. Leaving the test in place with the current behavior. - EXPECT_THROW(stan::math::beta_binomial_lcdf(n, N, alpha, beta), - std::domain_error); + // log(15 / 26), computed to more digits than the Mathematica value + // -0.5500463 + EXPECT_NEAR(-0.55004633691927198, + stan::math::beta_binomial_lcdf(n, N, alpha, beta), 1e-12); } diff --git a/test/unit/math/rev/prob/beta_binomial_cdfs_test.cpp b/test/unit/math/rev/prob/beta_binomial_cdfs_test.cpp new file mode 100644 index 00000000000..1074e6584f8 --- /dev/null +++ b/test/unit/math/rev/prob/beta_binomial_cdfs_test.cpp @@ -0,0 +1,121 @@ +#include +#include +#include +#include +#include + +namespace { + +// Values and shape gradients of beta_binomial_cdf, beta_binomial_lcdf and +// beta_binomial_lccdf at points where the 3F2 series of these functions +// ends at a zero numerator parameter. Reference values: sums of the pmf and +// of its closed-form shape derivatives (mpmath, 50 digits). +struct beta_binomial_cdfs_point { + int n; + int N; + double alpha; + double beta; + // rows: cdf, lcdf, lccdf; columns: value, d/dalpha, d/dbeta + double ref[3][3]; +}; + +const std::vector beta_binomial_cdfs_points = { + // large shapes: the 3F2 terms grow before the series ends + {3, + 10, + 298.30333970, + 146.75306521, + {{0.019577975830469777, -0.00025213357314530317, 0.0005048815375647927}, + {-3.9333500266686489, -0.012878429074005716, 0.025788239904711231}, + {-0.019772163107297, 0.00025716841006185449, -0.00051496348013239942}}}, + // beta = 1: a denominator parameter of 3F2 equals the numerator + // parameter that ends the series + {3, + 10, + 2.5, + 1.0, + {{0.10178821044348033, -0.078609438975422664, 0.12549085581810137}, + {-2.2848609925400993, -0.77228432087497911, 1.232862384271726}, + {-0.1073493926558798, 0.08751771006505609, -0.13971187784126151}}}, + {8, + 10, + 3.0, + 1.0, + {{0.57692307692307692, -0.092455621301775148, 0.3774924861463323}, + {-0.55004633691927198, -0.16025641025641026, 0.65432030932030932}, + {-0.8602012652231115, 0.21853146853146853, -0.89225496725496725}}}, + // alpha + beta = 1: the sums of the 3F2 numerator and denominator + // parameters are equal + {3, + 10, + 0.5, + 0.5, + {{0.4080352783203125, -0.61252655150398375, 0.49626754881843688}, + {-0.89640164184463105, -1.5011607673372379, 1.2162368677074559}, + {-0.52430823763107531, 1.0347348905623168, -0.83833973654762416}}}, + // N - n = 1: a 3F2 numerator parameter is zero + {0, + 1, + 100.0, + 0.5, + {{0.0049751243781094527, -0.000049503725155317938, 0.0099007450310635875}, + {-5.3033049080590758, -0.0099502487562189055, 1.9900497512437811}, + {-0.0049875415110390736, 0.000049751243781094527, + -0.0099502487562189055}}}, +}; + +template +void expect_beta_binomial_cdfs(const char* name, int row, const F& f) { + using stan::math::var; + const double rel_tol = 1e-12; + for (const auto& p : beta_binomial_cdfs_points) { + const double* ref = p.ref[row]; + const double value_d = f(p.n, p.N, p.alpha, p.beta); + EXPECT_NEAR(ref[0], value_d, rel_tol * std::fabs(ref[0])) + << name << "(" << p.n << " | " << p.N << ", " << p.alpha << ", " + << p.beta << "), double arguments"; + + var alpha = p.alpha; + var beta = p.beta; + var value = f(p.n, p.N, alpha, beta); + std::vector vars = {alpha, beta}; + std::vector grad; + value.grad(vars, grad); + EXPECT_NEAR(ref[0], value.val(), rel_tol * std::fabs(ref[0])) + << name << "(" << p.n << " | " << p.N << ", " << p.alpha << ", " + << p.beta << ")"; + EXPECT_NEAR(ref[1], grad[0], rel_tol * std::fabs(ref[1])) + << "d/dalpha " << name << "(" << p.n << " | " << p.N << ", " << p.alpha + << ", " << p.beta << ")"; + EXPECT_NEAR(ref[2], grad[1], rel_tol * std::fabs(ref[2])) + << "d/dbeta " << name << "(" << p.n << " | " << p.N << ", " << p.alpha + << ", " << p.beta << ")"; + stan::math::recover_memory(); + } +} + +} // namespace + +TEST_F(AgradRev, ProbDistributionsBetaBinomial_cdf_terminating_series) { + expect_beta_binomial_cdfs( + "beta_binomial_cdf", 0, + [](int n, int N, const auto& alpha, const auto& beta) { + return stan::math::beta_binomial_cdf(n, N, alpha, beta); + }); +} + +TEST_F(AgradRev, ProbDistributionsBetaBinomial_lcdf_terminating_series) { + expect_beta_binomial_cdfs( + "beta_binomial_lcdf", 1, + [](int n, int N, const auto& alpha, const auto& beta) { + return stan::math::beta_binomial_lcdf(n, N, alpha, beta); + }); +} + +TEST_F(AgradRev, ProbDistributionsBetaBinomial_lccdf_terminating_series) { + expect_beta_binomial_cdfs( + "beta_binomial_lccdf", 2, + [](int n, int N, const auto& alpha, const auto& beta) { + return stan::math::beta_binomial_lccdf(n, N, alpha, beta); + }); +} From 4e31f27804c5f34452a09d55ab5e3a805171f224 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Fri, 9 Oct 2026 12:08:07 +0300 Subject: [PATCH 2/5] Fix the pole check, the step limit and the sign of z in the 3F2 series --- stan/math/prim/err/check_3F2_converges.hpp | 33 ++++++++++--------- stan/math/prim/fun/hypergeometric_3F2.hpp | 5 +-- .../prim/err/check_3F2_converges_test.cpp | 27 +++++++++++++++ test/unit/math/prim/fun/grad_F32_test.cpp | 13 ++++++++ .../math/prim/fun/hypergeometric_3F2_test.cpp | 28 ++++++++++++++++ 5 files changed, 88 insertions(+), 18 deletions(-) diff --git a/stan/math/prim/err/check_3F2_converges.hpp b/stan/math/prim/err/check_3F2_converges.hpp index 98260c9c133..8c74f7766f1 100644 --- a/stan/math/prim/err/check_3F2_converges.hpp +++ b/stan/math/prim/err/check_3F2_converges.hpp @@ -50,28 +50,29 @@ inline void check_3F2_converges(const char* function, const T_a1& a1, check_not_nan("check_3F2_converges", "b2", b2); check_not_nan("check_3F2_converges", "z", z); - int num_terms = 0; + // A numerator parameter a = -m (m a nonnegative integer) makes the terms + // zero from term m + 1 on, so the series is a polynomial that ends at the + // smallest such m. bool is_polynomial = false; - - if (is_nonpositive_integer(a1) && fabs(a1) >= num_terms) { - is_polynomial = true; - num_terms = floor(fabs(value_of_rec(a1))); - } - if (is_nonpositive_integer(a2) && fabs(a2) >= num_terms) { - is_polynomial = true; - num_terms = floor(fabs(value_of_rec(a2))); - } - if (is_nonpositive_integer(a3) && fabs(a3) >= num_terms) { - is_polynomial = true; - num_terms = floor(fabs(value_of_rec(a3))); - } + double num_terms = 0; + auto add_end = [&](const auto& a) { + if (is_nonpositive_integer(a)) { + const double m = fabs(value_of_rec(a)); + num_terms = is_polynomial ? std::fmin(num_terms, m) : m; + is_polynomial = true; + } + }; + add_end(a1); + add_end(a2); + add_end(a3); // A denominator parameter b = -m (m a nonnegative integer) is a pole from // term m + 1 on, because (b)_k is zero for k > m. A polynomial ends at - // term num_terms, so there b is a pole only if m < num_terms. + // term num_terms, so there b is a pole only if m < num_terms; in an + // infinite series b is always a pole. auto is_pole = [&](const auto& b) { return is_nonpositive_integer(b) - && (is_polynomial ? fabs(b) < num_terms : fabs(b) <= num_terms); + && (!is_polynomial || fabs(value_of_rec(b)) < num_terms); }; bool is_undefined = is_pole(b1) || is_pole(b2); diff --git a/stan/math/prim/fun/hypergeometric_3F2.hpp b/stan/math/prim/fun/hypergeometric_3F2.hpp index a2b2e5f61ed..e57d60f4fca 100644 --- a/stan/math/prim/fun/hypergeometric_3F2.hpp +++ b/stan/math/prim/fun/hypergeometric_3F2.hpp @@ -68,9 +68,10 @@ inline return_type_t hypergeometric_3F2_infsum( b_array += 1.0; a_signs = sign(value_of_rec(a_array)); b_signs = sign(value_of_rec(b_array)); - t_sign = a_signs.prod() * b_signs.prod() * t_sign; + t_sign = z_sign * a_signs.prod() * b_signs.prod() * t_sign; } - if (k == max_steps) { + // The loop ends with k = max_steps + 1 when the steps run out + if (k > max_steps) { throw_domain_error("hypergeometric_3F2", "k (internal counter)", max_steps, "exceeded iterations, hypergeometric function did not ", "converge."); diff --git a/test/unit/math/prim/err/check_3F2_converges_test.cpp b/test/unit/math/prim/err/check_3F2_converges_test.cpp index 1f1e5900e49..9f67e0ebc9e 100644 --- a/test/unit/math/prim/err/check_3F2_converges_test.cpp +++ b/test/unit/math/prim/err/check_3F2_converges_test.cpp @@ -115,3 +115,30 @@ TEST(passesOnConvergentArgs, Check3F2Converges) { EXPECT_NO_THROW(check_3F2_converges(function, a1, a2, a3, b1, b2, z)); EXPECT_NO_THROW(check_3F2_converges(function, a1, a2, a3, b1, b2, z)); } + +TEST(Check3F2Converges, polynomial_ends_at_smallest_numerator) { + using stan::math::check_3F2_converges; + const char* function = "check_3F2_converges"; + // a1 = -1 (or a3 = -1) ends the series at k = 1, before (b1)_k is zero + // from k = 3 on + EXPECT_NO_THROW( + check_3F2_converges(function, -1.0, -5.0, 1.0, -2.0, 1.0, 1.0)); + EXPECT_NO_THROW( + check_3F2_converges(function, -5.0, 1.0, -1.0, -2.0, 1.0, 1.0)); + // a1 = -3 ends the series at k = 3, and (b1)_3 is zero + EXPECT_THROW(check_3F2_converges(function, -3.0, -5.0, 1.0, -2.0, 1.0, 1.0), + std::domain_error); +} + +TEST(Check3F2Converges, infinite_series_pole) { + using stan::math::check_3F2_converges; + const char* function = "check_3F2_converges"; + // no numerator parameter ends the series, so a non-positive integer + // denominator parameter is a pole + EXPECT_THROW(check_3F2_converges(function, 1.0, 1.0, 1.0, -1.0, 1.0, 0.5), + std::domain_error); + EXPECT_THROW(check_3F2_converges(function, 1.0, 1.0, 1.0, 1.0, -2.0, 0.5), + std::domain_error); + EXPECT_THROW(check_3F2_converges(function, 1.0, 1.0, 1.0, 0.0, 1.0, 0.5), + std::domain_error); +} diff --git a/test/unit/math/prim/fun/grad_F32_test.cpp b/test/unit/math/prim/fun/grad_F32_test.cpp index f91aa0cd725..bad3056ef77 100644 --- a/test/unit/math/prim/fun/grad_F32_test.cpp +++ b/test/unit/math/prim/fun/grad_F32_test.cpp @@ -148,6 +148,19 @@ TEST(MathPrimScalFun, grad_F32_polynomial_large_parameters) { } } +// a1 = -1 ends the series at k = 1, before (b1)_k is zero from k = 3 on; the +// larger non-positive integer a2 = -5 does not matter. The series is +// 1 + a1 a2 a3 z / (b1 b2). The derivative with respect to a1 is not +// checked: a1 is the integer that ends the series. +TEST(MathPrimScalFun, grad_F32_polynomial_smallest_numerator_ends) { + std::vector g = {0.0, 0.25, -1.25, -0.625, 1.25, -2.5}; + double g_calc[6]; + stan::math::grad_F32(g_calc, -1.0, -5.0, 1.0, -2.0, 1.0, 0.5); + for (int i : {1, 2, 3, 4, 5}) { + EXPECT_NEAR(g[i], g_calc[i], 1e-15); + } +} + // // m = { // {1.0, 1.0, 1.0, 1.0, 1.0, 0.6, 1e-10}, diff --git a/test/unit/math/prim/fun/hypergeometric_3F2_test.cpp b/test/unit/math/prim/fun/hypergeometric_3F2_test.cpp index 087027e6a01..aa493a1cf97 100644 --- a/test/unit/math/prim/fun/hypergeometric_3F2_test.cpp +++ b/test/unit/math/prim/fun/hypergeometric_3F2_test.cpp @@ -96,3 +96,31 @@ TEST(MathPrimScalFun, F32_polynomial_equal_parameter_sums) { stan::math::hypergeometric_3F2({1.0, 4.5, -6.0}, {5.0, -5.5}, 1.0), 1e-12); } +// a1 = -1 ends the series at k = 1, before (b1)_k is zero from k = 3 on; the +// larger non-positive integer a2 = -5 does not matter +TEST(MathPrimScalFun, F32_polynomial_smallest_numerator_ends) { + EXPECT_NEAR( + -0.25, + stan::math::hypergeometric_3F2({-1.0, -5.0, 1.0}, {-2.0, 1.0}, 0.5), + 1e-15); +} +// the sign of z enters every term of the sum +TEST(MathPrimScalFun, F32_infsum_negative_z) { + Eigen::VectorXd a(3); + a << 1.0, 2.0, -3.0; + Eigen::VectorXd b(2); + b << 4.0, 5.0; + EXPECT_NEAR(204.0 / 175.0, + stan::math::internal::hypergeometric_3F2_infsum(a, b, -0.5), + 1e-15); +} +// the sum throws when the steps run out before the end of the series +TEST(MathPrimScalFun, F32_infsum_max_steps) { + Eigen::VectorXd a(3); + a << 1.0, 4.5, -6.0; + Eigen::VectorXd b(2); + b << 5.0, -5.5; + EXPECT_THROW( + stan::math::internal::hypergeometric_3F2_infsum(a, b, 1.0, 1e-6, 3), + std::domain_error); +} From 289d6462d43c6fd26348e7b24a5520c492fce4ca Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Fri, 9 Oct 2026 12:08:24 +0300 Subject: [PATCH 3/5] Stop a terminating 3F2 series only when a bound of its tail is negligible --- stan/math/prim/fun.hpp | 1 + stan/math/prim/fun/grad_F32.hpp | 117 ++++++++++++++-- stan/math/prim/fun/hypergeometric_3F2.hpp | 65 ++++++++- .../fun/hypergeometric_3F2_tail_bound.hpp | 130 ++++++++++++++++++ test/unit/math/prim/fun/grad_F32_test.cpp | 26 ++++ .../hypergeometric_3F2_tail_bound_test.cpp | 91 ++++++++++++ .../math/prim/fun/hypergeometric_3F2_test.cpp | 17 +++ .../math/rev/prob/beta_binomial_cdfs_test.cpp | 16 +++ 8 files changed, 451 insertions(+), 12 deletions(-) create mode 100644 stan/math/prim/fun/hypergeometric_3F2_tail_bound.hpp create mode 100644 test/unit/math/prim/fun/hypergeometric_3F2_tail_bound_test.cpp diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp index 6e767e25c45..a13ca9b260d 100644 --- a/stan/math/prim/fun.hpp +++ b/stan/math/prim/fun.hpp @@ -128,6 +128,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/fun/grad_F32.hpp b/stan/math/prim/fun/grad_F32.hpp index e9d7536eb8e..5cdf4b603e9 100644 --- a/stan/math/prim/fun/grad_F32.hpp +++ b/stan/math/prim/fun/grad_F32.hpp @@ -6,8 +6,12 @@ #include #include #include +#include #include +#include #include +#include +#include #include namespace stan { @@ -21,6 +25,12 @@ namespace math { * to within precision or throwing when the * function takes max_steps steps. * + * A terminating series (a1, a2 or a3 a non-positive integer) is a + * polynomial. Its sum stops at the first zero term, or earlier when bounds + * of the sums of all remaining terms of the series and of every requested + * gradient are below 1e-17 times the absolute values of their partial sums; + * precision does not apply to it. + * * This power-series representation converges for all gradients * under the same conditions as the 3F2 function itself. * @@ -45,7 +55,9 @@ namespace math { * @param[in] b1 b1 see generalized hypergeometric function definition. * @param[in] b2 b2 see generalized hypergeometric function definition. * @param[in] z z see generalized hypergeometric function definition. - * @param[in] precision precision of the infinite sum + * @param[in] precision precision of the infinite sum: the sum stops after a + * term below precision in absolute value. Not used for a + * terminating series. * @param[in] max_steps number of steps to take */ template a_val{value_of_rec(a1), value_of_rec(a2), + value_of_rec(a3)}; + const std::array b_val{value_of_rec(b1), value_of_rec(b2)}; + const double abs_z = std::fabs(value_of_rec(z)); + // index of the last term that can be nonzero + double last = max_steps + 1.0; + for (double a_i : a_val) { + if (a_i <= 0.0 && a_i == std::floor(a_i)) { + last = std::fmin(last, -a_i); + } + } + auto is_negligible = [](double increment, const auto& sum) { + return increment <= polynomial_precision * std::fabs(value_of_rec(sum)); + }; + double f_sum = 1.0; std::array term{0}; for (int k = 0; k <= max_steps; ++k) { // A numerator parameter that has reached zero ends the series (a @@ -92,13 +126,22 @@ inline void grad_F32(T1* g, const T2& a1, const T3& a2, const T4& a3, log_t_new += log(fabs(p)) + log_z; log_t_new_sign = p >= 0.0 ? log_t_new_sign : -log_t_new_sign; + bool negligible = false; + double t_new = 0.0; + if (is_polynomial) { + t_new = std::exp(value_of_rec(log_t_new)); + f_sum += value_of_rec(log_t_new_sign) * t_new; + negligible = is_negligible(t_new, f_sum); + } if constexpr (grad_a1) { term[0] = log_g_old_sign[0] * log_t_old_sign * exp(log_g_old[0] - log_t_old) + inv(a1 + k); log_g_old[0] = log_t_new + log(fabs(term[0])); log_g_old_sign[0] = term[0] >= 0.0 ? log_t_new_sign : -log_t_new_sign; - g[0] += log_g_old_sign[0] * exp(log_g_old[0]); + const T1 increment = exp(log_g_old[0]); + g[0] += log_g_old_sign[0] * increment; + negligible = negligible && is_negligible(value_of_rec(increment), g[0]); } if constexpr (grad_a2) { @@ -107,7 +150,9 @@ inline void grad_F32(T1* g, const T2& a1, const T3& a2, const T4& a3, + inv(a2 + k); log_g_old[1] = log_t_new + log(fabs(term[1])); log_g_old_sign[1] = term[1] >= 0.0 ? log_t_new_sign : -log_t_new_sign; - g[1] += log_g_old_sign[1] * exp(log_g_old[1]); + const T1 increment = exp(log_g_old[1]); + g[1] += log_g_old_sign[1] * increment; + negligible = negligible && is_negligible(value_of_rec(increment), g[1]); } if constexpr (grad_a3) { @@ -116,7 +161,9 @@ inline void grad_F32(T1* g, const T2& a1, const T3& a2, const T4& a3, + inv(a3 + k); log_g_old[2] = log_t_new + log(fabs(term[2])); log_g_old_sign[2] = term[2] >= 0.0 ? log_t_new_sign : -log_t_new_sign; - g[2] += log_g_old_sign[2] * exp(log_g_old[2]); + const T1 increment = exp(log_g_old[2]); + g[2] += log_g_old_sign[2] * increment; + negligible = negligible && is_negligible(value_of_rec(increment), g[2]); } if constexpr (grad_b1) { @@ -125,7 +172,9 @@ inline void grad_F32(T1* g, const T2& a1, const T3& a2, const T4& a3, - inv(b1 + k); log_g_old[3] = log_t_new + log(fabs(term[3])); log_g_old_sign[3] = term[3] >= 0.0 ? log_t_new_sign : -log_t_new_sign; - g[3] += log_g_old_sign[3] * exp(log_g_old[3]); + const T1 increment = exp(log_g_old[3]); + g[3] += log_g_old_sign[3] * increment; + negligible = negligible && is_negligible(value_of_rec(increment), g[3]); } if constexpr (grad_b2) { @@ -134,7 +183,9 @@ inline void grad_F32(T1* g, const T2& a1, const T3& a2, const T4& a3, - inv(b2 + k); log_g_old[4] = log_t_new + log(fabs(term[4])); log_g_old_sign[4] = term[4] >= 0.0 ? log_t_new_sign : -log_t_new_sign; - g[4] += log_g_old_sign[4] * exp(log_g_old[4]); + const T1 increment = exp(log_g_old[4]); + g[4] += log_g_old_sign[4] * increment; + negligible = negligible && is_negligible(value_of_rec(increment), g[4]); } if constexpr (grad_z) { @@ -143,10 +194,60 @@ inline void grad_F32(T1* g, const T2& a1, const T3& a2, const T4& a3, + inv(z); log_g_old[5] = log_t_new + log(fabs(term[5])); log_g_old_sign[5] = term[5] >= 0.0 ? log_t_new_sign : -log_t_new_sign; - g[5] += log_g_old_sign[5] * exp(log_g_old[5]); + const T1 increment = exp(log_g_old[5]); + g[5] += log_g_old_sign[5] * increment; + negligible = negligible && is_negligible(value_of_rec(increment), g[5]); } - if (log_t_new <= log(precision)) { + if (is_polynomial) { + if (negligible) { + // The remaining terms t_{k + 2}, ..., t_last have the ratios r_j for + // j = k + 1, ..., last - 1. A gradient term is t_j s_j, where s_j + // (term[i] for j = k + 1) grows by at most c per step. + const double n_left = last - k - 1; + if (n_left <= 0) { + return; + } + const double lo = k + 1; + const double hi = last - 1; + const auto sums = internal::hypergeometric_3F2_tail_sums( + internal::hypergeometric_3F2_ratio_bound(a_val, b_val, abs_z, lo, + hi), + n_left); + auto tail_is_negligible = [&](const auto& s, double c, + const auto& sum) { + return is_negligible( + t_new + * (std::fabs(value_of_rec(s)) * sums.first + c * sums.second), + sum); + }; + auto growth = [&](double p) { + return 1.0 / internal::hypergeometric_3F2_min_abs(p, lo, hi); + }; + bool done = is_negligible(t_new * sums.first, f_sum); + if constexpr (grad_a1) { + done = done && tail_is_negligible(term[0], growth(a_val[0]), g[0]); + } + if constexpr (grad_a2) { + done = done && tail_is_negligible(term[1], growth(a_val[1]), g[1]); + } + if constexpr (grad_a3) { + done = done && tail_is_negligible(term[2], growth(a_val[2]), g[2]); + } + if constexpr (grad_b1) { + done = done && tail_is_negligible(term[3], growth(b_val[0]), g[3]); + } + if constexpr (grad_b2) { + done = done && tail_is_negligible(term[4], growth(b_val[1]), g[4]); + } + if constexpr (grad_z) { + done = done && tail_is_negligible(term[5], 1.0 / abs_z, g[5]); + } + if (done) { + return; + } + } + } else if (log_t_new <= log(precision)) { return; // implicit abs } diff --git a/stan/math/prim/fun/hypergeometric_3F2.hpp b/stan/math/prim/fun/hypergeometric_3F2.hpp index e57d60f4fca..eca870d93b7 100644 --- a/stan/math/prim/fun/hypergeometric_3F2.hpp +++ b/stan/math/prim/fun/hypergeometric_3F2.hpp @@ -8,19 +8,44 @@ #include #include #include +#include #include #include #include #include +#include +#include namespace stan { namespace math { namespace internal { +/** + * Sum the power series of the hypergeometric function 3F2 term by term. + * + * `hypergeometric_3F2` calls this function at z = 1 with sum(b) <= sum(a). + * There `check_3F2_converges` accepts only a terminating series (a + * polynomial), so the sum stops at the first zero term. It stops earlier + * when a bound of the sum of all remaining terms is below `precision` times + * the absolute value of the partial sum. + * + * @tparam Ta type of Eigen/Std vector 'a' arguments + * @tparam Tb type of Eigen/Std vector 'b' arguments + * @tparam Tz type of z argument + * @param[in] a numerator parameters + * @param[in] b denominator parameters + * @param[in] z argument + * @param[in] precision relative precision of the sum. The default 1e-17 is + * below half an ulp of the sum. + * @param[in] max_steps number of steps to take + * @return the sum of the series + * @throw std::domain_error if the sum overflows or needs more than + * max_steps steps + */ template * = nullptr, require_stan_scalar_t* = nullptr> inline return_type_t hypergeometric_3F2_infsum( - const Ta& a, const Tb& b, const Tz& z, double precision = 1e-6, + const Ta& a, const Tb& b, const Tz& z, double precision = 1e-17, int max_steps = 1e5) { using T_return = return_type_t; Eigen::Array, 3, 1> a_array = as_array_or_scalar(a); @@ -37,9 +62,25 @@ inline return_type_t hypergeometric_3F2_infsum( int z_sign = sign(value_of_rec(z)); int t_sign = z_sign * a_signs.prod() * b_signs.prod(); + // For the bound of the remaining terms: the parameters, and the index of + // the last term that can be nonzero (the first zero numerator ends the + // series, otherwise the steps end the sum) + const std::array a_val{value_of_rec(a_array[0]), + value_of_rec(a_array[1]), + value_of_rec(a_array[2])}; + const std::array b_val{value_of_rec(b_array[0]), + value_of_rec(b_array[1])}; + const double abs_z = std::fabs(value_of_rec(z)); + double last = max_steps + 1.0; + for (double a_i : a_val) { + if (a_i <= 0.0 && a_i == std::floor(a_i)) { + last = std::fmin(last, -a_i); + } + } + int k = 0; - const double log_precision = log(precision); - while (k <= max_steps && log_t >= log_precision) { + double abs_term = 1.0; + while (k <= max_steps) { // A numerator parameter that has reached zero makes this term and every // later term zero: the series is a polynomial and has ended. Without // this stop the sign below is 0 while the magnitude keeps growing, and @@ -47,6 +88,20 @@ inline return_type_t hypergeometric_3F2_infsum( if ((value_of_rec(a_array) == 0.0).any()) { return t_acc; } + // Stop when the last term t_k and a bound of the sum of the remaining + // terms t_{k + 1}, ..., t_last are negligible against the partial sum. + // The ratios of the remaining terms are r_j for j = k, ..., last - 1. + const double abs_sum = std::fabs(value_of_rec(t_acc)); + if (abs_term <= precision * abs_sum) { + const double ratio_bound + = hypergeometric_3F2_ratio_bound(a_val, b_val, abs_z, k, last - 1); + const double tail + = abs_term + * hypergeometric_3F2_tail_sums(ratio_bound, last - k).first; + if (tail <= precision * abs_sum) { + return t_acc; + } + } // Replace zero values with 1 prior to taking the log so that we accumulate // 0.0 rather than -inf const auto& abs_apk = math::fabs((a_array == 0).select(1.0, a_array)); @@ -57,7 +112,9 @@ inline return_type_t hypergeometric_3F2_infsum( } log_t += p + log_z; - t_acc += t_sign * exp(log_t); + const auto term = exp(log_t); + t_acc += t_sign * term; + abs_term = value_of_rec(term); if (is_inf(t_acc)) { throw_domain_error("hypergeometric_3F2", "sum (output)", t_acc, diff --git a/stan/math/prim/fun/hypergeometric_3F2_tail_bound.hpp b/stan/math/prim/fun/hypergeometric_3F2_tail_bound.hpp new file mode 100644 index 00000000000..acceb734224 --- /dev/null +++ b/stan/math/prim/fun/hypergeometric_3F2_tail_bound.hpp @@ -0,0 +1,130 @@ +#ifndef STAN_MATH_PRIM_FUN_HYPERGEOMETRIC_3F2_TAIL_BOUND_HPP +#define STAN_MATH_PRIM_FUN_HYPERGEOMETRIC_3F2_TAIL_BOUND_HPP + +#include +#include +#include +#include +#include + +namespace stan { +namespace math { +namespace internal { + +/** + * Return the smallest value of |p + j| over the integers j in [lo, hi]. + * + * @param p parameter + * @param lo smallest j + * @param hi largest j, at least lo + * @return smallest |p + j| + */ +inline double hypergeometric_3F2_min_abs(double p, double lo, double hi) { + if ((p + lo) * (p + hi) > 0.0) { + return std::fmin(std::fabs(p + lo), std::fabs(p + hi)); + } + // p + j has a zero in [lo, hi]: the smallest value is at one of the two + // integers next to -p + const double j = std::floor(-p); + double min_abs = std::numeric_limits::infinity(); + for (double i : {j, j + 1}) { + if (i >= lo && i <= hi) { + min_abs = std::fmin(min_abs, std::fabs(p + i)); + } + } + return min_abs; +} + +/** + * Return an upper bound of |x + j| / |y + j| over the integers j in + * [lo, hi]. + * + * If neither x + j nor y + j changes sign in [lo, hi], the function of j is + * a Moebius transformation without a pole there, so it is monotone and has + * its largest value at lo or hi. Otherwise the bound is the largest + * |x + j| (at lo or hi) divided by the smallest |y + j|. + * + * @param x numerator parameter + * @param y denominator parameter + * @param lo smallest j + * @param hi largest j, at least lo + * @return upper bound + */ +inline double hypergeometric_3F2_factor_bound(double x, double y, double lo, + double hi) { + if ((x + lo) * (x + hi) >= 0.0 && (y + lo) * (y + hi) > 0.0) { + return std::fmax(std::fabs((x + lo) / (y + lo)), + std::fabs((x + hi) / (y + hi))); + } + return std::fmax(std::fabs(x + lo), std::fabs(x + hi)) + / hypergeometric_3F2_min_abs(y, lo, hi); +} + +/** + * Return an upper bound of the absolute value of the term ratio + * r_j = z (a1 + j)(a2 + j)(a3 + j) / ((b1 + j)(b2 + j)(1 + j)) + * of the 3F2 series over the integers j in [lo, hi]. + * + * The ratio is the product of three factors |a_i + j| / |d_l + j|, with + * d = (b1, b2, 1). Every assignment of the numerator parameters to the + * denominator parameters gives a bound; the function returns the smallest + * of the six. + * + * @param a numerator parameters + * @param b denominator parameters + * @param abs_z absolute value of the argument z + * @param lo smallest j + * @param hi largest j, at least lo + * @return upper bound of |r_j| + */ +inline double hypergeometric_3F2_ratio_bound(const std::array& a, + const std::array& b, + double abs_z, double lo, + double hi) { + const std::array d{b[0], b[1], 1.0}; + std::array, 3> factor; + for (int i = 0; i < 3; ++i) { + for (int l = 0; l < 3; ++l) { + factor[i][l] = hypergeometric_3F2_factor_bound(a[i], d[l], lo, hi); + } + } + static constexpr int pairings[6][3] + = {{0, 1, 2}, {0, 2, 1}, {1, 0, 2}, {1, 2, 0}, {2, 0, 1}, {2, 1, 0}}; + double bound = std::numeric_limits::infinity(); + for (const auto& p : pairings) { + bound + = std::fmin(bound, factor[0][p[0]] * factor[1][p[1]] * factor[2][p[2]]); + } + return abs_z * bound; +} + +/** + * Return upper bounds of sum_{d = 1}^{n} rho^d and sum_{d = 1}^{n} d rho^d. + * + * With these sums, a term t of the series and a bound rho of the next n + * term ratios bound the sum of the next n terms by |t| times the first + * sum. A derivative term t s_j with |s_j| <= |s| + (j - k) c has the bound + * |t| (|s| times the first sum + c times the second sum). + * + * @param rho bound of the absolute term ratios + * @param n number of remaining terms + * @return the two bounds + */ +inline std::pair hypergeometric_3F2_tail_sums(double rho, + double n) { + if (n <= 0) { + return {0.0, 0.0}; + } + if (rho < 1.0) { + return { + std::fmin(n * rho, rho / (1.0 - rho)), + std::fmin(0.5 * n * (n + 1) * rho, rho / ((1.0 - rho) * (1.0 - rho)))}; + } + const double rho_n = std::pow(rho, n); + return {n * rho_n, n * n * rho_n}; +} + +} // namespace internal +} // namespace math +} // namespace stan +#endif diff --git a/test/unit/math/prim/fun/grad_F32_test.cpp b/test/unit/math/prim/fun/grad_F32_test.cpp index bad3056ef77..5a7029ea544 100644 --- a/test/unit/math/prim/fun/grad_F32_test.cpp +++ b/test/unit/math/prim/fun/grad_F32_test.cpp @@ -161,6 +161,32 @@ TEST(MathPrimScalFun, grad_F32_polynomial_smallest_numerator_ends) { } } +// terminate by zero numerator at k = 116; the terms fall below 1e-6 before +// the end of the series, so the sum must not stop at an absolute tolerance +TEST(MathPrimScalFun, grad_F32_polynomial_small_terms) { + std::vector g + = {4.4789990965505127, 4.2105041895611112, 0.0, + -2.8590213153488275, 0.10966400496285045, 12.978200050704774}; + double g_calc[6]; + stan::math::grad_F32(g_calc, 1.0, 1.1, -116.0, 2.0, -125.0, 1.0); + for (int i : {0, 1, 3, 4, 5}) { + EXPECT_NEAR(g[i], g_calc[i], 1e-13 * std::fabs(g[i])); + } +} + +// all terms are positive; the term after 1 is about 2e-21, and the sum is +// about 6e36. With only the derivatives in a2 and b2 requested, their first +// terms are small too, so the sum must not stop at a small term. +TEST(MathPrimScalFun, grad_F32_polynomial_small_term_before_large_terms) { + const double dF_da2 = 1.4427627588420411e+36; + const double dF_db2 = 4.7972504951249748e+36; + double g_calc[6]; + stan::math::grad_F32( + g_calc, 1e-23, 300.0, -100.0, 1.0, -150.5, 1.0); + EXPECT_NEAR(dF_da2, g_calc[1], 1e-12 * dF_da2); + EXPECT_NEAR(dF_db2, g_calc[4], 1e-12 * dF_db2); +} + // // m = { // {1.0, 1.0, 1.0, 1.0, 1.0, 0.6, 1e-10}, diff --git a/test/unit/math/prim/fun/hypergeometric_3F2_tail_bound_test.cpp b/test/unit/math/prim/fun/hypergeometric_3F2_tail_bound_test.cpp new file mode 100644 index 00000000000..8793a023255 --- /dev/null +++ b/test/unit/math/prim/fun/hypergeometric_3F2_tail_bound_test.cpp @@ -0,0 +1,91 @@ +#include +#include +#include +#include +#include + +namespace { + +// largest |r_j| over the integers j in [lo, hi], by evaluation of every +// ratio r_j = z (a1 + j)(a2 + j)(a3 + j) / ((b1 + j)(b2 + j)(1 + j)) +double max_abs_ratio(const std::array& a, + const std::array& b, double z, int lo, int hi) { + double max_r = 0.0; + for (int j = lo; j <= hi; ++j) { + const double r = z * (a[0] + j) * (a[1] + j) * (a[2] + j) + / ((b[0] + j) * (b[1] + j) * (1.0 + j)); + max_r = std::fmax(max_r, std::fabs(r)); + } + return max_r; +} + +} // namespace + +TEST(MathPrimScalFun, hypergeometric_3F2_min_abs) { + using stan::math::internal::hypergeometric_3F2_min_abs; + // no zero of p + j in [1, 10] + EXPECT_DOUBLE_EQ(3.5, hypergeometric_3F2_min_abs(2.5, 1, 10)); + EXPECT_DOUBLE_EQ(5.25, hypergeometric_3F2_min_abs(-15.25, 1, 10)); + // zero at j = 5.25, between the integers 5 and 6 + EXPECT_DOUBLE_EQ(0.25, hypergeometric_3F2_min_abs(-5.25, 1, 10)); + EXPECT_DOUBLE_EQ(0.0, hypergeometric_3F2_min_abs(-3.0, 1, 10)); +} + +TEST(MathPrimScalFun, hypergeometric_3F2_ratio_bound_is_upper_bound) { + using stan::math::internal::hypergeometric_3F2_ratio_bound; + struct series { + std::array a; + std::array b; + double z; + int lo; + int hi; + }; + const std::vector cases = { + // numerator and denominator parameters with zeros inside [lo, hi] + {{0.5, -2.5, -10.0}, {-3.5, 1.5}, 0.7, 0, 9}, + {{2.0, 3.0, -20.0}, {-7.25, 0.5}, -0.9, 0, 19}, + {{-4.5, 7.0, -30.0}, {-12.5, -40.5}, 1.0, 3, 29}, + // a tiny first term before large terms + {{1e-23, 300.0, -100.0}, {1.0, -150.5}, 1.0, 0, 99}, + {{1e-23, 300.0, -100.0}, {1.0, -150.5}, 1.0, 60, 99}, + // the series of beta_binomial_lcdf(3 | 10, 0.5, 0.5) + {{1.0, 4.5, -6.0}, {5.0, -5.5}, 1.0, 0, 5}, + }; + for (const auto& s : cases) { + const double bound + = hypergeometric_3F2_ratio_bound(s.a, s.b, std::fabs(s.z), s.lo, s.hi); + EXPECT_GE(bound * (1 + 1e-14), max_abs_ratio(s.a, s.b, s.z, s.lo, s.hi)); + } +} + +TEST(MathPrimScalFun, hypergeometric_3F2_ratio_bound_is_tight) { + using stan::math::internal::hypergeometric_3F2_ratio_bound; + // the series of beta_binomial_lcdf(5000 | 10000, 1000, 3000): every factor + // decreases in j, so the bound equals the ratio at j = lo + const std::array a{1.0, 6001.0, -4999.0}; + const std::array b{5002.0, -7998.0}; + for (int lo : {0, 100, 1000, 4000}) { + const double max_r = max_abs_ratio(a, b, 1.0, lo, 4998); + EXPECT_NEAR(max_r, hypergeometric_3F2_ratio_bound(a, b, 1.0, lo, 4998), + 1e-14 * max_r); + } +} + +TEST(MathPrimScalFun, hypergeometric_3F2_tail_sums) { + using stan::math::internal::hypergeometric_3F2_tail_sums; + for (double rho : {0.0, 0.3, 0.9, 1.0, 1.7}) { + for (int n : {1, 5, 40}) { + double sum1 = 0.0; + double sum2 = 0.0; + for (int d = 1; d <= n; ++d) { + sum1 += std::pow(rho, d); + sum2 += d * std::pow(rho, d); + } + const auto sums = hypergeometric_3F2_tail_sums(rho, n); + EXPECT_GE(sums.first, sum1 * (1 - 1e-15)); + EXPECT_GE(sums.second, sum2 * (1 - 1e-15)); + } + } + EXPECT_EQ(0.0, hypergeometric_3F2_tail_sums(0.5, 0).first); + EXPECT_EQ(0.0, hypergeometric_3F2_tail_sums(0.5, 0).second); +} diff --git a/test/unit/math/prim/fun/hypergeometric_3F2_test.cpp b/test/unit/math/prim/fun/hypergeometric_3F2_test.cpp index aa493a1cf97..d344d5fac8e 100644 --- a/test/unit/math/prim/fun/hypergeometric_3F2_test.cpp +++ b/test/unit/math/prim/fun/hypergeometric_3F2_test.cpp @@ -124,3 +124,20 @@ TEST(MathPrimScalFun, F32_infsum_max_steps) { stan::math::internal::hypergeometric_3F2_infsum(a, b, 1.0, 1e-6, 3), std::domain_error); } +// terminate by zero numerator at k = 116; the terms fall below 1e-6 before +// the end of the series, so the sum must not stop at an absolute tolerance +TEST(MathPrimScalFun, F32_polynomial_small_terms) { + EXPECT_NEAR( + 3.1748878200738662, + stan::math::hypergeometric_3F2({1.0, 1.1, -116.0}, {2.0, -125.0}, 1.0), + 1e-13); +} +// all terms are positive; the term after 1 is about 2e-21, and the sum is +// about 6e36, so the sum must not stop at a small term +TEST(MathPrimScalFun, F32_polynomial_small_term_before_large_terms) { + const double F = 5.7668140904936036e+36; + EXPECT_NEAR(F, + stan::math::hypergeometric_3F2({1e-23, 300.0, -100.0}, + {1.0, -150.5}, 1.0), + 1e-13 * F); +} diff --git a/test/unit/math/rev/prob/beta_binomial_cdfs_test.cpp b/test/unit/math/rev/prob/beta_binomial_cdfs_test.cpp index 1074e6584f8..dcbc1eee6d2 100644 --- a/test/unit/math/rev/prob/beta_binomial_cdfs_test.cpp +++ b/test/unit/math/rev/prob/beta_binomial_cdfs_test.cpp @@ -62,6 +62,22 @@ const std::vector beta_binomial_cdfs_points = { {-5.3033049080590758, -0.0099502487562189055, 1.9900497512437811}, {-0.0049875415110390736, 0.000049751243781094527, -0.0099502487562189055}}}, + // the terms of the 117-term 3F2 series fall below 1e-6 before the end + // of the series; the second point has the mirrored series + {0, + 117, + 0.1, + 10.0, + {{0.77231348114687892, -1.9916580324442808, 0.0074694843556493395}, + {-0.25836474771713782, -2.5788207522762463, 0.0096715705966406531}, + {-1.4797855134046491, 8.7473691568409674, -0.032806001836533192}}}, + {116, + 117, + 10.0, + 0.1, + {{0.22768651885312108, -0.0074694843556493395, 1.9916580324442808}, + {-1.4797855134046491, -0.032806001836533192, 8.7473691568409674}, + {-0.25836474771713782, 0.0096715705966406531, -2.5788207522762463}}}, }; template From a8525b6d94e144f1e4a9536913420240b9de2646 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Fri, 9 Oct 2026 12:08:39 +0300 Subject: [PATCH 4/5] Use the beta_binomial_lcdf mirror at alpha = 1 and request only the used 3F2 gradients --- stan/math/prim/prob/beta_binomial_lccdf.hpp | 3 ++- stan/math/prim/prob/beta_binomial_lcdf.hpp | 12 ++++----- .../math/rev/prob/beta_binomial_cdfs_test.cpp | 27 +++++++++++++++++++ 3 files changed, 35 insertions(+), 7 deletions(-) diff --git a/stan/math/prim/prob/beta_binomial_lccdf.hpp b/stan/math/prim/prob/beta_binomial_lccdf.hpp index 8f17cb47054..8df12ccc380 100644 --- a/stan/math/prim/prob/beta_binomial_lccdf.hpp +++ b/stan/math/prim/prob/beta_binomial_lccdf.hpp @@ -123,7 +123,8 @@ inline return_type_t beta_binomial_lccdf( T_partials_return dF[6]; if constexpr (is_any_autodiff_v) { - grad_F32(dF, one, mu, -N_dbl + n_dbl + 1, n_dbl + 2, 1 - nu, one); + grad_F32( + dF, one, mu, -N_dbl + n_dbl + 1, n_dbl + 2, 1 - nu, one); } if constexpr (is_autodiff_v) { partials<0>(ops_partials)[i] diff --git a/stan/math/prim/prob/beta_binomial_lcdf.hpp b/stan/math/prim/prob/beta_binomial_lcdf.hpp index d95431a2960..a0407e835f1 100644 --- a/stan/math/prim/prob/beta_binomial_lcdf.hpp +++ b/stan/math/prim/prob/beta_binomial_lcdf.hpp @@ -115,14 +115,12 @@ inline return_type_t beta_binomial_lcdf(const T_n& n, - lbeta(N_minus_n, n_dbl + 2); C = F * exp(C) / (N_dbl + 1); - if (C > 0.5 && alpha_dbl != 1.0) { + if (C > 0.5) { // The cdf is below 1/2, and log(1 - C) loses its digits; C can even // round to 1 or above. Use the mirror: X <= n if and only if // N - X > N - n - 1, and N - X is beta_binomial(N, beta, alpha). Its // ccdf is the same series with the shapes swapped and needs no - // complement. For alpha = 1 the mirrored series has the parameters - // -n and 1 - alpha - n = -n, whose terms are 0 / 0 at k = n + 1, so - // that case keeps log(1 - C). + // complement. const T_partials_return n_m = N_minus_n - 1; const T_partials_return mu_m = beta_dbl + N_minus_n; const T_partials_return nu_m = alpha_dbl + n_dbl; @@ -134,7 +132,8 @@ inline return_type_t beta_binomial_lcdf(const T_n& n, - lbeta(n_dbl + 1, N_minus_n + 1) + log(F_m) - log(N_dbl + 1); if constexpr (is_any_autodiff_v) { T_partials_return dF_m[6]; - grad_F32(dF_m, one, mu_m, -n_dbl, n_m + 2, 1 - nu_m, one); + grad_F32( + dF_m, one, mu_m, -n_dbl, n_m + 2, 1 - nu_m, one); const T_partials_return dpsi_total = digamma_diff(alpha_dbl + beta_dbl, N_dbl); if constexpr (is_autodiff_v) { @@ -161,7 +160,8 @@ inline return_type_t beta_binomial_lcdf(const T_n& n, T_partials_return dF[6]; if constexpr (is_any_autodiff_v) { - grad_F32(dF, one, mu, 1 - N_minus_n, n_dbl + 2, 1 - nu, one); + grad_F32( + dF, one, mu, 1 - N_minus_n, n_dbl + 2, 1 - nu, one); } if constexpr (is_autodiff_v) { const T_partials_return g diff --git a/test/unit/math/rev/prob/beta_binomial_cdfs_test.cpp b/test/unit/math/rev/prob/beta_binomial_cdfs_test.cpp index dcbc1eee6d2..d241715322f 100644 --- a/test/unit/math/rev/prob/beta_binomial_cdfs_test.cpp +++ b/test/unit/math/rev/prob/beta_binomial_cdfs_test.cpp @@ -1,6 +1,7 @@ #include #include #include +#include #include #include @@ -135,3 +136,29 @@ TEST_F(AgradRev, ProbDistributionsBetaBinomial_lccdf_terminating_series) { return stan::math::beta_binomial_lccdf(n, N, alpha, beta); }); } + +// alpha = 1 and a small cdf: beta_binomial_lcdf uses the mirrored series, +// which has a denominator parameter equal to the numerator parameter that +// ends the series +TEST_F(AgradRev, ProbDistributionsBetaBinomial_lcdf_mirror_alpha_one) { + using stan::math::var; + // n, N, alpha, beta, lcdf, d/dalpha, d/dbeta + const std::vector> points + = {{0, 117, 1.0, 0.1, -7.0656133635977173, -5.1910469886639266, + 9.9914602903501275}, + {0, 1000, 1.0, 0.01, -11.512935464920229, -7.4691506484692035, + 99.999000009999898}}; + for (const auto& p : points) { + var alpha = p[2]; + var beta = p[3]; + var lcdf = stan::math::beta_binomial_lcdf( + static_cast(p[0]), static_cast(p[1]), alpha, beta); + std::vector vars = {alpha, beta}; + std::vector grad; + lcdf.grad(vars, grad); + EXPECT_NEAR(p[4], lcdf.val(), 1e-12 * std::fabs(p[4])); + EXPECT_NEAR(p[5], grad[0], 1e-12 * std::fabs(p[5])); + EXPECT_NEAR(p[6], grad[1], 1e-12 * std::fabs(p[6])); + stan::math::recover_memory(); + } +} From 8d0cf249f7f708364ced0b7983f00515fb4b552d Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Fri, 9 Oct 2026 12:08:58 +0300 Subject: [PATCH 5/5] Compute beta_binomial_cdf from beta_binomial_lcdf --- stan/math/prim/prob/beta_binomial_cdf.hpp | 105 +----------------- .../math/rev/prob/beta_binomial_cdfs_test.cpp | 27 +++++ 2 files changed, 33 insertions(+), 99 deletions(-) diff --git a/stan/math/prim/prob/beta_binomial_cdf.hpp b/stan/math/prim/prob/beta_binomial_cdf.hpp index 0bab43798d7..7f1406c4ce7 100644 --- a/stan/math/prim/prob/beta_binomial_cdf.hpp +++ b/stan/math/prim/prob/beta_binomial_cdf.hpp @@ -3,21 +3,9 @@ #include #include -#include -#include -#include -#include #include -#include -#include -#include -#include -#include -#include #include -#include -#include -#include +#include namespace stan { namespace math { @@ -27,6 +15,10 @@ namespace math { * size, prior success, and prior failure parameters. Given containers of * matching sizes, returns the product of probabilities. * + * The CDF is the exponential of `beta_binomial_lcdf`, which computes a small + * probability from the mirrored series instead of the complement of the + * upper tail. + * * @tparam T_n type of success parameter * @tparam T_N type of population size parameter * @tparam T_size1 type of prior success parameter @@ -45,8 +37,6 @@ inline return_type_t beta_binomial_cdf(const T_n& n, const T_N& N, const T_size1& alpha, const T_size2& beta) { - using T_partials_return = partials_return_t; - using std::exp; using T_N_ref = ref_type_t; using T_alpha_ref = ref_type_t; using T_beta_ref = ref_type_t; @@ -68,90 +58,7 @@ inline return_type_t beta_binomial_cdf(const T_n& n, check_positive_finite(function, "Second prior sample size parameter", beta_ref); - T_partials_return P(1.0); - auto ops_partials = make_partials_propagator(alpha_ref, beta_ref); - - scalar_seq_view n_vec(n); - scalar_seq_view N_vec(N_ref); - scalar_seq_view alpha_vec(alpha_ref); - scalar_seq_view beta_vec(beta_ref); - size_t max_size_seq_view = max_size(n, N, alpha, beta); - - // Explicit return for extreme values - // The gradients are technically ill-defined, but treated as zero - for (size_t i = 0; i < stan::math::size(n); i++) { - if (n_vec.val(i) < 0) { - return ops_partials.build(0.0); - } - } - - for (size_t i = 0; i < max_size_seq_view; i++) { - // Explicit results for extreme values - // The gradients are technically ill-defined, but treated as zero - if (n_vec.val(i) >= N_vec.val(i)) { - continue; - } - - const T_partials_return n_dbl = n_vec.val(i); - const T_partials_return N_dbl = N_vec.val(i); - const T_partials_return alpha_dbl = alpha_vec.val(i); - const T_partials_return beta_dbl = beta_vec.val(i); - const T_partials_return N_minus_n = N_dbl - n_dbl; - const T_partials_return mu = alpha_dbl + n_dbl + 1; - const T_partials_return nu = beta_dbl + N_minus_n - 1; - const T_partials_return one = 1; - - const T_partials_return F = hypergeometric_3F2({one, mu, 1 - N_minus_n}, - {n_dbl + 2, 1 - nu}, one); - - // lbeta(nu, mu) - lbeta(alpha, beta) without the cancellation for large - // shapes - T_partials_return C - = internal::log_beta_ratio( - alpha_dbl, beta_dbl, n_dbl + 1, N_minus_n - 1, - internal::log_beta_ratio_denominator(alpha_dbl, beta_dbl)) - - lbeta(N_minus_n, n_dbl + 2); - C = F * exp(C) / (N_dbl + 1); - - const T_partials_return Pi = 1 - C; - - P *= Pi; - - // digamma(alpha + beta) - digamma(mu + nu), mu + nu = alpha + beta + N - T_partials_return digammaDiff - = is_constant_all::value - ? 0 - : -digamma_diff(alpha_dbl + beta_dbl, N_dbl); - - T_partials_return dF[6]; - if constexpr (is_any_autodiff_v) { - grad_F32(dF, one, mu, 1 - N_minus_n, n_dbl + 2, 1 - nu, one); - } - if constexpr (is_autodiff_v) { - const T_partials_return g - = -C * (digamma_diff(alpha_dbl, n_dbl + 1) + digammaDiff + dF[1] / F); - partials<0>(ops_partials)[i] += g / Pi; - } - if constexpr (is_autodiff_v) { - const T_partials_return g - = -C - * (digamma_diff(beta_dbl, N_minus_n - 1) + digammaDiff - dF[4] / F); - partials<1>(ops_partials)[i] += g / Pi; - } - } - - if constexpr (is_autodiff_v) { - for (size_t i = 0; i < stan::math::size(alpha); ++i) { - partials<0>(ops_partials)[i] *= P; - } - } - if constexpr (is_autodiff_v) { - for (size_t i = 0; i < stan::math::size(beta); ++i) { - partials<1>(ops_partials)[i] *= P; - } - } - - return ops_partials.build(P); + return exp(beta_binomial_lcdf(n, N_ref, alpha_ref, beta_ref)); } } // namespace math diff --git a/test/unit/math/rev/prob/beta_binomial_cdfs_test.cpp b/test/unit/math/rev/prob/beta_binomial_cdfs_test.cpp index d241715322f..0267707fded 100644 --- a/test/unit/math/rev/prob/beta_binomial_cdfs_test.cpp +++ b/test/unit/math/rev/prob/beta_binomial_cdfs_test.cpp @@ -162,3 +162,30 @@ TEST_F(AgradRev, ProbDistributionsBetaBinomial_lcdf_mirror_alpha_one) { stan::math::recover_memory(); } } + +// a cdf far below 1: beta_binomial_cdf must not form it as the complement +// of the upper tail, which rounds to 0 +TEST_F(AgradRev, ProbDistributionsBetaBinomial_cdf_small) { + using stan::math::var; + // n, N, alpha, beta, cdf, d/dalpha, d/dbeta + const std::vector> points + = {{0, 117, 300.0, 300.0, 8.1475432021262664e-32, -1.4525622365169557e-32, + 1.2342679116106887e-32}, + {1, 10, 1e4, 0.1, 5.2292970409674409e-32, -4.7040109793376876e-35, + 6.5758879077878958e-31}, + {4, 10, 1e4, 0.1, 3.133010009882612e-21, -1.8786628411897913e-24, + 3.8056360838721468e-20}}; + for (const auto& p : points) { + var alpha = p[2]; + var beta = p[3]; + var cdf = stan::math::beta_binomial_cdf( + static_cast(p[0]), static_cast(p[1]), alpha, beta); + std::vector vars = {alpha, beta}; + std::vector grad; + cdf.grad(vars, grad); + EXPECT_NEAR(p[4], cdf.val(), 1e-12 * std::fabs(p[4])); + EXPECT_NEAR(p[5], grad[0], 1e-12 * std::fabs(p[5])); + EXPECT_NEAR(p[6], grad[1], 1e-12 * std::fabs(p[6])); + stan::math::recover_memory(); + } +}