1 #ifndef STAN_MATH_PRIM_SCAL_FUN_INC_BETA_DDA_HPP 2 #define STAN_MATH_PRIM_SCAL_FUN_INC_BETA_DDA_HPP 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);
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);
52 return -
inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab);
54 double threshold = 1
e-10;
56 const T a_plus_b = a + b;
57 const T a_plus_1 = a + 1;
62 T prefactor =
pow(a_plus_1 / a_plus_b, 3);
64 T sum_numer = (digamma_ab - digamma_a) * prefactor;
65 T sum_denom = prefactor;
67 T summand = prefactor * z * a_plus_b / a_plus_1;
70 digamma_ab +=
inv(a_plus_b);
71 digamma_a +=
inv(a_plus_1);
73 while (
fabs(summand) > threshold) {
74 sum_numer += (digamma_ab - digamma_a) * summand;
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);
84 domain_error(
"inc_beta_dda",
"did not converge within 10000 iterations",
87 return inc_beta(a, b, z) * (
log(z) + sum_numer / sum_denom);
fvar< T > fabs(const fvar< T > &x)
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)
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)
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.
fvar< T > pow(const fvar< T > &x1, const fvar< T > &x2)
fvar< T > inv(const fvar< T > &x)