From b4848f9125de9f97ede01b7aad9fef4dc7ecbca9 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Thu, 13 Aug 2026 17:31:41 +0100 Subject: [PATCH 1/2] Non Central T: track cancellation errors, And switch to asymptotic approximation when all else fails. Fixes: https://github.com/boostorg/math/issues/1430 --- .../math/distributions/non_central_t.hpp | 21 +++++++++++++++---- test/test_nc_t.hpp | 5 +++++ 2 files changed, 22 insertions(+), 4 deletions(-) diff --git a/include/boost/math/distributions/non_central_t.hpp b/include/boost/math/distributions/non_central_t.hpp index 72d1d1678d..dd873fa0cd 100644 --- a/include/boost/math/distributions/non_central_t.hpp +++ b/include/boost/math/distributions/non_central_t.hpp @@ -249,8 +249,7 @@ namespace boost // t-Distribution". C. van Eeden. International Statistical Review, 29, 4-31. // "Continuous Univariate Distributions". N.L. Johnson, S. Kotz and // N. Balkrishnan. 1995. John Wiley and Sons New York. - T result = cdf(students_t_distribution(v), t - delta); - return invert ? 1 - result : result; + return invert ? cdf(complement(students_t_distribution(v), t - delta)) : cdf(students_t_distribution(v), t - delta); } // // x and y are the corresponding random @@ -274,10 +273,11 @@ namespace boost // // Calculate p: // + T c = 0; if(x != 0) { - result = non_central_beta_p(a, b, d2, x, y, pol); - result = non_central_t2_p(v, delta, x, y, pol, result); + c = non_central_beta_p(a, b, d2, x, y, pol); + result = non_central_t2_p(v, delta, x, y, pol, c); result /= 2; } else @@ -285,6 +285,19 @@ namespace boost if (invert) { result = cdf(complement(boost::math::normal_distribution(), -delta)) - result; + if ((x != 0) && (fabs(result / (c * tools::epsilon())) < 1000)) + { + // We've cancelled out most of the digits in the result, try A&S 26.7.9, + // this is only accurate to a couple of digits at best, but it's better + // than nothing when all else fails, see https://github.com/boostorg/math/issues/1430 + t = -t; + delta = -delta; + + T z = (t * (1 - 1 / (4 * v)) - delta) / sqrt(1 + t * t / (2 * v)); + + return cdf(boost::math::normal_distribution(), z); + + } invert = false; } else diff --git a/test/test_nc_t.hpp b/test/test_nc_t.hpp index 83dbf12c99..1950a008c2 100644 --- a/test/test_nc_t.hpp +++ b/test/test_nc_t.hpp @@ -344,6 +344,11 @@ void test_spots(RealType) distro1 d(8.0f, 8.5f); BOOST_CHECK_CLOSE(pdf(d, -1), static_cast(6.1747948083757028903541988987716621647020752431287e-20), 2e-5); // Can we do better on accuracy here? } + // https://github.com/boostorg/math/issues/1430 + { + distro1 d(1000.f, 23.f); + BOOST_CHECK_CLOSE_FRACTION(cdf(d, -1), static_cast(1.61471461239552e-127), 1e-3); + } } // template void test_spots(RealType) From 811faf72059d34dff98fea8d50c20678df46efda Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Thu, 13 Aug 2026 20:47:56 +0100 Subject: [PATCH 2/2] Correct unrelated failures: Use of find_degrees_of_freedom requires no internal promotion if we are to get the expected error. --- test/test_nc_t.hpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/test/test_nc_t.hpp b/test/test_nc_t.hpp index 1950a008c2..6d4db2166a 100644 --- a/test/test_nc_t.hpp +++ b/test/test_nc_t.hpp @@ -547,8 +547,10 @@ void quantile_sanity_check(T& data, const char* type_name, const char* test) BOOST_ERROR(e.what()); } // Code coverage: - BOOST_CHECK_THROW(boost::math::non_central_t_distribution::find_degrees_of_freedom(data[i][1], boost::math::tools::epsilon() / 2, data[i][3]), boost::math::evaluation_error); - BOOST_CHECK_THROW(boost::math::non_central_t_distribution::find_degrees_of_freedom(boost::math::complement(data[i][1], boost::math::tools::epsilon() / 2, data[i][3])), boost::math::evaluation_error); + using no_promote_policy = boost::math::policies::policy, boost::math::policies::promote_double >; + using no_promote_distro = boost::math::non_central_t_distribution; + BOOST_CHECK_THROW(no_promote_distro::find_degrees_of_freedom(data[i][1], boost::math::tools::epsilon() / 2, data[i][3]), boost::math::evaluation_error); + BOOST_CHECK_THROW(no_promote_distro::find_degrees_of_freedom(boost::math::complement(data[i][1], boost::math::tools::epsilon() / 2, data[i][3])), boost::math::evaluation_error); BOOST_CHECK_THROW(boost::math::non_central_t_distribution::find_degrees_of_freedom(data[i][1], data[i][2], value_type(0)), boost::math::evaluation_error); BOOST_CHECK_THROW(boost::math::non_central_t_distribution::find_degrees_of_freedom(boost::math::complement(data[i][1], data[i][2], value_type(0))), boost::math::evaluation_error); BOOST_CHECK_THROW(boost::math::non_central_t_distribution::find_degrees_of_freedom(data[i][1], data[i][2], value_type(1)), boost::math::evaluation_error);