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