Initial implementation of find_non_centrality for non_central_f_distribution

CMakeLists.txt

@@ -39,7 +39,6 @@ else()
       Boost::lexical_cast
       Boost::predef
       Boost::random
-      Boost::static_assert
       Boost::throw_exception
   )
 

build.jam

@@ -14,7 +14,6 @@ constant boost_dependencies :
     /boost/lexical_cast//boost_lexical_cast
     /boost/predef//boost_predef
     /boost/random//boost_random
-    /boost/static_assert//boost_static_assert
     /boost/throw_exception//boost_throw_exception ;
 
 project /boost/math

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: 
+      BOOST_MATH_GPU_ENABLED static RealType find_non_centrality(const RealType x, const RealType v1, const RealType v2, const RealType p);
+      template <class A, class B, class C, class D>
+      BOOST_MATH_GPU_ENABLED static RealType find_non_centrality(const complemented3_type<A,B,C,D>& 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<RealType, Policy>(v1, v2, lambda), x) == p`
+
+      template <class A, class B, class C, class D>
+      BOOST_MATH_GPU_ENABLED static RealType find_non_centrality(const complemented4_type<A,B,C,D>& 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<RealType, Policy>(v1, v2, lambda), x)) == q`.
+
 [h4 Member Functions]
 
       BOOST_MATH_GPU_ENABLED non_central_f_distribution(RealType v1, RealType v2, RealType lambda);

include/boost/math/distributions/non_central_f.hpp

@@ -17,11 +17,73 @@
 #include <boost/math/distributions/detail/generic_mode.hpp>
 #include <boost/math/special_functions/pow.hpp>
 #include <boost/math/policies/policy.hpp>
+#include <boost/math/distributions/complement.hpp> // complements
+#include <boost/math/distributions/fisher_f.hpp>
+#include <boost/math/special_functions/relative_difference.hpp> 
 
 namespace boost
 {
    namespace math
    {
+      namespace detail
+      {
+         template <class RealType, class Policy>
+         struct non_centrality_finder_f
+         {
+            BOOST_MATH_GPU_ENABLED non_centrality_finder_f(const RealType x_, const RealType v1_, const RealType v2_, const RealType p_, const bool c)
+               : x(x_), v1(v1_), v2(v2_), p(p_), comp(c) {}
+
+            BOOST_MATH_GPU_ENABLED RealType operator()(RealType nc) const
+            {
+               non_central_f_distribution<RealType, Policy> 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 <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)         {
+            constexpr auto function = "non_central_f<%1%>::find_non_centrality";
+
+            if ( p == 0 || q == 0) {
+               return policies::raise_domain_error<RealType>(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<RealType>::quiet_NaN()), Policy()); // LCOV_EXCL_LINE
+            }
+
+            // Check if nc = 0 (which is just the F-distribution)
+            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>()) <= 20 * p_q_precision * p){
+               return 0;
+            }
+
+            RealType guess = RealType(10);                       // Starting guess.
+            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>());
+
+            boost::math::pair<RealType, RealType> 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<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<> >
       class non_central_f_distribution
       {
@@ -58,6 +120,51 @@ namespace boost
          { // Private data getter function.
             return ncp;
          }
+         BOOST_MATH_GPU_ENABLED 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<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_non_centrality_f(
+               static_cast<eval_type>(x),
+               static_cast<eval_type>(v1),
+               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,
+               function);
+         }
+         template <class A, class B, class C, class D>
+         BOOST_MATH_GPU_ENABLED static RealType find_non_centrality(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_non_centrality_f(
+               static_cast<eval_type>(c.dist),
+               static_cast<eval_type>(c.param1),
+               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,
+               function);
+         }
       private:
          // Data member, initialized by constructor.
          RealType v1;   // alpha.
@@ -404,7 +511,6 @@ namespace boost
                Policy());
          return (x / (1 - x)) * (c.dist.degrees_of_freedom2() / c.dist.degrees_of_freedom1());
       } // quantile complement.
-
    } // namespace math
 } // namespace boost
 

test/CMakeLists.txt

