420 lines · plain
1// boost\math\distributions\non_central_f.hpp2 3// Copyright John Maddock 2008.4// Copyright Matt Borland 2024.5// Use, modification and distribution are subject to the6// Boost Software License, Version 1.0.7// (See accompanying file LICENSE_1_0.txt8// or copy at http://www.boost.org/LICENSE_1_0.txt)9 10#ifndef BOOST_MATH_SPECIAL_NON_CENTRAL_F_HPP11#define BOOST_MATH_SPECIAL_NON_CENTRAL_F_HPP12 13#include <boost/math/tools/config.hpp>14#include <boost/math/tools/tuple.hpp>15#include <boost/math/tools/promotion.hpp>16#include <boost/math/distributions/non_central_beta.hpp>17#include <boost/math/distributions/detail/generic_mode.hpp>18#include <boost/math/special_functions/pow.hpp>19#include <boost/math/policies/policy.hpp>20 21namespace boost22{23 namespace math24 {25 template <class RealType = double, class Policy = policies::policy<> >26 class non_central_f_distribution27 {28 public:29 typedef RealType value_type;30 typedef Policy policy_type;31 32 BOOST_MATH_GPU_ENABLED non_central_f_distribution(RealType v1_, RealType v2_, RealType lambda) : v1(v1_), v2(v2_), ncp(lambda)33 {34 constexpr auto function = "boost::math::non_central_f_distribution<%1%>::non_central_f_distribution(%1%,%1%)";35 RealType r;36 detail::check_df(37 function,38 v1, &r, Policy());39 detail::check_df(40 function,41 v2, &r, Policy());42 detail::check_non_centrality(43 function,44 lambda,45 &r,46 Policy());47 } // non_central_f_distribution constructor.48 49 BOOST_MATH_GPU_ENABLED RealType degrees_of_freedom1()const50 {51 return v1;52 }53 BOOST_MATH_GPU_ENABLED RealType degrees_of_freedom2()const54 {55 return v2;56 }57 BOOST_MATH_GPU_ENABLED RealType non_centrality() const58 { // Private data getter function.59 return ncp;60 }61 private:62 // Data member, initialized by constructor.63 RealType v1; // alpha.64 RealType v2; // beta.65 RealType ncp; // non-centrality parameter66 }; // template <class RealType, class Policy> class non_central_f_distribution67 68 typedef non_central_f_distribution<double> non_central_f; // Reserved name of type double.69 70 #ifdef __cpp_deduction_guides71 template <class RealType>72 non_central_f_distribution(RealType,RealType,RealType)->non_central_f_distribution<typename boost::math::tools::promote_args<RealType>::type>;73 #endif74 75 // Non-member functions to give properties of the distribution.76 77 template <class RealType, class Policy>78 BOOST_MATH_GPU_ENABLED inline const boost::math::pair<RealType, RealType> range(const non_central_f_distribution<RealType, Policy>& /* dist */)79 { // Range of permissible values for random variable k.80 using boost::math::tools::max_value;81 return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());82 }83 84 template <class RealType, class Policy>85 BOOST_MATH_GPU_ENABLED inline const boost::math::pair<RealType, RealType> support(const non_central_f_distribution<RealType, Policy>& /* dist */)86 { // Range of supported values for random variable k.87 // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.88 using boost::math::tools::max_value;89 return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());90 }91 92 template <class RealType, class Policy>93 BOOST_MATH_GPU_ENABLED inline RealType mean(const non_central_f_distribution<RealType, Policy>& dist)94 {95 constexpr auto function = "mean(non_central_f_distribution<%1%> const&)";96 RealType v1 = dist.degrees_of_freedom1();97 RealType v2 = dist.degrees_of_freedom2();98 RealType l = dist.non_centrality();99 RealType r;100 if(!detail::check_df(101 function,102 v1, &r, Policy())103 ||104 !detail::check_df(105 function,106 v2, &r, Policy())107 ||108 !detail::check_non_centrality(109 function,110 l,111 &r,112 Policy()))113 return r;114 if(v2 <= 2)115 return policies::raise_domain_error(116 function,117 "Second degrees of freedom parameter was %1%, but must be > 2 !",118 v2, Policy());119 return v2 * (v1 + l) / (v1 * (v2 - 2));120 } // mean121 122 template <class RealType, class Policy>123 BOOST_MATH_GPU_ENABLED inline RealType mode(const non_central_f_distribution<RealType, Policy>& dist)124 { // mode.125 constexpr auto function = "mode(non_central_chi_squared_distribution<%1%> const&)";126 127 RealType n = dist.degrees_of_freedom1();128 RealType m = dist.degrees_of_freedom2();129 RealType l = dist.non_centrality();130 RealType r;131 if(!detail::check_df(132 function,133 n, &r, Policy())134 ||135 !detail::check_df(136 function,137 m, &r, Policy())138 ||139 !detail::check_non_centrality(140 function,141 l,142 &r,143 Policy()))144 return r;145 RealType guess = m > 2 ? RealType(m * (n + l) / (n * (m - 2))) : RealType(1);146 return detail::generic_find_mode(147 dist,148 guess,149 function);150 }151 152 template <class RealType, class Policy>153 BOOST_MATH_GPU_ENABLED inline RealType variance(const non_central_f_distribution<RealType, Policy>& dist)154 { // variance.155 constexpr auto function = "variance(non_central_f_distribution<%1%> const&)";156 RealType n = dist.degrees_of_freedom1();157 RealType m = dist.degrees_of_freedom2();158 RealType l = dist.non_centrality();159 RealType r;160 if(!detail::check_df(161 function,162 n, &r, Policy())163 ||164 !detail::check_df(165 function,166 m, &r, Policy())167 ||168 !detail::check_non_centrality(169 function,170 l,171 &r,172 Policy()))173 return r;174 if(m <= 4)175 return policies::raise_domain_error(176 function,177 "Second degrees of freedom parameter was %1%, but must be > 4 !",178 m, Policy());179 RealType result = 2 * m * m * ((n + l) * (n + l)180 + (m - 2) * (n + 2 * l));181 result /= (m - 4) * (m - 2) * (m - 2) * n * n;182 return result;183 }184 185 // RealType standard_deviation(const non_central_f_distribution<RealType, Policy>& dist)186 // standard_deviation provided by derived accessors.187 188 template <class RealType, class Policy>189 BOOST_MATH_GPU_ENABLED inline RealType skewness(const non_central_f_distribution<RealType, Policy>& dist)190 { // skewness = sqrt(l).191 constexpr auto function = "skewness(non_central_f_distribution<%1%> const&)";192 BOOST_MATH_STD_USING193 RealType n = dist.degrees_of_freedom1();194 RealType m = dist.degrees_of_freedom2();195 RealType l = dist.non_centrality();196 RealType r;197 if(!detail::check_df(198 function,199 n, &r, Policy())200 ||201 !detail::check_df(202 function,203 m, &r, Policy())204 ||205 !detail::check_non_centrality(206 function,207 l,208 &r,209 Policy()))210 return r;211 if(m <= 6)212 return policies::raise_domain_error(213 function,214 "Second degrees of freedom parameter was %1%, but must be > 6 !",215 m, Policy());216 RealType result = 2 * constants::root_two<RealType>();217 result *= sqrt(m - 4);218 result *= (n * (m + n - 2) *(m + 2 * n - 2)219 + 3 * (m + n - 2) * (m + 2 * n - 2) * l220 + 6 * (m + n - 2) * l * l + 2 * l * l * l);221 result /= (m - 6) * pow(n * (m + n - 2) + 2 * (m + n - 2) * l + l * l, RealType(1.5f));222 return result;223 }224 225 template <class RealType, class Policy>226 BOOST_MATH_GPU_ENABLED inline RealType kurtosis_excess(const non_central_f_distribution<RealType, Policy>& dist)227 {228 constexpr auto function = "kurtosis_excess(non_central_f_distribution<%1%> const&)";229 BOOST_MATH_STD_USING230 RealType n = dist.degrees_of_freedom1();231 RealType m = dist.degrees_of_freedom2();232 RealType l = dist.non_centrality();233 RealType r;234 if(!detail::check_df(235 function,236 n, &r, Policy())237 ||238 !detail::check_df(239 function,240 m, &r, Policy())241 ||242 !detail::check_non_centrality(243 function,244 l,245 &r,246 Policy()))247 return r;248 if(m <= 8)249 return policies::raise_domain_error(250 function,251 "Second degrees of freedom parameter was %1%, but must be > 8 !",252 m, Policy());253 RealType l2 = l * l;254 RealType l3 = l2 * l;255 RealType l4 = l2 * l2;256 RealType result = (3 * (m - 4) * (n * (m + n - 2)257 * (4 * (m - 2) * (m - 2)258 + (m - 2) * (m + 10) * n259 + (10 + m) * n * n)260 + 4 * (m + n - 2) * (4 * (m - 2) * (m - 2)261 + (m - 2) * (10 + m) * n262 + (10 + m) * n * n) * l + 2 * (10 + m)263 * (m + n - 2) * (2 * m + 3 * n - 4) * l2264 + 4 * (10 + m) * (-2 + m + n) * l3265 + (10 + m) * l4))266 /267 ((-8 + m) * (-6 + m) * boost::math::pow<2>(n * (-2 + m + n)268 + 2 * (-2 + m + n) * l + l2));269 return result;270 } // kurtosis_excess271 272 template <class RealType, class Policy>273 BOOST_MATH_GPU_ENABLED inline RealType kurtosis(const non_central_f_distribution<RealType, Policy>& dist)274 {275 return kurtosis_excess(dist) + 3;276 }277 278 template <class RealType, class Policy>279 BOOST_MATH_GPU_ENABLED inline RealType pdf(const non_central_f_distribution<RealType, Policy>& dist, const RealType& x)280 { // Probability Density/Mass Function.281 typedef typename policies::evaluation<RealType, Policy>::type value_type;282 typedef typename policies::normalise<283 Policy,284 policies::promote_float<false>,285 policies::promote_double<false>,286 policies::discrete_quantile<>,287 policies::assert_undefined<> >::type forwarding_policy;288 289 value_type alpha = dist.degrees_of_freedom1() / 2;290 value_type beta = dist.degrees_of_freedom2() / 2;291 value_type y = x * alpha / beta;292 value_type r = pdf(boost::math::non_central_beta_distribution<value_type, forwarding_policy>(alpha, beta, dist.non_centrality()), y / (1 + y));293 return policies::checked_narrowing_cast<RealType, forwarding_policy>(294 r * (dist.degrees_of_freedom1() / dist.degrees_of_freedom2()) / ((1 + y) * (1 + y)),295 "pdf(non_central_f_distribution<%1%>, %1%)");296 } // pdf297 298 template <class RealType, class Policy>299 BOOST_MATH_GPU_ENABLED RealType cdf(const non_central_f_distribution<RealType, Policy>& dist, const RealType& x)300 {301 constexpr auto function = "cdf(const non_central_f_distribution<%1%>&, %1%)";302 RealType r;303 if(!detail::check_df(304 function,305 dist.degrees_of_freedom1(), &r, Policy())306 ||307 !detail::check_df(308 function,309 dist.degrees_of_freedom2(), &r, Policy())310 ||311 !detail::check_non_centrality(312 function,313 dist.non_centrality(),314 &r,315 Policy()))316 return r;317 318 if((x < 0) || !(boost::math::isfinite)(x))319 {320 return policies::raise_domain_error<RealType>(321 function, "Random Variable parameter was %1%, but must be > 0 !", x, Policy());322 }323 324 RealType alpha = dist.degrees_of_freedom1() / 2;325 RealType beta = dist.degrees_of_freedom2() / 2;326 RealType y = x * alpha / beta;327 RealType c = y / (1 + y);328 RealType cp = 1 / (1 + y);329 //330 // To ensure accuracy, we pass both x and 1-x to the331 // non-central beta cdf routine, this ensures accuracy332 // even when we compute x to be ~ 1:333 //334 r = detail::non_central_beta_cdf(c, cp, alpha, beta,335 dist.non_centrality(), false, Policy());336 return r;337 } // cdf338 339 template <class RealType, class Policy>340 BOOST_MATH_GPU_ENABLED RealType cdf(const complemented2_type<non_central_f_distribution<RealType, Policy>, RealType>& c)341 { // Complemented Cumulative Distribution Function342 constexpr auto function = "cdf(complement(const non_central_f_distribution<%1%>&, %1%))";343 RealType r;344 if(!detail::check_df(345 function,346 c.dist.degrees_of_freedom1(), &r, Policy())347 ||348 !detail::check_df(349 function,350 c.dist.degrees_of_freedom2(), &r, Policy())351 ||352 !detail::check_non_centrality(353 function,354 c.dist.non_centrality(),355 &r,356 Policy()))357 return r;358 359 if((c.param < 0) || !(boost::math::isfinite)(c.param))360 {361 return policies::raise_domain_error<RealType>(362 function, "Random Variable parameter was %1%, but must be > 0 !", c.param, Policy());363 }364 365 RealType alpha = c.dist.degrees_of_freedom1() / 2;366 RealType beta = c.dist.degrees_of_freedom2() / 2;367 RealType y = c.param * alpha / beta;368 RealType x = y / (1 + y);369 RealType cx = 1 / (1 + y);370 //371 // To ensure accuracy, we pass both x and 1-x to the372 // non-central beta cdf routine, this ensures accuracy373 // even when we compute x to be ~ 1:374 //375 r = detail::non_central_beta_cdf(x, cx, alpha, beta,376 c.dist.non_centrality(), true, Policy());377 return r;378 } // ccdf379 380 template <class RealType, class Policy>381 BOOST_MATH_GPU_ENABLED inline RealType quantile(const non_central_f_distribution<RealType, Policy>& dist, const RealType& p)382 { // Quantile (or Percent Point) function.383 RealType alpha = dist.degrees_of_freedom1() / 2;384 RealType beta = dist.degrees_of_freedom2() / 2;385 RealType x = quantile(boost::math::non_central_beta_distribution<RealType, Policy>(alpha, beta, dist.non_centrality()), p);386 if(x == 1)387 return policies::raise_overflow_error<RealType>(388 "quantile(const non_central_f_distribution<%1%>&, %1%)",389 "Result of non central F quantile is too large to represent.",390 Policy());391 return (x / (1 - x)) * (dist.degrees_of_freedom2() / dist.degrees_of_freedom1());392 } // quantile393 394 template <class RealType, class Policy>395 BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type<non_central_f_distribution<RealType, Policy>, RealType>& c)396 { // Quantile (or Percent Point) function.397 RealType alpha = c.dist.degrees_of_freedom1() / 2;398 RealType beta = c.dist.degrees_of_freedom2() / 2;399 RealType x = quantile(complement(boost::math::non_central_beta_distribution<RealType, Policy>(alpha, beta, c.dist.non_centrality()), c.param));400 if(x == 1)401 return policies::raise_overflow_error<RealType>(402 "quantile(complement(const non_central_f_distribution<%1%>&, %1%))",403 "Result of non central F quantile is too large to represent.",404 Policy());405 return (x / (1 - x)) * (c.dist.degrees_of_freedom2() / c.dist.degrees_of_freedom1());406 } // quantile complement.407 408 } // namespace math409} // namespace boost410 411// This include must be at the end, *after* the accessors412// for this distribution have been defined, in order to413// keep compilers that support two-phase lookup happy.414#include <boost/math/distributions/detail/derived_accessors.hpp>415 416#endif // BOOST_MATH_SPECIAL_NON_CENTRAL_F_HPP417 418 419 420