diff --git a/include/boost/math/distributions/fisher_f.hpp b/include/boost/math/distributions/fisher_f.hpp index 56b288d88e..464d1cda46 100644 --- a/include/boost/math/distributions/fisher_f.hpp +++ b/include/boost/math/distributions/fisher_f.hpp @@ -17,7 +17,84 @@ #include // error checks #include -namespace boost{ namespace math{ +namespace boost{ namespace math{ + namespace detail{ + template + struct fisher_degrees_of_freedom_finder + { + fisher_degrees_of_freedom_finder( + RealType x_, RealType v_, bool find_v1_, RealType p_, bool c) + : x(x_), v(v_), 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; + fisher_f_distribution d(v1, v2); + return comp ? + p - cdf(complement(d, x)) + : cdf(d, x) - p; + } + private: + RealType x; + RealType v; + bool find_v1; + RealType p; + bool comp; + }; + + template + inline RealType find_degrees_of_freedom_fisher_f( + const RealType x, const RealType v, const bool find_v1, const RealType p, const RealType q, const Policy& pol) + { + BOOST_MATH_STD_USING + using std::fabs; + const char* function = "fisher_f_distribution<%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 + } + fisher_degrees_of_freedom_finder f(x, v, 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(); + // + // 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; + } + } template > class fisher_f_distribution @@ -44,7 +121,92 @@ class fisher_f_distribution { return m_df2; } - + BOOST_MATH_GPU_ENABLED static RealType find_v1(const RealType x, const RealType v2, const RealType p) + { + constexpr auto function = "fisher_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_fisher_f( + static_cast(x), + static_cast(v2), + 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 complemented3_type& c) + { + constexpr auto function = "fisher_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_fisher_f( + static_cast(c.dist), + static_cast(c.param1), + true, + static_cast(1-c.param2), + static_cast(c.param2), + 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 p) + { + constexpr auto function = "fisher_f_distribution<%1%>::find_v2"; + 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_fisher_f( + static_cast(x), + static_cast(v2), + 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 complemented3_type& c) + { + constexpr auto function = "fisher_f_distribution<%1%>::find_v2"; + 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_fisher_f( + static_cast(c.dist), + static_cast(c.param1), + false, + static_cast(1-c.param2), + static_cast(c.param2), + forwarding_policy()); + return policies::checked_narrowing_cast( + result, + function); + } private: // // Data members: diff --git a/include/boost/math/distributions/non_central_f.hpp b/include/boost/math/distributions/non_central_f.hpp index 9bae5c41ff..6a53a3f1d3 100644 --- a/include/boost/math/distributions/non_central_f.hpp +++ b/include/boost/math/distributions/non_central_f.hpp @@ -287,7 +287,7 @@ namespace boost template BOOST_MATH_GPU_ENABLED static RealType find_v1(const complemented4_type& c) { - constexpr auto function = "non_central_f_distribution<%1%>::find_non_centrality"; + constexpr auto function = "non_central_f_distribution<%1%>::find_v1"; typedef typename policies::evaluation::type eval_type; typedef typename policies::normalise< Policy, @@ -307,9 +307,9 @@ namespace boost result, function); } - BOOST_MATH_GPU_ENABLED static RealType find_v2(const RealType x, const RealType v2, const RealType nc, const RealType p) + 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"; + constexpr auto function = "non_central_f_distribution<%1%>::find_v2"; typedef typename policies::evaluation::type eval_type; typedef typename policies::normalise< Policy, @@ -332,7 +332,7 @@ namespace boost template BOOST_MATH_GPU_ENABLED static RealType find_v2(const complemented4_type& c) { - constexpr auto function = "non_central_f_distribution<%1%>::find_non_centrality"; + constexpr auto function = "non_central_f_distribution<%1%>::find_v2"; typedef typename policies::evaluation::type eval_type; typedef typename policies::normalise< Policy, diff --git a/test/test_fisher_f.cpp b/test/test_fisher_f.cpp index 0739dfbda2..37cabdafce 100644 --- a/test/test_fisher_f.cpp +++ b/test/test_fisher_f.cpp @@ -413,7 +413,19 @@ void test_spots(RealType) BOOST_CHECK_CLOSE( kurtosis(dist2) , static_cast(6272) * 12 / 3456 + 3, tol2); - // special cases: + + if (!std::is_same::value) + { + RealType df1 = 10; + RealType df2 = 2; + RealType x = 6; + RealType P = cdf(fisher_f_distribution(df1, df2), x); + BOOST_CHECK_CLOSE(fisher_f_distribution::find_v2(x, df1, P), df2, tol2*10); + BOOST_CHECK_CLOSE(fisher_f_distribution::find_v1(x, df2, P), df1, tol2*10); + BOOST_CHECK_CLOSE(fisher_f_distribution::find_v1(boost::math::complement(x, df2, 1-P)), df1, tol2*10); + BOOST_CHECK_CLOSE(fisher_f_distribution::find_v2(boost::math::complement(x, df1, 1-P)), df2, tol2*10); + } + // special cases: BOOST_MATH_CHECK_THROW( pdf( fisher_f_distribution(static_cast(1), static_cast(1)),