diff --git a/src/stan/math/prim/scal/fun.hpp b/src/stan/math/prim/scal/fun.hpp index a190a3ea541..240c07320e6 100644 --- a/src/stan/math/prim/scal/fun.hpp +++ b/src/stan/math/prim/scal/fun.hpp @@ -28,6 +28,9 @@ #include #include #include +#include +#include +#include #include #include #include diff --git a/src/stan/math/prim/scal/fun/inc_beta_dda.hpp b/src/stan/math/prim/scal/fun/inc_beta_dda.hpp new file mode 100644 index 00000000000..7891feea2ff --- /dev/null +++ b/src/stan/math/prim/scal/fun/inc_beta_dda.hpp @@ -0,0 +1,94 @@ +#ifndef STAN_MATH_PRIM_SCAL_FUN_INC_BETA_DDA_HPP +#define STAN_MATH_PRIM_SCAL_FUN_INC_BETA_DDA_HPP + +#include +#include +#include + +namespace stan { + namespace math { + + template + 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 a. + * The power series used to compute the deriative tends to + * converge slowly when a and b are large, especially if z + * approaches 1. The implementation will throw an exception + * if the series have not converged within 100,000 iterations. + * The current implementation has been tested for values + * of a and b up to 12500 and z = 0.999. + * + * @tparam T scalar types of arguments + * @param a a + * @param b b + * @param z upper bound of the integral + * @param digamma_a value of digamma(a) + * @param digamma_ab value of digamma(b) + * @return partial derivative of the incomplete beta with respect to a + * + * @pre a >= 0 + * @pre b >= 0 + * @pre 0 <= z <= 1 + */ + template + T inc_beta_dda(T a, T b, T z, + T digamma_a, T digamma_ab) { + using std::log; + + if (b > a) + if ((0.1 < z && z <= 0.75 && b > 500) + || (0.01 < z && z <= 0.1 && b > 2500) + || (0.001 < z && z <= 0.01 && b > 1e5)) + return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab); + + if (z > 0.75 && a < 500) + return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab); + if (z > 0.9 && a < 2500) + return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab); + if (z > 0.99 && a < 1e5) + return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab); + if (z > 0.999) + return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab); + + double threshold = 1e-10; + + digamma_a += 1.0 / a; // Need digamma(a + 1), not digamma(a); + + // Common prefactor to regularize numerator and denomentator + T prefactor = (a + 1) / (a + b); + prefactor = prefactor * prefactor * prefactor; + + T sum_numer = (digamma_ab - digamma_a) * prefactor; + T sum_denom = prefactor; + + T summand = prefactor * z * (a + b) / (a + 1); + + T k = 1; + digamma_ab += 1.0 / (a + b); + digamma_a += 1.0 / (a + 1); + + while (fabs(summand) > threshold) { + sum_numer += (digamma_ab - digamma_a) * summand; + sum_denom += summand; + + summand *= (1 + (a + b) / k) * (1 + k) / (1 + (a + 1) / k); + digamma_ab += 1.0 / (a + b + k); + digamma_a += 1.0 / (a + 1 + k); + ++k; + summand *= z / k; + + if (k > 1e5) + throw std::domain_error("stan::math::inc_beta_dda did " + "not converge within 100000 iterations"); + } + return inc_beta(a, b, z) * (log(z) + sum_numer / sum_denom); + } + + } // math +} // stan + +#endif diff --git a/src/stan/math/prim/scal/fun/inc_beta_ddb.hpp b/src/stan/math/prim/scal/fun/inc_beta_ddb.hpp new file mode 100644 index 00000000000..52af39d89e1 --- /dev/null +++ b/src/stan/math/prim/scal/fun/inc_beta_ddb.hpp @@ -0,0 +1,89 @@ +#ifndef STAN_MATH_PRIM_SCAL_FUN_INC_BETA_DDB_HPP +#define STAN_MATH_PRIM_SCAL_FUN_INC_BETA_DDB_HPP + +#include +#include +#include + +namespace stan { + namespace math { + + template + 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 b. + * The power series used to compute the deriative tends to + * converge slowly when a and b are large, especailly if z + * approaches 1. The implementation will throw an exception + * if the series have not converged within 100,000 iterations. + * The current implementation has been tested for values + * of a and b up to 12500 and z = 0.999. + * + * @tparam T scalar types of arguments + * @param a a + * @param b b + * @param z upper bound of the integral + * @param digamma_a value of digamma(a) + * @param digamma_ab value of digamma(b) + * @return partial derivative of the incomplete beta with respect to b + * + * @pre a >= 0 + * @pre b >= 0 + * @pre 0 <= z <= 1 + */ + template + T inc_beta_ddb(T a, T b, T z, + T digamma_b, T digamma_ab) { + using std::log; + + if (b > a) + if ((0.1 < z && z <= 0.75 && b > 500) + || (0.01 < z && z <= 0.1 && b > 2500) + || (0.001 < z && z <= 0.01 && b > 1e5)) + return -inc_beta_dda(b, a, 1 - z, digamma_b, digamma_ab); + + if ((z > 0.75 && a < 500) + || (z > 0.9 && a < 2500) + || (z > 0.99 && a < 1e5) + || (z > 0.999)) + return -inc_beta_dda(b, a, 1 - z, digamma_b, digamma_ab); + + double threshold = 1e-10; + + // Common prefactor to regularize numerator and denomentator + T prefactor = (a + 1) / (a + b); + prefactor = prefactor * prefactor * prefactor; + + T sum_numer = digamma_ab * prefactor; + T sum_denom = prefactor; + + T summand = prefactor * z * (a + b) / (a + 1); + + T k = 1; + digamma_ab += 1.0 / (a + b); + + while (fabs(summand) > threshold) { + sum_numer += digamma_ab * summand; + sum_denom += summand; + + summand *= (1 + (a + b) / k) * (1 + k) / (1 + (a + 1) / k); + digamma_ab += 1.0 / (a + b + k); + ++k; + summand *= z / k; + + if (k > 1e5) + throw std::domain_error("stan::math::inc_beta_ddb did " + "not converge within 100000 iterations"); + } + + return inc_beta(a, b, z) + * (log(1 - z) - digamma_b + sum_numer / sum_denom); + } + + } // math +} // stan + +#endif diff --git a/src/stan/math/prim/scal/fun/inc_beta_ddz.hpp b/src/stan/math/prim/scal/fun/inc_beta_ddz.hpp new file mode 100644 index 00000000000..d3edb7635d3 --- /dev/null +++ b/src/stan/math/prim/scal/fun/inc_beta_ddz.hpp @@ -0,0 +1,44 @@ +#ifndef STAN_MATH_PRIM_SCAL_FUN_INC_BETA_DERIVATIVES_HPP +#define STAN_MATH_PRIM_SCAL_FUN_INC_BETA_DERIVATIVES_HPP + +#include +#include +#include +#include + +namespace stan { + namespace math { + + + /** + * Returns the partial derivative of the regularized + * incomplete beta function, I_{z}(a, b) with respect to z. + * + * @tparam T scalar types of arguments + * @param a a + * @param b b + * @param z upper bound of the integral + * @return partial derivative of the incomplete beta with respect to z + * + * @pre a > 0 + * @pre b > 0 + * @pre 0 < z <= 1 + */ + template + T inc_beta_ddz(T a, T b, T z) { + using std::exp; + using std::log; + return exp((b - 1) * log(1 - z) + (a - 1) * log(z) + + lgamma(a + b) - lgamma(a) - lgamma(b)); + } + + template <> + double inc_beta_ddz(double a, double b, double z) { + using boost::math::ibeta_derivative; + return ibeta_derivative(a, b, z); + } + + } // math +} // stan + +#endif diff --git a/src/stan/math/prim/scal/prob/neg_binomial_2_cdf.hpp b/src/stan/math/prim/scal/prob/neg_binomial_2_cdf.hpp index fccedaf622d..5c45bb7e5fb 100644 --- a/src/stan/math/prim/scal/prob/neg_binomial_2_cdf.hpp +++ b/src/stan/math/prim/scal/prob/neg_binomial_2_cdf.hpp @@ -2,27 +2,26 @@ #define STAN_MATH_PRIM_SCAL_PROB_NEG_BINOMIAL_2_CDF_HPP #include +#include +#include +#include + #include #include #include -#include -#include -#include -#include -#include -#include +#include + #include -#include -#include -#include -#include -#include -#include -#include +#include +#include +#include +#include +#include #include -namespace stan { +#include +namespace stan { namespace math { template ::type + T_partials_return; + using stan::math::check_positive_finite; + using stan::math::check_nonnegative; using stan::math::check_not_nan; using stan::math::check_consistent_sizes; - using stan::math::check_less; - static const char* function("stan::math::neg_binomial_2_cdf"); + // Ensure non-zero arugment lengths + if (!(stan::length(n) && stan::length(mu) && stan::length(phi))) + return 1.0; + + T_partials_return P(1.0); + + // Validate arguments check_positive_finite(function, "Location parameter", mu); check_positive_finite(function, "Precision parameter", phi); check_not_nan(function, "Random variable", n); @@ -51,38 +55,94 @@ namespace stan { "Location parameter", mu, "Precision Parameter", phi); + // Wrap arguments in vector views VectorView n_vec(n); VectorView mu_vec(mu); VectorView phi_vec(phi); + size_t size = max_size(n, mu, phi); + + // Compute vectorized CDF and gradient + using stan::math::value_of; + using stan::math::inc_beta; + using stan::math::inc_beta_ddz; + using stan::math::inc_beta_dda; + using stan::math::digamma; + + OperandsAndPartials + operands_and_partials(mu, phi); - size_t size_phi_mu = max_size(mu, phi); - size_t size_n = length(n); - - std::vector::type> - phi_mu(size_phi_mu); - std::vector::type> np1(size_n); - - for (size_t i = 0; i < size_phi_mu; i++) - phi_mu[i] = phi_vec[i]/(phi_vec[i]+mu_vec[i]); - - for (size_t i = 0; i < size_n; i++) - if (n_vec[i] < 0) - return 0.0; - else - np1[i] = n_vec[i] + 1.0; - - if (size_n == 1) { - if (size_phi_mu == 1) - return beta_cdf(phi_mu[0], phi, np1[0]); - else - return beta_cdf(phi_mu, phi, np1[0]); - } else { - if (size_phi_mu == 1) - return beta_cdf(phi_mu[0], phi, np1); - else - return beta_cdf(phi_mu, phi, np1); + // Explicit return for extreme values + // The gradients are technically ill-defined, but treated as zero + for (size_t i = 0; i < stan::length(n); i++) { + if (value_of(n_vec[i]) < 0) + return operands_and_partials.to_var(0.0, mu, phi); } + + // Cache a few expensive function calls if phi is a parameter + VectorBuilder::value, + T_partials_return, T_precision> + digamma_phi_vec(stan::length(phi)); + + VectorBuilder::value, + T_partials_return, T_precision> + digamma_sum_vec(stan::length(phi)); + + if (!is_constant_struct::value) { + for (size_t i = 0; i < stan::length(phi); i++) { + const T_partials_return n_dbl = value_of(n_vec[i]); + const T_partials_return phi_dbl = value_of(phi_vec[i]); + + digamma_phi_vec[i] = digamma(phi_dbl); + digamma_sum_vec[i] = digamma(n_dbl + phi_dbl + 1); + } + } + + for (size_t i = 0; i < size; i++) { + // Explicit results for extreme values + // The gradients are technically ill-defined, but treated as zero + if (value_of(n_vec[i]) == std::numeric_limits::max()) + return operands_and_partials.to_var(1.0, mu, phi); + + const T_partials_return n_dbl = value_of(n_vec[i]); + const T_partials_return mu_dbl = value_of(mu_vec[i]); + const T_partials_return phi_dbl = value_of(phi_vec[i]); + + const T_partials_return p_dbl = phi_dbl / (mu_dbl + phi_dbl); + const T_partials_return d_dbl = 1.0 / ((mu_dbl + phi_dbl) + * (mu_dbl + phi_dbl)); + + const T_partials_return P_i = + inc_beta(phi_dbl, n_dbl + 1.0, p_dbl); + + P *= P_i; + + if (!is_constant_struct::value) + operands_and_partials.d_x1[i] += + - inc_beta_ddz(phi_dbl, n_dbl + 1.0, p_dbl) * phi_dbl * d_dbl / P_i; + + if (!is_constant_struct::value) { + operands_and_partials.d_x2[i] + += inc_beta_dda(phi_dbl, n_dbl + 1, p_dbl, + digamma_phi_vec[i], + digamma_sum_vec[i]) / P_i + + inc_beta_ddz(phi_dbl, n_dbl + 1.0, p_dbl) + * mu_dbl * d_dbl / P_i; + } + } + + if (!is_constant_struct::value) { + for (size_t i = 0; i < stan::length(mu); ++i) + operands_and_partials.d_x1[i] *= P; + } + + if (!is_constant_struct::value) { + for (size_t i = 0; i < stan::length(phi); ++i) + operands_and_partials.d_x2[i] *= P; + } + + return operands_and_partials.to_var(P, mu, phi); } - } -} + + } // math +} // stan #endif diff --git a/src/stan/math/prim/scal/prob/neg_binomial_2_log_log.hpp b/src/stan/math/prim/scal/prob/neg_binomial_2_log_log.hpp index 1f537ddc72f..52d1196ea28 100644 --- a/src/stan/math/prim/scal/prob/neg_binomial_2_log_log.hpp +++ b/src/stan/math/prim/scal/prob/neg_binomial_2_log_log.hpp @@ -37,7 +37,7 @@ namespace stan { T_precision>::type T_partials_return; - static const char* function("stan::math::neg_binomial_log"); + static const char* function("stan::prob::neg_binomial_2_log_log"); using stan::math::check_finite; using stan::math::check_nonnegative; diff --git a/src/stan/math/prim/scal/prob/neg_binomial_cdf.hpp b/src/stan/math/prim/scal/prob/neg_binomial_cdf.hpp index 25b28df1735..625660f3583 100644 --- a/src/stan/math/prim/scal/prob/neg_binomial_cdf.hpp +++ b/src/stan/math/prim/scal/prob/neg_binomial_cdf.hpp @@ -1,29 +1,25 @@ #ifndef STAN_MATH_PRIM_SCAL_PROB_NEG_BINOMIAL_CDF_HPP #define STAN_MATH_PRIM_SCAL_PROB_NEG_BINOMIAL_CDF_HPP -#include -#include -#include #include +#include +#include +#include + #include -#include #include -#include +#include + #include -#include -#include #include -#include -#include -#include -#include #include -#include +#include +#include +#include #include #include namespace stan { - namespace math { // Negative Binomial CDF @@ -38,9 +34,7 @@ namespace stan { T_partials_return; using stan::math::check_positive_finite; - using stan::math::check_nonnegative; using stan::math::check_consistent_sizes; - using stan::math::include_summand; // Ensure non-zero arugment lengths if (!(stan::length(n) && stan::length(alpha) && stan::length(beta))) @@ -65,11 +59,9 @@ namespace stan { // Compute vectorized CDF and gradient using stan::math::value_of; using stan::math::inc_beta; + using stan::math::inc_beta_ddz; + using stan::math::inc_beta_dda; using stan::math::digamma; - using stan::math::lbeta; - using std::exp; - using std::pow; - using std::exp; OperandsAndPartials operands_and_partials(alpha, beta); @@ -84,29 +76,26 @@ namespace stan { // Cache a few expensive function calls if alpha is a parameter VectorBuilder::value, T_partials_return, T_shape> - digammaN_vec(stan::length(alpha)); - VectorBuilder::value, - T_partials_return, T_shape> - digammaAlpha_vec(stan::length(alpha)); + digamma_alpha_vec(stan::length(alpha)); + VectorBuilder::value, T_partials_return, T_shape> - digammaSum_vec(stan::length(alpha)); + digamma_sum_vec(stan::length(alpha)); + if (!is_constant_struct::value) { for (size_t i = 0; i < stan::length(alpha); i++) { const T_partials_return n_dbl = value_of(n_vec[i]); const T_partials_return alpha_dbl = value_of(alpha_vec[i]); - digammaN_vec[i] = digamma(n_dbl + 1); - digammaAlpha_vec[i] = digamma(alpha_dbl); - digammaSum_vec[i] = digamma(n_dbl + alpha_dbl + 1); + digamma_alpha_vec[i] = digamma(alpha_dbl); + digamma_sum_vec[i] = digamma(n_dbl + alpha_dbl + 1); } } for (size_t i = 0; i < size; i++) { // Explicit results for extreme values // The gradients are technically ill-defined, but treated as zero - if (value_of(n_vec[i]) - == std::numeric_limits::max()) + if (value_of(n_vec[i]) == std::numeric_limits::max()) return operands_and_partials.to_var(1.0, alpha, beta); const T_partials_return n_dbl = value_of(n_vec[i]); @@ -117,31 +106,21 @@ namespace stan { const T_partials_return d_dbl = 1.0 / ( (1.0 + beta_dbl) * (1.0 + beta_dbl) ); - const T_partials_return Pi = inc_beta(alpha_dbl, n_dbl + 1.0, p_dbl); - - const T_partials_return beta_func = exp(lbeta(n_dbl + 1, alpha_dbl)); + const T_partials_return P_i = + inc_beta(alpha_dbl, n_dbl + 1.0, p_dbl); - - P *= Pi; + P *= P_i; if (!is_constant_struct::value) { - T_partials_return g1 = 0; - T_partials_return g2 = 0; - - stan::math::grad_reg_inc_beta(g1, g2, alpha_dbl, - n_dbl + 1, p_dbl, - digammaAlpha_vec[i], - digammaN_vec[i], - digammaSum_vec[i], - beta_func); - operands_and_partials.d_x1[i] - += g1 / Pi; + += inc_beta_dda(alpha_dbl, n_dbl + 1, p_dbl, + digamma_alpha_vec[i], + digamma_sum_vec[i]) / P_i; } if (!is_constant_struct::value) - operands_and_partials.d_x2[i] += d_dbl * pow(1-p_dbl, n_dbl) - * pow(p_dbl, alpha_dbl-1) / beta_func / Pi; + operands_and_partials.d_x2[i] += + inc_beta_ddz(alpha_dbl, n_dbl + 1.0, p_dbl) * d_dbl / P_i; } if (!is_constant_struct::value) { @@ -156,6 +135,7 @@ namespace stan { return operands_and_partials.to_var(P, alpha, beta); } - } -} + + } // prob +} // stan #endif diff --git a/src/test/unit/math/prim/scal/fun/inc_beta_dda_test.cpp b/src/test/unit/math/prim/scal/fun/inc_beta_dda_test.cpp new file mode 100644 index 00000000000..472e75f1de4 --- /dev/null +++ b/src/test/unit/math/prim/scal/fun/inc_beta_dda_test.cpp @@ -0,0 +1,71 @@ +#include +#include +#include + +TEST(MathFunctions, inc_beta_dda) { + using stan::math::digamma; + using stan::math::inc_beta_dda; + + double small_a = 1.5; + double large_a = 15000; + + double small_b = 1.25; + double large_b = 12500; + + double small_z = 0.001; + double mid_z = 0.6; + double large_z = 0.999; + + EXPECT_FLOAT_EQ(-0.00028665637, inc_beta_dda(small_a, small_b, small_z, + digamma(small_a), + digamma(small_a + small_b))) + << "reasonable values for a, b, x"; + EXPECT_FLOAT_EQ(-0.23806757, inc_beta_dda(small_a, small_b, mid_z, + digamma(small_a), + digamma(small_a + small_b))) + << "reasonable values for a, b, x"; + EXPECT_FLOAT_EQ(-0.00022264493, inc_beta_dda(small_a, small_b, large_z, + digamma(small_a), + digamma(small_a + small_b))) + << "reasonable values for a, b, x"; + + + EXPECT_FLOAT_EQ(0.0, inc_beta_dda(large_a, small_b, small_z, + digamma(large_a), + digamma(large_a + small_b))) + << "reasonable values for a, b, x"; + EXPECT_FLOAT_EQ(0.0, inc_beta_dda(large_a, small_b, mid_z, + digamma(large_a), + digamma(large_a + small_b))) + << "reasonable values for a, b, x"; + EXPECT_FLOAT_EQ(0.00091716705, inc_beta_dda(large_a, small_b, large_z, + digamma(large_a), + digamma(large_a + small_b))) + << "reasonable values for a, b, x"; + + EXPECT_FLOAT_EQ(1.8226241, inc_beta_dda(small_a, large_b, small_z, + digamma(small_a), + digamma(small_a + large_b))) + << "reasonable values for a, b, x"; + EXPECT_FLOAT_EQ(0.0, inc_beta_dda(small_a, large_b, mid_z, + digamma(small_a), + digamma(small_a + large_b))) + << "reasonable values for a, b, x"; + EXPECT_FLOAT_EQ(0.0, inc_beta_dda(small_a, large_b, large_z, + digamma(small_a), + digamma(small_a + large_b))) + << "reasonable values for a, b, x"; + + EXPECT_FLOAT_EQ(0.0, inc_beta_dda(large_a, large_b, small_z, + digamma(large_a), + digamma(large_a + large_b))) + << "reasonable values for a, b, x"; + EXPECT_FLOAT_EQ(-4.0856207e-14, inc_beta_dda(large_a, large_b, mid_z, + digamma(large_a), + digamma(large_a + large_b))) + << "reasonable values for a, b, x"; + EXPECT_FLOAT_EQ(0.0, inc_beta_dda(large_a, large_b, large_z, + digamma(large_a), + digamma(large_a + large_b))) + << "reasonable values for a, b, x"; +} diff --git a/src/test/unit/math/prim/scal/fun/inc_beta_ddb_test.cpp b/src/test/unit/math/prim/scal/fun/inc_beta_ddb_test.cpp new file mode 100644 index 00000000000..6b4ed5e23f2 --- /dev/null +++ b/src/test/unit/math/prim/scal/fun/inc_beta_ddb_test.cpp @@ -0,0 +1,72 @@ +#include +#include +#include + +TEST(MathFunctions, inc_beta_ddb) { + using stan::math::digamma; + using stan::math::inc_beta_ddb; + + double small_a = 1.5; + double large_a = 15000; + + double small_b = 1.25; + double large_b = 12500; + + double small_z = 0.001; + double mid_z = 0.5; + double large_z = 0.999; + + EXPECT_FLOAT_EQ(3.2996082e-05, inc_beta_ddb(small_a, small_b, small_z, + digamma(small_a), + digamma(small_a + small_b))) + << "reasonable values for a, b, x"; + EXPECT_FLOAT_EQ(0.17945045, inc_beta_ddb(small_a, small_b, mid_z, + digamma(small_a), + digamma(small_a + small_b))) + << "reasonable values for a, b, x"; + EXPECT_FLOAT_EQ(0.0019721228, inc_beta_ddb(small_a, small_b, large_z, + digamma(small_a), + digamma(small_a + small_b))) + << "reasonable values for a, b, x"; + + EXPECT_FLOAT_EQ(0.0, inc_beta_ddb(large_a, small_b, small_z, + digamma(large_a), + digamma(large_a + small_b))) + << "reasonable values for a, b, x"; + EXPECT_FLOAT_EQ(0.0, inc_beta_ddb(large_a, small_b, mid_z, + digamma(large_a), + digamma(large_a + small_b))) + << "reasonable values for a, b, x"; + EXPECT_FLOAT_EQ(7.7076669, inc_beta_ddb(large_a, small_b, large_z, + digamma(large_a), + digamma(large_a + small_b))) + << "reasonable values for a, b, x"; + + EXPECT_FLOAT_EQ(9.3959293, inc_beta_ddb(small_a, large_b, small_z, + digamma(small_a), + digamma(small_a + large_b))) + << "reasonable values for a, b, x"; + EXPECT_FLOAT_EQ(0.0, inc_beta_ddb(small_a, large_b, mid_z, + digamma(small_a), + digamma(small_a + large_b))) + << "reasonable values for a, b, x"; + EXPECT_FLOAT_EQ(0.0, inc_beta_ddb(small_a, large_b, large_z, + digamma(small_a), + digamma(small_a + large_b))) + << "reasonable values for a, b, x"; + + EXPECT_FLOAT_EQ(0.0, inc_beta_ddb(large_a, large_b, small_z, + digamma(large_a), + digamma(large_a + large_b))) + << "reasonable values for a, b, x"; + EXPECT_FLOAT_EQ(0.0, inc_beta_ddb(large_a, large_b, mid_z, + digamma(large_a), + digamma(large_a + large_b))) + << "reasonable values for a, b, x"; + EXPECT_FLOAT_EQ(0.0, inc_beta_ddb(large_a, large_b, large_z, + digamma(large_a), + digamma(large_a + large_b))) + << "reasonable values for a, b, x"; + +} + diff --git a/src/test/unit/math/prim/scal/fun/inc_beta_ddz_test.cpp b/src/test/unit/math/prim/scal/fun/inc_beta_ddz_test.cpp new file mode 100644 index 00000000000..346247497c6 --- /dev/null +++ b/src/test/unit/math/prim/scal/fun/inc_beta_ddz_test.cpp @@ -0,0 +1,47 @@ +#include +#include +#include + +TEST(MathFunctions, inc_beta_ddz) { + using stan::math::inc_beta_ddz; + + double small_a = 1.5; + double large_a = 15000; + + double small_b = 1.25; + double large_b = 12500; + + double small_z = 0.001; + double mid_z = 0.5; + double large_z = 0.999; + + EXPECT_FLOAT_EQ(0.063300692, inc_beta_ddz(small_a, small_b, small_z)) + << "reasonable values for a, b, x"; + EXPECT_FLOAT_EQ(1.1905416, inc_beta_ddz(small_a, small_b, mid_z)) + << "reasonable values for a, b, x"; + EXPECT_FLOAT_EQ(0.35587692, inc_beta_ddz(small_a, small_b, large_z)) + << "reasonable values for a, b, x"; + + EXPECT_FLOAT_EQ(0.0, inc_beta_ddz(large_a, small_b, small_z)) + << "reasonable values for a, b, x"; + EXPECT_FLOAT_EQ(0.0, inc_beta_ddz(large_a, small_b, mid_z)) + << "reasonable values for a, b, x"; + EXPECT_FLOAT_EQ(0.009898182, inc_beta_ddz(large_a, small_b, large_z)) + << "reasonable values for a, b, x"; + + EXPECT_FLOAT_EQ(0.1848717, inc_beta_ddz(small_a, large_b, small_z)) + << "reasonable values for a, b, x"; + EXPECT_FLOAT_EQ(0.0, inc_beta_ddz(small_a, large_b, mid_z)) + << "reasonable values for a, b, x"; + EXPECT_FLOAT_EQ(0.0, inc_beta_ddz(small_a, large_b, large_z)) + << "reasonable values for a, b, x"; + + EXPECT_FLOAT_EQ(0.0, inc_beta_ddz(large_a, large_b, small_z)) + << "reasonable values for a, b, x"; + EXPECT_FLOAT_EQ(0.0, inc_beta_ddz(large_a, large_b, mid_z)) + << "reasonable values for a, b, x"; + EXPECT_FLOAT_EQ(0.0, inc_beta_ddz(large_a, large_b, large_z)) + << "reasonable values for a, b, x"; + +} +