Stan Math Library  2.9.0
reverse mode automatic differentiation
inc_beta_dda.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_SCAL_FUN_INC_BETA_DDA_HPP
2 #define STAN_MATH_PRIM_SCAL_FUN_INC_BETA_DDA_HPP
3 
6 #include <cmath>
7 
8 namespace stan {
9  namespace math {
10 
11  template <typename T>
12  T inc_beta_ddb(T a, T b, T z,
13  T digamma_b, T digamma_ab);
14 
37  template <typename T>
38  T inc_beta_dda(T a, T b, T z,
39  T digamma_a, T digamma_ab) {
40  using std::log;
41 
42  if (b > a)
43  if ((0.1 < z && z <= 0.75 && b > 500)
44  || (0.01 < z && z <= 0.1 && b > 2500)
45  || (0.001 < z && z <= 0.01 && b > 1e5))
46  return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab);
47 
48  if (z > 0.75 && a < 500)
49  return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab);
50  if (z > 0.9 && a < 2500)
51  return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab);
52  if (z > 0.99 && a < 1e5)
53  return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab);
54  if (z > 0.999)
55  return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab);
56 
57  double threshold = 1e-10;
58 
59  digamma_a += 1.0 / a; // Need digamma(a + 1), not digamma(a);
60 
61  // Common prefactor to regularize numerator and denomentator
62  T prefactor = (a + 1) / (a + b);
63  prefactor = prefactor * prefactor * prefactor;
64 
65  T sum_numer = (digamma_ab - digamma_a) * prefactor;
66  T sum_denom = prefactor;
67 
68  T summand = prefactor * z * (a + b) / (a + 1);
69 
70  T k = 1;
71  digamma_ab += 1.0 / (a + b);
72  digamma_a += 1.0 / (a + 1);
73 
74  while (fabs(summand) > threshold) {
75  sum_numer += (digamma_ab - digamma_a) * summand;
76  sum_denom += summand;
77 
78  summand *= (1 + (a + b) / k) * (1 + k) / (1 + (a + 1) / k);
79  digamma_ab += 1.0 / (a + b + k);
80  digamma_a += 1.0 / (a + 1 + k);
81  ++k;
82  summand *= z / k;
83 
84  if (k > 1e5)
85  throw std::domain_error("stan::math::inc_beta_dda did "
86  "not converge within 100000 iterations");
87  }
88  return inc_beta(a, b, z) * (log(z) + sum_numer / sum_denom);
89  }
90 
91  } // math
92 } // stan
93 
94 #endif
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:14
T inc_beta_dda(T a, T b, T z, T digamma_a, T digamma_ab)
Returns the partial derivative of the regularized incomplete beta function, I_{z}(a, b) with respect to a.
fvar< T > log(const fvar< T > &x)
Definition: log.hpp:15
T inc_beta_ddb(T a, T b, T z, T digamma_b, T digamma_ab)
Returns the partial derivative of the regularized incomplete beta function, I_{z}(a, b) with respect to b.
fvar< T > inc_beta(const fvar< T > &a, const fvar< T > &b, const fvar< T > &x)
Definition: inc_beta.hpp:20
void domain_error(const char *function, const char *name, const T &y, const char *msg1, const char *msg2)
Throw a domain error with a consistently formatted message.
double e()
Return the base of the natural logarithm.
Definition: constants.hpp:95

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