From 44cbb533cc0a5cc49ccb236bdf630f16be29ccd2 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Fri, 20 Feb 2026 11:02:57 -0800 Subject: [PATCH 1/5] 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 2/5] 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 3/5] 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 4/5] 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 5/5] 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);