math/distributions/students_t.hpp

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