Stan Math Library  2.20.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 
8 #include <cmath>
9 
10 namespace stan {
11 namespace math {
12 
35 template <typename T>
36 T inc_beta_dda(T a, T b, T z, T digamma_a, T digamma_ab) {
37  using std::log;
38  using std::pow;
39 
40  if (b > a)
41  if ((0.1 < z && z <= 0.75 && b > 500) || (0.01 < z && z <= 0.1 && b > 2500)
42  || (0.001 < z && z <= 0.01 && b > 1e5))
43  return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab);
44 
45  if (z > 0.75 && a < 500)
46  return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab);
47  if (z > 0.9 && a < 2500)
48  return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab);
49  if (z > 0.99 && a < 1e5)
50  return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab);
51  if (z > 0.999)
52  return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab);
53 
54  double threshold = 1e-10;
55 
56  const T a_plus_b = a + b;
57  const T a_plus_1 = a + 1;
58 
59  digamma_a += inv(a); // Need digamma(a + 1), not digamma(a);
60 
61  // Common prefactor to regularize numerator and denomentator
62  T prefactor = pow(a_plus_1 / a_plus_b, 3);
63 
64  T sum_numer = (digamma_ab - digamma_a) * prefactor;
65  T sum_denom = prefactor;
66 
67  T summand = prefactor * z * a_plus_b / a_plus_1;
68 
69  T k = 1;
70  digamma_ab += inv(a_plus_b);
71  digamma_a += inv(a_plus_1);
72 
73  while (fabs(summand) > threshold) {
74  sum_numer += (digamma_ab - digamma_a) * summand;
75  sum_denom += summand;
76 
77  summand *= (1 + (a_plus_b) / k) * (1 + k) / (1 + a_plus_1 / k);
78  digamma_ab += inv(a_plus_b + k);
79  digamma_a += inv(a_plus_1 + k);
80  ++k;
81  summand *= z / k;
82 
83  if (k > 1e5)
84  domain_error("inc_beta_dda", "did not converge within 10000 iterations",
85  "", "");
86  }
87  return inc_beta(a, b, z) * (log(z) + sum_numer / sum_denom);
88 }
89 
90 } // namespace math
91 } // namespace stan
92 #endif
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
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:12
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:18
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:87
fvar< T > pow(const fvar< T > &x1, const fvar< T > &x2)
Definition: pow.hpp:16
fvar< T > inv(const fvar< T > &x)
Definition: inv.hpp:12

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