1 #ifndef STAN_MATH_PRIM_FUN_LOWER_REG_INC_GAMMA_HPP
2 #define STAN_MATH_PRIM_FUN_LOWER_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/digamma.hpp>
7 #include <stan/math/prim/fun/exp.hpp>
8 #include <stan/math/prim/fun/gamma_p.hpp>
9 #include <stan/math/prim/fun/grad_reg_inc_gamma.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/lgamma.hpp>
13 #include <stan/math/prim/fun/log.hpp>
14 #include <stan/math/prim/fun/log1p.hpp>
15 #include <stan/math/prim/fun/sqrt.hpp>
16 #include <stan/math/prim/fun/tgamma.hpp>
17 #include <stan/math/prim/fun/value_of_rec.hpp>
18 #include <limits>
19 #include <cmath>
20
21 namespace stan {
22 namespace math {
23
24 /**
25 * Computes the gradient of the lower regularized incomplete
26 * gamma function.
27 *
28 * The lower incomplete gamma function
29 * derivative w.r.t its first parameter (a) seems to have no
30 * standard source. It also appears to have no widely known
31 * approximate implementation. Gautschi (1979) has a thorough
32 * discussion of the calculation of the lower regularized
33 * incomplete gamma function itself and some stability issues.
34 *
35 * Reference: Gautschi, Walter (1979) ACM Transactions on
36 * mathematical software. 5(4):466-481
37 *
38 * We implemented calculations for d(gamma_p)/da by taking
39 * derivatives of formulas suggested by Gauschi and others and
40 * testing them against an outside source (Mathematica). We
41 * took three implementations which can cover the range {a:[0,20],
42 * z:[0,30]} with absolute error < 1e-10 with the exception of
43 * values near (0,0) where the error is near 1e-5. Relative error
44 * is also <<1e-6 except for regions where the gradient approaches
45 * zero.
46 *
47 * Gautschi suggests calculating the lower incomplete gamma
48 * function for small to moderate values of $z$ using the
49 * approximation:
50 *
51 * \f[
52 * \frac{\gamma(a,z)}{\Gamma(a)}=z^a e^-z
53 * \sum_n=0^\infty \frac{z^n}{\Gamma(a+n+1)}
54 * \f]
55 *
56 * We write the derivative in the form:
57 *
58 * \f[
59 * \frac{d\gamma(a,z)\Gamma(a)}{da} = \frac{\log z}{e^z}
60 * \sum_n=0^\infty \frac{z^{a+n}}{\Gamma(a+n+1)}
61 * - \frac{1}{e^z}
62 * \sum_n=0^\infty \frac{z^{a+n}}{\Gamma(a+n+1)}\psi^0(a+n+1)
63 * \f]
64 *
65 * This calculation is sufficiently accurate for small $a$ and
66 * small $z$. For larger values and $a$ and $z$ we use it in its
67 * log form:
68 *
69 * \f[
70 * \frac{d \gamma(a,z)\Gamma(a)}{da} = \frac{\log z}{e^z}
71 * \sum_n=0^\infty \exp[(a+n)\log z - \log\Gamma(a+n+1)]
72 * - \sum_n=0^\infty \exp[(a+n)\log z - \log\Gamma(a+n+1) +
73 * \log\psi^0(a+n+1)]
74 * \f]
75 *
76 * For large $z$, Gauschi recommends using the upper incomplete
77 * Gamma instead and the negative of its derivative turns out to be
78 * more stable and accurate for larger $z$ and for some combinations
79 * of $a$ and $z$. This is a log-scale implementation of the
80 * derivative of the formulation suggested by Gauschi (1979). For
81 * some values it defers to the negative of the gradient
82 * for the gamma_q function. This is based on the suggestion by Gauschi
83 * (1979) that for large values of $z$ it is better to
84 * carry out calculations using the upper incomplete Gamma function.
85 *
86 * Branching for choice of implementation for the lower incomplete
87 * regularized gamma function gradient. The derivative based on
88 * Gautschi's formulation appears to be sufficiently accurate
89 * everywhere except for large z and small to moderate a. The
90 * intersection between the two regions is a radius 12 quarter circle
91 * centered at a=0, z=30 although both implementations are
92 * satisfactory near the intersection.
93 *
94 * Some limits that could be treated, e.g., infinite z should
95 * return tgamma(a) * digamma(a), throw instead to match the behavior of,
96 * e.g., boost::math::gamma_p
97 *
98 * @tparam T1 type of a
99 * @tparam T2 type of z
100 * @param[in] a shared with complete Gamma
101 * @param[in] z value to integrate up to
102 * @param[in] precision series terminates when increment falls below
103 * this value.
104 * @param[in] max_steps number of terms to sum before throwing
105 * @throw std::domain_error if the series does not converge to
106 * requested precision before max_steps.
107 *
108 */
109 template <typename T1, typename T2>
grad_reg_lower_inc_gamma(const T1 & a,const T2 & z,double precision=1e-10,int max_steps=1e5)110 return_type_t<T1, T2> grad_reg_lower_inc_gamma(const T1& a, const T2& z,
111 double precision = 1e-10,
112 int max_steps = 1e5) {
113 using std::exp;
114 using std::log;
115 using std::pow;
116 using std::sqrt;
117
118 using TP = return_type_t<T1, T2>;
119
120 if (is_any_nan(a, z)) {
121 return std::numeric_limits<TP>::quiet_NaN();
122 }
123
124 check_positive_finite("grad_reg_lower_inc_gamma", "a", a);
125
126 if (z == 0.0) {
127 return 0.0;
128 }
129 check_positive_finite("grad_reg_lower_inc_gamma", "z", z);
130
131 if ((a < 0.8 && z > 15.0) || (a < 12.0 && z > 30.0)
132 || a < sqrt(-756 - value_of_rec(z) * value_of_rec(z)
133 + 60 * value_of_rec(z))) {
134 T1 tg = tgamma(a);
135 T1 dig = digamma(a);
136 return -grad_reg_inc_gamma(a, z, tg, dig, max_steps, precision);
137 }
138
139 T2 log_z = log(z);
140 T2 emz = exp(-z);
141
142 int n = 0;
143 T1 a_plus_n = a;
144 TP sum_a = 0.0;
145 T1 lgamma_a_plus_1 = lgamma(a + 1);
146 T1 lgamma_a_plus_n_plus_1 = lgamma_a_plus_1;
147 TP term;
148 while (true) {
149 term = exp(a_plus_n * log_z - lgamma_a_plus_n_plus_1);
150 sum_a += term;
151 if (term <= precision) {
152 break;
153 }
154 if (n >= max_steps) {
155 throw_domain_error("grad_reg_lower_inc_gamma", "n (internal counter)",
156 max_steps, "exceeded ",
157 " iterations, gamma_p(a,z) gradient (a) "
158 "did not converge.");
159 }
160 ++n;
161 lgamma_a_plus_n_plus_1 += log1p(a_plus_n);
162 ++a_plus_n;
163 }
164
165 n = 1;
166 a_plus_n = a + 1;
167 TP sum_b = digamma(a + 1) * exp(a * log_z - lgamma_a_plus_1);
168 lgamma_a_plus_n_plus_1 = lgamma_a_plus_1 + log(a_plus_n);
169 while (true) {
170 term = exp(a_plus_n * log_z - lgamma_a_plus_n_plus_1)
171 * digamma(a_plus_n + 1);
172 sum_b += term;
173 if (term <= precision) {
174 return emz * (log_z * sum_a - sum_b);
175 }
176 if (n >= max_steps) {
177 throw_domain_error("grad_reg_lower_inc_gamma", "n (internal counter)",
178 max_steps, "exceeded ",
179 " iterations, gamma_p(a,z) gradient (a) "
180 "did not converge.");
181 }
182 ++n;
183 lgamma_a_plus_n_plus_1 += log1p(a_plus_n);
184 ++a_plus_n;
185 }
186 }
187
188 } // namespace math
189 } // namespace stan
190
191 #endif
192