Initial implementation of find_non_centrality for non_central_f_distribution

[JacobHass8]

See #1305 for for details.

[JacobHass8]

@dschmitz89 how does this first commit look in regards to #1338 (comment)

[dschmitz89]

Thanks @JacobHass8 for looking into this so quickly. This is a good start. I am not a boost maintainer (barely a boost developer), so can only comment from a high level.

The most important missing bit here is tests for the new functionality. I could not find the current noncentral F distribution tests after 2 minutes though, so we will have to rely on the boost devs to guide us a little here.

[jzmaddock]

Yup, tests needed.

Non central F is treated as a special case of the non-central beta (which is where most of the tests are). Non-central F tests are fairly lightweight, just there to make sure we have the transformation to non-central beta correct really: https://github.com/boostorg/math/blob/develop/test/test_nc_f.cpp. A similar level of sanity checking would probably suffice for this function I would guess?

[JacobHass8]

Yup, tests needed.

Non central F is treated as a special case of the non-central beta (which is where most of the tests are). Non-central F tests are fairly lightweight, just there to make sure we have the transformation to non-central beta correct really: https://github.com/boostorg/math/blob/develop/test/test_nc_f.cpp. A similar level of sanity checking would probably suffice for this function I would guess?

I've added a couple of tests and documentation. The tests and documentation copy what is done for the non_central_chi_squared distribution. Is there anything else I need to do before this can be merged?

[JacobHass8]

@jzmaddock could I get a review on this when you get the chance? Are there any test cases that you can think of that I should add?

[JacobHass8]

