math/distributions/chi_squared.hpp

0001 // Copyright John Maddock 2006, 2007.
0002 // Copyright Paul A. Bristow 2008, 2010.
0003
0004 // Use, modification and distribution are subject to the
0005 // Boost Software License, Version 1.0.
0006 // (See accompanying file LICENSE_1_0.txt
0007 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0008
0009 #ifndef BOOST_MATH_DISTRIBUTIONS_CHI_SQUARED_HPP
0010 #define BOOST_MATH_DISTRIBUTIONS_CHI_SQUARED_HPP
0011
0012 #include <boost/math/distributions/fwd.hpp>
0013 #include <boost/math/special_functions/gamma.hpp> // for incomplete beta.
0014 #include <boost/math/distributions/complement.hpp> // complements
0015 #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks
0016 #include <boost/math/special_functions/fpclassify.hpp>
0017
0018 #include <utility>
0019
0020 namespace boost{ namespace math{
0021
0022 template <class RealType = double, class Policy = policies::policy<> >
0023 class chi_squared_distribution
0024 {
0025 public:
0026    using value_type = RealType;
0027    using policy_type = Policy;
0028
0029    explicit chi_squared_distribution(RealType i) : m_df(i)
0030    {
0031       RealType result;
0032       detail::check_df(
0033          "boost::math::chi_squared_distribution<%1%>::chi_squared_distribution", m_df, &result, Policy());
0034    } // chi_squared_distribution
0035
0036    RealType degrees_of_freedom()const
0037    {
0038       return m_df;
0039    }
0040
0041    // Parameter estimation:
0042    static RealType find_degrees_of_freedom(
0043       RealType difference_from_variance,
0044       RealType alpha,
0045       RealType beta,
0046       RealType variance,
0047       RealType hint = 100);
0048
0049 private:
0050    //
0051    // Data member:
0052    //
0053    RealType m_df; // degrees of freedom is a positive real number.
0054 }; // class chi_squared_distribution
0055
0056 using chi_squared = chi_squared_distribution<double>;
0057
0058 #ifdef __cpp_deduction_guides
0059 template <class RealType>
0060 chi_squared_distribution(RealType)->chi_squared_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0061 #endif
0062
0063 #ifdef _MSC_VER
0064 #pragma warning(push)
0065 #pragma warning(disable:4127)
0066 #endif
0067
0068 template <class RealType, class Policy>
0069 inline std::pair<RealType, RealType> range(const chi_squared_distribution<RealType, Policy>& /*dist*/)
0070 { // Range of permissible values for random variable x.
0071   if (std::numeric_limits<RealType>::has_infinity)
0072   {
0073     return std::pair<RealType, RealType>(static_cast<RealType>(0), std::numeric_limits<RealType>::infinity()); // 0 to + infinity.
0074   }
0075   else
0076   {
0077     using boost::math::tools::max_value;
0078     return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // 0 to + max.
0079   }
0080 }
0081
0082 #ifdef _MSC_VER
0083 #pragma warning(pop)
0084 #endif
0085
0086 template <class RealType, class Policy>
0087 inline std::pair<RealType, RealType> support(const chi_squared_distribution<RealType, Policy>& /*dist*/)
0088 { // Range of supported values for random variable x.
0089    // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
0090    return std::pair<RealType, RealType>(static_cast<RealType>(0), tools::max_value<RealType>()); // 0 to + infinity.
0091 }
0092
0093 template <class RealType, class Policy>
0094 RealType pdf(const chi_squared_distribution<RealType, Policy>& dist, const RealType& chi_square)
0095 {
0096    BOOST_MATH_STD_USING  // for ADL of std functions
0097    RealType degrees_of_freedom = dist.degrees_of_freedom();
0098    // Error check:
0099    RealType error_result;
0100
0101    static const char* function = "boost::math::pdf(const chi_squared_distribution<%1%>&, %1%)";
0102
0103    if(false == detail::check_df(
0104          function, degrees_of_freedom, &error_result, Policy()))
0105       return error_result;
0106
0107    if((chi_square < 0) || !(boost::math::isfinite)(chi_square))
0108    {
0109       return policies::raise_domain_error<RealType>(
0110          function, "Chi Square parameter was %1%, but must be > 0 !", chi_square, Policy());
0111    }
0112
0113    if(chi_square == 0)
0114    {
0115       // Handle special cases:
0116       if(degrees_of_freedom < 2)
0117       {
0118          return policies::raise_overflow_error<RealType>(
0119             function, 0, Policy());
0120       }
0121       else if(degrees_of_freedom == 2)
0122       {
0123          return 0.5f;
0124       }
0125       else
0126       {
0127          return 0;
0128       }
0129    }
0130
0131    return gamma_p_derivative(degrees_of_freedom / 2, chi_square / 2, Policy()) / 2;
0132 } // pdf
0133
0134 template <class RealType, class Policy>
0135 inline RealType cdf(const chi_squared_distribution<RealType, Policy>& dist, const RealType& chi_square)
0136 {
0137    RealType degrees_of_freedom = dist.degrees_of_freedom();
0138    // Error check:
0139    RealType error_result;
0140    static const char* function = "boost::math::cdf(const chi_squared_distribution<%1%>&, %1%)";
0141
0142    if(false == detail::check_df(
0143          function, degrees_of_freedom, &error_result, Policy()))
0144       return error_result;
0145
0146    if((chi_square < 0) || !(boost::math::isfinite)(chi_square))
0147    {
0148       return policies::raise_domain_error<RealType>(
0149          function, "Chi Square parameter was %1%, but must be > 0 !", chi_square, Policy());
0150    }
0151
0152    return boost::math::gamma_p(degrees_of_freedom / 2, chi_square / 2, Policy());
0153 } // cdf
0154
0155 template <class RealType, class Policy>
0156 inline RealType quantile(const chi_squared_distribution<RealType, Policy>& dist, const RealType& p)
0157 {
0158    RealType degrees_of_freedom = dist.degrees_of_freedom();
0159    static const char* function = "boost::math::quantile(const chi_squared_distribution<%1%>&, %1%)";
0160    // Error check:
0161    RealType error_result;
0162    if(false ==
0163      (
0164        detail::check_df(function, degrees_of_freedom, &error_result, Policy())
0165        && detail::check_probability(function, p, &error_result, Policy()))
0166      )
0167      return error_result;
0168
0169    return 2 * boost::math::gamma_p_inv(degrees_of_freedom / 2, p, Policy());
0170 } // quantile
0171
0172 template <class RealType, class Policy>
0173 inline RealType cdf(const complemented2_type<chi_squared_distribution<RealType, Policy>, RealType>& c)
0174 {
0175    RealType const& degrees_of_freedom = c.dist.degrees_of_freedom();
0176    RealType const& chi_square = c.param;
0177    static const char* function = "boost::math::cdf(const chi_squared_distribution<%1%>&, %1%)";
0178    // Error check:
0179    RealType error_result;
0180    if(false == detail::check_df(
0181          function, degrees_of_freedom, &error_result, Policy()))
0182       return error_result;
0183
0184    if((chi_square < 0) || !(boost::math::isfinite)(chi_square))
0185    {
0186       return policies::raise_domain_error<RealType>(
0187          function, "Chi Square parameter was %1%, but must be > 0 !", chi_square, Policy());
0188    }
0189
0190    return boost::math::gamma_q(degrees_of_freedom / 2, chi_square / 2, Policy());
0191 }
0192
0193 template <class RealType, class Policy>
0194 inline RealType quantile(const complemented2_type<chi_squared_distribution<RealType, Policy>, RealType>& c)
0195 {
0196    RealType const& degrees_of_freedom = c.dist.degrees_of_freedom();
0197    RealType const& q = c.param;
0198    static const char* function = "boost::math::quantile(const chi_squared_distribution<%1%>&, %1%)";
0199    // Error check:
0200    RealType error_result;
0201    if(false == (
0202      detail::check_df(function, degrees_of_freedom, &error_result, Policy())
0203      && detail::check_probability(function, q, &error_result, Policy()))
0204      )
0205     return error_result;
0206
0207    return 2 * boost::math::gamma_q_inv(degrees_of_freedom / 2, q, Policy());
0208 }
0209
0210 template <class RealType, class Policy>
0211 inline RealType mean(const chi_squared_distribution<RealType, Policy>& dist)
0212 { // Mean of Chi-Squared distribution = v.
0213   return dist.degrees_of_freedom();
0214 } // mean
0215
0216 template <class RealType, class Policy>
0217 inline RealType variance(const chi_squared_distribution<RealType, Policy>& dist)
0218 { // Variance of Chi-Squared distribution = 2v.
0219   return 2 * dist.degrees_of_freedom();
0220 } // variance
0221
0222 template <class RealType, class Policy>
0223 inline RealType mode(const chi_squared_distribution<RealType, Policy>& dist)
0224 {
0225    RealType df = dist.degrees_of_freedom();
0226    static const char* function = "boost::math::mode(const chi_squared_distribution<%1%>&)";
0227
0228    if(df < 2)
0229       return policies::raise_domain_error<RealType>(
0230          function,
0231          "Chi-Squared distribution only has a mode for degrees of freedom >= 2, but got degrees of freedom = %1%.",
0232          df, Policy());
0233    return df - 2;
0234 }
0235
0236 template <class RealType, class Policy>
0237 inline RealType skewness(const chi_squared_distribution<RealType, Policy>& dist)
0238 {
0239    BOOST_MATH_STD_USING // For ADL
0240    RealType df = dist.degrees_of_freedom();
0241    return sqrt (8 / df);
0242 }
0243
0244 template <class RealType, class Policy>
0245 inline RealType kurtosis(const chi_squared_distribution<RealType, Policy>& dist)
0246 {
0247    RealType df = dist.degrees_of_freedom();
0248    return 3 + 12 / df;
0249 }
0250
0251 template <class RealType, class Policy>
0252 inline RealType kurtosis_excess(const chi_squared_distribution<RealType, Policy>& dist)
0253 {
0254    RealType df = dist.degrees_of_freedom();
0255    return 12 / df;
0256 }
0257
0258 //
0259 // Parameter estimation comes last:
0260 //
0261 namespace detail
0262 {
0263
0264 template <class RealType, class Policy>
0265 struct df_estimator
0266 {
0267    df_estimator(RealType a, RealType b, RealType variance, RealType delta)
0268       : alpha(a), beta(b), ratio(delta/variance)
0269    { // Constructor
0270    }
0271
0272    RealType operator()(const RealType& df)
0273    {
0274       if(df <= tools::min_value<RealType>())
0275          return 1;
0276       chi_squared_distribution<RealType, Policy> cs(df);
0277
0278       RealType result;
0279       if(ratio > 0)
0280       {
0281          RealType r = 1 + ratio;
0282          result = cdf(cs, quantile(complement(cs, alpha)) / r) - beta;
0283       }
0284       else
0285       { // ratio <= 0
0286          RealType r = 1 + ratio;
0287          result = cdf(complement(cs, quantile(cs, alpha) / r)) - beta;
0288       }
0289       return result;
0290    }
0291 private:
0292    RealType alpha;
0293    RealType beta;
0294    RealType ratio; // Difference from variance / variance, so fractional.
0295 };
0296
0297 } // namespace detail
0298
0299 template <class RealType, class Policy>
0300 RealType chi_squared_distribution<RealType, Policy>::find_degrees_of_freedom(
0301    RealType difference_from_variance,
0302    RealType alpha,
0303    RealType beta,
0304    RealType variance,
0305    RealType hint)
0306 {
0307    static const char* function = "boost::math::chi_squared_distribution<%1%>::find_degrees_of_freedom(%1%,%1%,%1%,%1%,%1%)";
0308    // Check for domain errors:
0309    RealType error_result;
0310    if(false ==
0311      detail::check_probability(function, alpha, &error_result, Policy())
0312      && detail::check_probability(function, beta, &error_result, Policy()))
0313    { // Either probability is outside 0 to 1.
0314       return error_result;
0315    }
0316
0317    if(hint <= 0)
0318    { // No hint given, so guess df = 1.
0319       hint = 1;
0320    }
0321
0322    detail::df_estimator<RealType, Policy> f(alpha, beta, variance, difference_from_variance);
0323    tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
0324    std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
0325    std::pair<RealType, RealType> r =
0326      tools::bracket_and_solve_root(f, hint, RealType(2), false, tol, max_iter, Policy());
0327    RealType result = r.first + (r.second - r.first) / 2;
0328    if(max_iter >= policies::get_max_root_iterations<Policy>())
0329    {
0330       policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:"
0331          " either there is no answer to how many degrees of freedom are required"
0332          " or the answer is infinite.  Current best guess is %1%", result, Policy());
0333    }
0334    return result;
0335 }
0336
0337 } // namespace math
0338 } // namespace boost
0339
0340 // This include must be at the end, *after* the accessors
0341 // for this distribution have been defined, in order to
0342 // keep compilers that support two-phase lookup happy.
0343 #include <boost/math/distributions/detail/derived_accessors.hpp>
0344
0345 #endif // BOOST_MATH_DISTRIBUTIONS_CHI_SQUARED_HPP