Stan Math Library  2.20.0
reverse mode automatic differentiation
grad_reg_lower_inc_gamma.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_SCAL_FUN_LOWER_REG_INC_GAMMA_HPP
2 #define STAN_MATH_PRIM_SCAL_FUN_LOWER_REG_INC_GAMMA_HPP
3 
15 #include <limits>
16 
17 namespace stan {
18 namespace math {
19 
105 template <typename T1, typename T2>
107  const T1& a, const T2& z, double precision = 1e-10, int max_steps = 1e5) {
108  using std::exp;
109  using std::log;
110  using std::pow;
111  typedef typename return_type<T1, T2>::type TP;
112 
113  if (is_any_nan(a, z))
114  return std::numeric_limits<TP>::quiet_NaN();
115 
116  check_positive_finite("grad_reg_lower_inc_gamma", "a", a);
117 
118  if (z == 0.0)
119  return 0.0;
120  check_positive_finite("grad_reg_lower_inc_gamma", "z", z);
121 
122  if ((a < 0.8 && z > 15.0) || (a < 12.0 && z > 30.0)
123  || a < sqrt(-756 - value_of_rec(z) * value_of_rec(z)
124  + 60 * value_of_rec(z))) {
125  T1 tg = tgamma(a);
126  T1 dig = digamma(a);
127  return -grad_reg_inc_gamma(a, z, tg, dig, max_steps, precision);
128  }
129 
130  T2 log_z = log(z);
131  T2 emz = exp(-z);
132 
133  int n = 0;
134  T1 a_plus_n = a;
135  TP sum_a = 0.0;
136  T1 lgamma_a_plus_1 = lgamma(a + 1);
137  T1 lgamma_a_plus_n_plus_1 = lgamma_a_plus_1;
138  TP term;
139  while (true) {
140  term = exp(a_plus_n * log_z - lgamma_a_plus_n_plus_1);
141  sum_a += term;
142  if (term <= precision)
143  break;
144  if (n >= max_steps)
145  domain_error("grad_reg_lower_inc_gamma", "n (internal counter)",
146  max_steps, "exceeded ",
147  " iterations, gamma_p(a,z) gradient (a) "
148  "did not converge.");
149  ++n;
150  lgamma_a_plus_n_plus_1 += log1p(a_plus_n);
151  ++a_plus_n;
152  }
153 
154  n = 1;
155  a_plus_n = a + 1;
156  TP sum_b = digamma(a + 1) * exp(a * log_z - lgamma_a_plus_1);
157  lgamma_a_plus_n_plus_1 = lgamma_a_plus_1 + log(a_plus_n);
158  while (true) {
159  term = exp(a_plus_n * log_z - lgamma_a_plus_n_plus_1)
160  * digamma(a_plus_n + 1);
161  sum_b += term;
162  if (term <= precision)
163  return emz * (log_z * sum_a - sum_b);
164  if (n >= max_steps)
165  domain_error("grad_reg_lower_inc_gamma", "n (internal counter)",
166  max_steps, "exceeded ",
167  " iterations, gamma_p(a,z) gradient (a) "
168  "did not converge.");
169  ++n;
170  lgamma_a_plus_n_plus_1 += log1p(a_plus_n);
171  ++a_plus_n;
172  }
173 }
174 
175 } // namespace math
176 } // namespace stan
177 
178 #endif
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 > sqrt(const fvar< T > &x)
Definition: sqrt.hpp:13
double value_of_rec(const fvar< T > &v)
Return the value of the specified variable.
fvar< T > log(const fvar< T > &x)
Definition: log.hpp:12
bool is_any_nan(const T &x)
Returns true if the input is NaN and false otherwise.
Definition: is_any_nan.hpp:21
void check_positive_finite(const char *function, const char *name, const T_y &y)
Check if y is positive and finite.
return_type< T1, T2 >::type 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)
boost::math::tools::promote_args< double, typename scalar_type< T >::type, typename return_type< Types_pack... >::type >::type type
Definition: return_type.hpp:36
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:11
void 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.
double e()
Return the base of the natural logarithm.
Definition: constants.hpp:87
fvar< T > log1p(const fvar< T > &x)
Definition: log1p.hpp:12
fvar< T > pow(const fvar< T > &x1, const fvar< T > &x2)
Definition: pow.hpp:16
fvar< T > tgamma(const fvar< T > &x)
Return the result of applying the gamma function to the specified argument.
Definition: tgamma.hpp:21
return_type< T1, T2 >::type 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.
fvar< T > digamma(const fvar< T > &x)
Return the derivative of the log gamma function at the specified argument.
Definition: digamma.hpp:23

     [ Stan Home Page ] © 2011–2018, Stan Development Team.