boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-02-24 04:02:18 +00:00
Code Issues Projects Releases Wiki Activity

Accepted modifications for nc close to 0

Browse Source
This commit is contained in:
Jacob Hass
2026-01-27 12:51:55 -08:00
parent 59a4e7a46d
commit 03bfa3362d
2 changed files with 21 additions and 19 deletions
Download Patch File Download Diff File Expand all files Collapse all files

22
include/boost/math/distributions/non_central_f.hpp
View File

@@ -46,8 +46,7 @@ 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 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) {
@@ -56,17 +55,20 @@ namespace boost
}
// Check if nc = 0 (which is just the F-distribution)
// I first tried comparing to boost::math::tools::epsilon<RealType>(),
// but this was never satisfied for floats or doubles. Only
// long doubles would correctly return 0.
fisher_f_distribution<RealType, Policy> dist(v1, v2);
if (boost::math::relative_difference(pdf(dist, x), p) < 1e-7){
non_centrality_finder_f<RealType, Policy> f(x, v1, v2, p < q ? p : q, p < q ? false : true);
// This occurs when the root finder would need to find a result smaller than
// tools::min_value (which it cannot do). Note that we have to add in a small
// amount of "tolerance" since the subtraction in our termination condition
// implies a small amount of wobble in the result which should be of the
// order p * eps (note not q * eps, since q is calculated as 1-p).
// Also note that p_q_precision is passed down from our caller as the
// epsilon of the original called values, and not after possible promotion.
if (f(tools::min_value<RealType>()) <= 3 * p_q_precision * p){
return 0;
}
non_centrality_finder_f<RealType, Policy> f(x, v1, v2, p < q ? p : q, p < q ? false : true);
RealType guess = RealType(10); // Starting guess.
RealType factor = RealType(2); // 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<Policy>();
tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
@@ -134,6 +136,7 @@ namespace boost
static_cast<eval_type>(v2),
static_cast<eval_type>(p),
static_cast<eval_type>(1-p),
static_cast<eval_type>(tools::epsilon<RealType>()),
forwarding_policy());
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
result,
@@ -156,6 +159,7 @@ namespace boost
static_cast<eval_type>(c.param2),
static_cast<eval_type>(1-c.param3),
static_cast<eval_type>(c.param3),
static_cast<eval_type>(tools::epsilon<RealType>()),
forwarding_policy());
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
result,

18
test/test_nc_f.cpp
View File

@@ -159,12 +159,10 @@ void test_spot(
}
template <class RealType> // Any floating-point type RealType.
void test_spots(RealType)
void test_spots(RealType, const char* name = nullptr)
{
RealType tolerance = boost::math::tools::epsilon<RealType>() * 10000;
cout << "Tolerance = " << (tolerance / 100) << "%." << endl;
std::cout << "Testing spot values with type " << name << " (Tolerance = " << (tolerance / 100) << "%)." << std::endl;
//
// Spot tests from Mathematica computed values:
//
@@ -337,8 +335,8 @@ void test_spots(RealType)
boost::math::non_central_f_distribution<RealType> 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<RealType>(0), tolerance);
RealType P = cdf(dist_no_centrality, x);
BOOST_CHECK(dist.find_non_centrality(x, a, b, P) < tolerance);
}
// Case when P=1 or P=0
BOOST_MATH_CHECK_THROW(dist.find_non_centrality(x, a, b, 1), std::domain_error);
@@ -354,12 +352,12 @@ BOOST_AUTO_TEST_CASE( test_main )
// Basic sanity-check spot values.
expected_results();
// (Parameter value, arbitrarily zero, only communicates the floating point type).
test_spots(0.0F); // Test float.
test_spots(0.0); // Test double.
test_spots(0.0F, "float"); // Test float.
test_spots(0.0, "double"); // Test double.
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
test_spots(0.0L); // Test long double.
test_spots(0.0L, "long double"); // Test long double.
#if !BOOST_WORKAROUND(BOOST_BORLANDC, BOOST_TESTED_AT(0x582)) && !defined(BOOST_MATH_NO_REAL_CONCEPT_TESTS)
test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
test_spots(boost::math::concepts::real_concept(0.), "real_concept"); // Test real concept.
#endif
#endif