Merge pull request #1123 from boostorg/issue1120

Change skew normal quantile to use bracket_and_solve_root.

Author: jzmaddock

include/boost/math/distributions/skew_normal.hpp

@@ -27,6 +27,10 @@
 #include <utility>
 #include <algorithm> // std::lower_bound, std::distance
 
+#ifdef BOOST_MATH_INSTRUMENT_SKEW_NORMAL_ITERATIONS
+extern std::uintmax_t global_iter_count;
+#endif
+
 namespace boost{ namespace math{
 
   namespace detail
@@ -614,32 +618,6 @@ namespace boost{ namespace math{
     return factor * y*y / (x*x);
   }
 
-  namespace detail
-  {
-
-    template <class RealType, class Policy>
-    struct skew_normal_quantile_functor
-    {
-      skew_normal_quantile_functor(const boost::math::skew_normal_distribution<RealType, Policy> dist, RealType const& p)
-        : distribution(dist), prob(p)
-      {
-      }
-
-      boost::math::tuple<RealType, RealType> operator()(RealType const& x)
-      {
-        RealType c = cdf(distribution, x);
-        RealType fx = c - prob;  // Difference cdf - value - to minimize.
-        RealType dx = pdf(distribution, x); // pdf is 1st derivative.
-        // return both function evaluation difference f(x) and 1st derivative f'(x).
-        return boost::math::make_tuple(fx, dx);
-      }
-    private:
-      const boost::math::skew_normal_distribution<RealType, Policy> distribution;
-      RealType prob;
-    };
-
-  } // namespace detail
-
   template <class RealType, class Policy>
   inline RealType quantile(const skew_normal_distribution<RealType, Policy>& dist, const RealType& p)
   {
@@ -681,14 +659,58 @@ namespace boost{ namespace math{
 
     // refine the result by numerically searching the root of (p-cdf)
 
-    const RealType search_min = support(dist).first;
-    const RealType search_max = support(dist).second;
-
     const int get_digits = policies::digits<RealType, Policy>();// get digits from policy,
     std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); // and max iterations.
 
-    result = tools::newton_raphson_iterate(detail::skew_normal_quantile_functor<RealType, Policy>(dist, p), result,
-      search_min, search_max, get_digits, max_iter);
+    if (result == 0)
+       result = tools::min_value<RealType>(); // we need to be one side of zero or the other for the root finder to work.
+
+    auto fun = [&, dist, p](const RealType& x)->RealType { return cdf(dist, x) - p; };
+
+    RealType f_result = fun(result);
+
+    if (f_result == 0)
+       return result;
+
+    if (f_result * result > 0)
+    {
+       // If the root is in the direction of zero, we need to check that we're the correct side of it:
+       RealType f_zero = fun(static_cast<RealType>(0));
+       if (f_zero * f_result > 0)
+       {
+          // we're the wrong side of zero:
+          result = -result;
+          f_result = fun(result);
+       }
+    }
+
+    RealType scaling_factor = 1.25;
+    if (f_result * result > 0)
+    {
+       // We're heading towards zero... it's a long way down so use a larger scaling factor:
+       scaling_factor = 16;
+    }
+
+    auto p_result = tools::bracket_and_solve_root(fun, result, scaling_factor, true, tools::eps_tolerance<RealType>(get_digits), max_iter, Policy());
+
+#ifdef BOOST_MATH_INSTRUMENT_SKEW_NORMAL_ITERATIONS
+    global_iter_count += max_iter;
+#endif
+
+    result = (p_result.first + p_result.second) / 2;
+
+    //
+    // Try one last Newton step, just to close up the interval:
+    //
+    RealType step = fun(result) / pdf(dist, result);
+
+    if (result - step <= p_result.first)
+       result = p_result.first;
+    else if (result - step >= p_result.second)
+       result = p_result.second;
+    else
+       result -= step;
+
     if (max_iter >= policies::get_max_root_iterations<Policy>())
     {
        return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time: either there is no answer to quantile" // LCOV_EXCL_LINE

test/Jamfile.v2

@@ -874,6 +874,7 @@ test-suite distribution_tests :
    [ run scipy_issue_18302.cpp ../../test/build//boost_unit_test_framework  ]
    [ run scipy_issue_18511.cpp ../../test/build//boost_unit_test_framework ]
    [ compile scipy_issue_19762.cpp ]
+   [ run git_issue_1120.cpp ]
 ;
 
 test-suite new_floats : 

test/git_issue_1120.cpp

@@ -0,0 +1,60 @@
+//  Copyright John Maddock 2024.
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+//
+// The purpose of this test case is to probe the skew normal quantiles
+// for extreme values of skewness and ensure that our root finders don't
+// blow up, see https://github.com/boostorg/math/issues/1120 for original
+// test case.  We test both the maximum number of iterations taken, and the
+// overall total (ie average).  Any changes to the skew normal should
+// ideally NOT cause this test to fail, as this indicates that our root
+// finding has been made worse by the change!!
+//
+// Note that defining BOOST_MATH_INSTRUMENT_SKEW_NORMAL_ITERATIONS
+// causes the skew normal quantile to save the number of iterations
+// to a global variable "global_iter_count".
+//
+
+#define BOOST_MATH_INSTRUMENT_SKEW_NORMAL_ITERATIONS
+
+#include <random>
+#include <boost/math/distributions/skew_normal.hpp>
+#include "math_unit_test.hpp"
+
+std::uintmax_t global_iter_count;
+std::uintmax_t total_iter_count = 0;
+
+int main()
+{
+   using scipy_policy = boost::math::policies::policy<boost::math::policies::promote_double<false>>;
+
+   std::mt19937 gen;
+   std::uniform_real_distribution<double> location(-3, 3);
+   std::uniform_real_distribution<double> scale(0.001, 3);
+
+   for (unsigned skew = 50; skew < 2000; skew += 43)
+   {
+      constexpr double pn[] = { 0.0001, 0.001, 0.01, 0.1, 0.2, 0.3, 0.4 };
+      boost::math::skew_normal_distribution<double, scipy_policy> dist(location(gen), scale(gen), skew);
+      for (unsigned i = 0; i < 7; ++i)
+      {
+         global_iter_count = 0;
+         double x = quantile(dist, pn[i]);
+         total_iter_count += global_iter_count;
+         CHECK_LE(global_iter_count, static_cast<std::uintmax_t>(60));
+         double p = cdf(dist, x);
+         CHECK_ABSOLUTE_ERROR(p, pn[i], 45 * std::numeric_limits<double>::epsilon());
+
+         global_iter_count = 0;
+         x = quantile(complement(dist, pn[i]));
+         total_iter_count += global_iter_count;
+         CHECK_LE(global_iter_count, static_cast<std::uintmax_t>(60));
+         p = cdf(complement(dist, x));
+         CHECK_ABSOLUTE_ERROR(p, pn[i], 45 * std::numeric_limits<double>::epsilon());
+      }
+   }
+   CHECK_LE(total_iter_count, static_cast<std::uintmax_t>(10000));
+   return boost::math::test::report_errors();
+}