Fix skew normal by exp-sinh quadrature

include/boost/math/distributions/skew_normal.hpp

@@ -17,6 +17,7 @@
 #include <boost/math/special_functions/owens_t.hpp> // Owen's T function
 #include <boost/math/distributions/complement.hpp>
 #include <boost/math/distributions/normal.hpp>
+#include <boost/math/quadrature/exp_sinh.hpp>
 #include <boost/math/distributions/detail/common_error_handling.hpp>
 #include <boost/math/constants/constants.hpp>
 #include <boost/math/tools/tuple.hpp>
@@ -55,6 +56,22 @@ namespace boost{ namespace math{
       return true;
     }
 
+#ifndef BOOST_MATH_HAS_NVRTC
+    template <class RealType, class Policy>
+    inline RealType skew_normal_tail_integral(RealType x, RealType shape, bool upper)
+    {
+      normal_distribution<RealType, Policy> std_normal;
+      quadrature::exp_sinh<RealType, Policy> integrator;
+      const RealType direction = upper ? static_cast<RealType>(1) : static_cast<RealType>(-1);
+      const auto integrand = [&](RealType t)->RealType
+      {
+        const RealType z = x + direction * t;
+        return static_cast<RealType>(2) * pdf(std_normal, z) * cdf(std_normal, shape * z);
+      };
+      return integrator.integrate(integrand, policies::get_epsilon<RealType, Policy>() * 8);
+    }
+#endif
+
   } // namespace detail
 
   BOOST_MATH_EXPORT template <class RealType = double, class Policy = policies::policy<> >
