From 951847625768476ee9f922d0a9e817dbc6b55a80 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Thu, 22 Oct 2026 10:10:03 +0300 Subject: [PATCH 1/5] Fix discrete cdfs when one vector element has probability 0 or 1 --- stan/math/opencl/prim/bernoulli_lccdf.hpp | 18 +- stan/math/prim/prob/bernoulli_lccdf.hpp | 18 +- stan/math/prim/prob/beta_binomial_lccdf.hpp | 17 +- stan/math/prim/prob/beta_neg_binomial_cdf.hpp | 6 +- .../prim/prob/beta_neg_binomial_lccdf.hpp | 15 +- .../math/prim/prob/beta_neg_binomial_lcdf.hpp | 6 +- stan/math/prim/prob/binomial_lccdf.hpp | 22 +- stan/math/prim/prob/neg_binomial_2_cdf.hpp | 6 +- stan/math/prim/prob/neg_binomial_cdf.hpp | 6 +- stan/math/prim/prob/neg_binomial_lccdf.hpp | 15 +- stan/math/prim/prob/neg_binomial_lcdf.hpp | 6 +- stan/math/prim/prob/yule_simon_cdf.hpp | 12 +- stan/math/prim/prob/yule_simon_lccdf.hpp | 11 +- stan/math/prim/prob/yule_simon_lcdf.hpp | 12 +- .../prob/discrete_cdfs_edge_elements_test.cpp | 234 ++++++++++++++++++ 15 files changed, 315 insertions(+), 89 deletions(-) create mode 100644 test/unit/math/mix/prob/discrete_cdfs_edge_elements_test.cpp diff --git a/stan/math/opencl/prim/bernoulli_lccdf.hpp b/stan/math/opencl/prim/bernoulli_lccdf.hpp index 42a9f7cac6d..c1c930a7ce3 100644 --- a/stan/math/opencl/prim/bernoulli_lccdf.hpp +++ b/stan/math/opencl/prim/bernoulli_lccdf.hpp @@ -49,24 +49,22 @@ inline return_type_t bernoulli_lccdf(const T_n_cl& n, theta_val, "in the interval [0, 1]"); auto theta_bounded_expr = 0.0 <= theta_val && theta_val <= 1.0; - auto any_n_negative = colwise_max(cast(n < 0)); + // An element with probability 1 (n < 0) contributes nothing auto any_n_over_one = colwise_max(cast(n >= 1)); - auto P_expr = colwise_sum(log(theta_val)); - auto deriv = elt_divide(1.0, theta_val); + auto cond = n < 0; + auto P_expr = colwise_sum(select(cond, 0.0, log(theta_val))); + auto deriv = select(cond, 0.0, elt_divide(1.0, theta_val)); - matrix_cl any_n_negative_cl; matrix_cl any_n_over_one_cl; matrix_cl P_cl; matrix_cl deriv_cl; - results(check_theta_bounded, any_n_negative_cl, any_n_over_one_cl, P_cl, - deriv_cl) - = expressions(theta_bounded_expr, any_n_negative, any_n_over_one, P_expr, + results(check_theta_bounded, any_n_over_one_cl, P_cl, deriv_cl) + = expressions(theta_bounded_expr, any_n_over_one, P_expr, calc_if<(is_autodiff_v)>(deriv)); - if (from_matrix_cl(any_n_negative_cl).maxCoeff()) { - return 0.0; - } + // An element with probability 0 (n >= 1) makes the result -inf; the + // gradients are treated as zero if (from_matrix_cl(any_n_over_one_cl).maxCoeff()) { return NEGATIVE_INFTY; } diff --git a/stan/math/prim/prob/bernoulli_lccdf.hpp b/stan/math/prim/prob/bernoulli_lccdf.hpp index e8e8b4dce30..ca02ed259a8 100644 --- a/stan/math/prim/prob/bernoulli_lccdf.hpp +++ b/stan/math/prim/prob/bernoulli_lccdf.hpp @@ -9,6 +9,7 @@ #include #include #include +#include #include namespace stan { @@ -46,23 +47,18 @@ inline return_type_t bernoulli_lccdf(const T_n& n, auto ops_partials = make_partials_propagator(theta_ref); - // Explicit return for extreme values - // The gradients are technically ill-defined, but treated as zero - if (any(n_arr < 0)) { - return ops_partials.build(0.0); - } else if (any(n_arr >= 1)) { + // An element with probability 0 (n >= 1) makes the result -inf; the + // gradients are treated as zero + if (any(n_arr >= 1)) { return ops_partials.build(NEGATIVE_INFTY); } - size_t theta_size = math::size(theta_arr); - size_t n_size = math::size(n_arr); - double broadcast_n = theta_size == n_size ? 1 : n_size; - + // An element with probability 1 (n < 0) contributes nothing if constexpr (is_autodiff_v) { - partials<0>(ops_partials) = inv(theta_arr) * broadcast_n; + partials<0>(ops_partials) = select(n_arr == 0, inv(theta_arr), 0.0); } - return ops_partials.build(sum(log(theta_arr)) * broadcast_n); + return ops_partials.build(sum(select(n_arr == 0, log(theta_arr), 0.0))); } } // namespace math diff --git a/stan/math/prim/prob/beta_binomial_lccdf.hpp b/stan/math/prim/prob/beta_binomial_lccdf.hpp index d611d23b700..e8e7ba57420 100644 --- a/stan/math/prim/prob/beta_binomial_lccdf.hpp +++ b/stan/math/prim/prob/beta_binomial_lccdf.hpp @@ -76,19 +76,18 @@ inline return_type_t beta_binomial_lccdf( 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); + // An element with probability 0 (n >= N) makes the result -inf; the + // gradients are treated as zero + for (size_t i = 0; i < max_size_seq_view; i++) { + if (n_vec.val(i) >= N_vec.val(i)) { + return ops_partials.build(negative_infinity()); } } 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 neg infinity - if (n_vec.val(i) >= N_vec.val(i)) { - return ops_partials.build(negative_infinity()); + // An element with probability 1 (n < 0) contributes nothing + if (n_vec.val(i) < 0) { + continue; } const T_partials_return n_dbl = n_vec.val(i); diff --git a/stan/math/prim/prob/beta_neg_binomial_cdf.hpp b/stan/math/prim/prob/beta_neg_binomial_cdf.hpp index 00aa643ad8b..7fa3c3878bd 100644 --- a/stan/math/prim/prob/beta_neg_binomial_cdf.hpp +++ b/stan/math/prim/prob/beta_neg_binomial_cdf.hpp @@ -82,10 +82,10 @@ inline return_type_t beta_neg_binomial_cdf( T_partials_return cdf(1.0); auto ops_partials = make_partials_propagator(r_ref, alpha_ref, beta_ref); for (size_t i = 0; i < max_size_seq_view; i++) { - // Explicit return for extreme values - // The gradients are technically ill-defined, but treated as zero + // An element with probability 1 (n = INT_MAX, which stands for an + // infinite count) contributes nothing if (n_vec.val(i) == std::numeric_limits::max()) { - return 1.0; + continue; } auto n_dbl = n_vec.val(i); auto r_dbl = r_vec.val(i); diff --git a/stan/math/prim/prob/beta_neg_binomial_lccdf.hpp b/stan/math/prim/prob/beta_neg_binomial_lccdf.hpp index 49dc472e4e5..f94f6dca666 100644 --- a/stan/math/prim/prob/beta_neg_binomial_lccdf.hpp +++ b/stan/math/prim/prob/beta_neg_binomial_lccdf.hpp @@ -68,11 +68,11 @@ inline return_type_t beta_neg_binomial_lccdf( int size_n = stan::math::size(n); size_t max_size_seq_view = max_size(n, r, alpha, beta); - // Explicit return for extreme values - // The gradients are technically ill-defined, but treated as zero + // An element with probability 0 (n = INT_MAX, which stands for an infinite + // count) makes the result -inf; the gradients are treated as zero for (int i = 0; i < size_n; i++) { - if (n_vec.val(i) < 0) { - return 0.0; + if (n_vec.val(i) == std::numeric_limits::max()) { + return negative_infinity(); } } @@ -80,10 +80,9 @@ inline return_type_t beta_neg_binomial_lccdf( T_partials_return log_ccdf(0.0); auto ops_partials = make_partials_propagator(r_ref, alpha_ref, beta_ref); for (size_t i = 0; i < max_size_seq_view; i++) { - // Explicit return for extreme values - // The gradients are technically ill-defined, but treated as zero - if (n_vec.val(i) == std::numeric_limits::max()) { - return ops_partials.build(negative_infinity()); + // An element with probability 1 (n < 0) contributes nothing + if (n_vec.val(i) < 0) { + continue; } auto n_dbl = n_vec.val(i); auto r_dbl = r_vec.val(i); diff --git a/stan/math/prim/prob/beta_neg_binomial_lcdf.hpp b/stan/math/prim/prob/beta_neg_binomial_lcdf.hpp index 897ffa135a4..47e13fe6d83 100644 --- a/stan/math/prim/prob/beta_neg_binomial_lcdf.hpp +++ b/stan/math/prim/prob/beta_neg_binomial_lcdf.hpp @@ -80,10 +80,10 @@ inline return_type_t beta_neg_binomial_lcdf( T_partials_return log_cdf(0.0); auto ops_partials = make_partials_propagator(r_ref, alpha_ref, beta_ref); for (size_t i = 0; i < max_size_seq_view; i++) { - // Explicit return for extreme values - // The gradients are technically ill-defined, but treated as zero + // An element with probability 1 (n = INT_MAX, which stands for an + // infinite count) contributes nothing if (n_vec.val(i) == std::numeric_limits::max()) { - return 0.0; + continue; } auto n_dbl = n_vec.val(i); auto r_dbl = r_vec.val(i); diff --git a/stan/math/prim/prob/binomial_lccdf.hpp b/stan/math/prim/prob/binomial_lccdf.hpp index 5fe97006126..d1628f37754 100644 --- a/stan/math/prim/prob/binomial_lccdf.hpp +++ b/stan/math/prim/prob/binomial_lccdf.hpp @@ -69,25 +69,21 @@ inline return_type_t binomial_lccdf(const T_n& n, const T_N& N, scalar_seq_view theta_vec(theta_ref); size_t max_size_seq_view = max_size(n, N, theta); - // Explicit return for extreme values - // The gradients are technically ill-defined, - // but treated as negative infinity - for (size_t i = 0; i < stan::math::size(n); i++) { - if (n_vec.val(i) < 0) { - return ops_partials.build(0.0); + // An element with probability 0 (n >= N) makes the result -inf; the + // gradients are treated as zero + for (size_t i = 0; i < max_size_seq_view; i++) { + if (n_vec.val(i) >= N_vec.val(i)) { + return ops_partials.build(NEGATIVE_INFTY); } } for (size_t i = 0; i < max_size_seq_view; i++) { + // An element with probability 1 (n < 0) contributes nothing + if (n_vec.val(i) < 0) { + continue; + } const T_partials_return n_dbl = n_vec.val(i); const T_partials_return N_dbl = N_vec.val(i); - - // Explicit results for extreme values - // The gradients are technically ill-defined, but treated as zero - if (n_dbl >= N_dbl) { - return ops_partials.build(NEGATIVE_INFTY); - } - const T_partials_return theta_dbl = theta_vec.val(i); const T_partials_return Pi = 1.0 - inc_beta(N_dbl - n_dbl, n_dbl + 1, 1 - theta_dbl); diff --git a/stan/math/prim/prob/neg_binomial_2_cdf.hpp b/stan/math/prim/prob/neg_binomial_2_cdf.hpp index 47cab11d76d..c326a0c0cfb 100644 --- a/stan/math/prim/prob/neg_binomial_2_cdf.hpp +++ b/stan/math/prim/prob/neg_binomial_2_cdf.hpp @@ -78,10 +78,10 @@ inline return_type_t neg_binomial_2_cdf( } 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 + // An element with probability 1 (n = INT_MAX, which stands for an + // infinite count) contributes nothing if (n_vec.val(i) == std::numeric_limits::max()) { - return ops_partials.build(1.0); + continue; } const T_partials_return n_dbl_p1 = n_vec.val(i) + 1; diff --git a/stan/math/prim/prob/neg_binomial_cdf.hpp b/stan/math/prim/prob/neg_binomial_cdf.hpp index 31163102be1..1957810dd27 100644 --- a/stan/math/prim/prob/neg_binomial_cdf.hpp +++ b/stan/math/prim/prob/neg_binomial_cdf.hpp @@ -76,10 +76,10 @@ inline return_type_t neg_binomial_cdf( } 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 + // An element with probability 1 (n = INT_MAX, which stands for an + // infinite count) contributes nothing if (n_vec.val(i) == std::numeric_limits::max()) { - return ops_partials.build(1.0); + continue; } const T_partials_return n_dbl = n_vec.val(i); diff --git a/stan/math/prim/prob/neg_binomial_lccdf.hpp b/stan/math/prim/prob/neg_binomial_lccdf.hpp index cc24a6b79a7..ee5e328df10 100644 --- a/stan/math/prim/prob/neg_binomial_lccdf.hpp +++ b/stan/math/prim/prob/neg_binomial_lccdf.hpp @@ -57,11 +57,11 @@ inline return_type_t neg_binomial_lccdf( size_t size_n_alpha = max_size(n, alpha); size_t max_size_seq_view = max_size(n, alpha, beta_param); - // Explicit return for extreme values - // The gradients are technically ill-defined, but treated as zero + // An element with probability 0 (n = INT_MAX, which stands for an infinite + // count) makes the result -inf; the gradients are treated as zero for (size_t i = 0; i < size_n; i++) { - if (n_vec.val(i) < 0) { - return ops_partials.build(0.0); + if (n_vec.val(i) == std::numeric_limits::max()) { + return ops_partials.build(negative_infinity()); } } @@ -87,10 +87,9 @@ inline return_type_t neg_binomial_lccdf( } 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) == std::numeric_limits::max()) { - return ops_partials.build(negative_infinity()); + // An element with probability 1 (n < 0) contributes nothing + if (n_vec.val(i) < 0) { + continue; } const T_partials_return n_dbl = n_vec.val(i); diff --git a/stan/math/prim/prob/neg_binomial_lcdf.hpp b/stan/math/prim/prob/neg_binomial_lcdf.hpp index ba85ba41d0d..5b1f7b75202 100644 --- a/stan/math/prim/prob/neg_binomial_lcdf.hpp +++ b/stan/math/prim/prob/neg_binomial_lcdf.hpp @@ -87,10 +87,10 @@ inline return_type_t neg_binomial_lcdf( } 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 + // An element with probability 1 (n = INT_MAX, which stands for an + // infinite count) contributes nothing if (n_vec.val(i) == std::numeric_limits::max()) { - return ops_partials.build(0.0); + continue; } const T_partials_return n_dbl = n_vec.val(i); diff --git a/stan/math/prim/prob/yule_simon_cdf.hpp b/stan/math/prim/prob/yule_simon_cdf.hpp index 10a5ba2637e..4b19ff82fcf 100644 --- a/stan/math/prim/prob/yule_simon_cdf.hpp +++ b/stan/math/prim/prob/yule_simon_cdf.hpp @@ -52,20 +52,22 @@ inline return_type_t yule_simon_cdf(const T_n& n, scalar_seq_view alpha_vec(alpha_ref); const size_t max_size_seq_view = max_size(n_ref, alpha_ref); - // Explicit return for invalid or extreme values - // The gradients are technically ill-defined, but treated as zero + // An element with probability 0 (n < 1) makes the result 0; the gradients + // are treated as zero for (int i = 0; i < stan::math::size(n); i++) { if (n_vec.val(i) < 1.0) { return 0.0; } - if (n_vec.val(i) == std::numeric_limits::max()) { - return 1.0; - } } T_partials_return cdf(1.0); auto ops_partials = make_partials_propagator(alpha_ref); for (size_t i = 0; i < max_size_seq_view; i++) { + // An element with probability 1 (n = INT_MAX, which stands for an + // infinite count) contributes nothing + if (n_vec.val(i) == std::numeric_limits::max()) { + continue; + } auto np1 = n_vec.val(i) + 1.0; auto ap1 = alpha_vec.val(i) + 1.0; auto nap1 = n_vec.val(i) + ap1; diff --git a/stan/math/prim/prob/yule_simon_lccdf.hpp b/stan/math/prim/prob/yule_simon_lccdf.hpp index 335fa158e8c..9a289a77077 100644 --- a/stan/math/prim/prob/yule_simon_lccdf.hpp +++ b/stan/math/prim/prob/yule_simon_lccdf.hpp @@ -52,12 +52,9 @@ inline return_type_t yule_simon_lccdf(const T_n& n, scalar_seq_view alpha_vec(alpha_ref); const size_t max_size_seq_view = max_size(n_ref, alpha_ref); - // Explicit return for invalid or extreme values - // The gradients are technically ill-defined, but treated as zero + // An element with probability 0 (n = INT_MAX, which stands for an infinite + // count) makes the result -inf; the gradients are treated as zero for (int i = 0; i < stan::math::size(n); i++) { - if (n_vec.val(i) < 1.0) { - return 0.0; - } if (n_vec.val(i) == std::numeric_limits::max()) { return negative_infinity(); } @@ -66,6 +63,10 @@ inline return_type_t yule_simon_lccdf(const T_n& n, T_partials_return log_ccdf(0.0); auto ops_partials = make_partials_propagator(alpha_ref); for (size_t i = 0; i < max_size_seq_view; i++) { + // An element with probability 1 (n < 1) contributes nothing + if (n_vec.val(i) < 1.0) { + continue; + } auto np1 = n_vec.val(i) + 1.0; auto ap1 = alpha_vec.val(i) + 1.0; auto nap1 = n_vec.val(i) + ap1; diff --git a/stan/math/prim/prob/yule_simon_lcdf.hpp b/stan/math/prim/prob/yule_simon_lcdf.hpp index b84e54a6a82..2d3894659cd 100644 --- a/stan/math/prim/prob/yule_simon_lcdf.hpp +++ b/stan/math/prim/prob/yule_simon_lcdf.hpp @@ -52,20 +52,22 @@ inline return_type_t yule_simon_lcdf(const T_n& n, scalar_seq_view alpha_vec(alpha_ref); const size_t max_size_seq_view = max_size(n_ref, alpha_ref); - // Explicit return for invalid or extreme values - // The gradients are technically ill-defined, but treated as zero + // An element with probability 0 (n < 1) makes the result -inf; the + // gradients are treated as zero for (int i = 0; i < stan::math::size(n); i++) { if (n_vec.val(i) < 1.0) { return negative_infinity(); } - if (n_vec.val(i) == std::numeric_limits::max()) { - return 0.0; - } } T_partials_return log_cdf(0.0); auto ops_partials = make_partials_propagator(alpha_ref); for (size_t i = 0; i < max_size_seq_view; i++) { + // An element with probability 1 (n = INT_MAX, which stands for an + // infinite count) contributes nothing + if (n_vec.val(i) == std::numeric_limits::max()) { + continue; + } auto np1 = n_vec.val(i) + 1.0; auto ap1 = alpha_vec.val(i) + 1.0; auto nap1 = n_vec.val(i) + ap1; diff --git a/test/unit/math/mix/prob/discrete_cdfs_edge_elements_test.cpp b/test/unit/math/mix/prob/discrete_cdfs_edge_elements_test.cpp new file mode 100644 index 00000000000..eeb5b210c9c --- /dev/null +++ b/test/unit/math/mix/prob/discrete_cdfs_edge_elements_test.cpp @@ -0,0 +1,234 @@ +#include +#include +#include +#include +#include + +// In a vector argument of a discrete cdf, lcdf or lccdf, an element with +// probability 1 contributes nothing, and an element with probability 0 makes +// the result 0 (cdf) or -inf (lcdf, lccdf) with zero derivatives. The tests +// compare f({ordinary, edge}) and f({edge, ordinary}) with f({ordinary}) in +// reverse and forward mode. + +namespace { + +using stan::math::fvar; +using stan::math::var; + +struct value_and_grad { + double value; + std::vector grad; +}; + +template +value_and_grad eval_rev(const F& f, const std::vector& n, + const std::vector& p) { + std::vector pv(p.begin(), p.end()); + var y = f(n, pv); + value_and_grad res{y.val(), {}}; + y.grad(); + for (const auto& x : pv) { + res.grad.push_back(x.adj()); + } + stan::math::recover_memory(); + return res; +} + +template +value_and_grad eval_fwd(const F& f, const std::vector& n, + const std::vector& p) { + value_and_grad res{0, {}}; + for (size_t k = 0; k < p.size(); ++k) { + std::vector> pf; + for (size_t j = 0; j < p.size(); ++j) { + pf.emplace_back(p[j], j == k ? 1.0 : 0.0); + } + fvar y = f(n, pf); + res.value = y.val(); + res.grad.push_back(y.d_); + } + return res; +} + +void expect_same(double a, double b, const char* what) { + if (std::isinf(b)) { + EXPECT_EQ(a, b) << what; + } else { + EXPECT_NEAR(a, b, 1e-12 * std::fmax(1.0, std::fabs(b))) << what; + } +} + +template +void expect_probability_one(const F& f, int ordinary, int edge, + const std::vector& p) { + for (int mode = 0; mode < 2; ++mode) { + auto eval = [&](const std::vector& n) { + return mode == 0 ? eval_rev(f, n, p) : eval_fwd(f, n, p); + }; + const value_and_grad alone = eval({ordinary}); + for (const auto& n : + {std::vector{ordinary, edge}, std::vector{edge, ordinary}}) { + const value_and_grad both = eval(n); + expect_same(both.value, alone.value, + mode == 0 ? "rev value" : "fwd value"); + for (size_t k = 0; k < p.size(); ++k) { + expect_same(both.grad[k], alone.grad[k], + mode == 0 ? "rev gradient" : "fwd gradient"); + } + } + } +} + +template +void expect_probability_zero(const F& f, int ordinary, int edge, + const std::vector& p, double result) { + for (int mode = 0; mode < 2; ++mode) { + for (const auto& n : + {std::vector{ordinary, edge}, std::vector{edge, ordinary}}) { + const value_and_grad both + = mode == 0 ? eval_rev(f, n, p) : eval_fwd(f, n, p); + EXPECT_EQ(both.value, result) << (mode == 0 ? "rev" : "fwd"); + for (size_t k = 0; k < p.size(); ++k) { + EXPECT_EQ(both.grad[k], 0.0) << (mode == 0 ? "rev" : "fwd"); + } + } + } +} + +constexpr double NEG_INF = -std::numeric_limits::infinity(); +constexpr int INT_MAX_N = std::numeric_limits::max(); + +} // namespace + +TEST(ProbDiscreteCdfsEdgeElements, bernoulli) { + auto cdf = [](const auto& n, const auto& p) { + return stan::math::bernoulli_cdf(n, p[0]); + }; + auto lcdf = [](const auto& n, const auto& p) { + return stan::math::bernoulli_lcdf(n, p[0]); + }; + auto lccdf = [](const auto& n, const auto& p) { + return stan::math::bernoulli_lccdf(n, p[0]); + }; + const std::vector p{0.3}; + expect_probability_one(cdf, 0, 1, p); + expect_probability_one(lcdf, 0, 1, p); + expect_probability_one(lccdf, 0, -1, p); + expect_probability_zero(cdf, 0, -1, p, 0.0); + expect_probability_zero(lcdf, 0, -1, p, NEG_INF); + expect_probability_zero(lccdf, 0, 1, p, NEG_INF); +} + +TEST(ProbDiscreteCdfsEdgeElements, binomial) { + auto cdf = [](const auto& n, const auto& p) { + return stan::math::binomial_cdf(n, 10, p[0]); + }; + auto lcdf = [](const auto& n, const auto& p) { + return stan::math::binomial_lcdf(n, 10, p[0]); + }; + auto lccdf = [](const auto& n, const auto& p) { + return stan::math::binomial_lccdf(n, 10, p[0]); + }; + const std::vector p{0.3}; + expect_probability_one(cdf, 3, 10, p); + expect_probability_one(lcdf, 3, 10, p); + expect_probability_one(lccdf, 3, -1, p); + expect_probability_zero(cdf, 3, -1, p, 0.0); + expect_probability_zero(lcdf, 3, -1, p, NEG_INF); + expect_probability_zero(lccdf, 3, 10, p, NEG_INF); +} + +TEST(ProbDiscreteCdfsEdgeElements, beta_binomial) { + auto cdf = [](const auto& n, const auto& p) { + return stan::math::beta_binomial_cdf(n, 10, p[0], p[1]); + }; + auto lcdf = [](const auto& n, const auto& p) { + return stan::math::beta_binomial_lcdf(n, 10, p[0], p[1]); + }; + auto lccdf = [](const auto& n, const auto& p) { + return stan::math::beta_binomial_lccdf(n, 10, p[0], p[1]); + }; + const std::vector p{2.5, 1.5}; + expect_probability_one(cdf, 3, 10, p); + expect_probability_one(lcdf, 3, 10, p); + expect_probability_one(lccdf, 3, -1, p); + expect_probability_zero(cdf, 3, -1, p, 0.0); + expect_probability_zero(lcdf, 3, -1, p, NEG_INF); + expect_probability_zero(lccdf, 3, 10, p, NEG_INF); +} + +TEST(ProbDiscreteCdfsEdgeElements, neg_binomial) { + auto cdf = [](const auto& n, const auto& p) { + return stan::math::neg_binomial_cdf(n, p[0], p[1]); + }; + auto lcdf = [](const auto& n, const auto& p) { + return stan::math::neg_binomial_lcdf(n, p[0], p[1]); + }; + auto lccdf = [](const auto& n, const auto& p) { + return stan::math::neg_binomial_lccdf(n, p[0], p[1]); + }; + const std::vector p{2.5, 1.5}; + expect_probability_one(cdf, 3, INT_MAX_N, p); + expect_probability_one(lcdf, 3, INT_MAX_N, p); + expect_probability_one(lccdf, 3, -1, p); + expect_probability_zero(cdf, 3, -1, p, 0.0); + expect_probability_zero(lcdf, 3, -1, p, NEG_INF); + expect_probability_zero(lccdf, 3, INT_MAX_N, p, NEG_INF); +} + +TEST(ProbDiscreteCdfsEdgeElements, neg_binomial_2) { + auto cdf = [](const auto& n, const auto& p) { + return stan::math::neg_binomial_2_cdf(n, p[0], p[1]); + }; + auto lcdf = [](const auto& n, const auto& p) { + return stan::math::neg_binomial_2_lcdf(n, p[0], p[1]); + }; + auto lccdf = [](const auto& n, const auto& p) { + return stan::math::neg_binomial_2_lccdf(n, p[0], p[1]); + }; + const std::vector p{4.0, 2.0}; + expect_probability_one(cdf, 3, INT_MAX_N, p); + expect_probability_one(lcdf, 3, INT_MAX_N, p); + expect_probability_one(lccdf, 3, -1, p); + expect_probability_zero(cdf, 3, -1, p, 0.0); + expect_probability_zero(lcdf, 3, -1, p, NEG_INF); + expect_probability_zero(lccdf, 3, INT_MAX_N, p, NEG_INF); +} + +TEST(ProbDiscreteCdfsEdgeElements, beta_neg_binomial) { + auto cdf = [](const auto& n, const auto& p) { + return stan::math::beta_neg_binomial_cdf(n, p[0], p[1], p[2]); + }; + auto lcdf = [](const auto& n, const auto& p) { + return stan::math::beta_neg_binomial_lcdf(n, p[0], p[1], p[2]); + }; + auto lccdf = [](const auto& n, const auto& p) { + return stan::math::beta_neg_binomial_lccdf(n, p[0], p[1], p[2]); + }; + const std::vector p{3.0, 2.5, 1.5}; + expect_probability_one(cdf, 3, INT_MAX_N, p); + expect_probability_one(lcdf, 3, INT_MAX_N, p); + expect_probability_one(lccdf, 3, -1, p); + expect_probability_zero(cdf, 3, -1, p, 0.0); + expect_probability_zero(lcdf, 3, -1, p, NEG_INF); + expect_probability_zero(lccdf, 3, INT_MAX_N, p, NEG_INF); +} + +TEST(ProbDiscreteCdfsEdgeElements, yule_simon) { + auto cdf = [](const auto& n, const auto& p) { + return stan::math::yule_simon_cdf(n, p[0]); + }; + auto lcdf = [](const auto& n, const auto& p) { + return stan::math::yule_simon_lcdf(n, p[0]); + }; + auto lccdf = [](const auto& n, const auto& p) { + return stan::math::yule_simon_lccdf(n, p[0]); + }; + const std::vector p{1.5}; + expect_probability_one(cdf, 3, INT_MAX_N, p); + expect_probability_one(lcdf, 3, INT_MAX_N, p); + expect_probability_one(lccdf, 3, 0, p); + expect_probability_zero(cdf, 3, 0, p, 0.0); + expect_probability_zero(lcdf, 3, 0, p, NEG_INF); + expect_probability_zero(lccdf, 3, INT_MAX_N, p, NEG_INF); +} From 4072ed02fb23778d9282405acaf7cbda537b3841 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sat, 31 Oct 2026 17:50:32 +0200 Subject: [PATCH 2/5] Skip infinite counts in neg_binomial_2_lcdf before the beta lcdf call --- stan/math/prim/prob/neg_binomial_2_lcdf.hpp | 22 +++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/stan/math/prim/prob/neg_binomial_2_lcdf.hpp b/stan/math/prim/prob/neg_binomial_2_lcdf.hpp index 207e5532730..8e183f40c59 100644 --- a/stan/math/prim/prob/neg_binomial_2_lcdf.hpp +++ b/stan/math/prim/prob/neg_binomial_2_lcdf.hpp @@ -10,6 +10,8 @@ #include #include #include +#include +#include namespace stan { namespace math { @@ -43,10 +45,30 @@ inline return_type_t neg_binomial_2_lcdf( size_t size_n = stan::math::size(n); size_t size_phi_mu = max_size(mu, phi); + bool any_n_max = false; for (size_t i = 0; i < size_n; i++) { if (n_vec[i] < 0) { return LOG_ZERO; } + any_n_max = any_n_max || n_vec[i] == std::numeric_limits::max(); + } + + if (any_n_max) { + // An element with probability 1 (n = INT_MAX, which stands for an + // infinite count) contributes nothing + const size_t max_size_seq_view = max_size(n, mu, phi); + std::vector> phi_mu_k; + std::vector> phi_k; + std::vector np1_k; + for (size_t i = 0; i < max_size_seq_view; i++) { + if (n_vec[i] == std::numeric_limits::max()) { + continue; + } + phi_mu_k.push_back(phi_vec[i] / (phi_vec[i] + mu_vec[i])); + phi_k.push_back(phi_vec[i]); + np1_k.push_back(n_vec[i] + 1.0); + } + return beta_cdf_log(phi_mu_k, phi_k, np1_k); } VectorBuilder, T_location, From 3d3475fedba1bfd00d1a2466f021f87c8c1268ed Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sat, 31 Oct 2026 18:14:10 +0200 Subject: [PATCH 3/5] Apply the edge-element rule to poisson_binomial and test vector parameters --- stan/math/opencl/prim/bernoulli_lccdf.hpp | 2 +- stan/math/prim/prob/bernoulli_lccdf.hpp | 3 +- stan/math/prim/prob/neg_binomial_2_lcdf.hpp | 4 +- stan/math/prim/prob/neg_binomial_lcdf.hpp | 2 +- stan/math/prim/prob/poisson_binomial_cdf.hpp | 5 + .../math/prim/prob/poisson_binomial_lccdf.hpp | 15 +- stan/math/prim/prob/poisson_binomial_lcdf.hpp | 5 + .../prob/discrete_cdfs_edge_elements_test.cpp | 205 ++++++++++++++++++ .../prob/poisson_binomial_ccdf_log_test.cpp | 2 +- 9 files changed, 233 insertions(+), 10 deletions(-) diff --git a/stan/math/opencl/prim/bernoulli_lccdf.hpp b/stan/math/opencl/prim/bernoulli_lccdf.hpp index c1c930a7ce3..0770c4fc4c4 100644 --- a/stan/math/opencl/prim/bernoulli_lccdf.hpp +++ b/stan/math/opencl/prim/bernoulli_lccdf.hpp @@ -49,8 +49,8 @@ inline return_type_t bernoulli_lccdf(const T_n_cl& n, theta_val, "in the interval [0, 1]"); auto theta_bounded_expr = 0.0 <= theta_val && theta_val <= 1.0; - // An element with probability 1 (n < 0) contributes nothing auto any_n_over_one = colwise_max(cast(n >= 1)); + // An element with probability 1 (n < 0) contributes nothing auto cond = n < 0; auto P_expr = colwise_sum(select(cond, 0.0, log(theta_val))); auto deriv = select(cond, 0.0, elt_divide(1.0, theta_val)); diff --git a/stan/math/prim/prob/bernoulli_lccdf.hpp b/stan/math/prim/prob/bernoulli_lccdf.hpp index ca02ed259a8..5688fb2085e 100644 --- a/stan/math/prim/prob/bernoulli_lccdf.hpp +++ b/stan/math/prim/prob/bernoulli_lccdf.hpp @@ -53,7 +53,8 @@ inline return_type_t bernoulli_lccdf(const T_n& n, return ops_partials.build(NEGATIVE_INFTY); } - // An element with probability 1 (n < 0) contributes nothing + // Only the elements with n = 0 remain to contribute; an element with + // probability 1 (n < 0) contributes nothing if constexpr (is_autodiff_v) { partials<0>(ops_partials) = select(n_arr == 0, inv(theta_arr), 0.0); } diff --git a/stan/math/prim/prob/neg_binomial_2_lcdf.hpp b/stan/math/prim/prob/neg_binomial_2_lcdf.hpp index 8e183f40c59..f821516af82 100644 --- a/stan/math/prim/prob/neg_binomial_2_lcdf.hpp +++ b/stan/math/prim/prob/neg_binomial_2_lcdf.hpp @@ -55,7 +55,9 @@ inline return_type_t neg_binomial_2_lcdf( if (any_n_max) { // An element with probability 1 (n = INT_MAX, which stands for an - // infinite count) contributes nothing + // infinite count) contributes nothing. Remove such elements before the + // beta lcdf call: the beta lcdf gradient is not accurate when the + // second shape is 2^31. const size_t max_size_seq_view = max_size(n, mu, phi); std::vector> phi_mu_k; std::vector> phi_k; diff --git a/stan/math/prim/prob/neg_binomial_lcdf.hpp b/stan/math/prim/prob/neg_binomial_lcdf.hpp index 5b1f7b75202..fff1f7891df 100644 --- a/stan/math/prim/prob/neg_binomial_lcdf.hpp +++ b/stan/math/prim/prob/neg_binomial_lcdf.hpp @@ -74,7 +74,7 @@ inline return_type_t neg_binomial_lcdf( if constexpr (is_autodiff_v) { for (size_t i = 0; i < size_n; i++) { - digammaN_vec[i] = digamma(n_vec.val(i) + 1); + digammaN_vec[i] = digamma(n_vec.val(i) + 1.0); } for (size_t i = 0; i < size_alpha; i++) { digammaAlpha_vec[i] = digamma(alpha_vec.val(i)); diff --git a/stan/math/prim/prob/poisson_binomial_cdf.hpp b/stan/math/prim/prob/poisson_binomial_cdf.hpp index d520792d75a..839ba40f3ad 100644 --- a/stan/math/prim/prob/poisson_binomial_cdf.hpp +++ b/stan/math/prim/prob/poisson_binomial_cdf.hpp @@ -49,6 +49,11 @@ inline return_type_t poisson_binomial_cdf(const T_y& y, return_type_t lcdf = 0.0; for (size_t i = 0; i < max_sz; ++i) { + // An element with probability 1 (y equal to the number of trials) + // contributes nothing + if (y_vec[i] == static_cast(stan::math::size(theta_vec[i]))) { + continue; + } auto x = log_sum_exp(poisson_binomial_log_probs(y_vec[i], theta_vec[i])); lcdf += x; } diff --git a/stan/math/prim/prob/poisson_binomial_lccdf.hpp b/stan/math/prim/prob/poisson_binomial_lccdf.hpp index d63f5d2869e..c9cc805d928 100644 --- a/stan/math/prim/prob/poisson_binomial_lccdf.hpp +++ b/stan/math/prim/prob/poisson_binomial_lccdf.hpp @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -50,14 +51,18 @@ inline return_type_t poisson_binomial_lccdf(const T_y& y, 1.0); } + // An element with probability 0 (y equal to the number of trials) makes the + // result -inf; the gradients are treated as zero + for (size_t i = 0; i < max_sz; ++i) { + if (y_vec[i] == static_cast(stan::math::size(theta_vec[i]))) { + return negative_infinity(); + } + } + return_type_t lccdf = 0.0; for (size_t i = 0; i < max_sz; ++i) { if (stan::math::size(theta_vec[i]) == 1) { - if (y_vec[i] == 0) { - lccdf += log(theta_vec[i][0]); - } else { - lccdf -= stan::math::INFTY; - } + lccdf += log(theta_vec[i][0]); } else { auto x = log1m_exp( log_sum_exp(poisson_binomial_log_probs(y_vec[i], theta_vec[i]))); diff --git a/stan/math/prim/prob/poisson_binomial_lcdf.hpp b/stan/math/prim/prob/poisson_binomial_lcdf.hpp index 6ae587cff95..6872ff396dc 100644 --- a/stan/math/prim/prob/poisson_binomial_lcdf.hpp +++ b/stan/math/prim/prob/poisson_binomial_lcdf.hpp @@ -50,6 +50,11 @@ inline return_type_t poisson_binomial_lcdf(const T_y& y, return_type_t lcdf = 0.0; for (size_t i = 0; i < max_sz; ++i) { + // An element with probability 1 (y equal to the number of trials) + // contributes nothing + if (y_vec[i] == static_cast(stan::math::size(theta_vec[i]))) { + continue; + } auto x = log_sum_exp(poisson_binomial_log_probs(y_vec[i], theta_vec[i])); lcdf += x; } diff --git a/test/unit/math/mix/prob/discrete_cdfs_edge_elements_test.cpp b/test/unit/math/mix/prob/discrete_cdfs_edge_elements_test.cpp index eeb5b210c9c..aea780522e4 100644 --- a/test/unit/math/mix/prob/discrete_cdfs_edge_elements_test.cpp +++ b/test/unit/math/mix/prob/discrete_cdfs_edge_elements_test.cpp @@ -1,7 +1,9 @@ +#include #include #include #include #include +#include #include // In a vector argument of a discrete cdf, lcdf or lccdf, an element with @@ -95,6 +97,64 @@ void expect_probability_zero(const F& f, int ordinary, int edge, } } +// Parameter j of k, with p holding the k parameters one after the other, each +// as a vector of p.size() / k elements +template +std::vector part(const std::vector& p, size_t k, size_t j) { + const size_t m = p.size() / k; + return std::vector(p.begin() + j * m, p.begin() + (j + 1) * m); +} + +// Vector parameters: q_ord and q_edge are the parameters of the ordinary and +// of the edge element. The derivatives with respect to the parameters of the +// edge element must be exactly 0. +template +void expect_probability_one_vec(const F& f, int ordinary, int edge, + const std::vector& q_ord, + const std::vector& q_edge) { + const size_t k = q_ord.size(); + for (int mode = 0; mode < 2; ++mode) { + auto eval = [&](const std::vector& n, const std::vector& p) { + return mode == 0 ? eval_rev(f, n, p) : eval_fwd(f, n, p); + }; + const value_and_grad alone = eval({ordinary}, q_ord); + for (size_t pos_ord = 0; pos_ord < 2; ++pos_ord) { + const size_t pos_edge = 1 - pos_ord; + std::vector n(2); + n[pos_ord] = ordinary; + n[pos_edge] = edge; + std::vector p(2 * k); + for (size_t j = 0; j < k; ++j) { + p[2 * j + pos_ord] = q_ord[j]; + p[2 * j + pos_edge] = q_edge[j]; + } + const value_and_grad both = eval(n, p); + const char* what = mode == 0 ? "rev" : "fwd"; + expect_same(both.value, alone.value, what); + for (size_t j = 0; j < k; ++j) { + expect_same(both.grad[2 * j + pos_ord], alone.grad[j], what); + EXPECT_EQ(both.grad[2 * j + pos_edge], 0.0) << what; + } + } + } +} + +// A scalar edge count with vector parameters: the result is the given value +// with zero derivatives +template +void expect_scalar_edge(const F& f, int edge, const std::vector& p, + double result) { + auto g = [&](const std::vector& n, const auto& q) { return f(n[0], q); }; + for (int mode = 0; mode < 2; ++mode) { + const value_and_grad res + = mode == 0 ? eval_rev(g, {edge}, p) : eval_fwd(g, {edge}, p); + EXPECT_EQ(res.value, result) << (mode == 0 ? "rev" : "fwd"); + for (size_t k = 0; k < p.size(); ++k) { + EXPECT_EQ(res.grad[k], 0.0) << (mode == 0 ? "rev" : "fwd"); + } + } +} + constexpr double NEG_INF = -std::numeric_limits::infinity(); constexpr int INT_MAX_N = std::numeric_limits::max(); @@ -232,3 +292,148 @@ TEST(ProbDiscreteCdfsEdgeElements, yule_simon) { expect_probability_zero(lcdf, 3, 0, p, NEG_INF); expect_probability_zero(lccdf, 3, INT_MAX_N, p, NEG_INF); } + +TEST(ProbDiscreteCdfsEdgeElements, poisson_binomial) { + auto theta = [](const auto& y, const auto& p) { + using T = std::decay_t; + Eigen::Matrix th(3); + th << p[0], p[1], p[2]; + return std::vector>(y.size(), th); + }; + auto cdf = [&](const auto& y, const auto& p) { + return stan::math::poisson_binomial_cdf(y, theta(y, p)); + }; + auto lcdf = [&](const auto& y, const auto& p) { + return stan::math::poisson_binomial_lcdf(y, theta(y, p)); + }; + auto lccdf = [&](const auto& y, const auto& p) { + return stan::math::poisson_binomial_lccdf(y, theta(y, p)); + }; + const std::vector p{0.3, 0.6, 0.8}; + expect_probability_one(cdf, 1, 3, p); + expect_probability_one(lcdf, 1, 3, p); + expect_probability_zero(lccdf, 1, 3, p, NEG_INF); +} + +TEST(ProbDiscreteCdfsEdgeElements, vector_parameters) { + using stan::math::bernoulli_cdf; + using stan::math::bernoulli_lccdf; + using stan::math::bernoulli_lcdf; + const std::vector t_ord{0.3}, t_edge{0.6}; + auto b_cdf = [](const auto& n, const auto& p) { return bernoulli_cdf(n, p); }; + auto b_lcdf + = [](const auto& n, const auto& p) { return bernoulli_lcdf(n, p); }; + auto b_lccdf + = [](const auto& n, const auto& p) { return bernoulli_lccdf(n, p); }; + expect_probability_one_vec(b_cdf, 0, 1, t_ord, t_edge); + expect_probability_one_vec(b_lcdf, 0, 1, t_ord, t_edge); + expect_probability_one_vec(b_lccdf, 0, -1, t_ord, t_edge); + expect_probability_zero(b_lccdf, 0, 1, {0.3, 0.6}, NEG_INF); + + auto bin_lccdf = [](const auto& n, const auto& p) { + return stan::math::binomial_lccdf(n, 10, p); + }; + expect_probability_one_vec(bin_lccdf, 3, -1, t_ord, t_edge); + expect_probability_zero(bin_lccdf, 3, 10, {0.3, 0.6}, NEG_INF); + + const std::vector ab_ord{2.5, 1.5}, ab_edge{3.5, 0.8}; + const std::vector ab_both{2.5, 3.5, 1.5, 0.8}; + auto bb_lccdf = [](const auto& n, const auto& p) { + return stan::math::beta_binomial_lccdf(n, 10, part(p, 2, 0), part(p, 2, 1)); + }; + expect_probability_one_vec(bb_lccdf, 3, -1, ab_ord, ab_edge); + expect_probability_zero(bb_lccdf, 3, 10, ab_both, NEG_INF); + + auto nb_cdf = [](const auto& n, const auto& p) { + return stan::math::neg_binomial_cdf(n, part(p, 2, 0), part(p, 2, 1)); + }; + auto nb_lcdf = [](const auto& n, const auto& p) { + return stan::math::neg_binomial_lcdf(n, part(p, 2, 0), part(p, 2, 1)); + }; + auto nb_lccdf = [](const auto& n, const auto& p) { + return stan::math::neg_binomial_lccdf(n, part(p, 2, 0), part(p, 2, 1)); + }; + expect_probability_one_vec(nb_cdf, 3, INT_MAX_N, ab_ord, ab_edge); + expect_probability_one_vec(nb_lcdf, 3, INT_MAX_N, ab_ord, ab_edge); + expect_probability_one_vec(nb_lccdf, 3, -1, ab_ord, ab_edge); + expect_probability_zero(nb_lccdf, 3, INT_MAX_N, ab_both, NEG_INF); + + const std::vector mp_ord{4.0, 2.0}, mp_edge{1.5, 3.0}; + auto nb2_cdf = [](const auto& n, const auto& p) { + return stan::math::neg_binomial_2_cdf(n, part(p, 2, 0), part(p, 2, 1)); + }; + auto nb2_lcdf = [](const auto& n, const auto& p) { + return stan::math::neg_binomial_2_lcdf(n, part(p, 2, 0), part(p, 2, 1)); + }; + auto nb2_lccdf = [](const auto& n, const auto& p) { + return stan::math::neg_binomial_2_lccdf(n, part(p, 2, 0), part(p, 2, 1)); + }; + expect_probability_one_vec(nb2_cdf, 3, INT_MAX_N, mp_ord, mp_edge); + expect_probability_one_vec(nb2_lcdf, 3, INT_MAX_N, mp_ord, mp_edge); + expect_probability_one_vec(nb2_lccdf, 3, -1, mp_ord, mp_edge); + + const std::vector rab_ord{3.0, 2.5, 1.5}, rab_edge{2.0, 3.5, 1.2}; + const std::vector rab_both{3.0, 2.0, 2.5, 3.5, 1.5, 1.2}; + auto bnb_cdf = [](const auto& n, const auto& p) { + return stan::math::beta_neg_binomial_cdf(n, part(p, 3, 0), part(p, 3, 1), + part(p, 3, 2)); + }; + auto bnb_lcdf = [](const auto& n, const auto& p) { + return stan::math::beta_neg_binomial_lcdf(n, part(p, 3, 0), part(p, 3, 1), + part(p, 3, 2)); + }; + auto bnb_lccdf = [](const auto& n, const auto& p) { + return stan::math::beta_neg_binomial_lccdf(n, part(p, 3, 0), part(p, 3, 1), + part(p, 3, 2)); + }; + expect_probability_one_vec(bnb_cdf, 3, INT_MAX_N, rab_ord, rab_edge); + expect_probability_one_vec(bnb_lcdf, 3, INT_MAX_N, rab_ord, rab_edge); + expect_probability_one_vec(bnb_lccdf, 3, -1, rab_ord, rab_edge); + expect_probability_zero(bnb_lccdf, 3, INT_MAX_N, rab_both, NEG_INF); + + const std::vector a_ord{1.5}, a_edge{0.7}; + auto ys_cdf = [](const auto& n, const auto& p) { + return stan::math::yule_simon_cdf(n, p); + }; + auto ys_lcdf = [](const auto& n, const auto& p) { + return stan::math::yule_simon_lcdf(n, p); + }; + auto ys_lccdf = [](const auto& n, const auto& p) { + return stan::math::yule_simon_lccdf(n, p); + }; + expect_probability_one_vec(ys_cdf, 3, INT_MAX_N, a_ord, a_edge); + expect_probability_one_vec(ys_lcdf, 3, INT_MAX_N, a_ord, a_edge); + expect_probability_one_vec(ys_lccdf, 3, 0, a_ord, a_edge); + + // A scalar edge count with vector parameters + expect_scalar_edge(b_lccdf, -1, {0.3, 0.6}, 0.0); + expect_scalar_edge(bin_lccdf, -1, {0.3, 0.6}, 0.0); + expect_scalar_edge(bb_lccdf, -1, ab_both, 0.0); + expect_scalar_edge(nb_cdf, INT_MAX_N, ab_both, 1.0); + expect_scalar_edge(nb_lcdf, INT_MAX_N, ab_both, 0.0); + expect_scalar_edge(nb_lccdf, -1, ab_both, 0.0); + expect_scalar_edge(nb2_cdf, INT_MAX_N, {4.0, 1.5, 2.0, 3.0}, 1.0); + expect_scalar_edge(nb2_lcdf, INT_MAX_N, {4.0, 1.5, 2.0, 3.0}, 0.0); + expect_scalar_edge(nb2_lccdf, -1, {4.0, 1.5, 2.0, 3.0}, 0.0); + expect_scalar_edge(bnb_cdf, INT_MAX_N, rab_both, 1.0); + expect_scalar_edge(bnb_lcdf, INT_MAX_N, rab_both, 0.0); + expect_scalar_edge(bnb_lccdf, -1, rab_both, 0.0); + expect_scalar_edge(ys_cdf, INT_MAX_N, {1.5, 0.7}, 1.0); + expect_scalar_edge(ys_lcdf, INT_MAX_N, {1.5, 0.7}, 0.0); + expect_scalar_edge(ys_lccdf, 0, {1.5, 0.7}, 0.0); +} + +TEST(ProbDiscreteCdfsEdgeElements, neg_binomial_2_lcdf_expect_ad) { + // The infinite count is removed before the beta lcdf call + auto f = [](const auto& mu, const auto& phi) { + return stan::math::neg_binomial_2_lcdf(std::vector{3, INT_MAX_N}, mu, + phi); + }; + Eigen::VectorXd mu(2), phi(2); + mu << 4.0, 1.5; + phi << 2.0, 3.0; + stan::test::expect_ad(f, 4.0, 2.0); + stan::test::expect_ad(f, mu, 2.0); + stan::test::expect_ad(f, 4.0, phi); + stan::test::expect_ad(f, mu, phi); +} diff --git a/test/unit/math/prim/prob/poisson_binomial_ccdf_log_test.cpp b/test/unit/math/prim/prob/poisson_binomial_ccdf_log_test.cpp index a6ca3e9afcc..4a4d87e2e3b 100644 --- a/test/unit/math/prim/prob/poisson_binomial_ccdf_log_test.cpp +++ b/test/unit/math/prim/prob/poisson_binomial_ccdf_log_test.cpp @@ -38,7 +38,7 @@ TEST(ProbDistributionsPoissonBinomial, lccdf_works_on_scalar_arguments) { EXPECT_NEAR(-0.127833, poisson_binomial_lccdf(0, p), 1e-6); EXPECT_NEAR(-0.798508, poisson_binomial_lccdf(1, p), 1e-6); EXPECT_NEAR(-2.65926, poisson_binomial_lccdf(2, p), 1e-5); - EXPECT_NEAR(-36.7368, poisson_binomial_lccdf(3, p), 1e-4); + EXPECT_EQ(-stan::math::INFTY, poisson_binomial_lccdf(3, p)); } TEST(ProbDistributionsPoissonBinomial, lccdf_works_on_vectorial_y) { From a834863f5387919d372bb84574460eabbf1fac35 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sat, 31 Oct 2026 19:11:39 +0200 Subject: [PATCH 4/5] Test probability-0 edges and poisson_binomial theta per element with vector parameters --- .../prob/discrete_cdfs_edge_elements_test.cpp | 120 ++++++++++++++++++ 1 file changed, 120 insertions(+) diff --git a/test/unit/math/mix/prob/discrete_cdfs_edge_elements_test.cpp b/test/unit/math/mix/prob/discrete_cdfs_edge_elements_test.cpp index aea780522e4..7d78753ea4c 100644 --- a/test/unit/math/mix/prob/discrete_cdfs_edge_elements_test.cpp +++ b/test/unit/math/mix/prob/discrete_cdfs_edge_elements_test.cpp @@ -155,6 +155,23 @@ void expect_scalar_edge(const F& f, int edge, const std::vector& p, } } +// poisson_binomial: element e has the k probabilities p[j * m + e], +// j = 0, ..., k - 1, with m = p.size() / k elements +template +std::vector> make_theta(const std::vector& p, + size_t k) { + const size_t m = p.size() / k; + std::vector> theta; + for (size_t e = 0; e < m; ++e) { + Eigen::Matrix th(k); + for (size_t j = 0; j < k; ++j) { + th(j) = p[j * m + e]; + } + theta.push_back(th); + } + return theta; +} + constexpr double NEG_INF = -std::numeric_limits::infinity(); constexpr int INT_MAX_N = std::numeric_limits::max(); @@ -421,6 +438,109 @@ TEST(ProbDiscreteCdfsEdgeElements, vector_parameters) { expect_scalar_edge(ys_cdf, INT_MAX_N, {1.5, 0.7}, 1.0); expect_scalar_edge(ys_lcdf, INT_MAX_N, {1.5, 0.7}, 0.0); expect_scalar_edge(ys_lccdf, 0, {1.5, 0.7}, 0.0); + + // The binomial and beta_binomial cdf and lcdf with vector parameters + auto bin_cdf = [](const auto& n, const auto& p) { + return stan::math::binomial_cdf(n, 10, p); + }; + auto bin_lcdf = [](const auto& n, const auto& p) { + return stan::math::binomial_lcdf(n, 10, p); + }; + auto bb_cdf = [](const auto& n, const auto& p) { + return stan::math::beta_binomial_cdf(n, 10, part(p, 2, 0), part(p, 2, 1)); + }; + auto bb_lcdf = [](const auto& n, const auto& p) { + return stan::math::beta_binomial_lcdf(n, 10, part(p, 2, 0), part(p, 2, 1)); + }; + expect_probability_one_vec(bin_cdf, 3, 10, t_ord, t_edge); + expect_probability_one_vec(bin_lcdf, 3, 10, t_ord, t_edge); + expect_probability_one_vec(bb_cdf, 3, 10, ab_ord, ab_edge); + expect_probability_one_vec(bb_lcdf, 3, 10, ab_ord, ab_edge); + expect_scalar_edge(b_cdf, 1, {0.3, 0.6}, 1.0); + expect_scalar_edge(b_lcdf, 1, {0.3, 0.6}, 0.0); + expect_scalar_edge(bin_cdf, 10, {0.3, 0.6}, 1.0); + expect_scalar_edge(bin_lcdf, 10, {0.3, 0.6}, 0.0); + expect_scalar_edge(bb_cdf, 10, ab_both, 1.0); + expect_scalar_edge(bb_lcdf, 10, ab_both, 0.0); + + // Probability 0 with vector parameters + const std::vector t_both{0.3, 0.6}; + const std::vector mp_both{4.0, 1.5, 2.0, 3.0}; + const std::vector a_both{1.5, 0.7}; + expect_probability_zero(b_cdf, 0, -1, t_both, 0.0); + expect_probability_zero(b_lcdf, 0, -1, t_both, NEG_INF); + expect_probability_zero(bin_cdf, 3, -1, t_both, 0.0); + expect_probability_zero(bin_lcdf, 3, -1, t_both, NEG_INF); + expect_probability_zero(bb_cdf, 3, -1, ab_both, 0.0); + expect_probability_zero(bb_lcdf, 3, -1, ab_both, NEG_INF); + expect_probability_zero(nb_cdf, 3, -1, ab_both, 0.0); + expect_probability_zero(nb_lcdf, 3, -1, ab_both, NEG_INF); + expect_probability_zero(nb2_cdf, 3, -1, mp_both, 0.0); + expect_probability_zero(nb2_lcdf, 3, -1, mp_both, NEG_INF); + expect_probability_zero(nb2_lccdf, 3, INT_MAX_N, mp_both, NEG_INF); + expect_probability_zero(bnb_cdf, 3, -1, rab_both, 0.0); + expect_probability_zero(bnb_lcdf, 3, -1, rab_both, NEG_INF); + expect_probability_zero(ys_cdf, 3, 0, a_both, 0.0); + expect_probability_zero(ys_lcdf, 3, 0, a_both, NEG_INF); + expect_probability_zero(ys_lccdf, 3, INT_MAX_N, a_both, NEG_INF); + + // A scalar edge count with probability 0 and vector parameters + expect_scalar_edge(b_cdf, -1, t_both, 0.0); + expect_scalar_edge(b_lcdf, -1, t_both, NEG_INF); + expect_scalar_edge(b_lccdf, 1, t_both, NEG_INF); + expect_scalar_edge(bin_cdf, -1, t_both, 0.0); + expect_scalar_edge(bin_lcdf, -1, t_both, NEG_INF); + expect_scalar_edge(bin_lccdf, 10, t_both, NEG_INF); + expect_scalar_edge(bb_cdf, -1, ab_both, 0.0); + expect_scalar_edge(bb_lcdf, -1, ab_both, NEG_INF); + expect_scalar_edge(bb_lccdf, 10, ab_both, NEG_INF); + expect_scalar_edge(nb_cdf, -1, ab_both, 0.0); + expect_scalar_edge(nb_lcdf, -1, ab_both, NEG_INF); + expect_scalar_edge(nb_lccdf, INT_MAX_N, ab_both, NEG_INF); + expect_scalar_edge(nb2_cdf, -1, mp_both, 0.0); + expect_scalar_edge(nb2_lcdf, -1, mp_both, NEG_INF); + expect_scalar_edge(nb2_lccdf, INT_MAX_N, mp_both, NEG_INF); + expect_scalar_edge(bnb_cdf, -1, rab_both, 0.0); + expect_scalar_edge(bnb_lcdf, -1, rab_both, NEG_INF); + expect_scalar_edge(bnb_lccdf, INT_MAX_N, rab_both, NEG_INF); + expect_scalar_edge(ys_cdf, 0, a_both, 0.0); + expect_scalar_edge(ys_lcdf, 0, a_both, NEG_INF); + expect_scalar_edge(ys_lccdf, INT_MAX_N, a_both, NEG_INF); +} + +TEST(ProbDiscreteCdfsEdgeElements, poisson_binomial_theta_per_element) { + // Three trials per element + auto cdf = [](const auto& y, const auto& p) { + return stan::math::poisson_binomial_cdf(y, make_theta(p, 3)); + }; + auto lcdf = [](const auto& y, const auto& p) { + return stan::math::poisson_binomial_lcdf(y, make_theta(p, 3)); + }; + auto lccdf = [](const auto& y, const auto& p) { + return stan::math::poisson_binomial_lccdf(y, make_theta(p, 3)); + }; + const std::vector q_ord{0.3, 0.6, 0.8}, q_edge{0.2, 0.5, 0.9}; + const std::vector q_both{0.3, 0.2, 0.6, 0.5, 0.8, 0.9}; + expect_probability_one_vec(cdf, 1, 3, q_ord, q_edge); + expect_probability_one_vec(lcdf, 1, 3, q_ord, q_edge); + expect_probability_zero(lccdf, 1, 3, q_both, NEG_INF); + expect_scalar_edge(cdf, 3, q_both, 1.0); + expect_scalar_edge(lcdf, 3, q_both, 0.0); + expect_scalar_edge(lccdf, 3, q_both, NEG_INF); + + // One trial per element + auto cdf1 = [](const auto& y, const auto& p) { + return stan::math::poisson_binomial_cdf(y, make_theta(p, 1)); + }; + auto lcdf1 = [](const auto& y, const auto& p) { + return stan::math::poisson_binomial_lcdf(y, make_theta(p, 1)); + }; + auto lccdf1 = [](const auto& y, const auto& p) { + return stan::math::poisson_binomial_lccdf(y, make_theta(p, 1)); + }; + expect_probability_one_vec(cdf1, 0, 1, {0.3}, {0.6}); + expect_probability_one_vec(lcdf1, 0, 1, {0.3}, {0.6}); + expect_probability_zero(lccdf1, 0, 1, {0.3, 0.6}, NEG_INF); } TEST(ProbDiscreteCdfsEdgeElements, neg_binomial_2_lcdf_expect_ad) { From 1315143b74d4dc15f4c2fe7a53ed01c0fa721331 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sat, 31 Oct 2026 19:11:43 +0200 Subject: [PATCH 5/5] Test OpenCL bernoulli_lccdf broadcasting with negative counts --- .../math/opencl/rev/bernoulli_lccdf_test.cpp | 21 +++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/test/unit/math/opencl/rev/bernoulli_lccdf_test.cpp b/test/unit/math/opencl/rev/bernoulli_lccdf_test.cpp index 732b1e62f2c..d2aba277fc7 100644 --- a/test/unit/math/opencl/rev/bernoulli_lccdf_test.cpp +++ b/test/unit/math/opencl/rev/bernoulli_lccdf_test.cpp @@ -106,6 +106,27 @@ TEST(ProbDistributionsBernoulliLccdf, opencl_broadcast_theta) { bernoulli_lccdf_functor, n, theta_scal); } +TEST(ProbDistributionsBernoulliLccdf, opencl_broadcast_n_negative) { + int N = 3; + + int n_scal = -1; + Eigen::VectorXd theta(N); + theta << 0.3, 0.8, 1.0; + + stan::math::test::test_opencl_broadcasting_prim_rev<0>( + bernoulli_lccdf_functor, n_scal, theta); +} + +TEST(ProbDistributionsBernoulliLccdf, opencl_broadcast_theta_n_negative) { + int N = 3; + + std::vector n{0, -1, 0}; + double theta_scal = 0.4; + + stan::math::test::test_opencl_broadcasting_prim_rev<1>( + bernoulli_lccdf_functor, n, theta_scal); +} + TEST(ProbDistributionsBernoulliLccdf, opencl_matches_cpu_big) { int N = 153;