diff --git a/stan/math/prim/fun/inc_beta.hpp b/stan/math/prim/fun/inc_beta.hpp index f2c12f91328..de6f432a0a9 100644 --- a/stan/math/prim/fun/inc_beta.hpp +++ b/stan/math/prim/fun/inc_beta.hpp @@ -29,6 +29,69 @@ inline double inc_beta(double a, double b, double x) { return boost::math::ibeta(a, b, x, boost_policy_t<>()); } +namespace internal { + +/** + * Return the complement of the regularized incomplete beta function, + * 1 - I_x(a, b) = I_{1-x}(b, a), from both x and 1 - x. + * + * Boost forms 1 - x from its argument, so an argument close to 1 loses + * relative precision in 1 - x. This overload calls Boost with whichever + * of x and 1 - x is at most 1/2: ibetac(a, b, x) or ibeta(b, a, 1 - x). + * Both keep their relative precision when the complement is below eps. + * + * Below the mean a / (a + b), while I_x(a, b) is at most 1/2, the + * complement is formed as 1 - I_x(a, b), which then has full precision. + * This also avoids Boost's arcsine case a = b = 1/2, where ibetac + * evaluates asin(sqrt(1 - x)) and loses digits for small x. + * + * @param a first shape, a > 0 + * @param b second shape, b > 0 + * @param x argument, 0 <= x <= 1 + * @param one_m_x 1 - x, computed by the caller without cancellation + * @return 1 - I_x(a, b) + */ +inline double inc_beta_complement(double a, double b, double x, + double one_m_x) { + check_not_nan("inc_beta_complement", "a", a); + check_not_nan("inc_beta_complement", "b", b); + check_not_nan("inc_beta_complement", "x", x); + if (x > 0.5) { + return boost::math::ibeta(b, a, one_m_x, boost_policy_t<>()); + } + if (x < a / (a + b)) { + const double inc = boost::math::ibeta(a, b, x, boost_policy_t<>()); + if (inc <= 0.5) { + return 1.0 - inc; + } + } + return boost::math::ibetac(a, b, x, boost_policy_t<>()); +} + +/** + * Return the complement of the regularized incomplete beta function for + * autodiff arguments, by the symmetry relation through inc_beta. + * + * @tparam T1 type of a + * @tparam T2 type of b + * @tparam T3 type of x + * @tparam T4 type of one_m_x + * @param a first shape, a > 0 + * @param b second shape, b > 0 + * @param x argument, 0 <= x <= 1; not used + * @param one_m_x 1 - x + * @return 1 - I_x(a, b) + */ +template * = nullptr> +inline return_type_t inc_beta_complement(const T1& a, const T2& b, + const T3& /* x */, + const T4& one_m_x) { + return inc_beta(b, a, one_m_x); +} + +} // namespace internal + /** * Enables the vectorized application of the inc_beta function, when * any arguments are containers. diff --git a/stan/math/prim/prob/beta_lccdf.hpp b/stan/math/prim/prob/beta_lccdf.hpp index ed4234a6f5a..172e08dff1c 100644 --- a/stan/math/prim/prob/beta_lccdf.hpp +++ b/stan/math/prim/prob/beta_lccdf.hpp @@ -7,6 +7,7 @@ #include #include #include +#include #include #include #include @@ -73,38 +74,71 @@ inline return_type_t beta_lccdf( constexpr bool calc_beta = is_autodiff_v; for (size_t n = 0; n < N; n++) { const auto y_dbl = y_vec.val(n); + + // Explicit results for extreme values + // The gradients are technically ill-defined, but treated as zero + if (y_dbl >= 1.0) { + return ops_partials.build(negative_infinity()); + } + const auto alpha_dbl = alpha_vec.val(n); const auto beta_dbl = beta_vec.val(n); - const auto betafunc_dbl = beta(alpha_dbl, beta_dbl); - const auto Pn = 1.0 - inc_beta(alpha_dbl, beta_dbl, y_dbl); + // The complement 1 - I_y(a, b), evaluated from the smaller of y and + // 1 - y. Forming it as a difference loses every digit once I_y(a, b) + // rounds to 1, which happens for every complement below eps. + const auto Pn = internal::inc_beta_complement(alpha_dbl, beta_dbl, y_dbl, + 1.0 - y_dbl); const auto inv_Pn = is_constant_all::value ? 0 : inv(Pn); ccdf_log += log(Pn); if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[n] -= pow(1 - y_dbl, beta_dbl - 1) - * pow(y_dbl, alpha_dbl - 1) * inv_Pn - / betafunc_dbl; + // inc_beta_ddz is the same density, evaluated without the product of + // two powers over a beta function, each of which can underflow on + // its own. + partials<0>(ops_partials)[n] + -= inc_beta_ddz(alpha_dbl, beta_dbl, y_dbl) * inv_Pn; } + // The shape derivatives from the same argument as the value. For + // y <= 1/2 they are those of I_y(a, b), which are those of -Pn; on the + // reflected arguments they are those of Pn itself, with beta the first + // shape and alpha the second. + [[maybe_unused]] const bool reflect = y_dbl > 0.5; if constexpr (calc_alpha && calc_beta) { - const auto [g1, g2] = grad_reg_inc_beta( - alpha_dbl, beta_dbl, y_dbl, digamma(alpha_dbl), digamma(beta_dbl), - digamma(alpha_dbl + beta_dbl), 0.0); - partials<1>(ops_partials)[n] -= g1 * inv_Pn; - partials<2>(ops_partials)[n] -= g2 * inv_Pn; + if (reflect) { + const auto [g1, g2] = grad_reg_inc_beta( + beta_dbl, alpha_dbl, 1.0 - y_dbl, digamma(beta_dbl), + digamma(alpha_dbl), digamma(alpha_dbl + beta_dbl), 0.0); + partials<1>(ops_partials)[n] += g2 * inv_Pn; + partials<2>(ops_partials)[n] += g1 * inv_Pn; + } else { + const auto [g1, g2] = grad_reg_inc_beta( + alpha_dbl, beta_dbl, y_dbl, digamma(alpha_dbl), digamma(beta_dbl), + digamma(alpha_dbl + beta_dbl), 0.0); + partials<1>(ops_partials)[n] -= g1 * inv_Pn; + partials<2>(ops_partials)[n] -= g2 * inv_Pn; + } } else if constexpr (calc_alpha) { + const auto digamma_sum = digamma(alpha_dbl + beta_dbl); partials<1>(ops_partials)[n] - -= grad_reg_inc_beta(alpha_dbl, beta_dbl, y_dbl, - digamma(alpha_dbl), 0.0, - digamma(alpha_dbl + beta_dbl), 0.0) + += (reflect ? grad_reg_inc_beta( + beta_dbl, alpha_dbl, 1.0 - y_dbl, 0.0, + digamma(alpha_dbl), digamma_sum, 0.0) + : -grad_reg_inc_beta( + alpha_dbl, beta_dbl, y_dbl, digamma(alpha_dbl), + 0.0, digamma_sum, 0.0)) * inv_Pn; } else if constexpr (calc_beta) { + const auto digamma_sum = digamma(alpha_dbl + beta_dbl); partials<2>(ops_partials)[n] - -= grad_reg_inc_beta(alpha_dbl, beta_dbl, y_dbl, 0.0, - digamma(beta_dbl), - digamma(alpha_dbl + beta_dbl), 0.0) + += (reflect ? grad_reg_inc_beta( + beta_dbl, alpha_dbl, 1.0 - y_dbl, + digamma(beta_dbl), 0.0, digamma_sum, 0.0) + : -grad_reg_inc_beta( + alpha_dbl, beta_dbl, y_dbl, 0.0, digamma(beta_dbl), + digamma_sum, 0.0)) * inv_Pn; } } diff --git a/stan/math/prim/prob/beta_proportion_lccdf.hpp b/stan/math/prim/prob/beta_proportion_lccdf.hpp index 04e606fd3c8..42773d93274 100644 --- a/stan/math/prim/prob/beta_proportion_lccdf.hpp +++ b/stan/math/prim/prob/beta_proportion_lccdf.hpp @@ -7,6 +7,7 @@ #include #include #include +#include #include #include #include @@ -75,13 +76,22 @@ inline return_type_t beta_proportion_lccdf( for (size_t n = 0; n < N; n++) { const T_partials_return y_dbl = y_vec.val(n); + + // Explicit results for extreme values + // The gradients are technically ill-defined, but treated as zero + if (y_dbl >= 1.0) { + return ops_partials.build(negative_infinity()); + } + const T_partials_return mu_dbl = mu_vec.val(n); const T_partials_return kappa_dbl = kappa_vec.val(n); const T_partials_return mukappa_dbl = mu_dbl * kappa_dbl; const T_partials_return kappa_mukappa_dbl = kappa_dbl - mukappa_dbl; - const T_partials_return betafunc_dbl = beta(mukappa_dbl, kappa_mukappa_dbl); - const T_partials_return Pn - = 1 - inc_beta(mukappa_dbl, kappa_mukappa_dbl, y_dbl); + // The complement 1 - I_y(a, b), evaluated from the smaller of y and + // 1 - y. Forming it as a difference loses every digit once I_y(a, b) + // rounds to 1, which happens for every complement below eps. + const T_partials_return Pn = internal::inc_beta_complement( + mukappa_dbl, kappa_mukappa_dbl, y_dbl, 1.0 - y_dbl); ccdf_log += log(Pn); @@ -89,22 +99,41 @@ inline return_type_t beta_proportion_lccdf( = is_constant_all::value ? 0 : inv(Pn); if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[n] -= pow(1 - y_dbl, kappa_mukappa_dbl - 1) - * pow(y_dbl, mukappa_dbl - 1) * inv_Pn - / betafunc_dbl; + // inc_beta_ddz is the same density, evaluated without the product of + // two powers over a beta function, each of which can underflow on + // its own. + partials<0>(ops_partials)[n] + -= inc_beta_ddz(mukappa_dbl, kappa_mukappa_dbl, y_dbl) * inv_Pn; } // the mu and kappa partials each need both shape partials if constexpr (is_any_autodiff_v) { - const auto [g1, g2] = grad_reg_inc_beta( - mukappa_dbl, kappa_mukappa_dbl, y_dbl, digamma(mukappa_dbl), - digamma(kappa_mukappa_dbl), digamma(kappa_dbl), 0.0); + // The derivatives with respect to a = mu kappa and b = kappa - mu kappa, + // from the same argument as the value + T_partials_return dPn_da = 0; + T_partials_return dPn_db = 0; + if (y_dbl <= 0.5) { + // the derivatives of I_y(a, b), which are those of -Pn + const auto [g1, g2] = grad_reg_inc_beta( + mukappa_dbl, kappa_mukappa_dbl, y_dbl, digamma(mukappa_dbl), + digamma(kappa_mukappa_dbl), digamma(kappa_dbl), 0.0); + dPn_da = -g1; + dPn_db = -g2; + } else { + // on the reflected arguments, the derivatives of Pn itself + const auto [g1, g2] = grad_reg_inc_beta( + kappa_mukappa_dbl, mukappa_dbl, 1 - y_dbl, + digamma(kappa_mukappa_dbl), digamma(mukappa_dbl), + digamma(kappa_dbl), 0.0); + dPn_da = g2; + dPn_db = g1; + } if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] -= kappa_dbl * (g1 - g2) * inv_Pn; + partials<1>(ops_partials)[n] += kappa_dbl * (dPn_da - dPn_db) * inv_Pn; } if constexpr (is_autodiff_v) { partials<2>(ops_partials)[n] - -= (g1 * mu_dbl + g2 * (1 - mu_dbl)) * inv_Pn; + += (dPn_da * mu_dbl + dPn_db * (1 - mu_dbl)) * inv_Pn; } } } diff --git a/stan/math/prim/prob/binomial_lccdf.hpp b/stan/math/prim/prob/binomial_lccdf.hpp index 5fe97006126..cd77ba7a3f7 100644 --- a/stan/math/prim/prob/binomial_lccdf.hpp +++ b/stan/math/prim/prob/binomial_lccdf.hpp @@ -6,6 +6,7 @@ #include #include #include +#include #include #include #include @@ -89,16 +90,20 @@ inline return_type_t binomial_lccdf(const T_n& n, const T_N& N, } const T_partials_return theta_dbl = theta_vec.val(i); - const T_partials_return Pi - = 1.0 - inc_beta(N_dbl - n_dbl, n_dbl + 1, 1 - theta_dbl); + // The complement by the symmetry relation I_z(a, b) = I_{1-z}(b, a), + // which also removes the 1 - theta argument. Forming it as + // 1 - I_{1-theta}(N - n, n + 1) loses every digit once that rounds to + // 1, which happens for every complement below eps. + const T_partials_return Pi = inc_beta(n_dbl + 1, N_dbl - n_dbl, theta_dbl); P += log(Pi); if constexpr (is_autodiff_v) { - const T_partials_return denom = beta(N_dbl - n_dbl, n_dbl + 1) * Pi; - partials<0>(ops_partials)[i] += pow(theta_dbl, n_dbl) - * pow(1 - theta_dbl, N_dbl - n_dbl - 1) - / denom; + // inc_beta_ddz is the same density, evaluated without the product of + // two powers over a beta function, each of which can underflow on + // its own. + partials<0>(ops_partials)[i] + += inc_beta_ddz(n_dbl + 1, N_dbl - n_dbl, theta_dbl) / Pi; } } diff --git a/stan/math/prim/prob/neg_binomial_lccdf.hpp b/stan/math/prim/prob/neg_binomial_lccdf.hpp index 3de84a34589..efa46930fe9 100644 --- a/stan/math/prim/prob/neg_binomial_lccdf.hpp +++ b/stan/math/prim/prob/neg_binomial_lccdf.hpp @@ -8,6 +8,7 @@ #include #include #include +#include #include #include #include @@ -76,23 +77,47 @@ inline return_type_t neg_binomial_lccdf( const auto inv_beta_p1 = inv(beta_dbl + 1); const auto p_dbl = beta_dbl * inv_beta_p1; const auto d_dbl = square(inv_beta_p1); - const auto Pi = 1.0 - inc_beta(alpha_dbl, n_dbl + 1.0, p_dbl); - const auto beta_func = beta(n_dbl + 1, alpha_dbl); + // The complement 1 - I_p(alpha, n + 1), evaluated from the smaller of + // p and 1 - p = 1 / (beta + 1), both formed without cancellation. + // Forming it as a difference loses every digit once I_p(alpha, n + 1) + // rounds to 1, which happens for every complement below eps. + const auto Pi = internal::inc_beta_complement(alpha_dbl, n_dbl + 1.0, p_dbl, + inv_beta_p1); P += log(Pi); - // alpha is the first shape of inc_beta, so only its partial is needed + // Only the alpha partial is needed, from the same argument as the value if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[i] - -= grad_reg_inc_beta( - alpha_dbl, n_dbl + 1, p_dbl, digamma(alpha_dbl), 0.0, - digamma(n_dbl + alpha_dbl + 1), 0.0) - / Pi; + const auto digamma_alpha = digamma(alpha_dbl); + const auto digamma_sum = digamma(n_dbl + alpha_dbl + 1); + if (p_dbl <= 0.5) { + // alpha is the first shape of I_p(alpha, n + 1), whose derivative + // is that of -Pi + partials<0>(ops_partials)[i] + -= grad_reg_inc_beta(alpha_dbl, n_dbl + 1, p_dbl, + digamma_alpha, 0.0, digamma_sum, + 0.0) + / Pi; + } else { + // on the reflected arguments alpha is the second shape, and the + // derivative is that of Pi itself + partials<0>(ops_partials)[i] + += grad_reg_inc_beta(n_dbl + 1, alpha_dbl, + inv_beta_p1, 0.0, digamma_alpha, + digamma_sum, 0.0) + / Pi; + } } if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[i] -= d_dbl * pow(1 - p_dbl, n_dbl) - * pow(p_dbl, alpha_dbl - 1) - / (beta_func * Pi); + // inc_beta_ddz is the same density, evaluated without the product of + // two powers over a beta function, each of which can underflow on + // its own. The density at p equals the reflected density at 1 - p; + // it is evaluated at whichever of the two is at most 1/2, because + // inc_beta_ddz forms the complement of its argument internally. + const T_partials_return density + = p_dbl <= 0.5 ? inc_beta_ddz(alpha_dbl, n_dbl + 1.0, p_dbl) + : inc_beta_ddz(n_dbl + 1.0, alpha_dbl, inv_beta_p1); + partials<1>(ops_partials)[i] -= d_dbl * density / Pi; } } diff --git a/stan/math/prim/prob/student_t_cdf.hpp b/stan/math/prim/prob/student_t_cdf.hpp index 78cfa8f5ee3..cde85c298cc 100644 --- a/stan/math/prim/prob/student_t_cdf.hpp +++ b/stan/math/prim/prob/student_t_cdf.hpp @@ -72,15 +72,18 @@ inline return_type_t student_t_cdf( const auto nu_dbl = nu_vec.val(n); const auto q = nu_dbl / (t * t); const auto r = 1.0 / (1.0 + q); + // 1 - r = 1 / (1 + 1 / q) without the subtraction, which rounds to 0 + // once q is below eps; equal to 1 at t = 0 + const auto one_m_r = 1.0 / (1.0 + t * t / nu_dbl); const auto J = 2 * r * r * q / t; const auto betaNuHalf = beta(0.5, 0.5 * nu_dbl); double zJacobian = t > 0 ? -0.5 : 0.5; if (q < 2) { - auto z = inc_beta(0.5 * nu_dbl, 0.5, 1.0 - r); + auto z = inc_beta(0.5 * nu_dbl, 0.5, one_m_r); const auto Pn = t > 0 ? 1.0 - 0.5 * z : 0.5 * z; const auto d_ibeta - = pow(r, -0.5) * pow(1.0 - r, 0.5 * nu_dbl - 1) / betaNuHalf; + = pow(r, -0.5) * pow(one_m_r, 0.5 * nu_dbl - 1) / betaNuHalf; P *= Pn; @@ -91,7 +94,7 @@ inline return_type_t student_t_cdf( if constexpr (is_autodiff_v) { // nu / 2 is the first shape of inc_beta here and 1 / 2 the second const auto g1 = grad_reg_inc_beta( - 0.5 * nu_dbl, 0.5, 1.0 - r, digamma(0.5 * nu_dbl), 0.0, + 0.5 * nu_dbl, 0.5, one_m_r, digamma(0.5 * nu_dbl), 0.0, digamma(0.5 + 0.5 * nu_dbl), 0.0); partials<1>(ops_partials)[n] @@ -108,13 +111,16 @@ inline return_type_t student_t_cdf( } } else { - auto z = 1.0 - inc_beta(0.5, 0.5 * nu_dbl, r); - + // Same quantity as in the other branch, 1 - I_r(1/2, nu/2), with the + // argument r <= 1/3. Forming it as a difference loses every digit once + // I_r(1/2, nu/2) rounds to 1, which happens for every z below eps. + T_partials_return z = internal::inc_beta_complement( + (T_partials_return)0.5, 0.5 * nu_dbl, r, one_m_r); zJacobian *= -1; const auto Pn = t > 0 ? 1.0 - 0.5 * z : 0.5 * z; - auto d_ibeta = pow(1.0 - r, 0.5 * nu_dbl - 1) * pow(r, -0.5) / betaNuHalf; + auto d_ibeta = pow(one_m_r, 0.5 * nu_dbl - 1) * pow(r, -0.5) / betaNuHalf; P *= Pn; diff --git a/stan/math/prim/prob/student_t_lccdf.hpp b/stan/math/prim/prob/student_t_lccdf.hpp index 224a3165f50..f1f8ae552b4 100644 --- a/stan/math/prim/prob/student_t_lccdf.hpp +++ b/stan/math/prim/prob/student_t_lccdf.hpp @@ -74,15 +74,18 @@ inline return_type_t student_t_lccdf( const auto nu_dbl = nu_vec.val(n); const auto q = nu_dbl / (t * t); const auto r = 1.0 / (1.0 + q); + // 1 - r = 1 / (1 + 1 / q) without the subtraction, which rounds to 0 + // once q is below eps; equal to 1 at t = 0 + const auto one_m_r = 1.0 / (1.0 + t * t / nu_dbl); const auto J = 2 * r * r * q / t; const auto betaNuHalf = beta(0.5, 0.5 * nu_dbl); auto zJacobian = t > 0 ? -0.5 : 0.5; if (q < 2) { - auto z = inc_beta(0.5 * nu_dbl, 0.5, 1.0 - r); + auto z = inc_beta(0.5 * nu_dbl, 0.5, one_m_r); const auto Pn = t > 0 ? 0.5 * z : 1.0 - 0.5 * z; const auto d_ibeta - = pow(r, -0.5) * pow(1.0 - r, 0.5 * nu_dbl - 1) / betaNuHalf; + = pow(r, -0.5) * pow(one_m_r, 0.5 * nu_dbl - 1) / betaNuHalf; P += log(Pn); @@ -94,7 +97,7 @@ inline return_type_t student_t_lccdf( if constexpr (is_autodiff_v) { // nu / 2 is the first shape of inc_beta here and 1 / 2 the second const auto g1 = grad_reg_inc_beta( - 0.5 * nu_dbl, 0.5, 1.0 - r, digamma(0.5 * nu_dbl), 0.0, + 0.5 * nu_dbl, 0.5, one_m_r, digamma(0.5 * nu_dbl), 0.0, digamma(0.5 + 0.5 * nu_dbl), 0.0); partials<1>(ops_partials)[n] @@ -111,12 +114,17 @@ inline return_type_t student_t_lccdf( } } else { - auto z = 1.0 - inc_beta(0.5, 0.5 * nu_dbl, r); + // Same quantity as in the other branch, 1 - I_r(1/2, nu/2), with the + // argument r <= 1/3. Forming it as a difference loses every digit once + // I_r(1/2, nu/2) rounds to 1, which happens for every z below eps, + // that is for every tail probability below about e^-37. + T_partials_return z = internal::inc_beta_complement( + (T_partials_return)0.5, 0.5 * nu_dbl, r, one_m_r); zJacobian *= -1; const auto Pn = t > 0 ? 0.5 * z : 1.0 - 0.5 * z; - auto d_ibeta = pow(1.0 - r, 0.5 * nu_dbl - 1) * pow(r, -0.5) / betaNuHalf; + auto d_ibeta = pow(one_m_r, 0.5 * nu_dbl - 1) * pow(r, -0.5) / betaNuHalf; P += log(Pn); diff --git a/stan/math/prim/prob/student_t_lcdf.hpp b/stan/math/prim/prob/student_t_lcdf.hpp index be791924977..f1d57d340c3 100644 --- a/stan/math/prim/prob/student_t_lcdf.hpp +++ b/stan/math/prim/prob/student_t_lcdf.hpp @@ -74,15 +74,18 @@ inline return_type_t student_t_lcdf( const auto nu_dbl = nu_vec.val(n); const auto q = nu_dbl / (t * t); const auto r = 1.0 / (1.0 + q); + // 1 - r = 1 / (1 + 1 / q) without the subtraction, which rounds to 0 + // once q is below eps; equal to 1 at t = 0 + const auto one_m_r = 1.0 / (1.0 + t * t / nu_dbl); const auto J = 2 * r * r * q / t; const auto betaNuHalf = beta(0.5, 0.5 * nu_dbl); auto zJacobian = t > 0 ? -0.5 : 0.5; if (q < 2) { - auto z = inc_beta(0.5 * nu_dbl, 0.5, 1.0 - r); + auto z = inc_beta(0.5 * nu_dbl, 0.5, one_m_r); const auto Pn = t > 0 ? 1.0 - 0.5 * z : 0.5 * z; const auto d_ibeta - = pow(r, -0.5) * pow(1.0 - r, 0.5 * nu_dbl - 1) / betaNuHalf; + = pow(r, -0.5) * pow(one_m_r, 0.5 * nu_dbl - 1) / betaNuHalf; P += log(Pn); @@ -94,7 +97,7 @@ inline return_type_t student_t_lcdf( if constexpr (is_autodiff_v) { // nu / 2 is the first shape of inc_beta here and 1 / 2 the second const auto g1 = grad_reg_inc_beta( - 0.5 * nu_dbl, 0.5, 1.0 - r, digamma(0.5 * nu_dbl), 0.0, + 0.5 * nu_dbl, 0.5, one_m_r, digamma(0.5 * nu_dbl), 0.0, digamma(0.5 + 0.5 * nu_dbl), 0.0); partials<1>(ops_partials)[n] @@ -111,12 +114,17 @@ inline return_type_t student_t_lcdf( } } else { - auto z = 1.0 - inc_beta(0.5, 0.5 * nu_dbl, r); + // Same quantity as in the other branch, 1 - I_r(1/2, nu/2), with the + // argument r <= 1/3. Forming it as a difference loses every digit once + // I_r(1/2, nu/2) rounds to 1, which happens for every z below eps, + // that is for every tail probability below about e^-37. + T_partials_return z = internal::inc_beta_complement( + (T_partials_return)0.5, 0.5 * nu_dbl, r, one_m_r); zJacobian *= -1; const auto Pn = t > 0 ? 1.0 - 0.5 * z : 0.5 * z; - auto d_ibeta = pow(1.0 - r, 0.5 * nu_dbl - 1) * pow(r, -0.5) / betaNuHalf; + auto d_ibeta = pow(one_m_r, 0.5 * nu_dbl - 1) * pow(r, -0.5) / betaNuHalf; P += log(Pn); diff --git a/test/unit/math/rev/prob/beta_ccdf_log_test.cpp b/test/unit/math/rev/prob/beta_ccdf_log_test.cpp new file mode 100644 index 00000000000..b10b0f4d04a --- /dev/null +++ b/test/unit/math/rev/prob/beta_ccdf_log_test.cpp @@ -0,0 +1,196 @@ +#include +#include +#include +#include +#include +#include + +namespace { + +// beta_lccdf and its shape gradients against a high-precision reference +// (mpmath, 50 digits). The complement 1 - I_y(a, b) is below eps at these +// points, so the value and both gradients are the quantities that a +// linear-space complement cannot represent. +void expect_beta_lccdf(double y, double alpha, double beta, double expected, + double dalpha, double dbeta, double rtol) { + using stan::math::var; + var a = alpha; + var b = beta; + var lp = stan::math::beta_lccdf(y, a, b); + std::vector x{a, b}; + std::vector g; + lp.grad(x, g); + EXPECT_TRUE(std::isfinite(lp.val())) + << "y=" << y << " a=" << alpha << " b=" << beta; + EXPECT_TRUE(std::isfinite(g[0]) && std::isfinite(g[1])) + << "y=" << y << " a=" << alpha << " b=" << beta; + EXPECT_NEAR(expected, lp.val(), rtol * std::fabs(expected)); + EXPECT_NEAR(dalpha, g[0], rtol * std::fabs(dalpha)); + EXPECT_NEAR(dbeta, g[1], rtol * std::fabs(dbeta)); +} + +// neg_binomial_lccdf and its alpha and beta gradients against references +// (mpmath, 60 digits, two routes). A value close to 0 is the log of a +// probability that rounds to 1, so the value also has an absolute floor. +void expect_neg_binomial_lccdf(int n, double alpha_in, double beta_in, + double expected, double dalpha, double dbeta, + double rtol) { + using stan::math::var; + var alpha = alpha_in; + var beta = beta_in; + var lp = stan::math::neg_binomial_lccdf(n, alpha, beta); + std::vector x{alpha, beta}; + std::vector g; + lp.grad(x, g); + EXPECT_NEAR(expected, lp.val(), std::max(rtol * std::fabs(expected), 1e-15)) + << "n=" << n << " beta=" << beta_in; + EXPECT_NEAR(dalpha, g[0], rtol * std::fabs(dalpha)) + << "n=" << n << " beta=" << beta_in; + EXPECT_NEAR(dbeta, g[1], rtol * std::fabs(dbeta)) + << "n=" << n << " beta=" << beta_in; +} + +// beta_proportion_lccdf and its mu and kappa gradients, as above +void expect_beta_proportion_lccdf(double y, double mu_in, double kappa_in, + double expected, double dmu, double dkappa, + double rtol) { + using stan::math::var; + var mu = mu_in; + var kappa = kappa_in; + var lp = stan::math::beta_proportion_lccdf(y, mu, kappa); + std::vector x{mu, kappa}; + std::vector g; + lp.grad(x, g); + EXPECT_NEAR(expected, lp.val(), std::max(rtol * std::fabs(expected), 1e-15)) + << "y=" << y; + EXPECT_NEAR(dmu, g[0], rtol * std::fabs(dmu)) << "y=" << y; + EXPECT_NEAR(dkappa, g[1], rtol * std::fabs(dkappa)) << "y=" << y; +} + +} // namespace + +TEST_F(AgradRev, ProbDistributionsBeta_lccdf_deep_tail) { + // moderate shapes, y in the upper tail: the true value is -54, which the + // complement cannot reach in linear space because eps is e^-36.7 + expect_beta_lccdf(0.9, 50, 50, -54.088884708, 0.595231340835, -1.62672723211, + 1e-10); + expect_beta_lccdf(0.8, 50, 50, -25.04429524, 0.481321, -0.937243, 1e-6); + expect_beta_lccdf(0.758, 352, 590, -316.195574221, 0.708857759363, + -0.952707627587, 1e-10); + // still inside the representable range of the old form, unchanged + expect_beta_lccdf(0.99, 5, 5, -18.2230299025, 0.737273747257, -4.06065905449, + 1e-10); +} + +TEST_F(AgradRev, ProbDistributionsBeta_proportion_lccdf_deep_tail) { + using stan::math::var; + var mu = 0.4; + var kappa = 900; + var lp = stan::math::beta_proportion_lccdf(0.758, mu, kappa); + lp.grad(); + EXPECT_TRUE(std::isfinite(lp.val())); + EXPECT_TRUE(std::isfinite(mu.adj()) && std::isfinite(kappa.adj())); + EXPECT_NEAR(-264.2052063, lp.val(), 1e-7 * 264.0); +} + +TEST_F(AgradRev, ProbDistributionsNegBinomial_lccdf_deep_tail) { + using stan::math::var; + { + std::vector n{50}; + var alpha = 2; + var beta = 99; + var lp = stan::math::neg_binomial_lccdf(n, alpha, beta); + std::vector x{alpha, beta}; + std::vector g; + lp.grad(x, g); + EXPECT_NEAR(-230.9222919, lp.val(), 1e-9 * 231.0); + EXPECT_NEAR(3.528187864, g[0], 1e-9 * 3.53); + EXPECT_NEAR(-0.5099009516, g[1], 1e-9 * 0.51); + } + { + std::vector n{50}; + var alpha = 2; + var beta = 999; + var lp = stan::math::neg_binomial_lccdf(n, alpha, beta); + std::vector x{alpha, beta}; + std::vector g; + lp.grad(x, g); + EXPECT_NEAR(-348.3452568, lp.val(), 1e-9 * 348.0); + EXPECT_NEAR(3.5370627, g[0], 1e-9 * 3.54); + EXPECT_NEAR(-0.05099901827, g[1], 1e-9 * 0.051); + } +} + +TEST_F(AgradRev, ProbDistributionsNegBinomial_lccdf_large_beta) { + // stan-dev/math#2031: at large beta, p = beta / (beta + 1) rounds towards + // 1. The complement is formed from 1 / (beta + 1); the density factor of + // the beta partial must be too, or it is 0 for n >= 1 once p rounds to 1 + // (beta = 1e16, 1e18) and loses digits before that (beta = 1e8). The last + // two rows are controls on both sides of beta = 1. + const double rtol = 1e-12; + // the reported case, log1m_exp(-alpha * log1p(1 / beta)) + expect_neg_binomial_lccdf(0, 1.0, 1e18, -41.446531673892821, 1.0, + -1.0000000000000001e-18, rtol); + expect_neg_binomial_lccdf(1, 2.0, 1e18, -81.794451059117534, + 0.83333333333333337, -2.0000000000000001e-18, rtol); + expect_neg_binomial_lccdf(1, 2.0, 1e16, -72.584110687141347, + 0.83333333333333326, -1.9999999999999997e-16, rtol); + expect_neg_binomial_lccdf(3, 2.0, 1e8, -72.073285111375355, + 1.2833333253333334, -3.9999999520000002e-08, rtol); + expect_neg_binomial_lccdf(5, 0.5, 1e-3, -0.08928900804336741, + 0.48141786839236278, -46.496417635081194, rtol); + expect_neg_binomial_lccdf(5, 3.0, 2.0, -3.9290859049832054, + 0.8864463205559171, -1.7364341085271318, rtol); +} + +TEST_F(AgradRev, ProbDistributionsBeta_lccdf_small_y) { + // The complement 1 - I_y(a, b) and its shape gradients must be evaluated + // from y, not from 1 - y, which is close to 1 for small y: from 1 - y the + // alpha gradient at y = 1e-10 had a relative error of 4e-8. + expect_beta_lccdf(1e-10, 0.5, 3, -1.8750175782197275e-05, + 0.00041174242508107103, -3.3820441414792342e-06, 1e-11); + expect_beta_lccdf(1e-6, 0.5, 0.5, -0.00063682260695097634, + 0.0091917773883985632, -0.00088310453739762977, 1e-11); + // y <= 1/2 with the complement far below eps + expect_beta_lccdf(0.3, 1, 200, -71.334988787746468, 4.6855356557158476, + -0.35667494393873234, 1e-12); + expect_beta_lccdf(0.45, 50, 450, -150.59662983495602, 1.5162204045412182, + -0.49492426185967869, 1e-12); +} + +TEST_F(AgradRev, ProbDistributionsBeta_proportion_lccdf_small_y) { + // as for beta_lccdf; the value -8.4e-23 is the log of a probability that + // rounds to 1 + expect_beta_proportion_lccdf(1e-8, 0.3, 10, -8.3999996220000247e-23, + 1.4955371148778298e-20, 4.1682780319251551e-22, + 1e-11); + expect_beta_proportion_lccdf(0.4, 0.02, 500, -215.48431197855072, + 1770.6978109683873, -0.42187260226452955, 1e-12); +} + +TEST_F(AgradRev, ProbDistributionsNegBinomial_lccdf_small_beta) { + // beta < 1, so p = beta / (beta + 1) < 1/2: the complement and the alpha + // gradient must be evaluated from p, not from 1 / (beta + 1), which is + // close to 1: from 1 / (beta + 1) the alpha gradient at beta = 1e-9 had a + // relative error of 4e-8. The last two rows have the complement far below + // eps. + expect_neg_binomial_lccdf(5, 2.0, 1e-6, -2.0999888000598496e-11, + 2.6717432945696996e-10, -4.1999664002393978e-05, + 1e-11); + expect_neg_binomial_lccdf(5, 0.5, 1e-9, -8.56075085052458e-05, + 0.0016237738030832265, -42805.586280793519, 1e-11); + expect_neg_binomial_lccdf(100, 1.0, 0.5, -40.951975918924603, + 4.1179071920897643, -67.333333333333329, 1e-12); + expect_neg_binomial_lccdf(300, 2.0, 0.2, -50.943700316621204, + 3.5140890478701055, -246.74809989142236, 1e-12); +} + +TEST_F(AgradRev, ProbDistributionsBinomial_lccdf_deep_tail) { + using stan::math::var; + std::vector n{500}; + var theta = 0.2; + var lp = stan::math::binomial_lccdf(n, 1000, theta); + lp.grad(); + EXPECT_NEAR(-227.9265046, lp.val(), 1e-9 * 228.0); + EXPECT_NEAR(1883.30957, theta.adj(), 1e-8 * 1884.0); +} diff --git a/test/unit/math/rev/prob/student_t_ccdf_log_test.cpp b/test/unit/math/rev/prob/student_t_ccdf_log_test.cpp new file mode 100644 index 00000000000..9b0154e94f7 --- /dev/null +++ b/test/unit/math/rev/prob/student_t_ccdf_log_test.cpp @@ -0,0 +1,169 @@ +#include +#include +#include +#include +#include + +namespace { + +// student_t_lccdf and its nu and sigma gradients against a high-precision +// reference (mpmath, 50 digits, the tail as 0.5 I_{nu/(nu+t^2)}(nu/2, 1/2) +// differentiated numerically). +void expect_student_t_lccdf(double y, double nu_in, double sigma_in, + double expected, double dnu, double dsigma, + double rtol) { + using stan::math::var; + var nu = nu_in; + var sigma = sigma_in; + var lp = stan::math::student_t_lccdf(y, nu, 0.0, sigma); + std::vector x{nu, sigma}; + std::vector g; + lp.grad(x, g); + EXPECT_TRUE(std::isfinite(lp.val())) << "y=" << y << " nu=" << nu_in; + EXPECT_TRUE(std::isfinite(g[0]) && std::isfinite(g[1])) + << "y=" << y << " nu=" << nu_in; + EXPECT_NEAR(expected, lp.val(), rtol * std::fabs(expected)); + EXPECT_NEAR(dnu, g[0], rtol * std::fabs(dnu)); + EXPECT_NEAR(dsigma, g[1], rtol * std::fabs(dsigma)); +} + +// f(y, nu, sigma) with location 0: the value, then the y, nu and sigma +// gradients, against references (mpmath, 60 digits, two routes) +template +void expect_student_t_fn(F f, double y_in, double nu_in, double sigma_in, + const std::vector& expected, double rtol) { + using stan::math::var; + var y = y_in; + var nu = nu_in; + var sigma = sigma_in; + var p = f(y, nu, sigma); + std::vector x{y, nu, sigma}; + std::vector g; + p.grad(x, g); + EXPECT_NEAR(expected[0], p.val(), rtol * std::fabs(expected[0])) + << "y=" << y_in << " nu=" << nu_in; + for (int i = 0; i < 3; ++i) { + EXPECT_NEAR(expected[i + 1], g[i], rtol * std::fabs(expected[i + 1])) + << "gradient " << i << " y=" << y_in << " nu=" << nu_in; + } +} + +const auto student_t_cdf_fn = [](auto& y, auto& nu, auto& sigma) { + return stan::math::student_t_cdf(y, nu, 0.0, sigma); +}; +const auto student_t_lcdf_fn = [](auto& y, auto& nu, auto& sigma) { + return stan::math::student_t_lcdf(y, nu, 0.0, sigma); +}; +const auto student_t_lccdf_fn = [](auto& y, auto& nu, auto& sigma) { + return stan::math::student_t_lccdf(y, nu, 0.0, sigma); +}; + +} // namespace + +TEST_F(AgradRev, ProbDistributionsStudentT_lccdf_deep_tail) { + // q = nu / t^2 >= 2, the branch that formed the complement in linear + // space; the true value is -50.8, below eps = e^-36.7 + expect_student_t_lccdf(10, 1000, 1, -50.838951454, -0.0022457674682, + 91.799148523, 1e-9); + expect_student_t_lccdf(8, 200, 1, -30.638906854, -0.018204242375, + 49.213616295, 1e-9); + // q < 2, the other branch, unchanged + expect_student_t_lccdf(20, 50, 1, -57.754062809, -0.66295856977, 44.550802193, + 1e-9); + expect_student_t_lccdf(50, 1000, 1, -630.58671311, -0.26959552076, + 714.57063175, 1e-9); + expect_student_t_lccdf(4, 4, 1, -4.8202159939, -0.49244669515, 3.3270509831, + 1e-9); +} + +TEST_F(AgradRev, ProbDistributionsStudentT_lccdf_large_t) { + // stan-dev/math#2923: q = nu / t^2 < 2, where 1 - r with r = 1 / (1 + q) + // was formed by subtraction. It loses digits from |t| = 10 on and is 0 + // once q is below eps, which made the tail probability 0. + expect_student_t_lccdf(1e8, 0.1, 1, -2.7158255239161981, -19.691030420120541, + 0.10000000000000001, 1e-12); + expect_student_t_lccdf(1e29, 0.1, 1, -7.5512542192036944, -68.045317372995498, + 0.10000000000000001, 1e-12); + // Cauchy, t = 1e10 with sigma = 2.5 + expect_student_t_lccdf(2.5e10, 1, 2.5, -24.170580815789858, + -22.83270374938051, 0.40000000000000002, 1e-12); + // finite before the change, but with a relative error of 3e-5 + expect_student_t_lccdf(1e7, 4, 1, -63.373770315165238, -15.034762317625022, + 3.9999999999998668, 1e-12); +} + +TEST_F(AgradRev, ProbDistributionsStudentT_cdf_lcdf_large_t) { + // stan-dev/math#2923: the reported case nu = 0.1 at y = -1e8, where the + // cdf was exactly 0, its mirror at y = 1e8, where it was exactly 1, and + // the Cauchy lcdf at y = -1e10. Value, then the y, nu and sigma gradients. + expect_student_t_fn(student_t_cdf_fn, -1e8, 0.1, 1.0, + {0.066150321787786709, 6.6150321787786714e-11, + -1.3025679986240706, 0.0066150321787786714}, + 1e-12); + expect_student_t_fn(student_t_cdf_fn, 1e8, 0.1, 1.0, + {0.93384967821221332, 6.6150321787786714e-11, + 1.3025679986240706, -0.0066150321787786714}, + 1e-12); + expect_student_t_fn(student_t_lcdf_fn, -1e10, 1.0, 1.0, + {-24.170580815789858, 1e-10, -22.83270374938051, 1.0}, + 1e-12); +} + +TEST_F(AgradRev, ProbDistributionsStudentT_cdf_lccdf_near_zero) { + // q >= 2 with |t| small. The complement 1 - I_r(1/2, nu/2) and its nu + // gradient must be evaluated from r, not from 1 - r, which is close to 1 + // here: from 1 - r the nu gradient had a relative error of 1e-3 at + // t = 1e-7. nu = 1 is Boost's arcsine case a = b = 1/2. + expect_student_t_fn(student_t_cdf_fn, -1e-7, 1.0, 1.0, + {0.49999996816901138, 0.31830988618378747, + -6.1480657060756583e-09, 3.1830988618378744e-08}, + 1e-11); + expect_student_t_fn(student_t_lccdf_fn, 1e-3, 1.0, 1.0, + {-0.69378400284838371, -0.63702467811884134, + -0.00012303970872281812, 0.00063702467811884139}, + 1e-11); + expect_student_t_fn(student_t_cdf_fn, -1e-5, 30.0, 1.0, + {0.49999604367815115, 0.39563218487365676, + -1.0983690983393141e-09, 3.9563218487365682e-06}, + 1e-11); +} + +TEST_F(AgradRev, ProbDistributionsStudentT_lccdf_branch_seam) { + // the two branches meet at q = 2; the value and the gradients must be + // continuous across it + using stan::math::var; + double v[3], d[3]; + int k = 0; + for (double y : {1.9999, 2.0, 2.0001}) { + var nu = 8; + var lp = stan::math::student_t_lccdf(y, nu, 0.0, 1.0); + lp.grad(); + v[k] = lp.val(); + d[k] = nu.adj(); + ++k; + } + EXPECT_NEAR(v[1], 0.5 * (v[0] + v[2]), 1e-8); + EXPECT_NEAR(d[1], 0.5 * (d[0] + d[2]), 1e-8); +} + +TEST_F(AgradRev, ProbDistributionsStudentT_lcdf_cdf_deep_tail) { + using stan::math::var; + { + var nu = 1000; + var sigma = 1; + var lp = stan::math::student_t_lcdf(-10.0, nu, 0.0, sigma); + lp.grad(); + EXPECT_TRUE(std::isfinite(lp.val())); + EXPECT_NEAR(-50.838951454, lp.val(), 1e-9 * 50.9); + EXPECT_NEAR(-0.0022457674682, nu.adj(), 1e-9 * 0.00225); + } + { + // the cdf itself: the value is representable and the gradients were NaN + var nu = 1000; + var sigma = 1; + var p = stan::math::student_t_cdf(-10.0, nu, 0.0, sigma); + p.grad(); + EXPECT_TRUE(std::isfinite(nu.adj()) && std::isfinite(sigma.adj())); + EXPECT_NEAR(8.33535e-23, p.val(), 1e-5 * 8.3e-23); + } +}