| Stan Math Library 2.11.0 reverse mode automatic differentiation

- Main Page
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

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>
5 #include <stan/math/prim/scal/fun/gamma_p.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
core.hpp
stan::math::fabs
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:14
stan
Definition: log_sum_exp.hpp:8
stan::math::log
fvar< T > log(const fvar< T > &x)
Definition: log.hpp:15
stan::math::var
Independent (input) and dependent (output) variables for gradients.
Definition: var.hpp:31
stan::math::exp
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
stan::math::var::vi_
vari * vi_
Pointer to the implementation of this variable.
Definition: var.hpp:43
stan::math::gamma_p
fvar< T > gamma_p(const fvar< T > &x1, const fvar< T > &x2)
Definition: gamma_p.hpp:15
gamma_p.hpp
stan::math::e
double e()
Return the base of the natural logarithm.
Definition: constants.hpp:95
stan::math::pow
fvar< T > pow(const fvar< T > &x1, const fvar< T > &x2)
Definition: pow.hpp:18
stan::math::tgamma
fvar< T > tgamma(const fvar< T > &x)
Definition: tgamma.hpp:15
stan::math::digamma
fvar< T > digamma(const fvar< T > &x)
Definition: digamma.hpp:16

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