From b172d6165358be8556fd10380c3e870f4271d583 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sat, 19 Sep 2026 00:11:25 +0300 Subject: [PATCH 1/9] Use Eigen igamma_der_a for the incomplete gamma gradient roots --- stan/math/fwd/fun.hpp | 1 + stan/math/fwd/fun/Eigen_SpecialFunctions.hpp | 65 +++++++ stan/math/prim/fun/grad_reg_inc_gamma.hpp | 111 +++-------- .../prim/fun/grad_reg_lower_inc_gamma.hpp | 174 +++--------------- stan/math/rev/fun.hpp | 1 + stan/math/rev/fun/Eigen_SpecialFunctions.hpp | 62 +++++++ .../math/mix/fun/grad_reg_inc_gamma_test.cpp | 67 ++++++- .../fun/grad_reg_inc_gamma_accuracy_test.cpp | 106 +++++++++++ 8 files changed, 349 insertions(+), 238 deletions(-) create mode 100644 stan/math/fwd/fun/Eigen_SpecialFunctions.hpp create mode 100644 stan/math/rev/fun/Eigen_SpecialFunctions.hpp create mode 100644 test/unit/math/prim/fun/grad_reg_inc_gamma_accuracy_test.cpp diff --git a/stan/math/fwd/fun.hpp b/stan/math/fwd/fun.hpp index 21c6d70eee0..7a5c2ba86fc 100644 --- a/stan/math/fwd/fun.hpp +++ b/stan/math/fwd/fun.hpp @@ -3,6 +3,7 @@ #include #include +#include #include #include diff --git a/stan/math/fwd/fun/Eigen_SpecialFunctions.hpp b/stan/math/fwd/fun/Eigen_SpecialFunctions.hpp new file mode 100644 index 00000000000..73a41c05159 --- /dev/null +++ b/stan/math/fwd/fun/Eigen_SpecialFunctions.hpp @@ -0,0 +1,65 @@ +#ifndef STAN_MATH_FWD_FUN_EIGEN_SPECIALFUNCTIONS_HPP +#define STAN_MATH_FWD_FUN_EIGEN_SPECIALFUNCTIONS_HPP + +#include +#include +#include +#include +#include + +namespace Eigen { +namespace internal { + +/** + * Support for Eigen's incomplete gamma routines with Stan's forward-mode + * scalars. + * + * `Eigen::internal::igamma_generic_impl` is written for a general scalar + * type and propagates the derivative through the Cephes power series and + * the Cephes continued fraction. Only three of its helpers are restricted + * to `float` and `double`. Supplying those three here lets + * `Eigen::numext::igamma_der_a` serve `fvar` at every autodiff order, + * so `grad_reg_inc_gamma` and `grad_reg_lower_inc_gamma` can use one + * algorithm for every scalar type. + * + * This depends on the names `lgamma_impl`, `digamma_impl` and + * `cephes_helper` in `Eigen::internal`. An Eigen upgrade that renames them + * fails at compile time, not silently. + */ +template +struct lgamma_impl> { + EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE stan::math::fvar run( + const stan::math::fvar& x) { + return stan::math::lgamma(x); + } +}; + +template +struct digamma_impl> { + EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE stan::math::fvar run( + const stan::math::fvar& x) { + return stan::math::digamma(x); + } +}; + +/** + * The Cephes constants are properties of the underlying floating point + * format, not of the autodiff type, so they carry a zero tangent. + */ +template +struct cephes_helper> { + EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE stan::math::fvar machep() { + return stan::math::fvar(cephes_helper::machep()); + } + EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE stan::math::fvar big() { + return stan::math::fvar(cephes_helper::big()); + } + EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE stan::math::fvar biginv() { + return stan::math::fvar(cephes_helper::biginv()); + } +}; + +} // namespace internal +} // namespace Eigen + +#endif diff --git a/stan/math/prim/fun/grad_reg_inc_gamma.hpp b/stan/math/prim/fun/grad_reg_inc_gamma.hpp index 79f47def83e..4aa6441f18b 100644 --- a/stan/math/prim/fun/grad_reg_inc_gamma.hpp +++ b/stan/math/prim/fun/grad_reg_inc_gamma.hpp @@ -2,115 +2,56 @@ #define STAN_MATH_PRIM_FUN_GRAD_REG_INC_GAMMA_HPP #include -#include -#include -#include -#include -#include +#include #include -#include -#include -#include -#include +#include #include namespace stan { namespace math { /** - * Gradient of the regularized incomplete gamma functions igamma(a, z) + * Gradient of the regularized incomplete gamma function igamma(a, z) with + * respect to the shape parameter, d/da Q(a, z), where Q is the upper + * regularized incomplete gamma function. * - * For small z, the gradient is computed via the series expansion; - * for large z, the series is numerically inaccurate due to cancellation - * and the asymptotic expansion is used. + * The work is done by `Eigen::numext::igamma_der_a`, which differentiates + * the Cephes power series and the Cephes continued fraction term by term + * and stops at machine epsilon. It returns d/da of the lower P(a, z), so + * the sign is flipped here. + * + * Eigen's implementation is generic in the scalar type. Stan supplies the + * three Cephes helpers it needs for autodiff scalars in + * `stan/math/fwd/fun/Eigen_SpecialFunctions.hpp` and + * `stan/math/rev/core/Eigen_SpecialFunctions.hpp`, so every autodiff order + * uses this same algorithm. + * + * `g`, `dig`, `precision` and `max_steps` are accepted and ignored. They + * belonged to the previous hand-written series. The signature is unchanged + * so that the 27 call sites do not change; removing the arguments is + * proposed separately. * * @tparam T1 type of the shape parameter * @tparam T2 type of the location parameter * @param a shape parameter, a > 0 * @param z location z >= 0 - * @param g stan::math::tgamma(a) (precomputed value) - * @param dig boost::math::digamma(a) (precomputed value) - * @param precision required precision; applies to series expansion only - * @param max_steps number of steps to take. - * @throw throws std::domain_error if not converged after max_steps - * or increment overflows to inf. - * - * For the asymptotic expansion, the gradient is given by: - \f[ - \begin{array}{rcl} - \Gamma(a, z) & = & z^{a-1}e^{-z} \sum_{k=0}^N \frac{(a-1)_k}{z^k} \qquad , z - \gg a\\ - Q(a, z) & = & \frac{z^{a-1}e^{-z}}{\Gamma(a)} \sum_{k=0}^N - \frac{(a-1)_k}{z^k}\\ - (a)_k & = & (a)_{k-1}(a-k)\\ - \frac{d}{da} (a)_k & = & (a)_{k-1} + (a-k)\frac{d}{da} (a)_{k-1}\\ - \frac{d}{da}Q(a, z) & = & (log(z) - \psi(a)) Q(a, z)\\ - && + \frac{z^{a-1}e^{-z}}{\Gamma(a)} \sum_{k=0}^N \left(\frac{d}{da} - (a-1)_k\right) \frac{1}{z^k} \end{array} \f] + * @param g ignored; previously stan::math::tgamma(a) + * @param dig ignored; previously stan::math::digamma(a) + * @param precision ignored; previously the series tolerance + * @param max_steps ignored; previously the series iteration limit + * @return d/da of the upper regularized incomplete gamma function */ template inline return_type_t grad_reg_inc_gamma(T1 a, T2 z, T1 g, T1 dig, double precision = 1e-6, int max_steps = 1e5) { - using std::exp; - using std::fabs; - using std::log; using TP = return_type_t; if (is_any_nan(a, z, g, dig)) { return std::numeric_limits::quiet_NaN(); } - T2 l = log(z); - if (z >= a && z >= 8) { - // asymptotic expansion http://dlmf.nist.gov/8.11#E2 - TP S = 0; - T1 a_minus_one_minus_k = a - 1; - T1 fac = a_minus_one_minus_k; // falling_factorial(a-1, k) - T1 dfac = 1; // d/da[falling_factorial(a-1, k)] - T2 zpow = z; // z ** k - TP delta = dfac / zpow; - - for (int k = 1; k < 10; ++k) { - a_minus_one_minus_k -= 1; - - S += delta; - - zpow *= z; - dfac = a_minus_one_minus_k * dfac + fac; - fac *= a_minus_one_minus_k; - delta = dfac / zpow; - - if (is_inf(delta)) { - throw_domain_error("grad_reg_inc_gamma", "is not converging", "", ""); - } - } - - return gamma_q(a, z) * (l - dig) + exp(-z + (a - 1) * l) * S / g; - } else { - // gradient of series expansion http://dlmf.nist.gov/8.7#E3 - TP S = 0; - TP log_s = 0.0; - double s_sign = 1.0; - T2 log_z = log(z); - TP log_delta = log_s - multiply_log(2, a); - for (int k = 1; k <= max_steps; ++k) { - S += s_sign >= 0.0 ? exp(log_delta) : -exp(log_delta); - log_s += log_z - log(k); - s_sign = -s_sign; - log_delta = log_s - multiply_log(2, k + a); - if (is_inf(log_delta)) { - throw_domain_error("grad_reg_inc_gamma", "is not converging", "", ""); - } - if (log_delta <= log(precision)) { - return gamma_p(a, z) * (dig - l) + exp(a * l) * S / g; - } - } - throw_domain_error( - "grad_reg_inc_gamma", "k (internal counter)", max_steps, "exceeded ", - " iterations, gamma function gradient did not converge."); - return INFTY; - } + return -Eigen::numext::igamma_der_a(TP(a), TP(z)); } } // namespace math diff --git a/stan/math/prim/fun/grad_reg_lower_inc_gamma.hpp b/stan/math/prim/fun/grad_reg_lower_inc_gamma.hpp index fb6c3620904..7f154f5eed1 100644 --- a/stan/math/prim/fun/grad_reg_lower_inc_gamma.hpp +++ b/stan/math/prim/fun/grad_reg_lower_inc_gamma.hpp @@ -3,118 +3,49 @@ #include #include -#include -#include -#include -#include +#include #include -#include -#include -#include -#include -#include -#include -#include +#include #include -#include namespace stan { namespace math { /** - * Computes the gradient of the lower regularized incomplete - * gamma function. + * Computes the gradient of the lower regularized incomplete gamma + * function, d/da P(a, z). * - * The lower incomplete gamma function - * derivative w.r.t its first parameter (a) seems to have no - * standard source. It also appears to have no widely known - * approximate implementation. Gautschi (1979) has a thorough - * discussion of the calculation of the lower regularized - * incomplete gamma function itself and some stability issues. + * The work is done by `Eigen::numext::igamma_der_a`, which returns exactly + * this quantity. It differentiates the Cephes power series and the Cephes + * continued fraction term by term and stops at machine epsilon. Eigen's + * implementation is generic in the scalar type; Stan supplies the three + * Cephes helpers it needs for autodiff scalars in + * `stan/math/fwd/fun/Eigen_SpecialFunctions.hpp` and + * `stan/math/rev/core/Eigen_SpecialFunctions.hpp`. * - * Reference: Gautschi, Walter (1979) ACM Transactions on - * mathematical software. 5(4):466-481 + * The previous implementation used a Gautschi series whose two sums cancel + * completely for z well above a. It returned about 1e-14 of rounding noise, + * with arbitrary sign, where the true value can be 1e-148. * - * We implemented calculations for d(gamma_p)/da by taking - * derivatives of formulas suggested by Gauschi and others and - * testing them against an outside source (Mathematica). We - * took three implementations which can cover the range {a:[0,20], - * z:[0,30]} with absolute error < 1e-10 with the exception of - * values near (0,0) where the error is near 1e-5. Relative error - * is also <<1e-6 except for regions where the gradient approaches - * zero. + * `precision` and `max_steps` are accepted and ignored. They belonged to + * the previous series. The signature is unchanged so that the call sites do + * not change; removing the arguments is proposed separately. * - * Gautschi suggests calculating the lower incomplete gamma - * function for small to moderate values of $z$ using the - * approximation: - * - * \f[ - * \frac{\gamma(a,z)}{\Gamma(a)}=z^a e^-z - * \sum_n=0^\infty \frac{z^n}{\Gamma(a+n+1)} - * \f] - * - * We write the derivative in the form: - * - * \f[ - * \frac{d\gamma(a,z)\Gamma(a)}{da} = \frac{\log z}{e^z} - * \sum_n=0^\infty \frac{z^{a+n}}{\Gamma(a+n+1)} - * - \frac{1}{e^z} - * \sum_n=0^\infty \frac{z^{a+n}}{\Gamma(a+n+1)}\psi^0(a+n+1) - * \f] - * - * This calculation is sufficiently accurate for small $a$ and - * small $z$. For larger values and $a$ and $z$ we use it in its - * log form: - * - * \f[ - * \frac{d \gamma(a,z)\Gamma(a)}{da} = \frac{\log z}{e^z} - * \sum_n=0^\infty \exp[(a+n)\log z - \log\Gamma(a+n+1)] - * - \sum_n=0^\infty \exp[(a+n)\log z - \log\Gamma(a+n+1) + - * \log\psi^0(a+n+1)] - * \f] - * - * For large $z$, Gauschi recommends using the upper incomplete - * Gamma instead and the negative of its derivative turns out to be - * more stable and accurate for larger $z$ and for some combinations - * of $a$ and $z$. This is a log-scale implementation of the - * derivative of the formulation suggested by Gauschi (1979). For - * some values it defers to the negative of the gradient - * for the gamma_q function. This is based on the suggestion by Gauschi - * (1979) that for large values of $z$ it is better to - * carry out calculations using the upper incomplete Gamma function. - * - * Branching for choice of implementation for the lower incomplete - * regularized gamma function gradient. The derivative based on - * Gautschi's formulation appears to be sufficiently accurate - * everywhere except for large z and small to moderate a. The - * intersection between the two regions is a radius 12 quarter circle - * centered at a=0, z=30 although both implementations are - * satisfactory near the intersection. - * - * Some limits that could be treated, e.g., infinite z should - * return tgamma(a) * digamma(a), throw instead to match the behavior of, - * e.g., boost::math::gamma_p + * Infinite z still throws, to match the behaviour of, for example, + * boost::math::gamma_p. * * @tparam T1 type of a * @tparam T2 type of z - * @param[in] a shared with complete Gamma - * @param[in] z value to integrate up to - * @param[in] precision series terminates when increment falls below - * this value. - * @param[in] max_steps number of terms to sum before throwing - * @throw std::domain_error if the series does not converge to - * requested precision before max_steps. - * + * @param[in] a shared with complete Gamma, a > 0 + * @param[in] z value to integrate up to, z >= 0 + * @param[in] precision ignored; previously the series tolerance + * @param[in] max_steps ignored; previously the series iteration limit + * @return d/da of the lower regularized incomplete gamma function */ template inline return_type_t grad_reg_lower_inc_gamma(const T1& a, const T2& z, double precision = 1e-10, int max_steps = 1e5) { - using std::exp; - using std::log; - using std::pow; - using std::sqrt; - using TP = return_type_t; if (is_any_nan(a, z)) { @@ -128,62 +59,7 @@ inline return_type_t grad_reg_lower_inc_gamma(const T1& a, const T2& z, } check_positive_finite("grad_reg_lower_inc_gamma", "z", z); - if ((a < 0.8 && z > 15.0) || (a < 12.0 && z > 30.0) - || a < sqrt(-756 - value_of_rec(z) * value_of_rec(z) - + 60 * value_of_rec(z))) { - T1 tg = tgamma(a); - T1 dig = digamma(a); - return -grad_reg_inc_gamma(a, z, tg, dig, precision, max_steps); - } - - T2 log_z = log(z); - T2 emz = exp(-z); - - int n = 0; - T1 a_plus_n = a; - TP sum_a = 0.0; - T1 lgamma_a_plus_1 = lgamma(a + 1); - T1 lgamma_a_plus_n_plus_1 = lgamma_a_plus_1; - TP term; - while (true) { - term = exp(a_plus_n * log_z - lgamma_a_plus_n_plus_1); - sum_a += term; - if (term <= precision) { - break; - } - if (n >= max_steps) { - throw_domain_error("grad_reg_lower_inc_gamma", "n (internal counter)", - max_steps, "exceeded ", - " iterations, gamma_p(a,z) gradient (a) " - "did not converge."); - } - ++n; - lgamma_a_plus_n_plus_1 += log1p(a_plus_n); - ++a_plus_n; - } - - n = 1; - a_plus_n = a + 1; - TP digammap1 = digamma(a_plus_n); - TP sum_b = digammap1 * exp(a * log_z - lgamma_a_plus_1); - lgamma_a_plus_n_plus_1 = lgamma_a_plus_1 + log(a_plus_n); - while (true) { - digammap1 += 1 / a_plus_n; - term = exp(a_plus_n * log_z - lgamma_a_plus_n_plus_1) * digammap1; - sum_b += term; - if (term <= precision) { - return emz * (log_z * sum_a - sum_b); - } - if (n >= max_steps) { - throw_domain_error("grad_reg_lower_inc_gamma", "n (internal counter)", - max_steps, "exceeded ", - " iterations, gamma_p(a,z) gradient (a) " - "did not converge."); - } - ++n; - lgamma_a_plus_n_plus_1 += log1p(a_plus_n); - ++a_plus_n; - } + return Eigen::numext::igamma_der_a(TP(a), TP(z)); } } // namespace math diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp index b9a8bfa3e2c..5aeea8212bb 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -3,6 +3,7 @@ #include +#include #include #include #include diff --git a/stan/math/rev/fun/Eigen_SpecialFunctions.hpp b/stan/math/rev/fun/Eigen_SpecialFunctions.hpp new file mode 100644 index 00000000000..458104c7933 --- /dev/null +++ b/stan/math/rev/fun/Eigen_SpecialFunctions.hpp @@ -0,0 +1,62 @@ +#ifndef STAN_MATH_REV_FUN_EIGEN_SPECIALFUNCTIONS_HPP +#define STAN_MATH_REV_FUN_EIGEN_SPECIALFUNCTIONS_HPP + +#include +#include +#include +#include +#include + +namespace Eigen { +namespace internal { + +/** + * Support for Eigen's incomplete gamma routines with Stan's reverse-mode + * scalar. See `stan/math/fwd/fun/Eigen_SpecialFunctions.hpp` for the + * reasoning; this is the same three helpers for `var`. + * + * No function inside Stan Math reaches `igamma_der_a` with `var`. Every + * distribution calls the gradient roots with `T_partials_return`, which is + * `double` for reverse mode. These specializations exist so that a direct + * call such as `grad_reg_inc_gamma(var, var, var, var)` keeps working after + * the hand-written series was removed. + * + * A `var` call builds one autodiff node per operation inside the series, so + * it is far more expensive than the `double` call. That was equally true of + * the previous implementation. + */ +template <> +struct lgamma_impl { + static EIGEN_STRONG_INLINE stan::math::var run(const stan::math::var& x) { + return stan::math::lgamma(x); + } +}; + +template <> +struct digamma_impl { + static EIGEN_STRONG_INLINE stan::math::var run(const stan::math::var& x) { + return stan::math::digamma(x); + } +}; + +/** + * The Cephes constants are properties of the underlying floating point + * format, so they are constants with no adjoint. + */ +template <> +struct cephes_helper { + static EIGEN_STRONG_INLINE stan::math::var machep() { + return stan::math::var(cephes_helper::machep()); + } + static EIGEN_STRONG_INLINE stan::math::var big() { + return stan::math::var(cephes_helper::big()); + } + static EIGEN_STRONG_INLINE stan::math::var biginv() { + return stan::math::var(cephes_helper::biginv()); + } +}; + +} // namespace internal +} // namespace Eigen + +#endif diff --git a/test/unit/math/mix/fun/grad_reg_inc_gamma_test.cpp b/test/unit/math/mix/fun/grad_reg_inc_gamma_test.cpp index 87bd0b1240a..af11db17f01 100644 --- a/test/unit/math/mix/fun/grad_reg_inc_gamma_test.cpp +++ b/test/unit/math/mix/fun/grad_reg_inc_gamma_test.cpp @@ -8,7 +8,14 @@ TEST(ProbInternalMath, gradRegIncGamma_typical) { double g = 1.77245; double dig = -1.96351; - EXPECT_FLOAT_EQ(0.38984156, stan::math::grad_reg_inc_gamma(a, b, g, dig)); + // The expected value was 0.38984156, which is 1.1e-5 away from the true + // d/da Q(0.5, 1.0). That constant encoded two errors of the old + // implementation: the rounded support values that this test passes in + // (g = 1.77245 is 2.2e-6 away from tgamma(0.5)), and the 1e-6 series + // truncation. The arithmetic path no longer uses g or dig, so the result + // is now the true value, from mpmath at 80 digits. + EXPECT_FLOAT_EQ(0.38983726432851057, + stan::math::grad_reg_inc_gamma(a, b, g, dig)); } TEST(ProbInternalMath, gradRegIncGamma_infLoopInVersion2_0_1) { @@ -30,6 +37,8 @@ TEST(ProbInternalMath, gradRegIncGamma_largeZ) { stan::math::grad_reg_inc_gamma(a, z, g, dig)); } +// Every scalar type now uses the same algorithm, so the fvar cases expect +// the same value as the arithmetic case. TEST(ProbInternalMath, gradRegIncGamma_fd) { using stan::math::fvar; @@ -38,7 +47,7 @@ TEST(ProbInternalMath, gradRegIncGamma_fd) { fvar g = 1.77245; fvar dig = -1.96351; - EXPECT_FLOAT_EQ(0.38984156, + EXPECT_FLOAT_EQ(0.38983726432851057, stan::math::grad_reg_inc_gamma(a, b, g, dig).val()); } TEST(ProbInternalMath, gradRegIncGamma_ffd) { @@ -49,7 +58,7 @@ TEST(ProbInternalMath, gradRegIncGamma_ffd) { fvar > g = 1.77245; fvar > dig = -1.96351; - EXPECT_FLOAT_EQ(0.38984156, + EXPECT_FLOAT_EQ(0.38983726432851057, stan::math::grad_reg_inc_gamma(a, b, g, dig).val_.val_); } @@ -63,7 +72,7 @@ TEST(ProbInternalMath, gradRegIncGamma_fv) { fvar g = 1.77245; fvar dig = digamma(a); - EXPECT_FLOAT_EQ(0.38984156, + EXPECT_FLOAT_EQ(0.38983726432851057, stan::math::grad_reg_inc_gamma(a, b, g, dig).val_.val()); } @@ -108,3 +117,53 @@ TEST(ProbInternalMath, gradRegIncGamma_fv_2ndderiv) { z.d_.grad(y1, grad1); EXPECT_NEAR(-0.546236927878295422, grad1[0], 1e-6); } + +namespace { +struct gamma_cdf_shape_functor { + double y_; + explicit gamma_cdf_shape_functor(double y) : y_(y) {} + template + T operator()(const Eigen::Matrix& x) const { + return stan::math::gamma_cdf(y_, x(0), 1.0); + } +}; +} // namespace + +/** + * The first derivative must not depend on how it is obtained. + * + * Reverse mode reaches the gradient root with `double` partials; `hessian` + * reaches it with `fvar`. While those were two different algorithms, + * the two routes disagreed by 4.4e-02 at alpha = 20, 5.8e-01 at alpha = 100 + * and 7.2e-01 at alpha = 171, and `hessian` returned NaN at alpha = 180, + * all with `z = alpha`. This test locks that down. It asserts agreement + * between the two routes, so it needs no external reference value. + */ +TEST(ProbInternalMath, gradRegIncGamma_gradient_matches_hessian) { + using stan::math::var; + + for (double alpha : {5.0, 20.0, 100.0, 171.0, 180.0, 300.0}) { + const double y = alpha; + + double rev_grad; + { + stan::math::nested_rev_autodiff nested; + var a = alpha; + var p = stan::math::gamma_cdf(y, a, 1.0); + p.grad(); + rev_grad = a.adj(); + } + + Eigen::Matrix x(1); + x(0) = alpha; + double fx; + Eigen::Matrix grad(1); + Eigen::Matrix hess(1, 1); + stan::math::hessian(gamma_cdf_shape_functor(y), x, fx, grad, hess); + + ASSERT_TRUE(std::isfinite(rev_grad)) << "reverse mode, alpha = " << alpha; + ASSERT_TRUE(std::isfinite(grad(0))) << "hessian(), alpha = " << alpha; + EXPECT_NEAR(rev_grad, grad(0), 1e-12 * std::fabs(rev_grad)) + << "alpha = " << alpha; + } +} diff --git a/test/unit/math/prim/fun/grad_reg_inc_gamma_accuracy_test.cpp b/test/unit/math/prim/fun/grad_reg_inc_gamma_accuracy_test.cpp new file mode 100644 index 00000000000..346e05c6b59 --- /dev/null +++ b/test/unit/math/prim/fun/grad_reg_inc_gamma_accuracy_test.cpp @@ -0,0 +1,106 @@ +#include +#include +#include +#include + +/** + * Fixed-reference accuracy tests for the incomplete gamma gradient roots. + * + * Every expected value comes from mpmath at 80 digits, through two + * independent routes that agree to at least 13 digits: numerical + * differentiation of the regularized function, and an explicit integral + * identity evaluated in density form. A third route, the Gautschi series, + * confirms the deep-tail values. The generator is + * `inc_gamma_beta_work/ref_mpmath.py` in the audit workspace. + * + * The references reproduce two constants that were already in + * grad_reg_inc_gamma_test.cpp, which is the check that the two agree: + * (a, z) = (1.1, 0.2) gives 0.31416364892410886, and (9, 10) gives + * 0.120855166827777. + * + * These tests fail on develop at 5252d51d47. They exist because + * `expect_ad` compares against finite differences and cannot see an error + * of this size, and because the existing suites assert absolute tolerances + * on a grid that stops short of every failing region. + */ + +namespace { + +/** Relative comparison, needed because the values span 1e-1 to 1e-148. */ +void expect_rel(double expected, double actual, double tol, + const std::string& label) { + ASSERT_TRUE(std::isfinite(actual)) + << label << ": result is not finite, got " << actual; + const double err = std::fabs(actual - expected) / std::fabs(expected); + EXPECT_LT(err, tol) << label << ": relative error " << err << ", expected " + << expected << ", got " << actual; +} + +double upper(double a, double z) { + return stan::math::grad_reg_inc_gamma(a, z, stan::math::tgamma(a), + stan::math::digamma(a)); +} + +} // namespace + +/** + * D1. `z >= a && z >= 8` selects a fixed 10-term asymptotic expansion that + * needs z >> a. On develop the error is 4.4e-2 at a = 20, 5.8e-1 at + * a = 100 and 7.2e-1 at a = 171. + */ +TEST(MathPrimScalFun, grad_reg_inc_gamma_diagonal_z_equals_a) { + expect_rel(8.95792055391225334e-02, upper(20.0, 20.0), 1e-10, "a=z=20"); + expect_rel(3.99274978577860498e-02, upper(100.0, 100.0), 1e-10, "a=z=100"); + expect_rel(3.05227525380366363e-02, upper(171.0, 171.0), 1e-10, "a=z=171"); +} + +/** + * D2. Above a = 171.62, tgamma(a) is inf, so the correction term is + * inf/inf. On develop every one of these returns NaN. + */ +TEST(MathPrimScalFun, grad_reg_inc_gamma_above_tgamma_overflow) { + expect_rel(2.97491740225716214e-02, upper(180.0, 180.0), 1e-10, "a=z=180"); + expect_rel(2.30393429308686434e-02, upper(300.0, 300.0), 1e-10, "a=z=300"); + expect_rel(1.78442151466593514e-02, upper(500.0, 500.0), 1e-10, "a=z=500"); +} + +/** + * D3. `grad_reg_lower_inc_gamma` runs its Gautschi branch for every z > 42 + * with a >= 12, and there `emz * (log_z * sum_a - sum_b)` cancels + * completely. On develop it returns about 1e-14 of noise, with arbitrary + * sign, at all four points. + */ +TEST(MathPrimScalFun, grad_reg_lower_inc_gamma_far_upper_tail) { + expect_rel(-8.28941360740501778e-22, + stan::math::grad_reg_lower_inc_gamma(50.0, 150.0), 1e-10, + "a=50 z=150"); + expect_rel(-2.79418207415687613e-54, + stan::math::grad_reg_lower_inc_gamma(50.0, 250.0), 1e-10, + "a=50 z=250"); + expect_rel(-5.33817508735222058e-148, + stan::math::grad_reg_lower_inc_gamma(50.0, 500.0), 1e-10, + "a=50 z=500"); + expect_rel(-4.74445289261568464e-69, + stan::math::grad_reg_lower_inc_gamma(170.0, 510.0), 1e-10, + "a=170 z=510"); +} + +/** + * The two functions compute the same quantity with opposite sign. On + * develop they disagree by 72 % at a = z = 171. This test needs no external + * reference at all. + */ +TEST(MathPrimScalFun, grad_reg_inc_gamma_agrees_with_lower) { + for (double a : {5.0, 20.0, 50.0, 100.0, 171.0}) { + const double up = upper(a, a); + const double lo = stan::math::grad_reg_lower_inc_gamma(a, a); + expect_rel(-up, lo, 1e-10, "a=z=" + std::to_string(a)); + } +} + +/** Points that already worked, to show the change is not a regression. */ +TEST(MathPrimScalFun, grad_reg_inc_gamma_unchanged_points) { + expect_rel(3.14163648924108863e-01, upper(1.1, 0.2), 1e-10, "a=1.1 z=0.2"); + expect_rel(1.20855166827777014e-01, upper(9.0, 10.0), 1e-10, "a=9 z=10"); + expect_rel(3.32046796768609819e-11, upper(2.5, 30.0), 1e-10, "a=2.5 z=30"); +} From af3205c5acee2d34d66cb3ae15a94ab6df6a857e Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sat, 19 Sep 2026 19:54:14 +0300 Subject: [PATCH 2/9] Remove the duplicated incomplete gamma series from fwd gamma_q --- stan/math/fwd/fun/gamma_q.hpp | 52 ++++++++++++----------------------- 1 file changed, 18 insertions(+), 34 deletions(-) diff --git a/stan/math/fwd/fun/gamma_q.hpp b/stan/math/fwd/fun/gamma_q.hpp index 2e1951bcd82..9c03670ae03 100644 --- a/stan/math/fwd/fun/gamma_q.hpp +++ b/stan/math/fwd/fun/gamma_q.hpp @@ -5,37 +5,36 @@ #include #include #include -#include -#include #include #include #include +#include #include namespace stan { namespace math { +/* + * The derivative with respect to the shape parameter used to be an inlined + * copy of the series in `grad_reg_inc_gamma`, with the same hard-coded 1e-6 + * tolerance and without that function's second branch. The copy is now + * replaced by a call to the root itself. + * + * Keeping the copy hid the defect it shared with the root. In + * `test/prob/chi_square`, the generated `ffv` case compares the + * distribution's analytic partials against autodiff through + * `log(gamma_q(nu * 0.5, y * 0.5))`. Both routes used the same inaccurate + * series, so their errors cancelled and the comparison passed. Once the + * root became accurate and the copy did not, the same comparison failed by + * 2.4e-03 at third order. + */ + template inline fvar gamma_q(const fvar& x1, const fvar& x2) { T u = gamma_q(x1.val_, x2.val_); - T S = 0; - T s = 1; - T l = log(x2.val_); T g = tgamma(x1.val_); - T dig = digamma(x1.val_); - - int k = 0; - T delta = s / (x1.val_ * x1.val_); - - while (fabs(delta) > 1e-6) { - S += delta; - ++k; - s *= -x2.val_ / k; - delta = s / ((k + x1.val_) * (k + x1.val_)); - } - - T der1 = (1.0 - u) * (dig - l) + exp(x1.val_ * l) * S / g; + T der1 = grad_reg_inc_gamma(x1.val_, x2.val_, g, digamma(x1.val_)); T der2 = -exp(-x2.val_) * pow(x2.val_, x1.val_ - 1.0) / g; return fvar(u, x1.d_ * der1 + x2.d_ * der2); @@ -44,23 +43,8 @@ inline fvar gamma_q(const fvar& x1, const fvar& x2) { template inline fvar gamma_q(const fvar& x1, double x2) { T u = gamma_q(x1.val_, x2); - T S = 0; - double s = 1; - double l = log(x2); - T g = tgamma(x1.val_); - T dig = digamma(x1.val_); - - int k = 0; - T delta = s / (x1.val_ * x1.val_); - - while (fabs(delta) > 1e-6) { - S += delta; - ++k; - s *= -x2 / k; - delta = s / ((k + x1.val_) * (k + x1.val_)); - } - T der1 = (1.0 - u) * (dig - l) + exp(x1.val_ * l) * S / g; + T der1 = grad_reg_inc_gamma(x1.val_, x2, tgamma(x1.val_), digamma(x1.val_)); return fvar(u, x1.d_ * der1); } From f23eaea890b16900b7daf1159a563abaae4cbb74 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sun, 20 Sep 2026 11:16:50 +0300 Subject: [PATCH 3/9] Evaluate the fwd gamma_q second-argument derivative in log space --- stan/math/fwd/fun/gamma_q.hpp | 22 ++++++++++++++++------ 1 file changed, 16 insertions(+), 6 deletions(-) diff --git a/stan/math/fwd/fun/gamma_q.hpp b/stan/math/fwd/fun/gamma_q.hpp index 9c03670ae03..bf276ebac22 100644 --- a/stan/math/fwd/fun/gamma_q.hpp +++ b/stan/math/fwd/fun/gamma_q.hpp @@ -5,7 +5,8 @@ #include #include #include -#include +#include +#include #include #include #include @@ -15,6 +16,8 @@ namespace stan { namespace math { /* + * Two changes from the previous version of this file. + * * The derivative with respect to the shape parameter used to be an inlined * copy of the series in `grad_reg_inc_gamma`, with the same hard-coded 1e-6 * tolerance and without that function's second branch. The copy is now @@ -27,15 +30,21 @@ namespace math { * series, so their errors cancelled and the comparison passed. Once the * root became accurate and the copy did not, the same comparison failed by * 2.4e-03 at third order. + * + * The derivative with respect to the second argument used to be + * `-exp(-x2) * pow(x2, x1 - 1) / tgamma(x1)`. Both `pow` and `tgamma` + * overflow for moderate arguments, and the quotient then silently returns + * -0 rather than a small negative number. It is now evaluated in log space, + * which is the form `gamma_p` already uses for the same quantity. */ template inline fvar gamma_q(const fvar& x1, const fvar& x2) { T u = gamma_q(x1.val_, x2.val_); - T g = tgamma(x1.val_); - T der1 = grad_reg_inc_gamma(x1.val_, x2.val_, g, digamma(x1.val_)); - T der2 = -exp(-x2.val_) * pow(x2.val_, x1.val_ - 1.0) / g; + T der1 + = grad_reg_inc_gamma(x1.val_, x2.val_, tgamma(x1.val_), digamma(x1.val_)); + T der2 = -exp(-x2.val_ + (x1.val_ - 1.0) * log(x2.val_) - lgamma(x1.val_)); return fvar(u, x1.d_ * der1 + x2.d_ * der2); } @@ -52,8 +61,9 @@ inline fvar gamma_q(const fvar& x1, double x2) { template inline fvar gamma_q(double x1, const fvar& x2) { T u = gamma_q(x1, x2.val_); - double g = tgamma(x1); - T der2 = -exp(-x2.val_) * pow(x2.val_, x1 - 1.0) / g; + + T der2 = -exp(-x2.val_ + (x1 - 1.0) * log(x2.val_) - lgamma(x1)); + return fvar(u, x2.d_ * der2); } } // namespace math From 1a9df72569d6c4f5f3cafcd03a120fb171e28a6d Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sun, 20 Sep 2026 15:44:02 +0300 Subject: [PATCH 4/9] Test the forward-mode gamma_q derivatives at large z and large a --- test/unit/math/mix/fun/gamma_q_test.cpp | 56 +++++++++++++++++++++++++ 1 file changed, 56 insertions(+) diff --git a/test/unit/math/mix/fun/gamma_q_test.cpp b/test/unit/math/mix/fun/gamma_q_test.cpp index a14705608ad..0ef0bef08c2 100644 --- a/test/unit/math/mix/fun/gamma_q_test.cpp +++ b/test/unit/math/mix/fun/gamma_q_test.cpp @@ -1,4 +1,5 @@ #include +#include #include TEST(mathMixScalFun, gammaQ) { @@ -18,6 +19,61 @@ TEST(mathMixScalFun, gammaQ) { // stan::test::expect_value(f, 8.01006, 2.47579e+215); } +/** + * Forward-mode shape derivative at large z, fixed references from mpmath + * at 80 digits. + * + * The previous forward-mode gamma_q carried its own copy of the series in + * grad_reg_inc_gamma, without an iteration cap. For z above about 700 the + * alternating terms overflow before they decay and the loop never + * terminates; that is the failure the commented-out expect_value above + * ran into. Below that, the series terminates and returns garbage: at + * (a, z) = (1.5, 300) it gave 1.2e+112 where the derivative is 5.7e-129. + * + * These points are chosen so that the old code fails rather than hangs, + * which keeps the test usable in CI. The derivative now comes from + * grad_reg_inc_gamma, so forward mode and reverse mode share one algorithm. + */ +TEST(mathMixScalFun, gammaQ_fwd_shape_derivative_large_z) { + using stan::math::fvar; + struct Case { + double a; + double z; + double dqda; + }; + const Case cases[] = { + {1.5, 300.0, 5.71509832058708679e-129}, + {3.0, 650.0, 6.01811954875449219e-277}, + }; + for (const auto& c : cases) { + fvar a = c.a; + a.d_ = 1.0; + fvar z = c.z; + const fvar r = stan::math::gamma_q(a, z); + ASSERT_TRUE(std::isfinite(r.d_)) << "a=" << c.a << " z=" << c.z; + EXPECT_NEAR(c.dqda, r.d_, 1e-10 * std::fabs(c.dqda)) + << "a=" << c.a << " z=" << c.z; + } +} + +/** + * Forward-mode derivative with respect to the second argument, above the + * tgamma overflow limit. The previous form, -exp(-z) * pow(z, a-1) / + * tgamma(a), had both pow and tgamma overflow at (200, 150) and returned + * NaN. It is now evaluated in log space. Reference from mpmath, closed + * form -z^(a-1) e^(-z) / Gamma(a). + */ +TEST(mathMixScalFun, gammaQ_fwd_z_derivative_large_a) { + using stan::math::fvar; + fvar a = 200.0; + fvar z = 150.0; + z.d_ = 1.0; + const fvar r = stan::math::gamma_q(a, z); + const double expected = -2.00507038377120400e-05; + ASSERT_TRUE(std::isfinite(r.d_)); + EXPECT_NEAR(expected, r.d_, 1e-12 * std::fabs(expected)); +} + TEST(MathFunctions, gammaQ_vec) { auto f = [](const auto& x1, const auto& x2) { using stan::math::gamma_q; From ad3ba6982ea37e4f6260adfe1cb6dc5f2c3c2f2f Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sun, 20 Sep 2026 16:45:28 +0300 Subject: [PATCH 5/9] Keep code comments descriptive; move the rationale to the PR text --- stan/math/fwd/fun/Eigen_SpecialFunctions.hpp | 19 +++---- stan/math/fwd/fun/gamma_q.hpp | 24 ++------- stan/math/prim/fun/grad_reg_inc_gamma.hpp | 27 ++++------ .../prim/fun/grad_reg_lower_inc_gamma.hpp | 26 ++++------ stan/math/rev/fun/Eigen_SpecialFunctions.hpp | 19 +++---- test/unit/math/mix/fun/gamma_q_test.cpp | 21 ++------ .../math/mix/fun/grad_reg_inc_gamma_test.cpp | 22 +++----- .../fun/grad_reg_inc_gamma_accuracy_test.cpp | 52 ++++++------------- 8 files changed, 62 insertions(+), 148 deletions(-) diff --git a/stan/math/fwd/fun/Eigen_SpecialFunctions.hpp b/stan/math/fwd/fun/Eigen_SpecialFunctions.hpp index 73a41c05159..3018006b44b 100644 --- a/stan/math/fwd/fun/Eigen_SpecialFunctions.hpp +++ b/stan/math/fwd/fun/Eigen_SpecialFunctions.hpp @@ -11,20 +11,13 @@ namespace Eigen { namespace internal { /** - * Support for Eigen's incomplete gamma routines with Stan's forward-mode - * scalars. + * Specializations of the three `Eigen::internal` helpers that + * `Eigen::internal::igamma_generic_impl` restricts to `float` and + * `double`: `lgamma_impl`, `digamma_impl` and `cephes_helper`. With these, + * `Eigen::numext::igamma_der_a` accepts `fvar` at any autodiff order. * - * `Eigen::internal::igamma_generic_impl` is written for a general scalar - * type and propagates the derivative through the Cephes power series and - * the Cephes continued fraction. Only three of its helpers are restricted - * to `float` and `double`. Supplying those three here lets - * `Eigen::numext::igamma_der_a` serve `fvar` at every autodiff order, - * so `grad_reg_inc_gamma` and `grad_reg_lower_inc_gamma` can use one - * algorithm for every scalar type. - * - * This depends on the names `lgamma_impl`, `digamma_impl` and - * `cephes_helper` in `Eigen::internal`. An Eigen upgrade that renames them - * fails at compile time, not silently. + * These names are Eigen internals. An Eigen upgrade that renames them + * fails here at compile time. */ template struct lgamma_impl> { diff --git a/stan/math/fwd/fun/gamma_q.hpp b/stan/math/fwd/fun/gamma_q.hpp index bf276ebac22..92b57b50989 100644 --- a/stan/math/fwd/fun/gamma_q.hpp +++ b/stan/math/fwd/fun/gamma_q.hpp @@ -16,26 +16,10 @@ namespace stan { namespace math { /* - * Two changes from the previous version of this file. - * - * The derivative with respect to the shape parameter used to be an inlined - * copy of the series in `grad_reg_inc_gamma`, with the same hard-coded 1e-6 - * tolerance and without that function's second branch. The copy is now - * replaced by a call to the root itself. - * - * Keeping the copy hid the defect it shared with the root. In - * `test/prob/chi_square`, the generated `ffv` case compares the - * distribution's analytic partials against autodiff through - * `log(gamma_q(nu * 0.5, y * 0.5))`. Both routes used the same inaccurate - * series, so their errors cancelled and the comparison passed. Once the - * root became accurate and the copy did not, the same comparison failed by - * 2.4e-03 at third order. - * - * The derivative with respect to the second argument used to be - * `-exp(-x2) * pow(x2, x1 - 1) / tgamma(x1)`. Both `pow` and `tgamma` - * overflow for moderate arguments, and the quotient then silently returns - * -0 rather than a small negative number. It is now evaluated in log space, - * which is the form `gamma_p` already uses for the same quantity. + * The derivative with respect to the first argument is grad_reg_inc_gamma. + * The derivative with respect to the second argument is minus the gamma + * density, evaluated in log space so that neither pow nor tgamma can + * overflow; this is the same form that gamma_p uses. */ template diff --git a/stan/math/prim/fun/grad_reg_inc_gamma.hpp b/stan/math/prim/fun/grad_reg_inc_gamma.hpp index 4aa6441f18b..0dbd6e47ebe 100644 --- a/stan/math/prim/fun/grad_reg_inc_gamma.hpp +++ b/stan/math/prim/fun/grad_reg_inc_gamma.hpp @@ -15,30 +15,23 @@ namespace math { * respect to the shape parameter, d/da Q(a, z), where Q is the upper * regularized incomplete gamma function. * - * The work is done by `Eigen::numext::igamma_der_a`, which differentiates - * the Cephes power series and the Cephes continued fraction term by term - * and stops at machine epsilon. It returns d/da of the lower P(a, z), so - * the sign is flipped here. - * - * Eigen's implementation is generic in the scalar type. Stan supplies the - * three Cephes helpers it needs for autodiff scalars in + * Delegates to `Eigen::numext::igamma_der_a`, which returns d/da of the + * lower P(a, z); the sign is flipped here. For autodiff scalar types the + * Cephes helpers that Eigen needs are specialized in * `stan/math/fwd/fun/Eigen_SpecialFunctions.hpp` and - * `stan/math/rev/core/Eigen_SpecialFunctions.hpp`, so every autodiff order - * uses this same algorithm. + * `stan/math/rev/fun/Eigen_SpecialFunctions.hpp`. * - * `g`, `dig`, `precision` and `max_steps` are accepted and ignored. They - * belonged to the previous hand-written series. The signature is unchanged - * so that the 27 call sites do not change; removing the arguments is - * proposed separately. + * `g`, `dig`, `precision` and `max_steps` are accepted for signature + * compatibility and are not used. * * @tparam T1 type of the shape parameter * @tparam T2 type of the location parameter * @param a shape parameter, a > 0 * @param z location z >= 0 - * @param g ignored; previously stan::math::tgamma(a) - * @param dig ignored; previously stan::math::digamma(a) - * @param precision ignored; previously the series tolerance - * @param max_steps ignored; previously the series iteration limit + * @param g unused + * @param dig unused + * @param precision unused + * @param max_steps unused * @return d/da of the upper regularized incomplete gamma function */ template diff --git a/stan/math/prim/fun/grad_reg_lower_inc_gamma.hpp b/stan/math/prim/fun/grad_reg_lower_inc_gamma.hpp index 7f154f5eed1..ebde7f2aa1b 100644 --- a/stan/math/prim/fun/grad_reg_lower_inc_gamma.hpp +++ b/stan/math/prim/fun/grad_reg_lower_inc_gamma.hpp @@ -15,31 +15,23 @@ namespace math { * Computes the gradient of the lower regularized incomplete gamma * function, d/da P(a, z). * - * The work is done by `Eigen::numext::igamma_der_a`, which returns exactly - * this quantity. It differentiates the Cephes power series and the Cephes - * continued fraction term by term and stops at machine epsilon. Eigen's - * implementation is generic in the scalar type; Stan supplies the three - * Cephes helpers it needs for autodiff scalars in - * `stan/math/fwd/fun/Eigen_SpecialFunctions.hpp` and - * `stan/math/rev/core/Eigen_SpecialFunctions.hpp`. + * Delegates to `Eigen::numext::igamma_der_a`, which returns exactly this + * quantity. For autodiff scalar types the Cephes helpers that Eigen needs + * are specialized in `stan/math/fwd/fun/Eigen_SpecialFunctions.hpp` and + * `stan/math/rev/fun/Eigen_SpecialFunctions.hpp`. * - * The previous implementation used a Gautschi series whose two sums cancel - * completely for z well above a. It returned about 1e-14 of rounding noise, - * with arbitrary sign, where the true value can be 1e-148. + * `precision` and `max_steps` are accepted for signature compatibility and + * are not used. * - * `precision` and `max_steps` are accepted and ignored. They belonged to - * the previous series. The signature is unchanged so that the call sites do - * not change; removing the arguments is proposed separately. - * - * Infinite z still throws, to match the behaviour of, for example, + * Infinite z throws, to match the behaviour of, for example, * boost::math::gamma_p. * * @tparam T1 type of a * @tparam T2 type of z * @param[in] a shared with complete Gamma, a > 0 * @param[in] z value to integrate up to, z >= 0 - * @param[in] precision ignored; previously the series tolerance - * @param[in] max_steps ignored; previously the series iteration limit + * @param[in] precision unused + * @param[in] max_steps unused * @return d/da of the lower regularized incomplete gamma function */ template diff --git a/stan/math/rev/fun/Eigen_SpecialFunctions.hpp b/stan/math/rev/fun/Eigen_SpecialFunctions.hpp index 458104c7933..054a10057e9 100644 --- a/stan/math/rev/fun/Eigen_SpecialFunctions.hpp +++ b/stan/math/rev/fun/Eigen_SpecialFunctions.hpp @@ -11,19 +11,14 @@ namespace Eigen { namespace internal { /** - * Support for Eigen's incomplete gamma routines with Stan's reverse-mode - * scalar. See `stan/math/fwd/fun/Eigen_SpecialFunctions.hpp` for the - * reasoning; this is the same three helpers for `var`. + * The same three `Eigen::internal` helper specializations as + * `stan/math/fwd/fun/Eigen_SpecialFunctions.hpp`, for `var`, so that + * `Eigen::numext::igamma_der_a` accepts `var` arguments. * - * No function inside Stan Math reaches `igamma_der_a` with `var`. Every - * distribution calls the gradient roots with `T_partials_return`, which is - * `double` for reverse mode. These specializations exist so that a direct - * call such as `grad_reg_inc_gamma(var, var, var, var)` keeps working after - * the hand-written series was removed. - * - * A `var` call builds one autodiff node per operation inside the series, so - * it is far more expensive than the `double` call. That was equally true of - * the previous implementation. + * Inside Stan Math the gradient roots are called with `double` partials in + * reverse mode; this path serves direct calls with `var`. Such a call + * records one autodiff node per operation of the series, so it is far more + * expensive than the `double` call. */ template <> struct lgamma_impl { diff --git a/test/unit/math/mix/fun/gamma_q_test.cpp b/test/unit/math/mix/fun/gamma_q_test.cpp index 0ef0bef08c2..ee474c8f20d 100644 --- a/test/unit/math/mix/fun/gamma_q_test.cpp +++ b/test/unit/math/mix/fun/gamma_q_test.cpp @@ -20,19 +20,8 @@ TEST(mathMixScalFun, gammaQ) { } /** - * Forward-mode shape derivative at large z, fixed references from mpmath - * at 80 digits. - * - * The previous forward-mode gamma_q carried its own copy of the series in - * grad_reg_inc_gamma, without an iteration cap. For z above about 700 the - * alternating terms overflow before they decay and the loop never - * terminates; that is the failure the commented-out expect_value above - * ran into. Below that, the series terminates and returns garbage: at - * (a, z) = (1.5, 300) it gave 1.2e+112 where the derivative is 5.7e-129. - * - * These points are chosen so that the old code fails rather than hangs, - * which keeps the test usable in CI. The derivative now comes from - * grad_reg_inc_gamma, so forward mode and reverse mode share one algorithm. + * Forward-mode derivative with respect to the first argument at large z. + * References are d/da Q(a, z) from mpmath at 80 digits. */ TEST(mathMixScalFun, gammaQ_fwd_shape_derivative_large_z) { using stan::math::fvar; @@ -57,10 +46,8 @@ TEST(mathMixScalFun, gammaQ_fwd_shape_derivative_large_z) { } /** - * Forward-mode derivative with respect to the second argument, above the - * tgamma overflow limit. The previous form, -exp(-z) * pow(z, a-1) / - * tgamma(a), had both pow and tgamma overflow at (200, 150) and returned - * NaN. It is now evaluated in log space. Reference from mpmath, closed + * Forward-mode derivative with respect to the second argument, at an a + * above the range where tgamma(a) is finite. Reference from mpmath, closed * form -z^(a-1) e^(-z) / Gamma(a). */ TEST(mathMixScalFun, gammaQ_fwd_z_derivative_large_a) { diff --git a/test/unit/math/mix/fun/grad_reg_inc_gamma_test.cpp b/test/unit/math/mix/fun/grad_reg_inc_gamma_test.cpp index af11db17f01..8ab0ddda58a 100644 --- a/test/unit/math/mix/fun/grad_reg_inc_gamma_test.cpp +++ b/test/unit/math/mix/fun/grad_reg_inc_gamma_test.cpp @@ -8,12 +8,8 @@ TEST(ProbInternalMath, gradRegIncGamma_typical) { double g = 1.77245; double dig = -1.96351; - // The expected value was 0.38984156, which is 1.1e-5 away from the true - // d/da Q(0.5, 1.0). That constant encoded two errors of the old - // implementation: the rounded support values that this test passes in - // (g = 1.77245 is 2.2e-6 away from tgamma(0.5)), and the 1e-6 series - // truncation. The arithmetic path no longer uses g or dig, so the result - // is now the true value, from mpmath at 80 digits. + // d/da Q(0.5, 1.0) from mpmath at 80 digits. The g and dig arguments are + // accepted for signature compatibility and do not affect the result. EXPECT_FLOAT_EQ(0.38983726432851057, stan::math::grad_reg_inc_gamma(a, b, g, dig)); } @@ -37,8 +33,6 @@ TEST(ProbInternalMath, gradRegIncGamma_largeZ) { stan::math::grad_reg_inc_gamma(a, z, g, dig)); } -// Every scalar type now uses the same algorithm, so the fvar cases expect -// the same value as the arithmetic case. TEST(ProbInternalMath, gradRegIncGamma_fd) { using stan::math::fvar; @@ -130,14 +124,10 @@ struct gamma_cdf_shape_functor { } // namespace /** - * The first derivative must not depend on how it is obtained. - * - * Reverse mode reaches the gradient root with `double` partials; `hessian` - * reaches it with `fvar`. While those were two different algorithms, - * the two routes disagreed by 4.4e-02 at alpha = 20, 5.8e-01 at alpha = 100 - * and 7.2e-01 at alpha = 171, and `hessian` returned NaN at alpha = 180, - * all with `z = alpha`. This test locks that down. It asserts agreement - * between the two routes, so it needs no external reference value. + * The first derivative must not depend on how it is obtained. Reverse mode + * reaches the gradient root with `double` partials and `hessian` reaches + * it with `fvar`; the two must agree. This needs no external + * reference value. */ TEST(ProbInternalMath, gradRegIncGamma_gradient_matches_hessian) { using stan::math::var; diff --git a/test/unit/math/prim/fun/grad_reg_inc_gamma_accuracy_test.cpp b/test/unit/math/prim/fun/grad_reg_inc_gamma_accuracy_test.cpp index 346e05c6b59..22a60fd2486 100644 --- a/test/unit/math/prim/fun/grad_reg_inc_gamma_accuracy_test.cpp +++ b/test/unit/math/prim/fun/grad_reg_inc_gamma_accuracy_test.cpp @@ -4,29 +4,22 @@ #include /** - * Fixed-reference accuracy tests for the incomplete gamma gradient roots. + * Fixed-reference tests for grad_reg_inc_gamma and grad_reg_lower_inc_gamma. * - * Every expected value comes from mpmath at 80 digits, through two - * independent routes that agree to at least 13 digits: numerical - * differentiation of the regularized function, and an explicit integral - * identity evaluated in density form. A third route, the Gautschi series, - * confirms the deep-tail values. The generator is - * `inc_gamma_beta_work/ref_mpmath.py` in the audit workspace. + * Expected values are d/da of the regularized incomplete gamma function + * from mpmath at 80 digits, checked through two independent routes that + * agree to at least 13 digits. They reproduce the two constants already in + * grad_reg_inc_gamma_test.cpp: (a, z) = (1.1, 0.2) gives + * 0.31416364892410886 and (9, 10) gives 0.120855166827777. * - * The references reproduce two constants that were already in - * grad_reg_inc_gamma_test.cpp, which is the check that the two agree: - * (a, z) = (1.1, 0.2) gives 0.31416364892410886, and (9, 10) gives - * 0.120855166827777. - * - * These tests fail on develop at 5252d51d47. They exist because - * `expect_ad` compares against finite differences and cannot see an error - * of this size, and because the existing suites assert absolute tolerances - * on a grid that stops short of every failing region. + * These are absolute references rather than autodiff comparisons because + * expect_ad compares against finite differences, whose own error is far + * larger than the accuracy asserted here. */ namespace { -/** Relative comparison, needed because the values span 1e-1 to 1e-148. */ +/** Relative comparison, because the values span 1e-1 to 1e-148. */ void expect_rel(double expected, double actual, double tol, const std::string& label) { ASSERT_TRUE(std::isfinite(actual)) @@ -43,33 +36,21 @@ double upper(double a, double z) { } // namespace -/** - * D1. `z >= a && z >= 8` selects a fixed 10-term asymptotic expansion that - * needs z >> a. On develop the error is 4.4e-2 at a = 20, 5.8e-1 at - * a = 100 and 7.2e-1 at a = 171. - */ +/** z equal to a, moderate and large a. */ TEST(MathPrimScalFun, grad_reg_inc_gamma_diagonal_z_equals_a) { expect_rel(8.95792055391225334e-02, upper(20.0, 20.0), 1e-10, "a=z=20"); expect_rel(3.99274978577860498e-02, upper(100.0, 100.0), 1e-10, "a=z=100"); expect_rel(3.05227525380366363e-02, upper(171.0, 171.0), 1e-10, "a=z=171"); } -/** - * D2. Above a = 171.62, tgamma(a) is inf, so the correction term is - * inf/inf. On develop every one of these returns NaN. - */ +/** a above the range where tgamma(a) is finite. */ TEST(MathPrimScalFun, grad_reg_inc_gamma_above_tgamma_overflow) { expect_rel(2.97491740225716214e-02, upper(180.0, 180.0), 1e-10, "a=z=180"); expect_rel(2.30393429308686434e-02, upper(300.0, 300.0), 1e-10, "a=z=300"); expect_rel(1.78442151466593514e-02, upper(500.0, 500.0), 1e-10, "a=z=500"); } -/** - * D3. `grad_reg_lower_inc_gamma` runs its Gautschi branch for every z > 42 - * with a >= 12, and there `emz * (log_z * sum_a - sum_b)` cancels - * completely. On develop it returns about 1e-14 of noise, with arbitrary - * sign, at all four points. - */ +/** The far upper tail, z well above a, where the values are very small. */ TEST(MathPrimScalFun, grad_reg_lower_inc_gamma_far_upper_tail) { expect_rel(-8.28941360740501778e-22, stan::math::grad_reg_lower_inc_gamma(50.0, 150.0), 1e-10, @@ -86,9 +67,8 @@ TEST(MathPrimScalFun, grad_reg_lower_inc_gamma_far_upper_tail) { } /** - * The two functions compute the same quantity with opposite sign. On - * develop they disagree by 72 % at a = z = 171. This test needs no external - * reference at all. + * The two functions compute the same quantity with opposite sign, so they + * must agree with each other. This needs no external reference. */ TEST(MathPrimScalFun, grad_reg_inc_gamma_agrees_with_lower) { for (double a : {5.0, 20.0, 50.0, 100.0, 171.0}) { @@ -98,7 +78,7 @@ TEST(MathPrimScalFun, grad_reg_inc_gamma_agrees_with_lower) { } } -/** Points that already worked, to show the change is not a regression. */ +/** Points in the series region with small a. */ TEST(MathPrimScalFun, grad_reg_inc_gamma_unchanged_points) { expect_rel(3.14163648924108863e-01, upper(1.1, 0.2), 1e-10, "a=1.1 z=0.2"); expect_rel(1.20855166827777014e-01, upper(9.0, 10.0), 1e-10, "a=9 z=10"); From 9341c4ce3883c31a3e03a63ed4a9ba7f89ad6cf4 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Wed, 21 Oct 2026 17:38:23 +0300 Subject: [PATCH 6/9] Compute the gamma_lccdf shape derivative in log space where Q underflows --- stan/math/prim/fun/log_gamma_q_dgamma.hpp | 95 +++++++++++++++++++ stan/math/prim/prob/gamma_lccdf.hpp | 41 +++----- test/unit/math/mix/prob/gamma_lccdf_test.cpp | 62 ++++++++++++ .../rev/prob/gamma_lccdf_tail_reference.hpp | 90 ++++++++++++++++++ test/unit/math/rev/prob/gamma_lccdf_test.cpp | 23 +++++ 5 files changed, 283 insertions(+), 28 deletions(-) create mode 100644 test/unit/math/mix/prob/gamma_lccdf_test.cpp create mode 100644 test/unit/math/rev/prob/gamma_lccdf_tail_reference.hpp diff --git a/stan/math/prim/fun/log_gamma_q_dgamma.hpp b/stan/math/prim/fun/log_gamma_q_dgamma.hpp index 742b23784ba..7433bcbc784 100644 --- a/stan/math/prim/fun/log_gamma_q_dgamma.hpp +++ b/stan/math/prim/fun/log_gamma_q_dgamma.hpp @@ -16,7 +16,10 @@ #include #include #include +#include #include +#include +#include namespace stan { namespace math { @@ -65,12 +68,97 @@ inline return_type_t log_q_gamma_cf(const T_a& a, const T_z& z, f *= delta; const double delta_m1 = fabs(value_of_rec(delta) - 1.0); if (delta_m1 < precision) { + if constexpr (is_fvar::value) { + // The tangent converges more slowly than the value near z = a + 1. + const double tangent_step = fabs(value_of_rec(delta.d_)); + const double tangent_scale = fabs(value_of_rec(f.d_ / f.val_)); + if (tangent_step > precision * std::max(1.0, tangent_scale)) { + continue; + } + } break; } } return log_prefactor - log(f); } +/** + * Compute log(Q(a,z)) and d/da log(Q(a,z)) for z > a + 1 from the + * continued fraction of log_q_gamma_cf, with the derivative with respect + * to a carried through the recurrence. Both results stay in log space, so + * they are finite where Q(a,z) underflows. + * + * @tparam T_a Type of shape parameter a (double or autodiff types) + * @tparam T_z Type of value parameter z (double or autodiff types) + * @param a Shape parameter + * @param z Value at which to evaluate + * @param precision Convergence threshold for the value and the derivative + * @param max_steps Maximum number of continued fraction iterations + * @return pair of log(Q(a,z)) and d/da log(Q(a,z)) + */ +template +inline std::pair, return_type_t> +log_q_gamma_cf_and_da(const T_a& a, const T_z& z, + double precision = LOG_Q_GAMMA_CF_PRECISION, + int max_steps = 250) { + using T_return = return_type_t; + const T_return log_z = log(z); + const T_return log_prefactor = a * log_z - z - lgamma(a); + const T_return d_log_prefactor = log_z - digamma(a); + + // Each quantity x is carried with its derivative dx/da. + const T_return b_init = z + 1.0 - a; + T_return C = b_init; + T_return dC = -1.0; + if (fabs(value_of_rec(C)) < EPSILON) { + C = EPSILON; + dC = 0.0; + } + T_return D = 0.0; + T_return dD = 0.0; + T_return f = C; + T_return df = dC; + for (int i = 1; i <= max_steps; ++i) { + const T_a an = -i * (i - a); + const T_return b = b_init + 2.0 * i; + dD = -1.0 + i * D + an * dD; + D = b + an * D; + if (fabs(value_of_rec(D)) < EPSILON) { + D = EPSILON; + dD = 0.0; + } + dC = -1.0 + (i * C - an * dC) / (C * C); + C = b + an / C; + if (fabs(value_of_rec(C)) < EPSILON) { + C = EPSILON; + dC = 0.0; + } + dD = -dD / (D * D); + D = inv(D); + const T_return delta = C * D; + const T_return d_delta = dC * D + C * dD; + df = df * delta + f * d_delta; + f *= delta; + const double tangent_scale = fabs(value_of_rec(df / f)); + if (fabs(value_of_rec(delta) - 1.0) < precision + && fabs(value_of_rec(d_delta)) + <= precision * std::max(1.0, tangent_scale)) { + break; + } + } + return {log_prefactor - log(f), d_log_prefactor - df / f}; +} + +/** + * Lower limit of log(Q(a,z)) for d/da log(Q) = (dQ/da) / Q. Below it Q or + * dQ/da can be subnormal or zero, and log_q_gamma_cf_and_da must be used. + */ +inline double log_q_gamma_ratio_min() { + static const double log_min + = std::log(std::numeric_limits::min()) - LOG_EPSILON; + return log_min; +} + } // namespace internal /** @@ -102,6 +190,13 @@ log_gamma_q_dgamma(const T_a& a, const T_z& z, if (z_val > a_val + 1.0) { std::pair result{ internal::log_q_gamma_cf(a_val, z_val, precision, max_steps), 0.0}; + if (value_of_rec(result.first) < internal::log_q_gamma_ratio_min()) { + const std::pair log_q_and_da + = internal::log_q_gamma_cf_and_da(a_val, z_val, precision, max_steps); + result.first = log_q_and_da.first; + result.second = log_q_and_da.second; + return result; + } // For gradient, use: d/da log(Q) = (1/Q) * dQ/da // grad_reg_inc_gamma computes dQ/da const T_return Q_val = exp(result.first); diff --git a/stan/math/prim/prob/gamma_lccdf.hpp b/stan/math/prim/prob/gamma_lccdf.hpp index c5538fa5966..c91ac01971f 100644 --- a/stan/math/prim/prob/gamma_lccdf.hpp +++ b/stan/math/prim/prob/gamma_lccdf.hpp @@ -31,8 +31,7 @@ namespace internal { /** * Computes log q and d(log q) / d(alpha) using continued fraction. */ -template +template inline std::optional, return_type_t>> eval_q_cf(const T1& alpha, const T2& beta_y) { using scalar_t = return_type_t; @@ -51,17 +50,16 @@ eval_q_cf(const T1& alpha, const T2& beta_y) { return std::optional{std::nullopt}; } if constexpr (is_autodiff_v) { - if constexpr (!partials_fvar) { + if constexpr (is_autodiff_v) { + // The partial is differentiated again, so it is computed with the + // derivative recurrence in the autodiff type. + out.second = internal::log_q_gamma_cf_and_da(alpha, beta_y).second; + } else if (value_of_rec(out.first) < internal::log_q_gamma_ratio_min()) { + out.second = internal::log_q_gamma_cf_and_da(alpha, beta_y).second; + } else { out.second = grad_reg_inc_gamma(alpha, beta_y, tgamma(alpha), digamma(alpha)) / exp(out.first); - } else { - auto alpha_unit = alpha; - alpha_unit.d_ = 1; - auto beta_y_unit = beta_y; - beta_y_unit.d_ = 0; - auto log_Q_fvar = internal::log_q_gamma_cf(alpha_unit, beta_y_unit); - out.second = log_Q_fvar.d_; } } return std::optional{out}; @@ -71,7 +69,7 @@ eval_q_cf(const T1& alpha, const T2& beta_y) { /** * Computes log q and d(log q) / d(alpha) using log1m. */ -template +template inline std::optional, return_type_t>> eval_q_log1m(const T1& alpha, const T2& beta_y) { using scalar_t = return_type_t; @@ -81,16 +79,7 @@ eval_q_log1m(const T1& alpha, const T2& beta_y) { return std::optional{std::nullopt}; } if constexpr (is_autodiff_v) { - if constexpr (partials_fvar) { - auto alpha_unit = alpha; - alpha_unit.d_ = 1; - auto beta_unit = beta_y; - beta_unit.d_ = 0; - auto log_Q_fvar = log1m(gamma_p(alpha_unit, beta_unit)); - out.second = log_Q_fvar.d_; - } else { - out.second = -grad_reg_lower_inc_gamma(alpha, beta_y) / exp(out.first); - } + out.second = -grad_reg_lower_inc_gamma(alpha, beta_y) / exp(out.first); } return std::optional{out}; } @@ -131,7 +120,6 @@ inline return_type_t gamma_lccdf( constexpr bool is_shape_fvar = is_fvar_v>; constexpr bool is_beta_fvar = is_fvar_v>; constexpr bool any_fvar = is_y_fvar || is_shape_fvar || is_beta_fvar; - constexpr bool partials_fvar = is_fvar_v; for (size_t n = 0; n < N; n++) { // Explicit results for extreme values @@ -153,15 +141,12 @@ inline return_type_t gamma_lccdf( } std::optional> result; if (beta_y > alpha_val + 1.0) { - result = internal::eval_q_cf(alpha_val, - beta_y); + result = internal::eval_q_cf(alpha_val, beta_y); } else { - result - = internal::eval_q_log1m(alpha_val, beta_y); + result = internal::eval_q_log1m(alpha_val, beta_y); if (!result && beta_y > 0.0) { // Fallback to continued fraction if log1m fails - result = internal::eval_q_cf( - alpha_val, beta_y); + result = internal::eval_q_cf(alpha_val, beta_y); } } if (unlikely(!result)) { diff --git a/test/unit/math/mix/prob/gamma_lccdf_test.cpp b/test/unit/math/mix/prob/gamma_lccdf_test.cpp new file mode 100644 index 00000000000..acde3221dec --- /dev/null +++ b/test/unit/math/mix/prob/gamma_lccdf_test.cpp @@ -0,0 +1,62 @@ +#include +#include +#include +#include +#include + +TEST_F(AgradRev, mathMixScalFun_gamma_lccdf) { + auto f = [](const auto& y, const auto& alpha, const auto& beta) { + return stan::math::gamma_lccdf(y, alpha, beta); + }; + + stan::test::expect_ad(f, 0.8, 1.1, 2.3); + stan::test::expect_ad(f, 2.5, 3.0, 1.5); + stan::test::expect_ad(f, 6.5, 5.0, 1.0); +} + +TEST_F(AgradRev, mathMixScalFun_gamma_lccdf_upper_tail_reference) { + using stan::math::fvar; + using stan::math::gamma_lccdf; + using stan::math::var; + namespace ref = stan::test::gamma_lccdf_tail; + + for (const auto& r : ref::rows()) { + const std::string at = ref::label(r); + + fvar a1(r.alpha, 1.0); + ref::expect_rel(r.d_alpha, gamma_lccdf(r.y, a1, r.beta).d_, 1e-11, + "fvar d/dalpha " + at); + + var a2_val = r.alpha; + fvar a2(a2_val, 1.0); + fvar lp2 = gamma_lccdf(r.y, a2, r.beta); + ref::expect_rel(r.d_alpha, lp2.d_.val(), 1e-11, "fvar d/dalpha " + at); + lp2.d_.grad(); + ref::expect_rel(r.d2_alpha, a2_val.adj(), ref::d2_tol(r), + "fvar d2/dalpha2 " + at); + stan::math::recover_memory(); + + fvar> a3(fvar(r.alpha, 1.0), fvar(1.0, 0.0)); + fvar> lp3 = gamma_lccdf(r.y, a3, r.beta); + ref::expect_rel(r.value, lp3.val_.val_, 1e-10, + "fvar> value " + at); + ref::expect_rel(r.d_alpha, lp3.val_.d_, 1e-11, + "fvar> inner d/dalpha " + at); + ref::expect_rel(r.d_alpha, lp3.d_.val_, 1e-11, + "fvar> outer d/dalpha " + at); + ref::expect_rel(r.d2_alpha, lp3.d_.d_, ref::d2_tol(r), + "fvar> d2/dalpha2 " + at); + + var a4_val = r.alpha; + fvar> a4(fvar(a4_val, 1.0), fvar(1.0, 0.0)); + fvar> lp4 = gamma_lccdf(r.y, a4, r.beta); + ref::expect_rel(r.d_alpha, lp4.d_.val_.val(), 1e-11, + "fvar> d/dalpha " + at); + ref::expect_rel(r.d2_alpha, lp4.d_.d_.val(), ref::d2_tol(r), + "fvar> d2/dalpha2 " + at); + lp4.d_.d_.grad(); + EXPECT_TRUE(std::isfinite(a4_val.adj())) + << "fvar> d3/dalpha3 " << at; + stan::math::recover_memory(); + } +} diff --git a/test/unit/math/rev/prob/gamma_lccdf_tail_reference.hpp b/test/unit/math/rev/prob/gamma_lccdf_tail_reference.hpp new file mode 100644 index 00000000000..3ea31eb2549 --- /dev/null +++ b/test/unit/math/rev/prob/gamma_lccdf_tail_reference.hpp @@ -0,0 +1,90 @@ +#ifndef TEST_UNIT_MATH_REV_PROB_GAMMA_LCCDF_TAIL_REFERENCE_HPP +#define TEST_UNIT_MATH_REV_PROB_GAMMA_LCCDF_TAIL_REFERENCE_HPP + +#include +#include +#include +#include + +namespace stan { +namespace test { +namespace gamma_lccdf_tail { + +/** + * gamma_lccdf(y | alpha, beta), its first and second derivatives with + * respect to alpha, and its derivative with respect to y. + * + * Computed with mpmath at 50 digits by two routes: the regularized upper + * incomplete gamma function with numerical differentiation, and a + * quadrature of the density scaled to 1 at the lower limit, with the + * derivatives as ratios of positive integrals. The routes agree to 1e-47 + * where both converge; at alpha = 1e4 and 1e6 only the quadrature + * converges. + * + * Rows 0 and 1 are near y = alpha + 1. Rows 2 to 9 have + * log(Q) < log(DBL_MIN), where Q itself is subnormal or zero. + */ +struct Row { + double alpha; + double y; + double beta; + double value; + double d_alpha; + double d2_alpha; + double d_y; +}; + +inline const std::vector& rows() { + static const std::vector r = { + {5, 6.5, 1, -1.4975754049608005, 5.8767262159542573e-1, + -1.9055199558693436e-1, -4.9993873514327224e-1}, + {5, 7.5, 1, -2.0244848579180109, 6.9436601738250379e-1, + -1.9817346506198028e-1, -5.5213820632791284e-1}, + {150, 1200, 1, -7.4345557714901795e+2, 2.0837288947176418, + -6.6880371230448705e-3, -8.7595121899782028e-1}, + {300, 1500, 1, -7.2232706800763835e+2, 1.6119371090820549, + -3.3382049273506081e-3, -8.0083229548002851e-1}, + {150, 600, 2, -7.4345557714901795e+2, 2.0837288947176418, + -6.6880371230448705e-3, -1.7519024379956406}, + {10, 1000, 1, -9.5062299837015608e+2, 4.6570107370619179, + -1.0516532158355843e-1, -9.9100906329734302e-1}, + {20, 1000, 1, -9.0807337076913739e+2, 3.9382495777101035, + -5.1269788164528015e-2, -9.8101932786738336e-1}, + {200, 1600, 1, -9.8962693094236584e+2, 2.0826567480373627, + -5.0120131147584039e-3, -8.757136315228946e-1}, + {1000, 3000, 1, -9.0645450699007447e+2, 1.0996116242418308, + -1.0002512857967094e-3, -6.6716616809812125e-1}, + {1e6, 1.1e6, 1, -4.6953444144485314e+3, 9.5320677605584263e-2, + -9.9990056594105228e-7, -9.0919088810213022e-2}, + {700, 2100, 1, -6.3586021219348883e+2, 1.1000395060722886, + -1.4290853740612502e-3, -6.6737993613863034e-1}, + {1e4, 10100, 1, -1.84104688482372, 1.5226403888940706e-2, + -8.0334369099635953e-5, -1.5150648133166995e-2}, + {200, 200, 1, -7.1213311417597682e-1, 5.7524325573397242e-2, + -3.3091679741037347e-3, -5.7476404689539276e-2}, + }; + return r; +} + +/** Tolerance for second derivatives: the continued fraction stops on the + * convergence of the value and the first derivative only. */ +inline double d2_tol(const Row& r) { return r.y < r.alpha + 3 ? 5e-8 : 1e-9; } + +inline std::string label(const Row& r) { + return "alpha=" + std::to_string(r.alpha) + " y=" + std::to_string(r.y) + + " beta=" + std::to_string(r.beta); +} + +inline void expect_rel(double expected, double actual, double tol, + const std::string& what) { + ASSERT_TRUE(std::isfinite(actual)) << what << ": not finite, " << actual; + const double err = std::fabs(actual - expected) / std::fabs(expected); + EXPECT_LT(err, tol) << what << ": relative error " << err << ", expected " + << expected << ", got " << actual; +} + +} // namespace gamma_lccdf_tail +} // namespace test +} // namespace stan + +#endif diff --git a/test/unit/math/rev/prob/gamma_lccdf_test.cpp b/test/unit/math/rev/prob/gamma_lccdf_test.cpp index d069a365578..778d05d52f7 100644 --- a/test/unit/math/rev/prob/gamma_lccdf_test.cpp +++ b/test/unit/math/rev/prob/gamma_lccdf_test.cpp @@ -1,5 +1,6 @@ #include #include +#include #include #include #include @@ -480,3 +481,25 @@ TEST_F(AgradRev, } } } + +TEST_F(AgradRev, ProbDistributionsGamma_lccdf_upper_tail_reference) { + using stan::math::gamma_lccdf; + using stan::math::var; + namespace ref = stan::test::gamma_lccdf_tail; + + for (const auto& r : ref::rows()) { + const std::string at = ref::label(r); + var y = r.y; + var alpha = r.alpha; + var beta = r.beta; + var lp = gamma_lccdf(y, alpha, beta); + lp.grad(); + // The value is limited by the cancellation in + // alpha log(z) - z - lgamma(alpha) at large alpha. + ref::expect_rel(r.value, lp.val(), 1e-10, "value " + at); + ref::expect_rel(r.d_alpha, alpha.adj(), 1e-11, "d/dalpha " + at); + ref::expect_rel(r.d_y, y.adj(), 1e-8, "d/dy " + at); + ref::expect_rel(r.d_y * r.y / r.beta, beta.adj(), 1e-8, "d/dbeta " + at); + stan::math::recover_memory(); + } +} From 8d376d5763ed96bc20356721c3c5cbe520b78b08 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Wed, 21 Oct 2026 18:18:47 +0300 Subject: [PATCH 7/9] Use the log-space continued-fraction derivative for every z > a + 1 in gamma_lccdf --- stan/math/prim/fun/log_gamma_q_dgamma.hpp | 33 ++++------------------- stan/math/prim/prob/gamma_lccdf.hpp | 21 +++++++-------- 2 files changed, 15 insertions(+), 39 deletions(-) diff --git a/stan/math/prim/fun/log_gamma_q_dgamma.hpp b/stan/math/prim/fun/log_gamma_q_dgamma.hpp index 7433bcbc784..abfdfbeda89 100644 --- a/stan/math/prim/fun/log_gamma_q_dgamma.hpp +++ b/stan/math/prim/fun/log_gamma_q_dgamma.hpp @@ -18,7 +18,6 @@ #include #include #include -#include #include namespace stan { @@ -149,16 +148,6 @@ log_q_gamma_cf_and_da(const T_a& a, const T_z& z, return {log_prefactor - log(f), d_log_prefactor - df / f}; } -/** - * Lower limit of log(Q(a,z)) for d/da log(Q) = (dQ/da) / Q. Below it Q or - * dQ/da can be subnormal or zero, and log_q_gamma_cf_and_da must be used. - */ -inline double log_q_gamma_ratio_min() { - static const double log_min - = std::log(std::numeric_limits::min()) - LOG_EPSILON; - return log_min; -} - } // namespace internal /** @@ -186,24 +175,12 @@ log_gamma_q_dgamma(const T_a& a, const T_z& z, using T_return = return_type_t; const double a_val = value_of(a); const double z_val = value_of(z); - // For z > a + 1, use continued fraction for better numerical stability + // For z > a + 1, use the continued fraction, with the derivative carried + // through it in log space; it stays finite where Q underflows if (z_val > a_val + 1.0) { - std::pair result{ - internal::log_q_gamma_cf(a_val, z_val, precision, max_steps), 0.0}; - if (value_of_rec(result.first) < internal::log_q_gamma_ratio_min()) { - const std::pair log_q_and_da - = internal::log_q_gamma_cf_and_da(a_val, z_val, precision, max_steps); - result.first = log_q_and_da.first; - result.second = log_q_and_da.second; - return result; - } - // For gradient, use: d/da log(Q) = (1/Q) * dQ/da - // grad_reg_inc_gamma computes dQ/da - const T_return Q_val = exp(result.first); - const double dQ_da - = grad_reg_inc_gamma(a_val, z_val, tgamma(a_val), digamma(a_val)); - result.second = dQ_da / Q_val; - return result; + const std::pair log_q_and_da + = internal::log_q_gamma_cf_and_da(a_val, z_val, precision, max_steps); + return {log_q_and_da.first, log_q_and_da.second}; } else { // For z <= a + 1, use log1m(P(a,z)) for better numerical accuracy const double P_val = gamma_p(a_val, z_val); diff --git a/stan/math/prim/prob/gamma_lccdf.hpp b/stan/math/prim/prob/gamma_lccdf.hpp index c91ac01971f..f2db4ea25bb 100644 --- a/stan/math/prim/prob/gamma_lccdf.hpp +++ b/stan/math/prim/prob/gamma_lccdf.hpp @@ -44,23 +44,22 @@ eval_q_cf(const T1& alpha, const T2& beta_y) { } else { return std::optional{std::nullopt}; } + } else if constexpr (is_autodiff_v && !is_autodiff_v) { + // One continued fraction gives the value and the derivative. + const ret_t out = internal::log_q_gamma_cf_and_da(alpha, beta_y); + if (unlikely(!std::isfinite(out.first))) { + return std::optional{std::nullopt}; + } + return std::optional{out}; } else { ret_t out{internal::log_q_gamma_cf(alpha, beta_y), 0.0}; if (unlikely(!std::isfinite(value_of_rec(out.first)))) { return std::optional{std::nullopt}; } if constexpr (is_autodiff_v) { - if constexpr (is_autodiff_v) { - // The partial is differentiated again, so it is computed with the - // derivative recurrence in the autodiff type. - out.second = internal::log_q_gamma_cf_and_da(alpha, beta_y).second; - } else if (value_of_rec(out.first) < internal::log_q_gamma_ratio_min()) { - out.second = internal::log_q_gamma_cf_and_da(alpha, beta_y).second; - } else { - out.second - = grad_reg_inc_gamma(alpha, beta_y, tgamma(alpha), digamma(alpha)) - / exp(out.first); - } + // The partial is differentiated again, so it is computed with the + // derivative recurrence in the autodiff type. + out.second = internal::log_q_gamma_cf_and_da(alpha, beta_y).second; } return std::optional{out}; } From 7d3ea2f5c6cca180505f3d9e7ed2a9621e166a1a Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Wed, 21 Oct 2026 19:14:02 +0300 Subject: [PATCH 8/9] Compute the incomplete-gamma cdfs of six distributions in log space --- stan/math/prim/fun/log_gamma_q_dgamma.hpp | 155 ++++++++ stan/math/prim/prob/chi_square_cdf.hpp | 98 ++--- stan/math/prim/prob/chi_square_lccdf.hpp | 99 ++--- stan/math/prim/prob/chi_square_lcdf.hpp | 97 ++--- stan/math/prim/prob/gamma_cdf.hpp | 115 +++--- stan/math/prim/prob/gamma_lccdf.hpp | 177 +++------ stan/math/prim/prob/gamma_lcdf.hpp | 104 ++--- stan/math/prim/prob/inv_chi_square_cdf.hpp | 95 ++--- stan/math/prim/prob/inv_chi_square_lccdf.hpp | 101 ++--- stan/math/prim/prob/inv_chi_square_lcdf.hpp | 97 ++--- stan/math/prim/prob/inv_gamma_cdf.hpp | 119 +++--- stan/math/prim/prob/inv_gamma_lccdf.hpp | 115 +++--- stan/math/prim/prob/inv_gamma_lcdf.hpp | 111 +++--- stan/math/prim/prob/poisson_cdf.hpp | 61 ++- stan/math/prim/prob/poisson_lccdf.hpp | 55 ++- stan/math/prim/prob/poisson_lcdf.hpp | 58 ++- .../prim/prob/scaled_inv_chi_square_cdf.hpp | 121 +++--- .../prim/prob/scaled_inv_chi_square_lccdf.hpp | 116 +++--- .../prim/prob/scaled_inv_chi_square_lcdf.hpp | 112 +++--- .../mix/prob/incomplete_gamma_cdfs_test.cpp | 58 +++ .../prob/inv_incomplete_gamma_cdfs_test.cpp | 53 +++ test/unit/math/prim/prob/gamma_lccdf_test.cpp | 10 +- test/unit/math/rev/prob/gamma_lccdf_test.cpp | 11 +- .../prob/incomplete_gamma_cdfs_reference.hpp | 370 ++++++++++++++++++ .../incomplete_gamma_cdfs_reference_test.cpp | 125 ++++++ 25 files changed, 1574 insertions(+), 1059 deletions(-) create mode 100644 test/unit/math/mix/prob/incomplete_gamma_cdfs_test.cpp create mode 100644 test/unit/math/mix/prob/inv_incomplete_gamma_cdfs_test.cpp create mode 100644 test/unit/math/rev/prob/incomplete_gamma_cdfs_reference.hpp create mode 100644 test/unit/math/rev/prob/incomplete_gamma_cdfs_reference_test.cpp diff --git a/stan/math/prim/fun/log_gamma_q_dgamma.hpp b/stan/math/prim/fun/log_gamma_q_dgamma.hpp index abfdfbeda89..5107c2a2edc 100644 --- a/stan/math/prim/fun/log_gamma_q_dgamma.hpp +++ b/stan/math/prim/fun/log_gamma_q_dgamma.hpp @@ -9,15 +9,18 @@ #include #include #include +#include #include #include #include #include +#include #include #include #include #include #include +#include #include namespace stan { @@ -148,6 +151,158 @@ log_q_gamma_cf_and_da(const T_a& a, const T_z& z, return {log_prefactor - log(f), d_log_prefactor - df / f}; } +/** + * Compute log(P(a,z)) and d/da log(P(a,z)) for z < a + 1 from the power + * series P(a,z) = z^a exp(-z) / Gamma(a + 1) * sum_k t_k, with t_0 = 1 and + * t_k = t_{k-1} z / (a + k), in log space. The derivative of log(t_k) with + * respect to a is carried with every term. The series stops when a bound + * on the remaining terms is below machine epsilon for the sum and for its + * derivative. + * + * @tparam T_a Type of shape parameter a (double or autodiff types) + * @tparam T_z Type of value parameter z (double or autodiff types) + * @param a Shape parameter + * @param z Value at which to evaluate, 0 < z < a + 1 + * @param max_steps Maximum number of terms + * @return pair of log(P(a,z)) and d/da log(P(a,z)) + */ +template +inline std::pair, return_type_t> +log_p_gamma_series_and_da(const T_a& a, const T_z& z, int max_steps = 100000) { + using T_return = return_type_t; + const double a_val = value_of_rec(a); + const double z_val = value_of_rec(z); + T_return t = 1.0; + T_return d_log_t = 0.0; + T_return sum = 1.0; + T_return d_sum = 0.0; + for (int k = 1; k <= max_steps; ++k) { + t *= z / (a + k); + d_log_t -= inv(a + k); + sum += t; + d_sum += t * d_log_t; + // the remaining terms sum to at most t q / (1 - q) + const double q = z_val / (a_val + k + 1); + if (q < 1) { + const double tail = value_of_rec(t) * q / (1 - q); + if (tail <= EPSILON * value_of_rec(sum) + && tail * (fabs(value_of_rec(d_log_t)) + 1) + <= EPSILON * fabs(value_of_rec(d_sum))) { + break; + } + } + } + const T_return log_z = log(z); + return {a * log_z - z - lgamma(a + 1.0) + log(sum), + log_z - digamma(a + 1.0) + d_sum / sum}; +} + +/** + * Result of log_gamma_tail: a log probability of the regularized incomplete + * gamma function and its derivatives with respect to a and z. + */ +template +struct log_gamma_tail_result { + T log_prob; + T d_a; + T d_z; +}; + +/** + * Compute log(P(a,z)) (lower tail) or log(Q(a,z)) (upper tail) with its + * derivatives with respect to a and z. Every result stays in log space, so + * it is finite where the probability underflows. + * + * For z > a + 1, log(Q) and its derivative come from the continued + * fraction, and log(P) = log1m_exp(log(Q)). Otherwise P comes from gamma_p; + * where P is below DBL_MIN / EPSILON, log(P) and its derivative come from + * the power series in log space. + * + * @tparam Upper true for log(Q), false for log(P) + * @tparam NeedDa compute the derivative with respect to a + * @tparam NeedDz compute the derivative with respect to z + * @tparam T_a Type of shape parameter a (double or autodiff types) + * @tparam T_z Type of value parameter z (double or autodiff types) + * @param a Shape parameter, a > 0 + * @param z Value at which to evaluate, 0 <= z < inf + * @return log probability, d/da and d/dz of it (0 where not computed) + */ +template +inline log_gamma_tail_result> log_gamma_tail( + const T_a& a, const T_z& z) { + using T_return = return_type_t; + const double a_val = value_of_rec(a); + const double z_val = value_of_rec(z); + log_gamma_tail_result out{0.0, 0.0, 0.0}; + if (z_val == 0) { + if constexpr (Upper) { + // Q(a, 0) = 1; d/dz log(Q) at 0 is minus the density at 0 + if constexpr (NeedDz) { + out.d_z = a_val > 1 ? 0.0 : (a_val == 1 ? -1.0 : NEGATIVE_INFTY); + } + } else { + out.log_prob = NEGATIVE_INFTY; + } + return out; + } + if (z_val == INFTY) { + // P(a, inf) = 1, Q(a, inf) = 0; the derivatives are treated as zero + if constexpr (Upper) { + out.log_prob = NEGATIVE_INFTY; + } + return out; + } + if (z_val > a_val + 1.0) { + T_return log_q; + T_return d_a_log_q = 0.0; + if constexpr (NeedDa) { + std::pair cf = log_q_gamma_cf_and_da(a, z); + log_q = cf.first; + d_a_log_q = cf.second; + } else { + log_q = log_q_gamma_cf(a, z); + } + if constexpr (Upper) { + out.log_prob = log_q; + out.d_a = d_a_log_q; + } else { + out.log_prob = log1m_exp(log_q); + if constexpr (NeedDa) { + out.d_a = -exp(log_q - out.log_prob) * d_a_log_q; + } + } + } else { + const T_return p = gamma_p(a, z); + if constexpr (Upper) { + out.log_prob = log1m(p); + if constexpr (NeedDa) { + out.d_a = -grad_reg_lower_inc_gamma(a, z) / (1.0 - p); + } + } else if (value_of_rec(p) + >= std::numeric_limits::min() / EPSILON) { + out.log_prob = log(p); + if constexpr (NeedDa) { + out.d_a = grad_reg_lower_inc_gamma(a, z) / p; + } + } else { + std::pair series = log_p_gamma_series_and_da(a, z); + out.log_prob = series.first; + out.d_a = series.second; + } + } + if constexpr (NeedDz) { + // the density of the gamma distribution at z, divided by P or Q + const T_return log_density = (a - 1.0) * log(z) - z - lgamma(a); + const T_return ratio = exp(log_density - out.log_prob); + if constexpr (Upper) { + out.d_z = -ratio; + } else { + out.d_z = ratio; + } + } + return out; +} + } // namespace internal /** diff --git a/stan/math/prim/prob/chi_square_cdf.hpp b/stan/math/prim/prob/chi_square_cdf.hpp index 46e0ea4e623..16acc4393fc 100644 --- a/stan/math/prim/prob/chi_square_cdf.hpp +++ b/stan/math/prim/prob/chi_square_cdf.hpp @@ -4,16 +4,14 @@ #include #include #include -#include #include -#include -#include +#include #include #include #include #include -#include #include +#include #include #include @@ -21,30 +19,29 @@ namespace stan { namespace math { /** \ingroup prob_dists - * Returns the chi square cumulative distribution function for the given - * variate and degrees of freedom. If given containers of matching sizes, - * returns the product of probabilities. + * The cumulative distribution function of the chi-square distribution. + * It is computed from the log of the regularized lower + * incomplete gamma function, so it stays finite where the probability + * underflows. * - * @tparam T_y type of scalar parameter - * @tparam T_dof type of degrees of freedom parameter - * @param y scalar parameter + * @tparam T_y type of y + * @tparam T_dof type of nu + * @param y random variable * @param nu degrees of freedom parameter - * @return probability or product of probabilities - * @throw std::domain_error if y is negative or nu is nonpositive + * @return probability, or product of probabilities for containers + * @throw std::domain_error if a parameter is out of its support * @throw std::invalid_argument if container sizes mismatch */ template inline return_type_t chi_square_cdf(const T_y& y, const T_dof& nu) { using T_partials_return = partials_return_t; - using std::exp; - using std::pow; using T_y_ref = ref_type_t; - using T_nu_ref = ref_type_t; + using T_dof_ref = ref_type_t; static constexpr const char* function = "chi_square_cdf"; check_consistent_sizes(function, "Random variable", y, "Degrees of freedom parameter", nu); T_y_ref y_ref = y; - T_nu_ref nu_ref = nu; + T_dof_ref nu_ref = nu; check_not_nan(function, "Random variable", y_ref); check_nonnegative(function, "Random variable", y_ref); check_positive_finite(function, "Degrees of freedom parameter", nu_ref); @@ -53,74 +50,55 @@ inline return_type_t chi_square_cdf(const T_y& y, const T_dof& nu) { return 1.0; } - T_partials_return cdf(1.0); + T_partials_return log_P(0.0); auto ops_partials = make_partials_propagator(y_ref, nu_ref); scalar_seq_view y_vec(y_ref); - scalar_seq_view nu_vec(nu_ref); - size_t N = max_size(y, nu); + scalar_seq_view nu_vec(nu_ref); + const size_t N = max_size(y, nu); - // Explicit return for extreme values - // The gradients are technically ill-defined, but treated as zero - for (size_t i = 0; i < stan::math::size(y); i++) { - if (y_vec.val(i) == 0) { + // An element with probability 0 (z = 0) makes the + // result 0; the gradients are treated as zero + for (size_t n = 0; n < N; n++) { + const double y_rec = value_of_rec(y_vec.val(n)); + if (0.5 * y_rec == 0) { return ops_partials.build(0.0); } } - VectorBuilder, T_partials_return, T_dof> gamma_vec( - math::size(nu)); - VectorBuilder, T_partials_return, T_dof> digamma_vec( - math::size(nu)); - - if constexpr (is_autodiff_v) { - for (size_t i = 0; i < stan::math::size(nu); i++) { - const T_partials_return alpha_dbl = nu_vec.val(i) * 0.5; - gamma_vec[i] = tgamma(alpha_dbl); - digamma_vec[i] = digamma(alpha_dbl); - } - } - + constexpr bool need_d_a = is_autodiff_v; + constexpr bool need_d_z = is_autodiff_v; for (size_t n = 0; n < N; n++) { - // Explicit results for extreme values - // The gradients are technically ill-defined, but treated as zero - if (y_vec.val(n) == INFTY) { + const T_partials_return y_val = y_vec.val(n); + const T_partials_return nu_val = nu_vec.val(n); + const T_partials_return z = 0.5 * y_val; + const double z_rec = value_of_rec(z); + // An element with probability 1 contributes nothing + if (z_rec == INFTY) { continue; } - - const T_partials_return y_dbl = y_vec.val(n); - const T_partials_return alpha_dbl = nu_vec.val(n) * 0.5; - const T_partials_return beta_dbl = 0.5; - - const T_partials_return Pn = gamma_p(alpha_dbl, beta_dbl * y_dbl); - - cdf *= Pn; - + const auto tail + = internal::log_gamma_tail(0.5 * nu_val, z); + log_P += tail.log_prob; if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[n] += beta_dbl * exp(-beta_dbl * y_dbl) - * pow(beta_dbl * y_dbl, alpha_dbl - 1) - / tgamma(alpha_dbl) / Pn; + partials<0>(ops_partials)[n] += 0.5 * tail.d_z; } if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] - -= 0.5 - * grad_reg_inc_gamma(alpha_dbl, beta_dbl * y_dbl, gamma_vec[n], - digamma_vec[n]) - / Pn; + partials<1>(ops_partials)[n] += 0.5 * tail.d_a; } } - + const T_partials_return P = exp(log_P); if constexpr (is_autodiff_v) { for (size_t n = 0; n < stan::math::size(y); ++n) { - partials<0>(ops_partials)[n] *= cdf; + partials<0>(ops_partials)[n] *= P; } } if constexpr (is_autodiff_v) { for (size_t n = 0; n < stan::math::size(nu); ++n) { - partials<1>(ops_partials)[n] *= cdf; + partials<1>(ops_partials)[n] *= P; } } - return ops_partials.build(cdf); + return ops_partials.build(P); } } // namespace math diff --git a/stan/math/prim/prob/chi_square_lccdf.hpp b/stan/math/prim/prob/chi_square_lccdf.hpp index 29da6302dc0..7926c175c96 100644 --- a/stan/math/prim/prob/chi_square_lccdf.hpp +++ b/stan/math/prim/prob/chi_square_lccdf.hpp @@ -4,17 +4,14 @@ #include #include #include -#include #include -#include -#include -#include +#include #include #include #include #include -#include #include +#include #include #include @@ -22,32 +19,31 @@ namespace stan { namespace math { /** \ingroup prob_dists - * Returns the chi square log complementary cumulative distribution - * function for the given variate and degrees of freedom. If given - * containers of matching sizes, returns the log sum of probabilities. + * The log complementary cumulative distribution function of the + * chi-square distribution. + * It is computed from the log of the regularized upper + * incomplete gamma function, so it stays finite where the probability + * underflows. * - * @tparam T_y type of scalar parameter - * @tparam T_dof type of degrees of freedom parameter - * @param y scalar parameter + * @tparam T_y type of y + * @tparam T_dof type of nu + * @param y random variable * @param nu degrees of freedom parameter - * @return log probability or log sum of probabilities - * @throw std::domain_error if y is negative or nu is nonpositive + * @return log probability, or sum of log probabilities for containers + * @throw std::domain_error if a parameter is out of its support * @throw std::invalid_argument if container sizes mismatch */ template inline return_type_t chi_square_lccdf(const T_y& y, const T_dof& nu) { using T_partials_return = partials_return_t; - using std::exp; - using std::log; - using std::pow; using T_y_ref = ref_type_t; - using T_nu_ref = ref_type_t; + using T_dof_ref = ref_type_t; static constexpr const char* function = "chi_square_lccdf"; check_consistent_sizes(function, "Random variable", y, "Degrees of freedom parameter", nu); T_y_ref y_ref = y; - T_nu_ref nu_ref = nu; + T_dof_ref nu_ref = nu; check_not_nan(function, "Random variable", y_ref); check_nonnegative(function, "Random variable", y_ref); check_positive_finite(function, "Degrees of freedom parameter", nu_ref); @@ -56,63 +52,44 @@ inline return_type_t chi_square_lccdf(const T_y& y, return 0; } - T_partials_return ccdf_log(0.0); + T_partials_return log_P(0.0); auto ops_partials = make_partials_propagator(y_ref, nu_ref); scalar_seq_view y_vec(y_ref); - scalar_seq_view nu_vec(nu_ref); - size_t N = max_size(y, nu); + scalar_seq_view nu_vec(nu_ref); + const size_t N = max_size(y, nu); - // Explicit return for extreme values - // The gradients are technically ill-defined, but treated as zero - for (size_t i = 0; i < stan::math::size(y); i++) { - if (y_vec.val(i) == 0) { - return ops_partials.build(0.0); - } - } - - VectorBuilder, T_partials_return, T_dof> gamma_vec( - math::size(nu)); - VectorBuilder, T_partials_return, T_dof> digamma_vec( - math::size(nu)); - - if constexpr (is_autodiff_v) { - for (size_t i = 0; i < stan::math::size(nu); i++) { - const T_partials_return alpha_dbl = nu_vec.val(i) * 0.5; - gamma_vec[i] = tgamma(alpha_dbl); - digamma_vec[i] = digamma(alpha_dbl); + // An element with probability 0 (z = inf) makes the + // result -inf; the gradients are treated as zero + for (size_t n = 0; n < N; n++) { + const double y_rec = value_of_rec(y_vec.val(n)); + if (0.5 * y_rec == INFTY) { + return ops_partials.build(negative_infinity()); } } + constexpr bool need_d_a = is_autodiff_v; + constexpr bool need_d_z = is_autodiff_v; for (size_t n = 0; n < N; n++) { - // Explicit results for extreme values - // The gradients are technically ill-defined, but treated as zero - if (y_vec.val(n) == INFTY) { - return ops_partials.build(negative_infinity()); + const T_partials_return y_val = y_vec.val(n); + const T_partials_return nu_val = nu_vec.val(n); + const T_partials_return z = 0.5 * y_val; + const double z_rec = value_of_rec(z); + // An element with probability 1 contributes nothing + if (z_rec == 0) { + continue; } - - const T_partials_return y_dbl = y_vec.val(n); - const T_partials_return alpha_dbl = nu_vec.val(n) * 0.5; - const T_partials_return beta_dbl = 0.5; - - const T_partials_return Pn = gamma_q(alpha_dbl, beta_dbl * y_dbl); - - ccdf_log += log(Pn); - + const auto tail + = internal::log_gamma_tail(0.5 * nu_val, z); + log_P += tail.log_prob; if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[n] -= beta_dbl * exp(-beta_dbl * y_dbl) - * pow(beta_dbl * y_dbl, alpha_dbl - 1) - / tgamma(alpha_dbl) / Pn; + partials<0>(ops_partials)[n] += 0.5 * tail.d_z; } if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] - += 0.5 - * grad_reg_inc_gamma(alpha_dbl, beta_dbl * y_dbl, gamma_vec[n], - digamma_vec[n]) - / Pn; + partials<1>(ops_partials)[n] += 0.5 * tail.d_a; } } - return ops_partials.build(ccdf_log); + return ops_partials.build(log_P); } } // namespace math diff --git a/stan/math/prim/prob/chi_square_lcdf.hpp b/stan/math/prim/prob/chi_square_lcdf.hpp index c75e14f8b14..d537fa3ca0a 100644 --- a/stan/math/prim/prob/chi_square_lcdf.hpp +++ b/stan/math/prim/prob/chi_square_lcdf.hpp @@ -4,17 +4,14 @@ #include #include #include -#include #include -#include -#include -#include +#include #include #include #include #include -#include #include +#include #include #include @@ -22,32 +19,31 @@ namespace stan { namespace math { /** \ingroup prob_dists - * Returns the chi square log cumulative distribution function for the given - * variate and degrees of freedom. If given containers of matching sizes, - * returns the log sum of probabilities. + * The log cumulative distribution function of the chi-square + * distribution. + * It is computed from the log of the regularized lower + * incomplete gamma function, so it stays finite where the probability + * underflows. * - * @tparam T_y type of scalar parameter - * @tparam T_dof type of degrees of freedom parameter - * @param y scalar parameter + * @tparam T_y type of y + * @tparam T_dof type of nu + * @param y random variable * @param nu degrees of freedom parameter - * @return log probability or log sum of probabilities - * @throw std::domain_error if y is negative or nu is nonpositive + * @return log probability, or sum of log probabilities for containers + * @throw std::domain_error if a parameter is out of its support * @throw std::invalid_argument if container sizes mismatch */ template inline return_type_t chi_square_lcdf(const T_y& y, const T_dof& nu) { using T_partials_return = partials_return_t; - using std::exp; - using std::log; - using std::pow; using T_y_ref = ref_type_t; - using T_nu_ref = ref_type_t; + using T_dof_ref = ref_type_t; static constexpr const char* function = "chi_square_lcdf"; check_consistent_sizes(function, "Random variable", y, "Degrees of freedom parameter", nu); T_y_ref y_ref = y; - T_nu_ref nu_ref = nu; + T_dof_ref nu_ref = nu; check_not_nan(function, "Random variable", y_ref); check_nonnegative(function, "Random variable", y_ref); check_positive_finite(function, "Degrees of freedom parameter", nu_ref); @@ -56,63 +52,44 @@ inline return_type_t chi_square_lcdf(const T_y& y, return 0; } - T_partials_return cdf_log(0.0); + T_partials_return log_P(0.0); auto ops_partials = make_partials_propagator(y_ref, nu_ref); scalar_seq_view y_vec(y_ref); - scalar_seq_view nu_vec(nu_ref); - size_t N = max_size(y, nu); + scalar_seq_view nu_vec(nu_ref); + const size_t N = max_size(y, nu); - // Explicit return for extreme values - // The gradients are technically ill-defined, but treated as zero - for (size_t i = 0; i < stan::math::size(y); i++) { - if (y_vec.val(i) == 0) { + // An element with probability 0 (z = 0) makes the + // result -inf; the gradients are treated as zero + for (size_t n = 0; n < N; n++) { + const double y_rec = value_of_rec(y_vec.val(n)); + if (0.5 * y_rec == 0) { return ops_partials.build(negative_infinity()); } } - VectorBuilder, T_partials_return, T_dof> gamma_vec( - math::size(nu)); - VectorBuilder, T_partials_return, T_dof> digamma_vec( - math::size(nu)); - - if constexpr (is_autodiff_v) { - for (size_t i = 0; i < stan::math::size(nu); i++) { - const T_partials_return alpha_dbl = nu_vec.val(i) * 0.5; - gamma_vec[i] = tgamma(alpha_dbl); - digamma_vec[i] = digamma(alpha_dbl); - } - } - + constexpr bool need_d_a = is_autodiff_v; + constexpr bool need_d_z = is_autodiff_v; for (size_t n = 0; n < N; n++) { - // Explicit results for extreme values - // The gradients are technically ill-defined, but treated as zero - if (y_vec.val(n) == INFTY) { - return ops_partials.build(0.0); + const T_partials_return y_val = y_vec.val(n); + const T_partials_return nu_val = nu_vec.val(n); + const T_partials_return z = 0.5 * y_val; + const double z_rec = value_of_rec(z); + // An element with probability 1 contributes nothing + if (z_rec == INFTY) { + continue; } - - const T_partials_return y_dbl = y_vec.val(n); - const T_partials_return alpha_dbl = nu_vec.val(n) * 0.5; - const T_partials_return beta_dbl = 0.5; - - const T_partials_return Pn = gamma_p(alpha_dbl, beta_dbl * y_dbl); - - cdf_log += log(Pn); - + const auto tail + = internal::log_gamma_tail(0.5 * nu_val, z); + log_P += tail.log_prob; if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[n] += beta_dbl * exp(-beta_dbl * y_dbl) - * pow(beta_dbl * y_dbl, alpha_dbl - 1) - / tgamma(alpha_dbl) / Pn; + partials<0>(ops_partials)[n] += 0.5 * tail.d_z; } if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] - -= 0.5 - * grad_reg_inc_gamma(alpha_dbl, beta_dbl * y_dbl, gamma_vec[n], - digamma_vec[n]) - / Pn; + partials<1>(ops_partials)[n] += 0.5 * tail.d_a; } } - return ops_partials.build(cdf_log); + return ops_partials.build(log_P); } } // namespace math diff --git a/stan/math/prim/prob/gamma_cdf.hpp b/stan/math/prim/prob/gamma_cdf.hpp index 4112a61ac39..ffdf3983237 100644 --- a/stan/math/prim/prob/gamma_cdf.hpp +++ b/stan/math/prim/prob/gamma_cdf.hpp @@ -4,19 +4,14 @@ #include #include #include -#include #include -#include -#include +#include #include -#include #include #include #include -#include #include -#include -#include +#include #include #include @@ -24,33 +19,34 @@ namespace stan { namespace math { /** \ingroup prob_dists - * The cumulative density function for a gamma distribution for y - * with the specified shape and inverse scale parameters. + * The cumulative distribution function of the gamma distribution. + * It is computed from the log of the regularized lower + * incomplete gamma function, so it stays finite where the probability + * underflows. * - * @tparam T_y type of scalar - * @tparam T_shape type of shape - * @tparam T_inv_scale type of inverse scale - * @param y A scalar variable. - * @param alpha Shape parameter. - * @param beta Inverse scale parameter. - * @throw std::domain_error if alpha is not greater than 0. - * @throw std::domain_error if beta is not greater than 0. - * @throw std::domain_error if y is not greater than or equal to 0. + * @tparam T_y type of y + * @tparam T_shape type of alpha + * @tparam T_inv_scale type of beta + * @param y random variable + * @param alpha shape parameter + * @param beta inverse scale parameter + * @return probability, or product of probabilities for containers + * @throw std::domain_error if a parameter is out of its support + * @throw std::invalid_argument if container sizes mismatch */ template inline return_type_t gamma_cdf( const T_y& y, const T_shape& alpha, const T_inv_scale& beta) { using T_partials_return = partials_return_t; using T_y_ref = ref_type_t; - using T_alpha_ref = ref_type_t; - using T_beta_ref = ref_type_t; - using std::exp; + using T_shape_ref = ref_type_t; + using T_inv_scale_ref = ref_type_t; static constexpr const char* function = "gamma_cdf"; check_consistent_sizes(function, "Random variable", y, "Shape parameter", alpha, "Inverse scale parameter", beta); T_y_ref y_ref = y; - T_alpha_ref alpha_ref = alpha; - T_beta_ref beta_ref = beta; + T_shape_ref alpha_ref = alpha; + T_inv_scale_ref beta_ref = beta; check_positive_finite(function, "Shape parameter", alpha_ref); check_positive_finite(function, "Inverse scale parameter", beta_ref); check_nonnegative(function, "Random variable", y_ref); @@ -59,71 +55,50 @@ inline return_type_t gamma_cdf( return 1.0; } - T_partials_return P(1.0); + T_partials_return log_P(0.0); auto ops_partials = make_partials_propagator(y_ref, alpha_ref, beta_ref); scalar_seq_view y_vec(y_ref); - scalar_seq_view alpha_vec(alpha_ref); - scalar_seq_view beta_vec(beta_ref); - size_t N = max_size(y, alpha, beta); + scalar_seq_view alpha_vec(alpha_ref); + scalar_seq_view beta_vec(beta_ref); + const size_t N = max_size(y, alpha, beta); - // Explicit return for extreme values - // The gradients are technically ill-defined, but treated as zero - for (size_t i = 0; i < stan::math::size(y); i++) { - if (y_vec.val(i) == 0) { + // An element with probability 0 (z = 0) makes the + // result 0; the gradients are treated as zero + for (size_t n = 0; n < N; n++) { + const double y_rec = value_of_rec(y_vec.val(n)); + const double beta_rec = value_of_rec(beta_vec.val(n)); + if (beta_rec * y_rec == 0) { return ops_partials.build(0.0); } } - using std::exp; - using std::pow; - - VectorBuilder, T_partials_return, T_shape> gamma_vec( - math::size(alpha)); - VectorBuilder, T_partials_return, T_shape> digamma_vec( - math::size(alpha)); - - if constexpr (is_autodiff_v) { - for (size_t i = 0; i < stan::math::size(alpha); i++) { - const T_partials_return alpha_dbl = alpha_vec.val(i); - gamma_vec[i] = tgamma(alpha_dbl); - digamma_vec[i] = digamma(alpha_dbl); - } - } - + constexpr bool need_d_a = is_autodiff_v; + constexpr bool need_d_z = is_any_autodiff_v; for (size_t n = 0; n < N; n++) { - // Explicit results for extreme values - // The gradients are technically ill-defined, but treated as zero - if (y_vec.val(n) == INFTY) { + const T_partials_return y_val = y_vec.val(n); + const T_partials_return alpha_val = alpha_vec.val(n); + const T_partials_return beta_val = beta_vec.val(n); + const T_partials_return z = beta_val * y_val; + const double z_rec = value_of_rec(z); + // An element with probability 1 contributes nothing + if (z_rec == INFTY) { continue; } - - const T_partials_return y_dbl = y_vec.val(n); - const T_partials_return alpha_dbl = alpha_vec.val(n); - const T_partials_return beta_dbl = beta_vec.val(n); - - const T_partials_return Pn = gamma_p(alpha_dbl, beta_dbl * y_dbl); - - P *= Pn; - + const auto tail + = internal::log_gamma_tail(alpha_val, z); + log_P += tail.log_prob; if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[n] += beta_dbl * exp(-beta_dbl * y_dbl) - * pow(beta_dbl * y_dbl, alpha_dbl - 1) - / tgamma(alpha_dbl) / Pn; + partials<0>(ops_partials)[n] += beta_val * tail.d_z; } if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] - -= grad_reg_inc_gamma(alpha_dbl, beta_dbl * y_dbl, gamma_vec[n], - digamma_vec[n]) - / Pn; + partials<1>(ops_partials)[n] += tail.d_a; } if constexpr (is_autodiff_v) { - partials<2>(ops_partials)[n] += y_dbl * exp(-beta_dbl * y_dbl) - * pow(beta_dbl * y_dbl, alpha_dbl - 1) - / tgamma(alpha_dbl) / Pn; + partials<2>(ops_partials)[n] += y_val * tail.d_z; } } - + const T_partials_return P = exp(log_P); if constexpr (is_autodiff_v) { for (size_t n = 0; n < stan::math::size(y); ++n) { partials<0>(ops_partials)[n] *= P; diff --git a/stan/math/prim/prob/gamma_lccdf.hpp b/stan/math/prim/prob/gamma_lccdf.hpp index f2db4ea25bb..50a37c5d4fe 100644 --- a/stan/math/prim/prob/gamma_lccdf.hpp +++ b/stan/math/prim/prob/gamma_lccdf.hpp @@ -4,101 +4,50 @@ #include #include #include -#include #include -#include -#include -#include -#include -#include -#include -#include +#include #include #include #include #include -#include +#include #include -#include #include #include -#include namespace stan { namespace math { -namespace internal { - -/** - * Computes log q and d(log q) / d(alpha) using continued fraction. - */ -template -inline std::optional, return_type_t>> -eval_q_cf(const T1& alpha, const T2& beta_y) { - using scalar_t = return_type_t; - using ret_t = std::pair; - if constexpr (!any_fvar && is_autodiff_v) { - std::pair log_q_result - = log_gamma_q_dgamma(value_of(alpha), value_of(beta_y)); - if (likely(std::isfinite(log_q_result.first))) { - return std::optional{log_q_result}; - } else { - return std::optional{std::nullopt}; - } - } else if constexpr (is_autodiff_v && !is_autodiff_v) { - // One continued fraction gives the value and the derivative. - const ret_t out = internal::log_q_gamma_cf_and_da(alpha, beta_y); - if (unlikely(!std::isfinite(out.first))) { - return std::optional{std::nullopt}; - } - return std::optional{out}; - } else { - ret_t out{internal::log_q_gamma_cf(alpha, beta_y), 0.0}; - if (unlikely(!std::isfinite(value_of_rec(out.first)))) { - return std::optional{std::nullopt}; - } - if constexpr (is_autodiff_v) { - // The partial is differentiated again, so it is computed with the - // derivative recurrence in the autodiff type. - out.second = internal::log_q_gamma_cf_and_da(alpha, beta_y).second; - } - return std::optional{out}; - } -} -/** - * Computes log q and d(log q) / d(alpha) using log1m. +/** \ingroup prob_dists + * The log complementary cumulative distribution function of the + * gamma distribution. + * It is computed from the log of the regularized upper + * incomplete gamma function, so it stays finite where the probability + * underflows. + * + * @tparam T_y type of y + * @tparam T_shape type of alpha + * @tparam T_inv_scale type of beta + * @param y random variable + * @param alpha shape parameter + * @param beta inverse scale parameter + * @return log probability, or sum of log probabilities for containers + * @throw std::domain_error if a parameter is out of its support + * @throw std::invalid_argument if container sizes mismatch */ -template -inline std::optional, return_type_t>> -eval_q_log1m(const T1& alpha, const T2& beta_y) { - using scalar_t = return_type_t; - using ret_t = std::pair; - ret_t out{log1m(gamma_p(alpha, beta_y)), 0.0}; - if (unlikely(!std::isfinite(value_of_rec(out.first)))) { - return std::optional{std::nullopt}; - } - if constexpr (is_autodiff_v) { - out.second = -grad_reg_lower_inc_gamma(alpha, beta_y) / exp(out.first); - } - return std::optional{out}; -} -} // namespace internal - template inline return_type_t gamma_lccdf( const T_y& y, const T_shape& alpha, const T_inv_scale& beta) { - using std::exp; - using std::log; using T_partials_return = partials_return_t; using T_y_ref = ref_type_t; - using T_alpha_ref = ref_type_t; - using T_beta_ref = ref_type_t; + using T_shape_ref = ref_type_t; + using T_inv_scale_ref = ref_type_t; static constexpr const char* function = "gamma_lccdf"; check_consistent_sizes(function, "Random variable", y, "Shape parameter", alpha, "Inverse scale parameter", beta); T_y_ref y_ref = y; - T_alpha_ref alpha_ref = alpha; - T_beta_ref beta_ref = beta; + T_shape_ref alpha_ref = alpha; + T_inv_scale_ref beta_ref = beta; check_positive_finite(function, "Shape parameter", alpha_ref); check_positive_finite(function, "Inverse scale parameter", beta_ref); check_nonnegative(function, "Random variable", y_ref); @@ -107,78 +56,52 @@ inline return_type_t gamma_lccdf( return 0; } - T_partials_return P(0.0); + T_partials_return log_P(0.0); auto ops_partials = make_partials_propagator(y_ref, alpha_ref, beta_ref); scalar_seq_view y_vec(y_ref); - scalar_seq_view alpha_vec(alpha_ref); - scalar_seq_view beta_vec(beta_ref); + scalar_seq_view alpha_vec(alpha_ref); + scalar_seq_view beta_vec(beta_ref); const size_t N = max_size(y, alpha, beta); - constexpr bool is_y_fvar = is_fvar_v>; - constexpr bool is_shape_fvar = is_fvar_v>; - constexpr bool is_beta_fvar = is_fvar_v>; - constexpr bool any_fvar = is_y_fvar || is_shape_fvar || is_beta_fvar; - + // An element with probability 0 (z = inf) makes the + // result -inf; the gradients are treated as zero for (size_t n = 0; n < N; n++) { - // Explicit results for extreme values - // The gradients are technically ill-defined, but treated as zero - const T_partials_return y_val = y_vec.val(n); - if (y_val == 0.0) { - continue; - } - if (y_val == INFTY) { + const double y_rec = value_of_rec(y_vec.val(n)); + const double beta_rec = value_of_rec(beta_vec.val(n)); + if (beta_rec * y_rec == INFTY) { return ops_partials.build(negative_infinity()); } + } + constexpr bool need_d_a = is_autodiff_v; + constexpr bool need_d_z = is_any_autodiff_v; + for (size_t n = 0; n < N; n++) { + const T_partials_return y_val = y_vec.val(n); const T_partials_return alpha_val = alpha_vec.val(n); const T_partials_return beta_val = beta_vec.val(n); - - const T_partials_return beta_y = beta_val * y_val; - if (beta_y == INFTY) { - return ops_partials.build(negative_infinity()); - } - std::optional> result; - if (beta_y > alpha_val + 1.0) { - result = internal::eval_q_cf(alpha_val, beta_y); - } else { - result = internal::eval_q_log1m(alpha_val, beta_y); - if (!result && beta_y > 0.0) { - // Fallback to continued fraction if log1m fails - result = internal::eval_q_cf(alpha_val, beta_y); - } - } - if (unlikely(!result)) { - return ops_partials.build(negative_infinity()); + const T_partials_return z = beta_val * y_val; + const double z_rec = value_of_rec(z); + // An element with probability 1 contributes nothing + if (z_rec == 0) { + continue; } - - P += result->first; - - if constexpr (is_autodiff_v || is_autodiff_v) { - const T_partials_return log_y = log(y_val); - const T_partials_return alpha_minus_one = fma(alpha_val, log_y, -log_y); - - const T_partials_return log_pdf = alpha_val * log(beta_val) - - lgamma(alpha_val) + alpha_minus_one - - beta_y; - - const T_partials_return hazard = exp(log_pdf - result->first); // f/Q - - if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[n] -= hazard; - } - if constexpr (is_autodiff_v) { - partials<2>(ops_partials)[n] -= (y_val / beta_val) * hazard; - } + const auto tail + = internal::log_gamma_tail(alpha_val, z); + log_P += tail.log_prob; + if constexpr (is_autodiff_v) { + partials<0>(ops_partials)[n] += beta_val * tail.d_z; } if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] += result->second; + partials<1>(ops_partials)[n] += tail.d_a; + } + if constexpr (is_autodiff_v) { + partials<2>(ops_partials)[n] += y_val * tail.d_z; } } - return ops_partials.build(P); + return ops_partials.build(log_P); } } // namespace math } // namespace stan - #endif diff --git a/stan/math/prim/prob/gamma_lcdf.hpp b/stan/math/prim/prob/gamma_lcdf.hpp index acac32bf92f..090959ff986 100644 --- a/stan/math/prim/prob/gamma_lcdf.hpp +++ b/stan/math/prim/prob/gamma_lcdf.hpp @@ -4,39 +4,49 @@ #include #include #include -#include #include -#include -#include -#include +#include #include #include #include #include -#include #include +#include #include #include namespace stan { namespace math { +/** \ingroup prob_dists + * The log cumulative distribution function of the gamma distribution. + * It is computed from the log of the regularized lower + * incomplete gamma function, so it stays finite where the probability + * underflows. + * + * @tparam T_y type of y + * @tparam T_shape type of alpha + * @tparam T_inv_scale type of beta + * @param y random variable + * @param alpha shape parameter + * @param beta inverse scale parameter + * @return log probability, or sum of log probabilities for containers + * @throw std::domain_error if a parameter is out of its support + * @throw std::invalid_argument if container sizes mismatch + */ template inline return_type_t gamma_lcdf( const T_y& y, const T_shape& alpha, const T_inv_scale& beta) { using T_partials_return = partials_return_t; - using std::exp; - using std::log; - using std::pow; using T_y_ref = ref_type_t; - using T_alpha_ref = ref_type_t; - using T_beta_ref = ref_type_t; + using T_shape_ref = ref_type_t; + using T_inv_scale_ref = ref_type_t; static constexpr const char* function = "gamma_lcdf"; check_consistent_sizes(function, "Random variable", y, "Shape parameter", alpha, "Inverse scale parameter", beta); T_y_ref y_ref = y; - T_alpha_ref alpha_ref = alpha; - T_beta_ref beta_ref = beta; + T_shape_ref alpha_ref = alpha; + T_inv_scale_ref beta_ref = beta; check_positive_finite(function, "Shape parameter", alpha_ref); check_positive_finite(function, "Inverse scale parameter", beta_ref); check_nonnegative(function, "Random variable", y_ref); @@ -45,60 +55,50 @@ inline return_type_t gamma_lcdf( return 0; } - T_partials_return P(0.0); + T_partials_return log_P(0.0); auto ops_partials = make_partials_propagator(y_ref, alpha_ref, beta_ref); scalar_seq_view y_vec(y_ref); - scalar_seq_view alpha_vec(alpha_ref); - scalar_seq_view beta_vec(beta_ref); - size_t N = max_size(y, alpha, beta); + scalar_seq_view alpha_vec(alpha_ref); + scalar_seq_view beta_vec(beta_ref); + const size_t N = max_size(y, alpha, beta); - // Explicit return for extreme values - // The gradients are technically ill-defined - for (size_t i = 0; i < stan::math::size(y); i++) { - if (y_vec.val(i) == 0) { + // An element with probability 0 (z = 0) makes the + // result -inf; the gradients are treated as zero + for (size_t n = 0; n < N; n++) { + const double y_rec = value_of_rec(y_vec.val(n)); + const double beta_rec = value_of_rec(beta_vec.val(n)); + if (beta_rec * y_rec == 0) { return ops_partials.build(negative_infinity()); } } + constexpr bool need_d_a = is_autodiff_v; + constexpr bool need_d_z = is_any_autodiff_v; for (size_t n = 0; n < N; n++) { - // Explicit results for extreme values - // The gradients are technically ill-defined, but treated as zero - if (y_vec.val(n) == INFTY) { - return ops_partials.build(0.0); + const T_partials_return y_val = y_vec.val(n); + const T_partials_return alpha_val = alpha_vec.val(n); + const T_partials_return beta_val = beta_vec.val(n); + const T_partials_return z = beta_val * y_val; + const double z_rec = value_of_rec(z); + // An element with probability 1 contributes nothing + if (z_rec == INFTY) { + continue; } - - const T_partials_return y_dbl = y_vec.val(n); - const T_partials_return alpha_dbl = alpha_vec.val(n); - const T_partials_return beta_dbl = beta_vec.val(n); - const T_partials_return log_beta_dbl = log(beta_dbl); - const T_partials_return beta_y_dbl = beta_dbl * y_dbl; - - const T_partials_return Pn = gamma_p(alpha_dbl, beta_y_dbl); - const T_partials_return log_Pn = log(Pn); - - P += log_Pn; - - if constexpr (is_any_autodiff_v) { - const T_partials_return log_y_dbl = log(y_dbl); - const T_partials_return d_num - = (-beta_y_dbl) + (alpha_dbl - 1) * (log_beta_dbl + log_y_dbl); - const T_partials_return d_den = lgamma(alpha_dbl) + log_Pn; - const T_partials_return d = exp(d_num - d_den); - - if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[n] += beta_dbl * d; - } - if constexpr (is_autodiff_v) { - partials<2>(ops_partials)[n] += y_dbl * d; - } + const auto tail + = internal::log_gamma_tail(alpha_val, z); + log_P += tail.log_prob; + if constexpr (is_autodiff_v) { + partials<0>(ops_partials)[n] += beta_val * tail.d_z; } if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] - += grad_reg_lower_inc_gamma(alpha_dbl, beta_y_dbl) / Pn; + partials<1>(ops_partials)[n] += tail.d_a; + } + if constexpr (is_autodiff_v) { + partials<2>(ops_partials)[n] += y_val * tail.d_z; } } - return ops_partials.build(P); + return ops_partials.build(log_P); } } // namespace math diff --git a/stan/math/prim/prob/inv_chi_square_cdf.hpp b/stan/math/prim/prob/inv_chi_square_cdf.hpp index d017fc1390f..e5d8bd04bb3 100644 --- a/stan/math/prim/prob/inv_chi_square_cdf.hpp +++ b/stan/math/prim/prob/inv_chi_square_cdf.hpp @@ -4,16 +4,14 @@ #include #include #include -#include #include -#include -#include +#include #include #include #include #include -#include #include +#include #include #include @@ -21,32 +19,31 @@ namespace stan { namespace math { /** \ingroup prob_dists - * Returns the inverse chi square cumulative distribution function for the - * given variate and degrees of freedom. If given containers of matching - * sizes, returns the product of probabilities. + * The cumulative distribution function of the inverse chi-square + * distribution. + * It is computed from the log of the regularized upper + * incomplete gamma function, so it stays finite where the probability + * underflows. * - * @tparam T_y type of scalar parameter - * @tparam T_dof type of degrees of freedom parameter - * @param y scalar parameter + * @tparam T_y type of y + * @tparam T_dof type of nu + * @param y random variable * @param nu degrees of freedom parameter - * @return probability or product of probabilities - * @throw std::domain_error if y is negative or nu is nonpositive + * @return probability, or product of probabilities for containers + * @throw std::domain_error if a parameter is out of its support * @throw std::invalid_argument if container sizes mismatch */ template inline return_type_t inv_chi_square_cdf(const T_y& y, const T_dof& nu) { using T_partials_return = partials_return_t; - using std::exp; - using std::pow; using T_y_ref = ref_type_t; - using T_nu_ref = ref_type_t; + using T_dof_ref = ref_type_t; static constexpr const char* function = "inv_chi_square_cdf"; check_consistent_sizes(function, "Random variable", y, "Degrees of freedom parameter", nu); - T_y_ref y_ref = y; - T_nu_ref nu_ref = nu; + T_dof_ref nu_ref = nu; check_positive_finite(function, "Degrees of freedom parameter", nu_ref); check_nonnegative(function, "Random variable", y_ref); @@ -54,64 +51,44 @@ inline return_type_t inv_chi_square_cdf(const T_y& y, return 1.0; } - T_partials_return P(1.0); + T_partials_return log_P(0.0); auto ops_partials = make_partials_propagator(y_ref, nu_ref); scalar_seq_view y_vec(y_ref); - scalar_seq_view nu_vec(nu_ref); - size_t N = max_size(y, nu); + scalar_seq_view nu_vec(nu_ref); + const size_t N = max_size(y, nu); - // Explicit return for extreme values - // The gradients are technically ill-defined, but treated as zero - for (size_t i = 0; i < stan::math::size(y); i++) { - if (y_vec.val(i) == 0) { + // An element with probability 0 (z = inf) makes the + // result 0; the gradients are treated as zero + for (size_t n = 0; n < N; n++) { + const double y_rec = value_of_rec(y_vec.val(n)); + if (0.5 / y_rec == INFTY) { return ops_partials.build(0.0); } } - VectorBuilder, T_partials_return, T_dof> gamma_vec( - math::size(nu)); - VectorBuilder, T_partials_return, T_dof> digamma_vec( - math::size(nu)); - - if constexpr (is_autodiff_v) { - for (size_t i = 0; i < stan::math::size(nu); i++) { - const T_partials_return nu_dbl = nu_vec.val(i); - gamma_vec[i] = tgamma(0.5 * nu_dbl); - digamma_vec[i] = digamma(0.5 * nu_dbl); - } - } - + constexpr bool need_d_a = is_autodiff_v; + constexpr bool need_d_z = is_autodiff_v; for (size_t n = 0; n < N; n++) { - // Explicit results for extreme values - // The gradients are technically ill-defined, but treated as zero - if (y_vec.val(n) == INFTY) { + const T_partials_return y_val = y_vec.val(n); + const T_partials_return nu_val = nu_vec.val(n); + const T_partials_return z = 0.5 / y_val; + const double z_rec = value_of_rec(z); + // An element with probability 1 contributes nothing + if (z_rec == 0) { continue; } - - const T_partials_return y_dbl = y_vec.val(n); - const T_partials_return y_inv_dbl = 1.0 / y_dbl; - const T_partials_return nu_dbl = nu_vec.val(n); - - const T_partials_return Pn = gamma_q(0.5 * nu_dbl, 0.5 * y_inv_dbl); - - P *= Pn; - + const auto tail + = internal::log_gamma_tail(0.5 * nu_val, z); + log_P += tail.log_prob; if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[n] += 0.5 * y_inv_dbl * y_inv_dbl - * exp(-0.5 * y_inv_dbl) - * pow(0.5 * y_inv_dbl, 0.5 * nu_dbl - 1) - / tgamma(0.5 * nu_dbl) / Pn; + partials<0>(ops_partials)[n] += -z / y_val * tail.d_z; } if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] - += 0.5 - * grad_reg_inc_gamma(0.5 * nu_dbl, 0.5 * y_inv_dbl, gamma_vec[n], - digamma_vec[n]) - / Pn; + partials<1>(ops_partials)[n] += 0.5 * tail.d_a; } } - + const T_partials_return P = exp(log_P); if constexpr (is_autodiff_v) { for (size_t n = 0; n < stan::math::size(y); ++n) { partials<0>(ops_partials)[n] *= P; diff --git a/stan/math/prim/prob/inv_chi_square_lccdf.hpp b/stan/math/prim/prob/inv_chi_square_lccdf.hpp index 27f73e38512..e588d6391c0 100644 --- a/stan/math/prim/prob/inv_chi_square_lccdf.hpp +++ b/stan/math/prim/prob/inv_chi_square_lccdf.hpp @@ -4,17 +4,14 @@ #include #include #include -#include #include -#include -#include -#include +#include #include #include #include #include -#include #include +#include #include #include @@ -22,33 +19,31 @@ namespace stan { namespace math { /** \ingroup prob_dists - * Returns the inverse chi square log complementary cumulative distribution - * function for the given variate and degrees of freedom. If given - * containers of matching sizes, returns the log sum of probabilities. + * The log complementary cumulative distribution function of the + * inverse chi-square distribution. + * It is computed from the log of the regularized lower + * incomplete gamma function, so it stays finite where the probability + * underflows. * - * @tparam T_y type of scalar parameter - * @tparam T_dof type of degrees of freedom parameter - * @param y scalar parameter + * @tparam T_y type of y + * @tparam T_dof type of nu + * @param y random variable * @param nu degrees of freedom parameter - * @return log probability or log sum of probabilities - * @throw std::domain_error if y is negative or nu is nonpositive + * @return log probability, or sum of log probabilities for containers + * @throw std::domain_error if a parameter is out of its support * @throw std::invalid_argument if container sizes mismatch */ template inline return_type_t inv_chi_square_lccdf(const T_y& y, const T_dof& nu) { using T_partials_return = partials_return_t; - using std::exp; - using std::log; - using std::pow; using T_y_ref = ref_type_t; - using T_nu_ref = ref_type_t; + using T_dof_ref = ref_type_t; static constexpr const char* function = "inv_chi_square_lccdf"; check_consistent_sizes(function, "Random variable", y, "Degrees of freedom parameter", nu); - T_y_ref y_ref = y; - T_nu_ref nu_ref = nu; + T_dof_ref nu_ref = nu; check_positive_finite(function, "Degrees of freedom parameter", nu_ref); check_nonnegative(function, "Random variable", y_ref); @@ -56,64 +51,44 @@ inline return_type_t inv_chi_square_lccdf(const T_y& y, return 0; } - T_partials_return P(0.0); + T_partials_return log_P(0.0); auto ops_partials = make_partials_propagator(y_ref, nu_ref); scalar_seq_view y_vec(y_ref); - scalar_seq_view nu_vec(nu_ref); - size_t N = max_size(y, nu); - - // Explicit return for extreme values - // The gradients are technically ill-defined, but treated as zero - for (size_t i = 0; i < stan::math::size(y); i++) { - if (y_vec.val(i) == 0) { - return ops_partials.build(0.0); - } - } - - VectorBuilder, T_partials_return, T_dof> gamma_vec( - math::size(nu)); - VectorBuilder, T_partials_return, T_dof> digamma_vec( - math::size(nu)); + scalar_seq_view nu_vec(nu_ref); + const size_t N = max_size(y, nu); - if constexpr (is_autodiff_v) { - for (size_t i = 0; i < stan::math::size(nu); i++) { - const T_partials_return nu_dbl = nu_vec.val(i); - gamma_vec[i] = tgamma(0.5 * nu_dbl); - digamma_vec[i] = digamma(0.5 * nu_dbl); + // An element with probability 0 (z = 0) makes the + // result -inf; the gradients are treated as zero + for (size_t n = 0; n < N; n++) { + const double y_rec = value_of_rec(y_vec.val(n)); + if (0.5 / y_rec == 0) { + return ops_partials.build(negative_infinity()); } } + constexpr bool need_d_a = is_autodiff_v; + constexpr bool need_d_z = is_autodiff_v; for (size_t n = 0; n < N; n++) { - // Explicit results for extreme values - // The gradients are technically ill-defined, but treated as zero - if (y_vec.val(n) == INFTY) { - return ops_partials.build(negative_infinity()); + const T_partials_return y_val = y_vec.val(n); + const T_partials_return nu_val = nu_vec.val(n); + const T_partials_return z = 0.5 / y_val; + const double z_rec = value_of_rec(z); + // An element with probability 1 contributes nothing + if (z_rec == INFTY) { + continue; } - - const T_partials_return y_dbl = y_vec.val(n); - const T_partials_return y_inv_dbl = 1.0 / y_dbl; - const T_partials_return nu_dbl = nu_vec.val(n); - - const T_partials_return Pn = gamma_p(0.5 * nu_dbl, 0.5 * y_inv_dbl); - - P += log(Pn); - + const auto tail + = internal::log_gamma_tail(0.5 * nu_val, z); + log_P += tail.log_prob; if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[n] -= 0.5 * y_inv_dbl * y_inv_dbl - * exp(-0.5 * y_inv_dbl) - * pow(0.5 * y_inv_dbl, 0.5 * nu_dbl - 1) - / tgamma(0.5 * nu_dbl) / Pn; + partials<0>(ops_partials)[n] += -z / y_val * tail.d_z; } if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] - -= 0.5 - * grad_reg_inc_gamma(0.5 * nu_dbl, 0.5 * y_inv_dbl, gamma_vec[n], - digamma_vec[n]) - / Pn; + partials<1>(ops_partials)[n] += 0.5 * tail.d_a; } } - return ops_partials.build(P); + return ops_partials.build(log_P); } } // namespace math diff --git a/stan/math/prim/prob/inv_chi_square_lcdf.hpp b/stan/math/prim/prob/inv_chi_square_lcdf.hpp index 86280b86ad2..e48c8e7fa46 100644 --- a/stan/math/prim/prob/inv_chi_square_lcdf.hpp +++ b/stan/math/prim/prob/inv_chi_square_lcdf.hpp @@ -4,17 +4,14 @@ #include #include #include -#include #include -#include -#include -#include +#include #include #include #include #include -#include #include +#include #include #include @@ -22,33 +19,31 @@ namespace stan { namespace math { /** \ingroup prob_dists - * Returns the inverse chi square log cumulative distribution function for - * the given variate and degrees of freedom. If given containers of - * matching sizes, returns the log sum of probabilities. + * The log cumulative distribution function of the inverse chi-square + * distribution. + * It is computed from the log of the regularized upper + * incomplete gamma function, so it stays finite where the probability + * underflows. * - * @tparam T_y type of scalar parameter - * @tparam T_dof type of degrees of freedom parameter - * @param y scalar parameter + * @tparam T_y type of y + * @tparam T_dof type of nu + * @param y random variable * @param nu degrees of freedom parameter - * @return log probability or log sum of probabilities - * @throw std::domain_error if y is negative or nu is nonpositive + * @return log probability, or sum of log probabilities for containers + * @throw std::domain_error if a parameter is out of its support * @throw std::invalid_argument if container sizes mismatch */ template inline return_type_t inv_chi_square_lcdf(const T_y& y, const T_dof& nu) { using T_partials_return = partials_return_t; - using std::exp; - using std::log; - using std::pow; using T_y_ref = ref_type_t; - using T_nu_ref = ref_type_t; + using T_dof_ref = ref_type_t; static constexpr const char* function = "inv_chi_square_lcdf"; check_consistent_sizes(function, "Random variable", y, "Degrees of freedom parameter", nu); - T_y_ref y_ref = y; - T_nu_ref nu_ref = nu; + T_dof_ref nu_ref = nu; check_positive_finite(function, "Degrees of freedom parameter", nu_ref); check_nonnegative(function, "Random variable", y_ref); @@ -56,64 +51,44 @@ inline return_type_t inv_chi_square_lcdf(const T_y& y, return 0; } - T_partials_return P(0.0); + T_partials_return log_P(0.0); auto ops_partials = make_partials_propagator(y_ref, nu_ref); scalar_seq_view y_vec(y_ref); - scalar_seq_view nu_vec(nu_ref); - size_t N = max_size(y, nu); + scalar_seq_view nu_vec(nu_ref); + const size_t N = max_size(y, nu); - // Explicit return for extreme values - // The gradients are technically ill-defined, but treated as zero - for (size_t i = 0; i < stan::math::size(y); i++) { - if (y_vec.val(i) == 0) { + // An element with probability 0 (z = inf) makes the + // result -inf; the gradients are treated as zero + for (size_t n = 0; n < N; n++) { + const double y_rec = value_of_rec(y_vec.val(n)); + if (0.5 / y_rec == INFTY) { return ops_partials.build(negative_infinity()); } } - VectorBuilder, T_partials_return, T_dof> gamma_vec( - math::size(nu)); - VectorBuilder, T_partials_return, T_dof> digamma_vec( - math::size(nu)); - - if constexpr (is_autodiff_v) { - for (size_t i = 0; i < stan::math::size(nu); i++) { - const T_partials_return nu_dbl = nu_vec.val(i); - gamma_vec[i] = tgamma(0.5 * nu_dbl); - digamma_vec[i] = digamma(0.5 * nu_dbl); - } - } - + constexpr bool need_d_a = is_autodiff_v; + constexpr bool need_d_z = is_autodiff_v; for (size_t n = 0; n < N; n++) { - // Explicit results for extreme values - // The gradients are technically ill-defined, but treated as zero - if (y_vec.val(n) == INFTY) { + const T_partials_return y_val = y_vec.val(n); + const T_partials_return nu_val = nu_vec.val(n); + const T_partials_return z = 0.5 / y_val; + const double z_rec = value_of_rec(z); + // An element with probability 1 contributes nothing + if (z_rec == 0) { continue; } - - const T_partials_return y_dbl = y_vec.val(n); - const T_partials_return y_inv_dbl = 1.0 / y_dbl; - const T_partials_return nu_dbl = nu_vec.val(n); - - const T_partials_return Pn = gamma_q(0.5 * nu_dbl, 0.5 * y_inv_dbl); - - P += log(Pn); - + const auto tail + = internal::log_gamma_tail(0.5 * nu_val, z); + log_P += tail.log_prob; if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[n] += 0.5 * y_inv_dbl * y_inv_dbl - * exp(-0.5 * y_inv_dbl) - * pow(0.5 * y_inv_dbl, 0.5 * nu_dbl - 1) - / tgamma(0.5 * nu_dbl) / Pn; + partials<0>(ops_partials)[n] += -z / y_val * tail.d_z; } if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] - += 0.5 - * grad_reg_inc_gamma(0.5 * nu_dbl, 0.5 * y_inv_dbl, gamma_vec[n], - digamma_vec[n]) - / Pn; + partials<1>(ops_partials)[n] += 0.5 * tail.d_a; } } - return ops_partials.build(P); + return ops_partials.build(log_P); } } // namespace math diff --git a/stan/math/prim/prob/inv_gamma_cdf.hpp b/stan/math/prim/prob/inv_gamma_cdf.hpp index 9c801521d03..414673b1993 100644 --- a/stan/math/prim/prob/inv_gamma_cdf.hpp +++ b/stan/math/prim/prob/inv_gamma_cdf.hpp @@ -4,16 +4,14 @@ #include #include #include -#include #include -#include -#include +#include #include #include #include #include -#include #include +#include #include #include @@ -21,39 +19,36 @@ namespace stan { namespace math { /** \ingroup prob_dists - * The CDF of an inverse gamma density for y with the specified - * shape and scale parameters. y, shape, and scale parameters must - * be greater than 0. + * The cumulative distribution function of the inverse gamma + * distribution. + * It is computed from the log of the regularized upper + * incomplete gamma function, so it stays finite where the probability + * underflows. * - * @tparam T_y type of scalar - * @tparam T_shape type of shape - * @tparam T_scale type of scale - * @param y A scalar variable. - * @param alpha Shape parameter. - * @param beta Scale parameter. - * @throw std::domain_error if alpha is not greater than 0. - * @throw std::domain_error if beta is not greater than 0. - * @throw std::domain_error if y is not greater than 0. + * @tparam T_y type of y + * @tparam T_shape type of alpha + * @tparam T_scale type of beta + * @param y random variable + * @param alpha shape parameter + * @param beta scale parameter + * @return probability, or product of probabilities for containers + * @throw std::domain_error if a parameter is out of its support + * @throw std::invalid_argument if container sizes mismatch */ - template inline return_type_t inv_gamma_cdf(const T_y& y, const T_shape& alpha, const T_scale& beta) { using T_partials_return = partials_return_t; using T_y_ref = ref_type_t; - using T_alpha_ref = ref_type_t; - using T_beta_ref = ref_type_t; - using std::exp; - using std::pow; + using T_shape_ref = ref_type_t; + using T_scale_ref = ref_type_t; static constexpr const char* function = "inv_gamma_cdf"; check_consistent_sizes(function, "Random variable", y, "Shape parameter", - alpha, "Scale Parameter", beta); - + alpha, "Scale parameter", beta); T_y_ref y_ref = y; - T_alpha_ref alpha_ref = alpha; - T_beta_ref beta_ref = beta; - + T_shape_ref alpha_ref = alpha; + T_scale_ref beta_ref = beta; check_positive_finite(function, "Shape parameter", alpha_ref); check_positive_finite(function, "Scale parameter", beta_ref); check_nonnegative(function, "Random variable", y_ref); @@ -62,70 +57,50 @@ inline return_type_t inv_gamma_cdf(const T_y& y, return 1.0; } - T_partials_return P(1.0); + T_partials_return log_P(0.0); auto ops_partials = make_partials_propagator(y_ref, alpha_ref, beta_ref); scalar_seq_view y_vec(y_ref); - scalar_seq_view alpha_vec(alpha_ref); - scalar_seq_view beta_vec(beta_ref); - size_t N = max_size(y, alpha, beta); + scalar_seq_view alpha_vec(alpha_ref); + scalar_seq_view beta_vec(beta_ref); + const size_t N = max_size(y, alpha, beta); - // Explicit return for extreme values - // The gradients are technically ill-defined, but treated as zero - for (size_t i = 0; i < stan::math::size(y); i++) { - if (y_vec.val(i) == 0) { + // An element with probability 0 (z = inf) makes the + // result 0; the gradients are treated as zero + for (size_t n = 0; n < N; n++) { + const double y_rec = value_of_rec(y_vec.val(n)); + const double beta_rec = value_of_rec(beta_vec.val(n)); + if (beta_rec / y_rec == INFTY) { return ops_partials.build(0.0); } } - VectorBuilder, T_partials_return, T_shape> gamma_vec( - math::size(alpha)); - VectorBuilder, T_partials_return, T_shape> digamma_vec( - math::size(alpha)); - - if constexpr (is_autodiff_v) { - for (size_t i = 0; i < stan::math::size(alpha); i++) { - const T_partials_return alpha_dbl = alpha_vec.val(i); - gamma_vec[i] = tgamma(alpha_dbl); - digamma_vec[i] = digamma(alpha_dbl); - } - } - + constexpr bool need_d_a = is_autodiff_v; + constexpr bool need_d_z = is_any_autodiff_v; for (size_t n = 0; n < N; n++) { - // Explicit results for extreme values - // The gradients are technically ill-defined, but treated as zero - if (y_vec.val(n) == INFTY) { + const T_partials_return y_val = y_vec.val(n); + const T_partials_return alpha_val = alpha_vec.val(n); + const T_partials_return beta_val = beta_vec.val(n); + const T_partials_return z = beta_val / y_val; + const double z_rec = value_of_rec(z); + // An element with probability 1 contributes nothing + if (z_rec == 0) { continue; } - - const T_partials_return y_dbl = y_vec.val(n); - const T_partials_return y_inv_dbl = 1.0 / y_dbl; - const T_partials_return alpha_dbl = alpha_vec.val(n); - const T_partials_return beta_dbl = beta_vec.val(n); - - const T_partials_return Pn = gamma_q(alpha_dbl, beta_dbl * y_inv_dbl); - - P *= Pn; - + const auto tail + = internal::log_gamma_tail(alpha_val, z); + log_P += tail.log_prob; if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[n] += beta_dbl * y_inv_dbl * y_inv_dbl - * exp(-beta_dbl * y_inv_dbl) - * pow(beta_dbl * y_inv_dbl, alpha_dbl - 1) - / tgamma(alpha_dbl) / Pn; + partials<0>(ops_partials)[n] += -z / y_val * tail.d_z; } if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] - += grad_reg_inc_gamma(alpha_dbl, beta_dbl * y_inv_dbl, gamma_vec[n], - digamma_vec[n]) - / Pn; + partials<1>(ops_partials)[n] += tail.d_a; } if constexpr (is_autodiff_v) { - partials<2>(ops_partials)[n] += -y_inv_dbl * exp(-beta_dbl * y_inv_dbl) - * pow(beta_dbl * y_inv_dbl, alpha_dbl - 1) - / tgamma(alpha_dbl) / Pn; + partials<2>(ops_partials)[n] += tail.d_z / y_val; } } - + const T_partials_return P = exp(log_P); if constexpr (is_autodiff_v) { for (size_t n = 0; n < stan::math::size(y); ++n) { partials<0>(ops_partials)[n] *= P; diff --git a/stan/math/prim/prob/inv_gamma_lccdf.hpp b/stan/math/prim/prob/inv_gamma_lccdf.hpp index 55d4485c8eb..48a0a5c4007 100644 --- a/stan/math/prim/prob/inv_gamma_lccdf.hpp +++ b/stan/math/prim/prob/inv_gamma_lccdf.hpp @@ -4,41 +4,50 @@ #include #include #include -#include #include -#include -#include -#include +#include #include #include #include #include -#include #include +#include #include #include namespace stan { namespace math { +/** \ingroup prob_dists + * The log complementary cumulative distribution function of the + * inverse gamma distribution. + * It is computed from the log of the regularized lower + * incomplete gamma function, so it stays finite where the probability + * underflows. + * + * @tparam T_y type of y + * @tparam T_shape type of alpha + * @tparam T_scale type of beta + * @param y random variable + * @param alpha shape parameter + * @param beta scale parameter + * @return log probability, or sum of log probabilities for containers + * @throw std::domain_error if a parameter is out of its support + * @throw std::invalid_argument if container sizes mismatch + */ template inline return_type_t inv_gamma_lccdf( const T_y& y, const T_shape& alpha, const T_scale& beta) { using T_partials_return = partials_return_t; using T_y_ref = ref_type_t; - using T_alpha_ref = ref_type_t; - using T_beta_ref = ref_type_t; - using std::exp; - using std::log; - using std::pow; + using T_shape_ref = ref_type_t; + using T_scale_ref = ref_type_t; static constexpr const char* function = "inv_gamma_lccdf"; check_consistent_sizes(function, "Random variable", y, "Shape parameter", - alpha, "Scale Parameter", beta); - + alpha, "Scale parameter", beta); T_y_ref y_ref = y; - T_alpha_ref alpha_ref = alpha; - T_beta_ref beta_ref = beta; - + T_shape_ref alpha_ref = alpha; + T_scale_ref beta_ref = beta; check_positive_finite(function, "Shape parameter", alpha_ref); check_positive_finite(function, "Scale parameter", beta_ref); check_nonnegative(function, "Random variable", y_ref); @@ -47,70 +56,50 @@ inline return_type_t inv_gamma_lccdf( return 0; } - T_partials_return P(0.0); + T_partials_return log_P(0.0); auto ops_partials = make_partials_propagator(y_ref, alpha_ref, beta_ref); scalar_seq_view y_vec(y_ref); - scalar_seq_view alpha_vec(alpha_ref); - scalar_seq_view beta_vec(beta_ref); - size_t N = max_size(y, alpha, beta); + scalar_seq_view alpha_vec(alpha_ref); + scalar_seq_view beta_vec(beta_ref); + const size_t N = max_size(y, alpha, beta); - // Explicit return for extreme values - // The gradients are technically ill-defined, but treated as zero - for (size_t i = 0; i < stan::math::size(y); i++) { - if (y_vec.val(i) == 0) { - return ops_partials.build(0.0); - } - } - - VectorBuilder, T_partials_return, T_shape> gamma_vec( - math::size(alpha)); - VectorBuilder, T_partials_return, T_shape> digamma_vec( - math::size(alpha)); - - if constexpr (is_autodiff_v) { - for (size_t i = 0; i < stan::math::size(alpha); i++) { - const T_partials_return alpha_dbl = alpha_vec.val(i); - gamma_vec[i] = tgamma(alpha_dbl); - digamma_vec[i] = digamma(alpha_dbl); + // An element with probability 0 (z = 0) makes the + // result -inf; the gradients are treated as zero + for (size_t n = 0; n < N; n++) { + const double y_rec = value_of_rec(y_vec.val(n)); + const double beta_rec = value_of_rec(beta_vec.val(n)); + if (beta_rec / y_rec == 0) { + return ops_partials.build(negative_infinity()); } } + constexpr bool need_d_a = is_autodiff_v; + constexpr bool need_d_z = is_any_autodiff_v; for (size_t n = 0; n < N; n++) { - // Explicit results for extreme values - // The gradients are technically ill-defined, but treated as zero - if (y_vec.val(n) == INFTY) { - return ops_partials.build(negative_infinity()); + const T_partials_return y_val = y_vec.val(n); + const T_partials_return alpha_val = alpha_vec.val(n); + const T_partials_return beta_val = beta_vec.val(n); + const T_partials_return z = beta_val / y_val; + const double z_rec = value_of_rec(z); + // An element with probability 1 contributes nothing + if (z_rec == INFTY) { + continue; } - - const T_partials_return y_dbl = y_vec.val(n); - const T_partials_return y_inv_dbl = 1.0 / y_dbl; - const T_partials_return alpha_dbl = alpha_vec.val(n); - const T_partials_return beta_dbl = beta_vec.val(n); - - const T_partials_return Pn = gamma_p(alpha_dbl, beta_dbl * y_inv_dbl); - - P += log(Pn); - + const auto tail + = internal::log_gamma_tail(alpha_val, z); + log_P += tail.log_prob; if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[n] -= beta_dbl * y_inv_dbl * y_inv_dbl - * exp(-beta_dbl * y_inv_dbl) - * pow(beta_dbl * y_inv_dbl, alpha_dbl - 1) - / tgamma(alpha_dbl) / Pn; + partials<0>(ops_partials)[n] += -z / y_val * tail.d_z; } if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] - -= grad_reg_inc_gamma(alpha_dbl, beta_dbl * y_inv_dbl, gamma_vec[n], - digamma_vec[n]) - / Pn; + partials<1>(ops_partials)[n] += tail.d_a; } if constexpr (is_autodiff_v) { - partials<2>(ops_partials)[n] += y_inv_dbl * exp(-beta_dbl * y_inv_dbl) - * pow(beta_dbl * y_inv_dbl, alpha_dbl - 1) - / tgamma(alpha_dbl) / Pn; + partials<2>(ops_partials)[n] += tail.d_z / y_val; } } - return ops_partials.build(P); + return ops_partials.build(log_P); } } // namespace math diff --git a/stan/math/prim/prob/inv_gamma_lcdf.hpp b/stan/math/prim/prob/inv_gamma_lcdf.hpp index 5e0a78c192d..37b7ace53f4 100644 --- a/stan/math/prim/prob/inv_gamma_lcdf.hpp +++ b/stan/math/prim/prob/inv_gamma_lcdf.hpp @@ -4,41 +4,50 @@ #include #include #include -#include #include -#include -#include -#include +#include #include #include #include #include -#include #include +#include #include #include namespace stan { namespace math { +/** \ingroup prob_dists + * The log cumulative distribution function of the inverse gamma + * distribution. + * It is computed from the log of the regularized upper + * incomplete gamma function, so it stays finite where the probability + * underflows. + * + * @tparam T_y type of y + * @tparam T_shape type of alpha + * @tparam T_scale type of beta + * @param y random variable + * @param alpha shape parameter + * @param beta scale parameter + * @return log probability, or sum of log probabilities for containers + * @throw std::domain_error if a parameter is out of its support + * @throw std::invalid_argument if container sizes mismatch + */ template inline return_type_t inv_gamma_lcdf( const T_y& y, const T_shape& alpha, const T_scale& beta) { using T_partials_return = partials_return_t; using T_y_ref = ref_type_t; - using T_alpha_ref = ref_type_t; - using T_beta_ref = ref_type_t; - using std::exp; - using std::log; - using std::pow; + using T_shape_ref = ref_type_t; + using T_scale_ref = ref_type_t; static constexpr const char* function = "inv_gamma_lcdf"; check_consistent_sizes(function, "Random variable", y, "Shape parameter", - alpha, "Scale Parameter", beta); - + alpha, "Scale parameter", beta); T_y_ref y_ref = y; - T_alpha_ref alpha_ref = alpha; - T_beta_ref beta_ref = beta; - + T_shape_ref alpha_ref = alpha; + T_scale_ref beta_ref = beta; check_positive_finite(function, "Shape parameter", alpha_ref); check_positive_finite(function, "Scale parameter", beta_ref); check_nonnegative(function, "Random variable", y_ref); @@ -47,70 +56,50 @@ inline return_type_t inv_gamma_lcdf( return 0; } - T_partials_return P(0.0); + T_partials_return log_P(0.0); auto ops_partials = make_partials_propagator(y_ref, alpha_ref, beta_ref); scalar_seq_view y_vec(y_ref); - scalar_seq_view alpha_vec(alpha_ref); - scalar_seq_view beta_vec(beta_ref); - size_t N = max_size(y, alpha, beta); + scalar_seq_view alpha_vec(alpha_ref); + scalar_seq_view beta_vec(beta_ref); + const size_t N = max_size(y, alpha, beta); - // Explicit return for extreme values - // The gradients are technically ill-defined, but treated as zero - for (size_t i = 0; i < stan::math::size(y); i++) { - if (y_vec.val(i) == 0) { + // An element with probability 0 (z = inf) makes the + // result -inf; the gradients are treated as zero + for (size_t n = 0; n < N; n++) { + const double y_rec = value_of_rec(y_vec.val(n)); + const double beta_rec = value_of_rec(beta_vec.val(n)); + if (beta_rec / y_rec == INFTY) { return ops_partials.build(negative_infinity()); } } - VectorBuilder, T_partials_return, T_shape> gamma_vec( - math::size(alpha)); - VectorBuilder, T_partials_return, T_shape> digamma_vec( - math::size(alpha)); - - if constexpr (is_autodiff_v) { - for (size_t i = 0; i < stan::math::size(alpha); i++) { - const T_partials_return alpha_dbl = alpha_vec.val(i); - gamma_vec[i] = tgamma(alpha_dbl); - digamma_vec[i] = digamma(alpha_dbl); - } - } - + constexpr bool need_d_a = is_autodiff_v; + constexpr bool need_d_z = is_any_autodiff_v; for (size_t n = 0; n < N; n++) { - // Explicit results for extreme values - // The gradients are technically ill-defined, but treated as zero - if (y_vec.val(n) == INFTY) { + const T_partials_return y_val = y_vec.val(n); + const T_partials_return alpha_val = alpha_vec.val(n); + const T_partials_return beta_val = beta_vec.val(n); + const T_partials_return z = beta_val / y_val; + const double z_rec = value_of_rec(z); + // An element with probability 1 contributes nothing + if (z_rec == 0) { continue; } - - const T_partials_return y_dbl = y_vec.val(n); - const T_partials_return y_inv_dbl = 1.0 / y_dbl; - const T_partials_return alpha_dbl = alpha_vec.val(n); - const T_partials_return beta_dbl = beta_vec.val(n); - - const T_partials_return Pn = gamma_q(alpha_dbl, beta_dbl * y_inv_dbl); - - P += log(Pn); - + const auto tail + = internal::log_gamma_tail(alpha_val, z); + log_P += tail.log_prob; if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[n] += beta_dbl * y_inv_dbl * y_inv_dbl - * exp(-beta_dbl * y_inv_dbl) - * pow(beta_dbl * y_inv_dbl, alpha_dbl - 1) - / tgamma(alpha_dbl) / Pn; + partials<0>(ops_partials)[n] += -z / y_val * tail.d_z; } if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] - += grad_reg_inc_gamma(alpha_dbl, beta_dbl * y_inv_dbl, gamma_vec[n], - digamma_vec[n]) - / Pn; + partials<1>(ops_partials)[n] += tail.d_a; } if constexpr (is_autodiff_v) { - partials<2>(ops_partials)[n] += -y_inv_dbl * exp(-beta_dbl * y_inv_dbl) - * pow(beta_dbl * y_inv_dbl, alpha_dbl - 1) - / tgamma(alpha_dbl) / Pn; + partials<2>(ops_partials)[n] += tail.d_z / y_val; } } - return ops_partials.build(P); + return ops_partials.build(log_P); } } // namespace math diff --git a/stan/math/prim/prob/poisson_cdf.hpp b/stan/math/prim/prob/poisson_cdf.hpp index c6ba4b866e0..a23e0ee195c 100644 --- a/stan/math/prim/prob/poisson_cdf.hpp +++ b/stan/math/prim/prob/poisson_cdf.hpp @@ -3,18 +3,19 @@ #include #include -#include -#include #include +#include #include -#include +#include #include #include +#include #include #include -#include +#include #include #include +#include #include #include #include @@ -22,13 +23,25 @@ namespace stan { namespace math { -// Poisson CDF +/** \ingroup prob_dists + * The cumulative distribution function of the Poisson distribution, + * Q(n + 1, lambda) with the regularized upper incomplete gamma function + * Q. It is computed from log Q, so its gradient stays finite where the + * probability underflows. + * + * @tparam T_n type of n + * @tparam T_rate type of lambda + * @param n random variable + * @param lambda rate parameter + * @return probability, or product of probabilities for containers + * @throw std::domain_error if lambda is negative + * @throw std::invalid_argument if container sizes mismatch + */ template inline return_type_t poisson_cdf(const T_n& n, const T_rate& lambda) { using T_partials_return = partials_return_t; using T_n_ref = ref_type_if_not_constant_t; using T_lambda_ref = ref_type_if_not_constant_t; - using std::pow; static constexpr const char* function = "poisson_cdf"; check_consistent_sizes(function, "Random variable", n, "Rate parameter", lambda); @@ -52,16 +65,40 @@ inline return_type_t poisson_cdf(const T_n& n, const T_rate& lambda) { return ops_partials.build(0.0); } - const auto& Pi - = to_ref_if>(gamma_q(n_val + 1.0, lambda_val)); + scalar_seq_view n_vec(n_ref); + scalar_seq_view lambda_vec(lambda_ref); + const size_t N = max_size(n, lambda); - T_partials_return P = prod(Pi); + // An element with probability 0 (lambda = inf) makes the result 0; the + // gradients are treated as zero + for (size_t i = 0; i < N; i++) { + if (value_of_rec(lambda_vec.val(i)) == INFTY) { + return ops_partials.build(0.0); + } + } + T_partials_return log_P(0.0); + for (size_t i = 0; i < N; i++) { + // An element with probability 1 (n = INT_MAX, which stands for an + // infinite count) contributes nothing + if (n_vec.val(i) == std::numeric_limits::max()) { + continue; + } + const T_partials_return lambda_dbl = lambda_vec.val(i); + const auto tail + = internal::log_gamma_tail>( + n_vec.val(i) + 1.0, lambda_dbl); + log_P += tail.log_prob; + if constexpr (is_autodiff_v) { + partials<0>(ops_partials)[i] += tail.d_z; + } + } + const T_partials_return P = exp(log_P); if constexpr (is_autodiff_v) { - partials<0>(ops_partials) = -exp(-lambda_val) * pow(lambda_val, n_val) - / (tgamma(n_val + 1.0) * Pi) * P; + for (size_t i = 0; i < stan::math::size(lambda); ++i) { + partials<0>(ops_partials)[i] *= P; + } } - return ops_partials.build(P); } diff --git a/stan/math/prim/prob/poisson_lccdf.hpp b/stan/math/prim/prob/poisson_lccdf.hpp index c3a622cfb9d..b5248360af1 100644 --- a/stan/math/prim/prob/poisson_lccdf.hpp +++ b/stan/math/prim/prob/poisson_lccdf.hpp @@ -3,20 +3,17 @@ #include #include -#include -#include #include #include -#include -#include -#include -#include +#include #include #include -#include +#include #include +#include #include #include +#include #include #include #include @@ -24,6 +21,20 @@ namespace stan { namespace math { +/** \ingroup prob_dists + * The log complementary cumulative distribution function of the Poisson + * distribution, log P(n + 1, lambda) with the regularized lower incomplete + * gamma function P. It is computed in log space, so it stays finite where + * the probability underflows. + * + * @tparam T_n type of n + * @tparam T_rate type of lambda + * @param n random variable + * @param lambda rate parameter + * @return log probability, or sum of log probabilities for containers + * @throw std::domain_error if lambda is negative + * @throw std::invalid_argument if container sizes mismatch + */ template inline return_type_t poisson_lccdf(const T_n& n, const T_rate& lambda) { using T_partials_return = partials_return_t; @@ -52,16 +63,32 @@ inline return_type_t poisson_lccdf(const T_n& n, const T_rate& lambda) { return ops_partials.build(0.0); } - const auto& log_Pi - = to_ref_if>(log(gamma_p(n_val + 1.0, lambda_val))); - T_partials_return P = sum(log_Pi); + scalar_seq_view n_vec(n_ref); + scalar_seq_view lambda_vec(lambda_ref); + const size_t N = max_size(n, lambda); - if constexpr (is_autodiff_v) { - partials<0>(ops_partials) = exp(n_val * log(lambda_val) - lambda_val - - lgamma(n_val + 1.0) - log_Pi); + // An element with probability 0 (lambda = 0, or n = INT_MAX, which + // stands for an infinite count) makes the result -inf; the gradients are + // treated as zero + for (size_t i = 0; i < N; i++) { + if (value_of_rec(lambda_vec.val(i)) == 0 + || n_vec.val(i) == std::numeric_limits::max()) { + return ops_partials.build(negative_infinity()); + } } - return ops_partials.build(P); + T_partials_return log_P(0.0); + for (size_t i = 0; i < N; i++) { + const T_partials_return lambda_dbl = lambda_vec.val(i); + const auto tail + = internal::log_gamma_tail>( + n_vec.val(i) + 1.0, lambda_dbl); + log_P += tail.log_prob; + if constexpr (is_autodiff_v) { + partials<0>(ops_partials)[i] += tail.d_z; + } + } + return ops_partials.build(log_P); } } // namespace math diff --git a/stan/math/prim/prob/poisson_lcdf.hpp b/stan/math/prim/prob/poisson_lcdf.hpp index 31baff192c3..2889ca6765f 100644 --- a/stan/math/prim/prob/poisson_lcdf.hpp +++ b/stan/math/prim/prob/poisson_lcdf.hpp @@ -3,20 +3,17 @@ #include #include -#include -#include #include #include -#include -#include -#include -#include +#include #include #include -#include +#include #include +#include #include #include +#include #include #include #include @@ -24,6 +21,20 @@ namespace stan { namespace math { +/** \ingroup prob_dists + * The log cumulative distribution function of the Poisson distribution, + * log Q(n + 1, lambda) with the regularized upper incomplete gamma + * function Q. It is computed in log space, so it stays finite where the + * probability underflows. + * + * @tparam T_n type of n + * @tparam T_rate type of lambda + * @param n random variable + * @param lambda rate parameter + * @return log probability, or sum of log probabilities for containers + * @throw std::domain_error if lambda is negative + * @throw std::invalid_argument if container sizes mismatch + */ template inline return_type_t poisson_lcdf(const T_n& n, const T_rate& lambda) { using T_partials_return = partials_return_t; @@ -52,16 +63,35 @@ inline return_type_t poisson_lcdf(const T_n& n, const T_rate& lambda) { return ops_partials.build(negative_infinity()); } - const auto& log_Pi - = to_ref_if>(log(gamma_q(n_val + 1.0, lambda_val))); - T_partials_return P = sum(log_Pi); + scalar_seq_view n_vec(n_ref); + scalar_seq_view lambda_vec(lambda_ref); + const size_t N = max_size(n, lambda); - if constexpr (is_autodiff_v) { - partials<0>(ops_partials) = -exp(n_val * log(lambda_val) - lambda_val - - lgamma(n_val + 1.0) - log_Pi); + // An element with probability 0 (lambda = inf) makes the result -inf; + // the gradients are treated as zero + for (size_t i = 0; i < N; i++) { + if (value_of_rec(lambda_vec.val(i)) == INFTY) { + return ops_partials.build(negative_infinity()); + } } - return ops_partials.build(P); + T_partials_return log_P(0.0); + for (size_t i = 0; i < N; i++) { + // An element with probability 1 (n = INT_MAX, which stands for an + // infinite count) contributes nothing + if (n_vec.val(i) == std::numeric_limits::max()) { + continue; + } + const T_partials_return lambda_dbl = lambda_vec.val(i); + const auto tail + = internal::log_gamma_tail>( + n_vec.val(i) + 1.0, lambda_dbl); + log_P += tail.log_prob; + if constexpr (is_autodiff_v) { + partials<0>(ops_partials)[i] += tail.d_z; + } + } + return ops_partials.build(log_P); } } // namespace math diff --git a/stan/math/prim/prob/scaled_inv_chi_square_cdf.hpp b/stan/math/prim/prob/scaled_inv_chi_square_cdf.hpp index 1ba20033d68..48a62dee6dc 100644 --- a/stan/math/prim/prob/scaled_inv_chi_square_cdf.hpp +++ b/stan/math/prim/prob/scaled_inv_chi_square_cdf.hpp @@ -4,17 +4,14 @@ #include #include #include -#include #include -#include -#include +#include #include #include #include #include -#include -#include #include +#include #include #include @@ -22,35 +19,36 @@ namespace stan { namespace math { /** \ingroup prob_dists - * The CDF of a scaled inverse chi-squared density for y with the - * specified degrees of freedom parameter and scale parameter. + * The cumulative distribution function of the scaled inverse + * chi-square distribution. + * It is computed from the log of the regularized upper + * incomplete gamma function, so it stays finite where the probability + * underflows. * - * @tparam T_y type of scalar - * @tparam T_dof type of degrees of freedom - * @param y A scalar variable. - * @param nu Degrees of freedom. - * @param s Scale parameter. - * @throw std::domain_error if nu is not greater than 0 - * @throw std::domain_error if s is not greater than 0. - * @throw std::domain_error if y is not greater than 0. + * @tparam T_y type of y + * @tparam T_dof type of nu + * @tparam T_scale type of s + * @param y random variable + * @param nu degrees of freedom parameter + * @param s scale parameter + * @return probability, or product of probabilities for containers + * @throw std::domain_error if a parameter is out of its support + * @throw std::invalid_argument if container sizes mismatch */ template inline return_type_t scaled_inv_chi_square_cdf( const T_y& y, const T_dof& nu, const T_scale& s) { using T_partials_return = partials_return_t; - using std::exp; - using std::pow; using T_y_ref = ref_type_t; - using T_nu_ref = ref_type_t; - using T_s_ref = ref_type_t; + using T_dof_ref = ref_type_t; + using T_scale_ref = ref_type_t; static constexpr const char* function = "scaled_inv_chi_square_cdf"; check_consistent_sizes(function, "Random variable", y, "Degrees of freedom parameter", nu, "Scale parameter", s); T_y_ref y_ref = y; - T_nu_ref nu_ref = nu; - T_s_ref s_ref = s; - + T_dof_ref nu_ref = nu; + T_scale_ref s_ref = s; check_nonnegative(function, "Random variable", y_ref); check_positive_finite(function, "Degrees of freedom parameter", nu_ref); check_positive_finite(function, "Scale parameter", s_ref); @@ -59,77 +57,52 @@ inline return_type_t scaled_inv_chi_square_cdf( return 1.0; } - T_partials_return P(1.0); + T_partials_return log_P(0.0); auto ops_partials = make_partials_propagator(y_ref, nu_ref, s_ref); scalar_seq_view y_vec(y_ref); - scalar_seq_view nu_vec(nu_ref); - scalar_seq_view s_vec(s_ref); - size_t N = max_size(y, nu, s); + scalar_seq_view nu_vec(nu_ref); + scalar_seq_view s_vec(s_ref); + const size_t N = max_size(y, nu, s); - // Explicit return for extreme values - // The gradients are technically ill-defined, but treated as zero - for (size_t i = 0; i < stan::math::size(y); i++) { - if (y_vec.val(i) == 0) { + // An element with probability 0 (z = inf) makes the + // result 0; the gradients are treated as zero + for (size_t n = 0; n < N; n++) { + const double y_rec = value_of_rec(y_vec.val(n)); + const double nu_rec = value_of_rec(nu_vec.val(n)); + const double s_rec = value_of_rec(s_vec.val(n)); + if (0.5 * nu_rec * s_rec * s_rec / y_rec == INFTY) { return ops_partials.build(0.0); } } - VectorBuilder, T_partials_return, T_dof> gamma_vec( - math::size(nu)); - VectorBuilder, T_partials_return, T_dof> digamma_vec( - math::size(nu)); - - if constexpr (is_autodiff_v) { - for (size_t i = 0; i < stan::math::size(nu); i++) { - const T_partials_return half_nu_dbl = 0.5 * nu_vec.val(i); - gamma_vec[i] = tgamma(half_nu_dbl); - digamma_vec[i] = digamma(half_nu_dbl); - } - } - + constexpr bool need_d_a = is_autodiff_v; + constexpr bool need_d_z = is_any_autodiff_v; for (size_t n = 0; n < N; n++) { - // Explicit results for extreme values - // The gradients are technically ill-defined, but treated as zero - if (y_vec.val(n) == INFTY) { + const T_partials_return y_val = y_vec.val(n); + const T_partials_return nu_val = nu_vec.val(n); + const T_partials_return s_val = s_vec.val(n); + const T_partials_return z = 0.5 * nu_val * s_val * s_val / y_val; + const double z_rec = value_of_rec(z); + // An element with probability 1 contributes nothing + if (z_rec == 0) { continue; } - - const T_partials_return y_dbl = y_vec.val(n); - const T_partials_return y_inv_dbl = 1.0 / y_dbl; - const T_partials_return half_nu_dbl = 0.5 * nu_vec.val(n); - const T_partials_return s_dbl = s_vec.val(n); - const T_partials_return half_s2_overx_dbl = 0.5 * s_dbl * s_dbl * y_inv_dbl; - const T_partials_return half_nu_s2_overx_dbl - = 2.0 * half_nu_dbl * half_s2_overx_dbl; - - const T_partials_return Pn = gamma_q(half_nu_dbl, half_nu_s2_overx_dbl); - const T_partials_return gamma_p_deriv - = exp(-half_nu_s2_overx_dbl) - * pow(half_nu_s2_overx_dbl, half_nu_dbl - 1) / tgamma(half_nu_dbl); - - P *= Pn; - + const auto tail + = internal::log_gamma_tail(0.5 * nu_val, z); + log_P += tail.log_prob; if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[n] - += half_nu_s2_overx_dbl * y_inv_dbl * gamma_p_deriv / Pn; + partials<0>(ops_partials)[n] += -z / y_val * tail.d_z; } - if constexpr (is_autodiff_v) { partials<1>(ops_partials)[n] - += (0.5 - * grad_reg_inc_gamma(half_nu_dbl, half_nu_s2_overx_dbl, - gamma_vec[n], digamma_vec[n]) - - half_s2_overx_dbl * gamma_p_deriv) - / Pn; + += 0.5 * tail.d_a + 0.5 * s_val * s_val / y_val * tail.d_z; } - if constexpr (is_autodiff_v) { - partials<2>(ops_partials)[n] - += -2.0 * half_nu_dbl * s_dbl * y_inv_dbl * gamma_p_deriv / Pn; + partials<2>(ops_partials)[n] += nu_val * s_val / y_val * tail.d_z; } } - + const T_partials_return P = exp(log_P); if constexpr (is_autodiff_v) { for (size_t n = 0; n < stan::math::size(y); ++n) { partials<0>(ops_partials)[n] *= P; diff --git a/stan/math/prim/prob/scaled_inv_chi_square_lccdf.hpp b/stan/math/prim/prob/scaled_inv_chi_square_lccdf.hpp index 5d14e326e18..20fba1ae2aa 100644 --- a/stan/math/prim/prob/scaled_inv_chi_square_lccdf.hpp +++ b/stan/math/prim/prob/scaled_inv_chi_square_lccdf.hpp @@ -4,40 +4,51 @@ #include #include #include -#include #include -#include -#include -#include +#include #include #include #include #include -#include #include +#include #include #include namespace stan { namespace math { +/** \ingroup prob_dists + * The log complementary cumulative distribution function of the + * scaled inverse chi-square distribution. + * It is computed from the log of the regularized lower + * incomplete gamma function, so it stays finite where the probability + * underflows. + * + * @tparam T_y type of y + * @tparam T_dof type of nu + * @tparam T_scale type of s + * @param y random variable + * @param nu degrees of freedom parameter + * @param s scale parameter + * @return log probability, or sum of log probabilities for containers + * @throw std::domain_error if a parameter is out of its support + * @throw std::invalid_argument if container sizes mismatch + */ template inline return_type_t scaled_inv_chi_square_lccdf( const T_y& y, const T_dof& nu, const T_scale& s) { using T_partials_return = partials_return_t; - using std::exp; - using std::log; - using std::pow; using T_y_ref = ref_type_t; - using T_nu_ref = ref_type_t; - using T_s_ref = ref_type_t; + using T_dof_ref = ref_type_t; + using T_scale_ref = ref_type_t; static constexpr const char* function = "scaled_inv_chi_square_lccdf"; check_consistent_sizes(function, "Random variable", y, "Degrees of freedom parameter", nu, "Scale parameter", s); T_y_ref y_ref = y; - T_nu_ref nu_ref = nu; - T_s_ref s_ref = s; + T_dof_ref nu_ref = nu; + T_scale_ref s_ref = s; check_nonnegative(function, "Random variable", y_ref); check_positive_finite(function, "Degrees of freedom parameter", nu_ref); check_positive_finite(function, "Scale parameter", s_ref); @@ -46,75 +57,52 @@ inline return_type_t scaled_inv_chi_square_lccdf( return 0; } - T_partials_return P(0.0); + T_partials_return log_P(0.0); auto ops_partials = make_partials_propagator(y_ref, nu_ref, s_ref); scalar_seq_view y_vec(y_ref); - scalar_seq_view nu_vec(nu_ref); - scalar_seq_view s_vec(s_ref); - size_t N = max_size(y, nu, s); + scalar_seq_view nu_vec(nu_ref); + scalar_seq_view s_vec(s_ref); + const size_t N = max_size(y, nu, s); - // Explicit return for extreme values - // The gradients are technically ill-defined, but treated as zero - for (size_t i = 0; i < stan::math::size(y); i++) { - if (y_vec.val(i) == 0) { - return ops_partials.build(0.0); - } - } - - VectorBuilder, T_partials_return, T_dof> gamma_vec( - math::size(nu)); - VectorBuilder, T_partials_return, T_dof> digamma_vec( - math::size(nu)); - - if constexpr (is_autodiff_v) { - for (size_t i = 0; i < stan::math::size(nu); i++) { - const T_partials_return half_nu_dbl = 0.5 * nu_vec.val(i); - gamma_vec[i] = tgamma(half_nu_dbl); - digamma_vec[i] = digamma(half_nu_dbl); + // An element with probability 0 (z = 0) makes the + // result -inf; the gradients are treated as zero + for (size_t n = 0; n < N; n++) { + const double y_rec = value_of_rec(y_vec.val(n)); + const double nu_rec = value_of_rec(nu_vec.val(n)); + const double s_rec = value_of_rec(s_vec.val(n)); + if (0.5 * nu_rec * s_rec * s_rec / y_rec == 0) { + return ops_partials.build(negative_infinity()); } } + constexpr bool need_d_a = is_autodiff_v; + constexpr bool need_d_z = is_any_autodiff_v; for (size_t n = 0; n < N; n++) { - // Explicit results for extreme values - // The gradients are technically ill-defined, but treated as zero - if (y_vec.val(n) == INFTY) { - return ops_partials.build(negative_infinity()); + const T_partials_return y_val = y_vec.val(n); + const T_partials_return nu_val = nu_vec.val(n); + const T_partials_return s_val = s_vec.val(n); + const T_partials_return z = 0.5 * nu_val * s_val * s_val / y_val; + const double z_rec = value_of_rec(z); + // An element with probability 1 contributes nothing + if (z_rec == INFTY) { + continue; } - - const T_partials_return y_dbl = y_vec.val(n); - const T_partials_return y_inv_dbl = 1.0 / y_dbl; - const T_partials_return half_nu_dbl = 0.5 * nu_vec.val(n); - const T_partials_return s_dbl = s_vec.val(n); - const T_partials_return half_s2_overx_dbl = 0.5 * s_dbl * s_dbl * y_inv_dbl; - const T_partials_return half_nu_s2_overx_dbl - = 2.0 * half_nu_dbl * half_s2_overx_dbl; - - const T_partials_return Pn = gamma_p(half_nu_dbl, half_nu_s2_overx_dbl); - const T_partials_return gamma_p_deriv - = exp(-half_nu_s2_overx_dbl) - * pow(half_nu_s2_overx_dbl, half_nu_dbl - 1) / tgamma(half_nu_dbl); - - P += log(Pn); - + const auto tail + = internal::log_gamma_tail(0.5 * nu_val, z); + log_P += tail.log_prob; if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[n] - -= half_nu_s2_overx_dbl * y_inv_dbl * gamma_p_deriv / Pn; + partials<0>(ops_partials)[n] += -z / y_val * tail.d_z; } if constexpr (is_autodiff_v) { partials<1>(ops_partials)[n] - -= (0.5 - * grad_reg_inc_gamma(half_nu_dbl, half_nu_s2_overx_dbl, - gamma_vec[n], digamma_vec[n]) - - half_s2_overx_dbl * gamma_p_deriv) - / Pn; + += 0.5 * tail.d_a + 0.5 * s_val * s_val / y_val * tail.d_z; } if constexpr (is_autodiff_v) { - partials<2>(ops_partials)[n] - += 2.0 * half_nu_dbl * s_dbl * y_inv_dbl * gamma_p_deriv / Pn; + partials<2>(ops_partials)[n] += nu_val * s_val / y_val * tail.d_z; } } - return ops_partials.build(P); + return ops_partials.build(log_P); } } // namespace math diff --git a/stan/math/prim/prob/scaled_inv_chi_square_lcdf.hpp b/stan/math/prim/prob/scaled_inv_chi_square_lcdf.hpp index ad017867c8d..5b364a874b5 100644 --- a/stan/math/prim/prob/scaled_inv_chi_square_lcdf.hpp +++ b/stan/math/prim/prob/scaled_inv_chi_square_lcdf.hpp @@ -4,40 +4,51 @@ #include #include #include -#include #include -#include -#include -#include +#include #include #include #include #include -#include #include +#include #include #include namespace stan { namespace math { +/** \ingroup prob_dists + * The log cumulative distribution function of the scaled inverse + * chi-square distribution. + * It is computed from the log of the regularized upper + * incomplete gamma function, so it stays finite where the probability + * underflows. + * + * @tparam T_y type of y + * @tparam T_dof type of nu + * @tparam T_scale type of s + * @param y random variable + * @param nu degrees of freedom parameter + * @param s scale parameter + * @return log probability, or sum of log probabilities for containers + * @throw std::domain_error if a parameter is out of its support + * @throw std::invalid_argument if container sizes mismatch + */ template inline return_type_t scaled_inv_chi_square_lcdf( const T_y& y, const T_dof& nu, const T_scale& s) { using T_partials_return = partials_return_t; - using std::exp; - using std::log; - using std::pow; using T_y_ref = ref_type_t; - using T_nu_ref = ref_type_t; - using T_s_ref = ref_type_t; + using T_dof_ref = ref_type_t; + using T_scale_ref = ref_type_t; static constexpr const char* function = "scaled_inv_chi_square_lcdf"; check_consistent_sizes(function, "Random variable", y, "Degrees of freedom parameter", nu, "Scale parameter", s); T_y_ref y_ref = y; - T_nu_ref nu_ref = nu; - T_s_ref s_ref = s; + T_dof_ref nu_ref = nu; + T_scale_ref s_ref = s; check_nonnegative(function, "Random variable", y_ref); check_positive_finite(function, "Degrees of freedom parameter", nu_ref); check_positive_finite(function, "Scale parameter", s_ref); @@ -46,75 +57,52 @@ inline return_type_t scaled_inv_chi_square_lcdf( return 0; } - T_partials_return P(0.0); + T_partials_return log_P(0.0); auto ops_partials = make_partials_propagator(y_ref, nu_ref, s_ref); scalar_seq_view y_vec(y_ref); - scalar_seq_view nu_vec(nu_ref); - scalar_seq_view s_vec(s_ref); - size_t N = max_size(y, nu, s); + scalar_seq_view nu_vec(nu_ref); + scalar_seq_view s_vec(s_ref); + const size_t N = max_size(y, nu, s); - // Explicit return for extreme values - // The gradients are technically ill-defined, but treated as zero - for (size_t i = 0; i < stan::math::size(y); i++) { - if (y_vec.val(i) == 0) { + // An element with probability 0 (z = inf) makes the + // result -inf; the gradients are treated as zero + for (size_t n = 0; n < N; n++) { + const double y_rec = value_of_rec(y_vec.val(n)); + const double nu_rec = value_of_rec(nu_vec.val(n)); + const double s_rec = value_of_rec(s_vec.val(n)); + if (0.5 * nu_rec * s_rec * s_rec / y_rec == INFTY) { return ops_partials.build(negative_infinity()); } } - VectorBuilder, T_partials_return, T_dof> gamma_vec( - math::size(nu)); - VectorBuilder, T_partials_return, T_dof> digamma_vec( - math::size(nu)); - - if constexpr (is_autodiff_v) { - for (size_t i = 0; i < stan::math::size(nu); i++) { - const T_partials_return half_nu_dbl = 0.5 * nu_vec.val(i); - gamma_vec[i] = tgamma(half_nu_dbl); - digamma_vec[i] = digamma(half_nu_dbl); - } - } - + constexpr bool need_d_a = is_autodiff_v; + constexpr bool need_d_z = is_any_autodiff_v; for (size_t n = 0; n < N; n++) { - // Explicit results for extreme values - // The gradients are technically ill-defined, but treated as zero - if (y_vec.val(n) == INFTY) { + const T_partials_return y_val = y_vec.val(n); + const T_partials_return nu_val = nu_vec.val(n); + const T_partials_return s_val = s_vec.val(n); + const T_partials_return z = 0.5 * nu_val * s_val * s_val / y_val; + const double z_rec = value_of_rec(z); + // An element with probability 1 contributes nothing + if (z_rec == 0) { continue; } - - const T_partials_return y_dbl = y_vec.val(n); - const T_partials_return y_inv_dbl = 1.0 / y_dbl; - const T_partials_return half_nu_dbl = 0.5 * nu_vec.val(n); - const T_partials_return s_dbl = s_vec.val(n); - const T_partials_return half_s2_overx_dbl = 0.5 * s_dbl * s_dbl * y_inv_dbl; - const T_partials_return half_nu_s2_overx_dbl - = 2.0 * half_nu_dbl * half_s2_overx_dbl; - - const T_partials_return Pn = gamma_q(half_nu_dbl, half_nu_s2_overx_dbl); - const T_partials_return gamma_p_deriv - = exp(-half_nu_s2_overx_dbl) - * pow(half_nu_s2_overx_dbl, half_nu_dbl - 1) / tgamma(half_nu_dbl); - - P += log(Pn); - + const auto tail + = internal::log_gamma_tail(0.5 * nu_val, z); + log_P += tail.log_prob; if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[n] - += half_nu_s2_overx_dbl * y_inv_dbl * gamma_p_deriv / Pn; + partials<0>(ops_partials)[n] += -z / y_val * tail.d_z; } if constexpr (is_autodiff_v) { partials<1>(ops_partials)[n] - += (0.5 - * grad_reg_inc_gamma(half_nu_dbl, half_nu_s2_overx_dbl, - gamma_vec[n], digamma_vec[n]) - - half_s2_overx_dbl * gamma_p_deriv) - / Pn; + += 0.5 * tail.d_a + 0.5 * s_val * s_val / y_val * tail.d_z; } if constexpr (is_autodiff_v) { - partials<2>(ops_partials)[n] - += -2.0 * half_nu_dbl * s_dbl * y_inv_dbl * gamma_p_deriv / Pn; + partials<2>(ops_partials)[n] += nu_val * s_val / y_val * tail.d_z; } } - return ops_partials.build(P); + return ops_partials.build(log_P); } } // namespace math diff --git a/test/unit/math/mix/prob/incomplete_gamma_cdfs_test.cpp b/test/unit/math/mix/prob/incomplete_gamma_cdfs_test.cpp new file mode 100644 index 00000000000..9b463d11b91 --- /dev/null +++ b/test/unit/math/mix/prob/incomplete_gamma_cdfs_test.cpp @@ -0,0 +1,58 @@ +#include +#include + +// Each function at two points: z > a + 1 (continued fraction) and +// z < a + 1 (power series side) of the incomplete gamma function. + +TEST_F(AgradRev, mathMixScalFun_gamma_cdfs) { + auto check = [](const auto& f) { + stan::test::expect_ad(f, 2.5, 3.0, 1.5); + stan::test::expect_ad(f, 1.0, 4.0, 1.5); + }; + check([](const auto& y, const auto& a, const auto& b) { + return stan::math::gamma_cdf(y, a, b); + }); + check([](const auto& y, const auto& a, const auto& b) { + return stan::math::gamma_lcdf(y, a, b); + }); + check([](const auto& y, const auto& a, const auto& b) { + return stan::math::gamma_lccdf(y, a, b); + }); +} + +TEST_F(AgradRev, mathMixScalFun_chi_square_cdfs) { + auto check = [](const auto& f) { + stan::test::expect_ad(f, 5.0, 3.0); + stan::test::expect_ad(f, 2.0, 6.0); + }; + check([](const auto& y, const auto& nu) { + return stan::math::chi_square_cdf(y, nu); + }); + check([](const auto& y, const auto& nu) { + return stan::math::chi_square_lcdf(y, nu); + }); + check([](const auto& y, const auto& nu) { + return stan::math::chi_square_lccdf(y, nu); + }); +} + +TEST_F(AgradRev, mathMixScalFun_poisson_cdfs) { + // the count is bound into the functor; only lambda is differentiated + auto check = [](const auto& f) { + stan::test::expect_ad(f(3), 6.5); + stan::test::expect_ad(f(5), 2.0); + }; + check([](int n) { + return + [n](const auto& lambda) { return stan::math::poisson_cdf(n, lambda); }; + }); + check([](int n) { + return + [n](const auto& lambda) { return stan::math::poisson_lcdf(n, lambda); }; + }); + check([](int n) { + return [n](const auto& lambda) { + return stan::math::poisson_lccdf(n, lambda); + }; + }); +} diff --git a/test/unit/math/mix/prob/inv_incomplete_gamma_cdfs_test.cpp b/test/unit/math/mix/prob/inv_incomplete_gamma_cdfs_test.cpp new file mode 100644 index 00000000000..7e90d485ee2 --- /dev/null +++ b/test/unit/math/mix/prob/inv_incomplete_gamma_cdfs_test.cpp @@ -0,0 +1,53 @@ +#include +#include + +// Each function at two points: z > a + 1 (continued fraction) and +// z < a + 1 (power series side) of the incomplete gamma function. + +TEST_F(AgradRev, mathMixScalFun_inv_gamma_cdfs) { + auto check = [](const auto& f) { + stan::test::expect_ad(f, 0.4, 2.0, 1.5); + stan::test::expect_ad(f, 2.0, 3.0, 1.0); + }; + check([](const auto& y, const auto& a, const auto& b) { + return stan::math::inv_gamma_cdf(y, a, b); + }); + check([](const auto& y, const auto& a, const auto& b) { + return stan::math::inv_gamma_lcdf(y, a, b); + }); + check([](const auto& y, const auto& a, const auto& b) { + return stan::math::inv_gamma_lccdf(y, a, b); + }); +} + +TEST_F(AgradRev, mathMixScalFun_inv_chi_square_cdfs) { + auto check = [](const auto& f) { + stan::test::expect_ad(f, 0.1, 3.0); + stan::test::expect_ad(f, 1.0, 4.0); + }; + check([](const auto& y, const auto& nu) { + return stan::math::inv_chi_square_cdf(y, nu); + }); + check([](const auto& y, const auto& nu) { + return stan::math::inv_chi_square_lcdf(y, nu); + }); + check([](const auto& y, const auto& nu) { + return stan::math::inv_chi_square_lccdf(y, nu); + }); +} + +TEST_F(AgradRev, mathMixScalFun_scaled_inv_chi_square_cdfs) { + auto check = [](const auto& f) { + stan::test::expect_ad(f, 0.5, 3.0, 1.0); + stan::test::expect_ad(f, 2.0, 4.0, 0.8); + }; + check([](const auto& y, const auto& nu, const auto& s) { + return stan::math::scaled_inv_chi_square_cdf(y, nu, s); + }); + check([](const auto& y, const auto& nu, const auto& s) { + return stan::math::scaled_inv_chi_square_lcdf(y, nu, s); + }); + check([](const auto& y, const auto& nu, const auto& s) { + return stan::math::scaled_inv_chi_square_lccdf(y, nu, s); + }); +} diff --git a/test/unit/math/prim/prob/gamma_lccdf_test.cpp b/test/unit/math/prim/prob/gamma_lccdf_test.cpp index e0c84861e0c..9a7e3d3082c 100644 --- a/test/unit/math/prim/prob/gamma_lccdf_test.cpp +++ b/test/unit/math/prim/prob/gamma_lccdf_test.cpp @@ -91,24 +91,24 @@ TEST(ProbGamma, lccdf_alpha_gt_30_small_y_old_code_rounds_to_zero) { EXPECT_DOUBLE_EQ(new_val, expected); } -TEST(ProbGamma, lccdf_log1m_exp_lcdf_rounds_to_inf) { +TEST(ProbGamma, lccdf_far_tail_matches_log1m_exp_lcdf) { using stan::math::gamma_lccdf; using stan::math::gamma_lcdf; using stan::math::log1m_exp; - using stan::math::negative_infinity; double y = 20000.0; double alpha = 800.0; double beta = 0.1; + // lcdf is about -1e-204, so it does not round to 0, and + // log1m_exp(lcdf) recovers lccdf double lcdf = gamma_lcdf(y, alpha, beta); - double log1m_lcdf = log1m_exp(lcdf); double lccdf = gamma_lccdf(y, alpha, beta); - EXPECT_EQ(lcdf, 0.0); - EXPECT_EQ(log1m_lcdf, negative_infinity()); + EXPECT_LT(lcdf, 0.0); EXPECT_TRUE(std::isfinite(lccdf)); EXPECT_LT(lccdf, 0.0); + EXPECT_NEAR(log1m_exp(lcdf), lccdf, 1e-12 * std::fabs(lccdf)); } TEST(ProbGamma, lccdf_large_alpha_large_y) { diff --git a/test/unit/math/rev/prob/gamma_lccdf_test.cpp b/test/unit/math/rev/prob/gamma_lccdf_test.cpp index 778d05d52f7..d36a57e4290 100644 --- a/test/unit/math/rev/prob/gamma_lccdf_test.cpp +++ b/test/unit/math/rev/prob/gamma_lccdf_test.cpp @@ -275,27 +275,28 @@ TEST_F( EXPECT_LE(grads[0], 0.0); } -TEST_F(AgradRev, ProbDistributionsGammalccdf_log1m_exp_lcdf_rounds_to_inf) { +TEST_F(AgradRev, ProbDistributionsGammalccdf_far_tail_matches_log1m_exp_lcdf) { using stan::math::gamma_lccdf; using stan::math::gamma_lcdf; using stan::math::log1m_exp; - using stan::math::negative_infinity; using stan::math::var; double y_d = 20000.0; double alpha_d = 800.0; double beta_d = 0.1; + // lcdf is about -1e-204, so it does not round to 0, and + // log1m_exp(lcdf) recovers lccdf double lcdf = gamma_lcdf(y_d, alpha_d, beta_d); - double log1m_lcdf = log1m_exp(lcdf); var y_v = y_d; var alpha_v = alpha_d; var beta_v = beta_d; var lccdf_var = gamma_lccdf(y_v, alpha_v, beta_v); - EXPECT_EQ(lcdf, 0.0); - EXPECT_EQ(log1m_lcdf, negative_infinity()); + EXPECT_LT(lcdf, 0.0); + EXPECT_NEAR(log1m_exp(lcdf), lccdf_var.val(), + 1e-12 * std::fabs(lccdf_var.val())); EXPECT_TRUE(std::isfinite(lccdf_var.val())); EXPECT_LT(lccdf_var.val(), 0.0); } diff --git a/test/unit/math/rev/prob/incomplete_gamma_cdfs_reference.hpp b/test/unit/math/rev/prob/incomplete_gamma_cdfs_reference.hpp new file mode 100644 index 00000000000..2a33658e286 --- /dev/null +++ b/test/unit/math/rev/prob/incomplete_gamma_cdfs_reference.hpp @@ -0,0 +1,370 @@ +#ifndef TEST_UNIT_MATH_REV_PROB_INCOMPLETE_GAMMA_CDFS_REFERENCE_HPP +#define TEST_UNIT_MATH_REV_PROB_INCOMPLETE_GAMMA_CDFS_REFERENCE_HPP + +#include +#include + +namespace stan { +namespace test { +namespace incomplete_gamma_cdfs { + +/** + * Values and gradients of the cdf, lcdf and lccdf of chi_square, gamma, + * inv_chi_square, inv_gamma, scaled_inv_chi_square and poisson, computed + * with mpmath at 60 digits by two routes: the closed form in the + * regularized incomplete gamma function with numerical differentiation, + * and a quadrature of the scaled density with the chain rule. Only rows + * where the routes agree to 1e-25 are kept. Each row is one point (a, z) + * of the incomplete gamma function, mapped to the arguments of the + * function; the comment names the region. The gradient is with respect to + * (y, p1, p2); for poisson with respect to lambda only. + */ +struct Row { + std::string function; + double y; + double p1; + double p2; + double value; + double grad[3]; + int n_grad; +}; + +inline const std::vector& rows() { + // clang-format off + // NOLINTBEGIN(whitespace/line_length) + static const std::vector r = { + {"chi_square_lcdf", 2000.0, 20.0, 0.0, -1.4114888035653841e-413, {6.9939909853800936e-414, -3.2866592567233371e-413, 0}, 2}, // upper_deep + {"chi_square_lcdf", 2400.0, 300.0, 0.0, -1.3223466032724075e-323, {5.7915555953704621e-324, -1.3777059130352209e-323, 0}, 2}, // upper_deep + {"chi_square_lcdf", 3200.0, 400.0, 0.0, -1.6236213327686189e-430, {7.1091366676842467e-431, -1.6907229624739902e-430, 0}, 2}, // upper_deep + {"chi_square_lcdf", 600.0, 20.0, 0.0, -2.8784992655370182e-114, {1.3962194799621868e-114, -4.9732607083897553e-114, 0}, 2}, // upper_tail + {"chi_square_lcdf", 1200.0, 300.0, 0.0, -8.1414836560590812e-108, {3.0620693166254771e-108, -5.665823750109828e-108, 0}, 2}, // upper_tail + {"chi_square_lcdf", 15.0, 10.0, 0.0, -0.14163482986540183, {0.042005525960158956, -0.05282592190639068, 0}, 2}, // bulk + {"chi_square_lcdf", 120.0, 100.0, 0.0, -0.088182987887198306, {0.010590217406820947, -0.011623606434664016, 0}, 2}, // bulk + {"chi_square_lcdf", 400.0, 400.0, 0.0, -0.67451500659035485, {0.027677206992671042, -0.027700282830990659, 0}, 2}, // bulk_big + {"chi_square_lcdf", 600.0, 400.0, 0.0, -3.3711032560941012e-10, {5.7795826456264662e-11, -7.036087390357163e-11, 0}, 2}, // bulk_big + {"chi_square_lcdf", 800.0, 800.0, 0.0, -0.6799365665689943, {0.019681232662310234, -0.019689435223116438, 0}, 2}, // bulk_big + {"chi_square_lcdf", 2200.0, 2000.0, 0.0, -0.0010598847333686436, {5.2741012602257179e-5, -5.5303523228577231e-5, 0}, 2}, // bulk_big + {"chi_square_lcdf", 1800.0, 2000.0, 0.0, -7.5057699942338924, {0.059878855399343508, -0.056789264431900816, 0}, 2}, // bulk_big + {"chi_square_lcdf", 600.0, 800.0, 0.0, -17.629204302744157, {0.17136338373499565, -0.14795598637042751, 0}, 2}, // bulk_big + {"chi_square_lcdf", 100.0, 400.0, 0.0, -130.54189961219741, {1.5032969848784755, -0.69521496152576584, 0}, 2}, // lower_tail + {"chi_square_lcdf", 10.0, 100.0, 0.0, -72.902912338991312, {4.5108198474797784, -1.1573345022795921, 0}, 2}, // lower_tail + {"chi_square_lcdf", 200.0, 2000.0, 0.0, -1406.8527552662533, {4.5005548024741409, -1.1515979235241671, 0}, 2}, // lower_deep + {"chi_square_lcdf", 60.0, 1000.0, 0.0, -940.67002769935046, {7.8343946181427776, -1.4072687336610712, 0}, 2}, // lower_deep + {"chi_square_lcdf", 0.002, 200.0, 0.0, -1054.515893552739, {49999.504950543103, -5.7589586148753858, 0}, 2}, // lower_deep + {"chi_square_lccdf", 2000.0, 20.0, 0.0, -950.62299837015608, {-0.49550453164867151, 2.328505368530959, 0}, 2}, // upper_deep + {"chi_square_lccdf", 2400.0, 300.0, 0.0, -743.45557714901795, {-0.43797560949891014, 1.0418644473588209, 0}, 2}, // upper_deep + {"chi_square_lccdf", 3200.0, 400.0, 0.0, -989.62693094236584, {-0.4378568157614473, 1.0413283740186813, 0}, 2}, // upper_deep + {"chi_square_lccdf", 600.0, 20.0, 0.0, -261.43743153134228, {-0.48505118506664113, 1.727726933243436, 0}, 2}, // upper_tail + {"chi_square_lccdf", 1200.0, 300.0, 0.0, -246.58221761262203, {-0.37610703969744065, 0.69592030021373197, 0}, 2}, // upper_tail + {"chi_square_lccdf", 15.0, 10.0, 0.0, -2.0244848579180109, {-0.27606910316395642, 0.34718300869125189, 0}, 2}, // bulk + {"chi_square_lccdf", 120.0, 100.0, 0.0, -2.4721087206353381, {-0.11487636023370675, 0.12608594788084989, 0}, 2}, // bulk + {"chi_square_lccdf", 400.0, 400.0, 0.0, -0.71213311417597682, {-0.028738202344769638, 0.028762162786698621, 0}, 2}, // bulk_big + {"chi_square_lccdf", 600.0, 400.0, 0.0, -21.810610863624893, {-0.17144484180970326, 0.20871764685497929, 0}, 2}, // bulk_big + {"chi_square_lccdf", 800.0, 800.0, 0.0, -0.70653465390004754, {-0.020211739769021841, 0.02022016342963798, 0}, 2}, // bulk_big + {"chi_square_lccdf", 2200.0, 2000.0, 0.0, -6.8501250144362653, {-0.049734717543219674, 0.052151162277844894, 0}, 2}, // bulk_big + {"chi_square_lccdf", 1800.0, 2000.0, 0.0, -0.00055005351741434473, {-3.2945635131731474e-5, 3.1245727275430397e-5, 0}, 2}, // bulk_big + {"chi_square_lccdf", 600.0, 800.0, 0.0, -2.2066520386350113e-8, {-3.7813936423833186e-9, 3.2648738255477803e-9, 0}, 2}, // bulk_big + {"chi_square_lccdf", 100.0, 400.0, 0.0, -2.0247590148473564e-57, {-3.0438141221255432e-57, 1.4076427606060524e-57, 0}, 2}, // lower_tail + {"chi_square_lccdf", 10.0, 100.0, 0.0, -2.1810592140784888e-32, {-9.8383651913938939e-32, 2.5242150799678461e-32, 0}, 2}, // lower_tail + {"chi_square_lccdf", 200.0, 2000.0, 0.0, -1.0270971815582278e-611, {-4.6225071530695367e-611, 1.1828029815399796e-611, 0}, 2}, // lower_deep + {"chi_square_lccdf", 60.0, 1000.0, 0.0, -2.9661812029194365e-409, {-2.3238234052588304e-408, 4.1742140652417082e-409, 0}, 2}, // lower_deep + {"chi_square_lccdf", 0.002, 200.0, 0.0, -1.070449911926898e-458, {-5.3521965670697366e-454, 6.1646767420840073e-458, 0}, 2}, // lower_deep + {"chi_square_cdf", 2000.0, 20.0, 0.0, 1.0, {6.9939909853800936e-414, -3.2866592567233371e-413, 0}, 2}, // upper_deep + {"chi_square_cdf", 2400.0, 300.0, 0.0, 1.0, {5.7915555953704621e-324, -1.3777059130352209e-323, 0}, 2}, // upper_deep + {"chi_square_cdf", 3200.0, 400.0, 0.0, 1.0, {7.1091366676842467e-431, -1.6907229624739902e-430, 0}, 2}, // upper_deep + {"chi_square_cdf", 600.0, 20.0, 0.0, 1.0, {1.3962194799621868e-114, -4.9732607083897553e-114, 0}, 2}, // upper_tail + {"chi_square_cdf", 1200.0, 300.0, 0.0, 1.0, {3.0620693166254771e-108, -5.665823750109828e-108, 0}, 2}, // upper_tail + {"chi_square_cdf", 15.0, 10.0, 0.0, 0.86793814371227939, {0.036458198227518327, -0.045849632599322562, 0}, 2}, // bulk + {"chi_square_cdf", 120.0, 100.0, 0.0, 0.91559331890630817, {0.0096963323034505477, -0.010642496393174746, 0}, 2}, // bulk + {"chi_square_cdf", 400.0, 400.0, 0.0, 0.50940341800723633, {0.014098863842960411, -0.014110618753873806, 0}, 2}, // bulk_big + {"chi_square_cdf", 600.0, 400.0, 0.0, 0.99999999966288967, {5.7795826436781092e-11, -7.0360873879852253e-11, 0}, 2}, // bulk_big + {"chi_square_cdf", 800.0, 800.0, 0.0, 0.50664912983890547, {0.0099714794025165248, -0.0099756352228114391, 0}, 2}, // bulk_big + {"chi_square_cdf", 2200.0, 2000.0, 0.0, 0.99894067674607002, {5.2685142821171794e-5, -5.5244938920396943e-5, 0}, 2}, // bulk_big + {"chi_square_cdf", 1800.0, 2000.0, 0.0, 0.00054990226571178292, {3.2927518252327221e-5, -3.1228545179207826e-5, 0}, 2}, // bulk_big + {"chi_square_cdf", 600.0, 800.0, 0.0, 2.2066520142884454e-8, {3.7813935589411196e-9, -3.2648737535033762e-9, 0}, 2}, // bulk_big + {"chi_square_cdf", 100.0, 400.0, 0.0, 2.0247590148473564e-57, {3.0438141221255432e-57, -1.4076427606060524e-57, 0}, 2}, // lower_tail + {"chi_square_cdf", 10.0, 100.0, 0.0, 2.1810592140784888e-32, {9.8383651913938939e-32, -2.5242150799678461e-32, 0}, 2}, // lower_tail + {"chi_square_cdf", 200.0, 2000.0, 0.0, 1.0270971815582278e-611, {4.6225071530695367e-611, -1.1828029815399796e-611, 0}, 2}, // lower_deep + {"chi_square_cdf", 60.0, 1000.0, 0.0, 2.9661812029194365e-409, {2.3238234052588304e-408, -4.1742140652417082e-409, 0}, 2}, // lower_deep + {"chi_square_cdf", 0.002, 200.0, 0.0, 1.070449911926898e-458, {5.3521965670697366e-454, -6.1646767420840073e-458, 0}, 2}, // lower_deep + {"gamma_lcdf", 500.0, 10.0, 2.0, -1.4114888035653841e-413, {2.7975963941520375e-413, -6.5733185134466742e-413, 6.9939909853800936e-411}, 3}, // upper_deep + {"gamma_lcdf", 600.0, 150.0, 2.0, -1.3223466032724075e-323, {2.3166222381481848e-323, -2.7554118260704417e-323, 6.9498667144445545e-321}, 3}, // upper_deep + {"gamma_lcdf", 800.0, 200.0, 2.0, -1.6236213327686189e-430, {2.8436546670736987e-430, -3.3814459249479805e-430, 1.1374618668294795e-427}, 3}, // upper_deep + {"gamma_lcdf", 150.0, 10.0, 2.0, -2.8784992655370182e-114, {5.5848779198487471e-114, -9.9465214167795107e-114, 4.1886584398865603e-112}, 3}, // upper_tail + {"gamma_lcdf", 300.0, 150.0, 2.0, -8.1414836560590812e-108, {1.2248277266501908e-107, -1.1331647500219656e-107, 1.8372415899752863e-105}, 3}, // upper_tail + {"gamma_lcdf", 3.75, 5.0, 2.0, -0.14163482986540183, {0.16802210384063583, -0.10565184381278136, 0.31504144470119217}, 3}, // bulk + {"gamma_lcdf", 30.0, 50.0, 2.0, -0.088182987887198306, {0.04236086962728379, -0.023247212869328032, 0.63541304440925685}, 3}, // bulk + {"gamma_lcdf", 100.0, 200.0, 2.0, -0.67451500659035485, {0.11070882797068417, -0.055400565661981318, 5.5354413985342084}, 3}, // bulk_big + {"gamma_lcdf", 150.0, 200.0, 2.0, -3.3711032560941012e-10, {2.3118330582505865e-10, -1.4072174780714326e-10, 1.7338747936879398e-8}, 3}, // bulk_big + {"gamma_lcdf", 200.0, 400.0, 2.0, -0.6799365665689943, {0.078724930649240935, -0.039378870446232877, 7.8724930649240935}, 3}, // bulk_big + {"gamma_lcdf", 550.0, 1000.0, 2.0, -0.0010598847333686436, {0.00021096405040902872, -0.00011060704645715446, 0.058015113862482897}, 3}, // bulk_big + {"gamma_lcdf", 450.0, 1000.0, 2.0, -7.5057699942338924, {0.23951542159737403, -0.11357852886380163, 53.890969859409157}, 3}, // bulk_big + {"gamma_lcdf", 150.0, 400.0, 2.0, -17.629204302744157, {0.68545353493998259, -0.29591197274085501, 51.409015120498695}, 3}, // bulk_big + {"gamma_lcdf", 25.0, 200.0, 2.0, -130.54189961219741, {6.0131879395139018, -1.3904299230515317, 75.164849243923773}, 3}, // lower_tail + {"gamma_lcdf", 2.5, 50.0, 2.0, -72.902912338991312, {18.043279389919113, -2.3146690045591842, 22.554099237398892}, 3}, // lower_tail + {"gamma_lcdf", 50.0, 1000.0, 2.0, -1406.8527552662533, {18.002219209896564, -2.3031958470483342, 450.05548024741409}, 3}, // lower_deep + {"gamma_lcdf", 15.0, 500.0, 2.0, -940.67002769935046, {31.337578472571111, -2.8145374673221423, 235.03183854428333}, 3}, // lower_deep + {"gamma_lcdf", 0.0005, 100.0, 2.0, -1054.515893552739, {199998.01980217241, -11.517917229750772, 49.999504950543104}, 3}, // lower_deep + {"gamma_lccdf", 500.0, 10.0, 2.0, -950.62299837015608, {-1.982018126594686, 4.6570107370619179, -495.50453164867151}, 3}, // upper_deep + {"gamma_lccdf", 600.0, 150.0, 2.0, -743.45557714901795, {-1.7519024379956406, 2.0837288947176418, -525.57073139869217}, 3}, // upper_deep + {"gamma_lccdf", 800.0, 200.0, 2.0, -989.62693094236584, {-1.7514272630457892, 2.0826567480373627, -700.57090521831568}, 3}, // upper_deep + {"gamma_lccdf", 150.0, 10.0, 2.0, -261.43743153134228, {-1.9402047402665645, 3.455453866486872, -145.51535551999234}, 3}, // upper_tail + {"gamma_lccdf", 300.0, 150.0, 2.0, -246.58221761262203, {-1.5044281587897626, 1.3918406004274639, -225.66422381846439}, 3}, // upper_tail + {"gamma_lccdf", 3.75, 5.0, 2.0, -2.0244848579180109, {-1.1042764126558257, 0.69436601738250379, -2.0705182737296731}, 3}, // bulk + {"gamma_lccdf", 30.0, 50.0, 2.0, -2.4721087206353381, {-0.45950544093482701, 0.25217189576169979, -6.8925816140224052}, 3}, // bulk + {"gamma_lccdf", 100.0, 200.0, 2.0, -0.71213311417597682, {-0.11495280937907855, 0.057524325573397242, -5.7476404689539276}, 3}, // bulk_big + {"gamma_lccdf", 150.0, 200.0, 2.0, -21.810610863624893, {-0.68577936723881305, 0.41743529370995857, -51.433452542910979}, 3}, // bulk_big + {"gamma_lccdf", 200.0, 400.0, 2.0, -0.70653465390004754, {-0.080846959076087362, 0.040440326859275961, -8.0846959076087362}, 3}, // bulk_big + {"gamma_lccdf", 550.0, 1000.0, 2.0, -6.8501250144362653, {-0.1989388701728787, 0.10430232455568979, -54.708189297541642}, 3}, // bulk_big + {"gamma_lccdf", 450.0, 1000.0, 2.0, -0.00055005351741434473, {-0.0001317825405269259, 6.2491454550860795e-5, -0.029651071618558327}, 3}, // bulk_big + {"gamma_lccdf", 150.0, 400.0, 2.0, -2.2066520386350113e-8, {-1.5125574569533275e-8, 6.5297476510955605e-9, -1.1344180927149956e-6}, 3}, // bulk_big + {"gamma_lccdf", 25.0, 200.0, 2.0, -2.0247590148473564e-57, {-1.2175256488502173e-56, 2.8152855212121049e-57, -1.5219070610627716e-55}, 3}, // lower_tail + {"gamma_lccdf", 2.5, 50.0, 2.0, -2.1810592140784888e-32, {-3.9353460765575576e-31, 5.0484301599356922e-32, -4.919182595696947e-31}, 3}, // lower_tail + {"gamma_lccdf", 50.0, 1000.0, 2.0, -1.0270971815582278e-611, {-1.8490028612278147e-610, 2.3656059630799592e-611, -4.6225071530695367e-609}, 3}, // lower_deep + {"gamma_lccdf", 15.0, 500.0, 2.0, -2.9661812029194365e-409, {-9.2952936210353214e-408, 8.3484281304834163e-409, -6.9714702157764911e-407}, 3}, // lower_deep + {"gamma_lccdf", 0.0005, 100.0, 2.0, -1.070449911926898e-458, {-2.1408786268278946e-453, 1.2329353484168015e-457, -5.3521965670697367e-457}, 3}, // lower_deep + {"gamma_cdf", 500.0, 10.0, 2.0, 1.0, {2.7975963941520375e-413, -6.5733185134466742e-413, 6.9939909853800936e-411}, 3}, // upper_deep + {"gamma_cdf", 600.0, 150.0, 2.0, 1.0, {2.3166222381481848e-323, -2.7554118260704417e-323, 6.9498667144445545e-321}, 3}, // upper_deep + {"gamma_cdf", 800.0, 200.0, 2.0, 1.0, {2.8436546670736987e-430, -3.3814459249479805e-430, 1.1374618668294795e-427}, 3}, // upper_deep + {"gamma_cdf", 150.0, 10.0, 2.0, 1.0, {5.5848779198487471e-114, -9.9465214167795107e-114, 4.1886584398865603e-112}, 3}, // upper_tail + {"gamma_cdf", 300.0, 150.0, 2.0, 1.0, {1.2248277266501908e-107, -1.1331647500219656e-107, 1.8372415899752863e-105}, 3}, // upper_tail + {"gamma_cdf", 3.75, 5.0, 2.0, 0.86793814371227939, {0.14583279291007331, -0.091699265198645124, 0.27343648670638745}, 3}, // bulk + {"gamma_cdf", 30.0, 50.0, 2.0, 0.91559331890630817, {0.038785329213802191, -0.021284992786349492, 0.58177993820703286}, 3}, // bulk + {"gamma_cdf", 100.0, 200.0, 2.0, 0.50940341800723633, {0.056395455371841644, -0.028221237507747612, 2.8197727685920822}, 3}, // bulk_big + {"gamma_cdf", 150.0, 200.0, 2.0, 0.99999999966288967, {2.3118330574712437e-10, -1.4072174775970451e-10, 1.7338747931034327e-8}, 3}, // bulk_big + {"gamma_cdf", 200.0, 400.0, 2.0, 0.50664912983890547, {0.039885917610066099, -0.019951270445622878, 3.9885917610066099}, 3}, // bulk_big + {"gamma_cdf", 550.0, 1000.0, 2.0, 0.99894067674607002, {0.00021074057128468718, -0.00011048987784079389, 0.057953657103288974}, 3}, // bulk_big + {"gamma_cdf", 450.0, 1000.0, 2.0, 0.00054990226571178292, {0.00013171007300930889, -6.2457090358415652e-5, 0.029634766427094499}, 3}, // bulk_big + {"gamma_cdf", 150.0, 400.0, 2.0, 2.2066520142884454e-8, {1.5125574235764479e-8, -6.5297475070067525e-9, 1.1344180676823359e-6}, 3}, // bulk_big + {"gamma_cdf", 25.0, 200.0, 2.0, 2.0247590148473564e-57, {1.2175256488502173e-56, -2.8152855212121049e-57, 1.5219070610627716e-55}, 3}, // lower_tail + {"gamma_cdf", 2.5, 50.0, 2.0, 2.1810592140784888e-32, {3.9353460765575576e-31, -5.0484301599356922e-32, 4.919182595696947e-31}, 3}, // lower_tail + {"gamma_cdf", 50.0, 1000.0, 2.0, 1.0270971815582278e-611, {1.8490028612278147e-610, -2.3656059630799592e-611, 4.6225071530695367e-609}, 3}, // lower_deep + {"gamma_cdf", 15.0, 500.0, 2.0, 2.9661812029194365e-409, {9.2952936210353214e-408, -8.3484281304834163e-409, 6.9714702157764911e-407}, 3}, // lower_deep + {"gamma_cdf", 0.0005, 100.0, 2.0, 1.070449911926898e-458, {2.1408786268278946e-453, -1.2329353484168015e-457, 5.3521965670697367e-457}, 3}, // lower_deep + {"inv_gamma_lcdf", 0.002, 10.0, 2.0, -950.62299837015606, {495504.53164867149, 4.6570107370619179, -495.5045316486715}, 3}, // upper_deep + {"inv_gamma_lcdf", 0.0016666666666666668, 150.0, 2.0, -743.45557714901788, {630684.87767843052, 2.0837288947176418, -525.57073139869213}, 3}, // upper_deep + {"inv_gamma_lcdf", 0.00125, 200.0, 2.0, -989.62693094236581, {1120913.448349305, 2.0826567480373627, -700.57090521831566}, 3}, // upper_deep + {"inv_gamma_lcdf", 0.006666666666666667, 10.0, 2.0, -261.43743153134226, {43654.606655997696, 3.4554538664868719, -145.51535551999233}, 3}, // upper_tail + {"inv_gamma_lcdf", 0.0033333333333333335, 150.0, 2.0, -246.582217612622, {135398.53429107861, 1.3918406004274639, -225.66422381846437}, 3}, // upper_tail + {"inv_gamma_lcdf", 0.26666666666666666, 5.0, 2.0, -2.024484857918011, {15.528887052972549, 0.6943660173825038, -2.0705182737296732}, 3}, // bulk + {"inv_gamma_lcdf", 0.03333333333333333, 50.0, 2.0, -2.4721087206353383, {413.55489684134434, 0.2521718957616998, -6.8925816140224055}, 3}, // bulk + {"inv_gamma_lcdf", 0.01, 200.0, 2.0, -0.71213311417597658, {1149.5280937907852, 0.057524325573397229, -5.7476404689539262}, 3}, // bulk_big + {"inv_gamma_lcdf", 0.006666666666666667, 200.0, 2.0, -21.810610863624886, {15430.03576287329, 0.41743529370995851, -51.43345254291097}, 3}, // bulk_big + {"inv_gamma_lcdf", 0.005, 400.0, 2.0, -0.70653465390004721, {3233.8783630434933, 0.040440326859275947, -8.0846959076087335}, 3}, // bulk_big + {"inv_gamma_lcdf", 0.0018181818181818182, 1000.0, 2.0, -6.8501250144362654, {60179.008227295806, 0.10430232455568979, -54.708189297541642}, 3}, // bulk_big + {"inv_gamma_lcdf", 0.0022222222222222222, 1000.0, 2.0, -0.00055005351741434478, {26.685964456702496, 6.24914545508608e-5, -0.029651071618558329}, 3}, // bulk_big + {"inv_gamma_lcdf", 0.006666666666666667, 400.0, 2.0, -2.2066520386349967e-8, {0.00034032542781449647, 6.5297476510955188e-9, -1.1344180927149883e-6}, 3}, // bulk_big + {"inv_gamma_lcdf", 0.04, 200.0, 2.0, -2.0247590148473501e-57, {7.6095353053138342e-54, 2.8152855212120961e-57, -1.5219070610627669e-55}, 3}, // lower_tail + {"inv_gamma_lcdf", 0.4, 50.0, 2.0, -2.1810592140784833e-32, {2.4595912978484672e-30, 5.0484301599356797e-32, -4.9191825956969347e-31}, 3}, // lower_tail + {"inv_gamma_lcdf", 0.02, 1000.0, 2.0, -1.0270971815582085e-611, {4.62250715306945e-607, 2.3656059630799149e-611, -4.6225071530694501e-609}, 3}, // lower_deep + {"inv_gamma_lcdf", 0.06666666666666667, 500.0, 2.0, -2.9661812029194559e-409, {2.091441064732961e-405, 8.3484281304834707e-409, -6.9714702157765366e-407}, 3}, // lower_deep + {"inv_gamma_lcdf", 2000.0, 100.0, 2.0, -1.0704499119268958e-458, {5.3521965670697255e-460, 1.2329353484167989e-457, -5.3521965670697255e-457}, 3}, // lower_deep + {"inv_gamma_lccdf", 0.002, 10.0, 2.0, -1.4114888035654132e-413, {-6.9939909853802376e-408, -6.5733185134468098e-413, 6.9939909853802378e-411}, 3}, // upper_deep + {"inv_gamma_lccdf", 0.0016666666666666668, 150.0, 2.0, -1.3223466032724967e-323, {-8.3398400573340269e-318, -2.7554118260706276e-323, 6.9498667144450229e-321}, 3}, // upper_deep + {"inv_gamma_lccdf", 0.00125, 200.0, 2.0, -1.6236213327686662e-430, {-1.8199389869272201e-424, -3.3814459249480791e-430, 1.1374618668295126e-427}, 3}, // upper_deep + {"inv_gamma_lccdf", 0.006666666666666667, 10.0, 2.0, -2.8784992655370719e-114, {-1.2565975319659914e-109, -9.9465214167796963e-114, 4.1886584398866383e-112}, 3}, // upper_tail + {"inv_gamma_lccdf", 0.0033333333333333335, 150.0, 2.0, -8.141483656059317e-108, {-1.1023449539852035e-102, -1.1331647500219984e-107, 1.8372415899753393e-105}, 3}, // upper_tail + {"inv_gamma_lccdf", 0.26666666666666666, 5.0, 2.0, -0.14163482986540182, {-2.3628108352589412, -0.10565184381278135, 0.31504144470119216}, 3}, // bulk + {"inv_gamma_lccdf", 0.03333333333333333, 50.0, 2.0, -0.088182987887198289, {-38.124782664555405, -0.023247212869328028, 0.63541304440925675}, 3}, // bulk + {"inv_gamma_lccdf", 0.01, 200.0, 2.0, -0.67451500659035508, {-1107.0882797068419, -0.055400565661981331, 5.5354413985342096}, 3}, // bulk_big + {"inv_gamma_lccdf", 0.006666666666666667, 200.0, 2.0, -3.3711032560941234e-10, {-5.2016243810638526e-6, -1.4072174780714417e-10, 1.733874793687951e-8}, 3}, // bulk_big + {"inv_gamma_lccdf", 0.005, 400.0, 2.0, -0.67993656656899462, {-3148.9972259696384, -0.03937887044623289, 7.8724930649240961}, 3}, // bulk_big + {"inv_gamma_lccdf", 0.0018181818181818182, 1000.0, 2.0, -0.0010598847333686435, {-63.816625248731181, -0.00011060704645715445, 0.058015113862482892}, 3}, // bulk_big + {"inv_gamma_lccdf", 0.0022222222222222222, 1000.0, 2.0, -7.5057699942338923, {-48501.872873468241, -0.11357852886380163, 53.890969859409157}, 3}, // bulk_big + {"inv_gamma_lccdf", 0.006666666666666667, 400.0, 2.0, -17.629204302744164, {-15422.70453614961, -0.29591197274085507, 51.409015120498704}, 3}, // bulk_big + {"inv_gamma_lccdf", 0.04, 200.0, 2.0, -130.54189961219741, {-3758.2424621961886, -1.3904299230515317, 75.164849243923774}, 3}, // lower_tail + {"inv_gamma_lccdf", 0.4, 50.0, 2.0, -72.902912338991314, {-112.77049618699445, -2.3146690045591842, 22.554099237398892}, 3}, // lower_tail + {"inv_gamma_lccdf", 0.02, 1000.0, 2.0, -1406.8527552662534, {-45005.548024741408, -2.3031958470483343, 450.05548024741409}, 3}, // lower_deep + {"inv_gamma_lccdf", 0.06666666666666667, 500.0, 2.0, -940.67002769935045, {-7050.9551563285, -2.8145374673221423, 235.03183854428333}, 3}, // lower_deep + {"inv_gamma_lccdf", 2000.0, 100.0, 2.0, -1054.515893552739, {-0.049999504950543104, -11.517917229750772, 49.999504950543104}, 3}, // lower_deep + {"inv_gamma_cdf", 0.002, 10.0, 2.0, 1.4114888035654132e-413, {6.9939909853802376e-408, 6.5733185134468098e-413, -6.9939909853802378e-411}, 3}, // upper_deep + {"inv_gamma_cdf", 0.0016666666666666668, 150.0, 2.0, 1.3223466032724967e-323, {8.3398400573340269e-318, 2.7554118260706276e-323, -6.9498667144450229e-321}, 3}, // upper_deep + {"inv_gamma_cdf", 0.00125, 200.0, 2.0, 1.6236213327686662e-430, {1.8199389869272201e-424, 3.3814459249480791e-430, -1.1374618668295126e-427}, 3}, // upper_deep + {"inv_gamma_cdf", 0.006666666666666667, 10.0, 2.0, 2.8784992655370719e-114, {1.2565975319659914e-109, 9.9465214167796963e-114, -4.1886584398866383e-112}, 3}, // upper_tail + {"inv_gamma_cdf", 0.0033333333333333335, 150.0, 2.0, 8.141483656059317e-108, {1.1023449539852035e-102, 1.1331647500219984e-107, -1.8372415899753393e-105}, 3}, // upper_tail + {"inv_gamma_cdf", 0.26666666666666666, 5.0, 2.0, 0.1320618562877206, {2.0507736502979059, 0.091699265198645121, -0.27343648670638744}, 3}, // bulk + {"inv_gamma_cdf", 0.03333333333333333, 50.0, 2.0, 0.084406681093691813, {34.906796292421967, 0.021284992786349489, -0.58177993820703278}, 3}, // bulk + {"inv_gamma_cdf", 0.01, 200.0, 2.0, 0.49059658199276379, {563.95455371841642, 0.028221237507747613, -2.8197727685920822}, 3}, // bulk_big + {"inv_gamma_cdf", 0.006666666666666667, 200.0, 2.0, 3.3711032555259066e-10, {5.2016243793103313e-6, 1.4072174775970541e-10, -1.7338747931034439e-8}, 3}, // bulk_big + {"inv_gamma_cdf", 0.005, 400.0, 2.0, 0.49335087016109469, {1595.4367044026439, 0.019951270445622878, -3.9885917610066099}, 3}, // bulk_big + {"inv_gamma_cdf", 0.0018181818181818182, 1000.0, 2.0, 0.0010593232539299772, {63.749022813617866, 0.00011048987784079388, -0.057953657103288969}, 3}, // bulk_big + {"inv_gamma_cdf", 0.0022222222222222222, 1000.0, 2.0, 0.99945009773428822, {26.671289784385052, 6.2457090358415657e-5, -0.029634766427094502}, 3}, // bulk_big + {"inv_gamma_cdf", 0.006666666666666667, 400.0, 2.0, 0.99999997793347986, {0.00034032542030469856, 6.5297475070067108e-9, -1.1344180676823286e-6}, 3}, // bulk_big + {"inv_gamma_cdf", 0.04, 200.0, 2.0, 1.0, {7.6095353053138342e-54, 2.8152855212120961e-57, -1.5219070610627669e-55}, 3}, // lower_tail + {"inv_gamma_cdf", 0.4, 50.0, 2.0, 1.0, {2.4595912978484672e-30, 5.0484301599356797e-32, -4.9191825956969347e-31}, 3}, // lower_tail + {"inv_gamma_cdf", 0.02, 1000.0, 2.0, 1.0, {4.62250715306945e-607, 2.3656059630799149e-611, -4.6225071530694501e-609}, 3}, // lower_deep + {"inv_gamma_cdf", 0.06666666666666667, 500.0, 2.0, 1.0, {2.091441064732961e-405, 8.3484281304834707e-409, -6.9714702157765366e-407}, 3}, // lower_deep + {"inv_gamma_cdf", 2000.0, 100.0, 2.0, 1.0, {5.3521965670697255e-460, 1.2329353484167989e-457, -5.3521965670697255e-457}, 3}, // lower_deep + {"inv_chi_square_lcdf", 0.0005, 20.0, 0.0, -950.62299837015606, {1982018.126594686, 2.328505368530959, 0}, 2}, // upper_deep + {"inv_chi_square_lcdf", 0.0004166666666666667, 300.0, 0.0, -743.45557714901788, {2522739.5107137221, 1.0418644473588209, 0}, 2}, // upper_deep + {"inv_chi_square_lcdf", 0.0003125, 400.0, 0.0, -989.62693094236581, {4483653.7933972201, 1.0413283740186813, 0}, 2}, // upper_deep + {"inv_chi_square_lcdf", 0.0016666666666666668, 20.0, 0.0, -261.43743153134226, {174618.42662399079, 1.727726933243436, 0}, 2}, // upper_tail + {"inv_chi_square_lcdf", 0.0008333333333333334, 300.0, 0.0, -246.582217612622, {541594.13716431446, 0.69592030021373194, 0}, 2}, // upper_tail + {"inv_chi_square_lcdf", 0.06666666666666667, 10.0, 0.0, -2.024484857918011, {62.115548211890196, 0.3471830086912519, 0}, 2}, // bulk + {"inv_chi_square_lcdf", 0.008333333333333333, 100.0, 0.0, -2.4721087206353383, {1654.2195873653773, 0.1260859478808499, 0}, 2}, // bulk + {"inv_chi_square_lcdf", 0.0025, 400.0, 0.0, -0.71213311417597658, {4598.1123751631408, 0.028762162786698614, 0}, 2}, // bulk_big + {"inv_chi_square_lcdf", 0.0016666666666666668, 400.0, 0.0, -21.810610863624886, {61720.14305149316, 0.20871764685497926, 0}, 2}, // bulk_big + {"inv_chi_square_lcdf", 0.00125, 800.0, 0.0, -0.70653465390004721, {12935.513452173973, 0.020220163429637974, 0}, 2}, // bulk_big + {"inv_chi_square_lcdf", 0.00045454545454545455, 2000.0, 0.0, -6.8501250144362654, {240716.03290918323, 0.052151162277844894, 0}, 2}, // bulk_big + {"inv_chi_square_lcdf", 0.0005555555555555556, 2000.0, 0.0, -0.00055005351741434478, {106.74385782680999, 3.12457272754304e-5, 0}, 2}, // bulk_big + {"inv_chi_square_lcdf", 0.0016666666666666668, 800.0, 0.0, -2.2066520386349967e-8, {0.0013613017112579859, 3.2648738255477594e-9, 0}, 2}, // bulk_big + {"inv_chi_square_lcdf", 0.01, 400.0, 0.0, -2.0247590148473501e-57, {3.0438141221255337e-53, 1.407642760606048e-57, 0}, 2}, // lower_tail + {"inv_chi_square_lcdf", 0.1, 100.0, 0.0, -2.1810592140784833e-32, {9.8383651913938688e-30, 2.5242150799678398e-32, 0}, 2}, // lower_tail + {"inv_chi_square_lcdf", 0.005, 2000.0, 0.0, -1.0270971815582085e-611, {1.84900286122778e-606, 1.1828029815399574e-611, 0}, 2}, // lower_deep + {"inv_chi_square_lcdf", 0.016666666666666666, 1000.0, 0.0, -2.9661812029194559e-409, {8.365764258931844e-405, 4.1742140652417354e-409, 0}, 2}, // lower_deep + {"inv_chi_square_lcdf", 500.0, 200.0, 0.0, -1.0704499119268958e-458, {2.1408786268278902e-459, 6.1646767420839945e-458, 0}, 2}, // lower_deep + {"inv_chi_square_lccdf", 0.0005, 20.0, 0.0, -1.4114888035654132e-413, {-2.7975963941520951e-407, -3.2866592567234049e-413, 0}, 2}, // upper_deep + {"inv_chi_square_lccdf", 0.0004166666666666667, 300.0, 0.0, -1.3223466032724967e-323, {-3.3359360229336108e-317, -1.3777059130353138e-323, 0}, 2}, // upper_deep + {"inv_chi_square_lccdf", 0.0003125, 400.0, 0.0, -1.6236213327686662e-430, {-7.2797559477088806e-424, -1.6907229624740395e-430, 0}, 2}, // upper_deep + {"inv_chi_square_lccdf", 0.0016666666666666668, 20.0, 0.0, -2.8784992655370719e-114, {-5.0263901278639656e-109, -4.9732607083898482e-114, 0}, 2}, // upper_tail + {"inv_chi_square_lccdf", 0.0008333333333333334, 300.0, 0.0, -8.141483656059317e-108, {-4.4093798159408141e-102, -5.6658237501099918e-108, 0}, 2}, // upper_tail + {"inv_chi_square_lccdf", 0.06666666666666667, 10.0, 0.0, -0.14163482986540182, {-9.4512433410357649, -0.052825921906390677, 0}, 2}, // bulk + {"inv_chi_square_lccdf", 0.008333333333333333, 100.0, 0.0, -0.088182987887198289, {-152.49913065822162, -0.011623606434664014, 0}, 2}, // bulk + {"inv_chi_square_lccdf", 0.0025, 400.0, 0.0, -0.67451500659035508, {-4428.3531188273676, -0.027700282830990666, 0}, 2}, // bulk_big + {"inv_chi_square_lccdf", 0.0016666666666666668, 400.0, 0.0, -3.3711032560941234e-10, {-2.080649752425541e-5, -7.0360873903572084e-11, 0}, 2}, // bulk_big + {"inv_chi_square_lccdf", 0.00125, 800.0, 0.0, -0.67993656656899462, {-12595.988903878553, -0.019689435223116445, 0}, 2}, // bulk_big + {"inv_chi_square_lccdf", 0.00045454545454545455, 2000.0, 0.0, -0.0010598847333686435, {-255.26650099492473, -5.5303523228577226e-5, 0}, 2}, // bulk_big + {"inv_chi_square_lccdf", 0.0005555555555555556, 2000.0, 0.0, -7.5057699942338923, {-194007.49149387296, -0.056789264431900816, 0}, 2}, // bulk_big + {"inv_chi_square_lccdf", 0.0016666666666666668, 800.0, 0.0, -17.629204302744164, {-61690.818144598441, -0.14795598637042754, 0}, 2}, // bulk_big + {"inv_chi_square_lccdf", 0.01, 400.0, 0.0, -130.54189961219741, {-15032.969848784754, -0.69521496152576585, 0}, 2}, // lower_tail + {"inv_chi_square_lccdf", 0.1, 100.0, 0.0, -72.902912338991314, {-451.08198474797781, -1.1573345022795921, 0}, 2}, // lower_tail + {"inv_chi_square_lccdf", 0.005, 2000.0, 0.0, -1406.8527552662534, {-180022.19209896563, -1.1515979235241671, 0}, 2}, // lower_deep + {"inv_chi_square_lccdf", 0.016666666666666666, 1000.0, 0.0, -940.67002769935045, {-28203.820625314, -1.4072687336610712, 0}, 2}, // lower_deep + {"inv_chi_square_lccdf", 500.0, 200.0, 0.0, -1054.515893552739, {-0.19999801980217241, -5.7589586148753858, 0}, 2}, // lower_deep + {"inv_chi_square_cdf", 0.0005, 20.0, 0.0, 1.4114888035654132e-413, {2.7975963941520951e-407, 3.2866592567234049e-413, 0}, 2}, // upper_deep + {"inv_chi_square_cdf", 0.0004166666666666667, 300.0, 0.0, 1.3223466032724967e-323, {3.3359360229336108e-317, 1.3777059130353138e-323, 0}, 2}, // upper_deep + {"inv_chi_square_cdf", 0.0003125, 400.0, 0.0, 1.6236213327686662e-430, {7.2797559477088806e-424, 1.6907229624740395e-430, 0}, 2}, // upper_deep + {"inv_chi_square_cdf", 0.0016666666666666668, 20.0, 0.0, 2.8784992655370719e-114, {5.0263901278639656e-109, 4.9732607083898482e-114, 0}, 2}, // upper_tail + {"inv_chi_square_cdf", 0.0008333333333333334, 300.0, 0.0, 8.141483656059317e-108, {4.4093798159408141e-102, 5.6658237501099918e-108, 0}, 2}, // upper_tail + {"inv_chi_square_cdf", 0.06666666666666667, 10.0, 0.0, 0.1320618562877206, {8.2030946011916234, 0.04584963259932256, 0}, 2}, // bulk + {"inv_chi_square_cdf", 0.008333333333333333, 100.0, 0.0, 0.084406681093691813, {139.62718516968787, 0.010642496393174745, 0}, 2}, // bulk + {"inv_chi_square_cdf", 0.0025, 400.0, 0.0, 0.49059658199276379, {2255.8182148736657, 0.014110618753873806, 0}, 2}, // bulk_big + {"inv_chi_square_cdf", 0.0016666666666666668, 400.0, 0.0, 3.3711032555259066e-10, {2.0806497517241325e-5, 7.0360873879852707e-11, 0}, 2}, // bulk_big + {"inv_chi_square_cdf", 0.00125, 800.0, 0.0, 0.49335087016109469, {6381.7468176105757, 0.0099756352228114392, 0}, 2}, // bulk_big + {"inv_chi_square_cdf", 0.00045454545454545455, 2000.0, 0.0, 0.0010593232539299772, {254.99609125447146, 5.5244938920396938e-5, 0}, 2}, // bulk_big + {"inv_chi_square_cdf", 0.0005555555555555556, 2000.0, 0.0, 0.99945009773428822, {106.68515913754021, 3.1228545179207828e-5, 0}, 2}, // bulk_big + {"inv_chi_square_cdf", 0.0016666666666666668, 800.0, 0.0, 0.99999997793347986, {0.0013613016812187942, 3.2648737535033554e-9, 0}, 2}, // bulk_big + {"inv_chi_square_cdf", 0.01, 400.0, 0.0, 1.0, {3.0438141221255337e-53, 1.407642760606048e-57, 0}, 2}, // lower_tail + {"inv_chi_square_cdf", 0.1, 100.0, 0.0, 1.0, {9.8383651913938688e-30, 2.5242150799678398e-32, 0}, 2}, // lower_tail + {"inv_chi_square_cdf", 0.005, 2000.0, 0.0, 1.0, {1.84900286122778e-606, 1.1828029815399574e-611, 0}, 2}, // lower_deep + {"inv_chi_square_cdf", 0.016666666666666666, 1000.0, 0.0, 1.0, {8.365764258931844e-405, 4.1742140652417354e-409, 0}, 2}, // lower_deep + {"inv_chi_square_cdf", 500.0, 200.0, 0.0, 1.0, {2.1408786268278902e-459, 6.1646767420839945e-458, 0}, 2}, // lower_deep + {"scaled_inv_chi_square_cdf", 0.01, 20.0, 1.0, 1.4114888035654132e-413, {1.3987981970760475e-408, -6.6653250597078973e-412, -2.7975963941520951e-410}, 3}, // upper_deep + {"scaled_inv_chi_square_cdf", 0.125, 300.0, 1.0, 1.3223466032724075e-323, {1.1119786743111287e-319, -3.2555385632611488e-323, -2.7799466857778218e-320}, 3}, // upper_deep + {"scaled_inv_chi_square_cdf", 0.125, 400.0, 1.0, 1.6236213327686189e-430, {1.8199389869271671e-426, -3.9965863716734071e-430, -4.5498474673179179e-427}, 3}, // upper_deep + {"scaled_inv_chi_square_cdf", 0.03333333333333333, 20.0, 1.0, 2.8784992655370065e-114, {2.5131950639319261e-110, -3.6913323690475699e-113, -1.6754633759546174e-111}, 3}, // upper_tail + {"scaled_inv_chi_square_cdf", 0.25, 300.0, 1.0, 8.1414836560590812e-108, {1.469793271980229e-104, -6.5824535163920805e-108, -7.3489663599011451e-105}, 3}, // upper_tail + {"scaled_inv_chi_square_cdf", 0.6666666666666666, 10.0, 1.0, 0.13206185628772058, {0.82030946011916229, -0.0088376647419549284, -1.0937459468255497}, 3}, // bulk + {"scaled_inv_chi_square_cdf", 0.8333333333333334, 100.0, 1.0, 0.084406681093691881, {1.3962718516968794, -0.00099310237096591134, -2.3271197528281325}, 3}, // bulk + {"scaled_inv_chi_square_cdf", 1.0, 400.0, 1.0, 0.49059658199276367, {5.6395455371841644, 1.1754910913395298e-5, -11.279091074368329}, 3}, // bulk_big + {"scaled_inv_chi_square_cdf", 0.6666666666666666, 400.0, 1.0, 3.3711032555258651e-10, {5.2016243793102697e-8, -1.6332865775319296e-11, -6.9354991724136925e-8}, 3}, // bulk_big + {"scaled_inv_chi_square_cdf", 1.0, 800.0, 1.0, 0.49335087016109453, {7.9771835220132198, 4.1558202949142874e-6, -15.95436704402644}, 3}, // bulk_big + {"scaled_inv_chi_square_cdf", 0.9090909090909091, 2000.0, 1.0, 0.0010593232539299735, {0.12749804562723532, -2.708718182892023e-6, -0.23181462841315512}, 3}, // bulk_big + {"scaled_inv_chi_square_cdf", 1.1111111111111112, 2000.0, 1.0, 0.99945009773428822, {0.053342579568769859, 1.5937787521133202e-6, -0.11853906570837747}, 3}, // bulk_big + {"scaled_inv_chi_square_cdf", 1.3333333333333333, 800.0, 1.0, 0.99999997793347986, {1.7016271015235134e-6, 4.288285842975388e-10, -4.5376722707293688e-6}, 3}, // bulk_big + {"scaled_inv_chi_square_cdf", 4.0, 400.0, 1.0, 1.0, {7.6095353053138581e-56, 6.4668923007466662e-58, -6.0876282442510865e-55}, 3}, // lower_tail + {"scaled_inv_chi_square_cdf", 10.0, 100.0, 1.0, 1.0, {9.8383651913938939e-32, 1.5403785608284567e-32, -1.9676730382787788e-30}, 3}, // lower_tail + {"scaled_inv_chi_square_cdf", 10.0, 2000.0, 1.0, 1.0, {9.2450143061390734e-610, 7.2055226623302592e-612, -1.8490028612278147e-608}, 3}, // lower_deep + {"scaled_inv_chi_square_cdf", 16.666666666666668, 1000.0, 1.0, 1.0, {8.3657642589315093e-408, 2.7799200220863172e-409, -2.7885880863105033e-406}, 3}, // lower_deep + {"scaled_inv_chi_square_cdf", 100000.0, 200.0, 1.0, 1.0, {1.0704393134139451e-461, 5.6294570853770219e-458, -2.1408786268278902e-456}, 3}, // lower_deep + {"scaled_inv_chi_square_lcdf", 0.01, 20.0, 1.0, -950.62299837015606, {99100.906329734298, -47.221947796336191, -1982.018126594686}, 3}, // upper_deep + {"scaled_inv_chi_square_lcdf", 0.125, 300.0, 1.0, -743.45557714901795, {8409.1317023790747, -2.4619404286324602, -2102.2829255947687}, 3}, // upper_deep + {"scaled_inv_chi_square_lcdf", 0.125, 400.0, 1.0, -989.62693094236584, {11209.134483493051, -2.461526152072897, -2802.2836208732627}, 3}, // upper_deep + {"scaled_inv_chi_square_lcdf", 0.03333333333333333, 20.0, 1.0, -261.43743153134228, {8730.9213311995407, -12.823808618755798, -582.06142207996937}, 3}, // upper_tail + {"scaled_inv_chi_square_lcdf", 0.25, 300.0, 1.0, -246.58221761262203, {1805.3137905477151, -0.80850785857603064, -902.65689527385757}, 3}, // upper_tail + {"scaled_inv_chi_square_lcdf", 0.6666666666666666, 10.0, 1.0, -2.0244848579180111, {6.2115548211890203, -0.06692064605468275, -8.2820730949186933}, 3}, // bulk + {"scaled_inv_chi_square_lcdf", 0.8333333333333334, 100.0, 1.0, -2.4721087206353375, {16.542195873653769, -0.011765684399598205, -27.570326456089616}, 3}, // bulk + {"scaled_inv_chi_square_lcdf", 1.0, 400.0, 1.0, -0.71213311417597682, {11.495280937907855, 2.3960441928983279e-5, -22.99056187581571}, 3}, // bulk_big + {"scaled_inv_chi_square_lcdf", 0.6666666666666666, 400.0, 1.0, -21.810610863624899, {154.30035762873297, -0.048449615859575622, -205.73381017164395}, 3}, // bulk_big + {"scaled_inv_chi_square_lcdf", 1.0, 800.0, 1.0, -0.70653465390004754, {16.169391815217472, 8.4236606161397501e-6, -32.338783630434945}, 3}, // bulk_big + {"scaled_inv_chi_square_lcdf", 0.9090909090909091, 2000.0, 1.0, -6.8501250144362689, {120.35801645459165, -0.0025570270196967494, -218.83275719016664}, 3}, // bulk_big + {"scaled_inv_chi_square_lcdf", 1.1111111111111112, 2000.0, 1.0, -0.0005500535174143421, {0.053371928913404749, 1.5946556568720643e-6, -0.11860428647423278}, 3}, // bulk_big + {"scaled_inv_chi_square_lcdf", 1.3333333333333333, 800.0, 1.0, -2.2066520386350239e-8, {1.7016271390725029e-6, 4.288285937602936e-10, -4.5376723708600075e-6}, 3}, // bulk_big + {"scaled_inv_chi_square_lcdf", 4.0, 400.0, 1.0, -2.0247590148473564e-57, {7.6095353053138581e-56, 6.4668923007466662e-58, -6.0876282442510865e-55}, 3}, // lower_tail + {"scaled_inv_chi_square_lcdf", 10.0, 100.0, 1.0, -2.1810592140784888e-32, {9.8383651913938939e-32, 1.5403785608284567e-32, -1.9676730382787788e-30}, 3}, // lower_tail + {"scaled_inv_chi_square_lcdf", 10.0, 2000.0, 1.0, -1.0270971815582278e-611, {9.2450143061390734e-610, 7.2055226623302592e-612, -1.8490028612278147e-608}, 3}, // lower_deep + {"scaled_inv_chi_square_lcdf", 16.666666666666668, 1000.0, 1.0, -2.9661812029193374e-409, {8.3657642589315093e-408, 2.7799200220863172e-409, -2.7885880863105033e-406}, 3}, // lower_deep + {"scaled_inv_chi_square_lcdf", 100000.0, 200.0, 1.0, -1.0704499119268958e-458, {1.0704393134139451e-461, 5.6294570853770219e-458, -2.1408786268278902e-456}, 3}, // lower_deep + {"scaled_inv_chi_square_lccdf", 0.01, 20.0, 1.0, -1.4114888035654132e-413, {-1.3987981970760475e-408, 6.6653250597078973e-412, 2.7975963941520951e-410}, 3}, // upper_deep + {"scaled_inv_chi_square_lccdf", 0.125, 300.0, 1.0, -1.3223466032724075e-323, {-1.1119786743111287e-319, 3.2555385632611488e-323, 2.7799466857778218e-320}, 3}, // upper_deep + {"scaled_inv_chi_square_lccdf", 0.125, 400.0, 1.0, -1.6236213327686189e-430, {-1.8199389869271671e-426, 3.9965863716734071e-430, 4.5498474673179179e-427}, 3}, // upper_deep + {"scaled_inv_chi_square_lccdf", 0.03333333333333333, 20.0, 1.0, -2.8784992655370065e-114, {-2.5131950639319261e-110, 3.6913323690475699e-113, 1.6754633759546174e-111}, 3}, // upper_tail + {"scaled_inv_chi_square_lccdf", 0.25, 300.0, 1.0, -8.1414836560590812e-108, {-1.469793271980229e-104, 6.5824535163920805e-108, 7.3489663599011451e-105}, 3}, // upper_tail + {"scaled_inv_chi_square_lccdf", 0.6666666666666666, 10.0, 1.0, -0.14163482986540179, {-0.9451243341035764, 0.010182367033847754, 1.2601657788047685}, 3}, // bulk + {"scaled_inv_chi_square_lccdf", 0.8333333333333334, 100.0, 1.0, -0.088182987887198363, {-1.5249913065822171, 0.0010846544535211213, 2.5416521776370287}, 3}, // bulk + {"scaled_inv_chi_square_lccdf", 1.0, 400.0, 1.0, -0.67451500659035485, {-11.070882797068417, -2.3075838319617074e-5, 22.141765594136833}, 3}, // bulk_big + {"scaled_inv_chi_square_lccdf", 0.6666666666666666, 400.0, 1.0, -3.3711032560940819e-10, {-5.201624381063791e-8, 1.6332865780825274e-11, 6.9354991747517209e-8}, 3}, // bulk_big + {"scaled_inv_chi_square_lccdf", 1.0, 800.0, 1.0, -0.6799365665689943, {-15.744986129848187, -8.2025608062046313e-6, 31.489972259696374}, 3}, // bulk_big + {"scaled_inv_chi_square_lccdf", 0.9090909090909091, 2000.0, 1.0, -0.0010598847333686397, {-0.12763325049746195, 2.711590633905658e-6, 0.23206045544993081}, 3}, // bulk_big + {"scaled_inv_chi_square_lccdf", 1.1111111111111112, 2000.0, 1.0, -7.5057699942338972, {-97.003745746936512, -0.0028982945724916613, 215.5638794376367}, 3}, // bulk_big + {"scaled_inv_chi_square_lccdf", 1.3333333333333333, 800.0, 1.0, -17.629204302744151, {-77.113522680748034, -0.019433448569180762, 205.63606048199475}, 3}, // bulk_big + {"scaled_inv_chi_square_lccdf", 4.0, 400.0, 1.0, -130.54189961219741, {-37.582424621961887, -0.31939071530614697, 300.65939697569509}, 3}, // lower_tail + {"scaled_inv_chi_square_lccdf", 10.0, 100.0, 1.0, -72.902912338991312, {-4.5108198474797784, -0.70625251753161426, 90.216396949595567}, 3}, // lower_tail + {"scaled_inv_chi_square_lccdf", 10.0, 2000.0, 1.0, -1406.8527552662533, {-90.011096049482818, -0.70154244327675303, 1800.2219209896564}, 3}, // lower_deep + {"scaled_inv_chi_square_lccdf", 16.666666666666668, 1000.0, 1.0, -940.67002769935049, {-28.203820625313998, -0.93720505657250454, 940.12735417713332}, 3}, // lower_deep + {"scaled_inv_chi_square_lccdf", 100000.0, 200.0, 1.0, -1054.515893552739, {-0.00099999009901086207, -5.2589635653699547, 199.99801980217241}, 3}, // lower_deep + {"poisson_lcdf", 9.0, 1000.0, 0.0, -950.62299837015608, {-0.99100906329734302, 0, 0}, 1}, // upper_deep + {"poisson_lcdf", 149.0, 1200.0, 0.0, -743.45557714901795, {-0.87595121899782028, 0, 0}, 1}, // upper_deep + {"poisson_lcdf", 199.0, 1600.0, 0.0, -989.62693094236584, {-0.8757136315228946, 0, 0}, 1}, // upper_deep + {"poisson_lcdf", 9.0, 300.0, 0.0, -261.43743153134228, {-0.97010237013328227, 0, 0}, 1}, // upper_tail + {"poisson_lcdf", 149.0, 600.0, 0.0, -246.58221761262203, {-0.75221407939488131, 0, 0}, 1}, // upper_tail + {"poisson_lcdf", 4.0, 7.5, 0.0, -2.0244848579180109, {-0.55213820632791284, 0, 0}, 1}, // bulk + {"poisson_lcdf", 49.0, 60.0, 0.0, -2.4721087206353381, {-0.22975272046741351, 0, 0}, 1}, // bulk + {"poisson_lcdf", 199.0, 200.0, 0.0, -0.71213311417597682, {-0.057476404689539276, 0, 0}, 1}, // bulk_big + {"poisson_lcdf", 199.0, 300.0, 0.0, -21.810610863624893, {-0.34288968361940653, 0, 0}, 1}, // bulk_big + {"poisson_lcdf", 399.0, 400.0, 0.0, -0.70653465390004754, {-0.040423479538043681, 0, 0}, 1}, // bulk_big + {"poisson_lcdf", 999.0, 1100.0, 0.0, -6.8501250144362653, {-0.099469435086439349, 0, 0}, 1}, // bulk_big + {"poisson_lcdf", 999.0, 900.0, 0.0, -0.00055005351741434473, {-6.5891270263462948e-5, 0, 0}, 1}, // bulk_big + {"poisson_lcdf", 399.0, 300.0, 0.0, -2.2066520386350113e-8, {-7.5627872847666373e-9, 0, 0}, 1}, // bulk_big + {"poisson_lcdf", 199.0, 50.0, 0.0, -2.0247590148473564e-57, {-6.0876282442510865e-57, 0, 0}, 1}, // lower_tail + {"poisson_lcdf", 49.0, 5.0, 0.0, -2.1810592140784888e-32, {-1.9676730382787788e-31, 0, 0}, 1}, // lower_tail + {"poisson_lcdf", 999.0, 100.0, 0.0, -1.0270971815582278e-611, {-9.2450143061390734e-611, 0, 0}, 1}, // lower_deep + {"poisson_lcdf", 499.0, 30.0, 0.0, -2.9661812029194365e-409, {-4.6476468105176607e-408, 0, 0}, 1}, // lower_deep + {"poisson_lcdf", 99.0, 0.001, 0.0, -1.070449911926898e-458, {-1.0704393134139473e-453, 0, 0}, 1}, // lower_deep + {"poisson_lccdf", 9.0, 1000.0, 0.0, -1.4114888035653841e-413, {1.3987981970760187e-413, 0, 0}, 1}, // upper_deep + {"poisson_lccdf", 149.0, 1200.0, 0.0, -1.3223466032724075e-323, {1.1583111190740924e-323, 0, 0}, 1}, // upper_deep + {"poisson_lccdf", 199.0, 1600.0, 0.0, -1.6236213327686189e-430, {1.4218273335368493e-430, 0, 0}, 1}, // upper_deep + {"poisson_lccdf", 9.0, 300.0, 0.0, -2.8784992655370182e-114, {2.7924389599243735e-114, 0, 0}, 1}, // upper_tail + {"poisson_lccdf", 149.0, 600.0, 0.0, -8.1414836560590812e-108, {6.1241386332509542e-108, 0, 0}, 1}, // upper_tail + {"poisson_lccdf", 4.0, 7.5, 0.0, -0.14163482986540183, {0.084011051920317913, 0, 0}, 1}, // bulk + {"poisson_lccdf", 49.0, 60.0, 0.0, -0.088182987887198306, {0.021180434813641895, 0, 0}, 1}, // bulk + {"poisson_lccdf", 199.0, 200.0, 0.0, -0.67451500659035485, {0.055354413985342084, 0, 0}, 1}, // bulk_big + {"poisson_lccdf", 199.0, 300.0, 0.0, -3.3711032560941012e-10, {1.1559165291252932e-10, 0, 0}, 1}, // bulk_big + {"poisson_lccdf", 399.0, 400.0, 0.0, -0.6799365665689943, {0.039362465324620467, 0, 0}, 1}, // bulk_big + {"poisson_lccdf", 999.0, 1100.0, 0.0, -0.0010598847333686436, {0.00010548202520451436, 0, 0}, 1}, // bulk_big + {"poisson_lccdf", 999.0, 900.0, 0.0, -7.5057699942338924, {0.11975771079868702, 0, 0}, 1}, // bulk_big + {"poisson_lccdf", 399.0, 300.0, 0.0, -17.629204302744157, {0.3427267674699913, 0, 0}, 1}, // bulk_big + {"poisson_lccdf", 199.0, 50.0, 0.0, -130.54189961219741, {3.0065939697569509, 0, 0}, 1}, // lower_tail + {"poisson_lccdf", 49.0, 5.0, 0.0, -72.902912338991312, {9.0216396949595567, 0, 0}, 1}, // lower_tail + {"poisson_lccdf", 999.0, 100.0, 0.0, -1406.8527552662533, {9.0011096049482818, 0, 0}, 1}, // lower_deep + {"poisson_lccdf", 499.0, 30.0, 0.0, -940.67002769935046, {15.668789236285555, 0, 0}, 1}, // lower_deep + {"poisson_lccdf", 99.0, 0.001, 0.0, -1054.515893552739, {99999.009901086205, 0, 0}, 1}, // lower_deep + {"poisson_cdf", 9.0, 1000.0, 0.0, 1.4114888035653841e-413, {-1.3987981970760187e-413, 0, 0}, 1}, // upper_deep + {"poisson_cdf", 149.0, 1200.0, 0.0, 1.3223466032724075e-323, {-1.1583111190740924e-323, 0, 0}, 1}, // upper_deep + {"poisson_cdf", 199.0, 1600.0, 0.0, 1.6236213327686189e-430, {-1.4218273335368493e-430, 0, 0}, 1}, // upper_deep + {"poisson_cdf", 9.0, 300.0, 0.0, 2.8784992655370182e-114, {-2.7924389599243735e-114, 0, 0}, 1}, // upper_tail + {"poisson_cdf", 149.0, 600.0, 0.0, 8.1414836560590812e-108, {-6.1241386332509542e-108, 0, 0}, 1}, // upper_tail + {"poisson_cdf", 4.0, 7.5, 0.0, 0.13206185628772061, {-0.072916396455036654, 0, 0}, 1}, // bulk + {"poisson_cdf", 49.0, 60.0, 0.0, 0.08440668109369183, {-0.019392664606901095, 0, 0}, 1}, // bulk + {"poisson_cdf", 199.0, 200.0, 0.0, 0.49059658199276367, {-0.028197727685920822, 0, 0}, 1}, // bulk_big + {"poisson_cdf", 199.0, 300.0, 0.0, 3.3711032555258843e-10, {-1.1559165287356218e-10, 0, 0}, 1}, // bulk_big + {"poisson_cdf", 399.0, 400.0, 0.0, 0.49335087016109453, {-0.01994295880503305, 0, 0}, 1}, // bulk_big + {"poisson_cdf", 999.0, 1100.0, 0.0, 0.0010593232539299773, {-0.00010537028564234359, 0, 0}, 1}, // bulk_big + {"poisson_cdf", 999.0, 900.0, 0.0, 0.99945009773428822, {-6.5855036504654443e-5, 0, 0}, 1}, // bulk_big + {"poisson_cdf", 399.0, 300.0, 0.0, 0.99999997793347986, {-7.5627871178822393e-9, 0, 0}, 1}, // bulk_big + {"poisson_cdf", 199.0, 50.0, 0.0, 1.0, {-6.0876282442510865e-57, 0, 0}, 1}, // lower_tail + {"poisson_cdf", 49.0, 5.0, 0.0, 1.0, {-1.9676730382787788e-31, 0, 0}, 1}, // lower_tail + {"poisson_cdf", 999.0, 100.0, 0.0, 1.0, {-9.2450143061390734e-611, 0, 0}, 1}, // lower_deep + {"poisson_cdf", 499.0, 30.0, 0.0, 1.0, {-4.6476468105176607e-408, 0, 0}, 1}, // lower_deep + {"poisson_cdf", 99.0, 0.001, 0.0, 1.0, {-1.0704393134139473e-453, 0, 0}, 1}, // lower_deep + }; + // NOLINTEND + // clang-format on + return r; +} + +} // namespace incomplete_gamma_cdfs +} // namespace test +} // namespace stan + +#endif diff --git a/test/unit/math/rev/prob/incomplete_gamma_cdfs_reference_test.cpp b/test/unit/math/rev/prob/incomplete_gamma_cdfs_reference_test.cpp new file mode 100644 index 00000000000..09d185d9ff4 --- /dev/null +++ b/test/unit/math/rev/prob/incomplete_gamma_cdfs_reference_test.cpp @@ -0,0 +1,125 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + +using stan::math::var; +using F3 = std::function; +using F2 = std::function; +using FP = std::function; + +const std::map& three_arg() { + namespace sm = stan::math; + static const std::map f = { + {"gamma_cdf", + [](auto& y, auto& a, auto& b) { return sm::gamma_cdf(y, a, b); }}, + {"gamma_lcdf", + [](auto& y, auto& a, auto& b) { return sm::gamma_lcdf(y, a, b); }}, + {"gamma_lccdf", + [](auto& y, auto& a, auto& b) { return sm::gamma_lccdf(y, a, b); }}, + {"inv_gamma_cdf", + [](auto& y, auto& a, auto& b) { return sm::inv_gamma_cdf(y, a, b); }}, + {"inv_gamma_lcdf", + [](auto& y, auto& a, auto& b) { return sm::inv_gamma_lcdf(y, a, b); }}, + {"inv_gamma_lccdf", + [](auto& y, auto& a, auto& b) { return sm::inv_gamma_lccdf(y, a, b); }}, + {"scaled_inv_chi_square_cdf", + [](auto& y, auto& n, auto& s) { + return sm::scaled_inv_chi_square_cdf(y, n, s); + }}, + {"scaled_inv_chi_square_lcdf", + [](auto& y, auto& n, auto& s) { + return sm::scaled_inv_chi_square_lcdf(y, n, s); + }}, + {"scaled_inv_chi_square_lccdf", + [](auto& y, auto& n, auto& s) { + return sm::scaled_inv_chi_square_lccdf(y, n, s); + }}, + }; + return f; +} + +const std::map& two_arg() { + namespace sm = stan::math; + static const std::map f = { + {"chi_square_cdf", + [](auto& y, auto& n) { return sm::chi_square_cdf(y, n); }}, + {"chi_square_lcdf", + [](auto& y, auto& n) { return sm::chi_square_lcdf(y, n); }}, + {"chi_square_lccdf", + [](auto& y, auto& n) { return sm::chi_square_lccdf(y, n); }}, + {"inv_chi_square_cdf", + [](auto& y, auto& n) { return sm::inv_chi_square_cdf(y, n); }}, + {"inv_chi_square_lcdf", + [](auto& y, auto& n) { return sm::inv_chi_square_lcdf(y, n); }}, + {"inv_chi_square_lccdf", + [](auto& y, auto& n) { return sm::inv_chi_square_lccdf(y, n); }}, + }; + return f; +} + +const std::map& poisson() { + namespace sm = stan::math; + static const std::map f = { + {"poisson_cdf", [](int n, auto& l) { return sm::poisson_cdf(n, l); }}, + {"poisson_lcdf", [](int n, auto& l) { return sm::poisson_lcdf(n, l); }}, + {"poisson_lccdf", [](int n, auto& l) { return sm::poisson_lccdf(n, l); }}, + }; + return f; +} + +/** + * Relative error; a reference below 1e-300 (zero or subnormal, with few + * significant bits) needs an absolute error below 1e-300. + */ +void expect_close(double expected, double actual, double tol, + const std::string& what) { + ASSERT_TRUE(std::isfinite(actual)) << what << ": not finite, " << actual; + if (std::fabs(expected) < 1e-300) { + EXPECT_LE(std::fabs(actual - expected), 1e-300) + << what << ": expected " << expected << ", got " << actual; + return; + } + const double err = std::fabs(actual - expected) / std::fabs(expected); + EXPECT_LT(err, tol) << what << ": relative error " << err << ", expected " + << expected << ", got " << actual; +} + +} // namespace + +TEST_F(AgradRev, ProbDistributionsIncompleteGammaCdfs_reference) { + namespace ref = stan::test::incomplete_gamma_cdfs; + for (const auto& r : ref::rows()) { + const std::string at = r.function + " at (" + std::to_string(r.y) + ", " + + std::to_string(r.p1) + ", " + std::to_string(r.p2) + + ")"; + std::vector args; + var lp; + if (poisson().count(r.function)) { + args = {var(r.p1)}; + lp = poisson().at(r.function)(static_cast(r.y), args[0]); + } else if (two_arg().count(r.function)) { + args = {var(r.y), var(r.p1)}; + lp = two_arg().at(r.function)(args[0], args[1]); + } else { + ASSERT_TRUE(three_arg().count(r.function)) << r.function; + args = {var(r.y), var(r.p1), var(r.p2)}; + lp = three_arg().at(r.function)(args[0], args[1], args[2]); + } + ASSERT_EQ(static_cast(args.size()), r.n_grad) << at; + lp.grad(); + expect_close(r.value, lp.val(), 1e-12, "value of " + at); + for (int k = 0; k < r.n_grad; ++k) { + expect_close(r.grad[k], args[k].adj(), 1e-9, + "gradient " + std::to_string(k) + " of " + at); + } + stan::math::recover_memory(); + } +} From 22d781f290cb909b7cfc7cd086a144e043750f4d Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Wed, 21 Oct 2026 23:56:55 +0300 Subject: [PATCH 9/9] Return the reverse-mode gamma_p derivatives for z / a > 10 and keep other elements of poisson_lccdf for negative counts --- stan/math/prim/prob/poisson_lccdf.hpp | 22 ++++----- stan/math/rev/fun/gamma_p.hpp | 22 --------- test/unit/math/mix/fun/gamma_p_test.cpp | 48 +++++++++++++++++++ .../math/prim/prob/poisson_ccdf_log_test.cpp | 13 +++++ 4 files changed, 71 insertions(+), 34 deletions(-) diff --git a/stan/math/prim/prob/poisson_lccdf.hpp b/stan/math/prim/prob/poisson_lccdf.hpp index b5248360af1..11c55430a39 100644 --- a/stan/math/prim/prob/poisson_lccdf.hpp +++ b/stan/math/prim/prob/poisson_lccdf.hpp @@ -7,10 +7,8 @@ #include #include #include -#include #include #include -#include #include #include #include @@ -47,7 +45,6 @@ inline return_type_t poisson_lccdf(const T_n& n, const T_rate& lambda) { T_n_ref n_ref = n; T_lambda_ref lambda_ref = lambda; - decltype(auto) n_val = to_ref(as_value_column_array_or_scalar(n_ref)); decltype(auto) lambda_val = to_ref(as_value_column_array_or_scalar(lambda_ref)); @@ -59,26 +56,27 @@ inline return_type_t poisson_lccdf(const T_n& n, const T_rate& lambda) { auto ops_partials = make_partials_propagator(lambda_ref); - if (sum(promote_scalar(n_val < 0))) { - return ops_partials.build(0.0); - } - scalar_seq_view n_vec(n_ref); scalar_seq_view lambda_vec(lambda_ref); const size_t N = max_size(n, lambda); - // An element with probability 0 (lambda = 0, or n = INT_MAX, which - // stands for an infinite count) makes the result -inf; the gradients are - // treated as zero + // An element with probability 0 (n >= 0 with lambda = 0, or n = INT_MAX, + // which stands for an infinite count) makes the result -inf; the + // gradients are treated as zero for (size_t i = 0; i < N; i++) { - if (value_of_rec(lambda_vec.val(i)) == 0 - || n_vec.val(i) == std::numeric_limits::max()) { + if (n_vec.val(i) >= 0 + && (value_of_rec(lambda_vec.val(i)) == 0 + || n_vec.val(i) == std::numeric_limits::max())) { return ops_partials.build(negative_infinity()); } } T_partials_return log_P(0.0); for (size_t i = 0; i < N; i++) { + // An element with probability 1 (n < 0) contributes nothing + if (n_vec.val(i) < 0) { + continue; + } const T_partials_return lambda_dbl = lambda_vec.val(i); const auto tail = internal::log_gamma_tail>( diff --git a/stan/math/rev/fun/gamma_p.hpp b/stan/math/rev/fun/gamma_p.hpp index e25d60132b7..698bc63be4a 100644 --- a/stan/math/rev/fun/gamma_p.hpp +++ b/stan/math/rev/fun/gamma_p.hpp @@ -19,10 +19,6 @@ class gamma_p_vv_vari : public op_vv_vari { gamma_p_vv_vari(vari* avi, vari* bvi) : op_vv_vari(gamma_p(avi->val_, bvi->val_), avi, bvi) {} void chain() { - using std::exp; - using std::fabs; - using std::log; - if (is_inf(avi_->val_)) { avi_->adj_ = NOT_A_NUMBER; bvi_->adj_ = NOT_A_NUMBER; @@ -34,12 +30,6 @@ class gamma_p_vv_vari : public op_vv_vari { return; } - // return zero derivative as gamma_p is flat - // to machine precision for b / a > 10 - if (std::fabs(bvi_->val_ / avi_->val_) > 10) { - return; - } - avi_->adj_ += adj_ * grad_reg_lower_inc_gamma(avi_->val_, bvi_->val_); bvi_->adj_ += adj_ @@ -62,12 +52,6 @@ class gamma_p_vd_vari : public op_vd_vari { return; } - // return zero derivative as gamma_p is flat - // to machine precision for b / a > 10 - if (std::fabs(bd_ / avi_->val_) > 10) { - return; - } - avi_->adj_ += adj_ * grad_reg_lower_inc_gamma(avi_->val_, bd_); } }; @@ -86,12 +70,6 @@ class gamma_p_dv_vari : public op_dv_vari { return; } - // return zero derivative as gamma_p is flat to - // machine precision for b / a > 10 - if (std::fabs(bvi_->val_ / ad_) > 10) { - return; - } - bvi_->adj_ += adj_ * std::exp(-bvi_->val_ + (ad_ - 1.0) * std::log(bvi_->val_) - lgamma(ad_)); diff --git a/test/unit/math/mix/fun/gamma_p_test.cpp b/test/unit/math/mix/fun/gamma_p_test.cpp index dedf35b5f5f..af952b9caf9 100644 --- a/test/unit/math/mix/fun/gamma_p_test.cpp +++ b/test/unit/math/mix/fun/gamma_p_test.cpp @@ -105,3 +105,51 @@ TEST(mathMixScalFun, gammaP_vec) { in2 << 0.5, 3.4; stan::test::expect_ad_vectorized_binary(f, in1, in2); } + +// Reverse-mode derivatives far in the upper tail, z / a > 10, where both +// derivatives are about 1e-13. References from mpmath at 50 digits. +TEST(mathMixScalFun, gammaP_rev_far_upper_tail) { + using stan::math::gamma_p; + using stan::math::var; + const double d_a = -4.2728631691092249555e-13; + const double d_z = 1.5030507939978329636e-13; + { + var a = 3.0; + var z = 36.0; + var p = gamma_p(a, z); + p.grad(); + stan::test::expect_near_rel("vv d/da", d_a, a.adj(), 1e-10); + stan::test::expect_near_rel("vv d/dz", d_z, z.adj(), 1e-10); + stan::math::recover_memory(); + } + { + var a = 3.0; + var p = gamma_p(a, 36.0); + p.grad(); + stan::test::expect_near_rel("vd d/da", d_a, a.adj(), 1e-10); + stan::math::recover_memory(); + } + { + var z = 36.0; + var p = gamma_p(3.0, z); + p.grad(); + stan::test::expect_near_rel("dv d/dz", d_z, z.adj(), 1e-10); + stan::math::recover_memory(); + } + // log1m(gamma_p(alpha, beta * y)) = log Q at y = 12, alpha = 3, beta = 3. + // 1 - P is about 1.6e-13 and keeps only about three digits in double. + { + var y = 12.0; + var alpha = 3.0; + var beta = 3.0; + var lp = stan::math::log1m(gamma_p(alpha, beta * y)); + lp.grad(); + stan::test::expect_near_rel("log1m d/dy", -2.8379562043795620438, y.adj(), + 1e-2); + stan::test::expect_near_rel("log1m d/dalpha", 2.6892412395869051037, + alpha.adj(), 1e-2); + stan::test::expect_near_rel("log1m d/dbeta", -11.351824817518248175, + beta.adj(), 1e-2); + stan::math::recover_memory(); + } +} diff --git a/test/unit/math/prim/prob/poisson_ccdf_log_test.cpp b/test/unit/math/prim/prob/poisson_ccdf_log_test.cpp index ba8c758eb8d..3f32b58b676 100644 --- a/test/unit/math/prim/prob/poisson_ccdf_log_test.cpp +++ b/test/unit/math/prim/prob/poisson_ccdf_log_test.cpp @@ -1,5 +1,6 @@ #include #include +#include TEST(ProbPoisson, ccdf_log_matches_lccdf) { int y = 3; @@ -10,3 +11,15 @@ TEST(ProbPoisson, ccdf_log_matches_lccdf) { EXPECT_FLOAT_EQ((stan::math::poisson_lccdf(y, lambda)), (stan::math::poisson_ccdf_log(y, lambda))); } + +// A negative count has P(N > n) = 1 and contributes 0; the other elements +// of the vector are kept. +TEST(ProbPoisson, lccdf_negative_count_contributes_zero) { + using stan::math::poisson_lccdf; + EXPECT_EQ(0.0, poisson_lccdf(-1, 5.0)); + EXPECT_EQ(0.0, poisson_lccdf(-1, 0.0)); + std::vector n{-1, 1}; + EXPECT_DOUBLE_EQ(poisson_lccdf(1, 5.0), poisson_lccdf(n, 5.0)); + std::vector lambda{0.0, 5.0}; + EXPECT_DOUBLE_EQ(poisson_lccdf(1, 5.0), poisson_lccdf(n, lambda)); +}