1 // Copyright John Maddock 2006.
2 // Copyright Paul A. Bristow 2006, 2012.
3 // Copyright Thomas Mang 2012.
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)
9 #ifndef BOOST_STATS_STUDENTS_T_HPP
10 #define BOOST_STATS_STUDENTS_T_HPP
12 // http://en.wikipedia.org/wiki/Student%27s_t_distribution
13 // http://www.itl.nist.gov/div898/handbook/eda/section3/eda3664.htm
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>
24 # pragma warning(push)
25 # pragma warning(disable: 4702) // unreachable code (return after domain_error throw).
28 namespace boost{ namespace math{
30 template <class RealType = double, class Policy = policies::policy<> >
31 class students_t_distribution
34 typedef RealType value_type;
35 typedef Policy policy_type;
37 students_t_distribution(RealType df) : df_(df)
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
44 RealType degrees_of_freedom()const
49 // Parameter estimation:
50 static RealType find_degrees_of_freedom(
51 RealType difference_from_mean,
59 RealType df_; // degrees of freedom is a real number or +infinity.
62 typedef students_t_distribution<double> students_t; // Convenience typedef for double version.
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 // NOT including infinity.
68 using boost::math::tools::max_value;
69 return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>());
72 template <class RealType, class Policy>
73 inline const std::pair<RealType, RealType> support(const students_t_distribution<RealType, Policy>& /*dist*/)
74 { // Range of supported values for random variable x.
75 // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
76 using boost::math::tools::max_value;
77 return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>());
80 template <class RealType, class Policy>
81 inline RealType pdf(const students_t_distribution<RealType, Policy>& dist, const RealType& x)
83 BOOST_FPU_EXCEPTION_GUARD
84 BOOST_MATH_STD_USING // for ADL of std functions.
86 RealType error_result;
87 if(false == detail::check_x(
88 "boost::math::pdf(const students_t_distribution<%1%>&, %1%)", x, &error_result, Policy()))
90 RealType df = dist.degrees_of_freedom();
91 if(false == detail::check_df_gt0_to_inf( // Check that df > 0 or == +infinity.
92 "boost::math::pdf(const students_t_distribution<%1%>&, %1%)", df, &error_result, Policy()))
96 if ((boost::math::isinf)(x))
98 normal_distribution<RealType, Policy> n(0, 1);
102 RealType limit = policies::get_epsilon<RealType, Policy>();
103 // Use policies so that if policy requests lower precision,
104 // then get the normal distribution approximation earlier.
105 limit = static_cast<RealType>(1) / limit; // 1/eps
106 // for 64-bit double 1/eps = 4503599627370496
108 { // Special case for really big degrees_of_freedom > 1 / eps
109 // - use normal distribution which is much faster and more accurate.
110 normal_distribution<RealType, Policy> n(0, 1);
115 RealType basem1 = x * x / df;
118 result = exp(-boost::math::log1p(basem1, Policy()) * (1+df) / 2);
122 result = pow(1 / (1 + basem1), (df + 1) / 2);
124 result /= sqrt(df) * boost::math::beta(df / 2, RealType(0.5f), Policy());
129 template <class RealType, class Policy>
130 inline RealType cdf(const students_t_distribution<RealType, Policy>& dist, const RealType& x)
132 RealType error_result;
133 if(false == detail::check_x(
134 "boost::math::pdf(const students_t_distribution<%1%>&, %1%)", x, &error_result, Policy()))
136 RealType df = dist.degrees_of_freedom();
139 if(false == detail::check_df_gt0_to_inf( // Check that df > 0 or == +infinity.
140 "boost::math::cdf(const students_t_distribution<%1%>&, %1%)", df, &error_result, Policy()))
144 { // Special case with exact result.
145 return static_cast<RealType>(0.5);
147 if ((boost::math::isinf)(x))
149 normal_distribution<RealType, Policy> n(0, 1);
150 RealType result = cdf(n, x);
153 RealType limit = policies::get_epsilon<RealType, Policy>();
154 // Use policies so that if policy requests lower precision,
155 // then get the normal distribution approximation earlier.
156 limit = static_cast<RealType>(1) / limit; // 1/eps
157 // for 64-bit double 1/eps = 4503599627370496
159 { // Special case for really big degrees_of_freedom > 1 / eps (perhaps infinite?)
160 // - use normal distribution which is much faster and more accurate.
161 normal_distribution<RealType, Policy> n(0, 1);
162 RealType result = cdf(n, x);
168 // Calculate probability of Student's t using the incomplete beta function.
169 // probability = ibeta(degrees_of_freedom / 2, 1/2, degrees_of_freedom / (degrees_of_freedom + t*t))
171 // However when t is small compared to the degrees of freedom, that formula
172 // suffers from rounding error, use the identity formula to work around
175 // I[x](a,b) = 1 - I[1-x](b,a)
179 // x = df / (df + t^2)
183 // 1 - x = t^2 / (df + t^2)
186 RealType probability;
189 RealType z = x2 / (df + x2);
190 probability = ibetac(static_cast<RealType>(0.5), df / 2, z, Policy()) / 2;
194 RealType z = df / (df + x2);
195 probability = ibeta(df / 2, static_cast<RealType>(0.5), z, Policy()) / 2;
197 return (x > 0 ? 1 - probability : probability);
201 template <class RealType, class Policy>
202 inline RealType quantile(const students_t_distribution<RealType, Policy>& dist, const RealType& p)
204 BOOST_MATH_STD_USING // for ADL of std functions
206 // Obtain parameters:
207 RealType probability = p;
209 // Check for domain errors:
210 RealType df = dist.degrees_of_freedom();
211 static const char* function = "boost::math::quantile(const students_t_distribution<%1%>&, %1%)";
212 RealType error_result;
213 if(false == (detail::check_df_gt0_to_inf( // Check that df > 0 or == +infinity.
214 function, df, &error_result, Policy())
215 && detail::check_probability(function, probability, &error_result, Policy())))
217 // Special cases, regardless of degrees_of_freedom.
218 if (probability == 0)
219 return -policies::raise_overflow_error<RealType>(function, 0, Policy());
220 if (probability == 1)
221 return policies::raise_overflow_error<RealType>(function, 0, Policy());
222 if (probability == static_cast<RealType>(0.5))
226 // This next block is disabled in favour of a faster method than
227 // incomplete beta inverse, but code retained for future reference:
229 // Calculate quantile of Student's t using the incomplete beta function inverse:
231 probability = (probability > 0.5) ? 1 - probability : probability;
233 x = ibeta_inv(degrees_of_freedom / 2, RealType(0.5), 2 * probability, &y);
234 if(degrees_of_freedom * y > tools::max_value<RealType>() * x)
235 t = tools::overflow_error<RealType>(function);
237 t = sqrt(degrees_of_freedom * y / x);
239 // Figure out sign based on the size of p:
247 // Depending on how many digits RealType has, this may forward
248 // to the incomplete beta inverse as above. Otherwise uses a
249 // faster method that is accurate to ~15 digits everywhere
250 // and a couple of epsilon at double precision and in the central
251 // region where most use cases will occur...
253 return boost::math::detail::fast_students_t_quantile(df, probability, Policy());
256 template <class RealType, class Policy>
257 inline RealType cdf(const complemented2_type<students_t_distribution<RealType, Policy>, RealType>& c)
259 return cdf(c.dist, -c.param);
262 template <class RealType, class Policy>
263 inline RealType quantile(const complemented2_type<students_t_distribution<RealType, Policy>, RealType>& c)
265 return -quantile(c.dist, c.param);
269 // Parameter estimation follows:
273 // Functors for finding degrees of freedom:
275 template <class RealType, class Policy>
276 struct sample_size_func
278 sample_size_func(RealType a, RealType b, RealType s, RealType d)
279 : alpha(a), beta(b), ratio(s*s/(d*d)) {}
281 RealType operator()(const RealType& df)
283 if(df <= tools::min_value<RealType>())
287 students_t_distribution<RealType, Policy> t(df);
288 RealType qa = quantile(complement(t, alpha));
289 RealType qb = quantile(complement(t, beta));
296 RealType alpha, beta, ratio;
299 } // namespace detail
301 template <class RealType, class Policy>
302 RealType students_t_distribution<RealType, Policy>::find_degrees_of_freedom(
303 RealType difference_from_mean,
309 static const char* function = "boost::math::students_t_distribution<%1%>::find_degrees_of_freedom";
311 // Check for domain errors:
313 RealType error_result;
314 if(false == detail::check_probability(
315 function, alpha, &error_result, Policy())
316 && detail::check_probability(function, beta, &error_result, Policy()))
322 detail::sample_size_func<RealType, Policy> f(alpha, beta, sd, difference_from_mean);
323 tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
324 boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
325 std::pair<RealType, RealType> r = tools::bracket_and_solve_root(f, hint, RealType(2), false, tol, max_iter, Policy());
326 RealType result = r.first + (r.second - r.first) / 2;
327 if(max_iter >= policies::get_max_root_iterations<Policy>())
329 return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:"
330 " either there is no answer to how many degrees of freedom are required"
331 " or the answer is infinite. Current best guess is %1%", result, Policy());
336 template <class RealType, class Policy>
337 inline RealType mode(const students_t_distribution<RealType, Policy>& /*dist*/)
339 // Assume no checks on degrees of freedom are useful (unlike mean).
340 return 0; // Always zero by definition.
343 template <class RealType, class Policy>
344 inline RealType median(const students_t_distribution<RealType, Policy>& /*dist*/)
346 // Assume no checks on degrees of freedom are useful (unlike mean).
347 return 0; // Always zero by definition.
350 // See section 5.1 on moments at http://en.wikipedia.org/wiki/Student%27s_t-distribution
352 template <class RealType, class Policy>
353 inline RealType mean(const students_t_distribution<RealType, Policy>& dist)
354 { // Revised for https://svn.boost.org/trac/boost/ticket/7177
355 RealType df = dist.degrees_of_freedom();
356 if(((boost::math::isnan)(df)) || (df <= 1) )
357 { // mean is undefined for moment <= 1!
358 return policies::raise_domain_error<RealType>(
359 "boost::math::mean(students_t_distribution<%1%> const&, %1%)",
360 "Mean is undefined for degrees of freedom < 1 but got %1%.", df, Policy());
361 return std::numeric_limits<RealType>::quiet_NaN();
366 template <class RealType, class Policy>
367 inline RealType variance(const students_t_distribution<RealType, Policy>& dist)
368 { // http://en.wikipedia.org/wiki/Student%27s_t-distribution
369 // Revised for https://svn.boost.org/trac/boost/ticket/7177
370 RealType df = dist.degrees_of_freedom();
371 if ((boost::math::isnan)(df) || (df <= 2))
372 { // NaN or undefined for <= 2.
373 return policies::raise_domain_error<RealType>(
374 "boost::math::variance(students_t_distribution<%1%> const&, %1%)",
375 "variance is undefined for degrees of freedom <= 2, but got %1%.",
377 return std::numeric_limits<RealType>::quiet_NaN(); // Undefined.
379 if ((boost::math::isinf)(df))
383 RealType limit = policies::get_epsilon<RealType, Policy>();
384 // Use policies so that if policy requests lower precision,
385 // then get the normal distribution approximation earlier.
386 limit = static_cast<RealType>(1) / limit; // 1/eps
387 // for 64-bit double 1/eps = 4503599627370496
389 { // Special case for really big degrees_of_freedom > 1 / eps.
394 return df / (df - 2);
398 template <class RealType, class Policy>
399 inline RealType skewness(const students_t_distribution<RealType, Policy>& dist)
401 RealType df = dist.degrees_of_freedom();
402 if( ((boost::math::isnan)(df)) || (dist.degrees_of_freedom() <= 3))
403 { // Undefined for moment k = 3.
404 return policies::raise_domain_error<RealType>(
405 "boost::math::skewness(students_t_distribution<%1%> const&, %1%)",
406 "Skewness is undefined for degrees of freedom <= 3, but got %1%.",
407 dist.degrees_of_freedom(), Policy());
408 return std::numeric_limits<RealType>::quiet_NaN();
410 return 0; // For all valid df, including infinity.
413 template <class RealType, class Policy>
414 inline RealType kurtosis(const students_t_distribution<RealType, Policy>& dist)
416 RealType df = dist.degrees_of_freedom();
417 if(((boost::math::isnan)(df)) || (df <= 4))
418 { // Undefined or infinity for moment k = 4.
419 return policies::raise_domain_error<RealType>(
420 "boost::math::kurtosis(students_t_distribution<%1%> const&, %1%)",
421 "Kurtosis is undefined for degrees of freedom <= 4, but got %1%.",
423 return std::numeric_limits<RealType>::quiet_NaN(); // Undefined.
425 if ((boost::math::isinf)(df))
429 RealType limit = policies::get_epsilon<RealType, Policy>();
430 // Use policies so that if policy requests lower precision,
431 // then get the normal distribution approximation earlier.
432 limit = static_cast<RealType>(1) / limit; // 1/eps
433 // for 64-bit double 1/eps = 4503599627370496
435 { // Special case for really big degrees_of_freedom > 1 / eps.
440 //return 3 * (df - 2) / (df - 4); re-arranged to
441 return 6 / (df - 4) + 3;
445 template <class RealType, class Policy>
446 inline RealType kurtosis_excess(const students_t_distribution<RealType, Policy>& dist)
448 // see http://mathworld.wolfram.com/Kurtosis.html
450 RealType df = dist.degrees_of_freedom();
451 if(((boost::math::isnan)(df)) || (df <= 4))
452 { // Undefined or infinity for moment k = 4.
453 return policies::raise_domain_error<RealType>(
454 "boost::math::kurtosis_excess(students_t_distribution<%1%> const&, %1%)",
455 "Kurtosis_excess is undefined for degrees of freedom <= 4, but got %1%.",
457 return std::numeric_limits<RealType>::quiet_NaN(); // Undefined.
459 if ((boost::math::isinf)(df))
463 RealType limit = policies::get_epsilon<RealType, Policy>();
464 // Use policies so that if policy requests lower precision,
465 // then get the normal distribution approximation earlier.
466 limit = static_cast<RealType>(1) / limit; // 1/eps
467 // for 64-bit double 1/eps = 4503599627370496
469 { // Special case for really big degrees_of_freedom > 1 / eps.
482 # pragma warning(pop)
485 // This include must be at the end, *after* the accessors
486 // for this distribution have been defined, in order to
487 // keep compilers that support two-phase lookup happy.
488 #include <boost/math/distributions/detail/derived_accessors.hpp>
490 #endif // BOOST_STATS_STUDENTS_T_HPP