Stan Math Library  2.11.0
reverse mode automatic differentiation
gamma_p.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_REV_SCAL_FUN_GAMMA_P_HPP
2 #define STAN_MATH_REV_SCAL_FUN_GAMMA_P_HPP
3 
4 #include <stan/math/rev/core.hpp>
6 #include <boost/math/special_functions/gamma.hpp>
7 #include <boost/math/special_functions/digamma.hpp>
8 #include <valarray>
9 
10 namespace stan {
11  namespace math {
12 
13  namespace {
14  class gamma_p_vv_vari : public op_vv_vari {
15  public:
16  gamma_p_vv_vari(vari* avi, vari* bvi) :
17  op_vv_vari(stan::math::gamma_p(avi->val_, bvi->val_),
18  avi, bvi) {
19  }
20  void chain() {
21  // return zero derivative as gamma_p is flat
22  // to machine precision for b / a > 10
23  if (std::fabs(bvi_->val_ / avi_->val_) > 10 ) return;
24 
25  double u = stan::math::gamma_p(avi_->val_, bvi_->val_);
26 
27  double S = 0.0;
28  double s = 1.0;
29  double l = std::log(bvi_->val_);
30  double g = boost::math::tgamma(avi_->val_);
31  double dig = boost::math::digamma(avi_->val_);
32 
33  int k = 0;
34  double delta = s / (avi_->val_ * avi_->val_);
35 
36  while (std::fabs(delta) > 1e-6) {
37  S += delta;
38  ++k;
39  s *= -bvi_->val_ / k;
40  delta = s / ((k + avi_->val_) * (k + avi_->val_));
41  }
42 
43 
44  avi_->adj_ -= adj_ * ((u) * (dig - l)
45  + std::exp(avi_->val_ * l) * S / g);
46  bvi_->adj_ += adj_ * (std::exp(-bvi_->val_)
47  * std::pow(bvi_->val_, avi_->val_ - 1.0) / g);
48  }
49  };
50 
51  class gamma_p_vd_vari : public op_vd_vari {
52  public:
53  gamma_p_vd_vari(vari* avi, double b) :
54  op_vd_vari(stan::math::gamma_p(avi->val_, b),
55  avi, b) {
56  }
57  void chain() {
58  // return zero derivative as gamma_p is flat
59  // to machine precision for b / a > 10
60  if (std::fabs(bd_ / avi_->val_) > 10)
61  return;
62 
63  double u = stan::math::gamma_p(avi_->val_, bd_);
64 
65  double S = 0.0;
66  double s = 1.0;
67  double l = std::log(bd_);
68  double g = boost::math::tgamma(avi_->val_);
69  double dig = boost::math::digamma(avi_->val_);
70 
71  int k = 0;
72  double delta = s / (avi_->val_ * avi_->val_);
73 
74  while (std::fabs(delta) > 1e-6) {
75  S += delta;
76  ++k;
77  s *= -bd_ / k;
78  delta = s / ((k + avi_->val_) * (k + avi_->val_));
79  }
80 
81  avi_->adj_ -= adj_ * ((u) * (dig - l)
82  + std::exp(avi_->val_ * l) * S / g);
83  }
84  };
85 
86  class gamma_p_dv_vari : public op_dv_vari {
87  public:
88  gamma_p_dv_vari(double a, vari* bvi) :
89  op_dv_vari(stan::math::gamma_p(a, bvi->val_),
90  a, bvi) {
91  }
92  void chain() {
93  // return zero derivative as gamma_p is flat to
94  // machine precision for b / a > 10
95  if (std::fabs(bvi_->val_ / ad_) > 10 )
96  return;
97  bvi_->adj_ += adj_
98  * (std::exp(-bvi_->val_) * std::pow(bvi_->val_, ad_ - 1.0)
99  / boost::math::tgamma(ad_));
100  }
101  };
102  }
103 
104  inline var gamma_p(const stan::math::var& a,
105  const stan::math::var& b) {
106  return var(new gamma_p_vv_vari(a.vi_, b.vi_));
107  }
108 
109  inline var gamma_p(const stan::math::var& a,
110  const double& b) {
111  return var(new gamma_p_vd_vari(a.vi_, b));
112  }
113 
114  inline var gamma_p(const double& a,
115  const stan::math::var& b) {
116  return var(new gamma_p_dv_vari(a, b.vi_));
117  }
118 
119  }
120 }
121 #endif
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:14
fvar< T > log(const fvar< T > &x)
Definition: log.hpp:15
Independent (input) and dependent (output) variables for gradients.
Definition: var.hpp:31
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
vari * vi_
Pointer to the implementation of this variable.
Definition: var.hpp:43
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.