math/distributions/students_t.hpp

0001 //  Copyright John Maddock 2006.
0002 //  Copyright Paul A. Bristow 2006, 2012, 2017.
0003 //  Copyright Thomas Mang 2012.
0004 //  Copyright Matt Borland 2024.
0005 //  Use, modification and distribution are subject to the
0006 //  Boost Software License, Version 1.0. (See accompanying file
0007 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0008
0009 #ifndef BOOST_STATS_STUDENTS_T_HPP
0010 #define BOOST_STATS_STUDENTS_T_HPP
0011
0012 // http://en.wikipedia.org/wiki/Student%27s_t_distribution
0013 // http://www.itl.nist.gov/div898/handbook/eda/section3/eda3664.htm
0014
0015 #include <boost/math/tools/config.hpp>
0016 #include <boost/math/tools/tuple.hpp>
0017 #include <boost/math/tools/cstdint.hpp>
0018 #include <boost/math/tools/numeric_limits.hpp>
0019 #include <boost/math/distributions/fwd.hpp>
0020 #include <boost/math/special_functions/beta.hpp> // for ibeta(a, b, x).
0021 #include <boost/math/special_functions/digamma.hpp>
0022 #include <boost/math/distributions/complement.hpp>
0023 #include <boost/math/distributions/detail/common_error_handling.hpp>
0024 #include <boost/math/distributions/normal.hpp>
0025 #include <boost/math/policies/policy.hpp>
0026
0027 #ifdef _MSC_VER
0028 # pragma warning(push)
0029 # pragma warning(disable: 4702) // unreachable code (return after domain_error throw).
0030 #endif
0031
0032 namespace boost { namespace math {
0033
0034 template <class RealType = double, class Policy = policies::policy<> >
0035 class students_t_distribution
0036 {
0037 public:
0038    typedef RealType value_type;
0039    typedef Policy policy_type;
0040
0041    BOOST_MATH_GPU_ENABLED students_t_distribution(RealType df) : df_(df)
0042    { // Constructor.
0043       RealType result;
0044       detail::check_df_gt0_to_inf( // Checks that df > 0 or df == inf.
0045          "boost::math::students_t_distribution<%1%>::students_t_distribution", df_, &result, Policy());
0046    } // students_t_distribution
0047
0048    BOOST_MATH_GPU_ENABLED RealType degrees_of_freedom()const
0049    {
0050       return df_;
0051    }
0052
0053    // Parameter estimation:
0054    BOOST_MATH_GPU_ENABLED static RealType find_degrees_of_freedom(
0055       RealType difference_from_mean,
0056       RealType alpha,
0057       RealType beta,
0058       RealType sd,
0059       RealType hint = 100);
0060
0061 private:
0062    // Data member:
0063    RealType df_;  // degrees of freedom is a real number > 0 or +infinity.
0064 };
0065
0066 typedef students_t_distribution<double> students_t; // Convenience typedef for double version.
0067
0068 #ifdef __cpp_deduction_guides
0069 template <class RealType>
0070 students_t_distribution(RealType)->students_t_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0071 #endif
0072
0073 template <class RealType, class Policy>
0074 BOOST_MATH_GPU_ENABLED inline const boost::math::pair<RealType, RealType> range(const students_t_distribution<RealType, Policy>& /*dist*/)
0075 { // Range of permissible values for random variable x.
0076   // Now including infinity.
0077    using boost::math::tools::max_value;
0078    //return boost::math::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>());
0079    return boost::math::pair<RealType, RealType>(((::boost::math::numeric_limits<RealType>::is_specialized & ::boost::math::numeric_limits<RealType>::has_infinity) ? -boost::math::numeric_limits<RealType>::infinity() : -max_value<RealType>()), ((::boost::math::numeric_limits<RealType>::is_specialized & ::boost::math::numeric_limits<RealType>::has_infinity) ? +boost::math::numeric_limits<RealType>::infinity() : +max_value<RealType>()));
0080 }
0081
0082 template <class RealType, class Policy>
0083 BOOST_MATH_GPU_ENABLED inline const boost::math::pair<RealType, RealType> support(const students_t_distribution<RealType, Policy>& /*dist*/)
0084 { // Range of supported values for random variable x.
0085   // Now including infinity.
0086    // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
0087    using boost::math::tools::max_value;
0088    //return boost::math::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>());
0089    return boost::math::pair<RealType, RealType>(((::boost::math::numeric_limits<RealType>::is_specialized & ::boost::math::numeric_limits<RealType>::has_infinity) ? -boost::math::numeric_limits<RealType>::infinity() : -max_value<RealType>()), ((::boost::math::numeric_limits<RealType>::is_specialized & ::boost::math::numeric_limits<RealType>::has_infinity) ? +boost::math::numeric_limits<RealType>::infinity() : +max_value<RealType>()));
0090 }
0091
0092 template <class RealType, class Policy>
0093 BOOST_MATH_GPU_ENABLED inline RealType pdf(const students_t_distribution<RealType, Policy>& dist, const RealType& x)
0094 {
0095    BOOST_FPU_EXCEPTION_GUARD
0096    BOOST_MATH_STD_USING  // for ADL of std functions.
0097
0098    RealType error_result;
0099    if(false == detail::check_x_not_NaN(
0100       "boost::math::pdf(const students_t_distribution<%1%>&, %1%)", x, &error_result, Policy()))
0101       return error_result;
0102    RealType df = dist.degrees_of_freedom();
0103    if(false == detail::check_df_gt0_to_inf( // Check that df > 0 or == +infinity.
0104       "boost::math::pdf(const students_t_distribution<%1%>&, %1%)", df, &error_result, Policy()))
0105       return error_result;
0106
0107    RealType result;
0108    if ((boost::math::isinf)(x))
0109    { // - or +infinity.
0110      result = static_cast<RealType>(0);
0111      return result;
0112    }
0113    RealType limit = policies::get_epsilon<RealType, Policy>();
0114    // Use policies so that if policy requests lower precision,
0115    // then get the normal distribution approximation earlier.
0116    limit = static_cast<RealType>(1) / limit; // 1/eps
0117    // for 64-bit double 1/eps = 4503599627370496
0118    if (df > limit)
0119    { // Special case for really big degrees_of_freedom > 1 / eps
0120      // - use normal distribution which is much faster and more accurate.
0121      normal_distribution<RealType, Policy> n(0, 1);
0122      result = pdf(n, x);
0123    }
0124    else
0125    { //
0126      RealType basem1 = x * x / df;
0127      if(basem1 < 0.125)
0128      {
0129         result = exp(-boost::math::log1p(basem1, Policy()) * (1+df) / 2);
0130      }
0131      else
0132      {
0133         result = pow(1 / (1 + basem1), (df + 1) / 2);
0134      }
0135      result /= sqrt(df) * boost::math::beta(df / 2, RealType(0.5f), Policy());
0136    }
0137    return result;
0138 } // pdf
0139
0140 template <class RealType, class Policy>
0141 BOOST_MATH_GPU_ENABLED inline RealType cdf(const students_t_distribution<RealType, Policy>& dist, const RealType& x)
0142 {
0143    RealType error_result;
0144    // degrees_of_freedom > 0 or infinity check:
0145    RealType df = dist.degrees_of_freedom();
0146    if (false == detail::check_df_gt0_to_inf(  // Check that df > 0 or == +infinity.
0147      "boost::math::cdf(const students_t_distribution<%1%>&, %1%)", df, &error_result, Policy()))
0148    {
0149      return error_result;
0150    }
0151    // Check for bad x first.
0152    if(false == detail::check_x_not_NaN(
0153       "boost::math::cdf(const students_t_distribution<%1%>&, %1%)", x, &error_result, Policy()))
0154    {
0155       return error_result;
0156    }
0157    if (x == 0)
0158    { // Special case with exact result.
0159      return static_cast<RealType>(0.5);
0160    }
0161    if ((boost::math::isinf)(x))
0162    { // x == - or + infinity, regardless of df.
0163      return ((x < 0) ? static_cast<RealType>(0) : static_cast<RealType>(1));
0164    }
0165
0166    RealType limit = policies::get_epsilon<RealType, Policy>();
0167    // Use policies so that if policy requests lower precision,
0168    // then get the normal distribution approximation earlier.
0169    limit = static_cast<RealType>(1) / limit; // 1/eps
0170    // for 64-bit double 1/eps = 4503599627370496
0171    if (df > limit)
0172    { // Special case for really big degrees_of_freedom > 1 / eps (perhaps infinite?)
0173      // - use normal distribution which is much faster and more accurate.
0174      normal_distribution<RealType, Policy> n(0, 1);
0175      RealType result = cdf(n, x);
0176      return result;
0177    }
0178    else
0179    { // normal df case.
0180      //
0181      // Calculate probability of Student's t using the incomplete beta function.
0182      // probability = ibeta(degrees_of_freedom / 2, 1/2, degrees_of_freedom / (degrees_of_freedom + t*t))
0183      //
0184      // However when t is small compared to the degrees of freedom, that formula
0185      // suffers from rounding error, use the identity formula to work around
0186      // the problem:
0187      //
0188      // I[x](a,b) = 1 - I[1-x](b,a)
0189      //
0190      // and:
0191      //
0192      //     x = df / (df + t^2)
0193      //
0194      // so:
0195      //
0196      // 1 - x = t^2 / (df + t^2)
0197      //
0198      RealType x2 = x * x;
0199      RealType probability;
0200      if(df > 2 * x2)
0201      {
0202         RealType z = x2 / (df + x2);
0203         probability = ibetac(static_cast<RealType>(0.5), df / 2, z, Policy()) / 2;
0204      }
0205      else
0206      {
0207         RealType z = df / (df + x2);
0208         probability = ibeta(df / 2, static_cast<RealType>(0.5), z, Policy()) / 2;
0209      }
0210      return (x > 0 ? 1   - probability : probability);
0211   }
0212 } // cdf
0213
0214 template <class RealType, class Policy>
0215 BOOST_MATH_GPU_ENABLED inline RealType quantile(const students_t_distribution<RealType, Policy>& dist, const RealType& p)
0216 {
0217    BOOST_MATH_STD_USING // for ADL of std functions
0218    //
0219    // Obtain parameters:
0220    RealType probability = p;
0221
0222    // Check for domain errors:
0223    RealType df = dist.degrees_of_freedom();
0224    constexpr auto function = "boost::math::quantile(const students_t_distribution<%1%>&, %1%)";
0225    RealType error_result;
0226    if(false == (detail::check_df_gt0_to_inf( // Check that df > 0 or == +infinity.
0227       function, df, &error_result, Policy())
0228          && detail::check_probability(function, probability, &error_result, Policy())))
0229       return error_result;
0230    // Special cases, regardless of degrees_of_freedom.
0231    if (probability == 0)
0232       return -policies::raise_overflow_error<RealType>(function, 0, Policy());
0233    if (probability == 1)
0234      return policies::raise_overflow_error<RealType>(function, 0, Policy());
0235    if (probability == static_cast<RealType>(0.5))
0236      return 0;  //
0237    //
0238 #if 0
0239    // This next block is disabled in favour of a faster method than
0240    // incomplete beta inverse, but code retained for future reference:
0241    //
0242    // Calculate quantile of Student's t using the incomplete beta function inverse:
0243    probability = (probability > 0.5) ? 1 - probability : probability;
0244    RealType t, x, y;
0245    x = ibeta_inv(degrees_of_freedom / 2, RealType(0.5), 2 * probability, &y);
0246    if(degrees_of_freedom * y > tools::max_value<RealType>() * x)
0247       t = tools::overflow_error<RealType>(function);
0248    else
0249       t = sqrt(degrees_of_freedom * y / x);
0250    //
0251    // Figure out sign based on the size of p:
0252    //
0253    if(p < 0.5)
0254       t = -t;
0255
0256    return t;
0257 #endif
0258    //
0259    // Depending on how many digits RealType has, this may forward
0260    // to the incomplete beta inverse as above.  Otherwise uses a
0261    // faster method that is accurate to ~15 digits everywhere
0262    // and a couple of epsilon at double precision and in the central
0263    // region where most use cases will occur...
0264    //
0265    return boost::math::detail::fast_students_t_quantile(df, probability, Policy());
0266 } // quantile
0267
0268 template <class RealType, class Policy>
0269 BOOST_MATH_GPU_ENABLED inline RealType cdf(const complemented2_type<students_t_distribution<RealType, Policy>, RealType>& c)
0270 {
0271    return cdf(c.dist, -c.param);
0272 }
0273
0274 template <class RealType, class Policy>
0275 BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type<students_t_distribution<RealType, Policy>, RealType>& c)
0276 {
0277    return -quantile(c.dist, c.param);
0278 }
0279
0280 //
0281 // Parameter estimation follows:
0282 //
0283 namespace detail{
0284 //
0285 // Functors for finding degrees of freedom:
0286 //
0287 template <class RealType, class Policy>
0288 struct sample_size_func
0289 {
0290    BOOST_MATH_GPU_ENABLED sample_size_func(RealType a, RealType b, RealType s, RealType d)
0291       : alpha(a), beta(b), ratio(s*s/(d*d)) {}
0292
0293    BOOST_MATH_GPU_ENABLED RealType operator()(const RealType& df)
0294    {
0295       if(df <= tools::min_value<RealType>())
0296       { //
0297          return 1;
0298       }
0299       students_t_distribution<RealType, Policy> t(df);
0300       RealType qa = quantile(complement(t, alpha));
0301       RealType qb = quantile(complement(t, beta));
0302       qa += qb;
0303       qa *= qa;
0304       qa *= ratio;
0305       qa -= (df + 1);
0306       return qa;
0307    }
0308    RealType alpha, beta, ratio;
0309 };
0310
0311 }  // namespace detail
0312
0313 template <class RealType, class Policy>
0314 BOOST_MATH_GPU_ENABLED RealType students_t_distribution<RealType, Policy>::find_degrees_of_freedom(
0315       RealType difference_from_mean,
0316       RealType alpha,
0317       RealType beta,
0318       RealType sd,
0319       RealType hint)
0320 {
0321    constexpr auto function = "boost::math::students_t_distribution<%1%>::find_degrees_of_freedom";
0322    //
0323    // Check for domain errors:
0324    //
0325    RealType error_result;
0326    if(false == detail::check_probability(
0327       function, alpha, &error_result, Policy())
0328          && detail::check_probability(function, beta, &error_result, Policy()))
0329       return error_result;
0330
0331    if(hint <= 0)
0332       hint = 1;
0333
0334    detail::sample_size_func<RealType, Policy> f(alpha, beta, sd, difference_from_mean);
0335    tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
0336    boost::math::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
0337    boost::math::pair<RealType, RealType> r = tools::bracket_and_solve_root(f, hint, RealType(2), false, tol, max_iter, Policy());
0338    RealType result = r.first + (r.second - r.first) / 2;
0339    if(max_iter >= policies::get_max_root_iterations<Policy>())
0340    {
0341       return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time: either there is no answer to how many degrees of freedom are required" // LCOV_EXCL_LINE
0342          " or the answer is infinite.  Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE
0343    }
0344    return result;
0345 }
0346
0347 template <class RealType, class Policy>
0348 BOOST_MATH_GPU_ENABLED inline RealType mode(const students_t_distribution<RealType, Policy>& /*dist*/)
0349 {
0350   // Assume no checks on degrees of freedom are useful (unlike mean).
0351    return 0; // Always zero by definition.
0352 }
0353
0354 template <class RealType, class Policy>
0355 BOOST_MATH_GPU_ENABLED inline RealType median(const students_t_distribution<RealType, Policy>& /*dist*/)
0356 {
0357    // Assume no checks on degrees of freedom are useful (unlike mean).
0358    return 0; // Always zero by definition.
0359 }
0360
0361 // See section 5.1 on moments at  http://en.wikipedia.org/wiki/Student%27s_t-distribution
0362
0363 template <class RealType, class Policy>
0364 BOOST_MATH_GPU_ENABLED inline RealType mean(const students_t_distribution<RealType, Policy>& dist)
0365 {  // Revised for https://svn.boost.org/trac/boost/ticket/7177
0366    RealType df = dist.degrees_of_freedom();
0367    if(((boost::math::isnan)(df)) || (df <= 1) )
0368    { // mean is undefined for moment <= 1!
0369       return policies::raise_domain_error<RealType>(
0370       "boost::math::mean(students_t_distribution<%1%> const&, %1%)",
0371       "Mean is undefined for degrees of freedom < 1 but got %1%.", df, Policy());
0372       return boost::math::numeric_limits<RealType>::quiet_NaN();
0373    }
0374    return 0;
0375 } // mean
0376
0377 template <class RealType, class Policy>
0378 BOOST_MATH_GPU_ENABLED inline RealType variance(const students_t_distribution<RealType, Policy>& dist)
0379 { // http://en.wikipedia.org/wiki/Student%27s_t-distribution
0380   // Revised for https://svn.boost.org/trac/boost/ticket/7177
0381   RealType df = dist.degrees_of_freedom();
0382   if ((boost::math::isnan)(df) || (df <= 2))
0383   { // NaN or undefined for <= 2.
0384      return policies::raise_domain_error<RealType>(
0385       "boost::math::variance(students_t_distribution<%1%> const&, %1%)",
0386       "variance is undefined for degrees of freedom <= 2, but got %1%.",
0387       df, Policy());
0388     return boost::math::numeric_limits<RealType>::quiet_NaN(); // Undefined.
0389   }
0390   if ((boost::math::isinf)(df))
0391   { // +infinity.
0392     return 1;
0393   }
0394   RealType limit = policies::get_epsilon<RealType, Policy>();
0395   // Use policies so that if policy requests lower precision,
0396   // then get the normal distribution approximation earlier.
0397   limit = static_cast<RealType>(1) / limit; // 1/eps
0398   // for 64-bit double 1/eps = 4503599627370496
0399   if (df > limit)
0400   { // Special case for really big degrees_of_freedom > 1 / eps.
0401     return 1;
0402   }
0403   else
0404   {
0405     return df / (df - 2);
0406   }
0407 } // variance
0408
0409 template <class RealType, class Policy>
0410 BOOST_MATH_GPU_ENABLED inline RealType skewness(const students_t_distribution<RealType, Policy>& dist)
0411 {
0412     RealType df = dist.degrees_of_freedom();
0413    if( ((boost::math::isnan)(df)) || (dist.degrees_of_freedom() <= 3))
0414    { // Undefined for moment k = 3.
0415       return policies::raise_domain_error<RealType>(
0416          "boost::math::skewness(students_t_distribution<%1%> const&, %1%)",
0417          "Skewness is undefined for degrees of freedom <= 3, but got %1%.",
0418          dist.degrees_of_freedom(), Policy());
0419       return boost::math::numeric_limits<RealType>::quiet_NaN();
0420    }
0421    return 0; // For all valid df, including infinity.
0422 } // skewness
0423
0424 template <class RealType, class Policy>
0425 BOOST_MATH_GPU_ENABLED inline RealType kurtosis(const students_t_distribution<RealType, Policy>& dist)
0426 {
0427    RealType df = dist.degrees_of_freedom();
0428    if(((boost::math::isnan)(df)) || (df <= 4))
0429    { // Undefined or infinity for moment k = 4.
0430       return policies::raise_domain_error<RealType>(
0431        "boost::math::kurtosis(students_t_distribution<%1%> const&, %1%)",
0432        "Kurtosis is undefined for degrees of freedom <= 4, but got %1%.",
0433         df, Policy());
0434         return boost::math::numeric_limits<RealType>::quiet_NaN(); // Undefined.
0435    }
0436    if ((boost::math::isinf)(df))
0437    { // +infinity.
0438      return 3;
0439    }
0440    RealType limit = policies::get_epsilon<RealType, Policy>();
0441    // Use policies so that if policy requests lower precision,
0442    // then get the normal distribution approximation earlier.
0443    limit = static_cast<RealType>(1) / limit; // 1/eps
0444    // for 64-bit double 1/eps = 4503599627370496
0445    if (df > limit)
0446    { // Special case for really big degrees_of_freedom > 1 / eps.
0447      return 3;
0448    }
0449    else
0450    {
0451      //return 3 * (df - 2) / (df - 4); re-arranged to
0452      return 6 / (df - 4) + 3;
0453    }
0454 } // kurtosis
0455
0456 template <class RealType, class Policy>
0457 BOOST_MATH_GPU_ENABLED inline RealType kurtosis_excess(const students_t_distribution<RealType, Policy>& dist)
0458 {
0459    // see http://mathworld.wolfram.com/Kurtosis.html
0460
0461    RealType df = dist.degrees_of_freedom();
0462    if(((boost::math::isnan)(df)) || (df <= 4))
0463    { // Undefined or infinity for moment k = 4.
0464      return policies::raise_domain_error<RealType>(
0465        "boost::math::kurtosis_excess(students_t_distribution<%1%> const&, %1%)",
0466        "Kurtosis_excess is undefined for degrees of freedom <= 4, but got %1%.",
0467       df, Policy());
0468      return boost::math::numeric_limits<RealType>::quiet_NaN(); // Undefined.
0469    }
0470    if ((boost::math::isinf)(df))
0471    { // +infinity.
0472      return 0;
0473    }
0474    RealType limit = policies::get_epsilon<RealType, Policy>();
0475    // Use policies so that if policy requests lower precision,
0476    // then get the normal distribution approximation earlier.
0477    limit = static_cast<RealType>(1) / limit; // 1/eps
0478    // for 64-bit double 1/eps = 4503599627370496
0479    if (df > limit)
0480    { // Special case for really big degrees_of_freedom > 1 / eps.
0481      return 0;
0482    }
0483    else
0484    {
0485      return 6 / (df - 4);
0486    }
0487 }
0488
0489 template <class RealType, class Policy>
0490 BOOST_MATH_GPU_ENABLED inline RealType entropy(const students_t_distribution<RealType, Policy>& dist)
0491 {
0492    BOOST_MATH_STD_USING
0493    RealType v = dist.degrees_of_freedom();
0494    RealType vp1 = (v+1)/2;
0495    RealType vd2 = v/2;
0496
0497    return vp1*(digamma(vp1) - digamma(vd2)) + log(sqrt(v)*beta(vd2, RealType(1)/RealType(2)));
0498 }
0499
0500 } // namespace math
0501 } // namespace boost
0502
0503 #ifdef _MSC_VER
0504 # pragma warning(pop)
0505 #endif
0506
0507 // This include must be at the end, *after* the accessors
0508 // for this distribution have been defined, in order to
0509 // keep compilers that support two-phase lookup happy.
0510 #include <boost/math/distributions/detail/derived_accessors.hpp>
0511
0512 #endif // BOOST_STATS_STUDENTS_T_HPP