Automatic Differentiation
 
Loading...
Searching...
No Matches
grad_reg_inc_gamma.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_FUN_GRAD_REG_INC_GAMMA_HPP
2#define STAN_MATH_PRIM_FUN_GRAD_REG_INC_GAMMA_HPP
3
14#include <cmath>
15#include <limits>
16
17namespace stan {
18namespace math {
19
51template <typename T1, typename T2>
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 }
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
typename return_type< Ts... >::type return_type_t
Convenience type for the return type of the specified template parameters.
fvar< T > multiply_log(const fvar< T > &x1, const fvar< T > &x2)
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.
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 > gamma_p(const fvar< T > &x1, const fvar< T > &x2)
Definition gamma_p.hpp:16
int is_inf(const fvar< T > &x)
Returns 1 if the input's value is infinite and 0 otherwise.
Definition is_inf.hpp:21
bool is_any_nan(const T &x)
Returns true if the input is NaN and false otherwise.
fvar< T > gamma_q(const fvar< T > &x1, const fvar< T > &x2)
Definition gamma_q.hpp:14
static constexpr double INFTY
Positive infinity.
Definition constants.hpp:46
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 ...