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..6d4db2166a 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) @@ -542,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);