1 #ifndef STAN_MATH_PRIM_FUN_INC_BETA_DDA_HPP
2 #define STAN_MATH_PRIM_FUN_INC_BETA_DDA_HPP
3
4 #include <stan/math/prim/meta.hpp>
5 #include <stan/math/prim/err.hpp>
6 #include <stan/math/prim/fun/fabs.hpp>
7 #include <stan/math/prim/fun/inc_beta.hpp>
8 #include <stan/math/prim/fun/inc_beta_ddb.hpp>
9 #include <stan/math/prim/fun/inv.hpp>
10 #include <stan/math/prim/fun/log.hpp>
11 #include <cmath>
12
13 namespace stan {
14 namespace math {
15
16 /**
17 * Returns the partial derivative of the regularized
18 * incomplete beta function, I_{z}(a, b) with respect to a.
19 * The power series used to compute the derivative tends to
20 * converge slowly when a and b are large, especially if z
21 * approaches 1. The implementation will throw an exception
22 * if the series have not converged within 100,000 iterations.
23 * The current implementation has been tested for values
24 * of a and b up to 12500 and z = 0.999.
25 *
26 * @tparam T scalar types of arguments
27 * @param a first argument
28 * @param b second argument
29 * @param z upper bound of the integral
30 * @param digamma_a value of digamma(a)
31 * @param digamma_ab value of digamma(b)
32 * @return partial derivative of the incomplete beta with respect to a
33 *
34 * @pre a >= 0
35 * @pre b >= 0
36 * @pre 0 <= z <= 1
37 */
38 template <typename T>
inc_beta_dda(T a,T b,T z,T digamma_a,T digamma_ab)39 T inc_beta_dda(T a, T b, T z, T digamma_a, T digamma_ab) {
40 using std::fabs;
41 using std::log;
42 using std::pow;
43
44 if (b > a) {
45 if ((0.1 < z && z <= 0.75 && b > 500) || (0.01 < z && z <= 0.1 && b > 2500)
46 || (0.001 < z && z <= 0.01 && b > 1e5)) {
47 return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab);
48 }
49 }
50
51 if (z > 0.75 && a < 500) {
52 return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab);
53 }
54 if (z > 0.9 && a < 2500) {
55 return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab);
56 }
57 if (z > 0.99 && a < 1e5) {
58 return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab);
59 }
60 if (z > 0.999) {
61 return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab);
62 }
63
64 double threshold = 1e-10;
65
66 const T a_plus_b = a + b;
67 const T a_plus_1 = a + 1;
68
69 digamma_a += inv(a); // Need digamma(a + 1), not digamma(a);
70
71 // Common prefactor to regularize numerator and denominator
72 T prefactor = pow(a_plus_1 / a_plus_b, 3);
73
74 T sum_numer = (digamma_ab - digamma_a) * prefactor;
75 T sum_denom = prefactor;
76
77 T summand = prefactor * z * a_plus_b / a_plus_1;
78
79 T k = 1;
80 digamma_ab += inv(a_plus_b);
81 digamma_a += inv(a_plus_1);
82
83 while (fabs(summand) > threshold) {
84 sum_numer += (digamma_ab - digamma_a) * summand;
85 sum_denom += summand;
86
87 summand *= (1 + (a_plus_b) / k) * (1 + k) / (1 + a_plus_1 / k);
88 digamma_ab += inv(a_plus_b + k);
89 digamma_a += inv(a_plus_1 + k);
90 ++k;
91 summand *= z / k;
92
93 if (k > 1e5) {
94 throw_domain_error("inc_beta_dda",
95 "did not converge within 10000 iterations", "", "");
96 }
97 }
98 return inc_beta(a, b, z) * (log(z) + sum_numer / sum_denom);
99 }
100
101 } // namespace math
102 } // namespace stan
103 #endif
104