1 #ifndef STAN_MATH_REV_FUN_GRAD_INC_BETA_HPP
2 #define STAN_MATH_REV_FUN_GRAD_INC_BETA_HPP
3
4 #include <stan/math/rev/meta.hpp>
5 #include <stan/math/rev/core.hpp>
6 #include <stan/math/rev/fun/beta.hpp>
7 #include <stan/math/rev/fun/exp.hpp>
8 #include <stan/math/rev/fun/fabs.hpp>
9 #include <stan/math/rev/fun/floor.hpp>
10 #include <stan/math/rev/fun/value_of_rec.hpp>
11 #include <stan/math/rev/fun/inc_beta.hpp>
12 #include <stan/math/rev/fun/lgamma.hpp>
13 #include <stan/math/rev/fun/log.hpp>
14 #include <stan/math/rev/fun/log1m.hpp>
15 #include <stan/math/rev/fun/value_of.hpp>
16 #include <stan/math/prim/fun/grad_2F1.hpp>
17 #include <stan/math/prim/fun/value_of.hpp>
18 #include <cmath>
19
20 namespace stan {
21 namespace math {
22
23 /**
24 * Gradient of the incomplete beta function beta(a, b, z) with
25 * respect to the first two arguments.
26 *
27 * Uses the equivalence to a hypergeometric function. See
28 * http://dlmf.nist.gov/8.17#ii
29 *
30 * @param[out] g1 d/da
31 * @param[out] g2 d/db
32 * @param[in] a a
33 * @param[in] b b
34 * @param[in] z z
35 */
grad_inc_beta(var & g1,var & g2,const var & a,const var & b,const var & z)36 inline void grad_inc_beta(var& g1, var& g2, const var& a, const var& b,
37 const var& z) {
38 var c1 = log(z);
39 var c2 = log1m(z);
40 var c3 = beta(a, b) * inc_beta(a, b, z);
41 var C = exp(a * c1 + b * c2) / a;
42 var dF1 = 0;
43 var dF2 = 0;
44 if (value_of(value_of(C))) {
45 grad_2F1(dF1, dF2, a + b, var(1.0), a + 1, z);
46 }
47 g1 = (c1 - 1.0 / a) * c3 + C * (dF1 + dF2);
48 g2 = c2 * c3 + C * dF1;
49 }
50
51 } // namespace math
52 } // namespace stan
53 #endif
54