From 560ff14187ab8b37dd19aef92a5aab0ee53c622e Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Mon, 26 Oct 2026 15:02:32 -0400 Subject: [PATCH] Reverse-mode grad_reg_inc_beta with a Hessian-vector product reverse pass Add the rev specialization of grad_reg_inc_beta and use it in rev inc_beta. The prim GRAD_RIB overload now excludes var arguments, which go to the rev overload. Add the rev and mix tests. --- stan/math/prim/fun/grad_reg_inc_beta.hpp | 3 +- stan/math/rev/fun.hpp | 1 + stan/math/rev/fun/grad_reg_inc_beta.hpp | 833 ++++++++++++++++++ stan/math/rev/fun/inc_beta.hpp | 77 +- .../math/mix/fun/grad_reg_inc_beta_test.cpp | 214 ++--- test/unit/math/mix/fun/inc_beta_test.cpp | 23 + test/unit/math/rev/fun/inc_beta_test.cpp | 62 ++ 7 files changed, 1051 insertions(+), 162 deletions(-) create mode 100644 stan/math/rev/fun/grad_reg_inc_beta.hpp create mode 100644 test/unit/math/rev/fun/inc_beta_test.cpp diff --git a/stan/math/prim/fun/grad_reg_inc_beta.hpp b/stan/math/prim/fun/grad_reg_inc_beta.hpp index 92fab1fc056..c34e463f968 100644 --- a/stan/math/prim/fun/grad_reg_inc_beta.hpp +++ b/stan/math/prim/fun/grad_reg_inc_beta.hpp @@ -420,7 +420,8 @@ inline auto inc_beta_continued_fraction_kernel(const T1& a, const T2& b, * the pair (d/da, d/db) for BOTH */ template + typename T4, typename T5, typename T6, typename T7, + require_all_not_var_t* = nullptr> 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 */) { diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp index e5ef938a302..efb00634da1 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -71,6 +71,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/rev/fun/grad_reg_inc_beta.hpp b/stan/math/rev/fun/grad_reg_inc_beta.hpp new file mode 100644 index 00000000000..daf151d46f0 --- /dev/null +++ b/stan/math/rev/fun/grad_reg_inc_beta.hpp @@ -0,0 +1,833 @@ +#ifndef STAN_MATH_REV_FUN_GRAD_REG_INC_BETA_HPP +#define STAN_MATH_REV_FUN_GRAD_REG_INC_BETA_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +namespace internal { + +/** + * Hessian-vector product H v of the regularized incomplete beta function + * I_x(a, b), where H is the Hessian of I_x(a, b) in (a, b), computed from its + * power series. + * + * As in internal::inc_beta_power_series_kernel, 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. + * + * Write grad for the gradient in (a, b) and D_v for the derivative along v. + * Differentiating I = N S twice and applying the result to v gives + * + * H v = N (D_v grad S + grad S D_v log N + grad log N D_v S + * + S (grad log N D_v log N + D_v grad log N)). + * + * When only one shape j is requested, the product needs only the diagonal + * entry + * + * H_jj = N (d2S + 2 dlog_N dS + S (dlog_N^2 + d2log_N)), + * + * where d and d2 are the first and second derivatives in that shape. + * + * The derivatives of S are the sums of those of its terms. The a derivatives + * of a term are closed form. Its b derivatives come from those of P_n, which + * follow from the recurrence for P_n: + * + * dt/da = -t / (a + n), d2t/da2 = 2 t / (a + n)^2, + * dt/db = P'_n / (a + n), d2t/dadb = -(dt/db) / (a + n), + * d2t/db2 = P''_n / (a + n). + * + * For n = 0 the term is 1 / a, with a derivatives -1 / a^2 and 2 / a^3, and + * it does not depend on b. + * + * With psi the digamma function and psi1 the trigamma function, the + * derivatives of log N = a log x - log B(a, b) are + * + * grad log N = (log x - psi(a) + psi(a + b), psi(a + b) - psi(b)), + * Hess log N = [[psi1(a + b) - psi1(a), psi1(a + b)], + * [psi1(a + b), psi1(a + b) - psi1(b)]]. + * + * psi1 is taken from Boost. stan::math::trigamma(double) is accurate only to + * about 1e-9 for arguments below 5, and psi1(a + b) appears alone in H_ab, + * where that error does not cancel. + * + * The loop stops once the terms are below eps / a in magnitude, or after 600 + * terms. For integer b, t_n and its a derivatives are 0 for n >= b, but its b + * derivatives are not. G2 and BOTH therefore also + * wait for the b derivatives of the terms to be small. + * + * The variable names follow two conventions. For G1 and + * G2 only the requested shape varies, and dX and d2X are the + * first and second derivatives of X in that shape. For BOTH, + * X_a and X_b are the partials of X, and Dv_X is the derivative of X along v. + * + * @tparam grads G1 for H_aa v_a, G2 for H_bb v_b, + * BOTH for H v + * @param a first shape + * @param b second shape + * @param x upper limit, in (0, 1) + * @param v_a a component of v, not used for G2 + * @param v_b b component of v, not used for G1 + * @param digamma_a digamma(a), not used for G2 + * @param digamma_b digamma(b), not used for G1 + * @param digamma_ab digamma(a + b) + * @return H_aa v_a for G1, H_bb v_b for G2, + * ((H v)_a, (H v)_b) for BOTH + */ +template +inline auto inc_beta_hvp_power_series(double a, double b, double x, double v_a, + double v_b, double digamma_a, + double digamma_b, double digamma_ab) { + constexpr int max_iterations = 600; + const double tolerance = std::numeric_limits::epsilon() / a; + const double log_x = std::log(x); + const double N = std::exp(a * log_x - lbeta(a, b)); + const double trigamma_ab = boost::math::trigamma(a + b, boost_policy_t<>()); + // n = 0: P_0 = 1 and the term is 1 / a + double product = 1.0; + double S = 1.0 / a; + if constexpr (grads == GRAD_RIB::G1) { + double dS = -S / a; + double d2S = 2.0 * S / (a * a); + for (int n = 1; n <= max_iterations; ++n) { + product *= (n - b) * x / n; + const double inv_apn = 1.0 / (a + n); + const double term = product * inv_apn; + S += term; + dS -= term * inv_apn; + d2S += 2.0 * term * inv_apn * inv_apn; + // for integer b the term and its a derivatives are 0 from n = b on + if (std::fabs(term) <= tolerance) { + break; + } + } + const double dlog_N = log_x - digamma_a + digamma_ab; + const double d2log_N + = trigamma_ab - boost::math::trigamma(a, boost_policy_t<>()); + return v_a * N + * (d2S + 2.0 * dlog_N * dS + S * (dlog_N * dlog_N + d2log_N)); + } else if constexpr (grads == GRAD_RIB::G2) { + double dproduct = 0.0; + double d2product = 0.0; + double dS = 0.0; + double d2S = 0.0; + for (int n = 1; n <= max_iterations; ++n) { + // P_n = f P_{n-1} with f = (n - b) x / n and df/db = -x / n + const double x_div_n = x / n; + const double factor = (n - b) * x_div_n; + d2product = factor * d2product - 2.0 * x_div_n * dproduct; + dproduct = factor * dproduct - x_div_n * product; + product *= factor; + const double inv_apn = 1.0 / (a + n); + const double term = product * inv_apn; + const double dterm = dproduct * inv_apn; + const double d2term = d2product * inv_apn; + S += term; + dS += dterm; + d2S += d2term; + if (std::fabs(term) <= tolerance && std::fabs(dterm) <= tolerance + && std::fabs(d2term) <= tolerance) { + break; + } + } + const double dlog_N = digamma_ab - digamma_b; + const double d2log_N + = trigamma_ab - boost::math::trigamma(b, boost_policy_t<>()); + return v_b * N + * (d2S + 2.0 * dlog_N * dS + S * (dlog_N * dlog_N + d2log_N)); + } else { + double product_b = 0.0; + double product_bb = 0.0; + double S_a = -S / a; + double S_b = 0.0; + double Dv_S_a = 2.0 * v_a * S / (a * a); + double Dv_S_b = 0.0; + for (int n = 1; n <= max_iterations; ++n) { + // P_n = f P_{n-1} with f = (n - b) x / n and df/db = -x / n + const double x_div_n = x / n; + const double factor = (n - b) * x_div_n; + product_bb = factor * product_bb - 2.0 * x_div_n * product_b; + product_b = factor * product_b - x_div_n * product; + product *= factor; + const double inv_apn = 1.0 / (a + n); + const double term = product * inv_apn; + const double term_b = product_b * inv_apn; + const double term_aa = 2.0 * term * inv_apn * inv_apn; + const double term_ab = -term_b * inv_apn; + const double term_bb = product_bb * inv_apn; + S += term; + S_a -= term * inv_apn; + S_b += term_b; + Dv_S_a += v_a * term_aa + v_b * term_ab; + Dv_S_b += v_a * term_ab + v_b * term_bb; + if (std::fabs(term) <= tolerance && std::fabs(term_b) <= tolerance + && std::fabs(term_bb) <= tolerance) { + break; + } + } + const double log_N_a = log_x - digamma_a + digamma_ab; + const double log_N_b = digamma_ab - digamma_b; + const double Dv_log_N = log_N_a * v_a + log_N_b * v_b; + const double Dv_S = S_a * v_a + S_b * v_b; + const double Dv_log_N_a + = (trigamma_ab - boost::math::trigamma(a, boost_policy_t<>())) * v_a + + trigamma_ab * v_b; + const double Dv_log_N_b + = trigamma_ab * v_a + + (trigamma_ab - boost::math::trigamma(b, boost_policy_t<>())) * v_b; + return std::pair{N + * (Dv_S_a + S_a * Dv_log_N + log_N_a * Dv_S + + S * (log_N_a * Dv_log_N + Dv_log_N_a)), + N + * (Dv_S_b + S_b * Dv_log_N + log_N_b * Dv_S + + S * (log_N_b * Dv_log_N + Dv_log_N_b))}; + } +} + +/** + * Hessian-vector product H v of the regularized incomplete beta function + * I_x(a, b), where H is the Hessian of I_x(a, b) in (a, b), computed from its + * continued fraction. + * + * As in internal::inc_beta_continued_fraction_kernel, I_x(a, b) = Q h with + * Q = x^a (1 - x)^b / (a B(a, b)) and h the continued fraction + * https://dlmf.nist.gov/8.17.E22, which is evaluated with the modified Lentz + * method. Write grad for the gradient in (a, b), Hess for the Hessian, and + * D_v for the derivative along v. Differentiating I = Q h twice gives + * + * H = I (grad log I grad log I^T + Hess log Q + Hess log h), + * + * with grad log I = grad log Q + grad log h. With psi the digamma function + * and psi1 the trigamma function, + * + * grad log Q = (log x - 1 / a - psi(a) + psi(a + b), + * log(1 - x) - psi(b) + psi(a + b)), + * Hess log Q = [[psi1(a + b) - psi1(a) + 1 / a^2, psi1(a + b)], + * [psi1(a + b), psi1(a + b) - psi1(b)]], + * + * with psi1 taken from Boost, as in the power series. + * + * Each Lentz step multiplies h by C D, so log h is the sum of log C and + * log D over the steps, and its derivatives are the sums of theirs. With d + * the partial numerator of the step, the step sets C' = 1 + d / C and + * D' = 1 / (1 + d D). Their log derivatives are + * + * grad log C' = (grad d - d grad log C) / (C + d), + * grad log D' = -(grad d + d grad log D) D D', + * + * using C C' = C + d. Differentiating these once more along v gives + * + * D_v grad log C' = (Hess d v - D_v d grad log C - d D_v grad log C) + * / (C + d) - grad log C' D_v log(C C'), + * D_v log(C C') = (C D_v log C + D_v d) / (C + d), + * D_v grad log D' = -(Hess d v + D_v d grad log D + d D_v grad log D) D D' + * + grad log D' D_v log(D D'), + * D_v log(D D') = D_v log D - D D' (D_v d + d D_v log D). + * + * The code applies each of the two second-derivative formulas in three parts. + * It first computes the part that uses only values from before the step. It + * then overwrites grad log C (or grad log D) with its new value, and last it + * adds the final term, which uses that new value. + * + * G1 and G2 take v = (1, 0) and v = (0, 1), so D_v + * is the derivative in the requested shape, and the result is scaled by v_a + * or v_b at the end. + * + * The partial numerators https://dlmf.nist.gov/8.17.E23 are d_n = num / den, + * with den = (a + n - 1)(a + n) and + * + * n = 2m: num = m (b - m) x, num_b = m x, + * n = 2m + 1: num = -(a + m)(a + b + m) x, num_a = -(2 (a + m) + b) x, + * num_b = -(a + m) x, num_aa = -2 x, num_ab = -x. + * + * Every other derivative of num is 0, so d_bb = 0. The derivatives of d_n + * follow from the quotient rule, with den_a = 2 (a + n) - 1 and den_aa = 2. + * n = 1 is the odd case with m = 0. The first step holds C at 1, so + * D = 1 / (1 + d_1). + * + * The test for stopping runs after each odd n. The loop stops once C D is + * within 3 eps of 1 and the second derivatives of log h have changed by at + * most a relative 3 eps since the previous test, or after n = 601. + * + * The variable names follow the same two conventions as in the power series. + * For G1 and G2, dX and d2X are the first and + * second derivatives of X in the requested shape. For BOTH, + * X_a and X_b are the partials of X, and Dv_X is the derivative of X along v. + * + * @tparam grads G1 for H_aa v_a, G2 for H_bb v_b, + * BOTH for H v + * @param a first shape + * @param b second shape + * @param x upper limit, in (0, 1) + * @param v_a a component of v, not used for G2 + * @param v_b b component of v, not used for G1 + * @param digamma_a digamma(a), not used for G2 + * @param digamma_b digamma(b), not used for G1 + * @param digamma_ab digamma(a + b) + * @return H_aa v_a for G1, H_bb v_b for G2, + * ((H v)_a, (H v)_b) for BOTH + */ +template +inline auto inc_beta_hvp_continued_fraction(double a, double b, double x, + double v_a, double v_b, + double digamma_a, double digamma_b, + double digamma_ab) { + constexpr int max_iterations = 300; + const double tolerance = 3.0 * std::numeric_limits::epsilon(); + const double small = std::sqrt(std::numeric_limits::min()); + const double log_x = std::log(x); + const double log1m_x = log1p(-x); + const double Q + = std::exp(a * log_x + b * log1m_x - std::log(a) - lbeta(a, b)); + // boost's trigamma, as in the power series + const double trigamma_ab = boost::math::trigamma(a + b, boost_policy_t<>()); + // first step, n = 1, with C held at 1 + const double ap1 = a + 1.0; + const double d_1 = -(a + b) * x / ap1; + double D = 1.0 + d_1; + if (std::fabs(D) < small) { + D = small; + } + D = 1.0 / D; + double C = 1.0; + double h = D; + if constexpr (grads == GRAD_RIB::G1) { + const double dd_1 = -(x + d_1) / ap1; + const double d2d_1 = 2.0 * (x + d_1) / (ap1 * ap1); + double dlog_C = 0.0; + double d2log_C = 0.0; + double dlog_D = -dd_1 * D; + double d2log_D = -(d2d_1 + dlog_D * dd_1) * D; + double dlog_h = dlog_D; + double d2log_h = d2log_D; + double d2log_h_prev = d2log_h; + for (int n = 2; n <= 2 * max_iterations + 1; ++n) { + const int m = n / 2; + const bool odd = n % 2 == 1; + const double a_plus_n = a + n; + const double inv_den = 1.0 / (a_plus_n * (a_plus_n - 1.0)); + const double dden = 2.0 * a_plus_n - 1.0; + double num; + double dnum = 0.0; + double d2num = 0.0; + if (odd) { + num = -(a + m) * (a + b + m) * x; + dnum = -(2.0 * (a + m) + b) * x; + d2num = -2.0 * x; + } else { + num = m * (b - m) * x; + } + const double d_n = num * inv_den; + const double dd_n = (dnum - d_n * dden) * inv_den; + const double d2d_n = (d2num - 2.0 * dd_n * dden - 2.0 * d_n) * inv_den; + + // one step of the modified Lentz method + // C_new = 1 + d_n / C, with C C_new = C + d_n + double C_new = 1.0 + d_n / C; + double C_C_new = C + d_n; + if (std::fabs(C_new) < small) { + C_new = small; + C_C_new = C * small; + } + const double inv_C_C_new = 1.0 / C_C_new; + const double dlog_C_C_new = (C * dlog_C + dd_n) * inv_C_C_new; + d2log_C = (d2d_n - dd_n * dlog_C - d_n * d2log_C) * inv_C_C_new; + dlog_C = (dd_n - d_n * dlog_C) * inv_C_C_new; + d2log_C -= dlog_C * dlog_C_C_new; + + // D_new = 1 / (1 + d_n D) + double D_new = 1.0 + d_n * D; + if (std::fabs(D_new) < small) { + D_new = small; + } + D_new = 1.0 / D_new; + const double D_D_new = D * D_new; + const double dlog_D_D_new = dlog_D - D_D_new * (dd_n + d_n * dlog_D); + d2log_D = -(d2d_n + dd_n * dlog_D + d_n * d2log_D) * D_D_new; + dlog_D = -(dd_n + d_n * dlog_D) * D_D_new; + d2log_D += dlog_D * dlog_D_D_new; + + dlog_h += dlog_C + dlog_D; + d2log_h += d2log_C + d2log_D; + C = C_new; + D = D_new; + const double delta = C_new * D_new; + h *= delta; + if (odd) { + if (std::fabs(delta - 1.0) < tolerance + && std::fabs(d2log_h - d2log_h_prev) + <= tolerance * std::fabs(d2log_h)) { + break; + } + d2log_h_prev = d2log_h; + } + } + const double dlog_I = log_x - 1.0 / a + digamma_ab - digamma_a + dlog_h; + const double d2log_Q = trigamma_ab + - boost::math::trigamma(a, boost_policy_t<>()) + + 1.0 / (a * a); + return v_a * h * Q * (dlog_I * dlog_I + d2log_Q + d2log_h); + } else if constexpr (grads == GRAD_RIB::G2) { + // d_n is linear in b, so d2d_n = 0 + const double dd_1 = -x / ap1; + double dlog_C = 0.0; + double d2log_C = 0.0; + double dlog_D = -dd_1 * D; + double d2log_D = -dlog_D * dd_1 * D; + double dlog_h = dlog_D; + double d2log_h = d2log_D; + double d2log_h_prev = d2log_h; + for (int n = 2; n <= 2 * max_iterations + 1; ++n) { + const int m = n / 2; + const bool odd = n % 2 == 1; + const double a_plus_n = a + n; + const double inv_den = 1.0 / (a_plus_n * (a_plus_n - 1.0)); + double dnum; + double num; + if (odd) { + dnum = -(a + m) * x; + num = dnum * (a + b + m); + } else { + dnum = m * x; + num = dnum * (b - m); + } + const double d_n = num * inv_den; + const double dd_n = dnum * inv_den; + + // one step of the modified Lentz method, as for G1 with d2d_n = 0. + // C_new = 1 + d_n / C, with C C_new = C + d_n + double C_new = 1.0 + d_n / C; + double C_C_new = C + d_n; + if (std::fabs(C_new) < small) { + C_new = small; + C_C_new = C * small; + } + const double inv_C_C_new = 1.0 / C_C_new; + const double dlog_C_C_new = (C * dlog_C + dd_n) * inv_C_C_new; + d2log_C = -(dd_n * dlog_C + d_n * d2log_C) * inv_C_C_new; + dlog_C = (dd_n - d_n * dlog_C) * inv_C_C_new; + d2log_C -= dlog_C * dlog_C_C_new; + + // D_new = 1 / (1 + d_n D) + double D_new = 1.0 + d_n * D; + if (std::fabs(D_new) < small) { + D_new = small; + } + D_new = 1.0 / D_new; + const double D_D_new = D * D_new; + const double dlog_D_D_new = dlog_D - D_D_new * (dd_n + d_n * dlog_D); + d2log_D = -(dd_n * dlog_D + d_n * d2log_D) * D_D_new; + dlog_D = -(dd_n + d_n * dlog_D) * D_D_new; + d2log_D += dlog_D * dlog_D_D_new; + + dlog_h += dlog_C + dlog_D; + d2log_h += d2log_C + d2log_D; + C = C_new; + D = D_new; + const double delta = C_new * D_new; + h *= delta; + if (odd) { + if (std::fabs(delta - 1.0) < tolerance + && std::fabs(d2log_h - d2log_h_prev) + <= tolerance * std::fabs(d2log_h)) { + break; + } + d2log_h_prev = d2log_h; + } + } + const double dlog_I = log1m_x + digamma_ab - digamma_b + dlog_h; + const double d2log_Q + = trigamma_ab - boost::math::trigamma(b, boost_policy_t<>()); + return v_b * h * Q * (dlog_I * dlog_I + d2log_Q + d2log_h); + } else { + double log_C_a = 0.0; + double log_C_b = 0.0; + double Dv_log_C_a = 0.0; + double Dv_log_C_b = 0.0; + const double d_1_a = -(x + d_1) / ap1; + const double d_1_b = -x / ap1; + double log_D_a = -d_1_a * D; + double log_D_b = -d_1_b * D; + // D_v log D = -D D_v d_1 + const double d_1_aa = 2.0 * (x + d_1) / (ap1 * ap1); + const double d_1_ab = x / (ap1 * ap1); + const double Dv_d_1 = d_1_a * v_a + d_1_b * v_b; + double Dv_log_D_a + = -(d_1_aa * v_a + d_1_ab * v_b) * D - log_D_a * D * Dv_d_1; + double Dv_log_D_b = -d_1_ab * v_a * D - log_D_b * D * Dv_d_1; + double log_h_a = log_D_a; + double log_h_b = log_D_b; + double Dv_log_h_a = Dv_log_D_a; + double Dv_log_h_b = Dv_log_D_b; + double Dv_log_h_a_prev = Dv_log_h_a; + double Dv_log_h_b_prev = Dv_log_h_b; + for (int n = 2; n <= 2 * max_iterations + 1; ++n) { + const int m = n / 2; + const bool odd = n % 2 == 1; + const double a_plus_n = a + n; + const double inv_den = 1.0 / (a_plus_n * (a_plus_n - 1.0)); + const double den_a = 2.0 * a_plus_n - 1.0; + double num_b; + double num; + double num_a = 0.0; + double num_aa = 0.0; + double num_ab = 0.0; + if (odd) { + num_b = -(a + m) * x; + num = num_b * (a + b + m); + num_a = -(2.0 * (a + m) + b) * x; + num_aa = -2.0 * x; + num_ab = -x; + } else { + num_b = m * x; + num = num_b * (b - m); + } + const double d_n = num * inv_den; + const double d_n_a = (num_a - d_n * den_a) * inv_den; + const double d_n_b = num_b * inv_den; + const double d_n_aa + = (num_aa - 2.0 * d_n_a * den_a - 2.0 * d_n) * inv_den; + const double d_n_ab = (num_ab - d_n_b * den_a) * inv_den; + const double Dv_d_n = d_n_a * v_a + d_n_b * v_b; + const double Dv_d_n_a = d_n_aa * v_a + d_n_ab * v_b; + const double Dv_d_n_b = d_n_ab * v_a; + + // one step of the modified Lentz method + // C_new = 1 + d_n / C, with C C_new = C + d_n + double C_new = 1.0 + d_n / C; + double C_C_new = C + d_n; + if (std::fabs(C_new) < small) { + C_new = small; + C_C_new = C * small; + } + const double inv_C_C_new = 1.0 / C_C_new; + const double Dv_log_C = log_C_a * v_a + log_C_b * v_b; + const double Dv_log_C_C_new = (C * Dv_log_C + Dv_d_n) * inv_C_C_new; + Dv_log_C_a + = (Dv_d_n_a - Dv_d_n * log_C_a - d_n * Dv_log_C_a) * inv_C_C_new; + log_C_a = (d_n_a - d_n * log_C_a) * inv_C_C_new; + Dv_log_C_a -= log_C_a * Dv_log_C_C_new; + Dv_log_C_b + = (Dv_d_n_b - Dv_d_n * log_C_b - d_n * Dv_log_C_b) * inv_C_C_new; + log_C_b = (d_n_b - d_n * log_C_b) * inv_C_C_new; + Dv_log_C_b -= log_C_b * Dv_log_C_C_new; + + // D_new = 1 / (1 + d_n D) + double D_new = 1.0 + d_n * D; + if (std::fabs(D_new) < small) { + D_new = small; + } + D_new = 1.0 / D_new; + const double D_D_new = D * D_new; + const double Dv_log_D = log_D_a * v_a + log_D_b * v_b; + const double Dv_log_D_D_new + = Dv_log_D - D_D_new * (Dv_d_n + d_n * Dv_log_D); + Dv_log_D_a = -(Dv_d_n_a + Dv_d_n * log_D_a + d_n * Dv_log_D_a) * D_D_new; + log_D_a = -(d_n_a + d_n * log_D_a) * D_D_new; + Dv_log_D_a += log_D_a * Dv_log_D_D_new; + Dv_log_D_b = -(Dv_d_n_b + Dv_d_n * log_D_b + d_n * Dv_log_D_b) * D_D_new; + log_D_b = -(d_n_b + d_n * log_D_b) * D_D_new; + Dv_log_D_b += log_D_b * Dv_log_D_D_new; + + log_h_a += log_C_a + log_D_a; + log_h_b += log_C_b + log_D_b; + Dv_log_h_a += Dv_log_C_a + Dv_log_D_a; + Dv_log_h_b += Dv_log_C_b + Dv_log_D_b; + C = C_new; + D = D_new; + const double delta = C_new * D_new; + h *= delta; + if (odd) { + if (std::fabs(delta - 1.0) < tolerance + && std::fabs(Dv_log_h_a - Dv_log_h_a_prev) + <= tolerance * std::fabs(Dv_log_h_a) + && std::fabs(Dv_log_h_b - Dv_log_h_b_prev) + <= tolerance * std::fabs(Dv_log_h_b)) { + break; + } + Dv_log_h_a_prev = Dv_log_h_a; + Dv_log_h_b_prev = Dv_log_h_b; + } + } + const double I = h * Q; + const double log_I_a = log_x - 1.0 / a + digamma_ab - digamma_a + log_h_a; + const double log_I_b = log1m_x + digamma_ab - digamma_b + log_h_b; + const double Dv_log_I = log_I_a * v_a + log_I_b * v_b; + // the Hessian of log Q is that of -log B(a, b) plus diag(1 / a^2, 0) + const double Dv_log_Q_a + = (trigamma_ab - boost::math::trigamma(a, boost_policy_t<>()) + + 1.0 / (a * a)) + * v_a + + trigamma_ab * v_b; + const double Dv_log_Q_b + = trigamma_ab * v_a + + (trigamma_ab - boost::math::trigamma(b, boost_policy_t<>())) * v_b; + return std::pair{I * (log_I_a * Dv_log_I + Dv_log_Q_a + Dv_log_h_a), + I * (log_I_b * Dv_log_I + Dv_log_Q_b + Dv_log_h_b)}; + } +} + +/** + * Hessian-vector product H v of the regularized incomplete beta function + * I_z(a, b), where H is the Hessian of I_z(a, b) in (a, b). When v holds the + * adjoints of (d/da I_z(a, b), d/db I_z(a, b)), H v is what the reverse pass + * adds to the adjoints of a and b. + * + * The power series or the continued fraction, and whether to apply the + * symmetry relation, are chosen by the same rule as in grad_reg_inc_beta. + * Under the symmetry relation I_z(a, b) = 1 - I_{1-z}(b, a), let H' be the + * Hessian of I_{1-z}(b, a) in its own arguments (b, a), and let P swap the + * two components of a vector. Then H = -P H' P, so H v = -P H' (P v). In + * particular H_aa = -H'_bb and H_bb = -H'_aa. + * + * @tparam grads G1 for H_aa v_a, G2 for H_bb v_b, + * BOTH for H v + * @param a first shape + * @param b second shape + * @param z upper limit, in [0, 1] + * @param v_a a component of v, not used for G2 + * @param v_b b component of v, not used for G1 + * @param digamma_a digamma(a), not used for G2 + * @param digamma_b digamma(b), not used for G1 + * @param digamma_ab digamma(a + b) + * @return H_aa v_a for G1, H_bb v_b for G2, + * ((H v)_a, (H v)_b) for BOTH. The result is 0 at z = 0 and + * z = 1. + */ +template +inline std::conditional_t, + double> +inc_beta_hvp(double a, double b, double z, double v_a, double v_b, + double digamma_a, double digamma_b, double digamma_ab) { + if (z == 0.0 || z == 1.0) { + return {}; + } + // the symmetry relation swaps which shape is requested + constexpr GRAD_RIB swapped_grads + = grads == GRAD_RIB::G1 + ? GRAD_RIB::G2 + : grads == GRAD_RIB::G2 ? GRAD_RIB::G1 : GRAD_RIB::BOTH; + const double mean = a / (a + b); + // the same branch rule as grad_reg_inc_beta + const bool use_power_series + = (z < mean && b * z <= 1.0 && z <= 0.5) + || (z >= mean && a * (1.0 - z) <= 1.0 && z >= 0.5); + if (use_power_series) { + if (z >= mean) { + const auto hv = inc_beta_hvp_power_series( + b, a, 1.0 - z, v_b, v_a, digamma_b, digamma_a, digamma_ab); + if constexpr (grads != GRAD_RIB::BOTH) { + return -hv; + } else { + return {-hv.second, -hv.first}; + } + } else { + return inc_beta_hvp_power_series(a, b, z, v_a, v_b, digamma_a, + digamma_b, digamma_ab); + } + } else { + // the same threshold as grad_reg_inc_beta, so both take the same branch + if (z >= (a + 1.0) / (a + b + 2.0)) { + const auto hv = inc_beta_hvp_continued_fraction( + b, a, 1.0 - z, v_b, v_a, digamma_b, digamma_a, digamma_ab); + if constexpr (grads != GRAD_RIB::BOTH) { + return -hv; + } else { + return {-hv.second, -hv.first}; + } + } else { + return inc_beta_hvp_continued_fraction( + a, b, z, v_a, v_b, digamma_a, digamma_b, digamma_ab); + } + } +} + +} // 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, for var arguments. + * + * The values come from the prim overload, evaluated in double. The reverse + * pass does not differentiate through that series or continued fraction. + * Instead, with v the adjoints of the outputs, it adds H v to the adjoints of + * a and b, where H is the Hessian of I_z(a, b) in (a, b), computed directly + * by internal::inc_beta_hvp. To the adjoint of z it adds v times the + * closed-form mixed derivatives + * + * d/dz d/da I_z(a, b) = inc_beta_ddz(a, b, z) (log z - psi(a) + psi(a+b)), + * d/dz d/db I_z(a, b) = inc_beta_ddz(a, b, z) (log1p(-z) - psi(b) + + * psi(a+b)), + * + * where psi is the digamma function. + * + * The digamma arguments are passed only so that the values are cheap to + * compute, and they must equal digamma(a), digamma(b), and digamma(a + b). + * They get no adjoint. The a and b adjoints above are total derivatives, so + * they already include any dependence on a and b that runs through the + * digamma values. + * + * A one-sided request still needs the cross term of H when both shapes are + * vars. For G1 the adjoint of b is H_ab v_a, which the + * BOTH product gives with v_b = 0. That product needs + * digamma(b), which a G1 call does not receive, so the reverse + * pass computes it when b is a var. When only a is a var, H_aa v_a does not + * need digamma(b). G2 is the same with a and b swapped. + * + * For BOTH there are two outputs, and one reverse pass callback + * reads both of their adjoints. It runs after every use of the outputs has + * added into those adjoints. + * + * Each argument is a var or arithmetic. Calls that mix in + * fvar go to the prim overload. + * + * @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 + * @return d/da I_z(a, b) for G1, d/db for G2, and + * the pair (d/da, d/db) for BOTH + */ +template * = nullptr> +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 */) { + const double digamma_a = value_of(digammaA); + const double digamma_b = value_of(digammaB); + const double digamma_ab = value_of(digammaSum); + if constexpr (grads_requested == GRAD_RIB::G1) { + return make_callback_var( + grad_reg_inc_beta(value_of(a), value_of(b), + value_of(z), digamma_a, digamma_b, + digamma_ab, 0.0), + [a, b, z, digamma_a, digamma_ab](auto& vi) { + const double a_val = value_of(a); + const double b_val = value_of(b); + const double z_val = value_of(z); + if constexpr (is_var_v && is_var_v) { + const auto [hv_a, hv_b] = internal::inc_beta_hvp( + a_val, b_val, z_val, vi.adj(), 0.0, digamma_a, digamma(b_val), + digamma_ab); + a.adj() += hv_a; + b.adj() += hv_b; + } else if constexpr (is_var_v) { + // H_aa v_a does not use digamma(b) + a.adj() += internal::inc_beta_hvp( + a_val, b_val, z_val, vi.adj(), 0.0, digamma_a, 0.0, digamma_ab); + } else if constexpr (is_var_v) { + b.adj() += internal::inc_beta_hvp( + a_val, b_val, z_val, vi.adj(), 0.0, digamma_a, + digamma(b_val), digamma_ab) + .second; + } + if constexpr (is_var_v) { + if (z_val > 0.0 && z_val < 1.0) { + const double log_term = std::log(z_val) - digamma_a + digamma_ab; + z.adj() + += vi.adj() * inc_beta_ddz(a_val, b_val, z_val) * log_term; + } + } + }); + } else if constexpr (grads_requested == GRAD_RIB::G2) { + return make_callback_var( + grad_reg_inc_beta(value_of(a), value_of(b), + value_of(z), digamma_a, digamma_b, + digamma_ab, 0.0), + [a, b, z, digamma_b, digamma_ab](auto& vi) { + const double a_val = value_of(a); + const double b_val = value_of(b); + const double z_val = value_of(z); + if constexpr (is_var_v && is_var_v) { + const auto [hv_a, hv_b] = internal::inc_beta_hvp( + a_val, b_val, z_val, 0.0, vi.adj(), digamma(a_val), digamma_b, + digamma_ab); + a.adj() += hv_a; + b.adj() += hv_b; + } else if constexpr (is_var_v) { + a.adj() += internal::inc_beta_hvp( + a_val, b_val, z_val, 0.0, vi.adj(), digamma(a_val), + digamma_b, digamma_ab) + .first; + } else if constexpr (is_var_v) { + // H_bb v_b does not use digamma(a) + b.adj() += internal::inc_beta_hvp( + a_val, b_val, z_val, 0.0, vi.adj(), 0.0, digamma_b, digamma_ab); + } + if constexpr (is_var_v) { + if (z_val > 0.0 && z_val < 1.0) { + const double log_term = log1p(-z_val) - digamma_b + digamma_ab; + z.adj() + += vi.adj() * inc_beta_ddz(a_val, b_val, z_val) * log_term; + } + } + }); + } else { // both + const auto [g1_val, g2_val] = grad_reg_inc_beta( + value_of(a), value_of(b), value_of(z), digamma_a, digamma_b, digamma_ab, + 0.0); + var g1 = g1_val; + var g2 = g2_val; + reverse_pass_callback([a, b, z, g1, g2, digamma_a, digamma_b, + digamma_ab]() { + const double a_val = value_of(a); + const double b_val = value_of(b); + const double z_val = value_of(z); + const double v_a = g1.adj(); + const double v_b = g2.adj(); + if constexpr (is_var_v || is_var_v) { + const auto [hv_a, hv_b] = internal::inc_beta_hvp( + a_val, b_val, z_val, v_a, v_b, digamma_a, digamma_b, digamma_ab); + if constexpr (is_var_v) { + a.adj() += hv_a; + } + if constexpr (is_var_v) { + b.adj() += hv_b; + } + } + if constexpr (is_var_v) { + if (z_val > 0.0 && z_val < 1.0) { + z.adj() += inc_beta_ddz(a_val, b_val, z_val) + * (v_a * (std::log(z_val) - digamma_a + digamma_ab) + + v_b * (log1p(-z_val) - digamma_b + digamma_ab)); + } + } + }); + return std::pair{g1, g2}; + } +} + +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/rev/fun/inc_beta.hpp b/stan/math/rev/fun/inc_beta.hpp index be8a65b7b6a..037f3f07433 100644 --- a/stan/math/rev/fun/inc_beta.hpp +++ b/stan/math/rev/fun/inc_beta.hpp @@ -3,41 +3,60 @@ #include #include -#include -#include -#include +#include #include #include -#include +#include +#include namespace stan { namespace math { -namespace internal { - -class inc_beta_vvv_vari : public op_vvv_vari { - public: - inc_beta_vvv_vari(vari* avi, vari* bvi, vari* cvi) - : op_vvv_vari(inc_beta(avi->val_, bvi->val_, cvi->val_), avi, bvi, cvi) {} - void chain() { - double d_a; - double d_b; - const double beta_ab = beta(avi_->val_, bvi_->val_); - grad_reg_inc_beta(d_a, d_b, avi_->val_, bvi_->val_, cvi_->val_, - digamma(avi_->val_), digamma(bvi_->val_), - digamma(avi_->val_ + bvi_->val_), beta_ab); - - avi_->adj_ += adj_ * d_a; - bvi_->adj_ += adj_ * d_b; - cvi_->adj_ += adj_ * std::pow(1 - cvi_->val_, bvi_->val_ - 1) - * std::pow(cvi_->val_, avi_->val_ - 1) / beta_ab; - } -}; - -} // namespace internal - -inline var inc_beta(const var& a, const var& b, const var& c) { - return var(new internal::inc_beta_vvv_vari(a.vi_, b.vi_, c.vi_)); +/** + * The regularized incomplete beta function I_z(a, b) for var + * arguments. Only the partials of the var arguments are + * computed in the reverse pass. Each argument is a var or + * arithmetic. Mixing with fvar is not supported. + * + * @tparam T1 type of a + * @tparam T2 type of b + * @tparam T3 type of z + * @param a first shape, positive + * @param b second shape, positive + * @param z upper limit, in [0, 1] + * @return I_z(a, b) + */ +template * = nullptr, + require_all_var_or_arithmetic_t* = nullptr, + require_any_var_t* = nullptr> +inline var inc_beta(const T1& a, const T2& b, const T3& z) { + return make_callback_var( + inc_beta(value_of(a), value_of(b), value_of(z)), [a, b, z](auto& vi) { + const double a_val = value_of(a); + const double b_val = value_of(b); + const double z_val = value_of(z); + if constexpr (is_autodiff_v && is_autodiff_v) { + const auto [d_a, d_b] = grad_reg_inc_beta( + a_val, b_val, z_val, digamma(a_val), digamma(b_val), + digamma(a_val + b_val), 0.0); + a.adj() += vi.adj() * d_a; + b.adj() += vi.adj() * d_b; + } else if constexpr (is_autodiff_v) { + a.adj() += vi.adj() + * grad_reg_inc_beta( + a_val, b_val, z_val, digamma(a_val), 0.0, + digamma(a_val + b_val), 0.0); + } else if constexpr (is_autodiff_v) { + b.adj() += vi.adj() + * grad_reg_inc_beta( + a_val, b_val, z_val, 0.0, digamma(b_val), + digamma(a_val + b_val), 0.0); + } + if constexpr (is_autodiff_v) { + z.adj() += vi.adj() * inc_beta_ddz(a_val, b_val, z_val); + } + }); } } // namespace math diff --git a/test/unit/math/mix/fun/grad_reg_inc_beta_test.cpp b/test/unit/math/mix/fun/grad_reg_inc_beta_test.cpp index e6613d43ea1..7831ab106ee 100644 --- a/test/unit/math/mix/fun/grad_reg_inc_beta_test.cpp +++ b/test/unit/math/mix/fun/grad_reg_inc_beta_test.cpp @@ -1,145 +1,95 @@ -#include -#include -#include +#include +#include -TEST(ProbInternalMath, grad_reg_inc_beta_fv) { - using stan::math::digamma; - using stan::math::exp; - using stan::math::fvar; - using stan::math::lbeta; - using stan::math::var; +namespace { - fvar a = 1.0; - fvar b = 1.0; - fvar g = 0.4; - a.d_ = 1.0; - b.d_ = 1.0; - g.d_ = 1.0; - fvar dig_a = digamma(a); - fvar dig_b = digamma(b); - fvar dig_sum = digamma(a + b); - fvar beta_ab = exp(lbeta(a, b)); - fvar g_a; - fvar g_b; - - stan::math::grad_reg_inc_beta(g_a, g_b, a, b, g, dig_a, dig_b, dig_sum, - beta_ab); - EXPECT_FLOAT_EQ(-0.36651629442883944183907601651838247842001142107486495485, - g_a.val_.val()); - EXPECT_NEAR(0.306495375042422864944011633197968575202046200428315551199, - g_b.val_.val(), 1e-6); +/** + * Runs expect_ad at points in each region of grad_reg_inc_beta: + * (1, 1, 0.4), (2, 3, 0.25), (0.5, 0.5, 0.3): power series + * (10, 1, 0.5): power series, integer b + * (2, 3, 0.6): power series through the symmetry relation + * (20, 20, 0.3), (5, 30, 0.1): continued fraction + * (20, 20, 0.7), (5, 30, 0.3): continued fraction through the symmetry + * relation + * (10, 1, 0.7), (4, 4, 0.5): continued fraction, integer b, where the + * partial numerator d_{2b} is 0 but its b derivative is not + */ +template +void expect_ad_grad_reg_inc_beta(const F& f) { + stan::test::expect_ad(f, 1.0, 1.0, 0.4); + stan::test::expect_ad(f, 2.0, 3.0, 0.25); + stan::test::expect_ad(f, 0.5, 0.5, 0.3); + stan::test::expect_ad(f, 10.0, 1.0, 0.5); + stan::test::expect_ad(f, 2.0, 3.0, 0.6); + stan::test::expect_ad(f, 20.0, 20.0, 0.3); + stan::test::expect_ad(f, 5.0, 30.0, 0.1); + stan::test::expect_ad(f, 20.0, 20.0, 0.7); + stan::test::expect_ad(f, 5.0, 30.0, 0.3); + stan::test::expect_ad(f, 10.0, 1.0, 0.7); + stan::test::expect_ad(f, 4.0, 4.0, 0.5); } -TEST(ProbInternalMath, grad_reg_inc_beta_fv_1stDeriv1) { - using stan::math::digamma; - using stan::math::exp; - using stan::math::fvar; - using stan::math::lbeta; - using stan::math::var; - - fvar a = 1.0; - fvar b = 1.0; - fvar g = 0.4; - a.d_ = 1.0; - b.d_ = 0.0; - g.d_ = 0.0; - fvar dig_a = digamma(a); - fvar dig_b = digamma(b); - fvar dig_sum = digamma(a + b); - fvar beta_ab = exp(lbeta(a, b)); - fvar g_a; - fvar g_b; - stan::math::grad_reg_inc_beta(g_a, g_b, a, b, g, dig_a, dig_b, dig_sum, - beta_ab); +} // namespace - std::vector y1{a.val_}; - std::vector grad1; - g_a.val_.grad(y1, grad1); - EXPECT_FLOAT_EQ(0.33583548212738989400284958327902414335945097838423129, - grad1[0]); +TEST(mathMixScalFun, grad_reg_inc_beta_g1) { + auto f = [](const auto& a, const auto& b, const auto& z) { + using stan::math::digamma; + return stan::math::grad_reg_inc_beta( + a, b, z, digamma(a), 0.0, digamma(a + b), 0.0); + }; + expect_ad_grad_reg_inc_beta(f); } -TEST(ProbInternalMath, grad_reg_inc_beta_fv_1stDeriv2) { - using stan::math::digamma; - using stan::math::exp; - using stan::math::fvar; - using stan::math::lbeta; - using stan::math::var; - - fvar a = 1.0; - fvar b = 1.0; - fvar g = 0.4; - a.d_ = 0.0; - b.d_ = 1.0; - g.d_ = 0.0; - fvar dig_a = digamma(a); - fvar dig_b = digamma(b); - fvar dig_sum = digamma(a + b); - fvar beta_ab = exp(lbeta(a, b)); - fvar g_a; - fvar g_b; - - stan::math::grad_reg_inc_beta(g_a, g_b, a, b, g, dig_a, dig_b, dig_sum, - beta_ab); - std::vector y1{b.val_}; - std::vector grad1; - g_b.val_.grad(y1, grad1); - EXPECT_NEAR(-0.156565690737548079304827886, grad1[0], 1e-6); +TEST(mathMixScalFun, grad_reg_inc_beta_g2) { + auto f = [](const auto& a, const auto& b, const auto& z) { + using stan::math::digamma; + return stan::math::grad_reg_inc_beta( + a, b, z, 0.0, digamma(b), digamma(a + b), 0.0); + }; + expect_ad_grad_reg_inc_beta(f); } -TEST(ProbInternalMath, grad_reg_inc_beta_fv_2ndDeriv1) { - using stan::math::digamma; - using stan::math::exp; - using stan::math::fvar; - using stan::math::lbeta; - using stan::math::var; - fvar a = 1.0; - fvar b = 1.0; - fvar g = 0.4; - a.d_ = 1.0; - b.d_ = 0.0; - g.d_ = 0.0; - fvar dig_a = digamma(a); - fvar dig_b = digamma(b); - fvar dig_sum = digamma(a + b); - fvar beta_ab = exp(lbeta(a, b)); - fvar g_a; - fvar g_b; - - stan::math::grad_reg_inc_beta(g_a, g_b, a, b, g, dig_a, dig_b, dig_sum, - beta_ab); - - std::vector y1{a.val_}; - std::vector grad1; - g_a.d_.grad(y1, grad1); - EXPECT_FLOAT_EQ(-0.30772293970781581317390510390046098438962772318921, - grad1[0]); +TEST(mathMixScalFun, grad_reg_inc_beta_both) { + auto f = [](const auto& a, const auto& b, const auto& z) { + using stan::math::digamma; + auto grads = stan::math::grad_reg_inc_beta( + a, b, z, digamma(a), digamma(b), digamma(a + b), 0.0); + Eigen::Matrix ret(2); + ret << grads.first, grads.second; + return ret; + }; + expect_ad_grad_reg_inc_beta(f); } -TEST(ProbInternalMath, grad_reg_inc_beta_fv_2ndDeriv2) { - using stan::math::digamma; - using stan::math::exp; - using stan::math::fvar; - using stan::math::lbeta; - using stan::math::var; - fvar a = 1.0; - fvar b = 1.0; - fvar g = 0.4; - a.d_ = 0.0; - b.d_ = 1.0; - g.d_ = 0.0; - fvar dig_a = digamma(a); - fvar dig_b = digamma(b); - fvar dig_sum = digamma(a + b); - fvar beta_ab = exp(lbeta(a, b)); - fvar g_a; - fvar g_b; - - stan::math::grad_reg_inc_beta(g_a, g_b, a, b, g, dig_a, dig_b, dig_sum, - beta_ab); +// the student_t_cdf call: a = nu / 2, b = 1 / 2 fixed, z = 1 - r +TEST(mathMixScalFun, grad_reg_inc_beta_student) { + auto f = [](const auto& a, const auto& /* b */, const auto& z) { + using stan::math::digamma; + return stan::math::grad_reg_inc_beta( + 0.5 * a, 0.5, 1.0 - z, digamma(0.5 * a), 0.0, digamma(0.5 + 0.5 * a), + 0.0); + }; + expect_ad_grad_reg_inc_beta(f); +} - std::vector y1{b.val_}; - std::vector grad1; - g_b.d_.grad(y1, grad1); - EXPECT_NEAR(0.079977766631361187517939795, grad1[0], 1e-4); +/** + * The overload writing both partials through out-parameters takes every + * argument as the same type. + */ +TEST(mathMixScalFun, grad_reg_inc_beta_out_params) { + auto f = [](const auto& a, const auto& b, const auto& z) { + using stan::math::digamma; + using T = stan::return_type_t; + const T a_t = a; + const T b_t = b; + const T z_t = z; + T g1 = 0; + T g2 = 0; + stan::math::grad_reg_inc_beta(g1, g2, a_t, b_t, z_t, digamma(a_t), + digamma(b_t), digamma(a_t + b_t), T(0)); + Eigen::Matrix ret(2); + ret << g1, g2; + return ret; + }; + expect_ad_grad_reg_inc_beta(f); } diff --git a/test/unit/math/mix/fun/inc_beta_test.cpp b/test/unit/math/mix/fun/inc_beta_test.cpp index 4f00deab729..8805333dddd 100644 --- a/test/unit/math/mix/fun/inc_beta_test.cpp +++ b/test/unit/math/mix/fun/inc_beta_test.cpp @@ -28,6 +28,29 @@ TEST(mathMixScalFun, inc_beta) { stan::test::expect_ad(f, nan, nan, nan); } +/** + * At z = 1 with b < 1 the density is infinite. The forward and reverse + * sweeps of fvar must not throw there. + */ +TEST(mathMixScalFun, inc_beta_upper_bound_infinite_density) { + using stan::math::fvar; + using stan::math::var; + // grad() runs every node on the stack; start from an empty one + stan::math::recover_memory(); + var a_v = 1.05; + var b_v = 0.35; + fvar a(a_v, 1.0); + fvar b(b_v, 0.0); + fvar z(1.0, 0.0); + fvar f; + EXPECT_NO_THROW(f = stan::math::inc_beta(a, b, z)); + EXPECT_EQ(1.0, f.val_.val()); + EXPECT_NO_THROW(f.val_.grad()); + stan::math::set_zero_all_adjoints(); + EXPECT_NO_THROW(f.d_.grad()); + stan::math::recover_memory(); +} + TEST(mathMixScalFun, inc_beta_vec) { auto f = [](const auto& x1, const auto& x2, const auto& x3) { return stan::math::inc_beta(x1, x2, x3); diff --git a/test/unit/math/rev/fun/inc_beta_test.cpp b/test/unit/math/rev/fun/inc_beta_test.cpp new file mode 100644 index 00000000000..a66e629482a --- /dev/null +++ b/test/unit/math/rev/fun/inc_beta_test.cpp @@ -0,0 +1,62 @@ +#include +#include +#include +#include +#include + +namespace { + +/** + * Shape gradients of inc_beta against a high-precision reference (mpmath, + * 80 digits), at points where a finite-difference check is not informative. + */ +void expect_inc_beta_grads(double a, double b, double z, double dda, double ddb, + double rtol) { + using stan::math::var; + var av = a; + var bv = b; + var zv = z; + var f = stan::math::inc_beta(av, bv, zv); + std::vector x{av, bv}; + std::vector g; + f.grad(x, g); + EXPECT_TRUE(std::isfinite(g[0])) << "a=" << a << " b=" << b << " z=" << z; + EXPECT_TRUE(std::isfinite(g[1])) << "a=" << a << " b=" << b << " z=" << z; + EXPECT_NEAR(dda, g[0], rtol * std::fabs(dda)) + << "a=" << a << " b=" << b << " z=" << z; + EXPECT_NEAR(ddb, g[1], rtol * std::fabs(ddb)) + << "a=" << a << " b=" << b << " z=" << z; +} + +} // namespace + +TEST_F(AgradRev, inc_beta_shape_gradients) { + expect_inc_beta_grads(2, 3, 0.25, -0.22805360232434637, 0.10692153005229139, + 1e-12); + // a + b above the beta(a, b) underflow limit + expect_inc_beta_grads(600, 600, 0.5, -0.01152047147304755, + 0.01152047147304755, 1e-12); + // z^a (1 - z)^b underflows while I is 1 + expect_inc_beta_grads(352, 590, 0.758, -3.3772717797123242e-138, + 4.5390666074606649e-138, 1e-12); + // integer b + expect_inc_beta_grads(10, 1, 0.5, -0.00067690154351557159, + 0.00226574342615196, 1e-12); +} + +/** + * At z = 1, I is 1 for all shapes, so the shape gradients are 0. With b < 1 + * the density, which is the z gradient, is infinite. + */ +TEST_F(AgradRev, inc_beta_upper_bound_infinite_density) { + using stan::math::var; + var a = 1.05; + var b = 0.35; + var z = 1.0; + var f = stan::math::inc_beta(a, b, z); + EXPECT_NO_THROW(f.grad()); + EXPECT_EQ(1.0, f.val()); + EXPECT_EQ(0.0, a.adj()); + EXPECT_EQ(0.0, b.adj()); + EXPECT_TRUE(std::isinf(z.adj()) && z.adj() > 0); +}