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