@@ -16,7 +16,7 @@ if(HAVE_BOOST_TEST)
 
         enable_testing()
 
-        boost_test_jamfile(FILE cuda_jamfile LINK_LIBRARIES Boost::math Boost::assert Boost::concept_check Boost::config Boost::core Boost::integer Boost::lexical_cast Boost::multiprecision Boost::predef Boost::random Boost::static_assert Boost::throw_exception Boost::unit_test_framework ${CUDA_LIBRARIES} INCLUDE_DIRECTORIES ${CUDA_INCLUDE_DIRS}  )
+        boost_test_jamfile(FILE cuda_jamfile LINK_LIBRARIES Boost::math Boost::assert Boost::concept_check Boost::config Boost::core Boost::integer Boost::lexical_cast Boost::multiprecision Boost::predef Boost::random Boost::throw_exception Boost::unit_test_framework ${CUDA_LIBRARIES} INCLUDE_DIRECTORIES ${CUDA_INCLUDE_DIRS}  )
 
     elseif (BOOST_MATH_ENABLE_NVRTC)
 
@@ -29,9 +29,9 @@ if(HAVE_BOOST_TEST)
         set(CUDA_nvrtc_LIBRARY /usr/local/cuda/lib64/libnvrtc.so)
 
         if (BOOST_MATH_NVRTC_CI_RUN)
-            boost_test_jamfile(FILE nvrtc_jamfile LINK_LIBRARIES Boost::math Boost::assert Boost::concept_check Boost::config Boost::core Boost::integer Boost::lexical_cast Boost::multiprecision Boost::predef Boost::random Boost::static_assert Boost::throw_exception ${CUDA_nvrtc_LIBRARY} ${CUDA_LIBRARIES} ${CUDA_CUDA_LIBRARY} COMPILE_DEFINITIONS BOOST_MATH_NVRTC_CI_RUN=1 INCLUDE_DIRECTORIES ${CUDA_INCLUDE_DIRS} )
+            boost_test_jamfile(FILE nvrtc_jamfile LINK_LIBRARIES Boost::math Boost::assert Boost::concept_check Boost::config Boost::core Boost::integer Boost::lexical_cast Boost::multiprecision Boost::predef Boost::random Boost::throw_exception ${CUDA_nvrtc_LIBRARY} ${CUDA_LIBRARIES} ${CUDA_CUDA_LIBRARY} COMPILE_DEFINITIONS BOOST_MATH_NVRTC_CI_RUN=1 INCLUDE_DIRECTORIES ${CUDA_INCLUDE_DIRS} )
         else ()
-            boost_test_jamfile(FILE nvrtc_jamfile LINK_LIBRARIES Boost::math Boost::assert Boost::concept_check Boost::config Boost::core Boost::integer Boost::lexical_cast Boost::multiprecision Boost::predef Boost::random Boost::static_assert Boost::throw_exception ${CUDA_nvrtc_LIBRARY} ${CUDA_LIBRARIES} ${CUDA_CUDA_LIBRARY} INCLUDE_DIRECTORIES ${CUDA_INCLUDE_DIRS} )
+            boost_test_jamfile(FILE nvrtc_jamfile LINK_LIBRARIES Boost::math Boost::assert Boost::concept_check Boost::config Boost::core Boost::integer Boost::lexical_cast Boost::multiprecision Boost::predef Boost::random Boost::throw_exception ${CUDA_nvrtc_LIBRARY} ${CUDA_LIBRARIES} ${CUDA_CUDA_LIBRARY} INCLUDE_DIRECTORIES ${CUDA_INCLUDE_DIRS} )
         endif()
 
     elseif (BOOST_MATH_ENABLE_SYCL)
@@ -43,7 +43,7 @@ if(HAVE_BOOST_TEST)
 
         enable_testing()
 
