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]); + } +}