From 511ab3ef90a841919291753970f9046b5834644f Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Mon, 17 Apr 2023 14:30:41 +0200 Subject: [PATCH] Fix for scipy issue 18302 --- include/boost/math/distributions/beta.hpp | 42 +++++++++++++++++-- include/boost/math/special_functions/beta.hpp | 11 ++++- test/Jamfile.v2 | 1 + test/scipy_issue_18302.cpp | 35 ++++++++++++++++ test/test_beta_dist.cpp | 8 ++-- 5 files changed, 88 insertions(+), 9 deletions(-) create mode 100644 test/scipy_issue_18302.cpp diff --git a/include/boost/math/distributions/beta.hpp b/include/boost/math/distributions/beta.hpp index 1de399ead2..10cbde8cc2 100644 --- a/include/boost/math/distributions/beta.hpp +++ b/include/boost/math/distributions/beta.hpp @@ -390,13 +390,47 @@ namespace boost } using boost::math::beta; - // Corner case: check_x ensures x element of [0, 1], but PDF is 0 for x = 0 and x = 1. PDF EQN: + // Corner cases: check_x ensures x element of [0, 1], but PDF is 0 for x = 0 and x = 1. PDF EQN: // https://wikimedia.org/api/rest_v1/media/math/render/svg/125fdaa41844a8703d1a8610ac00fbf3edacc8e7 - if(x == 0 || x == 1) + if(x == 0) { - return RealType(0); + if (a == 1) + { + return 1 / beta(a, b); + } + else if (a < 1) + { + policies::raise_overflow_error(function, nullptr, Policy()); + } + else + { + return RealType(0); + } + } + else if (x == 1) + { + if (b == 1) + { + return 1 / beta(a, b); + } + else if (b < 1) + { + policies::raise_overflow_error(function, nullptr, Policy()); + } + else + { + return RealType(0); + } + } + + // Bounds checking + RealType return_val = ibeta_derivative(a, b, x, Policy()); + if (return_val > b) + { + return_val = b; } - return ibeta_derivative(a, b, x, Policy()); + + return return_val; } // pdf template diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp index f0212800c7..77ca35dcd8 100644 --- a/include/boost/math/special_functions/beta.hpp +++ b/include/boost/math/special_functions/beta.hpp @@ -1452,7 +1452,16 @@ T ibeta_derivative_imp(T a, T b, T x, const Policy& pol) // typedef typename lanczos::lanczos::type lanczos_type; T y = (1 - x) * x; - T f1 = ibeta_power_terms(a, b, x, 1 - x, lanczos_type(), true, pol, 1 / y, function); + T f1; + if (!(boost::math::isinf)(1 / y)) + { + f1 = ibeta_power_terms(a, b, x, 1 - x, lanczos_type(), true, pol, 1 / y, function); + } + else + { + return (a > 1) ? 0 : (a == 1) ? 1 / boost::math::beta(a, b, pol) : policies::raise_overflow_error(function, nullptr, pol); + } + return f1; } // diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 1608aa2993..48c3ff21e3 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -899,6 +899,7 @@ test-suite distribution_tests : [ run scipy_issue_17388.cpp ../../test/build//boost_unit_test_framework ] [ run scipy_issue_17916.cpp ../../test/build//boost_unit_test_framework ] [ run scipy_issue_17916_nct.cpp ../../test/build//boost_unit_test_framework ] + [ run scipy_issue_18302.cpp ../../test/build//boost_unit_test_framework ] ; test-suite mp : diff --git a/test/scipy_issue_18302.cpp b/test/scipy_issue_18302.cpp new file mode 100644 index 0000000000..739a3047a4 --- /dev/null +++ b/test/scipy_issue_18302.cpp @@ -0,0 +1,35 @@ +// Copyright Matt Borland, 2023 +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE or copy at http://www.boost.org/LICENSE_1_0.txt) +// +// See: https://github.com/scipy/scipy/issues/17916 + +#include +#include +#include +#include "math_unit_test.hpp" + +int main(void) +{ + using my_policy = boost::math::policies::policy>; + + auto dist = boost::math::beta_distribution(1, 5); + + // https://www.wolframalpha.com/input?i=PDF%28beta+distribution%281%2C+5%29%2C+0%29 + double test_pdf_spot = boost::math::pdf(dist, 0); + CHECK_ULP_CLOSE(test_pdf_spot, 5.0, 1); + + // https://www.wolframalpha.com/input?i=PDF%28beta+distribution%281%2C+5%29%2C+1e-30%29 + test_pdf_spot = boost::math::pdf(dist, 1e-30); + CHECK_ULP_CLOSE(test_pdf_spot, 5.0, 1); + + // Appox equal to 5 + test_pdf_spot = boost::math::pdf(dist, 1e-310); + CHECK_ULP_CLOSE(test_pdf_spot, 5.0, 1); + + test_pdf_spot = boost::math::pdf(dist, 1); + CHECK_ULP_CLOSE(test_pdf_spot, 0.0, 1); + + return boost::math::test::report_errors(); +} diff --git a/test/test_beta_dist.cpp b/test/test_beta_dist.cpp index a0041a4b1b..d441e5b9dd 100644 --- a/test/test_beta_dist.cpp +++ b/test/test_beta_dist.cpp @@ -197,11 +197,11 @@ void test_spots(RealType) BOOST_CHECK_EQUAL( // a = b = 1 is uniform distribution. pdf(beta_distribution(static_cast(1), static_cast(1)), static_cast(1)), // x - static_cast(0)); + static_cast(1)); BOOST_CHECK_EQUAL( pdf(beta_distribution(static_cast(1), static_cast(1)), static_cast(0)), // x - static_cast(0)); + static_cast(1)); BOOST_CHECK_CLOSE_FRACTION( pdf(beta_distribution(static_cast(1), static_cast(1)), static_cast(0.5)), // x @@ -561,8 +561,8 @@ BOOST_AUTO_TEST_CASE( test_main ) beta_distribution<> mybetaH3(0.5, 3.); // // Check a few values using double. - BOOST_CHECK_EQUAL(pdf(mybeta11, 1), 0); // is uniform unity over (0, 1) - BOOST_CHECK_EQUAL(pdf(mybeta11, 0), 0); // https://www.wolframalpha.com/input/?i=beta+distribution+pdf+alpha+%3D+1%2C+beta+%3D+1 + BOOST_CHECK_EQUAL(pdf(mybeta11, 1), 1); // is uniform unity over (0, 1) + BOOST_CHECK_EQUAL(pdf(mybeta11, 0), 1); // Although these next three have an exact result, internally they're // *not* treated as special cases, and may be out by a couple of eps: BOOST_CHECK_CLOSE_FRACTION(pdf(mybeta11, 0.5), 1.0, 5*std::numeric_limits::epsilon());