Stan  1.3
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
ibeta.hpp
Go to the documentation of this file.
1 #ifndef __STAN__AGRAD__REV__IBETA_HPP__
2 #define __STAN__AGRAD__REV__IBETA_HPP__
3 
4 #include <boost/math/special_functions/digamma.hpp>
5 #include <boost/math/special_functions/gamma.hpp>
6 #include <stan/agrad/rev/var.hpp>
7 #include <stan/agrad/rev/op/vvv_vari.hpp>
8 #include <stan/agrad/rev/op/vvd_vari.hpp>
9 #include <stan/agrad/rev/op/vdv_vari.hpp>
10 #include <stan/agrad/rev/op/vdd_vari.hpp>
11 #include <stan/agrad/rev/op/dvv_vari.hpp>
12 #include <stan/agrad/rev/op/dvd_vari.hpp>
13 #include <stan/agrad/rev/op/ddv_vari.hpp>
14 #include <stan/math/functions/ibeta.hpp>
15 
16 namespace stan {
17  namespace agrad {
18 
19  namespace {
25  double ibeta_hypergeometric_helper(double a, double b, double z, double precision=1e-8, double max_steps=1000) {
26  double val=0;
27  double diff=1;
28  double k=0;
29  double a_2 = a*a;
30  double bprod = 1;
31  while (std::abs(diff) > precision && (++k < max_steps) && !std::isnan(diff)) {
32  val += diff;
33  bprod *= b+k-1.0;
34  diff = a_2*std::pow(a+k,-2)*bprod*std::pow(z,k)/boost::math::tgamma(k+1);
35  }
36  return val;
37  }
38 
39  class ibeta_vvv_vari : public op_vvv_vari {
40  public:
41  ibeta_vvv_vari(vari* avi, vari* bvi, vari* xvi) :
42  op_vvv_vari(stan::math::ibeta(avi->val_,bvi->val_,xvi->val_),avi,bvi,xvi) {
43  }
44  void chain() {
45  double a = avi_->val_;
46  double b = bvi_->val_;
47  double c = cvi_->val_;
48 
49  using std::sin;
50  using std::pow;
51  using std::log;
52  using boost::math::constants::pi;
53  using boost::math::tgamma;
54  using boost::math::digamma;
55  using boost::math::ibeta;
56  using stan::agrad::ibeta_hypergeometric_helper;
57  avi_->adj_ += adj_ *
58  (log(c) - digamma(a) + digamma(a+b)) * val_ -
59  tgamma(a)*tgamma(a+b)/tgamma(b) * pow(c,a) / tgamma(1+a) / tgamma(1+a) * ibeta_hypergeometric_helper(a, 1-b, c);
60  bvi_->adj_ += adj_ *
61  (tgamma(b)*tgamma(a+b)/tgamma(a)*pow(1-c,b) * ibeta_hypergeometric_helper(b,1-a,1-c)/tgamma(b+1)/tgamma(b+1)
62  + ibeta(b, a, 1-c) * (digamma(b) - digamma(a+b) - log(1-c)));
63  cvi_->adj_ += adj_ *
64  boost::math::ibeta_derivative(a,b,c);
65  }
66  };
67  class ibeta_vvd_vari : public op_vvd_vari {
68  public:
69  ibeta_vvd_vari(vari* avi, vari* bvi, double x) :
70  op_vvd_vari(stan::math::ibeta(avi->val_,bvi->val_,x),avi,bvi,x) {
71  }
72  void chain() {
73  double a = avi_->val_;
74  double b = bvi_->val_;
75  double c = cd_;
76 
77  using std::sin;
78  using std::pow;
79  using std::log;
80  using boost::math::constants::pi;
81  using boost::math::tgamma;
82  using boost::math::digamma;
83  using boost::math::ibeta;
84  using stan::agrad::ibeta_hypergeometric_helper;
85  avi_->adj_ += adj_ *
86  (log(c) - digamma(a) + digamma(a+b)) * val_ -
87  tgamma(a)*tgamma(a+b)/tgamma(b) * pow(c,a) / tgamma(1+a) / tgamma(1+a) * ibeta_hypergeometric_helper(a, 1-b, c);
88  bvi_->adj_ += adj_ *
89  (tgamma(b)*tgamma(a+b)/tgamma(a)*pow(1-c,b) * ibeta_hypergeometric_helper(b,1-a,1-c)/tgamma(b+1)/tgamma(b+1)
90  + ibeta(b, a, 1-c) * (digamma(b) - digamma(a+b) - log(1-c)));
91  }
92  };
93  class ibeta_vdv_vari : public op_vdv_vari {
94  public:
95  ibeta_vdv_vari(vari* avi, double b, vari* xvi) :
96  op_vdv_vari(stan::math::ibeta(avi->val_,b,xvi->val_),avi,b,xvi) {
97  }
98  void chain() {
99  double a = avi_->val_;
100  double b = bd_;
101  double c = cvi_->val_;
102 
103  using std::sin;
104  using std::pow;
105  using std::log;
106  using boost::math::constants::pi;
107  using boost::math::tgamma;
108  using boost::math::digamma;
109  using boost::math::ibeta;
110  using stan::agrad::ibeta_hypergeometric_helper;
111  avi_->adj_ += adj_ *
112  (log(c) - digamma(a) + digamma(a+b)) * val_ -
113  tgamma(a)*tgamma(a+b)/tgamma(b) * pow(c,a) / tgamma(1+a) / tgamma(1+a) * ibeta_hypergeometric_helper(a, 1-b, c);
114  cvi_->adj_ += adj_ *
115  boost::math::ibeta_derivative(a,b,c);
116  }
117  };
118  class ibeta_vdd_vari : public op_vdd_vari {
119  public:
120  ibeta_vdd_vari(vari* avi, double b, double x) :
121  op_vdd_vari(stan::math::ibeta(avi->val_,b,x),avi,b,x) {
122  }
123  void chain() {
124  double a = avi_->val_;
125  double b = bd_;
126  double c = cd_;
127 
128  using std::sin;
129  using std::pow;
130  using std::log;
131  using boost::math::constants::pi;
132  using boost::math::tgamma;
133  using boost::math::digamma;
134  using boost::math::ibeta;
135  using stan::agrad::ibeta_hypergeometric_helper;
136  avi_->adj_ += adj_ *
137  (log(c) - digamma(a) + digamma(a+b)) * val_ -
138  tgamma(a)*tgamma(a+b)/tgamma(b) * pow(c,a) / tgamma(1+a) / tgamma(1+a) * ibeta_hypergeometric_helper(a, 1-b, c);
139  }
140  };
141  class ibeta_dvv_vari : public op_dvv_vari {
142  public:
143  ibeta_dvv_vari(double a, vari* bvi, vari* xvi) :
144  op_dvv_vari(stan::math::ibeta(a,bvi->val_,xvi->val_),a,bvi,xvi) {
145  }
146  void chain() {
147  double a = ad_;
148  double b = bvi_->val_;
149  double c = cvi_->val_;
150 
151  using std::sin;
152  using std::pow;
153  using std::log;
154  using boost::math::constants::pi;
155  using boost::math::tgamma;
156  using boost::math::digamma;
157  using boost::math::ibeta;
158  using stan::agrad::ibeta_hypergeometric_helper;
159  bvi_->adj_ += adj_ *
160  (tgamma(b)*tgamma(a+b)/tgamma(a)*pow(1-c,b) * ibeta_hypergeometric_helper(b,1-a,1-c)/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(stan::math::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;
179  using boost::math::constants::pi;
180  using boost::math::tgamma;
181  using boost::math::digamma;
182  using boost::math::ibeta;
183  using stan::agrad::ibeta_hypergeometric_helper;
184  bvi_->adj_ += adj_ *
185  (tgamma(b)*tgamma(a+b)/tgamma(a)*pow(1-c,b) * ibeta_hypergeometric_helper(b,1-a,1-c)/tgamma(b+1)/tgamma(b+1)
186  + ibeta(b, a, 1-c) * (digamma(b) - digamma(a+b) - log(1-c)));
187  }
188  };
189  class ibeta_ddv_vari : public op_ddv_vari {
190  public:
191  ibeta_ddv_vari(double a, double b, vari* xvi) :
192  op_ddv_vari(stan::math::ibeta(a,b,xvi->val_),a,b,xvi) {
193  }
194  void chain() {
195  double a = ad_;
196  double b = bd_;
197  double c = cvi_->val_;
198 
199  cvi_->adj_ += adj_ *
200  boost::math::ibeta_derivative(a,b,c);
201  }
202  };
203  }
204 
222  inline var ibeta(const var& a,
223  const var& b,
224  const var& x) {
225  return var(new ibeta_vvv_vari(a.vi_, b.vi_, x.vi_));
226  }
227 
228  }
229 }
230 #endif
     [ Stan Home Page ] © 2011–2013, Stan Development Team.