@mborland (#1305) if you're interested, here's my first attempt at the inverse of the non-central F distribution. It was pretty much copied from the scipy implementation.

[mborland]

This looks pretty good to me. I made a couple minor changes to get it running with CUDA. Since scipy already has an implementation, it's probably worth copying over any tests they have as well.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 95.30%. Comparing base (1d670c8) to head (e51d71a).

Additional details and impacted files

@@           Coverage Diff            @@
##           develop    #1345   +/-   ##
========================================
  Coverage    95.29%   95.30%           
========================================
  Files          814      814           
  Lines        67422    67461   +39     
========================================
+ Hits         64249    64292   +43     
+ Misses        3173     3169    -4     

Files with missing lines	Coverage Δ
include/boost/math/distributions/non_central_f.hpp	87.76% <100.00%> (+2.26%)	⬆️
test/test_nc_f.cpp	100.00% <100.00%> (ø)

... and 1 file with indirect coverage changes

---

Continue to review full report in Codecov by Sentry.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 1d670c8...e51d71a. Read the comment docs.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[dschmitz89]

This looks pretty good to me. I made a couple minor changes to get it running with CUDA. Since scipy already has an implementation, it's probably worth copying over any tests they have as well.

SciPy does not have tests for these functions at the moment yet. The functions were simply incorporated a long time ago with the cdflib library but actually never unit tested. One more plus point for the migration we are planning.

[JacobHass8]

SciPy does not have tests for these functions at the moment yet. The functions were simply incorporated a long time ago with the cdflib library but actually never unit tested. One more plus point for the migration we are planning.

There's actually a couple here. It's only really 3 spot checks.

I could add some edge case testing like what happens if $p=0,1$ or $nc = \infty$. I could also do more spot checks for small and large values of $nc$ if that would help.

[dschmitz89]

SciPy does not have tests for these functions at the moment yet. The functions were simply incorporated a long time ago with the cdflib library but actually never unit tested. One more plus point for the migration we are planning.

There's actually a couple here. It's only really 3 spot checks.

I could add some edge case testing like what happens if p = 0 , 1 or n c = ∞ . I could also do more spot checks for small and large values of n c if that would help.

Ah, thanks for spotting those. The file where the function is supposed to be tested, lists it as untested.

About tests here: as long as boost maintainers are happy, it should be good. In SciPy, we do not need to test the function very extensively if it is already thoroughly tested in boost (although one never knows in future, SciPy has demanding users .. ).

[mborland]

I could add some edge case testing like what happens if p = 0 , 1 or n c = ∞ . I could also do more spot checks for small and large values of n c if that would help.

I think adding all of these makes sense. As @dschmitz89 alluded to, the SciPy users are pretty adept at finding weakness in our edge cases when they exist.

[JacobHass8]

All the tests are passing except for g++-14 c++23 autodiff which appears to be unrelated to this PR. The error is below.

Details

unknown location(0): fatal error: in "test_autodiff_8/heuman_lambda_hpp<_Float32>": Throw location unknown (consider using BOOST_THROW_EXCEPTION)
Dynamic exception type: boost::wrapexcept<std::domain_error>
std::exception::what: Error in function boost::math::heuman_lambda<N5boost4math15differentiation11autodiff_v16detail4fvarIDF32_Lm5EEE>(N5boost4math15differentiation11autodiff_v16detail4fvarIDF32_Lm5EEE, N5boost4math15differentiation11autodiff_v16detail4fvarIDF32_Lm5EEE): When 1-k^2 == 1 then phi must be < Pi/2, but got phi = depth(1)(4.20396852,1,0,0,0,0)

test_autodiff_8.cpp(38): last checkpoint

*** 1 failure is detected in the test module "test_autodiff"

[JacobHass8]

I think adding all of these makes sense. As @dschmitz89 alluded to, the SciPy users are pretty adept at finding weakness in our edge cases when they exist.

I've added tests to make sure that $p=0, 1$ and $q=0,1$ throw errors. The spot tests already included have large values of nc. I had to do some special handling when $nc = 0$ and the non-central distribution reduces to the F-distribution. This is done in the following lines

fisher_f_distribution<RealType, Policy> dist(v1, v2); 
if (boost::math::relative_difference(pdf(dist, x), p) < 1e-7){
   return 0;
}

I first tried a relative difference of boost::math::tools::epsilon<RealType>(), but this only worked for long doubles. I've added a number of spot checks for when $nc=0$, and this condition works in those cases.

[mborland]

I think adding all of these makes sense. As @dschmitz89 alluded to, the SciPy users are pretty adept at finding weakness in our edge cases when they exist.

I've added tests to make sure that p = 0 , 1 and q = 0 , 1 throw errors. The spot tests already included have large values of nc. I had to do some special handling when n c = 0 and the non-central distribution reduces to the F-distribution. This is done in the following lines

fisher_f_distribution<RealType, Policy> dist(v1, v2); 
if (boost::math::relative_difference(pdf(dist, x), p) < 1e-7){
   return 0;
}

I first tried a relative difference of boost::math::tools::epsilon<RealType>(), but this only worked for long doubles. I've added a number of spot checks for when n c = 0 , and this condition works in those cases.

There is also boost::math::epsilon_difference which would allow you to set a certain eps threshold. That might give more consistency across type since the relative distance tolerance should be dependent on precision of the type. Then you could have something like:

fisher_f_distribution<RealType, Policy> dist(v1, v2); 
if (boost::math::epsilon_difference(pdf(dist, x), p) <= 1){
    return 0;
}

float_distance could also work but I believe it's more expensive to calculate:

fisher_f_distribution<RealType, Policy> dist(v1, v2); 
if (abs(boost::math::float_distance(pdf(dist, x), p)) <= 1){
    return 0;
}

See also: https://www.boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/float_comparison.html

[JacobHass8]

There is also boost::math::epsilon_difference which would allow you to set a certain eps threshold. That might give more consistency across type since the relative distance tolerance should be dependent on precision of the type. Then you could have something like:

fisher_f_distribution<RealType, Policy> dist(v1, v2); 
if (boost::math::epsilon_difference(pdf(dist, x), p) <= 1){
    return 0;
}

float_distance could also work but I believe it's more expensive to calculate:

fisher_f_distribution<RealType, Policy> dist(v1, v2); 
if (abs(boost::math::float_distance(pdf(dist, x), p)) <= 1){
    return 0;
}

Hmmm, neither of these conditions worked either. For floats, abs(boost::math::float_distance(pdf(dist, x), p)) and boost::math::epsilon_difference(pdf(dist, x), p) were both of order 1e8.

I also tried

tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
fisher_f_distribution<RealType, Policy> dist(v1, v2); 
if (tol(pdf(dist, x), p)){
    return 0;
}

and that didn't work either. The relative difference decreases with increased precision though.

[JacobHass8]

I know what is going wrong now. When calculating the pdf of dist for floats:

fisher_f_distribution<RealType, Policy> dist(v1, v2); 
if (abs(boost::math::float_distance(pdf(dist, x), p)) <= 1){
    return 0;
}

the template RealType is upcast to a long double (maybe a double, I'm not 100% sure). If I change the dist variable above to non_central_f_distribution<float, policies::policy<> > dist(v1, v2, 0);, and take the relative distance I get 0 as expected. I also think we should be comparing to non_central_f_distribution with nc=0 rather than fisher_f_distribution because they calculate the pdf differently.

I'm not sure if there is an elegant way to avoid this upcasting based on type because this happens before find_non_centrality_f is called.

[mborland]

I know what is going wrong now. When calculating the pdf of dist for floats:

fisher_f_distribution<RealType, Policy> dist(v1, v2); 
if (abs(boost::math::float_distance(pdf(dist, x), p)) <= 1){
    return 0;
}

the template RealType is upcast to a long double (maybe a double, I'm not 100% sure). If I change the dist variable above to non_central_f_distribution<float, policies::policy<> > dist(v1, v2, 0);, and take the relative distance I get 0 as expected. I also think we should be comparing to non_central_f_distribution with nc=0 rather than fisher_f_distribution because they calculate the pdf differently.

I'm not sure if there is an elegant way to avoid this upcasting based on type because this happens before find_non_centrality_f is called.

Instead of using policies::policy<> as your template parameter to non_central_f_distribution you should be able to use policies::policy<policies::promote_float<false>> in order to avoid this.

[JacobHass8]

Instead of using policies::policy<> as your template parameter to non_central_f_distribution you should be able to use policies::policy<policies::promote_float<false>> in order to avoid this.

Not quite, because RealType at this point has already been upcast to a long double from a float. If we just stick with a relative precision less than 1e-7, then I don't think we'll have to deal with all the up/down casting based on type/policy.

[mborland]

Instead of using policies::policy<> as your template parameter to non_central_f_distribution you should be able to use policies::policy<policies::promote_float<false>> in order to avoid this.

Not quite, because RealType at this point has already been upcast to a long double from a float. If we just stick with a relative precision less than 1e-7, then I don't think we'll have to deal with all the up/down casting based on type/policy.

Ah, you're correct. @jzmaddock do you have any final comments for this one before merge?