math/distributions/fisher_f.hpp

0001 // Copyright John Maddock 2006.
0002 // Copyright Matt Borland 2024.
0003 // Use, modification and distribution are subject to the
0004 // Boost Software License, Version 1.0.
0005 // (See accompanying file LICENSE_1_0.txt
0006 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0007
0008 #ifndef BOOST_MATH_DISTRIBUTIONS_FISHER_F_HPP
0009 #define BOOST_MATH_DISTRIBUTIONS_FISHER_F_HPP
0010
0011 #include <boost/math/tools/config.hpp>
0012 #include <boost/math/tools/tuple.hpp>
0013 #include <boost/math/tools/promotion.hpp>
0014 #include <boost/math/distributions/fwd.hpp>
0015 #include <boost/math/special_functions/beta.hpp> // for incomplete beta.
0016 #include <boost/math/distributions/complement.hpp> // complements
0017 #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks
0018 #include <boost/math/special_functions/fpclassify.hpp>
0019
0020 namespace boost{ namespace math{
0021
0022 template <class RealType = double, class Policy = policies::policy<> >
0023 class fisher_f_distribution
0024 {
0025 public:
0026    typedef RealType value_type;
0027    typedef Policy policy_type;
0028
0029    BOOST_MATH_GPU_ENABLED fisher_f_distribution(const RealType& i, const RealType& j) : m_df1(i), m_df2(j)
0030    {
0031       constexpr auto function = "fisher_f_distribution<%1%>::fisher_f_distribution";
0032       RealType result;
0033       detail::check_df(
0034          function, m_df1, &result, Policy());
0035       detail::check_df(
0036          function, m_df2, &result, Policy());
0037    } // fisher_f_distribution
0038
0039    BOOST_MATH_GPU_ENABLED RealType degrees_of_freedom1()const
0040    {
0041       return m_df1;
0042    }
0043    BOOST_MATH_GPU_ENABLED RealType degrees_of_freedom2()const
0044    {
0045       return m_df2;
0046    }
0047
0048 private:
0049    //
0050    // Data members:
0051    //
0052    RealType m_df1;  // degrees of freedom are a real number.
0053    RealType m_df2;  // degrees of freedom are a real number.
0054 };
0055
0056 typedef fisher_f_distribution<double> fisher_f;
0057
0058 #ifdef __cpp_deduction_guides
0059 template <class RealType>
0060 fisher_f_distribution(RealType,RealType)->fisher_f_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0061 #endif
0062
0063 template <class RealType, class Policy>
0064 BOOST_MATH_GPU_ENABLED inline const boost::math::pair<RealType, RealType> range(const fisher_f_distribution<RealType, Policy>& /*dist*/)
0065 { // Range of permissible values for random variable x.
0066    using boost::math::tools::max_value;
0067    return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());
0068 }
0069
0070 template <class RealType, class Policy>
0071 BOOST_MATH_GPU_ENABLED inline const boost::math::pair<RealType, RealType> support(const fisher_f_distribution<RealType, Policy>& /*dist*/)
0072 { // Range of supported values for random variable x.
0073    // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
0074    using boost::math::tools::max_value;
0075    return boost::math::pair<RealType, RealType>(static_cast<RealType>(0),  max_value<RealType>());
0076 }
0077
0078 template <class RealType, class Policy>
0079 BOOST_MATH_GPU_ENABLED RealType pdf(const fisher_f_distribution<RealType, Policy>& dist, const RealType& x)
0080 {
0081    BOOST_MATH_STD_USING  // for ADL of std functions
0082    RealType df1 = dist.degrees_of_freedom1();
0083    RealType df2 = dist.degrees_of_freedom2();
0084    // Error check:
0085    RealType error_result = 0;
0086    constexpr auto function = "boost::math::pdf(fisher_f_distribution<%1%> const&, %1%)";
0087    if(false == (detail::check_df(
0088          function, df1, &error_result, Policy())
0089          && detail::check_df(
0090          function, df2, &error_result, Policy())))
0091       return error_result;
0092
0093    if((x < 0) || !(boost::math::isfinite)(x))
0094    {
0095       return policies::raise_domain_error<RealType>(
0096          function, "Random variable parameter was %1%, but must be > 0 !", x, Policy());
0097    }
0098
0099    if(x == 0)
0100    {
0101       // special cases:
0102       if(df1 < 2)
0103          return policies::raise_overflow_error<RealType>(
0104             function, 0, Policy());
0105       else if(df1 == 2)
0106          return 1;
0107       else
0108          return 0;
0109    }
0110
0111    //
0112    // You reach this formula by direct differentiation of the
0113    // cdf expressed in terms of the incomplete beta.
0114    //
0115    // There are two versions so we don't pass a value of z
0116    // that is very close to 1 to ibeta_derivative: for some values
0117    // of df1 and df2, all the change takes place in this area.
0118    //
0119    RealType v1x = df1 * x;
0120    RealType result;
0121    if(v1x > df2)
0122    {
0123       result = (df2 * df1) / ((df2 + v1x) * (df2 + v1x));
0124       result *= ibeta_derivative(df2 / 2, df1 / 2, df2 / (df2 + v1x), Policy());
0125    }
0126    else
0127    {
0128       result = df2 + df1 * x;
0129       result = (result * df1 - x * df1 * df1) / (result * result);
0130       result *= ibeta_derivative(df1 / 2, df2 / 2, v1x / (df2 + v1x), Policy());
0131    }
0132    return result;
0133 } // pdf
0134
0135 template <class RealType, class Policy>
0136 BOOST_MATH_GPU_ENABLED inline RealType cdf(const fisher_f_distribution<RealType, Policy>& dist, const RealType& x)
0137 {
0138    constexpr auto function = "boost::math::cdf(fisher_f_distribution<%1%> const&, %1%)";
0139    RealType df1 = dist.degrees_of_freedom1();
0140    RealType df2 = dist.degrees_of_freedom2();
0141    // Error check:
0142    RealType error_result = 0;
0143    if(false == detail::check_df(
0144          function, df1, &error_result, Policy())
0145          && detail::check_df(
0146          function, df2, &error_result, Policy()))
0147       return error_result;
0148
0149    if((x < 0) || !(boost::math::isfinite)(x))
0150    {
0151       return policies::raise_domain_error<RealType>(
0152          function, "Random Variable parameter was %1%, but must be > 0 !", x, Policy());
0153    }
0154
0155    RealType v1x = df1 * x;
0156    //
0157    // There are two equivalent formulas used here, the aim is
0158    // to prevent the final argument to the incomplete beta
0159    // from being too close to 1: for some values of df1 and df2
0160    // the rate of change can be arbitrarily large in this area,
0161    // whilst the value we're passing will have lost information
0162    // content as a result of being 0.999999something.  Better
0163    // to switch things around so we're passing 1-z instead.
0164    //
0165    return v1x > df2
0166       ? boost::math::ibetac(df2 / 2, df1 / 2, df2 / (df2 + v1x), Policy())
0167       : boost::math::ibeta(df1 / 2, df2 / 2, v1x / (df2 + v1x), Policy());
0168 } // cdf
0169
0170 template <class RealType, class Policy>
0171 BOOST_MATH_GPU_ENABLED inline RealType quantile(const fisher_f_distribution<RealType, Policy>& dist, const RealType& p)
0172 {
0173    constexpr auto function = "boost::math::quantile(fisher_f_distribution<%1%> const&, %1%)";
0174    RealType df1 = dist.degrees_of_freedom1();
0175    RealType df2 = dist.degrees_of_freedom2();
0176    // Error check:
0177    RealType error_result = 0;
0178    if(false == (detail::check_df(
0179             function, df1, &error_result, Policy())
0180          && detail::check_df(
0181             function, df2, &error_result, Policy())
0182          && detail::check_probability(
0183             function, p, &error_result, Policy())))
0184       return error_result;
0185
0186    // With optimizations turned on, gcc wrongly warns about y being used
0187    // uninitialized unless we initialize it to something:
0188    RealType x, y(0);
0189
0190    x = boost::math::ibeta_inv(df1 / 2, df2 / 2, p, &y, Policy());
0191
0192    return df2 * x / (df1 * y);
0193 } // quantile
0194
0195 template <class RealType, class Policy>
0196 BOOST_MATH_GPU_ENABLED inline RealType cdf(const complemented2_type<fisher_f_distribution<RealType, Policy>, RealType>& c)
0197 {
0198    constexpr auto function = "boost::math::cdf(fisher_f_distribution<%1%> const&, %1%)";
0199    RealType df1 = c.dist.degrees_of_freedom1();
0200    RealType df2 = c.dist.degrees_of_freedom2();
0201    RealType x = c.param;
0202    // Error check:
0203    RealType error_result = 0;
0204    if(false == detail::check_df(
0205          function, df1, &error_result, Policy())
0206          && detail::check_df(
0207          function, df2, &error_result, Policy()))
0208       return error_result;
0209
0210    if((x < 0) || !(boost::math::isfinite)(x))
0211    {
0212       return policies::raise_domain_error<RealType>(
0213          function, "Random Variable parameter was %1%, but must be > 0 !", x, Policy());
0214    }
0215
0216    RealType v1x = df1 * x;
0217    //
0218    // There are two equivalent formulas used here, the aim is
0219    // to prevent the final argument to the incomplete beta
0220    // from being too close to 1: for some values of df1 and df2
0221    // the rate of change can be arbitrarily large in this area,
0222    // whilst the value we're passing will have lost information
0223    // content as a result of being 0.999999something.  Better
0224    // to switch things around so we're passing 1-z instead.
0225    //
0226    return v1x > df2
0227       ? boost::math::ibeta(df2 / 2, df1 / 2, df2 / (df2 + v1x), Policy())
0228       : boost::math::ibetac(df1 / 2, df2 / 2, v1x / (df2 + v1x), Policy());
0229 }
0230
0231 template <class RealType, class Policy>
0232 BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type<fisher_f_distribution<RealType, Policy>, RealType>& c)
0233 {
0234    constexpr auto function = "boost::math::quantile(fisher_f_distribution<%1%> const&, %1%)";
0235    RealType df1 = c.dist.degrees_of_freedom1();
0236    RealType df2 = c.dist.degrees_of_freedom2();
0237    RealType p = c.param;
0238    // Error check:
0239    RealType error_result = 0;
0240    if(false == (detail::check_df(
0241             function, df1, &error_result, Policy())
0242          && detail::check_df(
0243             function, df2, &error_result, Policy())
0244          && detail::check_probability(
0245             function, p, &error_result, Policy())))
0246       return error_result;
0247
0248    RealType x, y;
0249
0250    x = boost::math::ibetac_inv(df1 / 2, df2 / 2, p, &y, Policy());
0251
0252    return df2 * x / (df1 * y);
0253 }
0254
0255 template <class RealType, class Policy>
0256 BOOST_MATH_GPU_ENABLED inline RealType mean(const fisher_f_distribution<RealType, Policy>& dist)
0257 { // Mean of F distribution = v.
0258    constexpr auto function = "boost::math::mean(fisher_f_distribution<%1%> const&)";
0259    RealType df1 = dist.degrees_of_freedom1();
0260    RealType df2 = dist.degrees_of_freedom2();
0261    // Error check:
0262    RealType error_result = 0;
0263    if(false == detail::check_df(
0264             function, df1, &error_result, Policy())
0265          && detail::check_df(
0266             function, df2, &error_result, Policy()))
0267       return error_result;
0268    if(df2 <= 2)
0269    {
0270       return policies::raise_domain_error<RealType>(
0271          function, "Second degree of freedom was %1% but must be > 2 in order for the distribution to have a mean.", df2, Policy());
0272    }
0273    return df2 / (df2 - 2);
0274 } // mean
0275
0276 template <class RealType, class Policy>
0277 BOOST_MATH_GPU_ENABLED inline RealType variance(const fisher_f_distribution<RealType, Policy>& dist)
0278 { // Variance of F distribution.
0279    constexpr auto function = "boost::math::variance(fisher_f_distribution<%1%> const&)";
0280    RealType df1 = dist.degrees_of_freedom1();
0281    RealType df2 = dist.degrees_of_freedom2();
0282    // Error check:
0283    RealType error_result = 0;
0284    if(false == detail::check_df(
0285             function, df1, &error_result, Policy())
0286          && detail::check_df(
0287             function, df2, &error_result, Policy()))
0288       return error_result;
0289    if(df2 <= 4)
0290    {
0291       return policies::raise_domain_error<RealType>(
0292          function, "Second degree of freedom was %1% but must be > 4 in order for the distribution to have a valid variance.", df2, Policy());
0293    }
0294    return 2 * df2 * df2 * (df1 + df2 - 2) / (df1 * (df2 - 2) * (df2 - 2) * (df2 - 4));
0295 } // variance
0296
0297 template <class RealType, class Policy>
0298 BOOST_MATH_GPU_ENABLED inline RealType mode(const fisher_f_distribution<RealType, Policy>& dist)
0299 {
0300    constexpr auto function = "boost::math::mode(fisher_f_distribution<%1%> const&)";
0301    RealType df1 = dist.degrees_of_freedom1();
0302    RealType df2 = dist.degrees_of_freedom2();
0303    // Error check:
0304    RealType error_result = 0;
0305    if(false == detail::check_df(
0306             function, df1, &error_result, Policy())
0307          && detail::check_df(
0308             function, df2, &error_result, Policy()))
0309       return error_result;
0310    if(df1 <= 2)
0311    {
0312       return policies::raise_domain_error<RealType>(
0313          function, "First degree of freedom was %1% but must be > 2 in order for the distribution to have a mode.", df1, Policy());
0314    }
0315    return df2 * (df1 - 2) / (df1 * (df2 + 2));
0316 }
0317
0318 //template <class RealType, class Policy>
0319 //inline RealType median(const fisher_f_distribution<RealType, Policy>& dist)
0320 //{ // Median of Fisher F distribution is not defined.
0321 //  return tools::domain_error<RealType>(BOOST_CURRENT_FUNCTION, "Median is not implemented, result is %1%!", boost::math::numeric_limits<RealType>::quiet_NaN());
0322 //  } // median
0323
0324 // Now implemented via quantile(half) in derived accessors.
0325
0326 template <class RealType, class Policy>
0327 BOOST_MATH_GPU_ENABLED inline RealType skewness(const fisher_f_distribution<RealType, Policy>& dist)
0328 {
0329    constexpr auto function = "boost::math::skewness(fisher_f_distribution<%1%> const&)";
0330    BOOST_MATH_STD_USING // ADL of std names
0331    // See http://mathworld.wolfram.com/F-Distribution.html
0332    RealType df1 = dist.degrees_of_freedom1();
0333    RealType df2 = dist.degrees_of_freedom2();
0334    // Error check:
0335    RealType error_result = 0;
0336    if(false == detail::check_df(
0337             function, df1, &error_result, Policy())
0338          && detail::check_df(
0339             function, df2, &error_result, Policy()))
0340       return error_result;
0341    if(df2 <= 6)
0342    {
0343       return policies::raise_domain_error<RealType>(
0344          function, "Second degree of freedom was %1% but must be > 6 in order for the distribution to have a skewness.", df2, Policy());
0345    }
0346    return 2 * (df2 + 2 * df1 - 2) * sqrt((2 * df2 - 8) / (df1 * (df2 + df1 - 2))) / (df2 - 6);
0347 }
0348
0349 template <class RealType, class Policy>
0350 BOOST_MATH_GPU_ENABLED RealType kurtosis_excess(const fisher_f_distribution<RealType, Policy>& dist);
0351
0352 template <class RealType, class Policy>
0353 BOOST_MATH_GPU_ENABLED inline RealType kurtosis(const fisher_f_distribution<RealType, Policy>& dist)
0354 {
0355    return 3 + kurtosis_excess(dist);
0356 }
0357
0358 template <class RealType, class Policy>
0359 BOOST_MATH_GPU_ENABLED inline RealType kurtosis_excess(const fisher_f_distribution<RealType, Policy>& dist)
0360 {
0361    constexpr auto function = "boost::math::kurtosis_excess(fisher_f_distribution<%1%> const&)";
0362    // See http://mathworld.wolfram.com/F-Distribution.html
0363    RealType df1 = dist.degrees_of_freedom1();
0364    RealType df2 = dist.degrees_of_freedom2();
0365    // Error check:
0366    RealType error_result = 0;
0367    if(false == detail::check_df(
0368             function, df1, &error_result, Policy())
0369          && detail::check_df(
0370             function, df2, &error_result, Policy()))
0371       return error_result;
0372    if(df2 <= 8)
0373    {
0374       return policies::raise_domain_error<RealType>(
0375          function, "Second degree of freedom was %1% but must be > 8 in order for the distribution to have a kurtosis.", df2, Policy());
0376    }
0377    RealType df2_2 = df2 * df2;
0378    RealType df1_2 = df1 * df1;
0379    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;
0380    n *= 12;
0381    RealType d = df1 * (df2 - 6) * (df2 - 8) * (df1 + df2 - 2);
0382    return n / d;
0383 }
0384
0385 } // namespace math
0386 } // namespace boost
0387
0388 // This include must be at the end, *after* the accessors
0389 // for this distribution have been defined, in order to
0390 // keep compilers that support two-phase lookup happy.
0391 #include <boost/math/distributions/detail/derived_accessors.hpp>
0392
0393 #endif // BOOST_MATH_DISTRIBUTIONS_FISHER_F_HPP