Automatic Differentiation
 
Loading...
Searching...
No Matches
grad_reg_lower_inc_gamma.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_FUN_LOWER_REG_INC_GAMMA_HPP
2#define STAN_MATH_PRIM_FUN_LOWER_REG_INC_GAMMA_HPP
3
18#include <limits>
19#include <cmath>
20
21namespace stan {
22namespace math {
23
109template <typename T1, typename T2>
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 digammap1 = digamma(a_plus_n);
168 TP sum_b = digammap1 * exp(a * log_z - lgamma_a_plus_1);
169 lgamma_a_plus_n_plus_1 = lgamma_a_plus_1 + log(a_plus_n);
170 while (true) {
171 digammap1 += 1 / a_plus_n;
172 term = exp(a_plus_n * log_z - lgamma_a_plus_n_plus_1) * digammap1;
173 sum_b += term;
174 if (term <= precision) {
175 return emz * (log_z * sum_a - sum_b);
176 }
177 if (n >= max_steps) {
178 throw_domain_error("grad_reg_lower_inc_gamma", "n (internal counter)",
179 max_steps, "exceeded ",
180 " iterations, gamma_p(a,z) gradient (a) "
181 "did not converge.");
182 }
183 ++n;
184 lgamma_a_plus_n_plus_1 += log1p(a_plus_n);
185 ++a_plus_n;
186 }
187}
188
189} // namespace math
190} // namespace stan
191
192#endif
typename return_type< Ts... >::type return_type_t
Convenience type for the return type of the specified template parameters.
double value_of_rec(const fvar< T > &v)
Return the value of the specified variable.
static constexpr double e()
Return the base of the natural logarithm.
Definition constants.hpp:20
fvar< T > log(const fvar< T > &x)
Definition log.hpp:15
void throw_domain_error(const char *function, const char *name, const T &y, const char *msg1, const char *msg2)
Throw a domain error with a consistently formatted message.
fvar< T > sqrt(const fvar< T > &x)
Definition sqrt.hpp:17
fvar< T > log1p(const fvar< T > &x)
Definition log1p.hpp:12
return_type_t< T1, T2 > grad_reg_inc_gamma(T1 a, T2 z, T1 g, T1 dig, double precision=1e-6, int max_steps=1e5)
Gradient of the regularized incomplete gamma functions igamma(a, z)
fvar< T > lgamma(const fvar< T > &x)
Return the natural logarithm of the gamma function applied to the specified argument.
Definition lgamma.hpp:21
fvar< T > tgamma(const fvar< T > &x)
Return the result of applying the gamma function to the specified argument.
Definition tgamma.hpp:21
bool is_any_nan(const T &x)
Returns true if the input is NaN and false otherwise.
return_type_t< T1, T2 > grad_reg_lower_inc_gamma(const T1 &a, const T2 &z, double precision=1e-10, int max_steps=1e5)
Computes the gradient of the lower regularized incomplete gamma function.
void check_positive_finite(const char *function, const char *name, const T_y &y)
Check if y is positive and finite.
fvar< T > digamma(const fvar< T > &x)
Return the derivative of the log gamma function at the specified argument.
Definition digamma.hpp:23
fvar< T > exp(const fvar< T > &x)
Definition exp.hpp:13
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...