diff --git a/include/boost/math/distributions/non_central_f.hpp b/include/boost/math/distributions/non_central_f.hpp index 2acaab255..02a75a28d 100644 --- a/include/boost/math/distributions/non_central_f.hpp +++ b/include/boost/math/distributions/non_central_f.hpp @@ -18,6 +18,8 @@ #include #include #include // complements +#include +#include namespace boost { @@ -53,9 +55,18 @@ namespace boost RealType(boost::math::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE } + // Check if nc = 0 (which is just the F-distribution) + // I first tried comparing to boost::math::tools::epsilon(), + // but this was never satisfied for floats or doubles. Only + // long doubles would correctly return 0. + fisher_f_distribution dist(v1, v2); + if (boost::math::relative_difference(pdf(dist, x), p) < 1e-7){ + return 0; + } + non_centrality_finder_f f(x, v1, v2, p < q ? p : q, p < q ? false : true); RealType guess = RealType(10); // Starting guess. - RealType factor = 1; // How big steps to take when searching. + RealType factor = RealType(2); // How big steps to take when searching. boost::math::uintmax_t max_iter = policies::get_max_root_iterations(); tools::eps_tolerance tol(policies::digits()); diff --git a/test/test_nc_f.cpp b/test/test_nc_f.cpp index fdd03ab28..3ccb9e3b4 100644 --- a/test/test_nc_f.cpp +++ b/test/test_nc_f.cpp @@ -327,6 +327,25 @@ void test_spots(RealType) { BOOST_CHECK_CLOSE(cdf(boost::math::non_central_f_distribution(static_cast(1e-100L), 3.f, 1.5f), static_cast(1e100L)), static_cast(0.6118152873453990639132215575213809716459L), tolerance); } + + // Check find_non_centrality_f edge case handling + // Case when nc=0 + RealType a = 5; + RealType b = 2; + RealType nc = 0; + RealType x_vals[] = { 0.25, 1.25, 10, 100}; + boost::math::non_central_f_distribution dist_no_centrality(a, b, nc); + for (RealType x : x_vals) + { + RealType P = pdf(dist_no_centrality, x); + BOOST_CHECK_CLOSE(dist.find_non_centrality(x, a, b, P), static_cast(0), tolerance); + } + // Case when P=1 or P=0 + BOOST_MATH_CHECK_THROW(dist.find_non_centrality(x, a, b, 1), std::domain_error); + BOOST_MATH_CHECK_THROW(dist.find_non_centrality(x, a, b, 0), std::domain_error); + // Case when Q=1 or Q=0 + BOOST_MATH_CHECK_THROW(dist.find_non_centrality(boost::math::complement(x, a, b, 1)), std::domain_error); + BOOST_MATH_CHECK_THROW(dist.find_non_centrality(boost::math::complement(x, a, b, 0)), std::domain_error); } // template void test_spots(RealType) BOOST_AUTO_TEST_CASE( test_main )