]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // Copyright John Maddock 2006. |
2 | // Copyright Paul A. Bristow 2006, 2012. | |
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 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 | // NOT including infinity. | |
68 | using boost::math::tools::max_value; | |
69 | return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>()); | |
70 | } | |
71 | ||
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>()); | |
78 | } | |
79 | ||
80 | template <class RealType, class Policy> | |
81 | inline RealType pdf(const students_t_distribution<RealType, Policy>& dist, const RealType& x) | |
82 | { | |
83 | BOOST_FPU_EXCEPTION_GUARD | |
84 | BOOST_MATH_STD_USING // for ADL of std functions. | |
85 | ||
86 | RealType error_result; | |
87 | if(false == detail::check_x( | |
88 | "boost::math::pdf(const students_t_distribution<%1%>&, %1%)", x, &error_result, Policy())) | |
89 | return error_result; | |
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())) | |
93 | return error_result; | |
94 | ||
95 | RealType result; | |
96 | if ((boost::math::isinf)(x)) | |
97 | { // +infinity. | |
98 | normal_distribution<RealType, Policy> n(0, 1); | |
99 | result = pdf(n, x); | |
100 | return result; | |
101 | } | |
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 | |
107 | if (df > limit) | |
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); | |
111 | result = pdf(n, x); | |
112 | } | |
113 | else | |
114 | { // | |
115 | RealType basem1 = x * x / df; | |
116 | if(basem1 < 0.125) | |
117 | { | |
118 | result = exp(-boost::math::log1p(basem1, Policy()) * (1+df) / 2); | |
119 | } | |
120 | else | |
121 | { | |
122 | result = pow(1 / (1 + basem1), (df + 1) / 2); | |
123 | } | |
124 | result /= sqrt(df) * boost::math::beta(df / 2, RealType(0.5f), Policy()); | |
125 | } | |
126 | return result; | |
127 | ||
128 | ||
129 | template <class RealType, class Policy> | |
130 | inline RealType cdf(const students_t_distribution<RealType, Policy>& dist, const RealType& x) | |
131 | { | |
132 | RealType error_result; | |
133 | if(false == detail::check_x( | |
134 | "boost::math::pdf(const students_t_distribution<%1%>&, %1%)", x, &error_result, Policy())) | |
135 | return error_result; | |
136 | RealType df = dist.degrees_of_freedom(); | |
137 | // Error check: | |
138 | ||
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())) | |
141 | return error_result; | |
142 | ||
143 | if (x == 0) | |
144 | { // Special case with exact result. | |
145 | return static_cast<RealType>(0.5); | |
146 | } | |
147 | if ((boost::math::isinf)(x)) | |
148 | { // +infinity. | |
149 | normal_distribution<RealType, Policy> n(0, 1); | |
150 | RealType result = cdf(n, x); | |
151 | return result; | |
152 | } | |
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 | |
158 | if (df > limit) | |
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); | |
163 | return result; | |
164 | } | |
165 | else | |
166 | { // normal df case. | |
167 | // | |
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)) | |
170 | // | |
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 | |
173 | // the problem: | |
174 | // | |
175 | // I[x](a,b) = 1 - I[1-x](b,a) | |
176 | // | |
177 | // and: | |
178 | // | |
179 | // x = df / (df + t^2) | |
180 | // | |
181 | // so: | |
182 | // | |
183 | // 1 - x = t^2 / (df + t^2) | |
184 | // | |
185 | RealType x2 = x * x; | |
186 | RealType probability; | |
187 | if(df > 2 * x2) | |
188 | { | |
189 | RealType z = x2 / (df + x2); | |
190 | probability = ibetac(static_cast<RealType>(0.5), df / 2, z, Policy()) / 2; | |
191 | } | |
192 | else | |
193 | { | |
194 | RealType z = df / (df + x2); | |
195 | probability = ibeta(df / 2, static_cast<RealType>(0.5), z, Policy()) / 2; | |
196 | } | |
197 | return (x > 0 ? 1 - probability : probability); | |
198 | } | |
199 | } // cdf | |
200 | ||
201 | template <class RealType, class Policy> | |
202 | inline RealType quantile(const students_t_distribution<RealType, Policy>& dist, const RealType& p) | |
203 | { | |
204 | BOOST_MATH_STD_USING // for ADL of std functions | |
205 | // | |
206 | // Obtain parameters: | |
207 | RealType probability = p; | |
208 | ||
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()))) | |
216 | return error_result; | |
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)) | |
223 | return 0; // | |
224 | // | |
225 | #if 0 | |
226 | // This next block is disabled in favour of a faster method than | |
227 | // incomplete beta inverse, but code retained for future reference: | |
228 | // | |
229 | // Calculate quantile of Student's t using the incomplete beta function inverse: | |
230 | // | |
231 | probability = (probability > 0.5) ? 1 - probability : probability; | |
232 | RealType t, x, y; | |
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); | |
236 | else | |
237 | t = sqrt(degrees_of_freedom * y / x); | |
238 | // | |
239 | // Figure out sign based on the size of p: | |
240 | // | |
241 | if(p < 0.5) | |
242 | t = -t; | |
243 | ||
244 | return t; | |
245 | #endif | |
246 | // | |
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... | |
252 | // | |
253 | return boost::math::detail::fast_students_t_quantile(df, probability, Policy()); | |
254 | } // quantile | |
255 | ||
256 | template <class RealType, class Policy> | |
257 | inline RealType cdf(const complemented2_type<students_t_distribution<RealType, Policy>, RealType>& c) | |
258 | { | |
259 | return cdf(c.dist, -c.param); | |
260 | } | |
261 | ||
262 | template <class RealType, class Policy> | |
263 | inline RealType quantile(const complemented2_type<students_t_distribution<RealType, Policy>, RealType>& c) | |
264 | { | |
265 | return -quantile(c.dist, c.param); | |
266 | } | |
267 | ||
268 | // | |
269 | // Parameter estimation follows: | |
270 | // | |
271 | namespace detail{ | |
272 | // | |
273 | // Functors for finding degrees of freedom: | |
274 | // | |
275 | template <class RealType, class Policy> | |
276 | struct sample_size_func | |
277 | { | |
278 | sample_size_func(RealType a, RealType b, RealType s, RealType d) | |
279 | : alpha(a), beta(b), ratio(s*s/(d*d)) {} | |
280 | ||
281 | RealType operator()(const RealType& df) | |
282 | { | |
283 | if(df <= tools::min_value<RealType>()) | |
284 | { // | |
285 | return 1; | |
286 | } | |
287 | students_t_distribution<RealType, Policy> t(df); | |
288 | RealType qa = quantile(complement(t, alpha)); | |
289 | RealType qb = quantile(complement(t, beta)); | |
290 | qa += qb; | |
291 | qa *= qa; | |
292 | qa *= ratio; | |
293 | qa -= (df + 1); | |
294 | return qa; | |
295 | } | |
296 | RealType alpha, beta, ratio; | |
297 | }; | |
298 | ||
299 | } // namespace detail | |
300 | ||
301 | template <class RealType, class Policy> | |
302 | RealType students_t_distribution<RealType, Policy>::find_degrees_of_freedom( | |
303 | RealType difference_from_mean, | |
304 | RealType alpha, | |
305 | RealType beta, | |
306 | RealType sd, | |
307 | RealType hint) | |
308 | { | |
309 | static const char* function = "boost::math::students_t_distribution<%1%>::find_degrees_of_freedom"; | |
310 | // | |
311 | // Check for domain errors: | |
312 | // | |
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())) | |
317 | return error_result; | |
318 | ||
319 | if(hint <= 0) | |
320 | hint = 1; | |
321 | ||
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>()) | |
328 | { | |
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()); | |
332 | } | |
333 | return result; | |
334 | } | |
335 | ||
336 | template <class RealType, class Policy> | |
337 | inline RealType mode(const students_t_distribution<RealType, Policy>& /*dist*/) | |
338 | { | |
339 | // Assume no checks on degrees of freedom are useful (unlike mean). | |
340 | return 0; // Always zero by definition. | |
341 | } | |
342 | ||
343 | template <class RealType, class Policy> | |
344 | inline RealType median(const students_t_distribution<RealType, Policy>& /*dist*/) | |
345 | { | |
346 | // Assume no checks on degrees of freedom are useful (unlike mean). | |
347 | return 0; // Always zero by definition. | |
348 | } | |
349 | ||
350 | // See section 5.1 on moments at http://en.wikipedia.org/wiki/Student%27s_t-distribution | |
351 | ||
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(); | |
362 | } | |
363 | return 0; | |
364 | } // mean | |
365 | ||
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%.", | |
376 | df, Policy()); | |
377 | return std::numeric_limits<RealType>::quiet_NaN(); // Undefined. | |
378 | } | |
379 | if ((boost::math::isinf)(df)) | |
380 | { // +infinity. | |
381 | return 1; | |
382 | } | |
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 | |
388 | if (df > limit) | |
389 | { // Special case for really big degrees_of_freedom > 1 / eps. | |
390 | return 1; | |
391 | } | |
392 | else | |
393 | { | |
394 | return df / (df - 2); | |
395 | } | |
396 | } // variance | |
397 | ||
398 | template <class RealType, class Policy> | |
399 | inline RealType skewness(const students_t_distribution<RealType, Policy>& dist) | |
400 | { | |
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(); | |
409 | } | |
410 | return 0; // For all valid df, including infinity. | |
411 | } // skewness | |
412 | ||
413 | template <class RealType, class Policy> | |
414 | inline RealType kurtosis(const students_t_distribution<RealType, Policy>& dist) | |
415 | { | |
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%.", | |
422 | df, Policy()); | |
423 | return std::numeric_limits<RealType>::quiet_NaN(); // Undefined. | |
424 | } | |
425 | if ((boost::math::isinf)(df)) | |
426 | { // +infinity. | |
427 | return 3; | |
428 | } | |
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 | |
434 | if (df > limit) | |
435 | { // Special case for really big degrees_of_freedom > 1 / eps. | |
436 | return 3; | |
437 | } | |
438 | else | |
439 | { | |
440 | //return 3 * (df - 2) / (df - 4); re-arranged to | |
441 | return 6 / (df - 4) + 3; | |
442 | } | |
443 | } // kurtosis | |
444 | ||
445 | template <class RealType, class Policy> | |
446 | inline RealType kurtosis_excess(const students_t_distribution<RealType, Policy>& dist) | |
447 | { | |
448 | // see http://mathworld.wolfram.com/Kurtosis.html | |
449 | ||
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%.", | |
456 | df, Policy()); | |
457 | return std::numeric_limits<RealType>::quiet_NaN(); // Undefined. | |
458 | } | |
459 | if ((boost::math::isinf)(df)) | |
460 | { // +infinity. | |
461 | return 0; | |
462 | } | |
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 | |
468 | if (df > limit) | |
469 | { // Special case for really big degrees_of_freedom > 1 / eps. | |
470 | return 0; | |
471 | } | |
472 | else | |
473 | { | |
474 | return 6 / (df - 4); | |
475 | } | |
476 | } | |
477 | ||
478 | } // namespace math | |
479 | } // namespace boost | |
480 | ||
481 | #ifdef BOOST_MSVC | |
482 | # pragma warning(pop) | |
483 | #endif | |
484 | ||
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> | |
489 | ||
490 | #endif // BOOST_STATS_STUDENTS_T_HPP |