diff --git a/stan/math/prim/prob.hpp b/stan/math/prim/prob.hpp index 4fae64c1ace..a58ee4292fa 100644 --- a/stan/math/prim/prob.hpp +++ b/stan/math/prim/prob.hpp @@ -17,6 +17,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/prob/beta_binomial_lccdf.hpp b/stan/math/prim/prob/beta_binomial_lccdf.hpp index 8df12ccc380..a69a2f6a593 100644 --- a/stan/math/prim/prob/beta_binomial_lccdf.hpp +++ b/stan/math/prim/prob/beta_binomial_lccdf.hpp @@ -3,30 +3,25 @@ #include #include -#include #include -#include -#include -#include -#include -#include -#include -#include #include #include #include #include -#include #include -#include +#include namespace stan { namespace math { /** \ingroup prob_dists - * Returns the log CCDF of the Beta-Binomial distribution with given population - * size, prior success, and prior failure parameters. Given containers of - * matching sizes, returns the log sum of probabilities. + * Returns the log CCDF of the Beta-Binomial distribution with given + * population size, 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 only where that loses no digits; see + * internal::beta_binomial_log_cdfs(). * * @tparam T_n type of success parameter * @tparam T_N type of population size parameter @@ -45,8 +40,6 @@ template inline return_type_t beta_binomial_lccdf( const T_n& n, const T_N& N, const T_size1& alpha, const T_size2& beta) { using T_partials_return = partials_return_t; - using std::exp; - using std::log; using T_N_ref = ref_type_t; using T_alpha_ref = ref_type_t; using T_beta_ref = ref_type_t; @@ -77,62 +70,25 @@ 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); - } - } - + constexpr bool any_autodiff = is_any_autodiff_v; 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 + // The gradients are technically ill-defined, but treated as zero + if (n_vec.val(i) < 0) { + continue; // a term of log(1) = 0 + } if (n_vec.val(i) >= N_vec.val(i)) { return ops_partials.build(negative_infinity()); } - - const T_partials_return n_dbl = n_vec.val(i); - const T_partials_return N_dbl = N_vec.val(i); - const T_partials_return alpha_dbl = alpha_vec.val(i); - const T_partials_return beta_dbl = beta_vec.val(i); - const T_partials_return mu = alpha_dbl + n_dbl + 1; - const T_partials_return nu = beta_dbl + N_dbl - n_dbl - 1; - const T_partials_return one = 1; - - const T_partials_return F = hypergeometric_3F2( - {one, mu, -N_dbl + n_dbl + 1}, {n_dbl + 2, 1 - nu}, one); - // lbeta(nu, mu) - lbeta(alpha, beta) without the cancellation for large - // shapes - T_partials_return C - = internal::log_beta_ratio( - alpha_dbl, beta_dbl, n_dbl + 1, N_dbl - n_dbl - 1, - internal::log_beta_ratio_denominator(alpha_dbl, beta_dbl)) - - lbeta(N_dbl - n_dbl, n_dbl + 2); - C = F * exp(C) / (N_dbl + 1); - - const T_partials_return Pi = C; - - P += log(Pi); - - // digamma(alpha + beta) - digamma(mu + nu), mu + nu = alpha + beta + N - T_partials_return digammaDiff - = is_constant_all::value - ? 0 - : -digamma_diff(alpha_dbl + beta_dbl, N_dbl); - - T_partials_return dF[6]; - if constexpr (is_any_autodiff_v) { - grad_F32( - dF, one, mu, -N_dbl + n_dbl + 1, n_dbl + 2, 1 - nu, one); - } + const auto res = internal::beta_binomial_log_cdfs( + n_vec.val(i), N_vec.val(i), T_partials_return(alpha_vec.val(i)), + T_partials_return(beta_vec.val(i))); + P += res.lccdf; if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[i] - += digamma_diff(alpha_dbl, n_dbl + 1) + digammaDiff + dF[1] / F; + partials<0>(ops_partials)[i] += res.dlccdf[0]; } if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[i] += digamma_diff(beta_dbl, N_dbl - n_dbl - 1) - + digammaDiff - dF[4] / F; + partials<1>(ops_partials)[i] += res.dlccdf[1]; } } diff --git a/stan/math/prim/prob/beta_binomial_lcdf.hpp b/stan/math/prim/prob/beta_binomial_lcdf.hpp index a0407e835f1..0e3b579cb5d 100644 --- a/stan/math/prim/prob/beta_binomial_lcdf.hpp +++ b/stan/math/prim/prob/beta_binomial_lcdf.hpp @@ -3,22 +3,13 @@ #include #include -#include #include -#include -#include -#include -#include -#include -#include -#include #include #include #include #include -#include #include -#include +#include namespace stan { namespace math { @@ -28,6 +19,10 @@ namespace math { * size, 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 only where that loses no digits; see + * internal::beta_binomial_log_cdfs(). + * * @tparam T_n type of success parameter * @tparam T_N type of population size parameter * @tparam T_size1 type of prior success parameter @@ -47,8 +42,6 @@ inline return_type_t beta_binomial_lcdf(const T_n& n, const T_size1& alpha, const T_size2& beta) { using T_partials_return = partials_return_t; - using std::exp; - using std::log; using T_N_ref = ref_type_t; using T_alpha_ref = ref_type_t; using T_beta_ref = ref_type_t; @@ -87,92 +80,22 @@ inline return_type_t beta_binomial_lcdf(const T_n& n, } } + constexpr bool any_autodiff = is_any_autodiff_v; for (size_t i = 0; i < max_size_seq_view; i++) { - 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) { - continue; - } - - const T_partials_return alpha_dbl = alpha_vec.val(i); - const T_partials_return beta_dbl = beta_vec.val(i); - const T_partials_return N_minus_n = N_dbl - n_dbl; - const T_partials_return mu = alpha_dbl + n_dbl + 1; - const T_partials_return nu = beta_dbl + N_minus_n - 1; - const T_partials_return one = 1; - - const T_partials_return F = hypergeometric_3F2({one, mu, 1 - N_minus_n}, - {n_dbl + 2, 1 - nu}, one); - // lbeta(nu, mu) - lbeta(alpha, beta) without the cancellation for large - // shapes - T_partials_return C - = internal::log_beta_ratio( - alpha_dbl, beta_dbl, n_dbl + 1, N_minus_n - 1, - internal::log_beta_ratio_denominator(alpha_dbl, beta_dbl)) - - lbeta(N_minus_n, n_dbl + 2); - C = F * exp(C) / (N_dbl + 1); - - if (C > 0.5) { - // The cdf is below 1/2, and log(1 - C) loses its digits; C can even - // round to 1 or above. Use the mirror: X <= n if and only if - // N - X > N - n - 1, and N - X is beta_binomial(N, beta, alpha). Its - // ccdf is the same series with the shapes swapped and needs no - // complement. - const T_partials_return n_m = N_minus_n - 1; - const T_partials_return mu_m = beta_dbl + N_minus_n; - const T_partials_return nu_m = alpha_dbl + n_dbl; - const T_partials_return F_m - = hypergeometric_3F2({one, mu_m, -n_dbl}, {n_m + 2, 1 - nu_m}, one); - P += internal::log_beta_ratio( - beta_dbl, alpha_dbl, N_minus_n, n_dbl, - internal::log_beta_ratio_denominator(beta_dbl, alpha_dbl)) - - lbeta(n_dbl + 1, N_minus_n + 1) + log(F_m) - log(N_dbl + 1); - if constexpr (is_any_autodiff_v) { - T_partials_return dF_m[6]; - grad_F32( - dF_m, one, mu_m, -n_dbl, n_m + 2, 1 - nu_m, one); - const T_partials_return dpsi_total - = digamma_diff(alpha_dbl + beta_dbl, N_dbl); - if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[i] - += digamma_diff(alpha_dbl, n_dbl) - dpsi_total - dF_m[4] / F_m; - } - if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[i] - += digamma_diff(beta_dbl, N_minus_n) - dpsi_total + dF_m[1] / F_m; - } - } + if (n_vec.val(i) >= N_vec.val(i)) { continue; } - - const T_partials_return Pi = 1 - C; - - P += log(Pi); - - // digamma(alpha + beta) - digamma(mu + nu), mu + nu = alpha + beta + N - T_partials_return digammaDiff - = is_constant_all::value - ? 0 - : -digamma_diff(alpha_dbl + beta_dbl, N_dbl); - - T_partials_return dF[6]; - if constexpr (is_any_autodiff_v) { - grad_F32( - dF, one, mu, 1 - N_minus_n, n_dbl + 2, 1 - nu, one); - } + const auto res = internal::beta_binomial_log_cdfs( + n_vec.val(i), N_vec.val(i), T_partials_return(alpha_vec.val(i)), + T_partials_return(beta_vec.val(i))); + P += res.lcdf; if constexpr (is_autodiff_v) { - const T_partials_return g - = -C * (digamma_diff(alpha_dbl, n_dbl + 1) + digammaDiff + dF[1] / F); - partials<0>(ops_partials)[i] += g / Pi; + partials<0>(ops_partials)[i] += res.dlcdf[0]; } if constexpr (is_autodiff_v) { - const T_partials_return g - = -C - * (digamma_diff(beta_dbl, N_minus_n - 1) + digammaDiff - dF[4] / F); - partials<1>(ops_partials)[i] += g / Pi; + partials<1>(ops_partials)[i] += res.dlcdf[1]; } } diff --git a/stan/math/prim/prob/beta_binomial_log_cdfs.hpp b/stan/math/prim/prob/beta_binomial_log_cdfs.hpp new file mode 100644 index 00000000000..9bd54e12649 --- /dev/null +++ b/stan/math/prim/prob/beta_binomial_log_cdfs.hpp @@ -0,0 +1,272 @@ +#ifndef STAN_MATH_PRIM_PROB_BETA_BINOMIAL_LOG_CDFS_HPP +#define STAN_MATH_PRIM_PROB_BETA_BINOMIAL_LOG_CDFS_HPP + +#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 binomial distribution at one + * point, and their derivatives with respect to (alpha, beta). + */ +template +struct beta_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 bb_cdfs_tol = 1e-17; + +template +inline double bb_cdfs_max_abs(const std::array& g) { + return std::max(std::fabs(value_of_rec(g[0])), std::fabs(value_of_rec(g[1]))); +} + +/** + * Returns log pmf(k) = log choose(N, k) + lbeta(k + alpha, N - k + beta) + * - lbeta(alpha, beta), without the cancellation of the two lbeta at large + * shapes. + */ +template +inline T bb_cdfs_log_pmf(int k, int N, const T& alpha, const T& beta) { + return binomial_coefficient_log(N, k) + + log_beta_ratio(alpha, beta, T(k), T(N - k), + log_beta_ratio_denominator(alpha, beta)); +} + +/** + * Writes d log pmf(k) / d(alpha, beta) to g, as differences of digamma + * functions that do not cancel at large arguments. + */ +template +inline void bb_cdfs_dlog_pmf(int k, int N, const T& alpha, const T& beta, + std::array& g) { + const T d_total = digamma_diff(T(alpha + beta), T(N)); + g[0] = digamma_diff(alpha, T(k)) - d_total; + g[1] = digamma_diff(beta, T(N - k)) - d_total; +} + +/** + * Returns whether A + B i, with i from lo to hi, has the sign sign (or is + * 0) at both ends; then it has that sign on the whole range. pmf(i + 1) / + * pmf(i) - 1 has the sign of A + B i, A = N alpha - N + 1 - beta, + * B = 2 - alpha - beta. + */ +inline bool bb_cdfs_same_sign(double a, double b, double lo, double hi, + double sign) { + return sign * (a + b * lo) >= 0.0 && sign * (a + b * hi) >= 0.0; +} + +/** + * Computes log cdf(n) as the sum of the pmf from k = n down to 0 in log + * scale, and, if Grad, its gradient as sum pmf(k) d log pmf(k) / cdf(n). + * All terms are positive. The pmf is unimodal for alpha + beta > 2 and + * U-shaped for alpha + beta < 2. Where it rises on [0, j], the j + 1 + * remaining terms are each at most pmf(j), and the sum stops when they are + * negligible. + */ +template +inline void bb_cdfs_lower_sum(int n, int N, const T& alpha, const T& beta, + T& log_cdf, std::array& dlog) { + const double a_lin = value_of_rec(N * alpha - N + 1.0 - beta); + const double b_lin = value_of_rec(2.0 - alpha - beta); + const double smaller = std::min(value_of_rec(alpha), value_of_rec(beta)); + T log_scale = bb_cdfs_log_pmf(n, N, alpha, beta); + std::array g{0, 0}; + std::array d{0, 0}; + if (Grad) { + bb_cdfs_dlog_pmf(n, N, 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;; --k) { + sum += t; + if (Grad) { + d[0] += t * g[0]; + d[1] += t * g[1]; + } + if (k == 0) { + break; + } + const double j = k - 1.0; + // pmf(j) / pmf(j + 1) + const T q = (j + 1.0) * (N - j - 1.0 + beta) / ((N - j) * (j + alpha)); + if (j == 0.0 || bb_cdfs_same_sign(a_lin, b_lin, 0.0, j - 1.0, 1.0)) { + // |d log pmf| changes by at most 2 (log1p(N / smaller) + 1 / smaller) + // over the remaining terms + const double rest = value_of_rec((j + 1.0) * t * q); + const double g_bound + = Grad ? bb_cdfs_max_abs(g) + + 2.0 * (std::log1p(N / smaller) + 1.0 / smaller) + 1.0 + : 1.0; + if (rest * g_bound <= bb_cdfs_tol * value_of_rec(sum)) { + break; + } + } + if (Grad) { + g[0] -= 1.0 / (j + alpha); + g[1] += 1.0 / (N - j - 1.0 + beta); + } + t *= q; + if (value_of_rec(t) > big || value_of_rec(sum) > big) { + t /= big; + sum /= big; + if (Grad) { + d[0] /= big; + d[1] /= big; + } + log_scale += log_big; + } + } + log_cdf = log_scale + log(sum); + if (Grad) { + dlog[0] = d[0] / sum; + dlog[1] = d[1] / sum; + } +} + +/** + * Computes log ccdf(n) as the sum of the pmf from k = n + 1 up to N in log + * scale, and, if Grad, its gradient. Where the pmf falls on [k + 1, N], the + * N - k remaining terms are each at most pmf(k + 1), and the sum stops when + * they are negligible. + */ +template +inline void bb_cdfs_upper_sum(int n, int N, const T& alpha, const T& beta, + T& log_ccdf, std::array& dlog) { + const double a_lin = value_of_rec(N * alpha - N + 1.0 - beta); + const double b_lin = value_of_rec(2.0 - alpha - beta); + const double smaller = std::min(value_of_rec(alpha), value_of_rec(beta)); + T log_scale = bb_cdfs_log_pmf(n + 1, N, alpha, beta); + std::array g{0, 0}; + std::array d{0, 0}; + if (Grad) { + bb_cdfs_dlog_pmf(n + 1, N, 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 + 1;; ++k) { + sum += t; + if (Grad) { + d[0] += t * g[0]; + d[1] += t * g[1]; + } + if (k == N) { + break; + } + // pmf(k + 1) / pmf(k) + const T rho = (N - k) * (k + alpha) / ((k + 1.0) * (N - k - 1.0 + beta)); + if (k + 1 == N || bb_cdfs_same_sign(a_lin, b_lin, k + 1.0, N - 1.0, -1.0)) { + const double rest = value_of_rec((N - k) * t * rho); + const double g_bound + = Grad ? bb_cdfs_max_abs(g) + + 2.0 * (std::log1p(N / smaller) + 1.0 / smaller) + 1.0 + : 1.0; + if (rest * g_bound <= bb_cdfs_tol * value_of_rec(sum)) { + break; + } + } + if (Grad) { + g[0] += 1.0 / (k + alpha); + g[1] -= 1.0 / (N - k - 1.0 + beta); + } + t *= rho; + if (value_of_rec(t) > big || value_of_rec(sum) > big) { + t /= big; + sum /= big; + if (Grad) { + d[0] /= big; + d[1] /= big; + } + log_scale += log_big; + } + } + log_ccdf = log_scale + log(sum); + if (Grad) { + dlog[0] = d[0] / sum; + dlog[1] = d[1] / sum; + } +} + +/** + * Returns the log cdf and the log ccdf of the beta binomial distribution at + * 0 <= n < N, and, if Grad, their derivatives with respect to (alpha, + * beta). The tail with fewer terms (n + 1 below, N - n above) is summed + * first. If it is at most 1/2, the other tail is its complement, which then + * loses no digits; else the other tail is summed too, and if that one is at + * most 1/2, the first tail is replaced by its complement. A sum near 1 + * whose first term is tiny keeps the rounding error of the log of that term + * (about 1e-13 for a log of -500); the complement of a small tail does not. + * + * @tparam Grad whether to compute the derivatives + * @tparam T type of the parameters + * @param n outcome, 0 <= n < N + * @param N number of trials + * @param alpha first shape, positive + * @param beta second shape, positive + * @return log cdf, log ccdf and, if Grad, their derivatives + */ +template +inline beta_binomial_log_cdfs_t beta_binomial_log_cdfs(int n, int N, + const T& alpha, + const T& beta) { + beta_binomial_log_cdfs_t out{0, 0, {0, 0}, {0, 0}}; + const double log_half = -LOG_TWO; + auto complement + = [](const T& l, const std::array& dl, T& c, std::array& dc) { + c = log1m_exp(l); + if (Grad) { + const T ratio = exp(l - c); + dc[0] = -dl[0] * ratio; + dc[1] = -dl[1] * ratio; + } + }; + if (n + 1 <= N - n) { + bb_cdfs_lower_sum(n, N, alpha, beta, out.lcdf, out.dlcdf); + if (value_of_rec(out.lcdf) <= log_half) { + complement(out.lcdf, out.dlcdf, out.lccdf, out.dlccdf); + } else { + bb_cdfs_upper_sum(n, N, alpha, beta, out.lccdf, out.dlccdf); + if (value_of_rec(out.lccdf) <= log_half) { + complement(out.lccdf, out.dlccdf, out.lcdf, out.dlcdf); + } + } + } else { + bb_cdfs_upper_sum(n, N, alpha, beta, out.lccdf, out.dlccdf); + if (value_of_rec(out.lccdf) <= log_half) { + complement(out.lccdf, out.dlccdf, out.lcdf, out.dlcdf); + } else { + bb_cdfs_lower_sum(n, N, alpha, beta, out.lcdf, out.dlcdf); + if (value_of_rec(out.lcdf) <= log_half) { + complement(out.lcdf, out.dlcdf, out.lccdf, out.dlccdf); + } + } + } + return out; +} + +} // namespace internal +} // namespace math +} // namespace stan +#endif diff --git a/test/unit/math/rev/prob/beta_binomial_log_cdfs_test.cpp b/test/unit/math/rev/prob/beta_binomial_log_cdfs_test.cpp new file mode 100644 index 00000000000..ae7886ac534 --- /dev/null +++ b/test/unit/math/rev/prob/beta_binomial_log_cdfs_test.cpp @@ -0,0 +1,157 @@ +#include +#include +#include +#include +#include +#include +#include + +// Values and gradients of beta_binomial_{cdf,lcdf,lccdf} at points that +// cover the routes of internal::beta_binomial_log_cdfs(): the lower or the +// upper tail summed first, the other one as its complement or summed too, +// then the first one as the complement of a small second one, the unimodal +// and the U-shaped pmf, an early stop, tails near 1e-276, and 200 000 +// trials. +// +// References: mpmath at 50 digits; both tails are direct sums of the pmf, +// with the derivatives as sums of pmf * d log pmf; the larger tail is then +// taken as 1 minus the smaller one, so that its log keeps its digits. The +// shapes 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_binomial_log_cdfs_test_internal { + +struct TestCase { + int n; + int N; + double alpha; + double beta; + double lcdf; + double lccdf; + std::vector dlcdf; // d/d(alpha, beta) + std::vector dlccdf; // d/d(alpha, beta) +}; + +// clang-format off +// NOLINTBEGIN(whitespace/line_length) +const std::vector test_cases = { + {2, 10, 0x1.8000000000000p+1, 0x1.0000000000000p+2, -1.4198170531585341, -0.27675300191959068, {-0.51879717504717504, 0.30789835164835166}, {0.16541359204402684, -0.098170488931358491}}, // lower tail first, upper as its complement + {8, 10, 0x1.8000000000000p+1, 0x1.0000000000000p+2, -0.036367644170874847, -3.3322045101752038, {-0.022532560032560031, 0.027348937071159294}, {0.60837912087912083, -0.73842130092130087}}, // upper tail first, lower as its complement + {8, 10, 0x1.0000000000000p-1, 0x1.e000000000000p+4, -7.8114337516020798e-09, -18.667677314898757, {-3.0283500032506201e-08, 1.994662475956812e-09}, {3.8768170962233981, -0.25535164626559514}}, // small upper tail + {1, 10, 0x1.e000000000000p+4, 0x1.0000000000000p-1, -18.667677314898757, -7.8114337516020798e-09, {-0.25535164626559514, 3.8768170962233981}, {1.994662475956812e-09, -3.0283500032506201e-08}}, // small lower tail + {58, 117, 0x1.0000000000000p+1, 0x1.0000000000000p+1, -0.69314718055994529, -0.69314718055994529, {-0.43689703406552266, 0.43689703406552266}, {0.43689703406552266, -0.43689703406552266}}, // lower tail above 1/2: both tails summed + {58, 117, 0x1.4000000000000p+3, 0x1.3880000000000p+13, -4.4746105258856488e-112, -256.39111110071553, {-8.7861915440757089e-112, 2.6123844465240352e-114}, {1.9635656540938118, -0.0058382387280666677}}, // both tails summed, lcdf as the complement of the upper tail + {59, 117, 0x1.3880000000000p+13, 0x1.4000000000000p+3, -251.39436575466911, -6.6193321069869591e-110, {-0.005738828270701637, 1.9488631096182969}, {3.7987210228739794e-112, -1.2900172153618836e-109}}, // both tails summed, lccdf as the complement of the lower tail + {60, 117, 0x1.3333333333333p-2, 0x1.999999999999ap-3, -0.91507251589898242, -0.51163859300372261, {-2.1340239197768334, 3.1499047202585775}, {1.4255752818494947, -2.1042061749013872}}, // U-shaped pmf: no early stop + {5, 117, 0x1.f400000000000p+9, 0x1.7700000000000p+11, -19.695873917018684, -2.7937699528797345e-09, {-0.023991525922721486, 0.0078730012808713602}, {6.7026804340263269e-11, -2.1995354448207051e-11}}, // large shapes: early stop of the lower sum + {3, 10, 0x1.2a4da7abab141p+8, 0x1.258191c361692p+7, -3.933350026668649, -0.019772163107297, {-0.012878429074005715, 0.025788239904711231}, {0.00025716841006185448, -0.00051496348013239945}}, // the 3F2 series once gave NaN here + {0, 117, 0x1.3880000000000p+13, 0x1.0000000000000p+0, -634.8223474979892, -1.9959860090434558e-276, {-0.011631506014216147, 5.3320255103901211}, {2.3216323268480239e-278, -1.0642648318601472e-275}}, // cdf about 2e-276 + {116, 117, 0x1.0000000000000p+0, 0x1.3880000000000p+13, -1.9959860090434558e-276, -634.8223474979892, {-1.0642648318601472e-275, 2.3216323268480239e-278}, {5.3320255103901211, -0.011631506014216147}}, // ccdf about 2e-276 + {0, 1, 0x1.0000000000000p+1, 0x1.8000000000000p+1, -0.51082562376599072, -0.91629073187415511, {-0.20000000000000001, 0.13333333333333333}, {0.29999999999999999, -0.20000000000000001}}, // N = 1 + {1000, 10000, 0x1.9000000000000p+6, 0x1.e000000000000p+4, -162.91049486808714, -1.773662707038273e-71, {-1.9664909412961633, 1.3654243090365596}, {3.4878916463055947e-71, -2.421802176221648e-71}}, // lccdf gradient about 1e-71 + {50000, 200000, 0x1.0000000000000p-1, 0x1.0000000000000p-1, -1.0986053971070815, -0.40546855390648867, {-1.9383781898621288, 1.2922610759880457}, {0.96919911383880364, -0.64613721731210594}}, // 200 000 trials, U-shaped +}; +// NOLINTEND +// clang-format on + +} // namespace beta_binomial_log_cdfs_test_internal + +TEST_F(AgradRev, ProbDistributionsBetaBinomialCdfs_values_and_gradients) { + using beta_binomial_log_cdfs_test_internal::test_cases; + using stan::math::var; + // the values relative to the reference, so that an lcdf or lccdf near 0 + // keeps its digits; the gradients relative to max(1, |reference|) + constexpr double tol_value = 1e-11; + constexpr double tol_grad = 1e-10; + for (const auto& c : test_cases) { + const std::vector x{c.alpha, c.beta}; + for (int which = 0; which < 3; ++which) { + var alpha = c.alpha; + var beta = c.beta; + var lp; + if (which == 0) { + lp = stan::math::beta_binomial_cdf(c.n, c.N, alpha, beta); + } else if (which == 1) { + lp = stan::math::beta_binomial_lcdf(c.n, c.N, alpha, beta); + } else { + lp = stan::math::beta_binomial_lccdf(c.n, c.N, alpha, beta); + } + lp.grad(); + const std::vector grads{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; + EXPECT_NEAR(lp.val(), ref, tol_value * std::fabs(ref)) + << "function " << which << ", n = " << c.n << ", N = " << c.N + << ", alpha = " << c.alpha << ", beta = " << c.beta; + for (int i = 0; i < 2; ++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 + << ", N = " << c.N << ", alpha = " << c.alpha + << ", beta = " << c.beta; + } + stan::math::recover_memory(); + } + } +} + +TEST_F(AgradRev, ProbDistributionsBetaBinomialCdfs_vectors_and_edges) { + using stan::math::beta_binomial_cdf; + using stan::math::beta_binomial_lccdf; + using stan::math::beta_binomial_lcdf; + using stan::math::var; + const std::vector ns{3, 40, 0}; + const std::vector Ns{10, 117, 5}; + const std::vector alphas{2.5, 30.0, 0.5}; + 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_binomial_lcdf(ns[i], Ns[i], alphas[i], beta); + lccdf_sum += beta_binomial_lccdf(ns[i], Ns[i], alphas[i], beta); + } + EXPECT_NEAR(beta_binomial_lcdf(ns, Ns, alphas, beta), lcdf_sum, 1e-14); + EXPECT_NEAR(beta_binomial_lccdf(ns, Ns, alphas, beta), lccdf_sum, 1e-14); + EXPECT_NEAR(beta_binomial_cdf(ns, Ns, alphas, 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 av(alphas.begin(), alphas.end()); + var lp = beta_binomial_lccdf(ns, Ns, av, beta); + lp.grad(); + std::vector adj_vector; + for (const var& v : av) { + adj_vector.push_back(v.adj()); + } + for (size_t i = 0; i < ns.size(); ++i) { + var ai = alphas[i]; + var lpi = beta_binomial_lccdf(ns[i], Ns[i], ai, beta); + stan::math::set_zero_all_adjoints(); + lpi.grad(); + EXPECT_NEAR(adj_vector[i], ai.adj(), + 1e-14 * std::max(1.0, std::fabs(ai.adj()))); + } + stan::math::recover_memory(); + + // n < 0 is a term of log(1) = 0 in the lccdf and makes the lcdf -inf; + // n >= N is a term of 0 in the lcdf and makes the lccdf -inf + const double inf = std::numeric_limits::infinity(); + EXPECT_DOUBLE_EQ(beta_binomial_lccdf(std::vector{-1, 7}, + std::vector{20, 20}, 2.5, 4.0), + beta_binomial_lccdf(7, 20, 2.5, 4.0)); + EXPECT_EQ(beta_binomial_lcdf(std::vector{-1, 7}, + std::vector{20, 20}, 2.5, 4.0), + -inf); + EXPECT_DOUBLE_EQ(beta_binomial_lcdf(std::vector{20, 7}, + std::vector{20, 20}, 2.5, 4.0), + beta_binomial_lcdf(7, 20, 2.5, 4.0)); + EXPECT_EQ(beta_binomial_lccdf(std::vector{20, 7}, + std::vector{20, 20}, 2.5, 4.0), + -inf); +}