From 2764cf584d1cdf5bbf4bc6146a303d3dd57fa4d6 Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Mon, 26 Oct 2026 15:02:32 -0400 Subject: [PATCH 1/2] Compute the incomplete beta shape gradients by continued fraction Port the Boik and Robinson-Cox continued fraction (as in TensorFlow Probability) to prim grad_reg_inc_beta, route inc_beta_dda, inc_beta_ddb and fwd inc_beta through it, and test against fixed mpmath references. --- stan/math/fwd/fun/inc_beta.hpp | 87 +-- stan/math/prim/fun.hpp | 3 + stan/math/prim/fun/grad_reg_inc_beta.hpp | 500 +++++++++++++++++- stan/math/prim/fun/inc_beta_dda.hpp | 135 +++-- stan/math/prim/fun/inc_beta_ddb.hpp | 88 +-- stan/math/prim/fun/inc_beta_ddz.hpp | 53 +- .../math/prim/fun/grad_reg_inc_beta_test.cpp | 100 +++- test/unit/math/prim/fun/inc_beta_dda_test.cpp | 130 +++-- test/unit/math/prim/fun/inc_beta_ddb_test.cpp | 98 ++-- test/unit/math/prim/fun/inc_beta_ddz_test.cpp | 48 ++ 10 files changed, 927 insertions(+), 315 deletions(-) diff --git a/stan/math/fwd/fun/inc_beta.hpp b/stan/math/fwd/fun/inc_beta.hpp index d0b89e35e06..d0c75a711aa 100644 --- a/stan/math/fwd/fun/inc_beta.hpp +++ b/stan/math/fwd/fun/inc_beta.hpp @@ -3,54 +3,59 @@ #include #include -#include #include -#include -#include -#include #include #include -#include +#include +#include namespace stan { namespace math { -template -inline fvar inc_beta(const fvar& a, const fvar& b, const fvar& x) { - T d_a; - T d_b; - const T beta_ab = beta(a.val_, b.val_); - grad_reg_inc_beta(d_a, d_b, a.val_, b.val_, x.val_, digamma(a.val_), - digamma(b.val_), digamma(a.val_ + b.val_), beta_ab); - T d_x = pow((1 - x.val_), b.val_ - 1) * pow(x.val_, a.val_ - 1) / beta_ab; - return fvar(inc_beta(a.val_, b.val_, x.val_), - a.d_ * d_a + b.d_ * d_b + x.d_ * d_x); -} - -template -inline fvar inc_beta(double a, const fvar& b, const fvar& x) { - return inc_beta(fvar(a), b, x); -} -template -inline fvar inc_beta(const fvar& a, double b, const fvar& x) { - return inc_beta(a, fvar(b), x); -} -template -inline fvar inc_beta(const fvar& a, const fvar& b, double x) { - return inc_beta(a, b, fvar(x)); -} - -template -inline fvar inc_beta(double a, double b, const fvar& x) { - return inc_beta(fvar(a), fvar(b), x); -} -template -inline fvar inc_beta(const fvar& a, double b, double x) { - return inc_beta(a, fvar(b), fvar(x)); -} -template -inline fvar inc_beta(double a, const fvar& b, double x) { - return inc_beta(fvar(a), b, fvar(x)); +/** + * The regularized incomplete beta function I_x(a, b) for fvar + * arguments. Only the partials of the fvar arguments are + * computed. Each argument is an fvar or arithmetic. Mixing with + * var is not supported. + * + * @tparam T1 type of a + * @tparam T2 type of b + * @tparam T3 type of x + * @param a first shape, positive + * @param b second shape, positive + * @param x upper limit, in [0, 1] + * @return I_x(a, b) + */ +template * = nullptr, + require_all_not_var_t* = nullptr, + require_any_fvar_t* = nullptr> +inline fvar> inc_beta(const T1& a, const T2& b, + const T3& x) { + using T_return = partials_return_t; + const auto a_val = value_of(a); + const auto b_val = value_of(b); + const auto x_val = value_of(x); + T_return d(0); + if constexpr (is_autodiff_v && is_autodiff_v) { + const auto [d_a, d_b] = grad_reg_inc_beta( + a_val, b_val, x_val, digamma(a_val), digamma(b_val), + digamma(a_val + b_val), 0.0); + d += a.d_ * d_a + b.d_ * d_b; + } else if constexpr (is_autodiff_v) { + d += a.d_ + * grad_reg_inc_beta(a_val, b_val, x_val, digamma(a_val), + 0.0, digamma(a_val + b_val), 0.0); + } else if constexpr (is_autodiff_v) { + d += b.d_ + * grad_reg_inc_beta(a_val, b_val, x_val, 0.0, + digamma(b_val), + digamma(a_val + b_val), 0.0); + } + if constexpr (is_autodiff_v) { + d += x.d_ * inc_beta_ddz(a_val, b_val, x_val); + } + return fvar(inc_beta(a_val, b_val, x_val), d); } } // namespace math diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp index 8ef4204c81b..00cab158ab8 100644 --- a/stan/math/prim/fun.hpp +++ b/stan/math/prim/fun.hpp @@ -134,6 +134,9 @@ #include #include #include +#include +#include +#include #include #include #include diff --git a/stan/math/prim/fun/grad_reg_inc_beta.hpp b/stan/math/prim/fun/grad_reg_inc_beta.hpp index f4fa1b58178..92fab1fc056 100644 --- a/stan/math/prim/fun/grad_reg_inc_beta.hpp +++ b/stan/math/prim/fun/grad_reg_inc_beta.hpp @@ -2,26 +2,484 @@ #define STAN_MATH_PRIM_FUN_GRAD_REG_INC_BETA_HPP #include -#include -#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include +#include +#include +#include namespace stan { namespace math { +/** + * Selects which partial derivatives of the regularized incomplete beta + * function grad_reg_inc_beta computes: G1 for the + * derivative with respect to a, G2 for the derivative with + * respect to b, BOTH for both. + */ +enum class GRAD_RIB : uint8_t { G1 = 0, G2 = 1, BOTH = 2 }; + +namespace internal { + +/** + * Partial derivatives of the regularized incomplete beta function I_x(a, b) + * with respect to a and b, computed from its power series. + * + * The series is + * + * I_x(a, b) = N S, N = x^a / B(a, b), S = sum_{n >= 0} t_n, + * t_n = P_n / (a + n), P_0 = 1, P_n = P_{n-1} (n - b) x / n, + * + * where S is 2F1(a, 1 - b; a + 1; x) / a. The derivatives of each term are + * + * dt_n/da = -t_n / (a + n), dt_n/db = (dP_n/db) / (a + n), + * + * and dP_n/db follows from the recurrence for P_n by the product rule. With + * log N = a log x - log B(a, b) and psi the digamma function, + * + * dI/da = N (dS/da + S (log x - psi(a) + psi(a + b))), + * dI/db = N (dS/db + S (psi(a + b) - psi(b))). + * + * Only the sums needed for the requested partials are accumulated. + * + * The loop stops once a term is below eps / a in magnitude, or after 600 + * terms. For a positive integer b, P_n is exactly 0 for n >= b but dP_n/db is + * not, so the sum for dS/db keeps changing after S has converged. The stop + * test therefore also requires dt_n/db to be small whenever that sum + * matters. That is the case when d/db is requested, and also when b is an + * autodiff type, because autodiff then differentiates the d/da result in b + * and needs the same terms. + * + * @tparam grads_requested which partials to compute + * @tparam flip_sign true when the caller evaluates its I_z(a, b) through the + * symmetry relation I_z(a, b) = 1 - I_{1-z}(b, a) and so passes + * (b, a, 1 - z). The results are then negated, and for BOTH + * they are returned in the caller's order of a and b. + * @tparam T1 type of a + * @tparam T2 type of b + * @tparam T3 type of x + * @tparam T4 type of digamma_a + * @tparam T5 type of digamma_b + * @tparam T6 type of digamma_ab + * @param a first shape + * @param b second shape + * @param x upper limit, in (0, 1) + * @param digamma_a digamma(a) + * @param digamma_b digamma(b) + * @param digamma_ab digamma(a + b) + * @return d/da I_x(a, b) for G1, d/db for G2, and + * the pair (d/da, d/db) for BOTH. With flip_sign + * these are negated and the pair is returned as (-d/db, -d/da), which are + * the caller's (d/da, d/db). + */ +template +inline auto inc_beta_power_series_kernel(const T1& a, const T2& b, const T3& x, + const T4& digamma_a, + const T5& digamma_b, + const T6& digamma_ab) { + using ret_t = return_type_t; + constexpr int max_iterations = 600; + const double tolerance + = std::numeric_limits::epsilon() / value_of_rec(a); + + // series for 2F1(a, 1 - b; a + 1; x) / a and its a and b derivatives + ret_t product = 1.0; + ret_t series_sum = inv(a); + ret_t da = -series_sum / a; + ret_t db = 0.0; + ret_t product_grad_b = 0.0; + for (int n = 1; n <= max_iterations; ++n) { + const auto x_div_n = x / n; + const auto factor = (n - b) * x_div_n; + const auto apn = a + n; + const auto new_product = product * factor; + const auto term = new_product / apn; + bool db_converged = true; + if constexpr (grads_requested != GRAD_RIB::G1 || is_autodiff_v) { + product_grad_b = factor * product_grad_b - product * x_div_n; + const auto db_increment = product_grad_b / apn; + if constexpr (grads_requested != GRAD_RIB::G1) { + db += db_increment; + } + db_converged = fabs(value_of_rec(db_increment)) <= tolerance; + } + if constexpr (grads_requested != GRAD_RIB::G2) { + da -= term / apn; + } + product = new_product; + series_sum += term; + if (fabs(value_of_rec(term)) <= tolerance && db_converged) { + break; + } + } + + // the symmetry relation negates the partials + constexpr double sign = flip_sign ? -1.0 : 1.0; + const auto log_x = log(x); + const auto normalize = sign * exp(a * log_x - lbeta(a, b)); + if constexpr (grads_requested == GRAD_RIB::G1) { + return normalize * (da + series_sum * (digamma_ab - digamma_a + log_x)); + } else if constexpr (grads_requested == GRAD_RIB::G2) { + return normalize * (db + series_sum * (digamma_ab - digamma_b)); + } else { + const auto grad_a + = normalize * (da + series_sum * (digamma_ab - digamma_a + log_x)); + const auto grad_b + = normalize * (db + series_sum * (digamma_ab - digamma_b)); + if constexpr (flip_sign) { + // the caller swapped a and b, so its (d/da, d/db) are (grad_b, grad_a) + return std::pair{grad_b, grad_a}; + } else { + return std::pair{grad_a, grad_b}; + } + } +} + +/** + * Partial derivatives of the regularized incomplete beta function I_x(a, b) + * with respect to a and b, computed from its continued fraction. + * + * The continued fraction https://dlmf.nist.gov/8.17.E22 is + * + * I_x(a, b) = Q h, Q = x^a (1 - x)^b / (a B(a, b)), + * h = 1 / (1 + d_1 / (1 + d_2 / (1 + ...))), + * + * with the partial numerators https://dlmf.nist.gov/8.17.E23 + * + * d_1 = -(a + b) x / (a + 1), + * d_{2m} = m (b - m) x / ((a + 2m - 1)(a + 2m)), + * d_{2m+1} = -(a + m)(a + b + m) x / ((a + 2m)(a + 2m + 1)). + * + * Each d_n is num / den with den = (a + n - 1)(a + n), so its a derivative is + * (num_a - d_n den_a) / den, with den_a = 2 (a + n) - 1. The even numerators + * do not depend on a, and the odd ones have num_a = -(2 (a + m) + b) x. Every + * numerator is linear in b, so the b derivative is num_b / den. + * + * h is evaluated with the modified Lentz method. With h_n = A_n / B_n the + * n-th convergent, it keeps the ratios C = A_n / A_{n-1} and + * D = B_{n-1} / B_n, and step n updates them as + * + * C <- 1 + d_n / C, D <- 1 / (1 + d_n D), h <- h C D. + * + * The first step holds C at 1, so D = 1 / (1 + d_1) and h = D. When + * 1 + d_n / C or 1 + d_n D is smaller in magnitude than sqrt(DBL_MIN), it is + * replaced by that value, so that nothing divides by 0. + * + * The derivatives are carried through the same steps, following Thompson and + * Barnett (1986), appendix. They are kept as logarithmic derivatives + * c = (dC/da) / C, e = (dD/da) / D, and r = (dh/da) / h, and likewise for b. + * Differentiating the updates gives + * + * c <- (dd_n/da - d_n c) / (C C_new), + * e <- -(dd_n/da + d_n e) D D_new, + * r <- r + c + e, + * + * where C, D, c, and e on the right are the values before the step and + * C_new, D_new are the values after it. At the end, with psi the digamma + * function, + * + * dI/da = I (r_a + log x - 1 / a + psi(a + b) - psi(a)), + * dI/db = I (r_b + log(1 - x) + psi(a + b) - psi(b)). + * + * Only the requested derivatives are carried. + * + * The test for stopping runs after each odd n. The loop stops once C D is + * within 3 eps of 1 and each carried derivative dh = h r has changed by at + * most a relative 3 eps since the previous test, or after n = 601. For + * integer b the partial numerator d_{2b} is 0, so from that step on h and its + * a derivative no longer change. The b derivative of d_{2b} is b x / den, + * which is not 0, so dh/db is still changing. The b derivative is therefore + * tested whenever it matters, which is when d/db is requested and also when b + * is an autodiff type, as in the power series. + * + * @tparam grads_requested which partials to compute + * @tparam flip_sign true when the caller evaluates its I_z(a, b) through the + * symmetry relation I_z(a, b) = 1 - I_{1-z}(b, a) and so passes + * (b, a, 1 - z). The results are then negated, and for BOTH + * they are returned in the caller's order of a and b. + * @tparam T1 type of a + * @tparam T2 type of b + * @tparam T3 type of x + * @tparam T4 type of digamma_a + * @tparam T5 type of digamma_b + * @tparam T6 type of digamma_ab + * @param a first shape + * @param b second shape + * @param x upper limit, in (0, 1) + * @param digamma_a digamma(a) + * @param digamma_b digamma(b) + * @param digamma_ab digamma(a + b) + * @return d/da I_x(a, b) for G1, d/db for G2, and + * the pair (d/da, d/db) for BOTH. With flip_sign + * these are negated and the pair is returned as (-d/db, -d/da), which are + * the caller's (d/da, d/db). + */ +template +inline auto inc_beta_continued_fraction_kernel(const T1& a, const T2& b, + const T3& x, const T4& digamma_a, + const T5& digamma_b, + const T6& digamma_ab) { + using ret_t = return_type_t; + constexpr int max_iterations = 300; + const double tolerance = 3.0 * std::numeric_limits::epsilon(); + const double small = std::sqrt(std::numeric_limits::min()); + + // first step, with C held at 1 + const auto ap1 = a + 1.0; + const auto d_1 = -(a + b) * x / ap1; + ret_t C = 1.0; + ret_t D = 1.0 + d_1; + if (fabs(value_of_rec(D)) < small) { + D = small; + } + D = inv(D); + ret_t h = D; + ret_t c_a = 0.0; + ret_t e_a = 0.0; + ret_t r_a = 0.0; + ret_t c_b = 0.0; + ret_t e_b = 0.0; + ret_t r_b = 0.0; + if constexpr (grads_requested != GRAD_RIB::G2) { + // d d_1 / da = -(x + d_1) / (a + 1) + e_a = (x + d_1) / ap1 * D; + r_a = e_a; + } + if constexpr (grads_requested != GRAD_RIB::G1 || is_autodiff_v) { + // d d_1 / db = -x / (a + 1) + e_b = x / ap1 * D; + r_b = e_b; + } + + // the stop test compares dh = h r after each odd step + double dh_a_prev = value_of_rec(h) * value_of_rec(r_a); + double dh_b_prev = value_of_rec(h) * value_of_rec(r_b); + for (int n = 2; n <= 2 * max_iterations + 1; ++n) { + // the partial numerator d_n and its b derivative dd_b + const int m = n / 2; + const bool odd = n % 2 == 1; + const auto a_plus_n = a + n; + const auto denominator = a_plus_n * (a_plus_n - 1.0); + ret_t dd_b; + ret_t d; + if (odd) { + dd_b = -(a + m) * x / denominator; + d = dd_b * (a + b + m); + } else { + dd_b = m * x / denominator; + d = dd_b * (b - m); + } + // one step of the modified Lentz method + auto C_new = 1.0 + d / C; + // C C_new = C + d, which does not wait on the division for C_new + ret_t C_C_new = C + d; + if (fabs(value_of_rec(C_new)) < small) { + C_new = small; + C_C_new = C * small; + } + auto D_new = 1.0 + d * D; + if (fabs(value_of_rec(D_new)) < small) { + D_new = small; + } + D_new = inv(D_new); + // c_new = (dd - d c) / (C C_new) and e_new = -(dd + d e) D D_new + const auto D_D_new = D * D_new; + if constexpr (grads_requested != GRAD_RIB::G2) { + auto dd_a = -d * (2.0 * a_plus_n - 1.0) / denominator; + if (odd) { + dd_a -= x * (2.0 * (a + m) + b) / denominator; + } + c_a = (dd_a - d * c_a) / C_C_new; + e_a = -(dd_a + d * e_a) * D_D_new; + r_a += c_a + e_a; + } + if constexpr (grads_requested != GRAD_RIB::G1 || is_autodiff_v) { + c_b = (dd_b - d * c_b) / C_C_new; + e_b = -(dd_b + d * e_b) * D_D_new; + r_b += c_b + e_b; + } + C = C_new; + D = D_new; + const auto delta = C_new * D_new; + h *= delta; + if (odd) { + const double h_d = value_of_rec(h); + const double dh_a_now = h_d * value_of_rec(r_a); + const double dh_b_now = h_d * value_of_rec(r_b); + if (fabs(value_of_rec(delta) - 1.0) < tolerance + && (grads_requested == GRAD_RIB::G2 + || fabs(dh_a_now - dh_a_prev) <= tolerance * fabs(dh_a_now)) + && ((grads_requested == GRAD_RIB::G1 && !is_autodiff_v) + || fabs(dh_b_now - dh_b_prev) <= tolerance * fabs(dh_b_now))) { + break; + } + dh_a_prev = dh_a_now; + dh_b_prev = dh_b_now; + } + } + + // the symmetry relation negates the partials + constexpr double sign = flip_sign ? -1.0 : 1.0; + const auto log_x = log(x); + const auto log1m_x = log1p(-x); + // h times the normalization is I_x(a, b), shared by both partials + const auto prefactor + = sign * h * exp(a * log_x + b * log1m_x - log(a) - lbeta(a, b)); + if constexpr (grads_requested == GRAD_RIB::G1) { + return ret_t(prefactor * (r_a + log_x - inv(a) + digamma_ab - digamma_a)); + } else if constexpr (grads_requested == GRAD_RIB::G2) { + return ret_t(prefactor * (r_b + log1m_x + digamma_ab - digamma_b)); + } else { + const ret_t grad_a + = prefactor * (r_a + log_x - inv(a) + digamma_ab - digamma_a); + const ret_t grad_b = prefactor * (r_b + log1m_x + digamma_ab - digamma_b); + if constexpr (flip_sign) { + // the caller swapped a and b, so its (d/da, d/db) are (grad_b, grad_a) + return std::pair{grad_b, grad_a}; + } else { + return std::pair{grad_a, grad_b}; + } + } +} + +} // namespace internal + +/** + * Computes the requested partial derivatives of the regularized incomplete + * beta function I_z(a, b) = inc_beta(a, b, z) with respect to + * a and b. + * + * The method is that of Boik and Robinson-Cox (1998), "Derivatives of the + * Incomplete Beta Function", Journal of Statistical Software 3(1), as + * implemented in TensorFlow Probability. I_z(a, b) is expanded either as a + * power series or as a continued fraction, and the expansion is + * differentiated term by term. Only the requested derivatives are carried + * through it. The normalizing factor is evaluated in log space, so it stays + * finite where B(a, b) underflows. + * + * Either expansion can instead be applied to the right side of the symmetry + * relation I_z(a, b) = 1 - I_{1-z}(b, a), which converges faster for some + * arguments. The expansion is then evaluated for I_{1-z}(b, a), and the + * partials of I_z(a, b) are its partials negated, with the roles of a and b + * swapped. + * + * With mean = a / (a + b), the branch is chosen as follows. + * + * - The power series for I_z(a, b) if z < mean, b z <= 1, and z <= 0.5. + * - The power series for I_{1-z}(b, a) if z >= mean, a (1 - z) <= 1, and + * z >= 0.5. + * - Otherwise, the continued fraction for I_z(a, b) if + * z < (a + 1) / (a + b + 2), and for I_{1-z}(b, a) if not. + * + * The terms of the power series shrink by a factor of about z per term + * (1 - z for the series of I_{1-z}(b, a)), so the series is limited to + * z <= 0.5. The condition b z <= 1 does not bound z by itself when b < 1. + * Without the limit, the series would take most of its 600 terms at z = 0.95, + * where the continued fraction needs tens of steps. + * + * The continued fraction converges rapidly for z < (a + 1) / (a + b + 2) + * (Numerical Recipes, 3rd ed., section 6.4). When z is at or above that + * point, 1 - z is at or below (b + 1) / (a + b + 2), the matching point for + * the fraction of I_{1-z}(b, a). The mode (a - 1) / (a + b - 2) is not used + * as the threshold because it exceeds 1 for b < 1. The fraction of I_z(a, b) + * would then be used for z near 1, where it needs hundreds of steps. + * + * The partials are 0 at z = 0 and z = 1. + * + * @tparam grads_requested G1 for the derivative with respect to + * a, G2 for the derivative with respect to b, + * BOTH for both + * @tparam T1 type of a + * @tparam T2 type of b + * @tparam T3 type of z + * @tparam T4 type of digammaA + * @tparam T5 type of digammaB + * @tparam T6 type of digammaSum + * @tparam T7 type of betaAB + * @param a a + * @param b b + * @param z z + * @param digammaA the value of digamma(a), not used for + * G2 + * @param digammaB the value of digamma(b), not used for + * G1 + * @param digammaSum the value of digamma(a + b) + * @param betaAB accepted for interface compatibility and not used. The + * normalization is evaluated in log space. + * @return d/da I_z(a, b) for G1, d/db for G2, and + * the pair (d/da, d/db) for BOTH + */ +template +inline auto grad_reg_inc_beta(const T1& a, const T2& b, const T3& z, + const T4& digammaA, const T5& digammaB, + const T6& digammaSum, const T7& /* betaAB */) { + // T7 is left out so that every branch returns the same type + using ret_t = return_type_t; + // the partials requested after the symmetry relation swaps a and b + constexpr GRAD_RIB swapped_grads + = grads_requested == GRAD_RIB::G1 + ? GRAD_RIB::G2 + : grads_requested == GRAD_RIB::G2 ? GRAD_RIB::G1 : GRAD_RIB::BOTH; + + const double z_d = value_of_rec(z); + if (z_d == 0.0 || z_d == 1.0) { + if constexpr (grads_requested == GRAD_RIB::BOTH) { + return std::pair{0.0, 0.0}; + } else { + return ret_t{0.0}; + } + } + + const double a_d = value_of_rec(a); + const double b_d = value_of_rec(b); + const double mean = a_d / (a_d + b_d); + const bool use_power_series + = (z_d < mean && b_d * z_d <= 1.0 && z_d <= 0.5) + || (z_d >= mean && a_d * (1.0 - z_d) <= 1.0 && z_d >= 0.5); + if (use_power_series) { + if (z_d >= mean) { + return internal::inc_beta_power_series_kernel( + b, a, 1.0 - z, digammaB, digammaA, digammaSum); + } else { + return internal::inc_beta_power_series_kernel( + a, b, z, digammaA, digammaB, digammaSum); + } + } else { + if (z_d >= (a_d + 1.0) / (a_d + b_d + 2.0)) { + return internal::inc_beta_continued_fraction_kernel( + b, a, 1.0 - z, digammaB, digammaA, digammaSum); + } else { + return internal::inc_beta_continued_fraction_kernel( + a, b, z, digammaA, digammaB, digammaSum); + } + } +} + /** * Computes the gradients of the regularized incomplete beta * function. Specifically, this function computes gradients of - * ibeta(a, b, z), with respect to the arguments + * inc_beta(a, b, z) with respect to the arguments * a and b. * - * Uses the equivalence to a hypergeometric function. See - * http://dlmf.nist.gov/8.17#ii + * Writes both partials. To compute only one, use the overload templated on + * GRAD_RIB. * * @tparam T type of arguments - * @param[out] g1 partial derivative of ibeta(a, b, z) + * @param[out] g1 partial derivative of inc_beta(a, b, z) * with respect to a - * @param[out] g2 partial derivative of ibeta(a, b, + * @param[out] g2 partial derivative of inc_beta(a, b, * z) with respect to b * @param[in] a a * @param[in] b b @@ -29,32 +487,18 @@ namespace math { * @param[in] digammaA the value of digamma(a) * @param[in] digammaB the value of digamma(b) * @param[in] digammaSum the value of digamma(a + b) - * @param[in] betaAB the value of beta(a, b) + * @param[in] betaAB the value of beta(a, b), accepted for + * interface compatibility and not used. The normalization is evaluated in + * log space. */ template inline void grad_reg_inc_beta(T& g1, T& g2, const T& a, const T& b, const T& z, const T& digammaA, const T& digammaB, const T& digammaSum, const T& betaAB) { - using std::exp; - - T c1 = log(z); - T c2 = log1m(z); - T c3 = betaAB * inc_beta(a, b, z); - T C = exp(a * c1 + b * c2) / a; - T dF1 = 0; - T dF2 = 0; - T dF3 = 0; - T dFz = 0; - if (value_of_rec(C)) { - std::forward_as_tuple(dF1, dF2, dF3, dFz) - = grad_2F1(a + b, 1.0, a + 1, z); - } - - T dBda = (c1 - 1.0 / a) * c3 + C * (dF1 + dF3); - T dBdb = c2 * c3 + C * dF1; - - g1 = (dBda - c3 * (digammaA - digammaSum)) / betaAB; - g2 = (dBdb - c3 * (digammaB - digammaSum)) / betaAB; + auto grads = grad_reg_inc_beta(a, b, z, digammaA, digammaB, + digammaSum, betaAB); + g1 = grads.first; + g2 = grads.second; } } // namespace math diff --git a/stan/math/prim/fun/inc_beta_dda.hpp b/stan/math/prim/fun/inc_beta_dda.hpp index 4ed2c0e4758..bf5d08258db 100644 --- a/stan/math/prim/fun/inc_beta_dda.hpp +++ b/stan/math/prim/fun/inc_beta_dda.hpp @@ -3,99 +3,84 @@ #include #include -#include -#include -#include -#include -#include -#include +#include +#include +#include +#include namespace stan { namespace math { /** - * Returns the partial derivative of the regularized - * incomplete beta function, I_{z}(a, b) with respect to a. - * The power series used to compute the derivative tends to - * converge slowly when a and b are large, especially if z - * approaches 1. The implementation will throw an exception - * if the series have not converged within 100,000 iterations. - * The current implementation has been tested for values - * of a and b up to 12500 and z = 0.999. + * Returns the partial derivative of the regularized incomplete beta function + * I_z(a, b) with respect to a. * - * @tparam T scalar types of arguments + * This calls grad_reg_inc_beta, which computes + * only the derivative with respect to a. + * + * @tparam T1 type of a + * @tparam T2 type of b + * @tparam T3 type of z + * @tparam T4 type of digamma_a + * @tparam T5 type of digamma_ab * @param a first argument * @param b second argument * @param z upper bound of the integral * @param digamma_a value of digamma(a) - * @param digamma_ab value of digamma(b) + * @param digamma_ab value of digamma(a + b) * @return partial derivative of the incomplete beta with respect to a * * @pre a >= 0 * @pre b >= 0 * @pre 0 <= z <= 1 */ -template -inline T inc_beta_dda(T a, T b, T z, T digamma_a, T digamma_ab) { - using std::fabs; - using std::log; - using std::pow; - - if (b > a) { - if ((0.1 < z && z <= 0.75 && b > 500) || (0.01 < z && z <= 0.1 && b > 2500) - || (0.001 < z && z <= 0.01 && b > 1e5)) { - return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab); - } - } - - if (z > 0.75 && a < 500) { - return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab); - } - if (z > 0.9 && a < 2500) { - return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab); - } - if (z > 0.99 && a < 1e5) { - return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab); - } - if (z > 0.999) { - return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab); - } - - double threshold = 1e-10; - - const T a_plus_b = a + b; - const T a_plus_1 = a + 1; - - digamma_a += inv(a); // Need digamma(a + 1), not digamma(a); - - // Common prefactor to regularize numerator and denominator - T prefactor = pow(a_plus_1 / a_plus_b, 3); - - T sum_numer = (digamma_ab - digamma_a) * prefactor; - T sum_denom = prefactor; - - T summand = prefactor * z * a_plus_b / a_plus_1; - - T k = 1; - digamma_ab += inv(a_plus_b); - digamma_a += inv(a_plus_1); - - while (fabs(summand) > threshold) { - sum_numer += (digamma_ab - digamma_a) * summand; - sum_denom += summand; - - summand *= (1 + (a_plus_b) / k) * (1 + k) / (1 + a_plus_1 / k); - digamma_ab += inv(a_plus_b + k); - digamma_a += inv(a_plus_1 + k); - ++k; - summand *= z / k; +template * = nullptr> +inline auto inc_beta_dda(const T1& a, const T2& b, const T3& z, + const T4& digamma_a, const T5& digamma_ab) { + return grad_reg_inc_beta(a, b, z, digamma_a, 0.0, digamma_ab, + 0.0); +} - if (k > 1e5) { - throw_domain_error("inc_beta_dda", - "did not converge within 10000 iterations", "", ""); - } +/** + * Enables the vectorized application of the inc_beta_dda function, when + * any arguments are vectors. Scalar arguments are broadcast against the + * vectors, which must all have the same size. + * + * @tparam T1 type of a + * @tparam T2 type of b + * @tparam T3 type of z + * @tparam T4 type of digamma_a + * @tparam T5 type of digamma_ab + * @param a first argument + * @param b second argument + * @param z upper bound of the integral + * @param digamma_a value of digamma(a) + * @param digamma_ab value of digamma(a + b) + * @return column array of the partial derivatives of the incomplete beta + * with respect to a + * @throw std::invalid_argument if the vector arguments differ in size + */ +template * = nullptr, + require_all_not_var_matrix_t* = nullptr> +inline auto inc_beta_dda(const T1& a, const T2& b, const T3& z, + const T4& digamma_a, const T5& digamma_ab) { + check_consistent_sizes("inc_beta_dda", "a", a, "b", b, "z", z, "digamma_a", + digamma_a, "digamma_ab", digamma_ab); + using ret_t = return_type_t; + scalar_seq_view a_vec(a); + scalar_seq_view b_vec(b); + scalar_seq_view z_vec(z); + scalar_seq_view digamma_a_vec(digamma_a); + scalar_seq_view digamma_ab_vec(digamma_ab); + const size_t N = max_size(a, b, z, digamma_a, digamma_ab); + Eigen::Array result(N); + for (size_t i = 0; i < N; ++i) { + result.coeffRef(i) = inc_beta_dda(a_vec[i], b_vec[i], z_vec[i], + digamma_a_vec[i], digamma_ab_vec[i]); } - return inc_beta(a, b, z) * (log(z) + sum_numer / sum_denom); + return result; } } // namespace math diff --git a/stan/math/prim/fun/inc_beta_ddb.hpp b/stan/math/prim/fun/inc_beta_ddb.hpp index 5514743f3e9..1e60ec1df01 100644 --- a/stan/math/prim/fun/inc_beta_ddb.hpp +++ b/stan/math/prim/fun/inc_beta_ddb.hpp @@ -2,92 +2,40 @@ #define STAN_MATH_PRIM_FUN_INC_BETA_DDB_HPP #include -#include -#include -#include -#include -#include -#include -#include +#include namespace stan { namespace math { -template -inline T inc_beta_dda(T a, T b, T z, T digamma_a, T digamma_ab); - /** - * Returns the partial derivative of the regularized - * incomplete beta function, I_{z}(a, b) with respect to b. - * The power series used to compute the derivative tends to - * converge slowly when a and b are large, especially if z - * approaches 1. The implementation will throw an exception - * if the series have not converged within 100,000 iterations. - * The current implementation has been tested for values - * of a and b up to 12500 and z = 0.999. + * Returns the partial derivative of the regularized incomplete beta function + * I_z(a, b) with respect to b. + * + * This calls grad_reg_inc_beta, which computes + * only the derivative with respect to b. * - * @tparam T scalar types of arguments + * @tparam T1 type of a + * @tparam T2 type of b + * @tparam T3 type of z + * @tparam T4 type of digamma_b + * @tparam T5 type of digamma_ab * @param a first argument * @param b second argument * @param z upper bound of the integral * @param digamma_b value of digamma(b) - * @param digamma_ab value of digamma(b) + * @param digamma_ab value of digamma(a + b) * @return partial derivative of the incomplete beta with respect to b * * @pre a >= 0 * @pre b >= 0 * @pre 0 <= z <= 1 */ -template -inline T inc_beta_ddb(T a, T b, T z, T digamma_b, T digamma_ab) { - using std::fabs; - using std::log; - using std::pow; - - if (b > a) { - if ((0.1 < z && z <= 0.75 && b > 500) || (0.01 < z && z <= 0.1 && b > 2500) - || (0.001 < z && z <= 0.01 && b > 1e5)) { - return -inc_beta_dda(b, a, 1 - z, digamma_b, digamma_ab); - } - } - - if ((z > 0.75 && a < 500) || (z > 0.9 && a < 2500) || (z > 0.99 && a < 1e5) - || (z > 0.999)) { - return -inc_beta_dda(b, a, 1 - z, digamma_b, digamma_ab); - } - - double threshold = 1e-10; - - const T a_plus_b = a + b; - const T a_plus_1 = a + 1; - - // Common prefactor to regularize numerator and denominator - T prefactor = pow(a_plus_1 / a_plus_b, 3); - - T sum_numer = digamma_ab * prefactor; - T sum_denom = prefactor; - - T summand = prefactor * z * a_plus_b / a_plus_1; - - T k = 1; - digamma_ab += inv(a_plus_b); - - while (fabs(summand) > threshold) { - sum_numer += digamma_ab * summand; - sum_denom += summand; - - summand *= (1 + (a_plus_b) / k) * (1 + k) / (1 + a_plus_1 / k); - digamma_ab += inv(a_plus_b + k); - ++k; - summand *= z / k; - - if (k > 1e5) { - throw_domain_error("inc_beta_ddb", - "did not converge within 100000 iterations", "", ""); - } - } - - return inc_beta(a, b, z) * (log1m(z) - digamma_b + sum_numer / sum_denom); +template * = nullptr> +inline auto inc_beta_ddb(const T1& a, const T2& b, const T3& z, + const T4& digamma_b, const T5& digamma_ab) { + return grad_reg_inc_beta(a, b, z, 0.0, digamma_b, digamma_ab, + 0.0); } } // namespace math diff --git a/stan/math/prim/fun/inc_beta_ddz.hpp b/stan/math/prim/fun/inc_beta_ddz.hpp index 4543df47880..4b955b2a189 100644 --- a/stan/math/prim/fun/inc_beta_ddz.hpp +++ b/stan/math/prim/fun/inc_beta_ddz.hpp @@ -6,6 +6,7 @@ #include #include #include +#include #include #include @@ -26,18 +27,58 @@ namespace math { * @pre b > 0 * @pre 0 < z <= 1 */ -template -inline T inc_beta_ddz(T a, T b, T z) { - using std::exp; - using std::log; +template * = nullptr> +inline auto inc_beta_ddz(T1 a, T2 b, T3 z) { return exp((b - 1) * log1m(z) + (a - 1) * log(z) + lgamma(a + b) - lgamma(a) - lgamma(b)); } -template <> +/** + * Returns the partial derivative of the regularized + * incomplete beta function, I_{z}(a, b) with respect to z. + * + * This overload for double arguments is evaluated by Boost. Where the density + * is infinite (z = 0 with a < 1, or z = 1 with b < 1) it returns infinity + * instead of throwing. + * + * @param a first argument + * @param b second argument + * @param z upper bound of the integral + * @return partial derivative of the incomplete beta with respect to z + */ inline double inc_beta_ddz(double a, double b, double z) { using boost::math::ibeta_derivative; - return ibeta_derivative(a, b, z); + using overflow_ignore + = boost::math::policies::policy>; + return ibeta_derivative(a, b, z, overflow_ignore()); +} + +/** + * Enables the vectorized application of the inc_beta_ddz function, when + * any arguments are containers. + * + * @tparam T1 type of first input + * @tparam T2 type of second input + * @tparam T3 type of third input + * @param a First input + * @param b Second input + * @param z Third input + * @return Partial derivative of the incomplete beta function with respect + * to z, applied to the three inputs. + */ +template * = nullptr, + require_all_not_var_matrix_t* = nullptr> +inline auto inc_beta_ddz(T1&& a, T2&& b, T3&& z) { + return apply_scalar_ternary( + [](auto&& d, auto&& e, auto&& f) { + return inc_beta_ddz(std::forward(d), + std::forward(e), + std::forward(f)); + }, + std::forward(a), std::forward(b), std::forward(z)); } } // namespace math diff --git a/test/unit/math/prim/fun/grad_reg_inc_beta_test.cpp b/test/unit/math/prim/fun/grad_reg_inc_beta_test.cpp index 129028fbab7..19b8e5c95dd 100644 --- a/test/unit/math/prim/fun/grad_reg_inc_beta_test.cpp +++ b/test/unit/math/prim/fun/grad_reg_inc_beta_test.cpp @@ -1,8 +1,51 @@ #include #include +#include #include +namespace { + +/** + * Both shape derivatives of I_z(a, b) against a high-precision reference + * (mpmath, 80 digits, density-form quadrature checked against numerical + * differentiation at 110 digits). Each one-sided request is also checked. + * It carries only its own derivative through the expansion, with the other + * digamma unused. + */ +void expect_grads(double a, double b, double z, double dda, double ddb, + double rtol) { + using stan::math::digamma; + using stan::math::lbeta; + double g1 = 0; + double g2 = 0; + stan::math::grad_reg_inc_beta(g1, g2, a, b, z, digamma(a), digamma(b), + digamma(a + b), std::exp(lbeta(a, b))); + EXPECT_TRUE(std::isfinite(g1)) << "a=" << a << " b=" << b << " z=" << z; + EXPECT_TRUE(std::isfinite(g2)) << "a=" << a << " b=" << b << " z=" << z; + EXPECT_NEAR(dda, g1, rtol * std::fabs(dda)) + << "a=" << a << " b=" << b << " z=" << z; + EXPECT_NEAR(ddb, g2, rtol * std::fabs(ddb)) + << "a=" << a << " b=" << b << " z=" << z; + + using stan::math::GRAD_RIB; + const double only_g1 = stan::math::grad_reg_inc_beta( + a, b, z, digamma(a), 0.0, digamma(a + b), 0.0); + EXPECT_NEAR(dda, only_g1, rtol * std::fabs(dda)) + << "G1 a=" << a << " b=" << b << " z=" << z; + const double only_g2 = stan::math::grad_reg_inc_beta( + a, b, z, 0.0, digamma(b), digamma(a + b), 0.0); + EXPECT_NEAR(ddb, only_g2, rtol * std::fabs(ddb)) + << "G2 a=" << a << " b=" << b << " z=" << z; + const auto both = stan::math::grad_reg_inc_beta( + a, b, z, digamma(a), digamma(b), digamma(a + b), 0.0); + EXPECT_EQ(g1, both.first) << "BOTH a=" << a << " b=" << b << " z=" << z; + EXPECT_EQ(g2, both.second) << "BOTH a=" << a << " b=" << b << " z=" << z; +} + +} // namespace + TEST(grad_reg_inc_beta, 1) { + // at y = 1 the function is identically 1 and both derivatives are 0 double alpha = 1.0; double beta = 1.0; double y = 1.0; @@ -16,22 +59,51 @@ TEST(grad_reg_inc_beta, 1) { stan::math::grad_reg_inc_beta(g1, g2, alpha, beta, y, digamma_alpha, digamma_beta, digamma_sum, betafunc); EXPECT_FLOAT_EQ(0, g1); - EXPECT_FLOAT_EQ(-std::numeric_limits::infinity(), g2); + EXPECT_FLOAT_EQ(0, g2); } TEST(grad_reg_inc_beta, 2) { - double alpha = 1.0; - double beta = 1.0; - double y = 0.4; - double digamma_alpha = stan::math::digamma(alpha); - double digamma_beta = stan::math::digamma(beta); - double digamma_sum = stan::math::digamma(alpha + beta); - double betafunc = std::exp(stan::math::lbeta(alpha, beta)); + expect_grads(1.0, 1.0, 0.4, -0.36651629274966202, 0.30649537425959442, 1e-12); +} - double g1 = 0; - double g2 = 0; - stan::math::grad_reg_inc_beta(g1, g2, alpha, beta, y, digamma_alpha, - digamma_beta, digamma_sum, betafunc); - EXPECT_NEAR(-0.36651629, g1, 1e-6); - EXPECT_NEAR(0.30649537, g2, 1e-6); +TEST(grad_reg_inc_beta, moderate) { + expect_grads(2, 3, 0.25, -0.22805360232434637, 0.10692153005229139, 1e-12); + expect_grads(0.5, 0.5, 0.5, -0.58312180806163756, 0.58312180806163756, 1e-12); + expect_grads(20, 20, 0.5, -0.063740468830439859, 0.063740468830439859, 1e-12); + expect_grads(1.5, 1.25, 0.6, -0.23806756196382869, 0.32279574755409757, + 1e-12); +} + +TEST(grad_reg_inc_beta, beta_ab_underflow) { + // beta(a, b) underflows for a + b above about 1100 + expect_grads(600, 600, 0.5, -0.01152047147304755, 0.01152047147304755, 1e-12); +} + +TEST(grad_reg_inc_beta, prefactor_underflow) { + // z^a (1 - z)^b underflows while I is 1 and the derivatives are 1e-138 + expect_grads(352, 590, 0.758, -3.3772717797123242e-138, + 4.5390666074606649e-138, 1e-12); +} + +/** + * For integer b the power series of the value terminates but the derivative + * series does not. The stop test must check both. + */ +TEST(grad_reg_inc_beta, integer_b) { + expect_grads(10, 1, 0.5, -0.00067690154351557159, 0.00226574342615196, 1e-12); + expect_grads(1, 100, 0.01, -0.42882413065860014, 0.0036787479630394138, + 1e-12); + expect_grads(2500, 2, 0.999, -0.00020509953516671709, 0.24850317856904407, + 1e-12); +} + +TEST(grad_reg_inc_beta, tails) { + expect_grads(1.145, 6786, 1.8e-4, -0.38621356822721375, 5.8439109664263782e-5, + 1e-10); + expect_grads(446, 5, 0.76, -1.0644675765185877e-47, 1.2763067943072555e-46, + 1e-12); + expect_grads(15000, 1.25, 0.999, -6.5954322569556673e-10, + 2.008497927784773e-6, 1e-10); + expect_grads(15000, 12500, 0.5, -8.7210244281937301e-53, + 9.5528279245811057e-53, 1e-10); } diff --git a/test/unit/math/prim/fun/inc_beta_dda_test.cpp b/test/unit/math/prim/fun/inc_beta_dda_test.cpp index ffe17385a6d..934a73c17f0 100644 --- a/test/unit/math/prim/fun/inc_beta_dda_test.cpp +++ b/test/unit/math/prim/fun/inc_beta_dda_test.cpp @@ -1,5 +1,23 @@ #include #include +#include + +namespace { + +/** + * d/da of I_z(a, b) against a high-precision reference (mpmath, 80 digits, + * density-form quadrature checked against numerical differentiation). + */ +void expect_dda(double a, double b, double z, double expected, double rtol) { + using stan::math::digamma; + using stan::math::inc_beta_dda; + const double got = inc_beta_dda(a, b, z, digamma(a), digamma(a + b)); + EXPECT_TRUE(std::isfinite(got)) << "a=" << a << " b=" << b << " z=" << z; + EXPECT_NEAR(expected, got, rtol * std::fabs(expected)) + << "a=" << a << " b=" << b << " z=" << z; +} + +} // namespace TEST(MathFunctions, inc_beta_dda) { using stan::math::digamma; @@ -15,49 +33,87 @@ TEST(MathFunctions, inc_beta_dda) { double mid_z = 0.6; double large_z = 0.999; - EXPECT_FLOAT_EQ(-0.00028665637, - inc_beta_dda(small_a, small_b, small_z, digamma(small_a), - digamma(small_a + small_b))) - << "reasonable values for a, b, x"; - EXPECT_FLOAT_EQ(-0.23806757, - inc_beta_dda(small_a, small_b, mid_z, digamma(small_a), - digamma(small_a + small_b))) - << "reasonable values for a, b, x"; - EXPECT_FLOAT_EQ(-0.00022264493, - inc_beta_dda(small_a, small_b, large_z, digamma(small_a), - digamma(small_a + small_b))) - << "reasonable values for a, b, x"; + expect_dda(small_a, small_b, small_z, -0.00028665636570426377, 1e-10); + expect_dda(small_a, small_b, mid_z, -0.23806756196382869, 1e-10); + expect_dda(small_a, small_b, large_z, -0.00022264492614285471, 1e-10); + // I underflows to 0 in double; the derivative is below 1e-3000 EXPECT_FLOAT_EQ(0.0, inc_beta_dda(large_a, small_b, small_z, digamma(large_a), - digamma(large_a + small_b))) - << "reasonable values for a, b, x"; + digamma(large_a + small_b))); EXPECT_FLOAT_EQ(0.0, inc_beta_dda(large_a, small_b, mid_z, digamma(large_a), - digamma(large_a + small_b))) - << "reasonable values for a, b, x"; - EXPECT_FLOAT_EQ(0.00091716705, - inc_beta_dda(large_a, small_b, large_z, digamma(large_a), - digamma(large_a + small_b))) - << "reasonable values for a, b, x"; - - EXPECT_FLOAT_EQ(1.8226241, - inc_beta_dda(small_a, large_b, small_z, digamma(small_a), - digamma(small_a + large_b))) - << "reasonable values for a, b, x"; + digamma(large_a + small_b))); + expect_dda(large_a, small_b, large_z, -6.5954322569556673e-10, 1e-10); + + expect_dda(small_a, large_b, small_z, -3.937563719959171e-5, 1e-10); EXPECT_FLOAT_EQ(0.0, inc_beta_dda(small_a, large_b, mid_z, digamma(small_a), - digamma(small_a + large_b))) - << "reasonable values for a, b, x"; + digamma(small_a + large_b))); EXPECT_FLOAT_EQ(0.0, inc_beta_dda(small_a, large_b, large_z, digamma(small_a), - digamma(small_a + large_b))) - << "reasonable values for a, b, x"; + digamma(small_a + large_b))); EXPECT_FLOAT_EQ(0.0, inc_beta_dda(large_a, large_b, small_z, digamma(large_a), - digamma(large_a + large_b))) - << "reasonable values for a, b, x"; - EXPECT_FLOAT_EQ(-4.2743586e-14, - inc_beta_dda(large_a, large_b, mid_z, digamma(large_a), - digamma(large_a + large_b))) - << "reasonable values for a, b, x"; + digamma(large_a + large_b))); + expect_dda(large_a, large_b, mid_z, -1.755811372560583e-76, 1e-10); EXPECT_FLOAT_EQ(0.0, inc_beta_dda(large_a, large_b, large_z, digamma(large_a), - digamma(large_a + large_b))) - << "reasonable values for a, b, x"; + digamma(large_a + large_b))); +} + +TEST(MathFunctions, inc_beta_dda_regressions) { + // b >> a with small z + expect_dda(1.145, 6786, 1.8e-4, -0.38621356822721375, 1e-10); + // deep lower tail, I of order 1e-47 + expect_dda(446, 5, 0.76, -1.0644675765185877e-47, 1e-10); + // large a with z near 1 + expect_dda(2500, 2, 0.999, -0.00020509953516671709, 1e-10); + expect_dda(1, 100, 0.01, -0.42882413065860014, 1e-10); + // integer b: the value series terminates but the derivative does not + expect_dda(2, 3, 0.25, -0.22805360232434637, 1e-10); + expect_dda(10, 1, 0.5, -0.00067690154351557159, 1e-10); + // a + b above the beta(a, b) underflow limit + expect_dda(600, 600, 0.5, -0.01152047147304755, 1e-10); + expect_dda(0.5, 0.5, 0.5, -0.58312180806163756, 1e-10); + expect_dda(20, 20, 0.5, -0.063740468830439859, 1e-10); +} + +TEST(MathFunctions, inc_beta_dda_vectorized) { + using stan::math::digamma; + using stan::math::inc_beta_dda; + + Eigen::VectorXd a(4); + a << 1.145, 446, 2, 600; + Eigen::VectorXd b(4); + b << 6786, 5, 3, 600; + Eigen::VectorXd z(4); + z << 1.8e-4, 0.76, 0.25, 0.5; + Eigen::VectorXd digamma_a(4); + Eigen::VectorXd digamma_ab(4); + for (int i = 0; i < 4; ++i) { + digamma_a[i] = digamma(a[i]); + digamma_ab[i] = digamma(a[i] + b[i]); + } + + Eigen::ArrayXd all_eigen = inc_beta_dda(a, b, z, digamma_a, digamma_ab); + ASSERT_EQ(4, all_eigen.size()); + for (int i = 0; i < 4; ++i) { + EXPECT_EQ(inc_beta_dda(a[i], b[i], z[i], digamma_a[i], digamma_ab[i]), + all_eigen[i]); + } + + std::vector z_std{1.8e-4, 0.76, 0.25, 0.5}; + Eigen::ArrayXd std_z = inc_beta_dda(a, b, z_std, digamma_a, digamma_ab); + for (int i = 0; i < 4; ++i) { + EXPECT_EQ(all_eigen[i], std_z[i]); + } + + // scalar shapes broadcast against a vector of upper limits + Eigen::ArrayXd only_z = inc_beta_dda(2.0, 3.0, z, digamma(2.0), digamma(5.0)); + ASSERT_EQ(4, only_z.size()); + for (int i = 0; i < 4; ++i) { + EXPECT_EQ(inc_beta_dda(2.0, 3.0, z[i], digamma(2.0), digamma(5.0)), + only_z[i]); + } + + Eigen::VectorXd too_short(3); + too_short << 0.1, 0.2, 0.3; + EXPECT_THROW(inc_beta_dda(a, b, too_short, digamma_a, digamma_ab), + std::invalid_argument); } diff --git a/test/unit/math/prim/fun/inc_beta_ddb_test.cpp b/test/unit/math/prim/fun/inc_beta_ddb_test.cpp index d6116dd437a..ec6c28928cc 100644 --- a/test/unit/math/prim/fun/inc_beta_ddb_test.cpp +++ b/test/unit/math/prim/fun/inc_beta_ddb_test.cpp @@ -1,5 +1,24 @@ #include #include +#include + +namespace { + +/** + * d/db of I_z(a, b) against a high-precision reference (mpmath, 80 digits, + * density-form quadrature checked against numerical differentiation). + * The fourth argument of inc_beta_ddb is digamma(b). + */ +void expect_ddb(double a, double b, double z, double expected, double rtol) { + using stan::math::digamma; + using stan::math::inc_beta_ddb; + const double got = inc_beta_ddb(a, b, z, digamma(b), digamma(a + b)); + EXPECT_TRUE(std::isfinite(got)) << "a=" << a << " b=" << b << " z=" << z; + EXPECT_NEAR(expected, got, rtol * std::fabs(expected)) + << "a=" << a << " b=" << b << " z=" << z; +} + +} // namespace TEST(MathFunctions, inc_beta_ddb) { using stan::math::digamma; @@ -15,48 +34,39 @@ TEST(MathFunctions, inc_beta_ddb) { double mid_z = 0.5; double large_z = 0.999; - EXPECT_FLOAT_EQ(3.2996082e-05, - inc_beta_ddb(small_a, small_b, small_z, digamma(small_a), - digamma(small_a + small_b))) - << "reasonable values for a, b, x"; - EXPECT_FLOAT_EQ(0.17945045, - inc_beta_ddb(small_a, small_b, mid_z, digamma(small_a), - digamma(small_a + small_b))) - << "reasonable values for a, b, x"; - EXPECT_FLOAT_EQ(0.0019721228, - inc_beta_ddb(small_a, small_b, large_z, digamma(small_a), - digamma(small_a + small_b))) - << "reasonable values for a, b, x"; - - EXPECT_FLOAT_EQ(0.0, inc_beta_ddb(large_a, small_b, small_z, digamma(large_a), - digamma(large_a + small_b))) - << "reasonable values for a, b, x"; - EXPECT_FLOAT_EQ(0.0, inc_beta_ddb(large_a, small_b, mid_z, digamma(large_a), - digamma(large_a + small_b))) - << "reasonable values for a, b, x"; - EXPECT_FLOAT_EQ(7.7076669, - inc_beta_ddb(large_a, small_b, large_z, digamma(large_a), - digamma(large_a + small_b))) - << "reasonable values for a, b, x"; - - EXPECT_FLOAT_EQ(9.3959293, - inc_beta_ddb(small_a, large_b, small_z, digamma(small_a), - digamma(small_a + large_b))) - << "reasonable values for a, b, x"; - EXPECT_FLOAT_EQ(0.0, inc_beta_ddb(small_a, large_b, mid_z, digamma(small_a), - digamma(small_a + large_b))) - << "reasonable values for a, b, x"; - EXPECT_FLOAT_EQ(0.0, inc_beta_ddb(small_a, large_b, large_z, digamma(small_a), - digamma(small_a + large_b))) - << "reasonable values for a, b, x"; - - EXPECT_FLOAT_EQ(0.0, inc_beta_ddb(large_a, large_b, small_z, digamma(large_a), - digamma(large_a + large_b))) - << "reasonable values for a, b, x"; - EXPECT_FLOAT_EQ(0.0, inc_beta_ddb(large_a, large_b, mid_z, digamma(large_a), - digamma(large_a + large_b))) - << "reasonable values for a, b, x"; - EXPECT_FLOAT_EQ(0.0, inc_beta_ddb(large_a, large_b, large_z, digamma(large_a), - digamma(large_a + large_b))) - << "reasonable values for a, b, x"; + expect_ddb(small_a, small_b, small_z, 4.4135732827365698e-5, 1e-10); + expect_ddb(small_a, small_b, mid_z, 0.29301795311347039, 1e-10); + expect_ddb(small_a, small_b, large_z, 0.0018969609780283825, 1e-10); + + // I underflows to 0 in double; the derivative is below 1e-3000 + EXPECT_FLOAT_EQ(0.0, inc_beta_ddb(large_a, small_b, small_z, digamma(small_b), + digamma(large_a + small_b))); + EXPECT_FLOAT_EQ(0.0, inc_beta_ddb(large_a, small_b, mid_z, digamma(small_b), + digamma(large_a + small_b))); + expect_ddb(large_a, small_b, large_z, 2.008497927784773e-6, 1e-10); + + expect_ddb(small_a, large_b, small_z, 1.4782043391986094e-8, 1e-10); + EXPECT_FLOAT_EQ(0.0, inc_beta_ddb(small_a, large_b, mid_z, digamma(large_b), + digamma(small_a + large_b))); + EXPECT_FLOAT_EQ(0.0, inc_beta_ddb(small_a, large_b, large_z, digamma(large_b), + digamma(small_a + large_b))); + + EXPECT_FLOAT_EQ(0.0, inc_beta_ddb(large_a, large_b, small_z, digamma(large_b), + digamma(large_a + large_b))); + expect_ddb(large_a, large_b, mid_z, 9.5528279245811057e-53, 1e-10); + EXPECT_FLOAT_EQ(0.0, inc_beta_ddb(large_a, large_b, large_z, digamma(large_b), + digamma(large_a + large_b))); +} + +TEST(MathFunctions, inc_beta_ddb_regressions) { + expect_ddb(1.145, 6786, 1.8e-4, 5.8439109664263782e-5, 1e-10); + expect_ddb(446, 5, 0.76, 1.2763067943072555e-46, 1e-10); + expect_ddb(2500, 2, 0.999, 0.24850317856904407, 1e-10); + expect_ddb(1, 100, 0.01, 0.0036787479630394138, 1e-10); + // integer b: the value series terminates but the derivative does not + expect_ddb(2, 3, 0.25, 0.10692153005229139, 1e-10); + expect_ddb(10, 1, 0.5, 0.00226574342615196, 1e-10); + expect_ddb(600, 600, 0.5, 0.01152047147304755, 1e-10); + expect_ddb(0.5, 0.5, 0.5, 0.58312180806163756, 1e-10); + expect_ddb(20, 20, 0.5, 0.063740468830439859, 1e-10); } diff --git a/test/unit/math/prim/fun/inc_beta_ddz_test.cpp b/test/unit/math/prim/fun/inc_beta_ddz_test.cpp index f8aa0579778..d44b613c015 100644 --- a/test/unit/math/prim/fun/inc_beta_ddz_test.cpp +++ b/test/unit/math/prim/fun/inc_beta_ddz_test.cpp @@ -42,3 +42,51 @@ TEST(MathFunctions, inc_beta_ddz) { EXPECT_FLOAT_EQ(0.0, inc_beta_ddz(large_a, large_b, large_z)) << "reasonable values for a, b, x"; } + +TEST(MathFunctions, inc_beta_ddz_infinite_density) { + using stan::math::inc_beta_ddz; + // the density is infinite at z = 1 with b < 1 and at z = 0 with a < 1 + double at_one = 0; + double at_zero = 0; + EXPECT_NO_THROW(at_one = inc_beta_ddz(1.05, 0.35, 1.0)); + EXPECT_NO_THROW(at_zero = inc_beta_ddz(0.5, 2.0, 0.0)); + EXPECT_TRUE(std::isinf(at_one) && at_one > 0); + EXPECT_TRUE(std::isinf(at_zero) && at_zero > 0); +} + +TEST(MathFunctions, inc_beta_ddz_vectorized) { + using stan::math::inc_beta_ddz; + + Eigen::VectorXd a(3); + a << 1.5, 15.0, 0.75; + Eigen::VectorXd b(3); + b << 1.25, 4.0, 30.0; + Eigen::VectorXd z(3); + z << 0.001, 0.5, 0.999; + + Eigen::VectorXd all_eigen = inc_beta_ddz(a, b, z); + ASSERT_EQ(3, all_eigen.size()); + for (int i = 0; i < 3; ++i) { + EXPECT_EQ(inc_beta_ddz(a[i], b[i], z[i]), all_eigen[i]); + } + + std::vector a_std{1.5, 15.0, 0.75}; + std::vector b_std{1.25, 4.0, 30.0}; + std::vector z_std{0.001, 0.5, 0.999}; + std::vector all_std = inc_beta_ddz(a_std, b_std, z_std); + ASSERT_EQ(3, all_std.size()); + for (int i = 0; i < 3; ++i) { + EXPECT_EQ(inc_beta_ddz(a_std[i], b_std[i], z_std[i]), all_std[i]); + } + + Eigen::VectorXd scalar_a = inc_beta_ddz(2.5, b, z); + Eigen::VectorXd scalar_b = inc_beta_ddz(a, 3.5, z); + Eigen::VectorXd scalar_z = inc_beta_ddz(a, b, 0.3); + Eigen::VectorXd only_z = inc_beta_ddz(2.5, 3.5, z); + for (int i = 0; i < 3; ++i) { + EXPECT_EQ(inc_beta_ddz(2.5, b[i], z[i]), scalar_a[i]); + EXPECT_EQ(inc_beta_ddz(a[i], 3.5, z[i]), scalar_b[i]); + EXPECT_EQ(inc_beta_ddz(a[i], b[i], 0.3), scalar_z[i]); + EXPECT_EQ(inc_beta_ddz(2.5, 3.5, z[i]), only_z[i]); + } +} From 7308fa514b5b4b73b7e22ccc0f34e87e50b0a1e7 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Wed, 28 Oct 2026 23:26:59 +0200 Subject: [PATCH 2/2] Split the overload tests, add throw and lcdf tests, and restore the TFP notice --- licenses/tensorflow-probability-license.txt | 201 ++++++++++++++++++ stan/math/prim/fun/grad_reg_inc_beta.hpp | 19 +- .../fun/grad_reg_inc_beta_one_sided_test.cpp | 76 +++++++ .../math/prim/fun/grad_reg_inc_beta_test.cpp | 27 +-- test/unit/math/prim/fun/inc_beta_dda_test.cpp | 53 +---- .../prim/fun/inc_beta_dda_vectorized_test.cpp | 47 ++++ test/unit/math/prim/fun/inc_beta_ddb_test.cpp | 13 ++ test/unit/math/prim/fun/inc_beta_ddz_test.cpp | 37 ---- .../prim/fun/inc_beta_ddz_vectorized_test.cpp | 40 ++++ test/unit/math/rev/prob/beta_cdf_log_test.cpp | 46 ++++ .../rev/prob/neg_binomial_cdf_log_test.cpp | 10 + 11 files changed, 471 insertions(+), 98 deletions(-) create mode 100644 licenses/tensorflow-probability-license.txt create mode 100644 test/unit/math/prim/fun/grad_reg_inc_beta_one_sided_test.cpp create mode 100644 test/unit/math/prim/fun/inc_beta_dda_vectorized_test.cpp create mode 100644 test/unit/math/prim/fun/inc_beta_ddz_vectorized_test.cpp create mode 100644 test/unit/math/rev/prob/beta_cdf_log_test.cpp diff --git a/licenses/tensorflow-probability-license.txt b/licenses/tensorflow-probability-license.txt new file mode 100644 index 00000000000..261eeb9e9f8 --- /dev/null +++ b/licenses/tensorflow-probability-license.txt @@ -0,0 +1,201 @@ + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + + Copyright [yyyy] [name of copyright owner] + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. diff --git a/stan/math/prim/fun/grad_reg_inc_beta.hpp b/stan/math/prim/fun/grad_reg_inc_beta.hpp index 92fab1fc056..9fc4b7e1265 100644 --- a/stan/math/prim/fun/grad_reg_inc_beta.hpp +++ b/stan/math/prim/fun/grad_reg_inc_beta.hpp @@ -1,3 +1,10 @@ +// Original code derived from TensorFlow Probability and is distributed here +// under the Apache License 2.0 (licenses/tensorflow-probability-license.txt) +// Copyright 2020 The TensorFlow Probability Authors. +// tensorflow_probability/python/math/special.py, betainc partials +// Secondary code copyright by its author and is distributed here +// under the BSD-3 license (LICENSE.md) + #ifndef STAN_MATH_PRIM_FUN_GRAD_REG_INC_BETA_HPP #define STAN_MATH_PRIM_FUN_GRAD_REG_INC_BETA_HPP @@ -56,7 +63,8 @@ namespace internal { * test therefore also requires dt_n/db to be small whenever that sum * matters. That is the case when d/db is requested, and also when b is an * autodiff type, because autodiff then differentiates the d/da result in b - * and needs the same terms. + * and needs the same terms. This test of dt_n/db is a change from the + * TensorFlow Probability source, which tests t_n only. * * @tparam grads_requested which partials to compute * @tparam flip_sign true when the caller evaluates its I_z(a, b) through the @@ -117,6 +125,7 @@ inline auto inc_beta_power_series_kernel(const T1& a, const T2& b, const T3& x, } product = new_product; series_sum += term; + // changed from TFP, which tests the term only: db_converged added if (fabs(value_of_rec(term)) <= tolerance && db_converged) { break; } @@ -199,7 +208,8 @@ inline auto inc_beta_power_series_kernel(const T1& a, const T2& b, const T3& x, * a derivative no longer change. The b derivative of d_{2b} is b x / den, * which is not 0, so dh/db is still changing. The b derivative is therefore * tested whenever it matters, which is when d/db is requested and also when b - * is an autodiff type, as in the power series. + * is an autodiff type, as in the power series. The tests of the derivatives + * are a change from the TensorFlow Probability source, which tests C D only. * * @tparam grads_requested which partials to compute * @tparam flip_sign true when the caller evaluates its I_z(a, b) through the @@ -316,6 +326,7 @@ inline auto inc_beta_continued_fraction_kernel(const T1& a, const T2& b, const double h_d = value_of_rec(h); const double dh_a_now = h_d * value_of_rec(r_a); const double dh_b_now = h_d * value_of_rec(r_b); + // changed from TFP, which tests C D only: the derivative tests added if (fabs(value_of_rec(delta) - 1.0) < tolerance && (grads_requested == GRAD_RIB::G2 || fabs(dh_a_now - dh_a_prev) <= tolerance * fabs(dh_a_now)) @@ -394,6 +405,10 @@ inline auto inc_beta_continued_fraction_kernel(const T1& a, const T2& b, * as the threshold because it exceeds 1 for b < 1. The fraction of I_z(a, b) * would then be used for z near 1, where it needs hundreds of steps. * + * Both choices are changes from the TensorFlow Probability source, which + * uses the power series up to z = 0.95 (from z = 0.05 for I_{1-z}(b, a)) + * and switches the continued fraction at the mode. + * * The partials are 0 at z = 0 and z = 1. * * @tparam grads_requested G1 for the derivative with respect to diff --git a/test/unit/math/prim/fun/grad_reg_inc_beta_one_sided_test.cpp b/test/unit/math/prim/fun/grad_reg_inc_beta_one_sided_test.cpp new file mode 100644 index 00000000000..449935fa8df --- /dev/null +++ b/test/unit/math/prim/fun/grad_reg_inc_beta_one_sided_test.cpp @@ -0,0 +1,76 @@ +#include +#include +#include + +namespace { + +/** + * The GRAD_RIB requests of grad_reg_inc_beta against the reference values + * of grad_reg_inc_beta_test.cpp (mpmath, 80 digits). G1 and G2 carry only + * their own derivative, with the other digamma unused. The comparison of + * BOTH with the 9-argument form checks that the latter forwards to BOTH. + */ +void expect_one_sided(double a, double b, double z, double dda, double ddb, + double rtol) { + using stan::math::digamma; + using stan::math::GRAD_RIB; + const double only_g1 = stan::math::grad_reg_inc_beta( + a, b, z, digamma(a), 0.0, digamma(a + b), 0.0); + EXPECT_NEAR(dda, only_g1, rtol * std::fabs(dda)) + << "G1 a=" << a << " b=" << b << " z=" << z; + const double only_g2 = stan::math::grad_reg_inc_beta( + a, b, z, 0.0, digamma(b), digamma(a + b), 0.0); + EXPECT_NEAR(ddb, only_g2, rtol * std::fabs(ddb)) + << "G2 a=" << a << " b=" << b << " z=" << z; + + double g1 = 0; + double g2 = 0; + stan::math::grad_reg_inc_beta(g1, g2, a, b, z, digamma(a), digamma(b), + digamma(a + b), 0.0); + const auto both = stan::math::grad_reg_inc_beta( + a, b, z, digamma(a), digamma(b), digamma(a + b), 0.0); + EXPECT_EQ(g1, both.first) << "BOTH a=" << a << " b=" << b << " z=" << z; + EXPECT_EQ(g2, both.second) << "BOTH a=" << a << " b=" << b << " z=" << z; +} + +} // namespace + +TEST(grad_reg_inc_beta_one_sided, moderate) { + expect_one_sided(1.0, 1.0, 0.4, -0.36651629274966202, 0.30649537425959442, + 1e-12); + expect_one_sided(2, 3, 0.25, -0.22805360232434637, 0.10692153005229139, + 1e-12); + expect_one_sided(0.5, 0.5, 0.5, -0.58312180806163756, 0.58312180806163756, + 1e-12); + expect_one_sided(20, 20, 0.5, -0.063740468830439859, 0.063740468830439859, + 1e-12); + expect_one_sided(1.5, 1.25, 0.6, -0.23806756196382869, 0.32279574755409757, + 1e-12); +} + +TEST(grad_reg_inc_beta_one_sided, underflow) { + expect_one_sided(600, 600, 0.5, -0.01152047147304755, 0.01152047147304755, + 1e-12); + expect_one_sided(352, 590, 0.758, -3.3772717797123242e-138, + 4.5390666074606649e-138, 1e-12); +} + +TEST(grad_reg_inc_beta_one_sided, integer_b) { + expect_one_sided(10, 1, 0.5, -0.00067690154351557159, 0.00226574342615196, + 1e-12); + expect_one_sided(1, 100, 0.01, -0.42882413065860014, 0.0036787479630394138, + 1e-12); + expect_one_sided(2500, 2, 0.999, -0.00020509953516671709, 0.24850317856904407, + 1e-12); +} + +TEST(grad_reg_inc_beta_one_sided, tails) { + expect_one_sided(1.145, 6786, 1.8e-4, -0.38621356822721375, + 5.8439109664263782e-5, 1e-10); + expect_one_sided(446, 5, 0.76, -1.0644675765185877e-47, + 1.2763067943072555e-46, 1e-12); + expect_one_sided(15000, 1.25, 0.999, -6.5954322569556673e-10, + 2.008497927784773e-6, 1e-10); + expect_one_sided(15000, 12500, 0.5, -8.7210244281937301e-53, + 9.5528279245811057e-53, 1e-10); +} diff --git a/test/unit/math/prim/fun/grad_reg_inc_beta_test.cpp b/test/unit/math/prim/fun/grad_reg_inc_beta_test.cpp index 19b8e5c95dd..68da5cf915f 100644 --- a/test/unit/math/prim/fun/grad_reg_inc_beta_test.cpp +++ b/test/unit/math/prim/fun/grad_reg_inc_beta_test.cpp @@ -8,9 +8,8 @@ namespace { /** * Both shape derivatives of I_z(a, b) against a high-precision reference * (mpmath, 80 digits, density-form quadrature checked against numerical - * differentiation at 110 digits). Each one-sided request is also checked. - * It carries only its own derivative through the expansion, with the other - * digamma unused. + * differentiation at 110 digits). The one-sided GRAD_RIB requests are + * tested in grad_reg_inc_beta_one_sided_test.cpp. */ void expect_grads(double a, double b, double z, double dda, double ddb, double rtol) { @@ -26,20 +25,6 @@ void expect_grads(double a, double b, double z, double dda, double ddb, << "a=" << a << " b=" << b << " z=" << z; EXPECT_NEAR(ddb, g2, rtol * std::fabs(ddb)) << "a=" << a << " b=" << b << " z=" << z; - - using stan::math::GRAD_RIB; - const double only_g1 = stan::math::grad_reg_inc_beta( - a, b, z, digamma(a), 0.0, digamma(a + b), 0.0); - EXPECT_NEAR(dda, only_g1, rtol * std::fabs(dda)) - << "G1 a=" << a << " b=" << b << " z=" << z; - const double only_g2 = stan::math::grad_reg_inc_beta( - a, b, z, 0.0, digamma(b), digamma(a + b), 0.0); - EXPECT_NEAR(ddb, only_g2, rtol * std::fabs(ddb)) - << "G2 a=" << a << " b=" << b << " z=" << z; - const auto both = stan::math::grad_reg_inc_beta( - a, b, z, digamma(a), digamma(b), digamma(a + b), 0.0); - EXPECT_EQ(g1, both.first) << "BOTH a=" << a << " b=" << b << " z=" << z; - EXPECT_EQ(g2, both.second) << "BOTH a=" << a << " b=" << b << " z=" << z; } } // namespace @@ -107,3 +92,11 @@ TEST(grad_reg_inc_beta, tails) { expect_grads(15000, 12500, 0.5, -8.7210244281937301e-53, 9.5528279245811057e-53, 1e-10); } + +TEST(grad_reg_inc_beta, upper_tail_no_throw) { + // deep upper tail: I rounds to 1 and the derivatives are near the smallest + // normal double + EXPECT_NO_THROW(expect_grads(4.0845138072907865, 150.97164240761546, + 0.991606836273885, -1.1190795130668706e-307, + 1.4196439874655394e-307, 1e-10)); +} diff --git a/test/unit/math/prim/fun/inc_beta_dda_test.cpp b/test/unit/math/prim/fun/inc_beta_dda_test.cpp index 934a73c17f0..4061ac4039b 100644 --- a/test/unit/math/prim/fun/inc_beta_dda_test.cpp +++ b/test/unit/math/prim/fun/inc_beta_dda_test.cpp @@ -74,46 +74,15 @@ TEST(MathFunctions, inc_beta_dda_regressions) { expect_dda(20, 20, 0.5, -0.063740468830439859, 1e-10); } -TEST(MathFunctions, inc_beta_dda_vectorized) { - using stan::math::digamma; - using stan::math::inc_beta_dda; - - Eigen::VectorXd a(4); - a << 1.145, 446, 2, 600; - Eigen::VectorXd b(4); - b << 6786, 5, 3, 600; - Eigen::VectorXd z(4); - z << 1.8e-4, 0.76, 0.25, 0.5; - Eigen::VectorXd digamma_a(4); - Eigen::VectorXd digamma_ab(4); - for (int i = 0; i < 4; ++i) { - digamma_a[i] = digamma(a[i]); - digamma_ab[i] = digamma(a[i] + b[i]); - } - - Eigen::ArrayXd all_eigen = inc_beta_dda(a, b, z, digamma_a, digamma_ab); - ASSERT_EQ(4, all_eigen.size()); - for (int i = 0; i < 4; ++i) { - EXPECT_EQ(inc_beta_dda(a[i], b[i], z[i], digamma_a[i], digamma_ab[i]), - all_eigen[i]); - } - - std::vector z_std{1.8e-4, 0.76, 0.25, 0.5}; - Eigen::ArrayXd std_z = inc_beta_dda(a, b, z_std, digamma_a, digamma_ab); - for (int i = 0; i < 4; ++i) { - EXPECT_EQ(all_eigen[i], std_z[i]); - } - - // scalar shapes broadcast against a vector of upper limits - Eigen::ArrayXd only_z = inc_beta_dda(2.0, 3.0, z, digamma(2.0), digamma(5.0)); - ASSERT_EQ(4, only_z.size()); - for (int i = 0; i < 4; ++i) { - EXPECT_EQ(inc_beta_dda(2.0, 3.0, z[i], digamma(2.0), digamma(5.0)), - only_z[i]); - } - - Eigen::VectorXd too_short(3); - too_short << 0.1, 0.2, 0.3; - EXPECT_THROW(inc_beta_dda(a, b, too_short, digamma_a, digamma_ab), - std::invalid_argument); +TEST(MathFunctions, inc_beta_dda_upper_tail_no_throw) { + // deep upper tail of large shapes: I rounds to 1 and the derivative is + // subnormal, so the tolerance allows for its 12 significant digits + const double a = 936.5962745738543; + const double b = 1294.2043377235907; + const double z = 0.7877755442980134; + double got = 0; + EXPECT_NO_THROW(got + = stan::math::inc_beta_dda(a, b, z, stan::math::digamma(a), + stan::math::digamma(a + b))); + EXPECT_NEAR(-3.941621974286e-312, got, 1e-9 * 3.941621974286e-312); } diff --git a/test/unit/math/prim/fun/inc_beta_dda_vectorized_test.cpp b/test/unit/math/prim/fun/inc_beta_dda_vectorized_test.cpp new file mode 100644 index 00000000000..9fda8c01691 --- /dev/null +++ b/test/unit/math/prim/fun/inc_beta_dda_vectorized_test.cpp @@ -0,0 +1,47 @@ +#include +#include +#include + +TEST(MathFunctions, inc_beta_dda_vectorized) { + using stan::math::digamma; + using stan::math::inc_beta_dda; + + Eigen::VectorXd a(4); + a << 1.145, 446, 2, 600; + Eigen::VectorXd b(4); + b << 6786, 5, 3, 600; + Eigen::VectorXd z(4); + z << 1.8e-4, 0.76, 0.25, 0.5; + Eigen::VectorXd digamma_a(4); + Eigen::VectorXd digamma_ab(4); + for (int i = 0; i < 4; ++i) { + digamma_a[i] = digamma(a[i]); + digamma_ab[i] = digamma(a[i] + b[i]); + } + + Eigen::ArrayXd all_eigen = inc_beta_dda(a, b, z, digamma_a, digamma_ab); + ASSERT_EQ(4, all_eigen.size()); + for (int i = 0; i < 4; ++i) { + EXPECT_EQ(inc_beta_dda(a[i], b[i], z[i], digamma_a[i], digamma_ab[i]), + all_eigen[i]); + } + + std::vector z_std{1.8e-4, 0.76, 0.25, 0.5}; + Eigen::ArrayXd std_z = inc_beta_dda(a, b, z_std, digamma_a, digamma_ab); + for (int i = 0; i < 4; ++i) { + EXPECT_EQ(all_eigen[i], std_z[i]); + } + + // scalar shapes broadcast against a vector of upper limits + Eigen::ArrayXd only_z = inc_beta_dda(2.0, 3.0, z, digamma(2.0), digamma(5.0)); + ASSERT_EQ(4, only_z.size()); + for (int i = 0; i < 4; ++i) { + EXPECT_EQ(inc_beta_dda(2.0, 3.0, z[i], digamma(2.0), digamma(5.0)), + only_z[i]); + } + + Eigen::VectorXd too_short(3); + too_short << 0.1, 0.2, 0.3; + EXPECT_THROW(inc_beta_dda(a, b, too_short, digamma_a, digamma_ab), + std::invalid_argument); +} diff --git a/test/unit/math/prim/fun/inc_beta_ddb_test.cpp b/test/unit/math/prim/fun/inc_beta_ddb_test.cpp index ec6c28928cc..0d242fd97e6 100644 --- a/test/unit/math/prim/fun/inc_beta_ddb_test.cpp +++ b/test/unit/math/prim/fun/inc_beta_ddb_test.cpp @@ -70,3 +70,16 @@ TEST(MathFunctions, inc_beta_ddb_regressions) { expect_ddb(0.5, 0.5, 0.5, 0.58312180806163756, 1e-10); expect_ddb(20, 20, 0.5, 0.063740468830439859, 1e-10); } + +TEST(MathFunctions, inc_beta_ddb_upper_tail_no_throw) { + // deep upper tail of large shapes: I rounds to 1 and the derivative is + // subnormal, so the tolerance allows for its 12 significant digits + const double a = 936.5962745738543; + const double b = 1294.2043377235907; + const double z = 0.7877755442980134; + double got = 0; + EXPECT_NO_THROW(got + = stan::math::inc_beta_ddb(a, b, z, stan::math::digamma(b), + stan::math::digamma(a + b))); + EXPECT_NEAR(6.29796155252e-312, got, 1e-9 * 6.29796155252e-312); +} diff --git a/test/unit/math/prim/fun/inc_beta_ddz_test.cpp b/test/unit/math/prim/fun/inc_beta_ddz_test.cpp index d44b613c015..16d16ca7c28 100644 --- a/test/unit/math/prim/fun/inc_beta_ddz_test.cpp +++ b/test/unit/math/prim/fun/inc_beta_ddz_test.cpp @@ -53,40 +53,3 @@ TEST(MathFunctions, inc_beta_ddz_infinite_density) { EXPECT_TRUE(std::isinf(at_one) && at_one > 0); EXPECT_TRUE(std::isinf(at_zero) && at_zero > 0); } - -TEST(MathFunctions, inc_beta_ddz_vectorized) { - using stan::math::inc_beta_ddz; - - Eigen::VectorXd a(3); - a << 1.5, 15.0, 0.75; - Eigen::VectorXd b(3); - b << 1.25, 4.0, 30.0; - Eigen::VectorXd z(3); - z << 0.001, 0.5, 0.999; - - Eigen::VectorXd all_eigen = inc_beta_ddz(a, b, z); - ASSERT_EQ(3, all_eigen.size()); - for (int i = 0; i < 3; ++i) { - EXPECT_EQ(inc_beta_ddz(a[i], b[i], z[i]), all_eigen[i]); - } - - std::vector a_std{1.5, 15.0, 0.75}; - std::vector b_std{1.25, 4.0, 30.0}; - std::vector z_std{0.001, 0.5, 0.999}; - std::vector all_std = inc_beta_ddz(a_std, b_std, z_std); - ASSERT_EQ(3, all_std.size()); - for (int i = 0; i < 3; ++i) { - EXPECT_EQ(inc_beta_ddz(a_std[i], b_std[i], z_std[i]), all_std[i]); - } - - Eigen::VectorXd scalar_a = inc_beta_ddz(2.5, b, z); - Eigen::VectorXd scalar_b = inc_beta_ddz(a, 3.5, z); - Eigen::VectorXd scalar_z = inc_beta_ddz(a, b, 0.3); - Eigen::VectorXd only_z = inc_beta_ddz(2.5, 3.5, z); - for (int i = 0; i < 3; ++i) { - EXPECT_EQ(inc_beta_ddz(2.5, b[i], z[i]), scalar_a[i]); - EXPECT_EQ(inc_beta_ddz(a[i], 3.5, z[i]), scalar_b[i]); - EXPECT_EQ(inc_beta_ddz(a[i], b[i], 0.3), scalar_z[i]); - EXPECT_EQ(inc_beta_ddz(2.5, 3.5, z[i]), only_z[i]); - } -} diff --git a/test/unit/math/prim/fun/inc_beta_ddz_vectorized_test.cpp b/test/unit/math/prim/fun/inc_beta_ddz_vectorized_test.cpp new file mode 100644 index 00000000000..61c575351da --- /dev/null +++ b/test/unit/math/prim/fun/inc_beta_ddz_vectorized_test.cpp @@ -0,0 +1,40 @@ +#include +#include +#include + +TEST(MathFunctions, inc_beta_ddz_vectorized) { + using stan::math::inc_beta_ddz; + + Eigen::VectorXd a(3); + a << 1.5, 15.0, 0.75; + Eigen::VectorXd b(3); + b << 1.25, 4.0, 30.0; + Eigen::VectorXd z(3); + z << 0.001, 0.5, 0.999; + + Eigen::VectorXd all_eigen = inc_beta_ddz(a, b, z); + ASSERT_EQ(3, all_eigen.size()); + for (int i = 0; i < 3; ++i) { + EXPECT_EQ(inc_beta_ddz(a[i], b[i], z[i]), all_eigen[i]); + } + + std::vector a_std{1.5, 15.0, 0.75}; + std::vector b_std{1.25, 4.0, 30.0}; + std::vector z_std{0.001, 0.5, 0.999}; + std::vector all_std = inc_beta_ddz(a_std, b_std, z_std); + ASSERT_EQ(3, all_std.size()); + for (int i = 0; i < 3; ++i) { + EXPECT_EQ(inc_beta_ddz(a_std[i], b_std[i], z_std[i]), all_std[i]); + } + + Eigen::VectorXd scalar_a = inc_beta_ddz(2.5, b, z); + Eigen::VectorXd scalar_b = inc_beta_ddz(a, 3.5, z); + Eigen::VectorXd scalar_z = inc_beta_ddz(a, b, 0.3); + Eigen::VectorXd only_z = inc_beta_ddz(2.5, 3.5, z); + for (int i = 0; i < 3; ++i) { + EXPECT_EQ(inc_beta_ddz(2.5, b[i], z[i]), scalar_a[i]); + EXPECT_EQ(inc_beta_ddz(a[i], 3.5, z[i]), scalar_b[i]); + EXPECT_EQ(inc_beta_ddz(a[i], b[i], 0.3), scalar_z[i]); + EXPECT_EQ(inc_beta_ddz(2.5, 3.5, z[i]), only_z[i]); + } +} diff --git a/test/unit/math/rev/prob/beta_cdf_log_test.cpp b/test/unit/math/rev/prob/beta_cdf_log_test.cpp new file mode 100644 index 00000000000..ec46553b456 --- /dev/null +++ b/test/unit/math/rev/prob/beta_cdf_log_test.cpp @@ -0,0 +1,46 @@ +#include +#include +#include +#include +#include + +namespace { + +/** + * Shape gradients of beta_lcdf against a high-precision reference. The + * reference is (d/da I) / I and (d/db I) / I from mpmath at 80 digits. + */ +void expect_beta_lcdf_grads(double y, double alpha, double beta, double dalpha, + double dbeta, double rtol) { + using stan::math::var; + var a = alpha; + var b = beta; + var lp = stan::math::beta_lcdf(y, a, b); + std::vector x{a, b}; + std::vector g; + lp.grad(x, g); + EXPECT_TRUE(std::isfinite(g[0])) + << "y=" << y << " a=" << alpha << " b=" << beta; + EXPECT_TRUE(std::isfinite(g[1])) + << "y=" << y << " a=" << alpha << " b=" << beta; + EXPECT_NEAR(dalpha, g[0], rtol * std::fabs(dalpha)) + << "y=" << y << " a=" << alpha << " b=" << beta; + EXPECT_NEAR(dbeta, g[1], rtol * std::fabs(dbeta)) + << "y=" << y << " a=" << alpha << " b=" << beta; +} + +} // namespace + +TEST_F(AgradRev, ProbDistributionsBeta_lcdf_shape_gradients) { + expect_beta_lcdf_grads(0.25, 2, 3, -0.87136898798556226, 0.40853599542368053, + 1e-12); + // large shapes: beta(alpha, beta) underflows + expect_beta_lcdf_grads(0.5, 600, 600, -0.0230409429460951, 0.0230409429460951, + 1e-12); + // y far above the mean: the cdf is 1 and the gradients are 1e-138 + expect_beta_lcdf_grads(0.758, 352, 590, -3.3772717797123241e-138, + 4.5390666074606647e-138, 1e-12); + // integer beta + expect_beta_lcdf_grads(0.5, 10, 1, -0.69314718055994531, 2.3201212683796071, + 1e-12); +} diff --git a/test/unit/math/rev/prob/neg_binomial_cdf_log_test.cpp b/test/unit/math/rev/prob/neg_binomial_cdf_log_test.cpp index 287f3fc5137..eca1e6527d9 100644 --- a/test/unit/math/rev/prob/neg_binomial_cdf_log_test.cpp +++ b/test/unit/math/rev/prob/neg_binomial_cdf_log_test.cpp @@ -31,3 +31,13 @@ TEST_F(AgradRev, ProbDistributionsNegBinomial_derivatives_lcdf) { / (2 * epsilon); EXPECT_FLOAT_EQ(grad_diff2, gradients[1]); } + +TEST_F(AgradRev, ProbDistributionsNegBinomial_lcdf_alpha_lower_tail) { + // lower tail: I = 1e-262 and B(a, b) I = 2e-381 underflows; the reference + // is (d/da I) / I from mpmath at 80 digits + using stan::math::var; + var alpha = 313; + var lp = stan::math::neg_binomial_lcdf(132, alpha, 0.068); + lp.grad(); + EXPECT_NEAR(-2.402731289220054, alpha.adj(), 1e-12 * 2.402731289220054); +}