]>
Commit | Line | Data |
---|---|---|
7c673cae | 1 | // Copyright John Maddock 2006. |
b32b8144 | 2 | // Copyright Paul A. Bristow 2006, 2012, 2017. |
7c673cae FG |
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 | ||
b32b8144 | 28 | namespace boost { namespace math { |
7c673cae FG |
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: | |
b32b8144 | 59 | RealType df_; // degrees of freedom is a real number > 0 or +infinity. |
7c673cae FG |
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. | |
b32b8144 | 67 | // Now including infinity. |
7c673cae | 68 | using boost::math::tools::max_value; |
b32b8144 FG |
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>())); | |
7c673cae FG |
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. | |
b32b8144 | 76 | // Now including infinity. |
7c673cae FG |
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; | |
b32b8144 FG |
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>())); | |
7c673cae FG |
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; | |
b32b8144 | 90 | if(false == detail::check_x_not_NaN( |
7c673cae FG |
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)) | |
b32b8144 FG |
100 | { // - or +infinity. |
101 | result = static_cast<RealType>(0); | |
7c673cae FG |
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 | ||
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; | |
b32b8144 | 135 | // degrees_of_freedom > 0 or infinity check: |
7c673cae | 136 | RealType df = dist.degrees_of_freedom(); |
b32b8144 FG |
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 | { | |
7c673cae | 146 | return error_result; |
b32b8144 | 147 | } |
7c673cae FG |
148 | if (x == 0) |
149 | { // Special case with exact result. | |
150 | return static_cast<RealType>(0.5); | |
151 | } | |
152 | if ((boost::math::isinf)(x)) | |
b32b8144 FG |
153 | { // x == - or + infinity, regardless of df. |
154 | return ((x < 0) ? static_cast<RealType>(0) : static_cast<RealType>(1)); | |
7c673cae | 155 | } |
b32b8144 | 156 | |
7c673cae FG |
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: | |
7c673cae FG |
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 |