Stan Math Library  2.10.0
reverse mode automatic differentiation
gamma_p.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_FWD_SCAL_FUN_GAMMA_P_HPP
2 #define STAN_MATH_FWD_SCAL_FUN_GAMMA_P_HPP
3 
4 #include <stan/math/fwd/core.hpp>
5 
7 
8 namespace stan {
9 
10  namespace math {
11 
12  template <typename T>
13  inline
14  fvar<T>
15  gamma_p(const fvar<T>& x1, const fvar<T>& x2) {
16  using stan::math::gamma_p;
17  using std::log;
18  using std::exp;
19  using std::pow;
20  using std::fabs;
21  using boost::math::tgamma;
23 
24  T u = gamma_p(x1.val_, x2.val_);
25 
26  T S = 0;
27  T s = 1;
28  T l = log(x2.val_);
29  T g = tgamma(x1.val_);
30  T dig = digamma(x1.val_);
31 
32  int k = 0;
33  T delta = s / (x1.val_ * x1.val_);
34 
35  while (fabs(delta) > 1e-6) {
36  S += delta;
37  ++k;
38  s *= -x2.val_ / k;
39  delta = s / ((k + x1.val_) * (k + x1.val_));
40  }
41 
42  T der1 = u * (dig - l) + exp(x1.val_ * l) * S / g;
43  T der2 = exp(-x2.val_) * pow(x2.val_, x1.val_ - 1.0) / g;
44 
45  return fvar<T>(u, x1.d_ * -der1 + x2.d_ * der2);
46  }
47 
48  template <typename T>
49  inline
50  fvar<T>
51  gamma_p(const fvar<T>& x1, const double x2) {
52  using stan::math::gamma_p;
53  using std::log;
54  using std::exp;
55  using std::pow;
56  using std::fabs;
57  using boost::math::tgamma;
59 
60  T u = gamma_p(x1.val_, x2);
61 
62  T S = 0.0;
63  double s = 1.0;
64  double l = log(x2);
65  T g = tgamma(x1.val_);
66  T dig = digamma(x1.val_);
67 
68  int k = 0;
69  T delta = s / (x1.val_ * x1.val_);
70 
71  while (fabs(delta) > 1e-6) {
72  S += delta;
73  ++k;
74  s *= -x2 / k;
75  delta = s / ((k + x1.val_) * (k + x1.val_));
76  }
77 
78  T der1 = u * (dig - l) + exp(x1.val_ * l) * S / g;
79 
80  return fvar<T>(u, x1.d_ * -der1);
81  }
82 
83  template <typename T>
84  inline
85  fvar<T>
86  gamma_p(const double x1, const fvar<T>& x2) {
87  using stan::math::gamma_p;
88  using std::exp;
89  using std::pow;
90 
91  T u = gamma_p(x1, x2.val_);
92 
93  double g = boost::math::tgamma(x1);
94 
95  T der2 = exp(-x2.val_) * pow(x2.val_, x1 - 1.0) / g;
96 
97  return fvar<T>(u, x2.d_ * der2);
98  }
99  }
100 }
101 #endif
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:14
fvar< T > log(const fvar< T > &x)
Definition: log.hpp:15
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
fvar< T > gamma_p(const fvar< T > &x1, const fvar< T > &x2)
Definition: gamma_p.hpp:15
double e()
Return the base of the natural logarithm.
Definition: constants.hpp:95
fvar< T > pow(const fvar< T > &x1, const fvar< T > &x2)
Definition: pow.hpp:18
fvar< T > tgamma(const fvar< T > &x)
Definition: tgamma.hpp:15
fvar< T > digamma(const fvar< T > &x)
Definition: digamma.hpp:16

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