diff --git a/include/boost/math/distributions/skew_normal.hpp b/include/boost/math/distributions/skew_normal.hpp index c4e429ed2f..a42cbd5841 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,19 @@ 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 + // 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(8))) + { + result = detail::skew_normal_tail_integral(transformed_x, shape, false); + } +#endif return result; } // cdf @@ -261,7 +290,17 @@ 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 + // The reflected tail has the same cancellation as the lower CDF. + if((shape < 0) && (transformed_x > 0) + && (result < normal_cdf / static_cast(8))) + { + result = detail::skew_normal_tail_integral(transformed_x, shape, true); + } +#endif return result; } // cdf complement @@ -653,11 +692,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 +711,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 +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()) { 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..1697f544cd 100644 --- a/test/test_skew_normal.cpp +++ b/test/test_skew_normal.cpp @@ -496,6 +496,61 @@ 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); + } + + // 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 > no_promote_policy; + typedef skew_normal_distribution distribution; + const double tolerance = 32 * numeric_limits::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 %