From f3b127638d8e0a0f242c08e467d5dfa84b4b4746 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Thu, 30 Dec 2021 20:51:16 +0800 Subject: [PATCH 01/12] Initial implementation, start tests --- stan/math/fwd/fun.hpp | 1 + stan/math/fwd/fun/inc_beta_inv.hpp | 123 +++++++++++++++++++ stan/math/prim/fun.hpp | 1 + stan/math/prim/fun/F32.hpp | 17 +-- stan/math/prim/fun/inc_beta_inv.hpp | 33 +++++ stan/math/rev/fun.hpp | 1 + stan/math/rev/fun/inc_beta_inv.hpp | 92 ++++++++++++++ test/unit/math/fwd/fun/inc_beta_inv_test.cpp | 32 +++++ test/unit/math/mix/fun/inc_beta_inv_test.cpp | 79 ++++++++++++ test/unit/math/rev/fun/inc_beta_inv_test.cpp | 21 ++++ 10 files changed, 392 insertions(+), 8 deletions(-) create mode 100644 stan/math/fwd/fun/inc_beta_inv.hpp create mode 100644 stan/math/prim/fun/inc_beta_inv.hpp create mode 100644 stan/math/rev/fun/inc_beta_inv.hpp create mode 100644 test/unit/math/fwd/fun/inc_beta_inv_test.cpp create mode 100644 test/unit/math/mix/fun/inc_beta_inv_test.cpp create mode 100644 test/unit/math/rev/fun/inc_beta_inv_test.cpp diff --git a/stan/math/fwd/fun.hpp b/stan/math/fwd/fun.hpp index d33dc0fa4a7..0c2d4559867 100644 --- a/stan/math/fwd/fun.hpp +++ b/stan/math/fwd/fun.hpp @@ -43,6 +43,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/fwd/fun/inc_beta_inv.hpp b/stan/math/fwd/fun/inc_beta_inv.hpp new file mode 100644 index 00000000000..274ce485b87 --- /dev/null +++ b/stan/math/fwd/fun/inc_beta_inv.hpp @@ -0,0 +1,123 @@ +#ifndef STAN_MATH_FWD_FUN_INC_BETA_INV_HPP +#define STAN_MATH_FWD_FUN_INC_BETA_INV_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * The fused multiply-add operation (C99). + * + * This double-based operation delegates to fma. + * + * The function is defined by + * + * fma(a, b, c) = (a * b) + c. + * + * + \f[ + \mbox{fma}(x, y, z) = + \begin{cases} + x\cdot y+z & \mbox{if } -\infty\leq x, y, z \leq \infty \\[6pt] + \textrm{NaN} & \mbox{if } x = \textrm{NaN} + \end{cases} + \f] + + \f[ + \frac{\partial\, \mbox{fma}(x, y, z)}{\partial x} = + \begin{cases} + y & \mbox{if } -\infty\leq x, y, z \leq \infty \\[6pt] + \textrm{NaN} & \mbox{if } x = \textrm{NaN} + \end{cases} + \f] + + \f[ + \frac{\partial\, \mbox{fma}(x, y, z)}{\partial y} = + \begin{cases} + x & \mbox{if } -\infty\leq x, y, z \leq \infty \\[6pt] + \textrm{NaN} & \mbox{if } x = \textrm{NaN} + \end{cases} + \f] + + \f[ + \frac{\partial\, \mbox{fma}(x, y, z)}{\partial z} = + \begin{cases} + 1 & \mbox{if } -\infty\leq x, y, z \leq \infty \\[6pt] + \textrm{NaN} & \mbox{if } x = \textrm{NaN} + \end{cases} + \f] + * + * @param x1 First value. + * @param x2 Second value. + * @param x3 Third value. + * @return Product of the first two values plus the third. + */ +template * = nullptr, + require_any_fvar_t* = nullptr> +inline fvar> inc_beta_inv(const T1& a, + const T2& b, + const T3& p) { + using T_return = partials_return_t; + auto a_val = value_of(a); + auto b_val = value_of(b); + auto p_val = value_of(p); + T_return w = inc_beta_inv(a_val, b_val, p_val); + T_return log_w = log(w); + T_return log1m_w = log1m(w); + auto one_m_a = 1 - a_val; + auto one_m_b = 1 - b_val; + T_return one_m_w = 1 - w; + auto ap1 = a_val + 1; + auto bp1 = b_val + 1; + auto lbeta_ab = lbeta(a_val, b_val); + auto digamma_apb = digamma(a_val + b_val); + + T_return inv_d_(0); + + if(is_fvar::value) { + auto da1 = one_m_b * log1m_w + one_m_a * log_w; + auto da2 = a_val * log_w + 2 * lgamma(a_val) + + log(F32(a_val, a_val, one_m_b, ap1, ap1, w)) + - 2 * lgamma(ap1); + auto da3 = lbeta_ab + log(inc_beta(a_val, b_val, w)) + + log(log_w - digamma(a_val) + digamma_apb); + inv_d_ += forward_as>(a).d_ * exp(da1 + log_diff_exp(da2, da3)); + } + + if(is_fvar::value) { + auto db1 = (w - 1) * exp(-b_val * log1m_w + one_m_a * log_w); + auto db2 = 2 * lgamma(b_val) + + log(F32(b_val, b_val, one_m_a, bp1, bp1, one_m_w)) + - 2 * lgamma(bp1) + + b_val * log1m_w; + + auto db3 = log(inc_beta(b_val, a_val, one_m_w)) + lbeta_ab + + log(log1m_w - digamma(b_val) + digamma_apb); + + inv_d_ += forward_as>(b).d_ * db1 * exp(log_diff_exp(db2, db3)); + } + + if(is_fvar::value) { + inv_d_ += forward_as>(p).d_ * + exp(one_m_b * log1m_w + one_m_a * log_w + lbeta_ab); + } + + return fvar(w, inv_d_); +} + + +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp index f3e9f3bc21c..6c181e3ef6c 100644 --- a/stan/math/prim/fun.hpp +++ b/stan/math/prim/fun.hpp @@ -134,6 +134,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/fun/F32.hpp b/stan/math/prim/fun/F32.hpp index c76eeef69fc..ecccd9bcaa1 100644 --- a/stan/math/prim/fun/F32.hpp +++ b/stan/math/prim/fun/F32.hpp @@ -49,29 +49,30 @@ namespace math { * @param[in] precision precision of the infinite sum. defaults to 1e-6 * @param[in] max_steps number of steps to take. defaults to 1e5 */ -template -T F32(const T& a1, const T& a2, const T& a3, const T& b1, const T& b2, - const T& z, double precision = 1e-6, int max_steps = 1e5) { +template +return_type_t F32(const Ta1& a1, const Ta2& a2, const Ta3& a3, const Tb1& b1, const Tb2& b2, + const Tz& z, double precision = 1e-6, int max_steps = 1e5) { check_3F2_converges("F32", a1, a2, a3, b1, b2, z); + using T_return = return_type_t; using std::exp; using std::fabs; using std::log; - T t_acc = 1.0; - T log_t = 0.0; - T log_z = log(z); + T_return t_acc = 1.0; + T_return log_t = 0.0; + Tz log_z = log(z); double t_sign = 1.0; for (int k = 0; k <= max_steps; ++k) { - T p = (a1 + k) * (a2 + k) * (a3 + k) / ((b1 + k) * (b2 + k) * (k + 1)); + T_return p = (a1 + k) * (a2 + k) * (a3 + k) / ((b1 + k) * (b2 + k) * (k + 1)); if (p == 0.0) { return t_acc; } log_t += log(fabs(p)) + log_z; t_sign = p >= 0.0 ? t_sign : -t_sign; - T t_new = t_sign > 0.0 ? exp(log_t) : -exp(log_t); + T_return t_new = t_sign > 0.0 ? exp(log_t) : -exp(log_t); t_acc += t_new; if (fabs(t_new) <= precision) { diff --git a/stan/math/prim/fun/inc_beta_inv.hpp b/stan/math/prim/fun/inc_beta_inv.hpp new file mode 100644 index 00000000000..0ecf33cf291 --- /dev/null +++ b/stan/math/prim/fun/inc_beta_inv.hpp @@ -0,0 +1,33 @@ +#ifndef STAN_MATH_PRIM_FUN_INC_BETA_INV_HPP +#define STAN_MATH_PRIM_FUN_INC_BETA_INV_HPP + +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * The inverse of the normalized incomplete beta function of a, b, with probability p. + * + * Used to compute the cumulative density function for the beta + * distribution. + * + * @param a Shape parameter a >= 0; a and b can't both be 0 + * @param b Shape parameter b >= 0 + * @param p Random variate. 0 <= p <= 1 + * @throws if constraints are violated or if any argument is NaN + * @return The inverse of the normalized incomplete beta function. + */ +inline double inc_beta_inv(double a, double b, double p) { + check_not_nan("inc_beta", "a", a); + check_not_nan("inc_beta", "b", b); + check_not_nan("inc_beta", "p", p); + return boost::math::ibeta_inv(a, b, p, boost_policy_t<>()); +} + +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp index 7545389c9e3..a21bb201691 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -78,6 +78,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/rev/fun/inc_beta_inv.hpp b/stan/math/rev/fun/inc_beta_inv.hpp new file mode 100644 index 00000000000..b804543602f --- /dev/null +++ b/stan/math/rev/fun/inc_beta_inv.hpp @@ -0,0 +1,92 @@ +#ifndef STAN_MATH_REV_FUN_INC_BETA_INV_HPP +#define STAN_MATH_REV_FUN_INC_BETA_INV_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * The fused multiply-add function for three variables (C99). + * This function returns the product of the first two arguments + * plus the third argument. + * + * The partial derivatives are + * + * \f$\frac{\partial}{\partial x} (x * y) + z = y\f$, and + * + * \f$\frac{\partial}{\partial y} (x * y) + z = x\f$, and + * + * \f$\frac{\partial}{\partial z} (x * y) + z = 1\f$. + * + * @param x First multiplicand. + * @param y Second multiplicand. + * @param z Summand. + * @return Product of the multiplicands plus the summand, ($a * $b) + $c. + */ +template * = nullptr, + require_any_var_t* = nullptr> +inline var inc_beta_inv(const T1& a, const T2& b, const T3& p) { + double a_val = value_of(a); + double b_val = value_of(b); + double p_val = value_of(p); + double w = inc_beta_inv(a_val, b_val, p_val); + return make_callback_var(w, [a, b, p, a_val, b_val, p_val, w](auto& vi) { + double log_w = log(w); + double log1m_w = log1m(w); + double one_m_a = 1 - a_val; + double one_m_b = 1 - b_val; + double one_m_w = 1 - w; + double ap1 = a_val + 1; + double bp1 = b_val + 1; + double lbeta_ab = lbeta(a_val, b_val); + double digamma_apb = digamma(a_val + b_val); + + if(!is_constant_all::value) { + double da1 = one_m_b * log1m_w + one_m_a * log_w; + double da2 = a_val * log_w + 2 * lgamma(a_val) + + log(F32(a_val, a_val, one_m_b, ap1, ap1, w)) + - 2 * lgamma(ap1); + double da3 = lbeta_ab + log(inc_beta(a_val, b_val, w)) + + log(log_w - digamma(a_val) + digamma_apb); + + forward_as(a).adj() += vi.adj() * exp(da1 + log_diff_exp(da2, da3)); + } + + if(!is_constant_all::value) { + double db1 = (w - 1) * exp(-b_val * log1m_w + one_m_a * log_w); + double db2 = 2 * lgamma(b_val) + + log(F32(b_val, b_val, one_m_a, bp1, bp1, one_m_w)) + - 2 * lgamma(bp1) + + b_val * log1m_w; + + double db3 = log(inc_beta(b_val, a_val, one_m_w)) + lbeta_ab + + log(log1m_w - digamma(b_val) + digamma_apb); + + forward_as(b).adj() += vi.adj() * db1 * exp(log_diff_exp(db2, db3)); + } + + if(!is_constant_all::value) { + forward_as(p).adj() += vi.adj() * + exp(one_m_b * log1m_w + one_m_a * log_w + lbeta_ab); + } + }); +} + + +} // namespace math +} // namespace stan +#endif diff --git a/test/unit/math/fwd/fun/inc_beta_inv_test.cpp b/test/unit/math/fwd/fun/inc_beta_inv_test.cpp new file mode 100644 index 00000000000..18be7930306 --- /dev/null +++ b/test/unit/math/fwd/fun/inc_beta_inv_test.cpp @@ -0,0 +1,32 @@ +#include +#include + +TEST(AgradFwdMatrixIncBetaInv, fd_scalar) { + using stan::math::fvar; + using stan::math::inc_beta_inv; + fvar a = 1; + fvar b = 2; + fvar p = 0.5; + a.d_ = 1.0; + b.d_ = 1.0; + p.d_ = 1.0; + + fvar res = inc_beta_inv(a, b, p); + + EXPECT_FLOAT_EQ(res.d_, 0.287698278597 - 0.122532267934 + 0.707106781187); +} + +TEST(AgradFwdMatrixIncBetaInv, ffd_scalar) { + using stan::math::fvar; + using stan::math::inc_beta_inv; + fvar> a = 1; + fvar> b = 2; + fvar> p = 0.5; + a.val_.d_ = 1.0; + b.val_.d_ = 1.0; + p.val_.d_ = 1.0; + + fvar> res = inc_beta_inv(a, b, p); + + EXPECT_FLOAT_EQ(res.val_.d_, 0.287698278597 - 0.122532267934 + 0.707106781187); +} \ No newline at end of file diff --git a/test/unit/math/mix/fun/inc_beta_inv_test.cpp b/test/unit/math/mix/fun/inc_beta_inv_test.cpp new file mode 100644 index 00000000000..236fecf664f --- /dev/null +++ b/test/unit/math/mix/fun/inc_beta_inv_test.cpp @@ -0,0 +1,79 @@ +#include +#include +#include + +TEST(ProbInternalMath, inc_beta_inv_fv1) { + using stan::math::fvar; + using stan::math::var; + using stan::math::inc_beta_inv; + double a_d = 1; + double b_d = 2; + double p_d = 0.5; + fvar a_v = a_d; + fvar b_v = b_d; + fvar p_v = p_d; + a_v.d_ = 1.0; + b_v.d_ = 1.0; + p_v.d_ = 1.0; + + fvar res = inc_beta_inv(a_v, b_v, p_v); + res.val_.grad(); + + EXPECT_FLOAT_EQ(a_v.val_.adj(), 0.287698278597); + EXPECT_FLOAT_EQ(b_v.val_.adj(), -0.122532267934); + EXPECT_FLOAT_EQ(p_v.val_.adj(), 0.707106781187); + + a_v = a_d; + b_v = b_d; + p_v = p_d; + a_v.d_ = 1.0; + b_v.d_ = 1.0; + p_v.d_ = 1.0; + + res = inc_beta_inv(a_d, b_v, p_v); + res.val_.grad(); + + EXPECT_FLOAT_EQ(b_v.val_.adj(), -0.122532267934); + EXPECT_FLOAT_EQ(p_v.val_.adj(), 0.707106781187); + + b_v = b_d; + p_v = p_d; + b_v.d_ = 1.0; + p_v.d_ = 1.0; + + res = inc_beta_inv(a_v, b_d, p_v); + res.val_.grad(); + + EXPECT_FLOAT_EQ(a_v.val_.adj(), 0.287698278597); + EXPECT_FLOAT_EQ(p_v.val_.adj(), 0.707106781187); + + a_v = a_d; + p_v = p_d; + a_v.d_ = 1.0; + p_v.d_ = 1.0; + + res = inc_beta_inv(a_v, b_v, p_d); + res.val_.grad(); + + EXPECT_FLOAT_EQ(a_v.val_.adj(), 0.287698278597); + EXPECT_FLOAT_EQ(b_v.val_.adj(), -0.122532267934); +} + +TEST(ProbInternalMath, inc_beta_inv_fv2) { + using stan::math::fvar; + using stan::math::var; + using stan::math::inc_beta_inv; + fvar> a = 1; + fvar> b = 2; + fvar> p = 0.5; + a.d_ = 1.0; + b.d_ = 1.0; + p.d_ = 1.0; + + fvar> res = inc_beta_inv(a, b, p); + res.val_.val_.grad(); + + EXPECT_FLOAT_EQ(a.val_.val_.adj(), 0.287698278597); + EXPECT_FLOAT_EQ(b.val_.val_.adj(), -0.122532267934); + EXPECT_FLOAT_EQ(p.val_.val_.adj(), 0.707106781187); +} \ No newline at end of file diff --git a/test/unit/math/rev/fun/inc_beta_inv_test.cpp b/test/unit/math/rev/fun/inc_beta_inv_test.cpp new file mode 100644 index 00000000000..6ef31ea657d --- /dev/null +++ b/test/unit/math/rev/fun/inc_beta_inv_test.cpp @@ -0,0 +1,21 @@ +#include +#include +#include +#include + +TEST(inc_beta_inv, values) { + using stan::math::var; + using stan::math::inc_beta_inv; + + // check autodiffing works with var types with large values + var a = 1; + var b = 2; + var p = 0.5; + + var res = inc_beta_inv(a, b, p); + res.grad(); + + EXPECT_FLOAT_EQ(a.adj(), 0.287698278597); + EXPECT_FLOAT_EQ(b.adj(), -0.122532267934); + EXPECT_FLOAT_EQ(p.adj(), 0.707106781187); +} From 6aeb64bfe1cf1316249f3ba2c0155a792fb8385f Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 31 Dec 2021 12:45:44 +0800 Subject: [PATCH 02/12] Update tests --- stan/math/fwd/fun/inc_beta_inv.hpp | 19 ++++++++++--------- stan/math/rev/fun/inc_beta_inv.hpp | 14 +++++++------- test/unit/math/fwd/fun/inc_beta_inv_test.cpp | 14 +++++++------- test/unit/math/mix/fun/inc_beta_inv_test.cpp | 12 ++++++------ test/unit/math/rev/fun/inc_beta_inv_test.cpp | 11 +++++------ 5 files changed, 35 insertions(+), 35 deletions(-) diff --git a/stan/math/fwd/fun/inc_beta_inv.hpp b/stan/math/fwd/fun/inc_beta_inv.hpp index 274ce485b87..3724ebdc995 100644 --- a/stan/math/fwd/fun/inc_beta_inv.hpp +++ b/stan/math/fwd/fun/inc_beta_inv.hpp @@ -87,13 +87,13 @@ inline fvar> inc_beta_inv(const T1& a, T_return inv_d_(0); if(is_fvar::value) { - auto da1 = one_m_b * log1m_w + one_m_a * log_w; - auto da2 = a_val * log_w + 2 * lgamma(a_val) + auto da1 = exp(one_m_b * log1m_w + one_m_a * log_w); + auto da2 = exp(a_val * log_w + 2 * lgamma(a_val) + log(F32(a_val, a_val, one_m_b, ap1, ap1, w)) - - 2 * lgamma(ap1); - auto da3 = lbeta_ab + log(inc_beta(a_val, b_val, w)) - + log(log_w - digamma(a_val) + digamma_apb); - inv_d_ += forward_as>(a).d_ * exp(da1 + log_diff_exp(da2, da3)); + - 2 * lgamma(ap1)); + auto da3 = inc_beta(a_val, b_val, w) * exp(lbeta_ab) + * (log_w - digamma(a_val) + digamma_apb); + inv_d_ += forward_as>(a).d_ * da1 * (da2 - da3); } if(is_fvar::value) { @@ -103,10 +103,11 @@ inline fvar> inc_beta_inv(const T1& a, - 2 * lgamma(bp1) + b_val * log1m_w; - auto db3 = log(inc_beta(b_val, a_val, one_m_w)) + lbeta_ab - + log(log1m_w - digamma(b_val) + digamma_apb); + auto db3 = inc_beta(b_val, a_val, one_m_w) * exp(lbeta_ab) + * (log1m_w - digamma(b_val) + digamma_apb); - inv_d_ += forward_as>(b).d_ * db1 * exp(log_diff_exp(db2, db3)); + + inv_d_ += forward_as>(b).d_ * db1 * (exp(db2) - db3); } if(is_fvar::value) { diff --git a/stan/math/rev/fun/inc_beta_inv.hpp b/stan/math/rev/fun/inc_beta_inv.hpp index b804543602f..c3fe853d205 100644 --- a/stan/math/rev/fun/inc_beta_inv.hpp +++ b/stan/math/rev/fun/inc_beta_inv.hpp @@ -56,14 +56,14 @@ inline var inc_beta_inv(const T1& a, const T2& b, const T3& p) { double digamma_apb = digamma(a_val + b_val); if(!is_constant_all::value) { - double da1 = one_m_b * log1m_w + one_m_a * log_w; + double da1 = exp(one_m_b * log1m_w + one_m_a * log_w); double da2 = a_val * log_w + 2 * lgamma(a_val) + log(F32(a_val, a_val, one_m_b, ap1, ap1, w)) - 2 * lgamma(ap1); - double da3 = lbeta_ab + log(inc_beta(a_val, b_val, w)) - + log(log_w - digamma(a_val) + digamma_apb); + double da3 = inc_beta(a_val, b_val, w) * exp(lbeta_ab) + * (log_w - digamma(a_val) + digamma_apb); - forward_as(a).adj() += vi.adj() * exp(da1 + log_diff_exp(da2, da3)); + forward_as(a).adj() += vi.adj() * da1 * (exp(da2) - da3); } if(!is_constant_all::value) { @@ -73,10 +73,10 @@ inline var inc_beta_inv(const T1& a, const T2& b, const T3& p) { - 2 * lgamma(bp1) + b_val * log1m_w; - double db3 = log(inc_beta(b_val, a_val, one_m_w)) + lbeta_ab - + log(log1m_w - digamma(b_val) + digamma_apb); + double db3 = inc_beta(b_val, a_val, one_m_w) * exp(lbeta_ab) + * (log1m_w - digamma(b_val) + digamma_apb); - forward_as(b).adj() += vi.adj() * db1 * exp(log_diff_exp(db2, db3)); + forward_as(b).adj() += vi.adj() * db1 * (exp(db2) - db3); } if(!is_constant_all::value) { diff --git a/test/unit/math/fwd/fun/inc_beta_inv_test.cpp b/test/unit/math/fwd/fun/inc_beta_inv_test.cpp index 18be7930306..288e9870596 100644 --- a/test/unit/math/fwd/fun/inc_beta_inv_test.cpp +++ b/test/unit/math/fwd/fun/inc_beta_inv_test.cpp @@ -4,29 +4,29 @@ TEST(AgradFwdMatrixIncBetaInv, fd_scalar) { using stan::math::fvar; using stan::math::inc_beta_inv; - fvar a = 1; + fvar a = 6; fvar b = 2; - fvar p = 0.5; + fvar p = 0.9; a.d_ = 1.0; b.d_ = 1.0; p.d_ = 1.0; fvar res = inc_beta_inv(a, b, p); - EXPECT_FLOAT_EQ(res.d_, 0.287698278597 - 0.122532267934 + 0.707106781187); + EXPECT_FLOAT_EQ(res.d_, 0.0117172527399 - 0.0680999818473 + 0.455387298585); } TEST(AgradFwdMatrixIncBetaInv, ffd_scalar) { using stan::math::fvar; using stan::math::inc_beta_inv; - fvar> a = 1; - fvar> b = 2; - fvar> p = 0.5; + fvar> a = 7; + fvar> b = 4; + fvar> p = 0.15; a.val_.d_ = 1.0; b.val_.d_ = 1.0; p.val_.d_ = 1.0; fvar> res = inc_beta_inv(a, b, p); - EXPECT_FLOAT_EQ(res.val_.d_, 0.287698278597 - 0.122532267934 + 0.707106781187); + EXPECT_FLOAT_EQ(res.val_.d_, 0.0428905418857 - 0.0563420377808 + 0.664919819507); } \ No newline at end of file diff --git a/test/unit/math/mix/fun/inc_beta_inv_test.cpp b/test/unit/math/mix/fun/inc_beta_inv_test.cpp index 236fecf664f..15ec5e01749 100644 --- a/test/unit/math/mix/fun/inc_beta_inv_test.cpp +++ b/test/unit/math/mix/fun/inc_beta_inv_test.cpp @@ -63,9 +63,9 @@ TEST(ProbInternalMath, inc_beta_inv_fv2) { using stan::math::fvar; using stan::math::var; using stan::math::inc_beta_inv; - fvar> a = 1; - fvar> b = 2; - fvar> p = 0.5; + fvar> a = 2; + fvar> b = 5; + fvar> p = 0.1; a.d_ = 1.0; b.d_ = 1.0; p.d_ = 1.0; @@ -73,7 +73,7 @@ TEST(ProbInternalMath, inc_beta_inv_fv2) { fvar> res = inc_beta_inv(a, b, p); res.val_.val_.grad(); - EXPECT_FLOAT_EQ(a.val_.val_.adj(), 0.287698278597); - EXPECT_FLOAT_EQ(b.val_.val_.adj(), -0.122532267934); - EXPECT_FLOAT_EQ(p.val_.val_.adj(), 0.707106781187); + EXPECT_FLOAT_EQ(a.val_.val_.adj(), 0.0783025374798); + EXPECT_FLOAT_EQ(b.val_.val_.adj(), -0.0161882044585); + EXPECT_FLOAT_EQ(p.val_.val_.adj(), 0.530989359806); } \ No newline at end of file diff --git a/test/unit/math/rev/fun/inc_beta_inv_test.cpp b/test/unit/math/rev/fun/inc_beta_inv_test.cpp index 6ef31ea657d..0c1c24b4fca 100644 --- a/test/unit/math/rev/fun/inc_beta_inv_test.cpp +++ b/test/unit/math/rev/fun/inc_beta_inv_test.cpp @@ -7,15 +7,14 @@ TEST(inc_beta_inv, values) { using stan::math::var; using stan::math::inc_beta_inv; - // check autodiffing works with var types with large values - var a = 1; + var a = 25; var b = 2; - var p = 0.5; + var p = 0.7; var res = inc_beta_inv(a, b, p); res.grad(); - EXPECT_FLOAT_EQ(a.adj(), 0.287698278597); - EXPECT_FLOAT_EQ(b.adj(), -0.122532267934); - EXPECT_FLOAT_EQ(p.adj(), 0.707106781187); + EXPECT_FLOAT_EQ(a.adj(), 0.00161775443606); + EXPECT_FLOAT_EQ(b.adj(), -0.0289198999936); + EXPECT_FLOAT_EQ(p.adj(), 0.102597843424); } From a444d9c6e5b0a5cbded9527124b7aed83cfe79c7 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 31 Dec 2021 13:10:28 +0800 Subject: [PATCH 03/12] cpplint --- stan/math/fwd/fun/inc_beta_inv.hpp | 59 ++++------------ stan/math/prim/fun/F32.hpp | 9 ++- stan/math/rev/fun/inc_beta_inv.hpp | 34 ++++----- test/unit/math/fwd/fun/inc_beta_inv_test.cpp | 5 +- test/unit/math/mix/fun/inc_beta_inv_test.cpp | 2 +- test/unit/math/prim/fun/inc_beta_inv_test.cpp | 70 +++++++++++++++++++ 6 files changed, 106 insertions(+), 73 deletions(-) create mode 100644 test/unit/math/prim/fun/inc_beta_inv_test.cpp diff --git a/stan/math/fwd/fun/inc_beta_inv.hpp b/stan/math/fwd/fun/inc_beta_inv.hpp index 3724ebdc995..d82ced9f6e9 100644 --- a/stan/math/fwd/fun/inc_beta_inv.hpp +++ b/stan/math/fwd/fun/inc_beta_inv.hpp @@ -17,51 +17,16 @@ namespace stan { namespace math { /** - * The fused multiply-add operation (C99). + * The inverse of the normalized incomplete beta function of a, b, with probability p. * - * This double-based operation delegates to fma. + * Used to compute the cumulative density function for the beta + * distribution. * - * The function is defined by - * - * fma(a, b, c) = (a * b) + c. - * - * - \f[ - \mbox{fma}(x, y, z) = - \begin{cases} - x\cdot y+z & \mbox{if } -\infty\leq x, y, z \leq \infty \\[6pt] - \textrm{NaN} & \mbox{if } x = \textrm{NaN} - \end{cases} - \f] - - \f[ - \frac{\partial\, \mbox{fma}(x, y, z)}{\partial x} = - \begin{cases} - y & \mbox{if } -\infty\leq x, y, z \leq \infty \\[6pt] - \textrm{NaN} & \mbox{if } x = \textrm{NaN} - \end{cases} - \f] - - \f[ - \frac{\partial\, \mbox{fma}(x, y, z)}{\partial y} = - \begin{cases} - x & \mbox{if } -\infty\leq x, y, z \leq \infty \\[6pt] - \textrm{NaN} & \mbox{if } x = \textrm{NaN} - \end{cases} - \f] - - \f[ - \frac{\partial\, \mbox{fma}(x, y, z)}{\partial z} = - \begin{cases} - 1 & \mbox{if } -\infty\leq x, y, z \leq \infty \\[6pt] - \textrm{NaN} & \mbox{if } x = \textrm{NaN} - \end{cases} - \f] - * - * @param x1 First value. - * @param x2 Second value. - * @param x3 Third value. - * @return Product of the first two values plus the third. + * @param a Shape parameter a >= 0; a and b can't both be 0 + * @param b Shape parameter b >= 0 + * @param p Random variate. 0 <= p <= 1 + * @throws if constraints are violated or if any argument is NaN + * @return The inverse of the normalized incomplete beta function. */ template * = nullptr, @@ -86,7 +51,7 @@ inline fvar> inc_beta_inv(const T1& a, T_return inv_d_(0); - if(is_fvar::value) { + if (is_fvar::value) { auto da1 = exp(one_m_b * log1m_w + one_m_a * log_w); auto da2 = exp(a_val * log_w + 2 * lgamma(a_val) + log(F32(a_val, a_val, one_m_b, ap1, ap1, w)) @@ -96,9 +61,9 @@ inline fvar> inc_beta_inv(const T1& a, inv_d_ += forward_as>(a).d_ * da1 * (da2 - da3); } - if(is_fvar::value) { + if (is_fvar::value) { auto db1 = (w - 1) * exp(-b_val * log1m_w + one_m_a * log_w); - auto db2 = 2 * lgamma(b_val) + auto db2 = 2 * lgamma(b_val) + log(F32(b_val, b_val, one_m_a, bp1, bp1, one_m_w)) - 2 * lgamma(bp1) + b_val * log1m_w; @@ -110,7 +75,7 @@ inline fvar> inc_beta_inv(const T1& a, inv_d_ += forward_as>(b).d_ * db1 * (exp(db2) - db3); } - if(is_fvar::value) { + if (is_fvar::value) { inv_d_ += forward_as>(p).d_ * exp(one_m_b * log1m_w + one_m_a * log_w + lbeta_ab); } diff --git a/stan/math/prim/fun/F32.hpp b/stan/math/prim/fun/F32.hpp index ecccd9bcaa1..ac5688652d0 100644 --- a/stan/math/prim/fun/F32.hpp +++ b/stan/math/prim/fun/F32.hpp @@ -49,8 +49,10 @@ namespace math { * @param[in] precision precision of the infinite sum. defaults to 1e-6 * @param[in] max_steps number of steps to take. defaults to 1e5 */ -template -return_type_t F32(const Ta1& a1, const Ta2& a2, const Ta3& a3, const Tb1& b1, const Tb2& b2, +template +return_type_t + F32(const Ta1& a1, const Ta2& a2, const Ta3& a3, const Tb1& b1, const Tb2& b2, const Tz& z, double precision = 1e-6, int max_steps = 1e5) { check_3F2_converges("F32", a1, a2, a3, b1, b2, z); @@ -65,7 +67,8 @@ return_type_t F32(const Ta1& a1, const Ta2& a2, con double t_sign = 1.0; for (int k = 0; k <= max_steps; ++k) { - T_return p = (a1 + k) * (a2 + k) * (a3 + k) / ((b1 + k) * (b2 + k) * (k + 1)); + T_return p = (a1 + k) * (a2 + k) * (a3 + k) + / ((b1 + k) * (b2 + k) * (k + 1)); if (p == 0.0) { return t_acc; } diff --git a/stan/math/rev/fun/inc_beta_inv.hpp b/stan/math/rev/fun/inc_beta_inv.hpp index c3fe853d205..54d88aa6f13 100644 --- a/stan/math/rev/fun/inc_beta_inv.hpp +++ b/stan/math/rev/fun/inc_beta_inv.hpp @@ -19,22 +19,16 @@ namespace stan { namespace math { /** - * The fused multiply-add function for three variables (C99). - * This function returns the product of the first two arguments - * plus the third argument. + * The inverse of the normalized incomplete beta function of a, b, with probability p. * - * The partial derivatives are + * Used to compute the cumulative density function for the beta + * distribution. * - * \f$\frac{\partial}{\partial x} (x * y) + z = y\f$, and - * - * \f$\frac{\partial}{\partial y} (x * y) + z = x\f$, and - * - * \f$\frac{\partial}{\partial z} (x * y) + z = 1\f$. - * - * @param x First multiplicand. - * @param y Second multiplicand. - * @param z Summand. - * @return Product of the multiplicands plus the summand, ($a * $b) + $c. + * @param a Shape parameter a >= 0; a and b can't both be 0 + * @param b Shape parameter b >= 0 + * @param p Random variate. 0 <= p <= 1 + * @throws if constraints are violated or if any argument is NaN + * @return The inverse of the normalized incomplete beta function. */ template * = nullptr, @@ -55,7 +49,7 @@ inline var inc_beta_inv(const T1& a, const T2& b, const T3& p) { double lbeta_ab = lbeta(a_val, b_val); double digamma_apb = digamma(a_val + b_val); - if(!is_constant_all::value) { + if (!is_constant_all::value) { double da1 = exp(one_m_b * log1m_w + one_m_a * log_w); double da2 = a_val * log_w + 2 * lgamma(a_val) + log(F32(a_val, a_val, one_m_b, ap1, ap1, w)) @@ -66,9 +60,9 @@ inline var inc_beta_inv(const T1& a, const T2& b, const T3& p) { forward_as(a).adj() += vi.adj() * da1 * (exp(da2) - da3); } - if(!is_constant_all::value) { + if (!is_constant_all::value) { double db1 = (w - 1) * exp(-b_val * log1m_w + one_m_a * log_w); - double db2 = 2 * lgamma(b_val) + double db2 = 2 * lgamma(b_val) + log(F32(b_val, b_val, one_m_a, bp1, bp1, one_m_w)) - 2 * lgamma(bp1) + b_val * log1m_w; @@ -77,12 +71,12 @@ inline var inc_beta_inv(const T1& a, const T2& b, const T3& p) { * (log1m_w - digamma(b_val) + digamma_apb); forward_as(b).adj() += vi.adj() * db1 * (exp(db2) - db3); - } + } - if(!is_constant_all::value) { + if (!is_constant_all::value) { forward_as(p).adj() += vi.adj() * exp(one_m_b * log1m_w + one_m_a * log_w + lbeta_ab); - } + } }); } diff --git a/test/unit/math/fwd/fun/inc_beta_inv_test.cpp b/test/unit/math/fwd/fun/inc_beta_inv_test.cpp index 288e9870596..9ff78cc9c1e 100644 --- a/test/unit/math/fwd/fun/inc_beta_inv_test.cpp +++ b/test/unit/math/fwd/fun/inc_beta_inv_test.cpp @@ -28,5 +28,6 @@ TEST(AgradFwdMatrixIncBetaInv, ffd_scalar) { fvar> res = inc_beta_inv(a, b, p); - EXPECT_FLOAT_EQ(res.val_.d_, 0.0428905418857 - 0.0563420377808 + 0.664919819507); -} \ No newline at end of file + EXPECT_FLOAT_EQ(res.val_.d_, + 0.0428905418857 - 0.0563420377808 + 0.664919819507); +} diff --git a/test/unit/math/mix/fun/inc_beta_inv_test.cpp b/test/unit/math/mix/fun/inc_beta_inv_test.cpp index 15ec5e01749..19381218bec 100644 --- a/test/unit/math/mix/fun/inc_beta_inv_test.cpp +++ b/test/unit/math/mix/fun/inc_beta_inv_test.cpp @@ -76,4 +76,4 @@ TEST(ProbInternalMath, inc_beta_inv_fv2) { EXPECT_FLOAT_EQ(a.val_.val_.adj(), 0.0783025374798); EXPECT_FLOAT_EQ(b.val_.val_.adj(), -0.0161882044585); EXPECT_FLOAT_EQ(p.val_.val_.adj(), 0.530989359806); -} \ No newline at end of file +} diff --git a/test/unit/math/prim/fun/inc_beta_inv_test.cpp b/test/unit/math/prim/fun/inc_beta_inv_test.cpp new file mode 100644 index 00000000000..815e4ba0e91 --- /dev/null +++ b/test/unit/math/prim/fun/inc_beta_inv_test.cpp @@ -0,0 +1,70 @@ +#include +#include +#include + +TEST(MathFunctions, inc_beta_inv_inv) { + using stan::math::inc_beta_inv; + + EXPECT_FLOAT_EQ(0.146446609407, inc_beta_inv(0.5, 0.5, 0.25)); + EXPECT_FLOAT_EQ(0.5, inc_beta_inv(0.5, 0.5, 0.5)); + EXPECT_FLOAT_EQ(0.853553390593, inc_beta_inv(0.5, 0.5, 0.75)); + EXPECT_FLOAT_EQ(1.0, inc_beta_inv(0.5, 0.5, 1.0)); + + EXPECT_FLOAT_EQ(0.0, inc_beta_inv(0.1, 1.5, 0.0)); + EXPECT_FLOAT_EQ(5.33624995327e-7, inc_beta_inv(0.1, 1.5, 0.25)); + EXPECT_FLOAT_EQ(0.000546567646391, inc_beta_inv(0.1, 1.5, 0.5)); + EXPECT_FLOAT_EQ(0.0319736161201, inc_beta_inv(0.1, 1.5, 0.75)); + EXPECT_FLOAT_EQ(1.0, inc_beta_inv(0.1, 1.5, 1.0)); + + EXPECT_FLOAT_EQ(0.804032164227, inc_beta_inv(0.6, 0.3, 0.5)); +} + +TEST(MathFunctions, inc_beta_inv_a_boundary) { + double b = 0.5; + double p = 0.5; + + const double inf = std::numeric_limits::infinity(); + + EXPECT_THROW(stan::math::inc_beta_inv(0.0, b, p), std::domain_error); + EXPECT_NO_THROW(stan::math::inc_beta_inv(inf, b, p)); + EXPECT_THROW(stan::math::inc_beta_inv(-0.01, b, p), std::domain_error); +} + +TEST(MathFunctions, inc_beta_inv_b_boundary) { + double a = 0.5; + double p = 0.5; + + const double inf = std::numeric_limits::infinity(); + + EXPECT_THROW(stan::math::inc_beta_inv(a, 0.0, p), std::domain_error); + EXPECT_NO_THROW(stan::math::inc_beta_inv(inf, 1.0, p)); + EXPECT_THROW(stan::math::inc_beta_inv(a, -0.01, p), std::domain_error); +} + +TEST(MathFunctions, inc_beta_inv_x_boundary) { + double a = 0.5; + double b = 0.5; + + EXPECT_NO_THROW(stan::math::inc_beta_inv(a, b, 0.0)); + EXPECT_NO_THROW(stan::math::inc_beta_inv(a, b, 1.0)); + EXPECT_THROW(stan::math::inc_beta_inv(a, b, -0.01), std::domain_error); + EXPECT_THROW(stan::math::inc_beta_inv(a, b, 1.01), std::domain_error); +} + +TEST(MathFunctions, inc_beta_inv_a_b_boundary) { + double p = 0.5; + + EXPECT_THROW(stan::math::inc_beta_inv(0.0, 0.0, p), std::domain_error); +} + +TEST(MathFunctions, inc_beta_inv_nan) { + double nan = std::numeric_limits::quiet_NaN(); + + EXPECT_THROW(stan::math::inc_beta_inv(0.5, 0.0, nan), std::domain_error); + EXPECT_THROW(stan::math::inc_beta_inv(0.5, nan, 0.0), std::domain_error); + EXPECT_THROW(stan::math::inc_beta_inv(0.5, nan, nan), std::domain_error); + EXPECT_THROW(stan::math::inc_beta_inv(nan, 0.0, 0.0), std::domain_error); + EXPECT_THROW(stan::math::inc_beta_inv(nan, 0.0, nan), std::domain_error); + EXPECT_THROW(stan::math::inc_beta_inv(nan, nan, 0.0), std::domain_error); + EXPECT_THROW(stan::math::inc_beta_inv(nan, nan, nan), std::domain_error); +} From be10487b369912d466dd93d16f7ace0a9911e122 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Fri, 31 Dec 2021 05:19:53 +0000 Subject: [PATCH 04/12] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/fwd/fun/inc_beta_inv.hpp | 30 +++++++++----------- stan/math/prim/fun/F32.hpp | 16 ++++++----- stan/math/prim/fun/inc_beta_inv.hpp | 3 +- stan/math/rev/fun/inc_beta_inv.hpp | 11 ++++--- test/unit/math/mix/fun/inc_beta_inv_test.cpp | 4 +-- test/unit/math/rev/fun/inc_beta_inv_test.cpp | 2 +- 6 files changed, 33 insertions(+), 33 deletions(-) diff --git a/stan/math/fwd/fun/inc_beta_inv.hpp b/stan/math/fwd/fun/inc_beta_inv.hpp index d82ced9f6e9..31a4674870b 100644 --- a/stan/math/fwd/fun/inc_beta_inv.hpp +++ b/stan/math/fwd/fun/inc_beta_inv.hpp @@ -17,7 +17,8 @@ namespace stan { namespace math { /** - * The inverse of the normalized incomplete beta function of a, b, with probability p. + * The inverse of the normalized incomplete beta function of a, b, with + * probability p. * * Used to compute the cumulative density function for the beta * distribution. @@ -32,8 +33,8 @@ template * = nullptr, require_any_fvar_t* = nullptr> inline fvar> inc_beta_inv(const T1& a, - const T2& b, - const T3& p) { + const T2& b, + const T3& p) { using T_return = partials_return_t; auto a_val = value_of(a); auto b_val = value_of(b); @@ -53,37 +54,34 @@ inline fvar> inc_beta_inv(const T1& a, if (is_fvar::value) { auto da1 = exp(one_m_b * log1m_w + one_m_a * log_w); - auto da2 = exp(a_val * log_w + 2 * lgamma(a_val) - + log(F32(a_val, a_val, one_m_b, ap1, ap1, w)) - - 2 * lgamma(ap1)); - auto da3 = inc_beta(a_val, b_val, w) * exp(lbeta_ab) - * (log_w - digamma(a_val) + digamma_apb); + auto da2 + = exp(a_val * log_w + 2 * lgamma(a_val) + + log(F32(a_val, a_val, one_m_b, ap1, ap1, w)) - 2 * lgamma(ap1)); + auto da3 = inc_beta(a_val, b_val, w) * exp(lbeta_ab) + * (log_w - digamma(a_val) + digamma_apb); inv_d_ += forward_as>(a).d_ * da1 * (da2 - da3); } if (is_fvar::value) { auto db1 = (w - 1) * exp(-b_val * log1m_w + one_m_a * log_w); auto db2 = 2 * lgamma(b_val) - + log(F32(b_val, b_val, one_m_a, bp1, bp1, one_m_w)) - - 2 * lgamma(bp1) - + b_val * log1m_w; + + log(F32(b_val, b_val, one_m_a, bp1, bp1, one_m_w)) + - 2 * lgamma(bp1) + b_val * log1m_w; auto db3 = inc_beta(b_val, a_val, one_m_w) * exp(lbeta_ab) - * (log1m_w - digamma(b_val) + digamma_apb); - + * (log1m_w - digamma(b_val) + digamma_apb); inv_d_ += forward_as>(b).d_ * db1 * (exp(db2) - db3); } if (is_fvar::value) { - inv_d_ += forward_as>(p).d_ * - exp(one_m_b * log1m_w + one_m_a * log_w + lbeta_ab); + inv_d_ += forward_as>(p).d_ + * exp(one_m_b * log1m_w + one_m_a * log_w + lbeta_ab); } return fvar(w, inv_d_); } - } // namespace math } // namespace stan #endif diff --git a/stan/math/prim/fun/F32.hpp b/stan/math/prim/fun/F32.hpp index ac5688652d0..bb48039bdb3 100644 --- a/stan/math/prim/fun/F32.hpp +++ b/stan/math/prim/fun/F32.hpp @@ -49,11 +49,13 @@ namespace math { * @param[in] precision precision of the infinite sum. defaults to 1e-6 * @param[in] max_steps number of steps to take. defaults to 1e5 */ -template -return_type_t - F32(const Ta1& a1, const Ta2& a2, const Ta3& a3, const Tb1& b1, const Tb2& b2, - const Tz& z, double precision = 1e-6, int max_steps = 1e5) { +template +return_type_t F32(const Ta1& a1, const Ta2& a2, + const Ta3& a3, const Tb1& b1, + const Tb2& b2, const Tz& z, + double precision = 1e-6, + int max_steps = 1e5) { check_3F2_converges("F32", a1, a2, a3, b1, b2, z); using T_return = return_type_t; @@ -67,8 +69,8 @@ return_type_t double t_sign = 1.0; for (int k = 0; k <= max_steps; ++k) { - T_return p = (a1 + k) * (a2 + k) * (a3 + k) - / ((b1 + k) * (b2 + k) * (k + 1)); + T_return p + = (a1 + k) * (a2 + k) * (a3 + k) / ((b1 + k) * (b2 + k) * (k + 1)); if (p == 0.0) { return t_acc; } diff --git a/stan/math/prim/fun/inc_beta_inv.hpp b/stan/math/prim/fun/inc_beta_inv.hpp index 0ecf33cf291..9d0301e3b6d 100644 --- a/stan/math/prim/fun/inc_beta_inv.hpp +++ b/stan/math/prim/fun/inc_beta_inv.hpp @@ -10,7 +10,8 @@ namespace stan { namespace math { /** - * The inverse of the normalized incomplete beta function of a, b, with probability p. + * The inverse of the normalized incomplete beta function of a, b, with + * probability p. * * Used to compute the cumulative density function for the beta * distribution. diff --git a/stan/math/rev/fun/inc_beta_inv.hpp b/stan/math/rev/fun/inc_beta_inv.hpp index 54d88aa6f13..a26e5e08893 100644 --- a/stan/math/rev/fun/inc_beta_inv.hpp +++ b/stan/math/rev/fun/inc_beta_inv.hpp @@ -19,7 +19,8 @@ namespace stan { namespace math { /** - * The inverse of the normalized incomplete beta function of a, b, with probability p. + * The inverse of the normalized incomplete beta function of a, b, with + * probability p. * * Used to compute the cumulative density function for the beta * distribution. @@ -64,8 +65,7 @@ inline var inc_beta_inv(const T1& a, const T2& b, const T3& p) { double db1 = (w - 1) * exp(-b_val * log1m_w + one_m_a * log_w); double db2 = 2 * lgamma(b_val) + log(F32(b_val, b_val, one_m_a, bp1, bp1, one_m_w)) - - 2 * lgamma(bp1) - + b_val * log1m_w; + - 2 * lgamma(bp1) + b_val * log1m_w; double db3 = inc_beta(b_val, a_val, one_m_w) * exp(lbeta_ab) * (log1m_w - digamma(b_val) + digamma_apb); @@ -74,13 +74,12 @@ inline var inc_beta_inv(const T1& a, const T2& b, const T3& p) { } if (!is_constant_all::value) { - forward_as(p).adj() += vi.adj() * - exp(one_m_b * log1m_w + one_m_a * log_w + lbeta_ab); + forward_as(p).adj() + += vi.adj() * exp(one_m_b * log1m_w + one_m_a * log_w + lbeta_ab); } }); } - } // namespace math } // namespace stan #endif diff --git a/test/unit/math/mix/fun/inc_beta_inv_test.cpp b/test/unit/math/mix/fun/inc_beta_inv_test.cpp index 19381218bec..271b1178057 100644 --- a/test/unit/math/mix/fun/inc_beta_inv_test.cpp +++ b/test/unit/math/mix/fun/inc_beta_inv_test.cpp @@ -4,8 +4,8 @@ TEST(ProbInternalMath, inc_beta_inv_fv1) { using stan::math::fvar; - using stan::math::var; using stan::math::inc_beta_inv; + using stan::math::var; double a_d = 1; double b_d = 2; double p_d = 0.5; @@ -61,8 +61,8 @@ TEST(ProbInternalMath, inc_beta_inv_fv1) { TEST(ProbInternalMath, inc_beta_inv_fv2) { using stan::math::fvar; - using stan::math::var; using stan::math::inc_beta_inv; + using stan::math::var; fvar> a = 2; fvar> b = 5; fvar> p = 0.1; diff --git a/test/unit/math/rev/fun/inc_beta_inv_test.cpp b/test/unit/math/rev/fun/inc_beta_inv_test.cpp index 0c1c24b4fca..389bb24b13a 100644 --- a/test/unit/math/rev/fun/inc_beta_inv_test.cpp +++ b/test/unit/math/rev/fun/inc_beta_inv_test.cpp @@ -4,8 +4,8 @@ #include TEST(inc_beta_inv, values) { - using stan::math::var; using stan::math::inc_beta_inv; + using stan::math::var; var a = 25; var b = 2; From 2b887739f5acf351dd49d0fbb1bd11ad955022d0 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sun, 2 Jan 2022 14:37:58 +0800 Subject: [PATCH 05/12] Update doc & latex --- stan/math/prim/fun/inc_beta_inv.hpp | 2 +- stan/math/rev/fun/inc_beta_inv.hpp | 20 ++++++++++++++++++++ 2 files changed, 21 insertions(+), 1 deletion(-) diff --git a/stan/math/prim/fun/inc_beta_inv.hpp b/stan/math/prim/fun/inc_beta_inv.hpp index 9d0301e3b6d..d533c683640 100644 --- a/stan/math/prim/fun/inc_beta_inv.hpp +++ b/stan/math/prim/fun/inc_beta_inv.hpp @@ -13,7 +13,7 @@ namespace math { * The inverse of the normalized incomplete beta function of a, b, with * probability p. * - * Used to compute the cumulative density function for the beta + * Used to compute the inverse cumulative density function for the beta * distribution. * * @param a Shape parameter a >= 0; a and b can't both be 0 diff --git a/stan/math/rev/fun/inc_beta_inv.hpp b/stan/math/rev/fun/inc_beta_inv.hpp index a26e5e08893..3910c2e82dc 100644 --- a/stan/math/rev/fun/inc_beta_inv.hpp +++ b/stan/math/rev/fun/inc_beta_inv.hpp @@ -24,6 +24,26 @@ namespace math { * * Used to compute the cumulative density function for the beta * distribution. + * + \f[ + \frac{\partial }{\partial a} = + (1-w)^{1-b}w^{1-a} + \left( + w^a\Gamma(a)^2 {}_3\tilde{F}_2(a,a,1-b;a+1,a+1;w) + - B(a,b)I_w(a,b)\left(\log(w)-\psi(a) + \psi(a+b)\right) + \right)/;w=I_z^{-1}(a,b) + \f] + \f[ + \frac{\partial }{\partial b} = + (1-w)^{-b}w^{1-a}(w-1) + \left( + (1-w)^{b}\Gamma(b)^2 {}_3\tilde{F}_2(b,b,1-a;b+1,b+1;1-w) + - B_{1-w}(b,a)\left(\log(1-w)-\psi(b) + \psi(a+b)\right) + \right)/;w=I_z^{-1}(a,b) + \f] + \f[ + \frac{\partial }{\partial z} = (1-w)^{1-b}w^{1-a}B(a,b)/;w=I_z^{-1}(a,b) + \f] * * @param a Shape parameter a >= 0; a and b can't both be 0 * @param b Shape parameter b >= 0 From 39a74219bf9317332d5f12ddf3b1756fb98248b8 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sun, 2 Jan 2022 09:52:34 +0000 Subject: [PATCH 06/12] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/inc_beta_inv.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/rev/fun/inc_beta_inv.hpp b/stan/math/rev/fun/inc_beta_inv.hpp index 3910c2e82dc..fe701e2249d 100644 --- a/stan/math/rev/fun/inc_beta_inv.hpp +++ b/stan/math/rev/fun/inc_beta_inv.hpp @@ -24,7 +24,7 @@ namespace math { * * Used to compute the cumulative density function for the beta * distribution. - * + * \f[ \frac{\partial }{\partial a} = (1-w)^{1-b}w^{1-a} From 085f8ea88b3657aac89f3ae277a7e83cf771c262 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 3 Jan 2022 14:13:25 +0800 Subject: [PATCH 07/12] Missing doc updates --- stan/math/fwd/fun/inc_beta_inv.hpp | 2 +- stan/math/prim/fun/inc_beta_inv.hpp | 6 +++--- stan/math/rev/fun/inc_beta_inv.hpp | 2 +- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/stan/math/fwd/fun/inc_beta_inv.hpp b/stan/math/fwd/fun/inc_beta_inv.hpp index 31a4674870b..0ae0e002d6a 100644 --- a/stan/math/fwd/fun/inc_beta_inv.hpp +++ b/stan/math/fwd/fun/inc_beta_inv.hpp @@ -20,7 +20,7 @@ namespace math { * The inverse of the normalized incomplete beta function of a, b, with * probability p. * - * Used to compute the cumulative density function for the beta + * Used to compute the inverse cumulative density function for the beta * distribution. * * @param a Shape parameter a >= 0; a and b can't both be 0 diff --git a/stan/math/prim/fun/inc_beta_inv.hpp b/stan/math/prim/fun/inc_beta_inv.hpp index d533c683640..e351ea2fae5 100644 --- a/stan/math/prim/fun/inc_beta_inv.hpp +++ b/stan/math/prim/fun/inc_beta_inv.hpp @@ -23,9 +23,9 @@ namespace math { * @return The inverse of the normalized incomplete beta function. */ inline double inc_beta_inv(double a, double b, double p) { - check_not_nan("inc_beta", "a", a); - check_not_nan("inc_beta", "b", b); - check_not_nan("inc_beta", "p", p); + check_not_nan("inc_beta_inv", "a", a); + check_not_nan("inc_beta_inv", "b", b); + check_not_nan("inc_beta_inv", "p", p); return boost::math::ibeta_inv(a, b, p, boost_policy_t<>()); } diff --git a/stan/math/rev/fun/inc_beta_inv.hpp b/stan/math/rev/fun/inc_beta_inv.hpp index fe701e2249d..0d349c682db 100644 --- a/stan/math/rev/fun/inc_beta_inv.hpp +++ b/stan/math/rev/fun/inc_beta_inv.hpp @@ -22,7 +22,7 @@ namespace math { * The inverse of the normalized incomplete beta function of a, b, with * probability p. * - * Used to compute the cumulative density function for the beta + * Used to compute the inverse cumulative density function for the beta * distribution. * \f[ From 1eef220d57b6109bdabbfcccdf9c193f6640e4ce Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 5 Jan 2022 21:38:54 +0800 Subject: [PATCH 08/12] Update inv naming --- stan/math/fwd/fun.hpp | 2 +- .../{inc_beta_inv.hpp => inv_inc_beta.hpp} | 10 +-- stan/math/prim/fun.hpp | 2 +- .../{inc_beta_inv.hpp => inv_inc_beta.hpp} | 12 ++-- stan/math/rev/fun.hpp | 2 +- .../{inc_beta_inv.hpp => inv_inc_beta.hpp} | 10 +-- ...eta_inv_test.cpp => inv_inc_beta_test.cpp} | 18 ++--- test/unit/math/prim/fun/inc_beta_inv_test.cpp | 70 ------------------- test/unit/math/prim/fun/inv_inc_beta_test.cpp | 70 +++++++++++++++++++ 9 files changed, 98 insertions(+), 98 deletions(-) rename stan/math/fwd/fun/{inc_beta_inv.hpp => inv_inc_beta.hpp} (91%) rename stan/math/prim/fun/{inc_beta_inv.hpp => inv_inc_beta.hpp} (74%) rename stan/math/rev/fun/{inc_beta_inv.hpp => inv_inc_beta.hpp} (93%) rename test/unit/math/mix/fun/{inc_beta_inv_test.cpp => inv_inc_beta_test.cpp} (80%) delete mode 100644 test/unit/math/prim/fun/inc_beta_inv_test.cpp create mode 100644 test/unit/math/prim/fun/inv_inc_beta_test.cpp diff --git a/stan/math/fwd/fun.hpp b/stan/math/fwd/fun.hpp index 0c2d4559867..3185423a7ce 100644 --- a/stan/math/fwd/fun.hpp +++ b/stan/math/fwd/fun.hpp @@ -43,10 +43,10 @@ #include #include #include -#include #include #include #include +#include #include #include #include diff --git a/stan/math/fwd/fun/inc_beta_inv.hpp b/stan/math/fwd/fun/inv_inc_beta.hpp similarity index 91% rename from stan/math/fwd/fun/inc_beta_inv.hpp rename to stan/math/fwd/fun/inv_inc_beta.hpp index 0ae0e002d6a..a4d11b053ab 100644 --- a/stan/math/fwd/fun/inc_beta_inv.hpp +++ b/stan/math/fwd/fun/inv_inc_beta.hpp @@ -1,9 +1,9 @@ -#ifndef STAN_MATH_FWD_FUN_INC_BETA_INV_HPP -#define STAN_MATH_FWD_FUN_INC_BETA_INV_HPP +#ifndef STAN_MATH_FWD_FUN_INV_INC_BETA_HPP +#define STAN_MATH_FWD_FUN_INV_INC_BETA_HPP #include #include -#include +#include #include #include #include @@ -32,14 +32,14 @@ namespace math { template * = nullptr, require_any_fvar_t* = nullptr> -inline fvar> inc_beta_inv(const T1& a, +inline fvar> inv_inc_beta(const T1& a, const T2& b, const T3& p) { using T_return = partials_return_t; auto a_val = value_of(a); auto b_val = value_of(b); auto p_val = value_of(p); - T_return w = inc_beta_inv(a_val, b_val, p_val); + T_return w = inv_inc_beta(a_val, b_val, p_val); T_return log_w = log(w); T_return log1m_w = log1m(w); auto one_m_a = 1 - a_val; diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp index 6c181e3ef6c..2a8502450cd 100644 --- a/stan/math/prim/fun.hpp +++ b/stan/math/prim/fun.hpp @@ -134,12 +134,12 @@ #include #include #include -#include #include #include #include #include #include +#include #include #include #include diff --git a/stan/math/prim/fun/inc_beta_inv.hpp b/stan/math/prim/fun/inv_inc_beta.hpp similarity index 74% rename from stan/math/prim/fun/inc_beta_inv.hpp rename to stan/math/prim/fun/inv_inc_beta.hpp index e351ea2fae5..bc81bad8f97 100644 --- a/stan/math/prim/fun/inc_beta_inv.hpp +++ b/stan/math/prim/fun/inv_inc_beta.hpp @@ -1,5 +1,5 @@ -#ifndef STAN_MATH_PRIM_FUN_INC_BETA_INV_HPP -#define STAN_MATH_PRIM_FUN_INC_BETA_INV_HPP +#ifndef STAN_MATH_PRIM_FUN_INV_INC_BETA_HPP +#define STAN_MATH_PRIM_FUN_INV_INC_BETA_HPP #include #include @@ -22,10 +22,10 @@ namespace math { * @throws if constraints are violated or if any argument is NaN * @return The inverse of the normalized incomplete beta function. */ -inline double inc_beta_inv(double a, double b, double p) { - check_not_nan("inc_beta_inv", "a", a); - check_not_nan("inc_beta_inv", "b", b); - check_not_nan("inc_beta_inv", "p", p); +inline double inv_inc_beta(double a, double b, double p) { + check_not_nan("inv_inc_beta", "a", a); + check_not_nan("inv_inc_beta", "b", b); + check_not_nan("inv_inc_beta", "p", p); return boost::math::ibeta_inv(a, b, p, boost_policy_t<>()); } diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp index a21bb201691..8549ad318b2 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -78,7 +78,7 @@ #include #include #include -#include +#include #include #include #include diff --git a/stan/math/rev/fun/inc_beta_inv.hpp b/stan/math/rev/fun/inv_inc_beta.hpp similarity index 93% rename from stan/math/rev/fun/inc_beta_inv.hpp rename to stan/math/rev/fun/inv_inc_beta.hpp index 0d349c682db..398ffb79866 100644 --- a/stan/math/rev/fun/inc_beta_inv.hpp +++ b/stan/math/rev/fun/inv_inc_beta.hpp @@ -1,10 +1,10 @@ -#ifndef STAN_MATH_REV_FUN_INC_BETA_INV_HPP -#define STAN_MATH_REV_FUN_INC_BETA_INV_HPP +#ifndef STAN_MATH_REV_FUN_INV_INC_BETA_HPP +#define STAN_MATH_REV_FUN_INV_INC_BETA_HPP #include #include #include -#include +#include #include #include #include @@ -54,11 +54,11 @@ namespace math { template * = nullptr, require_any_var_t* = nullptr> -inline var inc_beta_inv(const T1& a, const T2& b, const T3& p) { +inline var inv_inc_beta(const T1& a, const T2& b, const T3& p) { double a_val = value_of(a); double b_val = value_of(b); double p_val = value_of(p); - double w = inc_beta_inv(a_val, b_val, p_val); + double w = inv_inc_beta(a_val, b_val, p_val); return make_callback_var(w, [a, b, p, a_val, b_val, p_val, w](auto& vi) { double log_w = log(w); double log1m_w = log1m(w); diff --git a/test/unit/math/mix/fun/inc_beta_inv_test.cpp b/test/unit/math/mix/fun/inv_inc_beta_test.cpp similarity index 80% rename from test/unit/math/mix/fun/inc_beta_inv_test.cpp rename to test/unit/math/mix/fun/inv_inc_beta_test.cpp index 271b1178057..9286c4990d4 100644 --- a/test/unit/math/mix/fun/inc_beta_inv_test.cpp +++ b/test/unit/math/mix/fun/inv_inc_beta_test.cpp @@ -2,9 +2,9 @@ #include #include -TEST(ProbInternalMath, inc_beta_inv_fv1) { +TEST(ProbInternalMath, inv_inc_beta_fv1) { using stan::math::fvar; - using stan::math::inc_beta_inv; + using stan::math::inv_inc_beta; using stan::math::var; double a_d = 1; double b_d = 2; @@ -16,7 +16,7 @@ TEST(ProbInternalMath, inc_beta_inv_fv1) { b_v.d_ = 1.0; p_v.d_ = 1.0; - fvar res = inc_beta_inv(a_v, b_v, p_v); + fvar res = inv_inc_beta(a_v, b_v, p_v); res.val_.grad(); EXPECT_FLOAT_EQ(a_v.val_.adj(), 0.287698278597); @@ -30,7 +30,7 @@ TEST(ProbInternalMath, inc_beta_inv_fv1) { b_v.d_ = 1.0; p_v.d_ = 1.0; - res = inc_beta_inv(a_d, b_v, p_v); + res = inv_inc_beta(a_d, b_v, p_v); res.val_.grad(); EXPECT_FLOAT_EQ(b_v.val_.adj(), -0.122532267934); @@ -41,7 +41,7 @@ TEST(ProbInternalMath, inc_beta_inv_fv1) { b_v.d_ = 1.0; p_v.d_ = 1.0; - res = inc_beta_inv(a_v, b_d, p_v); + res = inv_inc_beta(a_v, b_d, p_v); res.val_.grad(); EXPECT_FLOAT_EQ(a_v.val_.adj(), 0.287698278597); @@ -52,16 +52,16 @@ TEST(ProbInternalMath, inc_beta_inv_fv1) { a_v.d_ = 1.0; p_v.d_ = 1.0; - res = inc_beta_inv(a_v, b_v, p_d); + res = inv_inc_beta(a_v, b_v, p_d); res.val_.grad(); EXPECT_FLOAT_EQ(a_v.val_.adj(), 0.287698278597); EXPECT_FLOAT_EQ(b_v.val_.adj(), -0.122532267934); } -TEST(ProbInternalMath, inc_beta_inv_fv2) { +TEST(ProbInternalMath, inv_inc_beta_fv2) { using stan::math::fvar; - using stan::math::inc_beta_inv; + using stan::math::inv_inc_beta; using stan::math::var; fvar> a = 2; fvar> b = 5; @@ -70,7 +70,7 @@ TEST(ProbInternalMath, inc_beta_inv_fv2) { b.d_ = 1.0; p.d_ = 1.0; - fvar> res = inc_beta_inv(a, b, p); + fvar> res = inv_inc_beta(a, b, p); res.val_.val_.grad(); EXPECT_FLOAT_EQ(a.val_.val_.adj(), 0.0783025374798); diff --git a/test/unit/math/prim/fun/inc_beta_inv_test.cpp b/test/unit/math/prim/fun/inc_beta_inv_test.cpp deleted file mode 100644 index 815e4ba0e91..00000000000 --- a/test/unit/math/prim/fun/inc_beta_inv_test.cpp +++ /dev/null @@ -1,70 +0,0 @@ -#include -#include -#include - -TEST(MathFunctions, inc_beta_inv_inv) { - using stan::math::inc_beta_inv; - - EXPECT_FLOAT_EQ(0.146446609407, inc_beta_inv(0.5, 0.5, 0.25)); - EXPECT_FLOAT_EQ(0.5, inc_beta_inv(0.5, 0.5, 0.5)); - EXPECT_FLOAT_EQ(0.853553390593, inc_beta_inv(0.5, 0.5, 0.75)); - EXPECT_FLOAT_EQ(1.0, inc_beta_inv(0.5, 0.5, 1.0)); - - EXPECT_FLOAT_EQ(0.0, inc_beta_inv(0.1, 1.5, 0.0)); - EXPECT_FLOAT_EQ(5.33624995327e-7, inc_beta_inv(0.1, 1.5, 0.25)); - EXPECT_FLOAT_EQ(0.000546567646391, inc_beta_inv(0.1, 1.5, 0.5)); - EXPECT_FLOAT_EQ(0.0319736161201, inc_beta_inv(0.1, 1.5, 0.75)); - EXPECT_FLOAT_EQ(1.0, inc_beta_inv(0.1, 1.5, 1.0)); - - EXPECT_FLOAT_EQ(0.804032164227, inc_beta_inv(0.6, 0.3, 0.5)); -} - -TEST(MathFunctions, inc_beta_inv_a_boundary) { - double b = 0.5; - double p = 0.5; - - const double inf = std::numeric_limits::infinity(); - - EXPECT_THROW(stan::math::inc_beta_inv(0.0, b, p), std::domain_error); - EXPECT_NO_THROW(stan::math::inc_beta_inv(inf, b, p)); - EXPECT_THROW(stan::math::inc_beta_inv(-0.01, b, p), std::domain_error); -} - -TEST(MathFunctions, inc_beta_inv_b_boundary) { - double a = 0.5; - double p = 0.5; - - const double inf = std::numeric_limits::infinity(); - - EXPECT_THROW(stan::math::inc_beta_inv(a, 0.0, p), std::domain_error); - EXPECT_NO_THROW(stan::math::inc_beta_inv(inf, 1.0, p)); - EXPECT_THROW(stan::math::inc_beta_inv(a, -0.01, p), std::domain_error); -} - -TEST(MathFunctions, inc_beta_inv_x_boundary) { - double a = 0.5; - double b = 0.5; - - EXPECT_NO_THROW(stan::math::inc_beta_inv(a, b, 0.0)); - EXPECT_NO_THROW(stan::math::inc_beta_inv(a, b, 1.0)); - EXPECT_THROW(stan::math::inc_beta_inv(a, b, -0.01), std::domain_error); - EXPECT_THROW(stan::math::inc_beta_inv(a, b, 1.01), std::domain_error); -} - -TEST(MathFunctions, inc_beta_inv_a_b_boundary) { - double p = 0.5; - - EXPECT_THROW(stan::math::inc_beta_inv(0.0, 0.0, p), std::domain_error); -} - -TEST(MathFunctions, inc_beta_inv_nan) { - double nan = std::numeric_limits::quiet_NaN(); - - EXPECT_THROW(stan::math::inc_beta_inv(0.5, 0.0, nan), std::domain_error); - EXPECT_THROW(stan::math::inc_beta_inv(0.5, nan, 0.0), std::domain_error); - EXPECT_THROW(stan::math::inc_beta_inv(0.5, nan, nan), std::domain_error); - EXPECT_THROW(stan::math::inc_beta_inv(nan, 0.0, 0.0), std::domain_error); - EXPECT_THROW(stan::math::inc_beta_inv(nan, 0.0, nan), std::domain_error); - EXPECT_THROW(stan::math::inc_beta_inv(nan, nan, 0.0), std::domain_error); - EXPECT_THROW(stan::math::inc_beta_inv(nan, nan, nan), std::domain_error); -} diff --git a/test/unit/math/prim/fun/inv_inc_beta_test.cpp b/test/unit/math/prim/fun/inv_inc_beta_test.cpp new file mode 100644 index 00000000000..9be2807ac9e --- /dev/null +++ b/test/unit/math/prim/fun/inv_inc_beta_test.cpp @@ -0,0 +1,70 @@ +#include +#include +#include + +TEST(MathFunctions, inv_inc_beta_inv) { + using stan::math::inv_inc_beta; + + EXPECT_FLOAT_EQ(0.146446609407, inv_inc_beta(0.5, 0.5, 0.25)); + EXPECT_FLOAT_EQ(0.5, inv_inc_beta(0.5, 0.5, 0.5)); + EXPECT_FLOAT_EQ(0.853553390593, inv_inc_beta(0.5, 0.5, 0.75)); + EXPECT_FLOAT_EQ(1.0, inv_inc_beta(0.5, 0.5, 1.0)); + + EXPECT_FLOAT_EQ(0.0, inv_inc_beta(0.1, 1.5, 0.0)); + EXPECT_FLOAT_EQ(5.33624995327e-7, inv_inc_beta(0.1, 1.5, 0.25)); + EXPECT_FLOAT_EQ(0.000546567646391, inv_inc_beta(0.1, 1.5, 0.5)); + EXPECT_FLOAT_EQ(0.0319736161201, inv_inc_beta(0.1, 1.5, 0.75)); + EXPECT_FLOAT_EQ(1.0, inv_inc_beta(0.1, 1.5, 1.0)); + + EXPECT_FLOAT_EQ(0.804032164227, inv_inc_beta(0.6, 0.3, 0.5)); +} + +TEST(MathFunctions, inv_inc_beta_a_boundary) { + double b = 0.5; + double p = 0.5; + + const double inf = std::numeric_limits::infinity(); + + EXPECT_THROW(stan::math::inv_inc_beta(0.0, b, p), std::domain_error); + EXPECT_NO_THROW(stan::math::inv_inc_beta(inf, b, p)); + EXPECT_THROW(stan::math::inv_inc_beta(-0.01, b, p), std::domain_error); +} + +TEST(MathFunctions, inv_inc_beta_b_boundary) { + double a = 0.5; + double p = 0.5; + + const double inf = std::numeric_limits::infinity(); + + EXPECT_THROW(stan::math::inv_inc_beta(a, 0.0, p), std::domain_error); + EXPECT_NO_THROW(stan::math::inv_inc_beta(inf, 1.0, p)); + EXPECT_THROW(stan::math::inv_inc_beta(a, -0.01, p), std::domain_error); +} + +TEST(MathFunctions, inv_inc_beta_x_boundary) { + double a = 0.5; + double b = 0.5; + + EXPECT_NO_THROW(stan::math::inv_inc_beta(a, b, 0.0)); + EXPECT_NO_THROW(stan::math::inv_inc_beta(a, b, 1.0)); + EXPECT_THROW(stan::math::inv_inc_beta(a, b, -0.01), std::domain_error); + EXPECT_THROW(stan::math::inv_inc_beta(a, b, 1.01), std::domain_error); +} + +TEST(MathFunctions, inv_inc_beta_a_b_boundary) { + double p = 0.5; + + EXPECT_THROW(stan::math::inv_inc_beta(0.0, 0.0, p), std::domain_error); +} + +TEST(MathFunctions, inv_inc_beta_nan) { + double nan = std::numeric_limits::quiet_NaN(); + + EXPECT_THROW(stan::math::inv_inc_beta(0.5, 0.0, nan), std::domain_error); + EXPECT_THROW(stan::math::inv_inc_beta(0.5, nan, 0.0), std::domain_error); + EXPECT_THROW(stan::math::inv_inc_beta(0.5, nan, nan), std::domain_error); + EXPECT_THROW(stan::math::inv_inc_beta(nan, 0.0, 0.0), std::domain_error); + EXPECT_THROW(stan::math::inv_inc_beta(nan, 0.0, nan), std::domain_error); + EXPECT_THROW(stan::math::inv_inc_beta(nan, nan, 0.0), std::domain_error); + EXPECT_THROW(stan::math::inv_inc_beta(nan, nan, nan), std::domain_error); +} From 5803eea074c0542a590439c32af81c4657f1f984 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 7 Jan 2022 10:55:03 +0800 Subject: [PATCH 09/12] Missed test naming --- .../fun/{inc_beta_inv_test.cpp => inv_inc_beta_test.cpp} | 8 ++++---- .../fun/{inc_beta_inv_test.cpp => inv_inc_beta_test.cpp} | 6 +++--- 2 files changed, 7 insertions(+), 7 deletions(-) rename test/unit/math/fwd/fun/{inc_beta_inv_test.cpp => inv_inc_beta_test.cpp} (79%) rename test/unit/math/rev/fun/{inc_beta_inv_test.cpp => inv_inc_beta_test.cpp} (78%) diff --git a/test/unit/math/fwd/fun/inc_beta_inv_test.cpp b/test/unit/math/fwd/fun/inv_inc_beta_test.cpp similarity index 79% rename from test/unit/math/fwd/fun/inc_beta_inv_test.cpp rename to test/unit/math/fwd/fun/inv_inc_beta_test.cpp index 9ff78cc9c1e..7caf4800ac8 100644 --- a/test/unit/math/fwd/fun/inc_beta_inv_test.cpp +++ b/test/unit/math/fwd/fun/inv_inc_beta_test.cpp @@ -3,7 +3,7 @@ TEST(AgradFwdMatrixIncBetaInv, fd_scalar) { using stan::math::fvar; - using stan::math::inc_beta_inv; + using stan::math::inv_inc_beta; fvar a = 6; fvar b = 2; fvar p = 0.9; @@ -11,14 +11,14 @@ TEST(AgradFwdMatrixIncBetaInv, fd_scalar) { b.d_ = 1.0; p.d_ = 1.0; - fvar res = inc_beta_inv(a, b, p); + fvar res = inv_inc_beta(a, b, p); EXPECT_FLOAT_EQ(res.d_, 0.0117172527399 - 0.0680999818473 + 0.455387298585); } TEST(AgradFwdMatrixIncBetaInv, ffd_scalar) { using stan::math::fvar; - using stan::math::inc_beta_inv; + using stan::math::inv_inc_beta; fvar> a = 7; fvar> b = 4; fvar> p = 0.15; @@ -26,7 +26,7 @@ TEST(AgradFwdMatrixIncBetaInv, ffd_scalar) { b.val_.d_ = 1.0; p.val_.d_ = 1.0; - fvar> res = inc_beta_inv(a, b, p); + fvar> res = inv_inc_beta(a, b, p); EXPECT_FLOAT_EQ(res.val_.d_, 0.0428905418857 - 0.0563420377808 + 0.664919819507); diff --git a/test/unit/math/rev/fun/inc_beta_inv_test.cpp b/test/unit/math/rev/fun/inv_inc_beta_test.cpp similarity index 78% rename from test/unit/math/rev/fun/inc_beta_inv_test.cpp rename to test/unit/math/rev/fun/inv_inc_beta_test.cpp index 389bb24b13a..9f485a78e9c 100644 --- a/test/unit/math/rev/fun/inc_beta_inv_test.cpp +++ b/test/unit/math/rev/fun/inv_inc_beta_test.cpp @@ -3,15 +3,15 @@ #include #include -TEST(inc_beta_inv, values) { - using stan::math::inc_beta_inv; +TEST(inv_inc_beta, values) { + using stan::math::inv_inc_beta; using stan::math::var; var a = 25; var b = 2; var p = 0.7; - var res = inc_beta_inv(a, b, p); + var res = inv_inc_beta(a, b, p); res.grad(); EXPECT_FLOAT_EQ(a.adj(), 0.00161775443606); From b33a57bd26f48ae2e6f1b81652a1cf820a0e300c Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 23 Mar 2022 20:19:38 +0800 Subject: [PATCH 10/12] Trigger CI From 488befa7b9b37c028e0979b73e8c2da884e35e80 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 25 Mar 2022 22:00:26 +0800 Subject: [PATCH 11/12] Add constraint checks --- stan/math/prim/fun/inv_inc_beta.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/stan/math/prim/fun/inv_inc_beta.hpp b/stan/math/prim/fun/inv_inc_beta.hpp index bc81bad8f97..ee3cf4157d7 100644 --- a/stan/math/prim/fun/inv_inc_beta.hpp +++ b/stan/math/prim/fun/inv_inc_beta.hpp @@ -26,6 +26,9 @@ inline double inv_inc_beta(double a, double b, double p) { check_not_nan("inv_inc_beta", "a", a); check_not_nan("inv_inc_beta", "b", b); check_not_nan("inv_inc_beta", "p", p); + check_positive("inv_inc_beta","a", a); + check_positive("inv_inc_beta","b", b); + check_bounded("inv_inc_beta", "p", p, 0, 1); return boost::math::ibeta_inv(a, b, p, boost_policy_t<>()); } From 7715e903cb0617f0d82bec87761c36e83ba9a58d Mon Sep 17 00:00:00 2001 From: Stan BuildBot Date: Fri, 25 Mar 2022 10:01:29 -0400 Subject: [PATCH 12/12] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/fun/inv_inc_beta.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/fun/inv_inc_beta.hpp b/stan/math/prim/fun/inv_inc_beta.hpp index ee3cf4157d7..920750b8112 100644 --- a/stan/math/prim/fun/inv_inc_beta.hpp +++ b/stan/math/prim/fun/inv_inc_beta.hpp @@ -26,8 +26,8 @@ inline double inv_inc_beta(double a, double b, double p) { check_not_nan("inv_inc_beta", "a", a); check_not_nan("inv_inc_beta", "b", b); check_not_nan("inv_inc_beta", "p", p); - check_positive("inv_inc_beta","a", a); - check_positive("inv_inc_beta","b", b); + check_positive("inv_inc_beta", "a", a); + check_positive("inv_inc_beta", "b", b); check_bounded("inv_inc_beta", "p", p, 0, 1); return boost::math::ibeta_inv(a, b, p, boost_policy_t<>()); }