diff --git a/include/boost/math/distributions/non_central_t.hpp b/include/boost/math/distributions/non_central_t.hpp index 85581a4c80..c0c898387a 100644 --- a/include/boost/math/distributions/non_central_t.hpp +++ b/include/boost/math/distributions/non_central_t.hpp @@ -719,11 +719,6 @@ namespace boost return result; } -#if 0 - // - // This code is disabled, since there can be multiple answers to the - // question, and it's not clear how to find the "right" one. - // template struct t_degrees_of_freedom_finder { @@ -749,15 +744,21 @@ namespace boost inline RealType find_t_degrees_of_freedom( RealType delta, RealType x, RealType p, RealType q, const Policy& pol) { + using std::fabs; const char* function = "non_central_t<%1%>::find_degrees_of_freedom"; if((p == 0) || (q == 0)) { // - // Can't a thing if one of p and q is zero: + // Can't find a thing if one of p and q is zero: // return policies::raise_evaluation_error(function, "Can't find degrees of freedom when the probability is 0 or 1, only possible answer is %1%", // LCOV_EXCL_LINE RealType(std::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE } + if (fabs(x) < tools::epsilon()) + { + return policies::raise_evaluation_error(function, "Can't find degrees of freedom when the abscissa value is very close to zero as all degrees of freedom generate the same CDF at x=0: try again further out in the tails!!", + RealType(std::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE + } t_degrees_of_freedom_finder f(delta, x, p < q ? p : q, p < q ? false : true); tools::eps_tolerance tol(policies::digits()); std::uintmax_t max_iter = policies::get_max_root_iterations(); @@ -766,7 +767,7 @@ namespace boost // RealType guess = 200; std::pair ir = tools::bracket_and_solve_root( - f, guess, RealType(2), false, tol, max_iter, pol); + f, guess, RealType(2), x < 0 ? false : true, tol, max_iter, pol); RealType result = ir.first + (ir.second - ir.first) / 2; if(max_iter >= policies::get_max_root_iterations()) { @@ -819,9 +820,9 @@ namespace boost // RealType guess; if(f(0) < 0) - guess = 1; - else guess = -1; + else + guess = 1; std::pair ir = tools::bracket_and_solve_root( f, guess, RealType(2), false, tol, max_iter, pol); RealType result = ir.first + (ir.second - ir.first) / 2; @@ -832,7 +833,6 @@ namespace boost } return result; } -#endif } // namespace detail ====================================================================== template > @@ -864,26 +864,22 @@ namespace boost { // Private data getter function. return ncp; } -#if 0 - // - // This code is disabled, since there can be multiple answers to the - // question, and it's not clear how to find the "right" one. - // + static RealType find_degrees_of_freedom(RealType delta, RealType x, RealType p) { const char* function = "non_central_t<%1%>::find_degrees_of_freedom"; - typedef typename policies::evaluation::type value_type; + typedef typename policies::evaluation::type eval_type; typedef typename policies::normalise< Policy, policies::promote_float, policies::promote_double, policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; - value_type result = detail::find_t_degrees_of_freedom( - static_cast(delta), - static_cast(x), - static_cast(p), - static_cast(1-p), + eval_type result = detail::find_t_degrees_of_freedom( + static_cast(delta), + static_cast(x), + static_cast(p), + static_cast(1-p), forwarding_policy()); return policies::checked_narrowing_cast( result, @@ -893,18 +889,18 @@ namespace boost static RealType find_degrees_of_freedom(const complemented3_type& c) { const char* function = "non_central_t<%1%>::find_degrees_of_freedom"; - typedef typename policies::evaluation::type value_type; + typedef typename policies::evaluation::type eval_type; typedef typename policies::normalise< Policy, policies::promote_float, policies::promote_double, policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; - value_type result = detail::find_t_degrees_of_freedom( - static_cast(c.dist), - static_cast(c.param1), - static_cast(1-c.param2), - static_cast(c.param2), + eval_type result = detail::find_t_degrees_of_freedom( + static_cast(c.dist), + static_cast(c.param1), + static_cast(1-c.param2), + static_cast(c.param2), forwarding_policy()); return policies::checked_narrowing_cast( result, @@ -913,18 +909,18 @@ namespace boost static RealType find_non_centrality(RealType v, RealType x, RealType p) { const char* function = "non_central_t<%1%>::find_t_non_centrality"; - typedef typename policies::evaluation::type value_type; + typedef typename policies::evaluation::type eval_type; typedef typename policies::normalise< Policy, policies::promote_float, policies::promote_double, policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; - value_type result = detail::find_t_non_centrality( - static_cast(v), - static_cast(x), - static_cast(p), - static_cast(1-p), + eval_type result = detail::find_t_non_centrality( + static_cast(v), + static_cast(x), + static_cast(p), + static_cast(1-p), forwarding_policy()); return policies::checked_narrowing_cast( result, @@ -934,24 +930,23 @@ namespace boost static RealType find_non_centrality(const complemented3_type& c) { const char* function = "non_central_t<%1%>::find_t_non_centrality"; - typedef typename policies::evaluation::type value_type; + typedef typename policies::evaluation::type eval_type; typedef typename policies::normalise< Policy, policies::promote_float, policies::promote_double, policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; - value_type result = detail::find_t_non_centrality( - static_cast(c.dist), - static_cast(c.param1), - static_cast(1-c.param2), - static_cast(c.param2), + eval_type result = detail::find_t_non_centrality( + static_cast(c.dist), + static_cast(c.param1), + static_cast(1-c.param2), + static_cast(c.param2), forwarding_policy()); return policies::checked_narrowing_cast( result, function); } -#endif private: // Data member, initialized by constructor. RealType v; // degrees of freedom diff --git a/test/test_nc_t.hpp b/test/test_nc_t.hpp index 8e0edabdc1..58a48511dc 100644 --- a/test/test_nc_t.hpp +++ b/test/test_nc_t.hpp @@ -494,7 +494,6 @@ void quantile_sanity_check(T& data, const char* type_name, const char* test) } catch(const boost::math::evaluation_error&) {} #endif -#if 0 // // Sanity check degrees-of-freedom finder, don't bother at float // precision though as there's not enough data in the probability @@ -502,30 +501,46 @@ void quantile_sanity_check(T& data, const char* type_name, const char* test) // non-centrality parameter: // try{ + Real non_centrality_precision_multiplier = 1; + Real df_precision_multiplier = 1; + if (data[i][0] > 10000) + { + non_centrality_precision_multiplier = 500; + df_precision_multiplier = 500; + } + if (data[i][0] > 1000000000) + { + df_precision_multiplier *= 50; // very little precision left at this point + } if((data[i][3] < 0.99) && (data[i][3] != 0)) { - BOOST_CHECK_CLOSE_EX( - boost::math::non_central_t_distribution::find_degrees_of_freedom(data[i][1], data[i][2], data[i][3]), - data[i][0], precision, i); + if (data[i][0] < 1.0e12) // no precision above this + { + BOOST_CHECK_CLOSE_EX( + boost::math::non_central_t_distribution::find_degrees_of_freedom(data[i][1], data[i][2], data[i][3]), + data[i][0], precision * df_precision_multiplier, i); + } BOOST_CHECK_CLOSE_EX( boost::math::non_central_t_distribution::find_non_centrality(data[i][0], data[i][2], data[i][3]), - data[i][1], precision, i); + data[i][1], precision * non_centrality_precision_multiplier, i); } if((data[i][4] < 0.99) && (data[i][4] != 0)) { - BOOST_CHECK_CLOSE_EX( - boost::math::non_central_t_distribution::find_degrees_of_freedom(boost::math::complement(data[i][1], data[i][2], data[i][4])), - data[i][0], precision, i); + if (data[i][0] < 1.0e12) // no precision above this + { + BOOST_CHECK_CLOSE_EX( + boost::math::non_central_t_distribution::find_degrees_of_freedom(boost::math::complement(data[i][1], data[i][2], data[i][4])), + data[i][0], precision * df_precision_multiplier, i); + } BOOST_CHECK_CLOSE_EX( boost::math::non_central_t_distribution::find_non_centrality(boost::math::complement(data[i][0], data[i][2], data[i][4])), - data[i][1], precision, i); + data[i][1], precision * non_centrality_precision_multiplier, i); } } catch(const std::exception& e) { BOOST_ERROR(e.what()); } -#endif } } #endif