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