Skip to content

Find degrees of freedom for non_central_f distribution #1368

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
JacobHass8 wants to merge 3 commits into
base: develop
Choose a base branch
from
+178 −5
Open

Find degrees of freedom for non_central_f distribution #1368

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
44cbb53
Initial implementation/test of find_v1
JacobHass8 Feb 20, 2026
3eb7561
Added v2 finder
JacobHass8 Feb 21, 2026
4801654
Included edge cases for code coverage [ci skip]
JacobHass8 Feb 21, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
155 changes: 154 additions & 1 deletion include/boost/math/distributions/non_central_f.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,8 @@ namespace boost
};

template <class RealType, class Policy>
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) {
Expand Down Expand Up @@ -80,6 +81,68 @@ namespace boost
}
return result;
}

template <class RealType, class Policy>
struct f_degrees_of_freedom_finder
{
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& input_v)
{
RealType v1 = find_v1 ? input_v : v;
RealType v2 = find_v1 ? v : input_v;
non_central_f_distribution<RealType, Policy> d(v1, v2, nc);
return comp ?
p - cdf(complement(d, x))
: cdf(d, x) - p;
}
private:
RealType x;
RealType v;
RealType nc;
bool find_v1;
RealType p;
bool comp;
};

template <class RealType, class Policy>
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_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<RealType>(function, "Can't find degrees of freedom when the probability is 0 or 1, only possible answer is %1%",
RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy()); // LCOV_EXCL_LINE
}
if (x < tools::epsilon<RealType>())
{
return policies::raise_evaluation_error<RealType>(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<RealType>::quiet_NaN()), Policy()); // LCOV_EXCL_LINE
}
f_degrees_of_freedom_finder<RealType, Policy> f(x, v, nc, find_v1, p < q ? p : q, p < q ? false : true);
tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
//
// Pick an initial guess:
//
RealType guess = 1;
std::pair<RealType, RealType> 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<Policy>())
{
return policies::raise_evaluation_error<RealType>(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 RealType = double, class Policy = policies::policy<> >
Expand Down Expand Up @@ -163,6 +226,96 @@ 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<RealType, Policy>::type eval_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
eval_type result = detail::find_degrees_of_freedom_f(
static_cast<eval_type>(x),
static_cast<eval_type>(v2),
static_cast<eval_type>(nc),
true,
static_cast<eval_type>(p),
static_cast<eval_type>(1-p),
forwarding_policy());
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
result,
function);
}
template <class A, class B, class C, class D>
BOOST_MATH_GPU_ENABLED static RealType find_v1(const complemented4_type<A,B,C, D>& c)
{
constexpr auto function = "non_central_f_distribution<%1%>::find_non_centrality";
typedef typename policies::evaluation<RealType, Policy>::type eval_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
eval_type result = detail::find_degrees_of_freedom_f(
static_cast<eval_type>(c.dist),
static_cast<eval_type>(c.param1),
static_cast<eval_type>(c.param2),
true,
static_cast<eval_type>(1-c.param3),
static_cast<eval_type>(c.param3),
forwarding_policy());
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
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<RealType, Policy>::type eval_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
eval_type result = detail::find_degrees_of_freedom_f(
static_cast<eval_type>(x),
static_cast<eval_type>(v2),
static_cast<eval_type>(nc),
false,
static_cast<eval_type>(p),
static_cast<eval_type>(1-p),
forwarding_policy());
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
result,
function);
}
template <class A, class B, class C, class D>
BOOST_MATH_GPU_ENABLED static RealType find_v2(const complemented4_type<A,B,C, D>& c)
{
constexpr auto function = "non_central_f_distribution<%1%>::find_non_centrality";
typedef typename policies::evaluation<RealType, Policy>::type eval_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
eval_type result = detail::find_degrees_of_freedom_f(
static_cast<eval_type>(c.dist),
static_cast<eval_type>(c.param1),
static_cast<eval_type>(c.param2),
false,
static_cast<eval_type>(1-c.param3),
static_cast<eval_type>(c.param3),
forwarding_policy());
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
result,
function);
}
private:
// Data member, initialized by constructor.
RealType v1; // alpha.
Expand Down
28 changes: 24 additions & 4 deletions test/test_nc_f.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -138,13 +138,21 @@ 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);
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<RealType>() > 50)
{
Expand Down Expand Up @@ -369,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<long double>() / 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 <class RealType>void test_spots(RealType)

BOOST_AUTO_TEST_CASE( test_main )
Expand Down