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