From 44cbb533cc0a5cc49ccb236bdf630f16be29ccd2 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Fri, 20 Feb 2026 11:02:57 -0800 Subject: [PATCH 01/11] Initial implementation/test of find_v1 --- .../math/distributions/non_central_f.hpp | 106 +++++++++++++++++- test/test_nc_f.cpp | 12 +- 2 files changed, 113 insertions(+), 5 deletions(-) diff --git a/include/boost/math/distributions/non_central_f.hpp b/include/boost/math/distributions/non_central_f.hpp index 783c321980..8efeb6216c 100644 --- a/include/boost/math/distributions/non_central_f.hpp +++ b/include/boost/math/distributions/non_central_f.hpp @@ -44,7 +44,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 +81,66 @@ namespace boost } return result; } + + // Find first degrees of freedom + template + struct f_v1_degrees_of_freedom_finder + { + f_v1_degrees_of_freedom_finder( + RealType x_, RealType v2_, RealType nc_, RealType p_, bool c) + : x(x_), v2(v2_), nc(nc_), p(p_), comp(c) {} + + RealType operator()(const RealType& v1) + { + non_central_f_distribution d(v1, v2, nc); + return comp ? + p - cdf(complement(d, x)) + : cdf(d, x) - p; + } + private: + RealType x; + RealType v2; + RealType nc; + RealType p; + bool comp; + }; + + template + inline RealType find_v1_degrees_of_freedom_f( + const RealType x, const RealType v2, const RealType nc, const RealType p, const RealType q, const Policy& pol) + { + using std::fabs; + const char* function = "non_central_f<%1%>::find_v1_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_evaluation_error(function, "Can't find v1 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 (fabs(x) < tools::epsilon()) + { + return policies::raise_evaluation_error(function, "Can't find v1 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_v1_degrees_of_freedom_finder f(x, v2, nc, p < q ? p : q, p < q ? false : true); + tools::eps_tolerance tol(policies::digits()); + std::uintmax_t max_iter = policies::get_max_root_iterations(); + // + // Pick an initial guess: + // + RealType guess = 200; + std::pair ir = tools::bracket_and_solve_root( + 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()) + { + 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 +224,49 @@ 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_v1_degrees_of_freedom_f( + static_cast(x), + static_cast(v2), + static_cast(nc), + 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_v1_degrees_of_freedom_f( + static_cast(c.dist), + static_cast(c.param1), + static_cast(c.param2), + 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/test/test_nc_f.cpp b/test/test_nc_f.cpp index 6b7c2347ab..96664e4d9a 100644 --- a/test/test_nc_f.cpp +++ b/test/test_nc_f.cpp @@ -138,13 +138,17 @@ 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); + BOOST_CHECK_CLOSE( + dist.find_v1(x, b, ncp, P), a, tol * 10); + BOOST_CHECK_CLOSE( + dist.find_v1(boost::math::complement(x, b, ncp, Q)), a, tol * 10); } if(boost::math::tools::digits() > 50) { From 3eb7561296c5a15cc2c891a68e5c7f7a624d31e4 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Fri, 20 Feb 2026 16:53:26 -0800 Subject: [PATCH 02/11] Added v2 finder --- .../math/distributions/non_central_f.hpp | 83 +++++++++++++++---- test/test_nc_f.cpp | 4 + 2 files changed, 70 insertions(+), 17 deletions(-) diff --git a/include/boost/math/distributions/non_central_f.hpp b/include/boost/math/distributions/non_central_f.hpp index 8efeb6216c..3112aeedb4 100644 --- a/include/boost/math/distributions/non_central_f.hpp +++ b/include/boost/math/distributions/non_central_f.hpp @@ -82,16 +82,17 @@ namespace boost return result; } - // Find first degrees of freedom template - struct f_v1_degrees_of_freedom_finder + struct f_degrees_of_freedom_finder { - f_v1_degrees_of_freedom_finder( - RealType x_, RealType v2_, RealType nc_, RealType p_, bool c) - : x(x_), v2(v2_), nc(nc_), p(p_), comp(c) {} + 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& v1) + 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)) @@ -99,40 +100,41 @@ namespace boost } private: RealType x; - RealType v2; + RealType v; RealType nc; + bool find_v1; RealType p; bool comp; }; template - inline RealType find_v1_degrees_of_freedom_f( - const RealType x, const RealType v2, const RealType nc, const RealType p, const RealType q, const Policy& pol) + 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) { using std::fabs; - const char* function = "non_central_f<%1%>::find_v1_degrees_of_freedom_f"; + 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_evaluation_error(function, "Can't find v1 degrees of freedom when the probability is 0 or 1, only possible answer is %1%", + return policies::raise_evaluation_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 (fabs(x) < tools::epsilon()) { - return policies::raise_evaluation_error(function, "Can't find v1 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!!", + 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_v1_degrees_of_freedom_finder f(x, v2, nc, p < q ? p : q, p < q ? false : true); + f_degrees_of_freedom_finder f(x, v, nc, find_v1, p < q ? p : q, p < q ? false : true); tools::eps_tolerance tol(policies::digits()); std::uintmax_t max_iter = policies::get_max_root_iterations(); // // Pick an initial guess: // - RealType guess = 200; + RealType guess = 1; std::pair ir = tools::bracket_and_solve_root( - f, guess, RealType(2), x < 0 ? false : true, tol, max_iter, pol); + 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()) { @@ -234,10 +236,11 @@ namespace boost policies::promote_double, policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; - eval_type result = detail::find_v1_degrees_of_freedom_f( + 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()); @@ -256,10 +259,56 @@ namespace boost policies::promote_double, policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; - eval_type result = detail::find_v1_degrees_of_freedom_f( + 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()); diff --git a/test/test_nc_f.cpp b/test/test_nc_f.cpp index 96664e4d9a..5a61d31947 100644 --- a/test/test_nc_f.cpp +++ b/test/test_nc_f.cpp @@ -149,6 +149,10 @@ void test_spot( dist.find_v1(x, b, ncp, P), a, tol * 10); BOOST_CHECK_CLOSE( dist.find_v1(boost::math::complement(x, b, ncp, Q)), a, tol * 10); + BOOST_CHECK_CLOSE( + dist.find_v2(x, a, ncp, P), b, tol * 10); + BOOST_CHECK_CLOSE( + dist.find_v2(boost::math::complement(x, a, ncp, Q)), b, tol * 10); } if(boost::math::tools::digits() > 50) { From 48016543b9a0a9545266d7242c3f35c586bf1952 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Sat, 21 Feb 2026 10:49:08 -0800 Subject: [PATCH 03/11] Included edge cases for code coverage [ci skip] --- include/boost/math/distributions/non_central_f.hpp | 4 ++-- test/test_nc_f.cpp | 12 ++++++++++++ 2 files changed, 14 insertions(+), 2 deletions(-) diff --git a/include/boost/math/distributions/non_central_f.hpp b/include/boost/math/distributions/non_central_f.hpp index 3112aeedb4..5831f6622c 100644 --- a/include/boost/math/distributions/non_central_f.hpp +++ b/include/boost/math/distributions/non_central_f.hpp @@ -118,10 +118,10 @@ namespace boost // // 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%", + 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 (fabs(x) < tools::epsilon()) + 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 diff --git a/test/test_nc_f.cpp b/test/test_nc_f.cpp index 5a61d31947..d96675320d 100644 --- a/test/test_nc_f.cpp +++ b/test/test_nc_f.cpp @@ -377,6 +377,18 @@ void test_spots(RealType, const char* name = nullptr) } } } + // 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); } // template void test_spots(RealType) BOOST_AUTO_TEST_CASE( test_main ) From 04eace43b14ae206f96fa8e5d2e94c16f21b4f3a Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Sun, 24 May 2026 11:11:05 -0700 Subject: [PATCH 04/11] Added checks for multiple degrees of freedom --- .../math/distributions/non_central_f.hpp | 15 +++++++++ test/test_nc_f.cpp | 33 ++++++++++++++----- 2 files changed, 39 insertions(+), 9 deletions(-) diff --git a/include/boost/math/distributions/non_central_f.hpp b/include/boost/math/distributions/non_central_f.hpp index 5831f6622c..ccfa18f65c 100644 --- a/include/boost/math/distributions/non_central_f.hpp +++ b/include/boost/math/distributions/non_central_f.hpp @@ -111,6 +111,7 @@ namespace boost 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)) @@ -127,6 +128,20 @@ namespace boost 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; + + if ((f(vLarge) < 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(); // diff --git a/test/test_nc_f.cpp b/test/test_nc_f.cpp index d96675320d..78d4eea701 100644 --- a/test/test_nc_f.cpp +++ b/test/test_nc_f.cpp @@ -145,14 +145,6 @@ void test_spot( 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); - BOOST_CHECK_CLOSE( - dist.find_v1(x, b, ncp, P), a, tol * 10); - BOOST_CHECK_CLOSE( - dist.find_v1(boost::math::complement(x, b, ncp, Q)), a, tol * 10); - BOOST_CHECK_CLOSE( - dist.find_v2(x, a, ncp, P), b, tol * 10); - BOOST_CHECK_CLOSE( - dist.find_v2(boost::math::complement(x, a, ncp, Q)), b, tol * 10); } if(boost::math::tools::digits() > 50) { @@ -377,8 +369,22 @@ void test_spots(RealType, const char* name = nullptr) } } } + + // Quick spot check for finding degrees of freedom + RealType v1 = 5; + RealType v2 = 2; + nc = 1; + x = 6; + boost::math::non_central_f_distribution ref(v1, v2, nc); + RealType P = cdf(ref, x); + BOOST_CHECK_CLOSE(ref.find_v1(x, v2, nc, 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 + // 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); @@ -389,6 +395,15 @@ void test_spots(RealType, const char* name = nullptr) 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 ) From d43d08c92d8cce7916d69cb9a3e3b2d9065ec52d Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Mon, 25 May 2026 12:18:04 -0700 Subject: [PATCH 05/11] Improved test coverage --- test/test_nc_f.cpp | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/test/test_nc_f.cpp b/test/test_nc_f.cpp index 78d4eea701..c836ddc5fb 100644 --- a/test/test_nc_f.cpp +++ b/test/test_nc_f.cpp @@ -370,14 +370,21 @@ void test_spots(RealType, const char* name = nullptr) } } - // Quick spot check for finding degrees of freedom - RealType v1 = 5; - RealType v2 = 2; - nc = 1; - x = 6; - boost::math::non_central_f_distribution ref(v1, v2, nc); - RealType P = cdf(ref, x); - BOOST_CHECK_CLOSE(ref.find_v1(x, v2, nc, P), v1, tolerance); + // 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); From bdce34ddf77b22dd4cefe632b7fc28e13f756757 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Tue, 2 Jun 2026 10:05:01 -0700 Subject: [PATCH 06/11] Added chi squared approximation --- .../math/distributions/non_central_f.hpp | 23 ++++++++++++++++++- 1 file changed, 22 insertions(+), 1 deletion(-) diff --git a/include/boost/math/distributions/non_central_f.hpp b/include/boost/math/distributions/non_central_f.hpp index ccfa18f65c..fcde87b2d2 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 @@ -107,6 +108,16 @@ namespace boost bool comp; }; + template + RealType large_v2_approximation(RealType x, RealType v1, RealType p, RealType q) + { // 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; + chi_squared_distribution d(v1); + 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) @@ -137,7 +148,17 @@ namespace boost RealType vLarge = sqrt(boost::math::tools::max_value()); RealType vSmall = 1 / vLarge; - if ((f(vLarge) < 0) == (f(vSmall) < 0)){ + // 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); + } + 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 } From a0bd79c0e2c9c29ce3ec29434c06bc36b4bc6edd Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Tue, 2 Jun 2026 13:18:04 -0700 Subject: [PATCH 07/11] Changed to non-central chi-squared dist --- include/boost/math/distributions/non_central_f.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/include/boost/math/distributions/non_central_f.hpp b/include/boost/math/distributions/non_central_f.hpp index fcde87b2d2..9bae5c41ff 100644 --- a/include/boost/math/distributions/non_central_f.hpp +++ b/include/boost/math/distributions/non_central_f.hpp @@ -14,7 +14,7 @@ #include #include #include -#include +#include #include #include #include @@ -109,12 +109,12 @@ namespace boost }; template - RealType large_v2_approximation(RealType x, RealType v1, RealType p, RealType q) + 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; - chi_squared_distribution d(v1); + non_central_chi_squared_distribution d(v1, nc); return comp ? pval - cdf(complement(d, x*v1)) : cdf(d, x*v1) - pval; } @@ -153,7 +153,7 @@ namespace boost RealType large_difference; if (find_v1) { - large_difference = large_v2_approximation(x, v, p, q); + large_difference = large_v2_approximation(x, v, p, q, nc); } else large_difference = f(vLarge); From 38ebc9a6cfa84e86fcaec62d0e261cec70f2041e Mon Sep 17 00:00:00 2001 From: dschmitz89 Date: Mon, 8 Jun 2026 21:16:39 +0200 Subject: [PATCH 08/11] ENH: add warm start for root finder --- .../boost/math/distributions/students_t.hpp | 29 +++++++++++++++++-- 1 file changed, 27 insertions(+), 2 deletions(-) diff --git a/include/boost/math/distributions/students_t.hpp b/include/boost/math/distributions/students_t.hpp index bc1c7139eb..50f591b690 100644 --- a/include/boost/math/distributions/students_t.hpp +++ b/include/boost/math/distributions/students_t.hpp @@ -25,6 +25,7 @@ #include #include #include +#include #ifdef _MSC_VER # pragma warning(push) @@ -58,7 +59,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, @@ -413,6 +414,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 +459,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()); From ce299579d420ce37f563d561bdb8af7715deb1a7 Mon Sep 17 00:00:00 2001 From: dschmitz89 Date: Mon, 8 Jun 2026 21:24:00 +0200 Subject: [PATCH 09/11] ENH: add initial guess for power equation inversion --- include/boost/math/distributions/students_t.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/include/boost/math/distributions/students_t.hpp b/include/boost/math/distributions/students_t.hpp index 50f591b690..94fe686c2d 100644 --- a/include/boost/math/distributions/students_t.hpp +++ b/include/boost/math/distributions/students_t.hpp @@ -25,7 +25,6 @@ #include #include #include -#include #ifdef _MSC_VER # pragma warning(push) From 4426160fb2afab7a08aae6aca9d62de2d9c3c927 Mon Sep 17 00:00:00 2001 From: Jacob Hass <72838561+JacobHass8@users.noreply.github.com> Date: Mon, 8 Jun 2026 17:21:48 -0700 Subject: [PATCH 10/11] Fix type for testing small x --- test/test_nc_f.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_nc_f.cpp b/test/test_nc_f.cpp index c836ddc5fb..7ad45c1ba3 100644 --- a/test/test_nc_f.cpp +++ b/test/test_nc_f.cpp @@ -399,7 +399,7 @@ void test_spots(RealType, const char* name = nullptr) 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; + 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); From 22da3792e592c3a67c5cba998e88df89807e6fcb Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Tue, 9 Jun 2026 17:51:20 +0100 Subject: [PATCH 11/11] students_t: Add exception handling guards. The degree of freedom finders can call the quantile with arguments which result in an infinite result, as far as our root finders are concerned, any large value as a result will do, but we need to not allow an exception to escape or the root finder terminates. --- .../boost/math/distributions/students_t.hpp | 27 +++++++++++++++++-- 1 file changed, 25 insertions(+), 2 deletions(-) diff --git a/include/boost/math/distributions/students_t.hpp b/include/boost/math/distributions/students_t.hpp index 94fe686c2d..559b155984 100644 --- a/include/boost/math/distributions/students_t.hpp +++ b/include/boost/math/distributions/students_t.hpp @@ -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;