-        boost_test_jamfile(FILE sycl_jamfile LINK_LIBRARIES Boost::math Boost::assert Boost::concept_check Boost::config Boost::core Boost::integer Boost::lexical_cast Boost::multiprecision Boost::predef Boost::random Boost::static_assert Boost::throw_exception Boost::unit_test_framework sycl COMPILE_DEFINITIONS BOOST_MATH_ENABLE_SYCL=1 COMPILE_OPTIONS -fsycl )
+        boost_test_jamfile(FILE sycl_jamfile LINK_LIBRARIES Boost::math Boost::assert Boost::concept_check Boost::config Boost::core Boost::integer Boost::lexical_cast Boost::multiprecision Boost::predef Boost::random Boost::throw_exception Boost::unit_test_framework sycl COMPILE_DEFINITIONS BOOST_MATH_ENABLE_SYCL=1 COMPILE_OPTIONS -fsycl )
     else()
 
         boost_test(SOURCES check_cmake_version.cpp ARGUMENTS ${PROJECT_VERSION} LINK_LIBRARIES Boost::core Boost::config)

test/cuda_jamfile

@@ -165,6 +165,8 @@ run test_nc_f_pdf_double.cu ;
 run test_nc_f_pdf_float.cu ;
 run test_nc_f_quan_double.cu ;
 run test_nc_f_quan_float.cu ;
+run test_nc_f_finder_double.cu ;
+run test_nc_f_finder_float.cu ;
 
 run test_nc_chi_squared_cdf_double.cu ;
 run test_nc_chi_squared_cdf_float.cu ;

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<RealType>() > 50)
    {
@@ -155,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:
    //
@@ -323,6 +325,50 @@ void test_spots(RealType)
    {
       BOOST_CHECK_CLOSE(cdf(boost::math::non_central_f_distribution<RealType>(static_cast<RealType>(1e-100L), 3.f, 1.5f), static_cast<RealType>(1e100L)), static_cast<RealType>(0.6118152873453990639132215575213809716459L), tolerance);
    }
+
+   // Check find_non_centrality_f edge case handling 
+   // Case when nc=0
+   RealType a = 5;
+   RealType b = 2;
+   RealType nc = 0;
+   RealType x_vals[] = { 0.25, 1.25, 10, 100};
+   boost::math::non_central_f_distribution<RealType> dist_no_centrality(a, b, nc);
+   for (RealType x : x_vals)
+   {
+      RealType P = cdf(dist_no_centrality, x);
+      BOOST_CHECK_LE(dist.find_non_centrality(x, a, b, P), boost::math::tools::min_value<RealType>() * 10);
+   }
+   // Case when P=1 or P=0 
+   BOOST_MATH_CHECK_THROW(dist.find_non_centrality(x, a, b, 1), std::domain_error);
+   BOOST_MATH_CHECK_THROW(dist.find_non_centrality(x, a, b, 0), std::domain_error);
+   // Case when Q=1 or Q=0
+   BOOST_MATH_CHECK_THROW(dist.find_non_centrality(boost::math::complement(x, a, b, 1)), std::domain_error);
+   BOOST_MATH_CHECK_THROW(dist.find_non_centrality(boost::math::complement(x, a, b, 0)), std::domain_error);
+   //
+   // Test non centrality finder over a grid of values:
+   //
+   RealType values[] = { 1.25, 3.5, 6.75, 8.25 };
+   for (RealType v1 : values)
+   {
+      for (RealType v2 : values)
+      {
+         for (RealType nc : values)
+         {
+            for (RealType x : values)
+            {
+               boost::math::non_central_f_distribution<RealType> ref(v1, v2, nc);
+               RealType P = cdf(ref, x);
+               RealType Q = cdf(complement(ref, x));
+
+               RealType nc1 = ref.find_non_centrality(x, v1, v2, P);
+               RealType nc2 = ref.find_non_centrality(boost::math::complement(x, v1, v2, Q));
+
+               BOOST_CHECK_CLOSE(nc1, nc, 2 * tolerance);
+               BOOST_CHECK_CLOSE(nc2, nc, tolerance);
+            }
+         }
+      }
+   }
 } // template <class RealType>void test_spots(RealType)
 
 BOOST_AUTO_TEST_CASE( test_main )
@@ -331,12 +377,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