diff --git a/include/boost/math/distributions/non_central_f.hpp b/include/boost/math/distributions/non_central_f.hpp index 783c321980..9bae5c41ff 100644 --- a/include/boost/math/distributions/non_central_f.hpp +++ b/include/boost/math/distributions/non_central_f.hpp @@ -14,6 +14,7 @@ #include #include #include +#include #include #include #include @@ -44,7 +45,8 @@ namespace boost }; template - BOOST_MATH_GPU_ENABLED RealType find_non_centrality_f(const RealType x, const RealType v1, const RealType v2, const RealType p, const RealType q, const RealType p_q_precision, const Policy& pol) { + BOOST_MATH_GPU_ENABLED RealType find_non_centrality_f(const RealType x, const RealType v1, const RealType v2, const RealType p, const RealType q, const RealType p_q_precision, const Policy& pol) + { constexpr auto function = "non_central_f<%1%>::find_non_centrality"; if ( p == 0 || q == 0) { @@ -80,6 +82,103 @@ namespace boost } return result; } + + template + struct f_degrees_of_freedom_finder + { + f_degrees_of_freedom_finder( + RealType x_, RealType v_, RealType nc_, bool find_v1_, RealType p_, bool c) + : x(x_), v(v_), nc(nc_), find_v1(find_v1_), p(p_), comp(c) {} + + RealType operator()(const RealType& input_v) + { + RealType v1 = find_v1 ? input_v : v; + RealType v2 = find_v1 ? v : input_v; + non_central_f_distribution d(v1, v2, nc); + return comp ? + p - cdf(complement(d, x)) + : cdf(d, x) - p; + } + private: + RealType x; + RealType v; + RealType nc; + bool find_v1; + RealType p; + bool comp; + }; + + template + RealType large_v2_approximation(RealType x, RealType v1, RealType p, RealType q, RealType nc) + { // For v2 -> inf approximate f_degreese_of_freedome_finder with chi squared distribution + // with degrees of freedom v1 at the cdf at x * v1 + bool comp = p < q ? false : true; + RealType pval = p < q ? p : q; + non_central_chi_squared_distribution d(v1, nc); + return comp ? pval - cdf(complement(d, x*v1)) : cdf(d, x*v1) - pval; + } + + template + inline RealType find_degrees_of_freedom_f( + const RealType x, const RealType v, const RealType nc, const bool find_v1, const RealType p, const RealType q, const Policy& pol) + { + BOOST_MATH_STD_USING + using std::fabs; + const char* function = "non_central_f<%1%>::find_degrees_of_freedom_f"; + if((p == 0) || (q == 0)) + { + // + // Can't find a thing if one of p and q is zero: + // + return policies::raise_domain_error(function, "Can't find degrees of freedom when the probability is 0 or 1, only possible answer is %1%", + RealType(std::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE + } + if (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 + } + f_degrees_of_freedom_finder f(x, v, nc, find_v1, p < q ? p : q, p < q ? false : true); + + // There are times when f has two roots - thus, two degrees of freedom can + // be found. We find this case by checking if the sign of f for large and + // small values of v have the same sign. If the sign is the same, then there + // are an even number of roots. If the signs differ, there is only one root + // and we can safely find the root. + RealType vLarge = sqrt(boost::math::tools::max_value()); + RealType vSmall = 1 / vLarge; + + // As v2 -> infinity, noncentral f converges to chi-squared distribution. + // Rather than evaluating f(vLarge), we can use the more stable chi-squared approximation + RealType large_difference; + if (find_v1) + { + large_difference = large_v2_approximation(x, v, p, q, nc); + } + else + large_difference = f(vLarge); + + if ((large_difference < 0) == (f(vSmall) < 0)){ + return policies::raise_evaluation_error(function, "Can't find degrees of freedom because two degrees of freedom can be found using the given parameters", + RealType(std::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE + } + + tools::eps_tolerance tol(policies::digits()); + std::uintmax_t max_iter = policies::get_max_root_iterations(); + // + // Pick an initial guess: + // + RealType guess = 1; + std::pair ir = tools::bracket_and_solve_root( + f, guess, RealType(2), f(guess) < 0 ? true : false, tol, max_iter, pol); + RealType result = ir.first + (ir.second - ir.first) / 2; + if(max_iter >= policies::get_max_root_iterations()) + { + return policies::raise_evaluation_error(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE + " or there is no answer to problem. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE + } + return result; + } } // namespace detail template > @@ -163,6 +262,96 @@ namespace boost result, function); } + BOOST_MATH_GPU_ENABLED static RealType find_v1(const RealType x, const RealType v2, const RealType nc, const RealType p) + { + constexpr auto function = "non_central_f_distribution<%1%>::find_v1"; + 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; + eval_type result = detail::find_degrees_of_freedom_f( + static_cast(x), + static_cast(v2), + static_cast(nc), + true, + static_cast(p), + static_cast(1-p), + forwarding_policy()); + return policies::checked_narrowing_cast( + result, + function); + } + template + BOOST_MATH_GPU_ENABLED static RealType find_v1(const complemented4_type& c) + { + constexpr auto function = "non_central_f_distribution<%1%>::find_non_centrality"; + 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; + eval_type result = detail::find_degrees_of_freedom_f( + static_cast(c.dist), + static_cast(c.param1), + static_cast(c.param2), + true, + static_cast(1-c.param3), + static_cast(c.param3), + forwarding_policy()); + return policies::checked_narrowing_cast( + result, + function); + } + BOOST_MATH_GPU_ENABLED static RealType find_v2(const RealType x, const RealType v2, const RealType nc, const RealType p) + { + constexpr auto function = "non_central_f_distribution<%1%>::find_v1"; + 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; + eval_type result = detail::find_degrees_of_freedom_f( + static_cast(x), + static_cast(v2), + static_cast(nc), + false, + static_cast(p), + static_cast(1-p), + forwarding_policy()); + return policies::checked_narrowing_cast( + result, + function); + } + template + BOOST_MATH_GPU_ENABLED static RealType find_v2(const complemented4_type& c) + { + constexpr auto function = "non_central_f_distribution<%1%>::find_non_centrality"; + 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; + eval_type result = detail::find_degrees_of_freedom_f( + static_cast(c.dist), + static_cast(c.param1), + static_cast(c.param2), + false, + static_cast(1-c.param3), + static_cast(c.param3), + forwarding_policy()); + return policies::checked_narrowing_cast( + result, + function); + } private: // Data member, initialized by constructor. RealType v1; // alpha. diff --git a/include/boost/math/distributions/students_t.hpp b/include/boost/math/distributions/students_t.hpp index bc1c7139eb..559b155984 100644 --- a/include/boost/math/distributions/students_t.hpp +++ b/include/boost/math/distributions/students_t.hpp @@ -58,7 +58,7 @@ class students_t_distribution RealType alpha, RealType beta, RealType sd, - RealType hint = 100); + RealType hint = 0); BOOST_MATH_GPU_ENABLED static RealType find_degrees_of_freedom( RealType t, @@ -315,8 +315,31 @@ struct sample_size_func return 1; } students_t_distribution t(df); - RealType qa = quantile(complement(t, alpha)); - RealType qb = quantile(complement(t, beta)); + RealType qa, qb; +#ifndef BOOST_MATH_NO_EXCEPTIONS + try { +#endif + qa = quantile(complement(t, alpha)); +#ifndef BOOST_MATH_NO_EXCEPTIONS + } + catch (const std::overflow_error&) + { + // Any arbitrary large value will do, must be positive since we calculate qa * qa below: + return tools::max_value(); + } +#endif +#ifndef BOOST_MATH_NO_EXCEPTIONS + try { +#endif + qb = quantile(complement(t, beta)); +#ifndef BOOST_MATH_NO_EXCEPTIONS + } + catch (const std::overflow_error&) + { + // Any arbitrary large value will do, will be negative if beta < 0: + return beta < 0 ? -tools::max_value() : tools::max_value(); + } +#endif qa += qb; qa *= qa; qa *= ratio; @@ -413,6 +436,30 @@ BOOST_MATH_GPU_ENABLED bool approximate_df_with_edgeworth_expansion(RealType x, return false; } +template +BOOST_MATH_GPU_ENABLED RealType calculate_hill_df_guess( + RealType difference_from_mean, + RealType alpha, + RealType beta, + RealType sd) + + // Hill-corrected normal approximation used as an initial guess + // for find_degrees_of_freedom(difference_from_mean, alpha, beta, sd). + // t_p,nu approx z_p * (1 + (z_p^2 + 1) / (4*nu)) + // nu + 1 = (sd/diff)^2 * (t_alpha + t_beta)^2 + +{ + BOOST_MATH_STD_USING + normal_distribution n(0, 1); + RealType za = quantile(complement(n, alpha)); + RealType zb = quantile(complement(n, beta)); + + RealType v_norm = pow((za + zb) * sd / difference_from_mean, 2) - 1; + RealType hint = v_norm + (za * za - za * zb + zb * zb + 1) / 2; + + return (hint <= 0) ? RealType(1) : hint; +} + } // namespace detail template @@ -434,7 +481,7 @@ BOOST_MATH_GPU_ENABLED RealType students_t_distribution::find_ return error_result; if(hint <= 0) - hint = 1; + hint = detail::calculate_hill_df_guess(difference_from_mean, alpha, beta, sd); detail::sample_size_func f(alpha, beta, sd, difference_from_mean); return detail::solve_for_degrees_of_freedom(f, hint, false, function, Policy()); diff --git a/test/test_nc_f.cpp b/test/test_nc_f.cpp index 6b7c2347ab..7ad45c1ba3 100644 --- a/test/test_nc_f.cpp +++ b/test/test_nc_f.cpp @@ -138,13 +138,13 @@ void test_spot( BOOST_CHECK_CLOSE( cdf(complement(dist, x)), Q, tol); BOOST_CHECK_CLOSE( - quantile(dist, P), x, tol * 10); + quantile(dist, P), x, tol * 10); BOOST_CHECK_CLOSE( - quantile(complement(dist, Q)), x, tol * 10); + quantile(complement(dist, Q)), x, tol * 10); BOOST_CHECK_CLOSE( - dist.find_non_centrality(x, a, b, P), ncp, tol * 10); + dist.find_non_centrality(x, a, b, P), ncp, tol * 10); BOOST_CHECK_CLOSE( - dist.find_non_centrality(boost::math::complement(x, a, b, Q)), ncp, tol * 10); + dist.find_non_centrality(boost::math::complement(x, a, b, Q)), ncp, tol * 10); } if(boost::math::tools::digits() > 50) { @@ -369,6 +369,48 @@ void test_spots(RealType, const char* name = nullptr) } } } + + // Quick spot check for finding degrees of freedom. When checking for two degrees + // of freedom for real_concept types, the cdf at large/small v2 can be greater than 1 + // or less than 0. + if (!std::is_same::value){ + RealType v1 = 10; + RealType v2 = 5; + nc = 1; + x = 6; + boost::math::non_central_f_distribution ref(v1, v2, nc); + RealType P = cdf(ref, x); + BOOST_CHECK_CLOSE(ref.find_v2(x, v1, nc, P), v2, tolerance); + BOOST_CHECK_CLOSE(ref.find_v2(boost::math::complement(x, v1, nc, 1-P)), v2, tolerance); + BOOST_CHECK_CLOSE(ref.find_v1(x, v2, nc, P), v1, tolerance); + BOOST_CHECK_CLOSE(ref.find_v1(boost::math::complement(x, v2, nc, 1-P)), v1, tolerance); + } + + // Check case where two degrees of freedom solve the inversion problem + BOOST_MATH_CHECK_THROW(dist.find_v1(RealType(1.5), RealType(2.0), RealType(1.0), RealType(0.49845842011686358665786775091245664L)), boost::math::evaluation_error); + BOOST_MATH_CHECK_THROW(dist.find_v1(RealType(3.51), RealType(5), RealType(0), RealType(0.85802971653663762108266155337333L)), boost::math::evaluation_error); + + // Check find_v1/v2 edge cases + // Case when P=1 or P=0 + nc = 2; + BOOST_MATH_CHECK_THROW(dist.find_v1(x, b, nc, 1), std::domain_error); + BOOST_MATH_CHECK_THROW(dist.find_v1(x, b, nc, 0), std::domain_error); + // Case when Q=1 or Q=0 + BOOST_MATH_CHECK_THROW(dist.find_v1(boost::math::complement(x, b, nc, 1)), std::domain_error); + BOOST_MATH_CHECK_THROW(dist.find_v1(boost::math::complement(x, b, nc, 0)), std::domain_error); + // Check very small values of x an evaluation error is thrown + x = boost::math::tools::epsilon() / 10; + BOOST_MATH_CHECK_THROW(dist.find_v1(boost::math::complement(x, b, nc, 0.5)), boost::math::evaluation_error); + BOOST_MATH_CHECK_THROW(dist.find_v1(x, b, nc, 0.5), boost::math::evaluation_error); + + BOOST_MATH_CHECK_THROW(dist.find_v2(x, b, nc, 1), std::domain_error); + BOOST_MATH_CHECK_THROW(dist.find_v2(x, b, nc, 0), std::domain_error); + // Case when Q=1 or Q=0 + BOOST_MATH_CHECK_THROW(dist.find_v2(boost::math::complement(x, b, nc, 1)), std::domain_error); + BOOST_MATH_CHECK_THROW(dist.find_v2(boost::math::complement(x, b, nc, 0)), std::domain_error); + // Check very small values of x an evaluation error is thrown + BOOST_MATH_CHECK_THROW(dist.find_v2(boost::math::complement(x, b, nc, 0.5)), boost::math::evaluation_error); + BOOST_MATH_CHECK_THROW(dist.find_v2(x, b, nc, 0.5), boost::math::evaluation_error); } // template void test_spots(RealType) BOOST_AUTO_TEST_CASE( test_main )