math/distributions/fisher_f.hpp

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