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