From 0b78413c024d86d7a2531a8b8165773865a22aa5 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sat, 10 Oct 2026 10:07:48 +0300 Subject: [PATCH] Compute the beta_neg_binomial cdfs from sums of the pmf instead of the 3F2 series --- stan/math/prim/prob.hpp | 1 + stan/math/prim/prob/beta_neg_binomial_cdf.hpp | 100 +--- .../prim/prob/beta_neg_binomial_lccdf.hpp | 84 +-- .../math/prim/prob/beta_neg_binomial_lcdf.hpp | 92 +--- .../prim/prob/beta_neg_binomial_log_cdfs.hpp | 498 ++++++++++++++++++ .../rev/prob/beta_neg_binomial_cdfs_test.cpp | 155 ++++++ test/unit/math/rev/prob/large_shapes_test.cpp | 8 +- 7 files changed, 736 insertions(+), 202 deletions(-) create mode 100644 stan/math/prim/prob/beta_neg_binomial_log_cdfs.hpp create mode 100644 test/unit/math/rev/prob/beta_neg_binomial_cdfs_test.cpp diff --git a/stan/math/prim/prob.hpp b/stan/math/prim/prob.hpp index 4fae64c1ace..7e3995723d0 100644 --- a/stan/math/prim/prob.hpp +++ b/stan/math/prim/prob.hpp @@ -28,6 +28,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/prob/beta_neg_binomial_cdf.hpp b/stan/math/prim/prob/beta_neg_binomial_cdf.hpp index e39c7385e8c..a5f13cd4fc2 100644 --- a/stan/math/prim/prob/beta_neg_binomial_cdf.hpp +++ b/stan/math/prim/prob/beta_neg_binomial_cdf.hpp @@ -4,20 +4,14 @@ #include #include #include -#include -#include #include -#include -#include -#include -#include -#include #include #include #include #include #include -#include +#include +#include namespace stan { namespace math { @@ -27,6 +21,9 @@ namespace math { * number of successes, prior success, and prior failure parameters. * Given containers of matching sizes, returns the product of probabilities. * + * The product is the exponential of the sum of the log cdfs of + * internal::beta_neg_binomial_log_cdfs(), so a small cdf keeps its digits. + * * @tparam T_n type of failure parameter * @tparam T_r type of number of successes parameter * @tparam T_alpha type of prior success parameter @@ -36,10 +33,12 @@ namespace math { * @param r Number of successes parameter * @param alpha prior success parameter * @param beta prior failure parameter - * @param precision precision for `grad_F32`, default \f$10^{-8}\f$ - * @param max_steps max iteration allowed for `grad_F32`, default \f$10^{8}\f$ - * @return probability or sum of probabilities - * @throw std::domain_error if r, alpha, or beta fails to be positive + * @param precision not used; the sums stop when their remaining terms are + * below 1e-17 of the partial sums + * @param max_steps largest number of terms of a sum, default \f$10^{8}\f$ + * @return probability or product of probabilities + * @throw std::domain_error if r, alpha, or beta fails to be positive, or if + * a sum needs more than max_steps terms * @throw std::invalid_argument if container sizes mismatch */ template @@ -80,72 +79,31 @@ inline return_type_t beta_neg_binomial_cdf( } using T_partials_return = partials_return_t; - T_partials_return cdf(1.0); + constexpr bool any_autodiff = is_any_autodiff_v; + 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 + // The largest int stands for infinity: its factor is 1 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); - auto alpha_dbl = alpha_vec.val(i); - auto beta_dbl = beta_vec.val(i); - auto b_plus_n = beta_dbl + n_dbl; - auto r_plus_n = r_dbl + n_dbl; - auto a_plus_r = alpha_dbl + r_dbl; - using a_t = return_type_t; - using b_t = return_type_t; - auto F = hypergeometric_3F2( - std::initializer_list{1.0, b_plus_n + 1.0, r_plus_n + 1.0}, - std::initializer_list{n_dbl + 2.0, a_plus_r + b_plus_n + 1.0}, - 1.0); - // C is the log pmf at n + 1; formed from lgamma and lbeta directly it - // cancels for large parameters - const T_partials_return k = n_dbl + 1.0; - const T_partials_return C = internal::beta_neg_binomial_log_pmf( - k, T_partials_return(r_dbl), T_partials_return(alpha_dbl), - T_partials_return(beta_dbl)); - auto ccdf = stan::math::exp(C + stan::math::log(F)); - cdf *= 1.0 - ccdf; - - if constexpr (is_any_autodiff_v) { - auto chain_rule_term = -ccdf / (1.0 - ccdf); - T_partials_return dF[6]; - grad_F32, is_autodiff_v, false, true, - false>(dF, 1.0, b_plus_n + 1.0, r_plus_n + 1.0, n_dbl + 2.0, - a_plus_r + b_plus_n + 1.0, 1.0, precision, max_steps); - // the partials of C as differences digamma(x + d) - digamma(x) - const T_partials_return alpha_plus_beta = alpha_dbl + beta_dbl; - const T_partials_return dpsi_r_plus_ab_k - = digamma_diff(r_dbl + alpha_plus_beta, k); - if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[i] - += (digamma_diff(r_dbl, k) - dpsi_r_plus_ab_k - - digamma_diff(a_plus_r, beta_dbl) + (dF[2] + dF[4]) / F) - * chain_rule_term; - } - if constexpr (is_any_autodiff_v) { - const T_partials_return dpsi_ab_r - = digamma_diff(alpha_plus_beta, r_dbl); - if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[i] - += (digamma_diff(alpha_dbl, r_dbl) - dpsi_ab_r - dpsi_r_plus_ab_k - + dF[4] / F) - * chain_rule_term; - } - if constexpr (is_autodiff_v) { - partials<2>(ops_partials)[i] - += (digamma_diff(beta_dbl, k) - dpsi_r_plus_ab_k - dpsi_ab_r - + (dF[1] + dF[4]) / F) - * chain_rule_term; - } - } + const auto res = internal::beta_neg_binomial_log_cdfs( + n_vec.val(i), T_partials_return(r_vec.val(i)), + T_partials_return(alpha_vec.val(i)), T_partials_return(beta_vec.val(i)), + max_steps, function); + log_cdf += res.lcdf; + if constexpr (is_autodiff_v) { + partials<0>(ops_partials)[i] += res.dlcdf[0]; + } + if constexpr (is_autodiff_v) { + partials<1>(ops_partials)[i] += res.dlcdf[1]; + } + if constexpr (is_autodiff_v) { + partials<2>(ops_partials)[i] += res.dlcdf[2]; } } + const T_partials_return cdf = exp(log_cdf); if constexpr (is_autodiff_v) { for (size_t i = 0; i < stan::math::size(r); ++i) { partials<0>(ops_partials)[i] *= cdf; diff --git a/stan/math/prim/prob/beta_neg_binomial_lccdf.hpp b/stan/math/prim/prob/beta_neg_binomial_lccdf.hpp index 259fddff2b8..0edce58b93d 100644 --- a/stan/math/prim/prob/beta_neg_binomial_lccdf.hpp +++ b/stan/math/prim/prob/beta_neg_binomial_lccdf.hpp @@ -4,18 +4,13 @@ #include #include #include -#include -#include -#include -#include -#include -#include #include #include #include #include #include -#include +#include +#include namespace stan { namespace math { @@ -25,6 +20,9 @@ namespace math { * number of successes, prior success, and prior failure parameters. * Given containers of matching sizes, returns the log sum of probabilities. * + * The lower or the upper tail is summed directly, and the other one is its + * complement; see internal::beta_neg_binomial_log_cdfs(). + * * @tparam T_n type of failure parameter * @tparam T_r type of number of successes parameter * @tparam T_alpha type of prior success parameter @@ -34,10 +32,12 @@ namespace math { * @param r Number of successes parameter * @param alpha prior success parameter * @param beta prior failure parameter - * @param precision precision for `grad_F32`, default \f$10^{-8}\f$ - * @param max_steps max iteration allowed for `grad_F32`, default \f$10^{8}\f$ + * @param precision not used; the sums stop when their remaining terms are + * below 1e-17 of the partial sums + * @param max_steps largest number of terms of a sum, default \f$10^{8}\f$ * @return log probability or log sum of probabilities - * @throw std::domain_error if r, alpha, or beta fails to be positive + * @throw std::domain_error if r, alpha, or beta fails to be positive, or if + * a sum needs more than max_steps terms * @throw std::invalid_argument if container sizes mismatch */ template @@ -78,6 +78,7 @@ inline return_type_t beta_neg_binomial_lccdf( } using T_partials_return = partials_return_t; + constexpr bool any_autodiff = is_any_autodiff_v; 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++) { @@ -86,56 +87,19 @@ inline return_type_t beta_neg_binomial_lccdf( if (n_vec.val(i) == std::numeric_limits::max()) { return ops_partials.build(negative_infinity()); } - auto n_dbl = n_vec.val(i); - auto r_dbl = r_vec.val(i); - auto alpha_dbl = alpha_vec.val(i); - auto beta_dbl = beta_vec.val(i); - auto b_plus_n = beta_dbl + n_dbl; - auto r_plus_n = r_dbl + n_dbl; - auto a_plus_r = alpha_dbl + r_dbl; - using a_t = return_type_t; - using b_t = return_type_t; - auto F = hypergeometric_3F2( - std::initializer_list{1.0, b_plus_n + 1.0, r_plus_n + 1.0}, - std::initializer_list{n_dbl + 2.0, a_plus_r + b_plus_n + 1.0}, - 1.0); - // C is the log pmf at n + 1; formed from lgamma and lbeta directly it - // cancels for large parameters - const T_partials_return k = n_dbl + 1.0; - const T_partials_return C = internal::beta_neg_binomial_log_pmf( - k, T_partials_return(r_dbl), T_partials_return(alpha_dbl), - T_partials_return(beta_dbl)); - log_ccdf += C + stan::math::log(F); - - if constexpr (is_any_autodiff_v) { - T_partials_return dF[6]; - grad_F32, is_autodiff_v, false, true, - false>(dF, 1.0, b_plus_n + 1.0, r_plus_n + 1.0, n_dbl + 2.0, - a_plus_r + b_plus_n + 1.0, 1.0, precision, max_steps); - // the partials of C as differences digamma(x + d) - digamma(x) - const T_partials_return alpha_plus_beta = alpha_dbl + beta_dbl; - const T_partials_return dpsi_r_plus_ab_k - = digamma_diff(r_dbl + alpha_plus_beta, k); - if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[i] - += digamma_diff(r_dbl, k) - dpsi_r_plus_ab_k - - digamma_diff(a_plus_r, beta_dbl) + (dF[2] + dF[4]) / F; - } - if constexpr (is_any_autodiff_v) { - const T_partials_return dpsi_ab_r - = digamma_diff(alpha_plus_beta, r_dbl); - if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[i] += digamma_diff(alpha_dbl, r_dbl) - - dpsi_ab_r - dpsi_r_plus_ab_k - + dF[4] / F; - } - if constexpr (is_autodiff_v) { - partials<2>(ops_partials)[i] += digamma_diff(beta_dbl, k) - - dpsi_r_plus_ab_k - dpsi_ab_r - + (dF[1] + dF[4]) / F; - } - } + const auto res = internal::beta_neg_binomial_log_cdfs( + n_vec.val(i), T_partials_return(r_vec.val(i)), + T_partials_return(alpha_vec.val(i)), T_partials_return(beta_vec.val(i)), + max_steps, function); + log_ccdf += res.lccdf; + if constexpr (is_autodiff_v) { + partials<0>(ops_partials)[i] += res.dlccdf[0]; + } + if constexpr (is_autodiff_v) { + partials<1>(ops_partials)[i] += res.dlccdf[1]; + } + if constexpr (is_autodiff_v) { + partials<2>(ops_partials)[i] += res.dlccdf[2]; } } diff --git a/stan/math/prim/prob/beta_neg_binomial_lcdf.hpp b/stan/math/prim/prob/beta_neg_binomial_lcdf.hpp index 6bd1cedc1f4..87bfc9502aa 100644 --- a/stan/math/prim/prob/beta_neg_binomial_lcdf.hpp +++ b/stan/math/prim/prob/beta_neg_binomial_lcdf.hpp @@ -4,18 +4,13 @@ #include #include #include -#include -#include -#include -#include -#include -#include #include #include #include #include #include -#include +#include +#include namespace stan { namespace math { @@ -25,6 +20,9 @@ namespace math { * number of successes, prior success, and prior failure parameters. * Given containers of matching sizes, returns the log sum of probabilities. * + * The lower or the upper tail is summed directly, and the other one is its + * complement; see internal::beta_neg_binomial_log_cdfs(). + * * @tparam T_n type of failure parameter * @tparam T_r type of number of successes parameter * @tparam T_alpha type of prior success parameter @@ -34,10 +32,12 @@ namespace math { * @param r Number of successes parameter * @param alpha prior success parameter * @param beta prior failure parameter - * @param precision precision for `grad_F32`, default \f$10^{-8}\f$ - * @param max_steps max iteration allowed for `grad_F32`, default \f$10^{8}\f$ + * @param precision not used; the sums stop when their remaining terms are + * below 1e-17 of the partial sums + * @param max_steps largest number of terms of a sum, default \f$10^{8}\f$ * @return log probability or log sum of probabilities - * @throw std::domain_error if r, alpha, or beta fails to be positive + * @throw std::domain_error if r, alpha, or beta fails to be positive, or if + * a sum needs more than max_steps terms * @throw std::invalid_argument if container sizes mismatch */ template @@ -78,67 +78,27 @@ inline return_type_t beta_neg_binomial_lcdf( } using T_partials_return = partials_return_t; + constexpr bool any_autodiff = is_any_autodiff_v; 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 + // The largest int stands for infinity: its term is log(1) = 0 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); - auto alpha_dbl = alpha_vec.val(i); - auto beta_dbl = beta_vec.val(i); - auto b_plus_n = beta_dbl + n_dbl; - auto r_plus_n = r_dbl + n_dbl; - auto a_plus_r = alpha_dbl + r_dbl; - using a_t = return_type_t; - using b_t = return_type_t; - auto F = hypergeometric_3F2( - std::initializer_list{1.0, b_plus_n + 1.0, r_plus_n + 1.0}, - std::initializer_list{n_dbl + 2.0, a_plus_r + b_plus_n + 1.0}, - 1.0); - // C is the log pmf at n + 1; formed from lgamma and lbeta directly it - // cancels for large parameters - const T_partials_return k = n_dbl + 1.0; - const T_partials_return C = internal::beta_neg_binomial_log_pmf( - k, T_partials_return(r_dbl), T_partials_return(alpha_dbl), - T_partials_return(beta_dbl)); - auto ccdf = stan::math::exp(C) * F; - log_cdf += log1m(ccdf); - - if constexpr (is_any_autodiff_v) { - auto chain_rule_term = -ccdf / (1.0 - ccdf); - T_partials_return dF[6]; - grad_F32, is_autodiff_v, false, true, - false>(dF, 1.0, b_plus_n + 1.0, r_plus_n + 1.0, n_dbl + 2.0, - a_plus_r + b_plus_n + 1.0, 1.0, precision, max_steps); - // the partials of C as differences digamma(x + d) - digamma(x) - const T_partials_return alpha_plus_beta = alpha_dbl + beta_dbl; - const T_partials_return dpsi_r_plus_ab_k - = digamma_diff(r_dbl + alpha_plus_beta, k); - if constexpr (is_autodiff_v) { - auto partial_lccdf = digamma_diff(r_dbl, k) - dpsi_r_plus_ab_k - - digamma_diff(a_plus_r, beta_dbl) - + (dF[2] + dF[4]) / F; - partials<0>(ops_partials)[i] += partial_lccdf * chain_rule_term; - } - if constexpr (is_any_autodiff_v) { - const T_partials_return dpsi_ab_r - = digamma_diff(alpha_plus_beta, r_dbl); - if constexpr (is_autodiff_v) { - auto partial_lccdf = digamma_diff(alpha_dbl, r_dbl) - dpsi_ab_r - - dpsi_r_plus_ab_k + dF[4] / F; - partials<1>(ops_partials)[i] += partial_lccdf * chain_rule_term; - } - if constexpr (is_autodiff_v) { - auto partial_lccdf = digamma_diff(beta_dbl, k) - dpsi_r_plus_ab_k - - dpsi_ab_r + (dF[1] + dF[4]) / F; - partials<2>(ops_partials)[i] += partial_lccdf * chain_rule_term; - } - } + const auto res = internal::beta_neg_binomial_log_cdfs( + n_vec.val(i), T_partials_return(r_vec.val(i)), + T_partials_return(alpha_vec.val(i)), T_partials_return(beta_vec.val(i)), + max_steps, function); + log_cdf += res.lcdf; + if constexpr (is_autodiff_v) { + partials<0>(ops_partials)[i] += res.dlcdf[0]; + } + if constexpr (is_autodiff_v) { + partials<1>(ops_partials)[i] += res.dlcdf[1]; + } + if constexpr (is_autodiff_v) { + partials<2>(ops_partials)[i] += res.dlcdf[2]; } } diff --git a/stan/math/prim/prob/beta_neg_binomial_log_cdfs.hpp b/stan/math/prim/prob/beta_neg_binomial_log_cdfs.hpp new file mode 100644 index 00000000000..60b6bd3d42a --- /dev/null +++ b/stan/math/prim/prob/beta_neg_binomial_log_cdfs.hpp @@ -0,0 +1,498 @@ +#ifndef STAN_MATH_PRIM_PROB_BETA_NEG_BINOMIAL_LOG_CDFS_HPP +#define STAN_MATH_PRIM_PROB_BETA_NEG_BINOMIAL_LOG_CDFS_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { +namespace internal { + +/** + * The log cdf and the log ccdf of the beta negative binomial distribution at + * one point, and their derivatives with respect to (r, alpha, beta). + */ +template +struct beta_neg_binomial_log_cdfs_t { + T lcdf; + T lccdf; + std::array dlcdf; + std::array dlccdf; +}; + +// A sum stops when a bound of its remaining terms is below this fraction of +// the partial sum. +constexpr double bnb_cdfs_tol = 1e-17; +// The complement 1 - p is used for a probability p of at least this size; +// it multiplies the relative error of p by at most (1 - p) / p < 100. +constexpr double bnb_cdfs_complement_min = 1e-2; +constexpr int bnb_cdfs_short_tail_steps = 2000; +// Up to this n the lower sum, with at most n + 1 terms, comes first. +constexpr int bnb_cdfs_lower_first_max_n = 100; +constexpr int bnb_cdfs_thomae_steps = 5000; +constexpr double bnb_cdfs_thomae_max_cancellation = 100; + +template +inline double bnb_cdfs_max_abs(const std::array& g) { + return std::max({std::fabs(value_of_rec(g[0])), std::fabs(value_of_rec(g[1])), + std::fabs(value_of_rec(g[2]))}); +} + +/** + * Writes d log pmf(k) / d(r, alpha, beta) to g, as differences of digamma + * functions that do not cancel at large arguments. + */ +template +inline void bnb_cdfs_dlog_pmf(const T& k, const T& r, const T& alpha, + const T& beta, std::array& g) { + const T dab = digamma_diff(T(alpha + beta), T(k + r)); + g[0] = digamma_diff(r, k) - digamma_diff(T(r + alpha), T(k + beta)); + g[1] = digamma_diff(alpha, r) - dab; + g[2] = digamma_diff(beta, k) - dab; +} + +/** + * Computes log cdf(n) as the sum of the pmf from k = n down to k = 0 in log + * scale, and its gradient as sum pmf(k) d log pmf(k) / cdf(n). All terms are + * positive. pmf(k + 1) / pmf(k) - 1 has the sign of + * r beta - s - k (1 + alpha), s = r + alpha + beta, so the pmf is unimodal: + * below the mode the terms fall as k falls, the j remaining terms are each + * below the current one, and the sum stops when they are negligible. + */ +template +inline void bnb_cdfs_lower_sum(int n, const T& r, const T& alpha, const T& beta, + int max_steps, const char* function, T& log_cdf, + std::array& dlog) { + const T s = r + alpha + beta; + const double k_peak = value_of_rec((r * beta - s) / (1.0 + alpha)); + const double smaller = std::min(value_of_rec(r), value_of_rec(beta)); + T log_scale = beta_neg_binomial_log_pmf(T(n), r, alpha, beta); + std::array g{0, 0, 0}; + std::array d{0, 0, 0}; + if (Grad) { + bnb_cdfs_dlog_pmf(T(n), r, alpha, beta, g); + } + T t = 1.0; + T sum = 0.0; + constexpr double big = 1e250; + const double log_big = std::log(big); + for (int k = n, steps = 0;; --k, ++steps) { + if (steps > max_steps) { + throw_domain_error(function, "number of terms of the lower sum", + max_steps, "exceeded ", ""); + } + sum += t; + if (Grad) { + for (int i = 0; i < 3; ++i) { + d[i] += t * g[i]; + } + } + if (k == 0) { + break; + } + const double j = k - 1.0; + // pmf(j) / pmf(j + 1) + const T q = (j + 1.0) * (j + s) / ((j + r) * (j + beta)); + if (j < k_peak) { + // the j + 1 terms below are each below t q; on the way down + // |d log pmf| changes by at most 3 (log1p(j / smaller) + 1 / smaller) + const double rest = value_of_rec((j + 1.0) * t * q); + const double g_bound + = Grad ? bnb_cdfs_max_abs(g) + + 3.0 * (std::log1p(j / smaller) + 1.0 / smaller) + 1.0 + : 1.0; + if (rest * g_bound <= bnb_cdfs_tol * value_of_rec(sum)) { + break; + } + } + if (Grad) { + const T js = j + s; + g[0] -= (alpha + beta) / ((j + r) * js); + g[1] += 1.0 / js; + g[2] -= (r + alpha) / ((j + beta) * js); + } + t *= q; + if (value_of_rec(t) > big || value_of_rec(sum) > big) { + t /= big; + sum /= big; + if (Grad) { + for (int i = 0; i < 3; ++i) { + d[i] /= big; + } + } + log_scale += log_big; + } + } + log_cdf = log_scale + log(sum); + if (Grad) { + for (int i = 0; i < 3; ++i) { + dlog[i] = d[i] / sum; + } + } +} + +/** + * Computes log ccdf(n) as the sum of the pmf over k > n, and its gradient. + * Past the peak of the terms the rest is about t / c for a geometric fall, + * c = 1 - pmf(k + 1) / pmf(k), and about t (k + s) / alpha for a fall like + * k^(-1 - alpha); the larger estimate is used. Returns false if the sum + * does not stop within max_terms terms. + */ +template +inline bool bnb_cdfs_tail_sum(int n, const T& r, const T& alpha, const T& beta, + int max_terms, T& log_ccdf, + std::array& dlog) { + const T s = r + alpha + beta; + const double k_peak = value_of_rec((r * beta - s) / (1.0 + alpha)); + const double alpha_dbl = value_of_rec(alpha); + double k = n + 1.0; + const T log_p = beta_neg_binomial_log_pmf(T(k), r, alpha, beta); + std::array g{0, 0, 0}; + std::array d{0, 0, 0}; + if (Grad) { + bnb_cdfs_dlog_pmf(T(k), r, alpha, beta, g); + } + T t = 1.0; + T sum = 0.0; + for (int terms = 1; terms <= max_terms; ++terms, k += 1.0) { + sum += t; + if (Grad) { + for (int i = 0; i < 3; ++i) { + d[i] += t * g[i]; + } + } + const T ks = k + s; + if (k > k_peak) { + const double c + = value_of_rec((k * (1.0 + alpha) + s - r * beta) / ((k + 1.0) * ks)); + const double rest + = value_of_rec(t) * std::max(1.0 / c, value_of_rec(ks) / alpha_dbl); + const double g_bound + = Grad ? bnb_cdfs_max_abs(g) + 1.0 + 1.0 / alpha_dbl : 1.0; + if (rest * g_bound <= bnb_cdfs_tol * value_of_rec(sum)) { + log_ccdf = log_p + log(sum); + if (Grad) { + for (int i = 0; i < 3; ++i) { + dlog[i] = d[i] / sum; + } + } + return true; + } + } + t *= (k + r) * (k + beta) / ((k + 1.0) * ks); + if (Grad) { + g[0] += (alpha + beta) / ((k + r) * ks); + g[1] -= 1.0 / ks; + g[2] += (r + alpha) / ((k + beta) * ks); + } + } + return false; +} + +/** + * Returns whether bnb_cdfs_thomae() takes x = r (else x = beta): the smaller + * of x (alpha + y) gives the smaller first terms. + */ +inline bool bnb_cdfs_thomae_x_is_r(double r, double alpha, double beta) { + return r * (alpha + beta) <= beta * (alpha + r); +} + +/** + * Returns an estimate of the number of terms of bnb_cdfs_tail_sum() for a + * geometric fall: log(1 / bnb_cdfs_tol) / c0, c0 = 1 - pmf(n + 2) / + * pmf(n + 1), plus the steps up to the peak of the terms. The estimate is + * optimistic where the terms later fall like k^(-1 - alpha); it only + * selects the method. + */ +inline double bnb_cdfs_tail_sum_terms(int n, double r, double alpha, + double beta) { + const double s = r + alpha + beta; + const double k_peak = (r * beta - s) / (1.0 + alpha); + const double k = n + 1.0; + const double log_tol = -std::log(bnb_cdfs_tol); + if (k <= k_peak) { + return k_peak - k + log_tol; + } + const double c0 = (k * (1.0 + alpha) + s - r * beta) / ((k + 1.0) * (k + s)); + return (c0 > 0.0) ? log_tol / c0 : INFTY; +} + +/** + * Returns an estimate of the number of terms of bnb_cdfs_thomae(), or + * infinity if it exceeds bnb_cdfs_thomae_steps. Over K terms the terms of + * the series fall by about x log(K) + log1p(K / alpha) + (n + 1) log1p(K / + * A); the series stops when the fall reaches log(1 / bnb_cdfs_tol) plus the + * log of the factor max(1, A / (x + n + 1)) of the rest estimate. The + * estimate only selects the method. + */ +inline double bnb_cdfs_thomae_terms(int n, double r, double alpha, + double beta) { + const bool x_is_r = bnb_cdfs_thomae_x_is_r(r, alpha, beta); + const double x = x_is_r ? r : beta; + const double y = x_is_r ? beta : r; + const double n1 = n + 1.0; + const double a = alpha + y + n1; + const double needed + = -std::log(bnb_cdfs_tol) + std::log(std::max(1.0, a / (x + n1))); + for (double k = 16; k <= bnb_cdfs_thomae_steps; k *= 4) { + const double fall + = x * std::log(k) + std::log1p(k / alpha) + n1 * std::log1p(k / a); + if (fall >= needed) { + return k; + } + } + return INFTY; +} + +/** + * Computes log ccdf(n) and its gradient by a Thomae transformation of the + * series in ccdf = pmf(n + 1) F, F = 3F2(1, beta + n + 1, r + n + 1; + * n + 2, alpha + r + beta + n + 1; 1), whose terms fall only like + * k^(-1 - alpha). With x one of r and beta, y the other one, and + * A = alpha + y + n + 1: + * + * F = (n + 1) / alpha * (A)_x / (n + 1)_x + * * 3F2(1 - x, alpha + y, alpha; alpha + 1, A; 1), + * + * where the new series has excess x + n + 1. For x > 1 its terms change + * sign, so the function returns false if the sum of their absolute values + * exceeds the sum by more than bnb_cdfs_thomae_max_cancellation, or if the + * series does not stop within bnb_cdfs_thomae_steps terms. For an integer + * x the factor 1 - x + k vanishes at k = x - 1: the later terms are 0, but + * their derivatives with respect to x are not. + */ +template +inline bool bnb_cdfs_thomae(int n, const T& r, const T& alpha, const T& beta, + T& log_ccdf, std::array& dlog) { + const bool x_is_r = bnb_cdfs_thomae_x_is_r( + value_of_rec(r), value_of_rec(alpha), value_of_rec(beta)); + const T& x = x_is_r ? r : beta; + const T& y = x_is_r ? beta : r; + const double n1 = n + 1.0; + const T a = alpha + y + n1; + const double x_plus_n1 = value_of_rec(x) + n1; + T t = 1.0; + T sum = 0.0; + double abs_sum = 0.0; + // derivatives of the sum and of log t with respect to (x, y, alpha) + std::array dsum{0, 0, 0}; + std::array h{0, 0, 0}; + bool polynomial = false; + T u = 0.0; + bool done = false; + for (int k = 0; k < bnb_cdfs_thomae_steps; ++k) { + const T a1 = 1.0 - x + k; + const T a2 = alpha + y + k; + const T a3 = alpha + k; + const T b1 = alpha + 1.0 + k; + const T b2 = a + k; + const T others = a2 * a3 / (b1 * b2 * (k + 1.0)); + if (polynomial) { + // d t_k / dx = -u_k after the zero factor + dsum[0] -= u; + const double abs_q = std::fabs(value_of_rec(a1 * others)); + if (abs_q < 1.0) { + const double rest + = std::fabs(value_of_rec(u)) * abs_q + * std::max(1.0 / (1.0 - abs_q), value_of_rec(b2) / x_plus_n1); + if (rest <= bnb_cdfs_tol + * (std::fabs(value_of_rec(sum)) + + std::fabs(value_of_rec(dsum[0])))) { + done = true; + break; + } + } + u *= a1 * others; + continue; + } + sum += t; + abs_sum += std::fabs(value_of_rec(t)); + if (Grad) { + for (int i = 0; i < 3; ++i) { + dsum[i] += t * h[i]; + } + } + if (value_of_rec(a1) == 0.0) { + if (!Grad) { + done = true; + break; + } + polynomial = true; + u = t * others; + continue; + } + const T q = a1 * others; + const double abs_q = std::fabs(value_of_rec(q)); + if (value_of_rec(a1) > 0.0 && abs_q < 1.0) { + const double rest + = std::fabs(value_of_rec(t)) * abs_q + * std::max(1.0 / (1.0 - abs_q), value_of_rec(b2) / x_plus_n1); + const double abs_partial = std::fabs(value_of_rec(sum)); + // no sign changes remain: the sum cannot grow by more than rest + if (abs_sum > bnb_cdfs_thomae_max_cancellation * (abs_partial + rest)) { + return false; + } + const double h_bound = Grad ? bnb_cdfs_max_abs(h) + 1.0 : 1.0; + if (rest * h_bound <= bnb_cdfs_tol * abs_partial) { + done = true; + break; + } + } + if (Grad) { + h[0] -= 1.0 / a1; + h[1] += n1 / (a2 * b2); + h[2] += n1 / (a2 * b2) + 1.0 / (a3 * b1); + } + t *= q; + } + if (!done || !(value_of_rec(sum) > 0.0) + || abs_sum > bnb_cdfs_thomae_max_cancellation * value_of_rec(sum)) { + return false; + } + log_ccdf = beta_neg_binomial_log_pmf(T(n1), r, alpha, beta) + std::log(n1) + - log_rising_factorial_ratio(T(n1), T(alpha + y), x) - log(alpha) + + log(sum); + if (Grad) { + std::array g; + bnb_cdfs_dlog_pmf(T(n1), r, alpha, beta, g); + const T d_y_prefactor = digamma_diff(a, x); + const T d_x = digamma_diff(T(x + n1), T(alpha + y)) + dsum[0] / sum; + const T d_y = d_y_prefactor + dsum[1] / sum; + dlog[0] = g[0] + (x_is_r ? d_x : d_y); + dlog[1] = g[1] + d_y_prefactor - 1.0 / alpha + dsum[2] / sum; + dlog[2] = g[2] + (x_is_r ? d_y : d_x); + } + return true; +} + +/** + * Returns the log cdf and the log ccdf of the beta negative binomial + * distribution at n >= 0, and, if Grad, their derivatives with respect to + * (r, alpha, beta). One of the two tails is summed directly and the other + * one is its complement, so that no result loses digits to cancellation: + * + * - past the mode and for n > bnb_cdfs_lower_first_max_n, the upper tail + * first: a short sum of the pmf over k > n, or a Thomae transformation + * of its series (fast for large n), the one with the smaller estimate + * of its number of terms first; + * - else the lower tail, the sum of the pmf from n down to 0; + * - the complement where the computed tail is at least 1/2, or the other + * tail is at least bnb_cdfs_complement_min; + * - else the upper tail by the short sum, the Thomae series or a long sum. + * + * @tparam Grad whether to compute the derivatives + * @tparam T type of the parameters + * @param n outcome, at least 0 + * @param r number of successes, positive + * @param alpha prior success, positive + * @param beta prior failure, positive + * @param max_steps largest number of terms of a sum + * @param function name of the calling function, for the error message + * @return log cdf, log ccdf and, if Grad, their derivatives + * @throw std::domain_error if a sum needs more than max_steps terms + */ +template +inline beta_neg_binomial_log_cdfs_t beta_neg_binomial_log_cdfs( + int n, const T& r, const T& alpha, const T& beta, int max_steps, + const char* function) { + beta_neg_binomial_log_cdfs_t out{0, 0, {0, 0, 0}, {0, 0, 0}}; + const double k_peak + = value_of_rec((r * beta - (r + alpha + beta)) / (1.0 + alpha)); + const double log_half = -LOG_TWO; + const double log_complement_min = std::log(bnb_cdfs_complement_min); + T lc = 0.0; + std::array d{0, 0, 0}; + bool have_upper = false; + auto from_upper = [&]() { + out.lccdf = lc; + out.lcdf = log1m_exp(lc); + if (Grad) { + const T ratio = exp(lc - out.lcdf); + for (int i = 0; i < 3; ++i) { + out.dlccdf[i] = d[i]; + out.dlcdf[i] = -d[i] * ratio; + } + } + return out; + }; + // the fast upper-tail methods, the one with the smaller estimate of its + // number of terms first; the tail-sum estimate is optimistic for a fall + // like k^(-1 - alpha), hence the factor 4. A method whose estimate exceeds + // its budget is not tried. + auto upper_fast = [&]() { + const double r_dbl = value_of_rec(r); + const double alpha_dbl = value_of_rec(alpha); + const double beta_dbl = value_of_rec(beta); + const double tail_terms + = bnb_cdfs_tail_sum_terms(n, r_dbl, alpha_dbl, beta_dbl); + const double thomae_terms + = bnb_cdfs_thomae_terms(n, r_dbl, alpha_dbl, beta_dbl); + auto tail = [&]() { + return tail_terms <= bnb_cdfs_short_tail_steps + && bnb_cdfs_tail_sum(n, r, alpha, beta, + bnb_cdfs_short_tail_steps, lc, d); + }; + auto thomae = [&]() { + return thomae_terms <= bnb_cdfs_thomae_steps + && bnb_cdfs_thomae(n, r, alpha, beta, lc, d); + }; + if (thomae_terms < 4.0 * tail_terms) { + return thomae() || tail(); + } + return tail() || thomae(); + }; + // the lower sum has at most n + 1 terms, and fewer below the mode + const bool tried_fast = n > k_peak && n > bnb_cdfs_lower_first_max_n; + if (tried_fast) { + have_upper = upper_fast(); + // the complement of the upper tail is accurate if the lower tail is + // not small + if (have_upper + && (value_of_rec(lc) <= log_half + || value_of_rec(log1m_exp(lc)) >= log_complement_min)) { + return from_upper(); + } + } + bnb_cdfs_lower_sum(n, r, alpha, beta, max_steps, function, out.lcdf, + out.dlcdf); + const double lcdf_dbl = value_of_rec(out.lcdf); + if (lcdf_dbl < log_half || -std::expm1(lcdf_dbl) >= bnb_cdfs_complement_min) { + out.lccdf = log1m_exp(out.lcdf); + if (Grad) { + const T ratio = exp(out.lcdf - out.lccdf); + for (int i = 0; i < 3; ++i) { + out.dlccdf[i] = -out.dlcdf[i] * ratio; + } + } + return out; + } + // the upper tail is below bnb_cdfs_complement_min + if (!have_upper && !tried_fast) { + have_upper = upper_fast(); + } + if (!have_upper) { + have_upper = bnb_cdfs_tail_sum(n, r, alpha, beta, max_steps, lc, d); + } + if (!have_upper) { + throw_domain_error(function, "number of terms of the upper sum", max_steps, + "exceeded ", ""); + } + return from_upper(); +} + +} // namespace internal +} // namespace math +} // namespace stan +#endif diff --git a/test/unit/math/rev/prob/beta_neg_binomial_cdfs_test.cpp b/test/unit/math/rev/prob/beta_neg_binomial_cdfs_test.cpp new file mode 100644 index 00000000000..06f4a0c0579 --- /dev/null +++ b/test/unit/math/rev/prob/beta_neg_binomial_cdfs_test.cpp @@ -0,0 +1,155 @@ +#include +#include +#include +#include +#include +#include +#include + +// Values and gradients of beta_neg_binomial_{cdf,lcdf,lccdf} at points that +// cover the routes of internal::beta_neg_binomial_log_cdfs(): the lower sum +// from n down, the complement of either tail, the short and the long tail +// sum, and the Thomae series with a non-integer and with an integer x. The +// points include small alpha, where the terms of the upper tail fall only +// like k^(-1 - alpha); shapes of 1e6; a cdf of 2.6e-35; and n = 20000. +// +// References: mpmath at 50 digits or more. The cdf is the finite sum of the +// pmf; the ccdf is 1 - cdf computed with extra digits, or the sum of the pmf +// over k > n; the gradients are the sums of pmf * d log pmf, checked against +// mp.diff. The parameters are written in hex so that they are exact. The +// partials are compared as x * d/dx, the gradient in log(x) that a sampler +// sees. + +namespace beta_neg_binomial_cdfs_test_internal { + +struct TestCase { + int n; + double r; + double alpha; + double beta; + double lcdf; + double lccdf; + std::vector dlcdf; // d/d(r, alpha, beta) + std::vector dlccdf; // d/d(r, alpha, beta) +}; + +// clang-format off +// NOLINTBEGIN(whitespace/line_length) +const std::vector test_cases = { + {2, 0x1.8000000000000p+1, 0x1.0000000000000p+2, 0x1.8000000000000p+1, -0.49721022174799306, -0.93706785981591312, {-0.19621186445029862, 0.15355088762562072, -0.19621186445029862}, {0.30461620945046358, -0.23838563217016256, 0.30461620945046358}}, // ordinary point + {20, 0x1.8000000000000p+1, 0x1.0000000000000p-1, 0x1.8000000000000p+1, -0.87140397613822818, -0.54191005068750497, {-0.1891559327062384, 1.6158672896523629, -0.1891559327062384}, {0.13605752717370073, -1.1622733927801345, 0.13605752717370073}}, // small alpha: the 3F2 value was wrong by 7e-3 + {20, 0x1.8000000000000p+1, 0x1.0000000000000p-2, 0x1.8000000000000p+1, -1.4760580874852567, -0.25946629095556611, {-0.23559870423820939, 3.6843570031945654, -0.23559870423820939}, {0.069793309685380625, -1.0914468742386048, 0.069793309685380625}}, // upper tail first, its complement for the cdf + {20000, 0x1.8000000000000p+1, 0x1.0000000000000p-2, 0x1.8000000000000p+1, -0.16335403587918723, -1.8924008435000104, {-0.016711519065541223, 1.3561824508346239, -0.016711519065541223}, {0.094174086515496705, -7.6424676269584797, 0.094174086515496705}}, // heavy tail, large n + {20000, 0x1.0000000000000p-1, 0x1.0000000000000p+0, 0x1.0000000000000p-1, -1.2499375029946549e-05, -11.289838162191222, {-2.499843759374463e-05, 0.00012816151869144557, -2.499843759374463e-05}, {1.999962502265497, -10.25337006183252, 1.999962502265497}}, // Thomae series, x = 0.5 + {0, 0x1.0000000000000p-1, 0x1.0000000000000p-2, 0x1.0000000000000p-1, -0.78318878541367354, -0.61054758643660245, {-0.85840734641020677, 2.2831853071795867, -0.85840734641020677}, {0.72229782263411269, -1.921162234855734, 0.72229782263411269}}, // n = 0, small shapes + {0, 0x1.e000000000000p+4, 0x1.4000000000000p+3, 0x1.2c00000000000p+8, -79.627954609150891, -2.6182957935133381e-35, {-2.1511469344641077, 1.3320590035711872, -0.09251578139693474}, {5.6323389697364848e-35, -3.4877244857620081e-35, 2.4223368126519377e-36}}, // cdf 2.6e-35: 1 - ccdf gave a NaN gradient + {0, 0x1.e848000000000p+19, 0x1.0000000000000p-2, 0x1.e848000000000p+19, -1386298.010794742, 0, {-0.69314730555993753, 17.349816535780612, -0.69314730555993753}, {0, 0, 0}}, // shapes of 1e6: the 3F2 value overflowed + {2, 0x1.e848000000000p+19, 0x1.0000000000000p+2, 0x1.e848000000000p+19, -1386223.7540788064, 0, {-0.69314443056525776, 11.866249958966133, -0.69314443056525776}, {0, 0, 0}}, // shapes of 1e6: lcdf -1.4e6, lower sum from n down + {200, 0x1.8000000000000p+1, 0x1.3880000000000p+13, 0x1.e000000000000p+4, 0, -890.84885316612292, {0, 0, 0}, {4.370163163488674, -0.019838211951244258, 2.0357763248986651}}, // tail sum: lccdf -891 + {20000, 0x1.8000000000000p+1, 0x1.9000000000000p+6, 0x1.d4c0000000000p+14, -5.8901075871069618e-20, -44.27842759632194, {-1.6515719040741391e-19, 2.8870435083041056e-20, -7.6037348214479634e-23}, {2.8039757842272959, -0.49015123503408398, 0.0012909330957030416}}, // Thomae series, integer x = 3 + {20000, 0x1.e000000000000p+4, 0x1.9000000000000p+6, 0x1.d4c0000000000p+14, -1.9902072448370877e-05, -10.824696639541852, {-1.1663144244205176e-05, 5.2405299430920191e-06, -1.5404356555183488e-08}, {0.58602078826087212, -0.26331317043269375, 0.00077399995936819138}}, // long tail sum (5.7e4 terms) + {2000, 0x1.e000000000000p+4, 0x1.0000000000000p+1, 0x1.2c00000000000p+8, -2.5815723689661887, -0.078669849505074954, {-0.11166547403255754, 1.1893396752518539, -0.010649489062478583}, {0.0091394941869035913, -0.09734399233418238, 0.00087162969774917404}}, // the 3F2 value was wrong by 5e-4 + {200, 0x1.2c00000000000p+8, 0x1.4000000000000p+3, 0x1.7700000000000p+11, -649.56328188146324, 0, {-1.9181174130266818, 3.3297281096790732, -0.089172183391557094}, {1.5174913349159807e-282, -2.6342670786201108e-282, 7.0547305755752653e-284}}, // mid shapes, alpha 10 + {20, 0x1.d4c0000000000p+14, 0x1.3880000000000p+13, 0x1.0000000000000p-1, -0.00054974363964707012, -7.5063333573606696, {-1.0207926638369315e-07, 3.0634751461336257e-07, -0.0020553852035363543}, {0.00018563418088173196, -0.0005571020634753558, 3.7377791021806801}}, // ccdf 5.5e-4: 1 - cdf would lose 3 digits +}; +// NOLINTEND +// clang-format on + +} // namespace beta_neg_binomial_cdfs_test_internal + +TEST_F(AgradRev, ProbDistributionsBetaNegBinomialCdfs_values_and_gradients) { + using beta_neg_binomial_cdfs_test_internal::test_cases; + using stan::math::var; + // tolerances times max(1, |reference|); relative for the cdf + constexpr double tol_value = 1e-11; + constexpr double tol_grad = 1e-10; + for (const auto& c : test_cases) { + const std::vector x{c.r, c.alpha, c.beta}; + for (int which = 0; which < 3; ++which) { + var r = c.r; + var alpha = c.alpha; + var beta = c.beta; + var lp; + if (which == 0) { + lp = stan::math::beta_neg_binomial_cdf(c.n, r, alpha, beta); + } else if (which == 1) { + lp = stan::math::beta_neg_binomial_lcdf(c.n, r, alpha, beta); + } else { + lp = stan::math::beta_neg_binomial_lccdf(c.n, r, alpha, beta); + } + lp.grad(); + const std::vector grads{r.adj(), alpha.adj(), beta.adj()}; + // the partials of the cdf are cdf * d lcdf + const double cdf = std::exp(c.lcdf); + const double ref = (which == 0) ? cdf : (which == 1) ? c.lcdf : c.lccdf; + const std::vector& dref = (which == 2) ? c.dlccdf : c.dlcdf; + const double scale = (which == 0) ? cdf : 1.0; + const double value_tol = (which == 0) + ? tol_value * cdf + : tol_value * std::max(1.0, std::fabs(ref)); + EXPECT_NEAR(lp.val(), ref, value_tol) + << "function " << which << ", n = " << c.n << ", r = " << c.r + << ", alpha = " << c.alpha << ", beta = " << c.beta; + for (int i = 0; i < 3; ++i) { + const double log_scale_ref = x[i] * dref[i]; + EXPECT_NEAR(x[i] * grads[i], scale * log_scale_ref, + tol_grad * scale * std::max(1.0, std::fabs(log_scale_ref))) + << "function " << which << ", gradient " << i << ", n = " << c.n + << ", r = " << c.r << ", alpha = " << c.alpha + << ", beta = " << c.beta; + } + stan::math::recover_memory(); + } + } +} + +TEST_F(AgradRev, ProbDistributionsBetaNegBinomialCdfs_vectors_and_largest_int) { + using stan::math::beta_neg_binomial_cdf; + using stan::math::beta_neg_binomial_lccdf; + using stan::math::beta_neg_binomial_lcdf; + using stan::math::var; + const std::vector ns{3, 40, 0}; + const std::vector rs{2.5, 30.0, 0.5}; + const double alpha = 3.0; + const double beta = 7.0; + + // the lcdf and lccdf of vectors are sums, the cdf is a product + double lcdf_sum = 0; + double lccdf_sum = 0; + for (size_t i = 0; i < ns.size(); ++i) { + lcdf_sum += beta_neg_binomial_lcdf(ns[i], rs[i], alpha, beta); + lccdf_sum += beta_neg_binomial_lccdf(ns[i], rs[i], alpha, beta); + } + EXPECT_NEAR(beta_neg_binomial_lcdf(ns, rs, alpha, beta), lcdf_sum, 1e-14); + EXPECT_NEAR(beta_neg_binomial_lccdf(ns, rs, alpha, beta), lccdf_sum, 1e-14); + EXPECT_NEAR(beta_neg_binomial_cdf(ns, rs, alpha, beta), std::exp(lcdf_sum), + 1e-14 * std::exp(lcdf_sum)); + + // the gradient of a vector argument holds the partials of its elements + std::vector rv(rs.begin(), rs.end()); + var lp = beta_neg_binomial_lccdf(ns, rv, alpha, beta); + lp.grad(); + std::vector adj_vector; + for (const var& v : rv) { + adj_vector.push_back(v.adj()); + } + for (size_t i = 0; i < ns.size(); ++i) { + var ri = rs[i]; + var lpi = beta_neg_binomial_lccdf(ns[i], ri, alpha, beta); + stan::math::set_zero_all_adjoints(); + lpi.grad(); + EXPECT_NEAR(adj_vector[i], ri.adj(), + 1e-14 * std::max(1.0, std::fabs(ri.adj()))); + } + stan::math::recover_memory(); + + // the largest int stands for infinity: a cdf factor of 1, an lcdf term of + // 0, and an lccdf of -inf + const std::vector n_inf{std::numeric_limits::max(), 5}; + EXPECT_DOUBLE_EQ(beta_neg_binomial_lcdf(n_inf, 3.0, 2.0, 4.0), + beta_neg_binomial_lcdf(5, 3.0, 2.0, 4.0)); + EXPECT_DOUBLE_EQ(beta_neg_binomial_cdf(n_inf, 3.0, 2.0, 4.0), + beta_neg_binomial_cdf(5, 3.0, 2.0, 4.0)); + EXPECT_EQ(beta_neg_binomial_lccdf(n_inf, 3.0, 2.0, 4.0), + -std::numeric_limits::infinity()); +} diff --git a/test/unit/math/rev/prob/large_shapes_test.cpp b/test/unit/math/rev/prob/large_shapes_test.cpp index 834e942cf11..672c8f13cc4 100644 --- a/test/unit/math/rev/prob/large_shapes_test.cpp +++ b/test/unit/math/rev/prob/large_shapes_test.cpp @@ -22,9 +22,8 @@ // // For each function the first rows are points where develop was wrong by // more than 100 times the tolerance, and the last row is a point at a small -// shape where develop was correct. The cdfs of beta_binomial and -// beta_neg_binomial are limited by the tolerance of their 3F2 series, so -// their tolerances are larger. +// shape where develop was correct. The cdfs of beta_binomial are limited by +// the tolerance of their 3F2 series, so their tolerances are larger. namespace large_shapes_test_internal { @@ -37,8 +36,7 @@ struct TestCase { // (value, gradient) tolerance factors, applied as tol * max(1, |reference|) std::pair tolerances(const std::string& tag) { - if (tag == "BBC" || tag == "BBLC" || tag == "BBLCC" || tag == "BNBC" - || tag == "BNBLCC") { + if (tag == "BBC" || tag == "BBLC" || tag == "BBLCC") { return {1e-6, 1e-4}; } return {1e-12, 1e-11};