@@ -218,7 +235,19 @@ namespace boost{ namespace math{
 
     normal_distribution<RealType, Policy> std_normal;
 
-    result = cdf(std_normal, transformed_x) - owens_t(transformed_x, shape)*static_cast<RealType>(2);
+    const RealType normal_cdf = cdf(std_normal, transformed_x);
+    result = normal_cdf - owens_t(transformed_x, shape)*static_cast<RealType>(2);
+
+#ifndef BOOST_MATH_HAS_NVRTC
+    // Subtraction magnifies the error in the normal CDF and Owen's T.
+    // Switch after losing three bits, rather than waiting until half the
+    // precision is lost: even moderate cancellation can spoil quantiles.
+    if((shape > 0) && (transformed_x < 0)
+      && (result < normal_cdf / static_cast<RealType>(8)))
+    {
+      result = detail::skew_normal_tail_integral<RealType, Policy>(transformed_x, shape, false);
+    }
+#endif
 
     return result;
   } // cdf
@@ -261,7 +290,17 @@ namespace boost{ namespace math{
 
     normal_distribution<RealType, Policy> std_normal;
 
-    result = cdf(complement(std_normal, transformed_x)) + owens_t(transformed_x, shape)*static_cast<RealType>(2);
+    const RealType normal_cdf = cdf(complement(std_normal, transformed_x));
+    result = normal_cdf + owens_t(transformed_x, shape)*static_cast<RealType>(2);
+
+#ifndef BOOST_MATH_HAS_NVRTC
+    // The reflected tail has the same cancellation as the lower CDF.
+    if((shape < 0) && (transformed_x > 0)
+      && (result < normal_cdf / static_cast<RealType>(8)))
+    {
+      result = detail::skew_normal_tail_integral<RealType, Policy>(transformed_x, shape, true);
+    }
+#endif
     return result;
   } // cdf complement
 
@@ -653,11 +692,16 @@ namespace boost{ namespace math{
       - x*(static_cast<RealType>(2)*x*x-static_cast<RealType>(5))*skew*skew/static_cast<RealType>(36);
     } // if(shape != 0)
 
-    result = standard_deviation(dist)*x+mean(dist);
-
     // handle special case of non-skew normal distribution.
     if(shape == 0)
-      return result;
+      return standard_deviation(dist)*x+mean(dist);
+
+    // Search in standardized coordinates.  bracket_and_solve_root expands
+    // multiplicatively about zero, so searching in the user's location/scale
+    // can turn a good initial estimate into a very wide bracket.
+    skew_normal_distribution<RealType, Policy> standard_dist(
+       static_cast<RealType>(0), static_cast<RealType>(1), shape);
+    result = standard_deviation(standard_dist)*x+mean(standard_dist);
 
     // refine the result by numerically searching the root of (p-cdf)
 
@@ -667,12 +711,12 @@ namespace boost{ namespace math{
     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; };
+    auto fun = [&, standard_dist, p](const RealType& x)->RealType { return cdf(standard_dist, x) - p; };
 
     RealType f_result = fun(result);
 
     if (f_result == 0)
-       return result;
+       return location + scale * result;
 
     if (f_result * result > 0)
     {
@@ -704,7 +748,7 @@ namespace boost{ namespace math{
     //
     // Try one last Newton step, just to close up the interval:
     //
-    RealType step = fun(result) / pdf(dist, result);
+    RealType step = fun(result) / pdf(standard_dist, result);
 
     if (result - step <= p_result.first)
        result = p_result.first;
@@ -713,6 +757,8 @@ namespace boost{ namespace math{
     else
        result -= step;
 
+    result = location + scale * result;
+
     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/test_skew_normal.cpp

@@ -496,6 +496,61 @@ BOOST_AUTO_TEST_CASE( test_main )
   BOOST_CHECK_CLOSE_FRACTION(mean(w01), static_cast<double>(0), tolfeweps); // Default mean == zero
   BOOST_CHECK_CLOSE_FRACTION(scale(w01), static_cast<double>(1), tolfeweps); // Default scale == unity
 
+  // https://github.com/boostorg/math/issues/1190
+  // Avoid cancellation in the extreme left tail (and its reflected upper tail).
+  {
+    const double tail_tolerance = 128 * numeric_limits<double>::epsilon();
+    boost::math::normal_distribution<double> std_normal;
+    const double normal_tail = cdf(std_normal, -8.0);
+    const double alpha_one_tail = normal_tail * normal_tail;
+
+    BOOST_CHECK_CLOSE_FRACTION(
+      cdf(skew_normal_distribution<double>(0, 1, 1), -8.0),
+      alpha_one_tail,
+      tail_tolerance);
+    BOOST_CHECK_CLOSE_FRACTION(
+      cdf(complement(skew_normal_distribution<double>(0, 1, -1), 8.0)),
+      alpha_one_tail,
+      tail_tolerance);
+    BOOST_CHECK_CLOSE_FRACTION(
+      cdf(skew_normal_distribution<double>(0, 1, 2), -6.0),
+      7.1180791906932412294852794865326109874556091859546e-43,
+      tail_tolerance);
+  }
+
+  // Moderate cancellation must also use the accurate tail calculation.
+  // Disable promotion so extended long double cannot mask the loss of bits.
+  {
+    typedef boost::math::policies::policy<
+      boost::math::policies::promote_double<false> > no_promote_policy;
+    typedef skew_normal_distribution<double, no_promote_policy> distribution;
+    const double tolerance = 32 * numeric_limits<double>::epsilon();
+    const distribution lower(0, 1, 4), upper(0, 1, -4);
+    // Normal CDF - 2 * Owen's T, evaluated with cpp_bin_float_100.
+    const double expected = 8.179690339064550124947048052432943616e-7;
+    BOOST_CHECK_CLOSE_FRACTION(cdf(lower, -1.0), expected, tolerance);
+    BOOST_CHECK_CLOSE_FRACTION(cdf(complement(upper, 1.0)), expected, tolerance);
+
+    // Regression from git_issue_184: increasing probabilities must not give
+    // decreasing quantiles after the location/scale transformation.
+    const distribution shifted(573.39724735636185, 77.0, 4.0);
+    const distribution reflected(-573.39724735636185, 77.0, -4.0);
+    const double probabilities[] = {
+      0.00285612015554148, 0.00285612015554149, 0.00285612015554150
+    };
+    double previous = quantile(shifted, probabilities[0]);
+    double previous_complement = quantile(complement(reflected, probabilities[0]));
+    for (unsigned i = 1; i < 3; ++i)
+    {
+      const double current = quantile(shifted, probabilities[i]);
+      const double current_complement = quantile(complement(reflected, probabilities[i]));
+      BOOST_CHECK_LE(previous, current);
+      BOOST_CHECK_GE(previous_complement, current_complement);
+      previous = current;
+      previous_complement = current_complement;
+    }
+  }
+
   // Basic sanity-check spot values for all floating-point types..
   // (Parameter value, arbitrarily zero, only communicates the floating point type).
   test_spots(0.0F); // Test float. OK at decdigits = 0 tolerance = 0.0001 %