ceph/src/boost/boost/math/distributions/students_t.hpp

   1 //  Copyright John Maddock 2006.
   2 //  Copyright Paul A. Bristow 2006, 2012, 2017.
   3 //  Copyright Thomas Mang 2012.
   4
   5 //  Use, modification and distribution are subject to the
   6 //  Boost Software License, Version 1.0. (See accompanying file
   7 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   8
   9 #ifndef BOOST_STATS_STUDENTS_T_HPP
  10 #define BOOST_STATS_STUDENTS_T_HPP
  11
  12 // http://en.wikipedia.org/wiki/Student%27s_t_distribution
  13 // http://www.itl.nist.gov/div898/handbook/eda/section3/eda3664.htm
  14
  15 #include <boost/math/distributions/fwd.hpp>
  16 #include <boost/math/special_functions/beta.hpp> // for ibeta(a, b, x).
  17 #include <boost/math/distributions/complement.hpp>
  18 #include <boost/math/distributions/detail/common_error_handling.hpp>
  19 #include <boost/math/distributions/normal.hpp>
  20
  21 #include <utility>
  22
  23 #ifdef BOOST_MSVC
  24 # pragma warning(push)
  25 # pragma warning(disable: 4702) // unreachable code (return after domain_error throw).
  26 #endif
  27
  28 namespace boost { namespace math {
  29
  30 template <class RealType = double, class Policy = policies::policy<> >
  31 class students_t_distribution
  32 {
  33 public:
  34    typedef RealType value_type;
  35    typedef Policy policy_type;
  36
  37    students_t_distribution(RealType df) : df_(df)
  38    { // Constructor.
  39       RealType result;
  40       detail::check_df_gt0_to_inf( // Checks that df > 0 or df == inf.
  41          "boost::math::students_t_distribution<%1%>::students_t_distribution", df_, &result, Policy());
  42    } // students_t_distribution
  43
  44    RealType degrees_of_freedom()const
  45    {
  46       return df_;
  47    }
  48
  49    // Parameter estimation:
  50    static RealType find_degrees_of_freedom(
  51       RealType difference_from_mean,
  52       RealType alpha,
  53       RealType beta,
  54       RealType sd,
  55       RealType hint = 100);
  56
  57 private:
  58    // Data member:
  59    RealType df_;  // degrees of freedom is a real number > 0 or +infinity.
  60 };
  61
  62 typedef students_t_distribution<double> students_t; // Convenience typedef for double version.
  63
  64 template <class RealType, class Policy>
  65 inline const std::pair<RealType, RealType> range(const students_t_distribution<RealType, Policy>& /*dist*/)
  66 { // Range of permissible values for random variable x.
  67   // Now including infinity.
  68    using boost::math::tools::max_value;
  69    //return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>());
  70    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>()));
  71 }
  72
  73 template <class RealType, class Policy>
  74 inline const std::pair<RealType, RealType> support(const students_t_distribution<RealType, Policy>& /*dist*/)
  75 { // Range of supported values for random variable x.
  76   // Now including infinity.
  77    // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
  78    using boost::math::tools::max_value;
  79    //return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>());
  80    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>()));
  81 }
  82
  83 template <class RealType, class Policy>
  84 inline RealType pdf(const students_t_distribution<RealType, Policy>& dist, const RealType& x)
  85 {
  86    BOOST_FPU_EXCEPTION_GUARD
  87    BOOST_MATH_STD_USING  // for ADL of std functions.
  88
  89    RealType error_result;
  90    if(false == detail::check_x_not_NaN(
  91       "boost::math::pdf(const students_t_distribution<%1%>&, %1%)", x, &error_result, Policy()))
  92       return error_result;
  93    RealType df = dist.degrees_of_freedom();
  94    if(false == detail::check_df_gt0_to_inf( // Check that df > 0 or == +infinity.
  95       "boost::math::pdf(const students_t_distribution<%1%>&, %1%)", df, &error_result, Policy()))
  96       return error_result;
  97
  98    RealType result;
  99    if ((boost::math::isinf)(x))
 100    { // - or +infinity.
 101      result = static_cast<RealType>(0);
 102      return result;
 103    }
 104    RealType limit = policies::get_epsilon<RealType, Policy>();
 105    // Use policies so that if policy requests lower precision,
 106    // then get the normal distribution approximation earlier.
 107    limit = static_cast<RealType>(1) / limit; // 1/eps
 108    // for 64-bit double 1/eps = 4503599627370496
 109    if (df > limit)
 110    { // Special case for really big degrees_of_freedom > 1 / eps
 111      // - use normal distribution which is much faster and more accurate.
 112      normal_distribution<RealType, Policy> n(0, 1);
 113      result = pdf(n, x);
 114    }
 115    else
 116    { //
 117      RealType basem1 = x * x / df;
 118      if(basem1 < 0.125)
 119      {
 120         result = exp(-boost::math::log1p(basem1, Policy()) * (1+df) / 2);
 121      }
 122      else
 123      {
 124         result = pow(1 / (1 + basem1), (df + 1) / 2);
 125      }
 126      result /= sqrt(df) * boost::math::beta(df / 2, RealType(0.5f), Policy());
 127    }
 128    return result;
 129 } // pdf
 130
 131 template <class RealType, class Policy>
 132 inline RealType cdf(const students_t_distribution<RealType, Policy>& dist, const RealType& x)
 133 {
 134    RealType error_result;
 135    // degrees_of_freedom > 0 or infinity check:
 136    RealType df = dist.degrees_of_freedom();
 137    if (false == detail::check_df_gt0_to_inf(  // Check that df > 0 or == +infinity.
 138      "boost::math::cdf(const students_t_distribution<%1%>&, %1%)", df, &error_result, Policy()))
 139    {
 140      return error_result;
 141    }
 142    // Check for bad x first.
 143    if(false == detail::check_x_not_NaN(
 144       "boost::math::cdf(const students_t_distribution<%1%>&, %1%)", x, &error_result, Policy()))
 145    {
 146       return error_result;
 147    }
 148    if (x == 0)
 149    { // Special case with exact result.
 150      return static_cast<RealType>(0.5);
 151    }
 152    if ((boost::math::isinf)(x))
 153    { // x == - or + infinity, regardless of df.
 154      return ((x < 0) ? static_cast<RealType>(0) : static_cast<RealType>(1));
 155    }
 156
 157    RealType limit = policies::get_epsilon<RealType, Policy>();
 158    // Use policies so that if policy requests lower precision,
 159    // then get the normal distribution approximation earlier.
 160    limit = static_cast<RealType>(1) / limit; // 1/eps
 161    // for 64-bit double 1/eps = 4503599627370496
 162    if (df > limit)
 163    { // Special case for really big degrees_of_freedom > 1 / eps (perhaps infinite?)
 164      // - use normal distribution which is much faster and more accurate.
 165      normal_distribution<RealType, Policy> n(0, 1);
 166      RealType result = cdf(n, x);
 167      return result;
 168    }
 169    else
 170    { // normal df case.
 171      //
 172      // Calculate probability of Student's t using the incomplete beta function.
 173      // probability = ibeta(degrees_of_freedom / 2, 1/2, degrees_of_freedom / (degrees_of_freedom + t*t))
 174      //
 175      // However when t is small compared to the degrees of freedom, that formula
 176      // suffers from rounding error, use the identity formula to work around
 177      // the problem:
 178      //
 179      // I[x](a,b) = 1 - I[1-x](b,a)
 180      //
 181      // and:
 182      //
 183      //     x = df / (df + t^2)
 184      //
 185      // so:
 186      //
 187      // 1 - x = t^2 / (df + t^2)
 188      //
 189      RealType x2 = x * x;
 190      RealType probability;
 191      if(df > 2 * x2)
 192      {
 193         RealType z = x2 / (df + x2);
 194         probability = ibetac(static_cast<RealType>(0.5), df / 2, z, Policy()) / 2;
 195      }
 196      else
 197      {
 198         RealType z = df / (df + x2);
 199         probability = ibeta(df / 2, static_cast<RealType>(0.5), z, Policy()) / 2;
 200      }
 201      return (x > 0 ? 1   - probability : probability);
 202   }
 203 } // cdf
 204
 205 template <class RealType, class Policy>
 206 inline RealType quantile(const students_t_distribution<RealType, Policy>& dist, const RealType& p)
 207 {
 208    BOOST_MATH_STD_USING // for ADL of std functions
 209    //
 210    // Obtain parameters:
 211    RealType probability = p;
 212
 213    // Check for domain errors:
 214    RealType df = dist.degrees_of_freedom();
 215    static const char* function = "boost::math::quantile(const students_t_distribution<%1%>&, %1%)";
 216    RealType error_result;
 217    if(false == (detail::check_df_gt0_to_inf( // Check that df > 0 or == +infinity.
 218       function, df, &error_result, Policy())
 219          && detail::check_probability(function, probability, &error_result, Policy())))
 220       return error_result;
 221    // Special cases, regardless of degrees_of_freedom.
 222    if (probability == 0)
 223       return -policies::raise_overflow_error<RealType>(function, 0, Policy());
 224    if (probability == 1)
 225      return policies::raise_overflow_error<RealType>(function, 0, Policy());
 226    if (probability == static_cast<RealType>(0.5))
 227      return 0;  //
 228    //
 229 #if 0
 230    // This next block is disabled in favour of a faster method than
 231    // incomplete beta inverse, but code retained for future reference:
 232    //
 233    // Calculate quantile of Student's t using the incomplete beta function inverse:
 234    probability = (probability > 0.5) ? 1 - probability : probability;
 235    RealType t, x, y;
 236    x = ibeta_inv(degrees_of_freedom / 2, RealType(0.5), 2 * probability, &y);
 237    if(degrees_of_freedom * y > tools::max_value<RealType>() * x)
 238       t = tools::overflow_error<RealType>(function);
 239    else
 240       t = sqrt(degrees_of_freedom * y / x);
 241    //
 242    // Figure out sign based on the size of p:
 243    //
 244    if(p < 0.5)
 245       t = -t;
 246
 247    return t;
 248 #endif
 249    //
 250    // Depending on how many digits RealType has, this may forward
 251    // to the incomplete beta inverse as above.  Otherwise uses a
 252    // faster method that is accurate to ~15 digits everywhere
 253    // and a couple of epsilon at double precision and in the central
 254    // region where most use cases will occur...
 255    //
 256    return boost::math::detail::fast_students_t_quantile(df, probability, Policy());
 257 } // quantile
 258
 259 template <class RealType, class Policy>
 260 inline RealType cdf(const complemented2_type<students_t_distribution<RealType, Policy>, RealType>& c)
 261 {
 262    return cdf(c.dist, -c.param);
 263 }
 264
 265 template <class RealType, class Policy>
 266 inline RealType quantile(const complemented2_type<students_t_distribution<RealType, Policy>, RealType>& c)
 267 {
 268    return -quantile(c.dist, c.param);
 269 }
 270
 271 //
 272 // Parameter estimation follows:
 273 //
 274 namespace detail{
 275 //
 276 // Functors for finding degrees of freedom:
 277 //
 278 template <class RealType, class Policy>
 279 struct sample_size_func
 280 {
 281    sample_size_func(RealType a, RealType b, RealType s, RealType d)
 282       : alpha(a), beta(b), ratio(s*s/(d*d)) {}
 283
 284    RealType operator()(const RealType& df)
 285    {
 286       if(df <= tools::min_value<RealType>())
 287       { //
 288          return 1;
 289       }
 290       students_t_distribution<RealType, Policy> t(df);
 291       RealType qa = quantile(complement(t, alpha));
 292       RealType qb = quantile(complement(t, beta));
 293       qa += qb;
 294       qa *= qa;
 295       qa *= ratio;
 296       qa -= (df + 1);
 297       return qa;
 298    }
 299    RealType alpha, beta, ratio;
 300 };
 301
 302 }  // namespace detail
 303
 304 template <class RealType, class Policy>
 305 RealType students_t_distribution<RealType, Policy>::find_degrees_of_freedom(
 306       RealType difference_from_mean,
 307       RealType alpha,
 308       RealType beta,
 309       RealType sd,
 310       RealType hint)
 311 {
 312    static const char* function = "boost::math::students_t_distribution<%1%>::find_degrees_of_freedom";
 313    //
 314    // Check for domain errors:
 315    //
 316    RealType error_result;
 317    if(false == detail::check_probability(
 318       function, alpha, &error_result, Policy())
 319          && detail::check_probability(function, beta, &error_result, Policy()))
 320       return error_result;
 321
 322    if(hint <= 0)
 323       hint = 1;
 324
 325    detail::sample_size_func<RealType, Policy> f(alpha, beta, sd, difference_from_mean);
 326    tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
 327    boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
 328    std::pair<RealType, RealType> r = tools::bracket_and_solve_root(f, hint, RealType(2), false, tol, max_iter, Policy());
 329    RealType result = r.first + (r.second - r.first) / 2;
 330    if(max_iter >= policies::get_max_root_iterations<Policy>())
 331    {
 332       return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:"
 333          " either there is no answer to how many degrees of freedom are required"
 334          " or the answer is infinite.  Current best guess is %1%", result, Policy());
 335    }
 336    return result;
 337 }
 338
 339 template <class RealType, class Policy>
 340 inline RealType mode(const students_t_distribution<RealType, Policy>& /*dist*/)
 341 {
 342   // Assume no checks on degrees of freedom are useful (unlike mean).
 343    return 0; // Always zero by definition.
 344 }
 345
 346 template <class RealType, class Policy>
 347 inline RealType median(const students_t_distribution<RealType, Policy>& /*dist*/)
 348 {
 349    // Assume no checks on degrees of freedom are useful (unlike mean).
 350    return 0; // Always zero by definition.
 351 }
 352
 353 // See section 5.1 on moments at  http://en.wikipedia.org/wiki/Student%27s_t-distribution
 354
 355 template <class RealType, class Policy>
 356 inline RealType mean(const students_t_distribution<RealType, Policy>& dist)
 357 {  // Revised for https://svn.boost.org/trac/boost/ticket/7177
 358    RealType df = dist.degrees_of_freedom();
 359    if(((boost::math::isnan)(df)) || (df <= 1) )
 360    { // mean is undefined for moment <= 1!
 361       return policies::raise_domain_error<RealType>(
 362       "boost::math::mean(students_t_distribution<%1%> const&, %1%)",
 363       "Mean is undefined for degrees of freedom < 1 but got %1%.", df, Policy());
 364       return std::numeric_limits<RealType>::quiet_NaN();
 365    }
 366    return 0;
 367 } // mean
 368
 369 template <class RealType, class Policy>
 370 inline RealType variance(const students_t_distribution<RealType, Policy>& dist)
 371 { // http://en.wikipedia.org/wiki/Student%27s_t-distribution
 372   // Revised for https://svn.boost.org/trac/boost/ticket/7177
 373   RealType df = dist.degrees_of_freedom();
 374   if ((boost::math::isnan)(df) || (df <= 2))
 375   { // NaN or undefined for <= 2.
 376      return policies::raise_domain_error<RealType>(
 377       "boost::math::variance(students_t_distribution<%1%> const&, %1%)",
 378       "variance is undefined for degrees of freedom <= 2, but got %1%.",
 379       df, Policy());
 380     return std::numeric_limits<RealType>::quiet_NaN(); // Undefined.
 381   }
 382   if ((boost::math::isinf)(df))
 383   { // +infinity.
 384     return 1;
 385   }
 386   RealType limit = policies::get_epsilon<RealType, Policy>();
 387   // Use policies so that if policy requests lower precision,
 388   // then get the normal distribution approximation earlier.
 389   limit = static_cast<RealType>(1) / limit; // 1/eps
 390   // for 64-bit double 1/eps = 4503599627370496
 391   if (df > limit)
 392   { // Special case for really big degrees_of_freedom > 1 / eps.
 393     return 1;
 394   }
 395   else
 396   {
 397     return df / (df - 2);
 398   }
 399 } // variance
 400
 401 template <class RealType, class Policy>
 402 inline RealType skewness(const students_t_distribution<RealType, Policy>& dist)
 403 {
 404     RealType df = dist.degrees_of_freedom();
 405    if( ((boost::math::isnan)(df)) || (dist.degrees_of_freedom() <= 3))
 406    { // Undefined for moment k = 3.
 407       return policies::raise_domain_error<RealType>(
 408          "boost::math::skewness(students_t_distribution<%1%> const&, %1%)",
 409          "Skewness is undefined for degrees of freedom <= 3, but got %1%.",
 410          dist.degrees_of_freedom(), Policy());
 411       return std::numeric_limits<RealType>::quiet_NaN();
 412    }
 413    return 0; // For all valid df, including infinity.
 414 } // skewness
 415
 416 template <class RealType, class Policy>
 417 inline RealType kurtosis(const students_t_distribution<RealType, Policy>& dist)
 418 {
 419    RealType df = dist.degrees_of_freedom();
 420    if(((boost::math::isnan)(df)) || (df <= 4))
 421    { // Undefined or infinity for moment k = 4.
 422       return policies::raise_domain_error<RealType>(
 423        "boost::math::kurtosis(students_t_distribution<%1%> const&, %1%)",
 424        "Kurtosis is undefined for degrees of freedom <= 4, but got %1%.",
 425         df, Policy());
 426         return std::numeric_limits<RealType>::quiet_NaN(); // Undefined.
 427    }
 428    if ((boost::math::isinf)(df))
 429    { // +infinity.
 430      return 3;
 431    }
 432    RealType limit = policies::get_epsilon<RealType, Policy>();
 433    // Use policies so that if policy requests lower precision,
 434    // then get the normal distribution approximation earlier.
 435    limit = static_cast<RealType>(1) / limit; // 1/eps
 436    // for 64-bit double 1/eps = 4503599627370496
 437    if (df > limit)
 438    { // Special case for really big degrees_of_freedom > 1 / eps.
 439      return 3;
 440    }
 441    else
 442    {
 443      //return 3 * (df - 2) / (df - 4); re-arranged to
 444      return 6 / (df - 4) + 3;
 445    }
 446 } // kurtosis
 447
 448 template <class RealType, class Policy>
 449 inline RealType kurtosis_excess(const students_t_distribution<RealType, Policy>& dist)
 450 {
 451    // see http://mathworld.wolfram.com/Kurtosis.html
 452
 453    RealType df = dist.degrees_of_freedom();
 454    if(((boost::math::isnan)(df)) || (df <= 4))
 455    { // Undefined or infinity for moment k = 4.
 456      return policies::raise_domain_error<RealType>(
 457        "boost::math::kurtosis_excess(students_t_distribution<%1%> const&, %1%)",
 458        "Kurtosis_excess is undefined for degrees of freedom <= 4, but got %1%.",
 459       df, Policy());
 460      return std::numeric_limits<RealType>::quiet_NaN(); // Undefined.
 461    }
 462    if ((boost::math::isinf)(df))
 463    { // +infinity.
 464      return 0;
 465    }
 466    RealType limit = policies::get_epsilon<RealType, Policy>();
 467    // Use policies so that if policy requests lower precision,
 468    // then get the normal distribution approximation earlier.
 469    limit = static_cast<RealType>(1) / limit; // 1/eps
 470    // for 64-bit double 1/eps = 4503599627370496
 471    if (df > limit)
 472    { // Special case for really big degrees_of_freedom > 1 / eps.
 473      return 0;
 474    }
 475    else
 476    {
 477      return 6 / (df - 4);
 478    }
 479 }
 480
 481 } // namespace math
 482 } // namespace boost
 483
 484 #ifdef BOOST_MSVC
 485 # pragma warning(pop)
 486 #endif
 487
 488 // This include must be at the end, *after* the accessors
 489 // for this distribution have been defined, in order to
 490 // keep compilers that support two-phase lookup happy.
 491 #include <boost/math/distributions/detail/derived_accessors.hpp>
 492
 493 #endif // BOOST_STATS_STUDENTS_T_HPP