Stan Math Library  2.15.0
reverse mode automatic differentiation
ibeta.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_REV_SCAL_FUN_IBETA_HPP
2 #define STAN_MATH_REV_SCAL_FUN_IBETA_HPP
3 
4 #include <stan/math/rev/core.hpp>
8 
9 namespace stan {
10  namespace math {
11 
12  namespace {
18  double ibeta_hypergeometric_helper(double a, double b, double z,
19  double precision = 1e-8,
20  double max_steps = 1000) {
21  double val = 0;
22  double diff = 1;
23  double k = 0;
24  double a_2 = a * a;
25  double bprod = 1;
26  while (std::abs(diff) > precision && ++k < max_steps) {
27  val += diff;
28  bprod *= b + k - 1.0;
29  diff = a_2 * std::pow(a + k, -2) * bprod * std::pow(z, k)
30  / tgamma(k + 1);
31  }
32  return val;
33  }
34 
35  class ibeta_vvv_vari : public op_vvv_vari {
36  public:
37  ibeta_vvv_vari(vari* avi, vari* bvi, vari* xvi) :
38  op_vvv_vari(ibeta(avi->val_, bvi->val_, xvi->val_),
39  avi, bvi, xvi) {
40  }
41  void chain() {
42  double a = avi_->val_;
43  double b = bvi_->val_;
44  double c = cvi_->val_;
45 
46  using std::sin;
47  using std::pow;
48  using std::log;
50  avi_->adj_ += adj_ *
51  (log(c) - digamma(a) + digamma(a+b)) * val_
52  - tgamma(a) * tgamma(a+b) / tgamma(b) * pow(c, a)
53  / tgamma(1+a) / tgamma(1+a)
54  * ibeta_hypergeometric_helper(a, 1-b, c);
55  bvi_->adj_ += adj_ *
56  (tgamma(b) * tgamma(a+b) / tgamma(a) * pow(1-c, b)
57  * ibeta_hypergeometric_helper(b, 1-a, 1-c)
58  / tgamma(b+1) / tgamma(b+1)
59  + ibeta(b, a, 1-c) * (digamma(b) - digamma(a+b) - log(1-c)));
60  cvi_->adj_ += adj_ *
61  boost::math::ibeta_derivative(a, b, c);
62  }
63  };
64  class ibeta_vvd_vari : public op_vvd_vari {
65  public:
66  ibeta_vvd_vari(vari* avi, vari* bvi, double x) :
67  op_vvd_vari(ibeta(avi->val_, bvi->val_, x), avi, bvi, x) {
68  }
69  void chain() {
70  double a = avi_->val_;
71  double b = bvi_->val_;
72  double c = cd_;
73 
74  using std::sin;
75  using std::pow;
76  using std::log;
78  avi_->adj_ += adj_ *
79  (log(c) - digamma(a) + digamma(a+b)) * val_ -
80  tgamma(a) * tgamma(a+b) / tgamma(b) * pow(c, a)
81  / tgamma(1+a) / tgamma(1+a)
82  * ibeta_hypergeometric_helper(a, 1-b, c);
83  bvi_->adj_ += adj_ *
84  (tgamma(b) * tgamma(a+b) / tgamma(a) * pow(1-c, b)
85  * ibeta_hypergeometric_helper(b, 1-a, 1-c)
86  / tgamma(b+1) / tgamma(b+1)
87  + ibeta(b, a, 1-c) * (digamma(b) - digamma(a+b) - log(1-c)));
88  }
89  };
90  class ibeta_vdv_vari : public op_vdv_vari {
91  public:
92  ibeta_vdv_vari(vari* avi, double b, vari* xvi) :
93  op_vdv_vari(ibeta(avi->val_, b, xvi->val_), avi, b, xvi) {
94  }
95  void chain() {
96  double a = avi_->val_;
97  double b = bd_;
98  double c = cvi_->val_;
99 
100  using std::sin;
101  using std::pow;
102  using std::log;
104  using boost::math::tgamma;
105  using boost::math::digamma;
106  using boost::math::ibeta;
107  avi_->adj_ += adj_ *
108  (log(c) - digamma(a) + digamma(a+b)) * val_
109  - tgamma(a) * tgamma(a+b) / tgamma(b) * pow(c, a)
110  / tgamma(1+a) / tgamma(1+a)
111  * ibeta_hypergeometric_helper(a, 1-b, c);
112  cvi_->adj_ += adj_ *
113  boost::math::ibeta_derivative(a, b, c);
114  }
115  };
116  class ibeta_vdd_vari : public op_vdd_vari {
117  public:
118  ibeta_vdd_vari(vari* avi, double b, double x) :
119  op_vdd_vari(ibeta(avi->val_, b, x), avi, b, x) {
120  }
121  void chain() {
122  double a = avi_->val_;
123  double b = bd_;
124  double c = cd_;
125 
126  using std::sin;
127  using std::pow;
128  using std::log;
130  using boost::math::tgamma;
131  using boost::math::digamma;
132  using boost::math::ibeta;
133  avi_->adj_ += adj_ *
134  (log(c) - digamma(a) + digamma(a+b)) * val_
135  - tgamma(a) * tgamma(a+b) / tgamma(b) * pow(c, a)
136  / tgamma(1+a) / tgamma(1+a)
137  * ibeta_hypergeometric_helper(a, 1-b, c);
138  }
139  };
140  class ibeta_dvv_vari : public op_dvv_vari {
141  public:
142  ibeta_dvv_vari(double a, vari* bvi, vari* xvi) :
143  op_dvv_vari(ibeta(a, bvi->val_, xvi->val_), a, bvi, xvi) {
144  }
145  void chain() {
146  double a = ad_;
147  double b = bvi_->val_;
148  double c = cvi_->val_;
149 
150  using std::sin;
151  using std::pow;
152  using std::log;
154  using boost::math::tgamma;
155  using boost::math::digamma;
156  using boost::math::ibeta;
157  bvi_->adj_ += adj_ *
158  (tgamma(b) * tgamma(a+b) / tgamma(a) * pow(1-c, b)
159  * ibeta_hypergeometric_helper(b, 1-a, 1-c)
160  / tgamma(b+1) / tgamma(b+1)
161  + ibeta(b, a, 1-c) * (digamma(b) - digamma(a+b) - log(1-c)));
162  cvi_->adj_ += adj_ *
163  boost::math::ibeta_derivative(a, b, c);
164  }
165  };
166  class ibeta_dvd_vari : public op_dvd_vari {
167  public:
168  ibeta_dvd_vari(double a, vari* bvi, double x) :
169  op_dvd_vari(ibeta(a, bvi->val_, x), a, bvi, x) {
170  }
171  void chain() {
172  double a = ad_;
173  double b = bvi_->val_;
174  double c = cd_;
175 
176  using std::sin;
177  using std::pow;
178  using std::log;
180  using boost::math::tgamma;
181  using boost::math::digamma;
182  using boost::math::ibeta;
183  bvi_->adj_ += adj_ *
184  (tgamma(b) * tgamma(a+b) / tgamma(a) * pow(1-c, b)
185  * ibeta_hypergeometric_helper(b, 1-a, 1-c)
186  / tgamma(b+1) / tgamma(b+1)
187  + ibeta(b, a, 1-c) * (digamma(b) - digamma(a+b) - log(1-c)));
188  }
189  };
190  class ibeta_ddv_vari : public op_ddv_vari {
191  public:
192  ibeta_ddv_vari(double a, double b, vari* xvi) :
193  op_ddv_vari(ibeta(a, b, xvi->val_), a, b, xvi) {
194  }
195  void chain() {
196  double a = ad_;
197  double b = bd_;
198  double c = cvi_->val_;
199 
200  cvi_->adj_ += adj_ *
201  boost::math::ibeta_derivative(a, b, c);
202  }
203  };
204  }
205 
224  inline var ibeta(const var& a,
225  const var& b,
226  const var& x) {
227  return var(new ibeta_vvv_vari(a.vi_, b.vi_, x.vi_));
228  }
229 
230  }
231 }
232 #endif
fvar< T > abs(const fvar< T > &x)
Definition: abs.hpp:15
fvar< T > log(const fvar< T > &x)
Definition: log.hpp:14
Independent (input) and dependent (output) variables for gradients.
Definition: var.hpp:30
fvar< T > sin(const fvar< T > &x)
Definition: sin.hpp:12
vari * vi_
Pointer to the implementation of this variable.
Definition: var.hpp:42
double ibeta(double a, double b, double x)
The normalized incomplete beta function of a, b, and x.
Definition: ibeta.hpp:23
double e()
Return the base of the natural logarithm.
Definition: constants.hpp:94
var ibeta(const var &a, const var &b, const var &x)
The normalized incomplete beta function of a, b, and x.
Definition: ibeta.hpp:224
double pi()
Return the value of pi.
Definition: constants.hpp:85
double adj_
The adjoint of this variable, which is the partial derivative of this variable with respect to the ro...
Definition: vari.hpp:44
fvar< T > pow(const fvar< T > &x1, const fvar< T > &x2)
Definition: pow.hpp:17
fvar< T > tgamma(const fvar< T > &x)
Return the result of applying the gamma function to the specified argument.
Definition: tgamma.hpp:20
fvar< T > digamma(const fvar< T > &x)
Return the derivative of the log gamma function at the specified argument.
Definition: digamma.hpp:22

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