394 lines · plain
1// Copyright John Maddock 2006.2// Copyright Matt Borland 2024.3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0.5// (See accompanying file LICENSE_1_0.txt6// or copy at http://www.boost.org/LICENSE_1_0.txt)7 8#ifndef BOOST_MATH_DISTRIBUTIONS_FISHER_F_HPP9#define BOOST_MATH_DISTRIBUTIONS_FISHER_F_HPP10 11#include <boost/math/tools/config.hpp>12#include <boost/math/tools/tuple.hpp>13#include <boost/math/tools/promotion.hpp>14#include <boost/math/distributions/fwd.hpp>15#include <boost/math/special_functions/beta.hpp> // for incomplete beta.16#include <boost/math/distributions/complement.hpp> // complements17#include <boost/math/distributions/detail/common_error_handling.hpp> // error checks18#include <boost/math/special_functions/fpclassify.hpp>19 20namespace boost{ namespace math{21 22template <class RealType = double, class Policy = policies::policy<> >23class fisher_f_distribution24{25public:26 typedef RealType value_type;27 typedef Policy policy_type;28 29 BOOST_MATH_GPU_ENABLED fisher_f_distribution(const RealType& i, const RealType& j) : m_df1(i), m_df2(j)30 {31 constexpr auto function = "fisher_f_distribution<%1%>::fisher_f_distribution";32 RealType result;33 detail::check_df(34 function, m_df1, &result, Policy());35 detail::check_df(36 function, m_df2, &result, Policy());37 } // fisher_f_distribution38 39 BOOST_MATH_GPU_ENABLED RealType degrees_of_freedom1()const40 {41 return m_df1;42 }43 BOOST_MATH_GPU_ENABLED RealType degrees_of_freedom2()const44 {45 return m_df2;46 }47 48private:49 //50 // Data members:51 //52 RealType m_df1; // degrees of freedom are a real number.53 RealType m_df2; // degrees of freedom are a real number.54};55 56typedef fisher_f_distribution<double> fisher_f;57 58#ifdef __cpp_deduction_guides59template <class RealType>60fisher_f_distribution(RealType,RealType)->fisher_f_distribution<typename boost::math::tools::promote_args<RealType>::type>;61#endif62 63template <class RealType, class Policy>64BOOST_MATH_GPU_ENABLED inline const boost::math::pair<RealType, RealType> range(const fisher_f_distribution<RealType, Policy>& /*dist*/)65{ // Range of permissible values for random variable x.66 using boost::math::tools::max_value;67 return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());68}69 70template <class RealType, class Policy>71BOOST_MATH_GPU_ENABLED inline const boost::math::pair<RealType, RealType> support(const fisher_f_distribution<RealType, Policy>& /*dist*/)72{ // Range of supported values for random variable x.73 // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.74 using boost::math::tools::max_value;75 return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());76}77 78template <class RealType, class Policy>79BOOST_MATH_GPU_ENABLED RealType pdf(const fisher_f_distribution<RealType, Policy>& dist, const RealType& x)80{81 BOOST_MATH_STD_USING // for ADL of std functions82 RealType df1 = dist.degrees_of_freedom1();83 RealType df2 = dist.degrees_of_freedom2();84 // Error check:85 RealType error_result = 0;86 constexpr auto function = "boost::math::pdf(fisher_f_distribution<%1%> const&, %1%)";87 if(false == (detail::check_df(88 function, df1, &error_result, Policy())89 && detail::check_df(90 function, df2, &error_result, Policy())))91 return error_result;92 93 if((x < 0) || !(boost::math::isfinite)(x))94 {95 return policies::raise_domain_error<RealType>(96 function, "Random variable parameter was %1%, but must be > 0 !", x, Policy());97 }98 99 if(x == 0)100 {101 // special cases:102 if(df1 < 2)103 return policies::raise_overflow_error<RealType>(104 function, 0, Policy());105 else if(df1 == 2)106 return 1;107 else108 return 0;109 }110 111 //112 // You reach this formula by direct differentiation of the113 // cdf expressed in terms of the incomplete beta.114 //115 // There are two versions so we don't pass a value of z116 // that is very close to 1 to ibeta_derivative: for some values117 // of df1 and df2, all the change takes place in this area.118 //119 RealType v1x = df1 * x;120 RealType result;121 if(v1x > df2)122 {123 result = (df2 * df1) / ((df2 + v1x) * (df2 + v1x));124 result *= ibeta_derivative(df2 / 2, df1 / 2, df2 / (df2 + v1x), Policy());125 }126 else127 {128 result = df2 + df1 * x;129 result = (result * df1 - x * df1 * df1) / (result * result);130 result *= ibeta_derivative(df1 / 2, df2 / 2, v1x / (df2 + v1x), Policy());131 }132 return result;133} // pdf134 135template <class RealType, class Policy>136BOOST_MATH_GPU_ENABLED inline RealType cdf(const fisher_f_distribution<RealType, Policy>& dist, const RealType& x)137{138 constexpr auto function = "boost::math::cdf(fisher_f_distribution<%1%> const&, %1%)";139 RealType df1 = dist.degrees_of_freedom1();140 RealType df2 = dist.degrees_of_freedom2();141 // Error check:142 RealType error_result = 0;143 if(false == detail::check_df(144 function, df1, &error_result, Policy())145 && detail::check_df(146 function, df2, &error_result, Policy()))147 return error_result;148 149 if((x < 0) || !(boost::math::isfinite)(x))150 {151 return policies::raise_domain_error<RealType>(152 function, "Random Variable parameter was %1%, but must be > 0 !", x, Policy());153 }154 155 RealType v1x = df1 * x;156 //157 // There are two equivalent formulas used here, the aim is158 // to prevent the final argument to the incomplete beta159 // from being too close to 1: for some values of df1 and df2160 // the rate of change can be arbitrarily large in this area,161 // whilst the value we're passing will have lost information162 // content as a result of being 0.999999something. Better163 // to switch things around so we're passing 1-z instead.164 //165 return v1x > df2166 ? boost::math::ibetac(df2 / 2, df1 / 2, df2 / (df2 + v1x), Policy())167 : boost::math::ibeta(df1 / 2, df2 / 2, v1x / (df2 + v1x), Policy());168} // cdf169 170template <class RealType, class Policy>171BOOST_MATH_GPU_ENABLED inline RealType quantile(const fisher_f_distribution<RealType, Policy>& dist, const RealType& p)172{173 constexpr auto function = "boost::math::quantile(fisher_f_distribution<%1%> const&, %1%)";174 RealType df1 = dist.degrees_of_freedom1();175 RealType df2 = dist.degrees_of_freedom2();176 // Error check:177 RealType error_result = 0;178 if(false == (detail::check_df(179 function, df1, &error_result, Policy())180 && detail::check_df(181 function, df2, &error_result, Policy())182 && detail::check_probability(183 function, p, &error_result, Policy())))184 return error_result;185 186 // With optimizations turned on, gcc wrongly warns about y being used187 // uninitialized unless we initialize it to something:188 RealType x, y(0);189 190 x = boost::math::ibeta_inv(df1 / 2, df2 / 2, p, &y, Policy());191 192 return df2 * x / (df1 * y);193} // quantile194 195template <class RealType, class Policy>196BOOST_MATH_GPU_ENABLED inline RealType cdf(const complemented2_type<fisher_f_distribution<RealType, Policy>, RealType>& c)197{198 constexpr auto function = "boost::math::cdf(fisher_f_distribution<%1%> const&, %1%)";199 RealType df1 = c.dist.degrees_of_freedom1();200 RealType df2 = c.dist.degrees_of_freedom2();201 RealType x = c.param;202 // Error check:203 RealType error_result = 0;204 if(false == detail::check_df(205 function, df1, &error_result, Policy())206 && detail::check_df(207 function, df2, &error_result, Policy()))208 return error_result;209 210 if((x < 0) || !(boost::math::isfinite)(x))211 {212 return policies::raise_domain_error<RealType>(213 function, "Random Variable parameter was %1%, but must be > 0 !", x, Policy());214 }215 216 RealType v1x = df1 * x;217 //218 // There are two equivalent formulas used here, the aim is219 // to prevent the final argument to the incomplete beta220 // from being too close to 1: for some values of df1 and df2221 // the rate of change can be arbitrarily large in this area,222 // whilst the value we're passing will have lost information223 // content as a result of being 0.999999something. Better224 // to switch things around so we're passing 1-z instead.225 //226 return v1x > df2227 ? boost::math::ibeta(df2 / 2, df1 / 2, df2 / (df2 + v1x), Policy())228 : boost::math::ibetac(df1 / 2, df2 / 2, v1x / (df2 + v1x), Policy());229}230 231template <class RealType, class Policy>232BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type<fisher_f_distribution<RealType, Policy>, RealType>& c)233{234 constexpr auto function = "boost::math::quantile(fisher_f_distribution<%1%> const&, %1%)";235 RealType df1 = c.dist.degrees_of_freedom1();236 RealType df2 = c.dist.degrees_of_freedom2();237 RealType p = c.param;238 // Error check:239 RealType error_result = 0;240 if(false == (detail::check_df(241 function, df1, &error_result, Policy())242 && detail::check_df(243 function, df2, &error_result, Policy())244 && detail::check_probability(245 function, p, &error_result, Policy())))246 return error_result;247 248 RealType x, y;249 250 x = boost::math::ibetac_inv(df1 / 2, df2 / 2, p, &y, Policy());251 252 return df2 * x / (df1 * y);253}254 255template <class RealType, class Policy>256BOOST_MATH_GPU_ENABLED inline RealType mean(const fisher_f_distribution<RealType, Policy>& dist)257{ // Mean of F distribution = v.258 constexpr auto function = "boost::math::mean(fisher_f_distribution<%1%> const&)";259 RealType df1 = dist.degrees_of_freedom1();260 RealType df2 = dist.degrees_of_freedom2();261 // Error check:262 RealType error_result = 0;263 if(false == detail::check_df(264 function, df1, &error_result, Policy())265 && detail::check_df(266 function, df2, &error_result, Policy()))267 return error_result;268 if(df2 <= 2)269 {270 return policies::raise_domain_error<RealType>(271 function, "Second degree of freedom was %1% but must be > 2 in order for the distribution to have a mean.", df2, Policy());272 }273 return df2 / (df2 - 2);274} // mean275 276template <class RealType, class Policy>277BOOST_MATH_GPU_ENABLED inline RealType variance(const fisher_f_distribution<RealType, Policy>& dist)278{ // Variance of F distribution.279 constexpr auto function = "boost::math::variance(fisher_f_distribution<%1%> const&)";280 RealType df1 = dist.degrees_of_freedom1();281 RealType df2 = dist.degrees_of_freedom2();282 // Error check:283 RealType error_result = 0;284 if(false == detail::check_df(285 function, df1, &error_result, Policy())286 && detail::check_df(287 function, df2, &error_result, Policy()))288 return error_result;289 if(df2 <= 4)290 {291 return policies::raise_domain_error<RealType>(292 function, "Second degree of freedom was %1% but must be > 4 in order for the distribution to have a valid variance.", df2, Policy());293 }294 return 2 * df2 * df2 * (df1 + df2 - 2) / (df1 * (df2 - 2) * (df2 - 2) * (df2 - 4));295} // variance296 297template <class RealType, class Policy>298BOOST_MATH_GPU_ENABLED inline RealType mode(const fisher_f_distribution<RealType, Policy>& dist)299{300 constexpr auto function = "boost::math::mode(fisher_f_distribution<%1%> const&)";301 RealType df1 = dist.degrees_of_freedom1();302 RealType df2 = dist.degrees_of_freedom2();303 // Error check:304 RealType error_result = 0;305 if(false == detail::check_df(306 function, df1, &error_result, Policy())307 && detail::check_df(308 function, df2, &error_result, Policy()))309 return error_result;310 if(df1 <= 2)311 {312 return policies::raise_domain_error<RealType>(313 function, "First degree of freedom was %1% but must be > 2 in order for the distribution to have a mode.", df1, Policy());314 }315 return df2 * (df1 - 2) / (df1 * (df2 + 2));316}317 318//template <class RealType, class Policy>319//inline RealType median(const fisher_f_distribution<RealType, Policy>& dist)320//{ // Median of Fisher F distribution is not defined.321// return tools::domain_error<RealType>(BOOST_CURRENT_FUNCTION, "Median is not implemented, result is %1%!", boost::math::numeric_limits<RealType>::quiet_NaN());322// } // median323 324// Now implemented via quantile(half) in derived accessors.325 326template <class RealType, class Policy>327BOOST_MATH_GPU_ENABLED inline RealType skewness(const fisher_f_distribution<RealType, Policy>& dist)328{329 constexpr auto function = "boost::math::skewness(fisher_f_distribution<%1%> const&)";330 BOOST_MATH_STD_USING // ADL of std names331 // See http://mathworld.wolfram.com/F-Distribution.html332 RealType df1 = dist.degrees_of_freedom1();333 RealType df2 = dist.degrees_of_freedom2();334 // Error check:335 RealType error_result = 0;336 if(false == detail::check_df(337 function, df1, &error_result, Policy())338 && detail::check_df(339 function, df2, &error_result, Policy()))340 return error_result;341 if(df2 <= 6)342 {343 return policies::raise_domain_error<RealType>(344 function, "Second degree of freedom was %1% but must be > 6 in order for the distribution to have a skewness.", df2, Policy());345 }346 return 2 * (df2 + 2 * df1 - 2) * sqrt((2 * df2 - 8) / (df1 * (df2 + df1 - 2))) / (df2 - 6);347}348 349template <class RealType, class Policy>350BOOST_MATH_GPU_ENABLED RealType kurtosis_excess(const fisher_f_distribution<RealType, Policy>& dist);351 352template <class RealType, class Policy>353BOOST_MATH_GPU_ENABLED inline RealType kurtosis(const fisher_f_distribution<RealType, Policy>& dist)354{355 return 3 + kurtosis_excess(dist);356}357 358template <class RealType, class Policy>359BOOST_MATH_GPU_ENABLED inline RealType kurtosis_excess(const fisher_f_distribution<RealType, Policy>& dist)360{361 constexpr auto function = "boost::math::kurtosis_excess(fisher_f_distribution<%1%> const&)";362 // See http://mathworld.wolfram.com/F-Distribution.html363 RealType df1 = dist.degrees_of_freedom1();364 RealType df2 = dist.degrees_of_freedom2();365 // Error check:366 RealType error_result = 0;367 if(false == detail::check_df(368 function, df1, &error_result, Policy())369 && detail::check_df(370 function, df2, &error_result, Policy()))371 return error_result;372 if(df2 <= 8)373 {374 return policies::raise_domain_error<RealType>(375 function, "Second degree of freedom was %1% but must be > 8 in order for the distribution to have a kurtosis.", df2, Policy());376 }377 RealType df2_2 = df2 * df2;378 RealType df1_2 = df1 * df1;379 RealType n = -16 + 20 * df2 - 8 * df2_2 + df2_2 * df2 + 44 * df1 - 32 * df2 * df1 + 5 * df2_2 * df1 - 22 * df1_2 + 5 * df2 * df1_2;380 n *= 12;381 RealType d = df1 * (df2 - 6) * (df2 - 8) * (df1 + df2 - 2);382 return n / d;383}384 385} // namespace math386} // namespace boost387 388// This include must be at the end, *after* the accessors389// for this distribution have been defined, in order to390// keep compilers that support two-phase lookup happy.391#include <boost/math/distributions/detail/derived_accessors.hpp>392 393#endif // BOOST_MATH_DISTRIBUTIONS_FISHER_F_HPP394