1 #ifndef STAN_MATH_PRIM_FUN_GRAD_REG_INC_GAMMA_HPP
2 #define STAN_MATH_PRIM_FUN_GRAD_REG_INC_GAMMA_HPP
3
4 #include <stan/math/prim/meta.hpp>
5 #include <stan/math/prim/err.hpp>
6 #include <stan/math/prim/fun/constants.hpp>
7 #include <stan/math/prim/fun/exp.hpp>
8 #include <stan/math/prim/fun/gamma_p.hpp>
9 #include <stan/math/prim/fun/gamma_q.hpp>
10 #include <stan/math/prim/fun/is_any_nan.hpp>
11 #include <stan/math/prim/fun/is_inf.hpp>
12 #include <stan/math/prim/fun/log.hpp>
13 #include <stan/math/prim/fun/multiply_log.hpp>
14 #include <cmath>
15 #include <limits>
16
17 namespace stan {
18 namespace math {
19
20 /**
21 * Gradient of the regularized incomplete gamma functions igamma(a, z)
22 *
23 * For small z, the gradient is computed via the series expansion;
24 * for large z, the series is numerically inaccurate due to cancellation
25 * and the asymptotic expansion is used.
26 *
27 * @tparam T1 type of the shape parameter
28 * @tparam T2 type of the location parameter
29 * @param a shape parameter, a > 0
30 * @param z location z >= 0
31 * @param g stan::math::tgamma(a) (precomputed value)
32 * @param dig boost::math::digamma(a) (precomputed value)
33 * @param precision required precision; applies to series expansion only
34 * @param max_steps number of steps to take.
35 * @throw throws std::domain_error if not converged after max_steps
36 * or increment overflows to inf.
37 *
38 * For the asymptotic expansion, the gradient is given by:
39 \f[
40 \begin{array}{rcl}
41 \Gamma(a, z) & = & z^{a-1}e^{-z} \sum_{k=0}^N \frac{(a-1)_k}{z^k} \qquad , z
42 \gg a\\
43 Q(a, z) & = & \frac{z^{a-1}e^{-z}}{\Gamma(a)} \sum_{k=0}^N
44 \frac{(a-1)_k}{z^k}\\
45 (a)_k & = & (a)_{k-1}(a-k)\\
46 \frac{d}{da} (a)_k & = & (a)_{k-1} + (a-k)\frac{d}{da} (a)_{k-1}\\
47 \frac{d}{da}Q(a, z) & = & (log(z) - \psi(a)) Q(a, z)\\
48 && + \frac{z^{a-1}e^{-z}}{\Gamma(a)} \sum_{k=0}^N \left(\frac{d}{da}
49 (a-1)_k\right) \frac{1}{z^k} \end{array} \f]
50 */
51 template <typename T1, typename T2>
grad_reg_inc_gamma(T1 a,T2 z,T1 g,T1 dig,double precision=1e-6,int max_steps=1e5)52 return_type_t<T1, T2> grad_reg_inc_gamma(T1 a, T2 z, T1 g, T1 dig,
53 double precision = 1e-6,
54 int max_steps = 1e5) {
55 using std::exp;
56 using std::fabs;
57 using std::log;
58 using TP = return_type_t<T1, T2>;
59
60 if (is_any_nan(a, z, g, dig)) {
61 return std::numeric_limits<TP>::quiet_NaN();
62 }
63
64 T2 l = log(z);
65 if (z >= a && z >= 8) {
66 // asymptotic expansion http://dlmf.nist.gov/8.11#E2
67 TP S = 0;
68 T1 a_minus_one_minus_k = a - 1;
69 T1 fac = a_minus_one_minus_k; // falling_factorial(a-1, k)
70 T1 dfac = 1; // d/da[falling_factorial(a-1, k)]
71 T2 zpow = z; // z ** k
72 TP delta = dfac / zpow;
73
74 for (int k = 1; k < 10; ++k) {
75 a_minus_one_minus_k -= 1;
76
77 S += delta;
78
79 zpow *= z;
80 dfac = a_minus_one_minus_k * dfac + fac;
81 fac *= a_minus_one_minus_k;
82 delta = dfac / zpow;
83
84 if (is_inf(delta)) {
85 throw_domain_error("grad_reg_inc_gamma", "is not converging", "", "");
86 }
87 }
88
89 return gamma_q(a, z) * (l - dig) + exp(-z + (a - 1) * l) * S / g;
90 } else {
91 // gradient of series expansion http://dlmf.nist.gov/8.7#E3
92 TP S = 0;
93 TP log_s = 0.0;
94 double s_sign = 1.0;
95 T2 log_z = log(z);
96 TP log_delta = log_s - multiply_log(2, a);
97 for (int k = 1; k <= max_steps; ++k) {
98 S += s_sign >= 0.0 ? exp(log_delta) : -exp(log_delta);
99 log_s += log_z - log(k);
100 s_sign = -s_sign;
101 log_delta = log_s - multiply_log(2, k + a);
102 if (is_inf(log_delta)) {
103 throw_domain_error("grad_reg_inc_gamma", "is not converging", "", "");
104 }
105 if (log_delta <= log(precision)) {
106 return gamma_p(a, z) * (dig - l) + exp(a * l) * S / g;
107 }
108 }
109 throw_domain_error(
110 "grad_reg_inc_gamma", "k (internal counter)", max_steps, "exceeded ",
111 " iterations, gamma function gradient did not converge.");
112 return INFTY;
113 }
114 }
115
116 } // namespace math
117 } // namespace stan
118 #endif
119