diff --git a/doc/distributions/nc_f.qbk b/doc/distributions/nc_f.qbk index d31c8116bc..8611f5c674 100644 --- a/doc/distributions/nc_f.qbk +++ b/doc/distributions/nc_f.qbk @@ -26,6 +26,11 @@ // Accessor to non-centrality parameter lambda: BOOST_MATH_GPU_ENABLED RealType non_centrality()const; + + // Parameter finders: + static RealType find_non_centrality(const RealType x, const RealType v1, const RealType v2, const RealType p); + template + static RealType find_non_centrality(const complemented3_type& c); }; }} // namespaces @@ -53,6 +58,20 @@ for different values of [lambda]: [graph nc_f_pdf] + BOOST_MATH_GPU_ENABLED static RealType find_non_centrality(const RealType x, const RealType v1, const RealType v2, const RealType p); + +This function returns the non centrality parameter /lambda/ such that: + +`cdf(non_central_chi_squared(v1, v2, lambda), x) == p` + + template + BOOST_MATH_GPU_ENABLED static RealType find_non_centrality(const complemented4_type& c); + +When called with argument `boost::math::complement(x, v1, v2, q)` +this function returns the non centrality parameter /lambda/ such that: + +`cdf(complement(non_central_chi_squared(v1, v2, lambda), x)) == q`. + [h4 Member Functions] BOOST_MATH_GPU_ENABLED non_central_f_distribution(RealType v1, RealType v2, RealType lambda); diff --git a/include/boost/math/distributions/non_central_f.hpp b/include/boost/math/distributions/non_central_f.hpp index dedd437144..dfd9baf0d6 100644 --- a/include/boost/math/distributions/non_central_f.hpp +++ b/include/boost/math/distributions/non_central_f.hpp @@ -17,11 +17,60 @@ #include #include #include +#include // complements namespace boost { namespace math { + namespace detail + { + template + struct non_centrality_finder_f + { + non_centrality_finder_f(const RealType x_, const RealType v1_, const RealType v2_, const RealType p_, bool c) + : x(x_), v1(v1_), v2(v2_), p(p_), comp(c) {} + + RealType operator()(RealType nc) const + { + non_central_f_distribution d(v1, v2, nc); + return comp ? + RealType(p - cdf(complement(d, x))) + : RealType(cdf(d, x) - p); + } + private: + RealType x, v1, v2, p; + bool comp; + }; + + template + inline RealType find_non_centrality_f(const RealType x, const RealType v1, const RealType v2, const RealType p, const RealType q, const Policy& pol) + { + constexpr auto function = "non_central_f<%1%>::find_non_centrality"; + + if ( p == 0 || q == 0) { + return policies::raise_domain_error(function, "Can't find non centrality parameter when the probability is <=0 or >=1, only possible answer is %1%", // LCOV_EXCL_LINE + RealType(boost::math::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE + } + + 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. + boost::math::uintmax_t max_iter = policies::get_max_root_iterations(); + tools::eps_tolerance tol(policies::digits()); + + std::pair result_bracket = tools::bracket_and_solve_root( + f, guess, factor, false, tol, max_iter, pol); + + RealType result = result_bracket.first + (result_bracket.second - result_bracket.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 > class non_central_f_distribution { @@ -58,6 +107,49 @@ namespace boost { // Private data getter function. return ncp; } + static RealType find_non_centrality(const RealType x, const RealType v1, const RealType v2, const RealType p) + { + 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_non_centrality_f( + static_cast(x), + static_cast(v1), + static_cast(v2), + static_cast(p), + static_cast(1-p), + forwarding_policy()); + return policies::checked_narrowing_cast( + result, + function); + } + template + static RealType find_non_centrality(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_non_centrality_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. @@ -404,7 +496,6 @@ namespace boost Policy()); return (x / (1 - x)) * (c.dist.degrees_of_freedom2() / c.dist.degrees_of_freedom1()); } // quantile complement. - } // namespace math } // namespace boost diff --git a/test/test_nc_f.cpp b/test/test_nc_f.cpp index 3a582a29a1..fdd03ab287 100644 --- a/test/test_nc_f.cpp +++ b/test/test_nc_f.cpp @@ -141,6 +141,10 @@ void test_spot( quantile(dist, P), x, tol * 10); BOOST_CHECK_CLOSE( quantile(complement(dist, Q)), x, tol * 10); + BOOST_CHECK_CLOSE( + 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); } if(boost::math::tools::digits() > 50) {