Fix for scipy issue 18302

boostorg · Apr 17, 2023 · 511ab3e · 511ab3e
1 parent 109a814
commit 511ab3e
Show file tree

Hide file tree

Showing 5 changed files with 88 additions and 9 deletions.
diff --git 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<RealType>(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<RealType>(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 <class RealType, class Policy>

diff --git 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<T, Policy>::type lanczos_type;
    T y = (1 - x) * x;
-   T f1 = ibeta_power_terms<T>(a, b, x, 1 - x, lanczos_type(), true, pol, 1 / y, function);
+   T f1;
+   if (!(boost::math::isinf)(1 / y))
+   {
+      f1 = ibeta_power_terms<T>(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<T>(function, nullptr, pol);
+   }
+
    return f1;
 }
 //

diff --git 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
@@ -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 <boost/math/distributions/beta.hpp>
+#include <boost/math/policies/policy.hpp>
+#include <limits>
+#include "math_unit_test.hpp"
+
+int main(void)
+{
+    using my_policy = boost::math::policies::policy<boost::math::policies::promote_double<false>>;
+
+    auto dist = boost::math::beta_distribution<double, my_policy>(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
@@ -197,11 +197,11 @@ void test_spots(RealType)
   BOOST_CHECK_EQUAL( // a = b = 1 is uniform distribution.
      pdf(beta_distribution<RealType>(static_cast<RealType>(1), static_cast<RealType>(1)),
      static_cast<RealType>(1)),  // x
-     static_cast<RealType>(0));
+     static_cast<RealType>(1));
   BOOST_CHECK_EQUAL(
      pdf(beta_distribution<RealType>(static_cast<RealType>(1), static_cast<RealType>(1)),
      static_cast<RealType>(0)),  // x
-     static_cast<RealType>(0));
+     static_cast<RealType>(1));
   BOOST_CHECK_CLOSE_FRACTION(
      pdf(beta_distribution<RealType>(static_cast<RealType>(1), static_cast<RealType>(1)),
      static_cast<RealType>(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<double>::epsilon());