From c5adf11ec58b008052251bf408e8e4b631e114fc Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Tue, 1 Sep 2026 08:21:11 -0700 Subject: [PATCH] Fix skew normal by exp-sinh quadrature Closes #1190 --- .../boost/math/distributions/skew_normal.hpp | 58 ++++++++++++++++--- test/test_skew_normal.cpp | 22 +++++++ 2 files changed, 72 insertions(+), 8 deletions(-) diff --git a/include/boost/math/distributions/skew_normal.hpp b/include/boost/math/distributions/skew_normal.hpp index c4e429ed2f..6878fdc6bf 100644 --- a/include/boost/math/distributions/skew_normal.hpp +++ b/include/boost/math/distributions/skew_normal.hpp @@ -17,6 +17,7 @@ #include // Owen's T function #include #include +#include #include #include #include @@ -55,6 +56,22 @@ namespace boost{ namespace math{ return true; } +#ifndef BOOST_MATH_HAS_NVRTC + template + inline RealType skew_normal_tail_integral(RealType x, RealType shape, bool upper) + { + normal_distribution std_normal; + quadrature::exp_sinh integrator; + const RealType direction = upper ? static_cast(1) : static_cast(-1); + const auto integrand = [&](RealType t)->RealType + { + const RealType z = x + direction * t; + return static_cast(2) * pdf(std_normal, z) * cdf(std_normal, shape * z); + }; + return integrator.integrate(integrand, policies::get_epsilon() * 8); + } +#endif + } // namespace detail BOOST_MATH_EXPORT template > @@ -218,7 +235,16 @@ namespace boost{ namespace math{ normal_distribution std_normal; - result = cdf(std_normal, transformed_x) - owens_t(transformed_x, shape)*static_cast(2); + const RealType normal_cdf = cdf(std_normal, transformed_x); + result = normal_cdf - owens_t(transformed_x, shape)*static_cast(2); + +#ifndef BOOST_MATH_HAS_NVRTC + if((shape > 0) && (transformed_x < 0) + && (result < normal_cdf * boost::math::tools::root_epsilon())) + { + result = detail::skew_normal_tail_integral(transformed_x, shape, false); + } +#endif return result; } // cdf @@ -261,7 +287,16 @@ namespace boost{ namespace math{ normal_distribution std_normal; - result = cdf(complement(std_normal, transformed_x)) + owens_t(transformed_x, shape)*static_cast(2); + const RealType normal_cdf = cdf(complement(std_normal, transformed_x)); + result = normal_cdf + owens_t(transformed_x, shape)*static_cast(2); + +#ifndef BOOST_MATH_HAS_NVRTC + if((shape < 0) && (transformed_x > 0) + && (result < normal_cdf * boost::math::tools::root_epsilon())) + { + result = detail::skew_normal_tail_integral(transformed_x, shape, true); + } +#endif return result; } // cdf complement @@ -653,11 +688,16 @@ namespace boost{ namespace math{ - x*(static_cast(2)*x*x-static_cast(5))*skew*skew/static_cast(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 standard_dist( + static_cast(0), static_cast(1), shape); + result = standard_deviation(standard_dist)*x+mean(standard_dist); // refine the result by numerically searching the root of (p-cdf) @@ -667,12 +707,12 @@ namespace boost{ namespace math{ if (result == 0) result = tools::min_value(); // 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 +744,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 +753,8 @@ namespace boost{ namespace math{ else result -= step; + result = location + scale * result; + if (max_iter >= policies::get_max_root_iterations()) { return policies::raise_evaluation_error(function, "Unable to locate solution in a reasonable time: either there is no answer to quantile" // LCOV_EXCL_LINE diff --git a/test/test_skew_normal.cpp b/test/test_skew_normal.cpp index 617a8f6fab..58780cd0e5 100644 --- a/test/test_skew_normal.cpp +++ b/test/test_skew_normal.cpp @@ -496,6 +496,28 @@ BOOST_AUTO_TEST_CASE( test_main ) BOOST_CHECK_CLOSE_FRACTION(mean(w01), static_cast(0), tolfeweps); // Default mean == zero BOOST_CHECK_CLOSE_FRACTION(scale(w01), static_cast(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::epsilon(); + boost::math::normal_distribution 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(0, 1, 1), -8.0), + alpha_one_tail, + tail_tolerance); + BOOST_CHECK_CLOSE_FRACTION( + cdf(complement(skew_normal_distribution(0, 1, -1), 8.0)), + alpha_one_tail, + tail_tolerance); + BOOST_CHECK_CLOSE_FRACTION( + cdf(skew_normal_distribution(0, 1, 2), -6.0), + 7.1180791906932412294852794865326109874556091859546e-43, + tail_tolerance); + } + // 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 %