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