]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // boost\math\distributions\non_central_t.hpp |
2 | ||
3 | // Copyright John Maddock 2008. | |
4 | ||
5 | // Use, modification and distribution are subject to the | |
6 | // Boost Software License, Version 1.0. | |
7 | // (See accompanying file LICENSE_1_0.txt | |
8 | // or copy at http://www.boost.org/LICENSE_1_0.txt) | |
9 | ||
10 | #ifndef BOOST_MATH_SPECIAL_NON_CENTRAL_T_HPP | |
11 | #define BOOST_MATH_SPECIAL_NON_CENTRAL_T_HPP | |
12 | ||
13 | #include <boost/math/distributions/fwd.hpp> | |
14 | #include <boost/math/distributions/non_central_beta.hpp> // for nc beta | |
15 | #include <boost/math/distributions/normal.hpp> // for normal CDF and quantile | |
16 | #include <boost/math/distributions/students_t.hpp> | |
17 | #include <boost/math/distributions/detail/generic_quantile.hpp> // quantile | |
18 | ||
19 | namespace boost | |
20 | { | |
21 | namespace math | |
22 | { | |
23 | ||
24 | template <class RealType, class Policy> | |
25 | class non_central_t_distribution; | |
26 | ||
27 | namespace detail{ | |
28 | ||
29 | template <class T, class Policy> | |
30 | T non_central_t2_p(T v, T delta, T x, T y, const Policy& pol, T init_val) | |
31 | { | |
32 | BOOST_MATH_STD_USING | |
33 | // | |
34 | // Variables come first: | |
35 | // | |
36 | boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); | |
37 | T errtol = policies::get_epsilon<T, Policy>(); | |
38 | T d2 = delta * delta / 2; | |
39 | // | |
40 | // k is the starting point for iteration, and is the | |
41 | // maximum of the poisson weighting term, we don't | |
42 | // ever allow k == 0 as this can lead to catastrophic | |
43 | // cancellation errors later (test case is v = 1621286869049072.3 | |
44 | // delta = 0.16212868690490723, x = 0.86987415482475994). | |
45 | // | |
46 | int k = itrunc(d2); | |
47 | T pois; | |
48 | if(k == 0) k = 1; | |
49 | // Starting Poisson weight: | |
50 | pois = gamma_p_derivative(T(k+1), d2, pol) | |
51 | * tgamma_delta_ratio(T(k + 1), T(0.5f)) | |
52 | * delta / constants::root_two<T>(); | |
53 | if(pois == 0) | |
54 | return init_val; | |
55 | T xterm, beta; | |
56 | // Recurrance & starting beta terms: | |
57 | beta = x < y | |
58 | ? detail::ibeta_imp(T(k + 1), T(v / 2), x, pol, false, true, &xterm) | |
59 | : detail::ibeta_imp(T(v / 2), T(k + 1), y, pol, true, true, &xterm); | |
60 | xterm *= y / (v / 2 + k); | |
61 | T poisf(pois), betaf(beta), xtermf(xterm); | |
62 | T sum = init_val; | |
63 | if((xterm == 0) && (beta == 0)) | |
64 | return init_val; | |
65 | ||
66 | // | |
67 | // Backwards recursion first, this is the stable | |
68 | // direction for recursion: | |
69 | // | |
70 | boost::uintmax_t count = 0; | |
71 | T last_term = 0; | |
72 | for(int i = k; i >= 0; --i) | |
73 | { | |
74 | T term = beta * pois; | |
75 | sum += term; | |
76 | // Don't terminate on first term in case we "fixed" k above: | |
77 | if((fabs(last_term) > fabs(term)) && fabs(term/sum) < errtol) | |
78 | break; | |
79 | last_term = term; | |
80 | pois *= (i + 0.5f) / d2; | |
81 | beta += xterm; | |
82 | xterm *= (i) / (x * (v / 2 + i - 1)); | |
83 | ++count; | |
84 | } | |
85 | last_term = 0; | |
86 | for(int i = k + 1; ; ++i) | |
87 | { | |
88 | poisf *= d2 / (i + 0.5f); | |
89 | xtermf *= (x * (v / 2 + i - 1)) / (i); | |
90 | betaf -= xtermf; | |
91 | T term = poisf * betaf; | |
92 | sum += term; | |
93 | if((fabs(last_term) >= fabs(term)) && (fabs(term/sum) < errtol)) | |
94 | break; | |
95 | last_term = term; | |
96 | ++count; | |
97 | if(count > max_iter) | |
98 | { | |
99 | return policies::raise_evaluation_error( | |
100 | "cdf(non_central_t_distribution<%1%>, %1%)", | |
101 | "Series did not converge, closest value was %1%", sum, pol); | |
102 | } | |
103 | } | |
104 | return sum; | |
105 | } | |
106 | ||
107 | template <class T, class Policy> | |
108 | T non_central_t2_q(T v, T delta, T x, T y, const Policy& pol, T init_val) | |
109 | { | |
110 | BOOST_MATH_STD_USING | |
111 | // | |
112 | // Variables come first: | |
113 | // | |
114 | boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); | |
115 | T errtol = boost::math::policies::get_epsilon<T, Policy>(); | |
116 | T d2 = delta * delta / 2; | |
117 | // | |
118 | // k is the starting point for iteration, and is the | |
119 | // maximum of the poisson weighting term, we don't allow | |
120 | // k == 0 as this can cause catastrophic cancellation errors | |
121 | // (test case is v = 561908036470413.25, delta = 0.056190803647041321, | |
122 | // x = 1.6155232703966216): | |
123 | // | |
124 | int k = itrunc(d2); | |
125 | if(k == 0) k = 1; | |
126 | // Starting Poisson weight: | |
127 | T pois; | |
128 | if((k < (int)(max_factorial<T>::value)) && (d2 < tools::log_max_value<T>()) && (log(d2) * k < tools::log_max_value<T>())) | |
129 | { | |
130 | // | |
131 | // For small k we can optimise this calculation by using | |
132 | // a simpler reduced formula: | |
133 | // | |
134 | pois = exp(-d2); | |
135 | pois *= pow(d2, static_cast<T>(k)); | |
136 | pois /= boost::math::tgamma(T(k + 1 + 0.5), pol); | |
137 | pois *= delta / constants::root_two<T>(); | |
138 | } | |
139 | else | |
140 | { | |
141 | pois = gamma_p_derivative(T(k+1), d2, pol) | |
142 | * tgamma_delta_ratio(T(k + 1), T(0.5f)) | |
143 | * delta / constants::root_two<T>(); | |
144 | } | |
145 | if(pois == 0) | |
146 | return init_val; | |
147 | // Recurance term: | |
148 | T xterm; | |
149 | T beta; | |
150 | // Starting beta term: | |
151 | if(k != 0) | |
152 | { | |
153 | beta = x < y | |
154 | ? detail::ibeta_imp(T(k + 1), T(v / 2), x, pol, true, true, &xterm) | |
155 | : detail::ibeta_imp(T(v / 2), T(k + 1), y, pol, false, true, &xterm); | |
156 | ||
157 | xterm *= y / (v / 2 + k); | |
158 | } | |
159 | else | |
160 | { | |
161 | beta = pow(y, v / 2); | |
162 | xterm = beta; | |
163 | } | |
164 | T poisf(pois), betaf(beta), xtermf(xterm); | |
165 | T sum = init_val; | |
166 | if((xterm == 0) && (beta == 0)) | |
167 | return init_val; | |
168 | ||
169 | // | |
170 | // Fused forward and backwards recursion: | |
171 | // | |
172 | boost::uintmax_t count = 0; | |
173 | T last_term = 0; | |
174 | for(int i = k + 1, j = k; ; ++i, --j) | |
175 | { | |
176 | poisf *= d2 / (i + 0.5f); | |
177 | xtermf *= (x * (v / 2 + i - 1)) / (i); | |
178 | betaf += xtermf; | |
179 | T term = poisf * betaf; | |
180 | ||
181 | if(j >= 0) | |
182 | { | |
183 | term += beta * pois; | |
184 | pois *= (j + 0.5f) / d2; | |
185 | beta -= xterm; | |
186 | xterm *= (j) / (x * (v / 2 + j - 1)); | |
187 | } | |
188 | ||
189 | sum += term; | |
190 | // Don't terminate on first term in case we "fixed" the value of k above: | |
191 | if((fabs(last_term) > fabs(term)) && fabs(term/sum) < errtol) | |
192 | break; | |
193 | last_term = term; | |
194 | if(count > max_iter) | |
195 | { | |
196 | return policies::raise_evaluation_error( | |
197 | "cdf(non_central_t_distribution<%1%>, %1%)", | |
198 | "Series did not converge, closest value was %1%", sum, pol); | |
199 | } | |
200 | ++count; | |
201 | } | |
202 | return sum; | |
203 | } | |
204 | ||
205 | template <class T, class Policy> | |
206 | T non_central_t_cdf(T v, T delta, T t, bool invert, const Policy& pol) | |
207 | { | |
208 | BOOST_MATH_STD_USING | |
209 | if ((boost::math::isinf)(v)) | |
210 | { // Infinite degrees of freedom, so use normal distribution located at delta. | |
211 | normal_distribution<T, Policy> n(delta, 1); | |
212 | return cdf(n, t); | |
213 | } | |
214 | // | |
215 | // Otherwise, for t < 0 we have to use the reflection formula: | |
216 | if(t < 0) | |
217 | { | |
218 | t = -t; | |
219 | delta = -delta; | |
220 | invert = !invert; | |
221 | } | |
222 | if(fabs(delta / (4 * v)) < policies::get_epsilon<T, Policy>()) | |
223 | { | |
224 | // Approximate with a Student's T centred on delta, | |
225 | // the crossover point is based on eq 2.6 from | |
226 | // "A Comparison of Approximations To Percentiles of the | |
227 | // Noncentral t-Distribution". H. Sahai and M. M. Ojeda, | |
228 | // Revista Investigacion Operacional Vol 21, No 2, 2000. | |
229 | // Original sources referenced in the above are: | |
230 | // "Some Approximations to the Percentage Points of the Noncentral | |
231 | // t-Distribution". C. van Eeden. International Statistical Review, 29, 4-31. | |
232 | // "Continuous Univariate Distributions". N.L. Johnson, S. Kotz and | |
233 | // N. Balkrishnan. 1995. John Wiley and Sons New York. | |
234 | T result = cdf(students_t_distribution<T, Policy>(v), t - delta); | |
235 | return invert ? 1 - result : result; | |
236 | } | |
237 | // | |
238 | // x and y are the corresponding random | |
239 | // variables for the noncentral beta distribution, | |
240 | // with y = 1 - x: | |
241 | // | |
242 | T x = t * t / (v + t * t); | |
243 | T y = v / (v + t * t); | |
244 | T d2 = delta * delta; | |
245 | T a = 0.5f; | |
246 | T b = v / 2; | |
247 | T c = a + b + d2 / 2; | |
248 | // | |
249 | // Crossover point for calculating p or q is the same | |
250 | // as for the noncentral beta: | |
251 | // | |
252 | T cross = 1 - (b / c) * (1 + d2 / (2 * c * c)); | |
253 | T result; | |
254 | if(x < cross) | |
255 | { | |
256 | // | |
257 | // Calculate p: | |
258 | // | |
259 | if(x != 0) | |
260 | { | |
261 | result = non_central_beta_p(a, b, d2, x, y, pol); | |
262 | result = non_central_t2_p(v, delta, x, y, pol, result); | |
263 | result /= 2; | |
264 | } | |
265 | else | |
266 | result = 0; | |
267 | result += cdf(boost::math::normal_distribution<T, Policy>(), -delta); | |
268 | } | |
269 | else | |
270 | { | |
271 | // | |
272 | // Calculate q: | |
273 | // | |
274 | invert = !invert; | |
275 | if(x != 0) | |
276 | { | |
277 | result = non_central_beta_q(a, b, d2, x, y, pol); | |
278 | result = non_central_t2_q(v, delta, x, y, pol, result); | |
279 | result /= 2; | |
280 | } | |
281 | else // x == 0 | |
282 | result = cdf(complement(boost::math::normal_distribution<T, Policy>(), -delta)); | |
283 | } | |
284 | if(invert) | |
285 | result = 1 - result; | |
286 | return result; | |
287 | } | |
288 | ||
289 | template <class T, class Policy> | |
290 | T non_central_t_quantile(const char* function, T v, T delta, T p, T q, const Policy&) | |
291 | { | |
292 | BOOST_MATH_STD_USING | |
293 | // static const char* function = "quantile(non_central_t_distribution<%1%>, %1%)"; | |
294 | // now passed as function | |
295 | typedef typename policies::evaluation<T, Policy>::type value_type; | |
296 | typedef typename policies::normalise< | |
297 | Policy, | |
298 | policies::promote_float<false>, | |
299 | policies::promote_double<false>, | |
300 | policies::discrete_quantile<>, | |
301 | policies::assert_undefined<> >::type forwarding_policy; | |
302 | ||
303 | T r; | |
304 | if(!detail::check_df_gt0_to_inf( | |
305 | function, | |
306 | v, &r, Policy()) | |
307 | || | |
308 | !detail::check_finite( | |
309 | function, | |
310 | delta, | |
311 | &r, | |
312 | Policy()) | |
313 | || | |
314 | !detail::check_probability( | |
315 | function, | |
316 | p, | |
317 | &r, | |
318 | Policy())) | |
319 | return r; | |
320 | ||
321 | ||
322 | value_type guess = 0; | |
323 | if ( ((boost::math::isinf)(v)) || (v > 1 / boost::math::tools::epsilon<T>()) ) | |
324 | { // Infinite or very large degrees of freedom, so use normal distribution located at delta. | |
325 | normal_distribution<T, Policy> n(delta, 1); | |
326 | if (p < q) | |
327 | { | |
328 | return quantile(n, p); | |
329 | } | |
330 | else | |
331 | { | |
332 | return quantile(complement(n, q)); | |
333 | } | |
334 | } | |
335 | else if(v > 3) | |
336 | { // Use normal distribution to calculate guess. | |
337 | value_type mean = (v > 1 / policies::get_epsilon<T, Policy>()) ? delta : delta * sqrt(v / 2) * tgamma_delta_ratio((v - 1) * 0.5f, T(0.5f)); | |
338 | value_type var = (v > 1 / policies::get_epsilon<T, Policy>()) ? value_type(1) : (((delta * delta + 1) * v) / (v - 2) - mean * mean); | |
339 | if(p < q) | |
340 | guess = quantile(normal_distribution<value_type, forwarding_policy>(mean, var), p); | |
341 | else | |
342 | guess = quantile(complement(normal_distribution<value_type, forwarding_policy>(mean, var), q)); | |
343 | } | |
344 | // | |
345 | // We *must* get the sign of the initial guess correct, | |
346 | // or our root-finder will fail, so double check it now: | |
347 | // | |
348 | value_type pzero = non_central_t_cdf( | |
349 | static_cast<value_type>(v), | |
350 | static_cast<value_type>(delta), | |
351 | static_cast<value_type>(0), | |
352 | !(p < q), | |
353 | forwarding_policy()); | |
354 | int s; | |
355 | if(p < q) | |
356 | s = boost::math::sign(p - pzero); | |
357 | else | |
358 | s = boost::math::sign(pzero - q); | |
359 | if(s != boost::math::sign(guess)) | |
360 | { | |
361 | guess = static_cast<T>(s); | |
362 | } | |
363 | ||
364 | value_type result = detail::generic_quantile( | |
365 | non_central_t_distribution<value_type, forwarding_policy>(v, delta), | |
366 | (p < q ? p : q), | |
367 | guess, | |
368 | (p >= q), | |
369 | function); | |
370 | return policies::checked_narrowing_cast<T, forwarding_policy>( | |
371 | result, | |
372 | function); | |
373 | } | |
374 | ||
375 | template <class T, class Policy> | |
376 | T non_central_t2_pdf(T n, T delta, T x, T y, const Policy& pol, T init_val) | |
377 | { | |
378 | BOOST_MATH_STD_USING | |
379 | // | |
380 | // Variables come first: | |
381 | // | |
382 | boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); | |
383 | T errtol = boost::math::policies::get_epsilon<T, Policy>(); | |
384 | T d2 = delta * delta / 2; | |
385 | // | |
386 | // k is the starting point for iteration, and is the | |
387 | // maximum of the poisson weighting term: | |
388 | // | |
389 | int k = itrunc(d2); | |
390 | T pois, xterm; | |
391 | if(k == 0) | |
392 | k = 1; | |
393 | // Starting Poisson weight: | |
394 | pois = gamma_p_derivative(T(k+1), d2, pol) | |
395 | * tgamma_delta_ratio(T(k + 1), T(0.5f)) | |
396 | * delta / constants::root_two<T>(); | |
397 | // Starting beta term: | |
398 | xterm = x < y | |
399 | ? ibeta_derivative(T(k + 1), n / 2, x, pol) | |
400 | : ibeta_derivative(n / 2, T(k + 1), y, pol); | |
401 | T poisf(pois), xtermf(xterm); | |
402 | T sum = init_val; | |
403 | if((pois == 0) || (xterm == 0)) | |
404 | return init_val; | |
405 | ||
406 | // | |
407 | // Backwards recursion first, this is the stable | |
408 | // direction for recursion: | |
409 | // | |
410 | boost::uintmax_t count = 0; | |
411 | for(int i = k; i >= 0; --i) | |
412 | { | |
413 | T term = xterm * pois; | |
414 | sum += term; | |
415 | if(((fabs(term/sum) < errtol) && (i != k)) || (term == 0)) | |
416 | break; | |
417 | pois *= (i + 0.5f) / d2; | |
418 | xterm *= (i) / (x * (n / 2 + i)); | |
419 | ++count; | |
420 | if(count > max_iter) | |
421 | { | |
422 | return policies::raise_evaluation_error( | |
423 | "pdf(non_central_t_distribution<%1%>, %1%)", | |
424 | "Series did not converge, closest value was %1%", sum, pol); | |
425 | } | |
426 | } | |
427 | for(int i = k + 1; ; ++i) | |
428 | { | |
429 | poisf *= d2 / (i + 0.5f); | |
430 | xtermf *= (x * (n / 2 + i)) / (i); | |
431 | T term = poisf * xtermf; | |
432 | sum += term; | |
433 | if((fabs(term/sum) < errtol) || (term == 0)) | |
434 | break; | |
435 | ++count; | |
436 | if(count > max_iter) | |
437 | { | |
438 | return policies::raise_evaluation_error( | |
439 | "pdf(non_central_t_distribution<%1%>, %1%)", | |
440 | "Series did not converge, closest value was %1%", sum, pol); | |
441 | } | |
442 | } | |
443 | return sum; | |
444 | } | |
445 | ||
446 | template <class T, class Policy> | |
447 | T non_central_t_pdf(T n, T delta, T t, const Policy& pol) | |
448 | { | |
449 | BOOST_MATH_STD_USING | |
450 | if ((boost::math::isinf)(n)) | |
451 | { // Infinite degrees of freedom, so use normal distribution located at delta. | |
452 | normal_distribution<T, Policy> norm(delta, 1); | |
453 | return pdf(norm, t); | |
454 | } | |
455 | // | |
456 | // Otherwise, for t < 0 we have to use the reflection formula: | |
457 | if(t < 0) | |
458 | { | |
459 | t = -t; | |
460 | delta = -delta; | |
461 | } | |
462 | if(t == 0) | |
463 | { | |
464 | // | |
465 | // Handle this as a special case, using the formula | |
466 | // from Weisstein, Eric W. | |
467 | // "Noncentral Student's t-Distribution." | |
468 | // From MathWorld--A Wolfram Web Resource. | |
469 | // http://mathworld.wolfram.com/NoncentralStudentst-Distribution.html | |
470 | // | |
471 | // The formula is simplified thanks to the relation | |
472 | // 1F1(a,b,0) = 1. | |
473 | // | |
474 | return tgamma_delta_ratio(n / 2 + 0.5f, T(0.5f)) | |
475 | * sqrt(n / constants::pi<T>()) | |
476 | * exp(-delta * delta / 2) / 2; | |
477 | } | |
478 | if(fabs(delta / (4 * n)) < policies::get_epsilon<T, Policy>()) | |
479 | { | |
480 | // Approximate with a Student's T centred on delta, | |
481 | // the crossover point is based on eq 2.6 from | |
482 | // "A Comparison of Approximations To Percentiles of the | |
483 | // Noncentral t-Distribution". H. Sahai and M. M. Ojeda, | |
484 | // Revista Investigacion Operacional Vol 21, No 2, 2000. | |
485 | // Original sources referenced in the above are: | |
486 | // "Some Approximations to the Percentage Points of the Noncentral | |
487 | // t-Distribution". C. van Eeden. International Statistical Review, 29, 4-31. | |
488 | // "Continuous Univariate Distributions". N.L. Johnson, S. Kotz and | |
489 | // N. Balkrishnan. 1995. John Wiley and Sons New York. | |
490 | return pdf(students_t_distribution<T, Policy>(n), t - delta); | |
491 | } | |
492 | // | |
493 | // x and y are the corresponding random | |
494 | // variables for the noncentral beta distribution, | |
495 | // with y = 1 - x: | |
496 | // | |
497 | T x = t * t / (n + t * t); | |
498 | T y = n / (n + t * t); | |
499 | T a = 0.5f; | |
500 | T b = n / 2; | |
501 | T d2 = delta * delta; | |
502 | // | |
503 | // Calculate pdf: | |
504 | // | |
505 | T dt = n * t / (n * n + 2 * n * t * t + t * t * t * t); | |
506 | T result = non_central_beta_pdf(a, b, d2, x, y, pol); | |
507 | T tol = tools::epsilon<T>() * result * 500; | |
508 | result = non_central_t2_pdf(n, delta, x, y, pol, result); | |
509 | if(result <= tol) | |
510 | result = 0; | |
511 | result *= dt; | |
512 | return result; | |
513 | } | |
514 | ||
515 | template <class T, class Policy> | |
516 | T mean(T v, T delta, const Policy& pol) | |
517 | { | |
518 | if ((boost::math::isinf)(v)) | |
519 | { | |
520 | return delta; | |
521 | } | |
522 | BOOST_MATH_STD_USING | |
523 | if (v > 1 / boost::math::tools::epsilon<T>() ) | |
524 | { | |
525 | //normal_distribution<T, Policy> n(delta, 1); | |
526 | //return boost::math::mean(n); | |
527 | return delta; | |
528 | } | |
529 | else | |
530 | { | |
531 | return delta * sqrt(v / 2) * tgamma_delta_ratio((v - 1) * 0.5f, T(0.5f), pol); | |
532 | } | |
533 | // Other moments use mean so using normal distribution is propagated. | |
534 | } | |
535 | ||
536 | template <class T, class Policy> | |
537 | T variance(T v, T delta, const Policy& pol) | |
538 | { | |
539 | if ((boost::math::isinf)(v)) | |
540 | { | |
541 | return 1; | |
542 | } | |
543 | if (delta == 0) | |
544 | { // == Student's t | |
545 | return v / (v - 2); | |
546 | } | |
547 | T result = ((delta * delta + 1) * v) / (v - 2); | |
548 | T m = mean(v, delta, pol); | |
549 | result -= m * m; | |
550 | return result; | |
551 | } | |
552 | ||
553 | template <class T, class Policy> | |
554 | T skewness(T v, T delta, const Policy& pol) | |
555 | { | |
556 | BOOST_MATH_STD_USING | |
557 | if ((boost::math::isinf)(v)) | |
558 | { | |
559 | return 0; | |
560 | } | |
561 | if(delta == 0) | |
562 | { // == Student's t | |
563 | return 0; | |
564 | } | |
565 | T mean = boost::math::detail::mean(v, delta, pol); | |
566 | T l2 = delta * delta; | |
567 | T var = ((l2 + 1) * v) / (v - 2) - mean * mean; | |
568 | T result = -2 * var; | |
569 | result += v * (l2 + 2 * v - 3) / ((v - 3) * (v - 2)); | |
570 | result *= mean; | |
571 | result /= pow(var, T(1.5f)); | |
572 | return result; | |
573 | } | |
574 | ||
575 | template <class T, class Policy> | |
576 | T kurtosis_excess(T v, T delta, const Policy& pol) | |
577 | { | |
578 | BOOST_MATH_STD_USING | |
579 | if ((boost::math::isinf)(v)) | |
580 | { | |
581 | return 3; | |
582 | } | |
583 | if (delta == 0) | |
584 | { // == Student's t | |
585 | return 3; | |
586 | } | |
587 | T mean = boost::math::detail::mean(v, delta, pol); | |
588 | T l2 = delta * delta; | |
589 | T var = ((l2 + 1) * v) / (v - 2) - mean * mean; | |
590 | T result = -3 * var; | |
591 | result += v * (l2 * (v + 1) + 3 * (3 * v - 5)) / ((v - 3) * (v - 2)); | |
592 | result *= -mean * mean; | |
593 | result += v * v * (l2 * l2 + 6 * l2 + 3) / ((v - 4) * (v - 2)); | |
594 | result /= var * var; | |
595 | return result; | |
596 | } | |
597 | ||
598 | #if 0 | |
599 | // | |
600 | // This code is disabled, since there can be multiple answers to the | |
601 | // question, and it's not clear how to find the "right" one. | |
602 | // | |
603 | template <class RealType, class Policy> | |
604 | struct t_degrees_of_freedom_finder | |
605 | { | |
606 | t_degrees_of_freedom_finder( | |
607 | RealType delta_, RealType x_, RealType p_, bool c) | |
608 | : delta(delta_), x(x_), p(p_), comp(c) {} | |
609 | ||
610 | RealType operator()(const RealType& v) | |
611 | { | |
612 | non_central_t_distribution<RealType, Policy> d(v, delta); | |
613 | return comp ? | |
614 | p - cdf(complement(d, x)) | |
615 | : cdf(d, x) - p; | |
616 | } | |
617 | private: | |
618 | RealType delta; | |
619 | RealType x; | |
620 | RealType p; | |
621 | bool comp; | |
622 | }; | |
623 | ||
624 | template <class RealType, class Policy> | |
625 | inline RealType find_t_degrees_of_freedom( | |
626 | RealType delta, RealType x, RealType p, RealType q, const Policy& pol) | |
627 | { | |
628 | const char* function = "non_central_t<%1%>::find_degrees_of_freedom"; | |
629 | if((p == 0) || (q == 0)) | |
630 | { | |
631 | // | |
632 | // Can't a thing if one of p and q is zero: | |
633 | // | |
634 | return policies::raise_evaluation_error<RealType>(function, | |
635 | "Can't find degrees of freedom when the probability is 0 or 1, only possible answer is %1%", | |
636 | RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy()); | |
637 | } | |
638 | t_degrees_of_freedom_finder<RealType, Policy> f(delta, x, p < q ? p : q, p < q ? false : true); | |
639 | tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>()); | |
640 | boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); | |
641 | // | |
642 | // Pick an initial guess: | |
643 | // | |
644 | RealType guess = 200; | |
645 | std::pair<RealType, RealType> ir = tools::bracket_and_solve_root( | |
646 | f, guess, RealType(2), false, tol, max_iter, pol); | |
647 | RealType result = ir.first + (ir.second - ir.first) / 2; | |
648 | if(max_iter >= policies::get_max_root_iterations<Policy>()) | |
649 | { | |
650 | return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:" | |
651 | " or there is no answer to problem. Current best guess is %1%", result, Policy()); | |
652 | } | |
653 | return result; | |
654 | } | |
655 | ||
656 | template <class RealType, class Policy> | |
657 | struct t_non_centrality_finder | |
658 | { | |
659 | t_non_centrality_finder( | |
660 | RealType v_, RealType x_, RealType p_, bool c) | |
661 | : v(v_), x(x_), p(p_), comp(c) {} | |
662 | ||
663 | RealType operator()(const RealType& delta) | |
664 | { | |
665 | non_central_t_distribution<RealType, Policy> d(v, delta); | |
666 | return comp ? | |
667 | p - cdf(complement(d, x)) | |
668 | : cdf(d, x) - p; | |
669 | } | |
670 | private: | |
671 | RealType v; | |
672 | RealType x; | |
673 | RealType p; | |
674 | bool comp; | |
675 | }; | |
676 | ||
677 | template <class RealType, class Policy> | |
678 | inline RealType find_t_non_centrality( | |
679 | RealType v, RealType x, RealType p, RealType q, const Policy& pol) | |
680 | { | |
681 | const char* function = "non_central_t<%1%>::find_t_non_centrality"; | |
682 | if((p == 0) || (q == 0)) | |
683 | { | |
684 | // | |
685 | // Can't do a thing if one of p and q is zero: | |
686 | // | |
687 | return policies::raise_evaluation_error<RealType>(function, | |
688 | "Can't find non-centrality parameter when the probability is 0 or 1, only possible answer is %1%", | |
689 | RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy()); | |
690 | } | |
691 | t_non_centrality_finder<RealType, Policy> f(v, x, p < q ? p : q, p < q ? false : true); | |
692 | tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>()); | |
693 | boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); | |
694 | // | |
695 | // Pick an initial guess that we know is the right side of | |
696 | // zero: | |
697 | // | |
698 | RealType guess; | |
699 | if(f(0) < 0) | |
700 | guess = 1; | |
701 | else | |
702 | guess = -1; | |
703 | std::pair<RealType, RealType> ir = tools::bracket_and_solve_root( | |
704 | f, guess, RealType(2), false, tol, max_iter, pol); | |
705 | RealType result = ir.first + (ir.second - ir.first) / 2; | |
706 | if(max_iter >= policies::get_max_root_iterations<Policy>()) | |
707 | { | |
708 | return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:" | |
709 | " or there is no answer to problem. Current best guess is %1%", result, Policy()); | |
710 | } | |
711 | return result; | |
712 | } | |
713 | #endif | |
714 | } // namespace detail ====================================================================== | |
715 | ||
716 | template <class RealType = double, class Policy = policies::policy<> > | |
717 | class non_central_t_distribution | |
718 | { | |
719 | public: | |
720 | typedef RealType value_type; | |
721 | typedef Policy policy_type; | |
722 | ||
723 | non_central_t_distribution(RealType v_, RealType lambda) : v(v_), ncp(lambda) | |
724 | { | |
725 | const char* function = "boost::math::non_central_t_distribution<%1%>::non_central_t_distribution(%1%,%1%)"; | |
726 | RealType r; | |
727 | detail::check_df_gt0_to_inf( | |
728 | function, | |
729 | v, &r, Policy()); | |
730 | detail::check_finite( | |
731 | function, | |
732 | lambda, | |
733 | &r, | |
734 | Policy()); | |
735 | } // non_central_t_distribution constructor. | |
736 | ||
737 | RealType degrees_of_freedom() const | |
738 | { // Private data getter function. | |
739 | return v; | |
740 | } | |
741 | RealType non_centrality() const | |
742 | { // Private data getter function. | |
743 | return ncp; | |
744 | } | |
745 | #if 0 | |
746 | // | |
747 | // This code is disabled, since there can be multiple answers to the | |
748 | // question, and it's not clear how to find the "right" one. | |
749 | // | |
750 | static RealType find_degrees_of_freedom(RealType delta, RealType x, RealType p) | |
751 | { | |
752 | const char* function = "non_central_t<%1%>::find_degrees_of_freedom"; | |
753 | typedef typename policies::evaluation<RealType, Policy>::type value_type; | |
754 | typedef typename policies::normalise< | |
755 | Policy, | |
756 | policies::promote_float<false>, | |
757 | policies::promote_double<false>, | |
758 | policies::discrete_quantile<>, | |
759 | policies::assert_undefined<> >::type forwarding_policy; | |
760 | value_type result = detail::find_t_degrees_of_freedom( | |
761 | static_cast<value_type>(delta), | |
762 | static_cast<value_type>(x), | |
763 | static_cast<value_type>(p), | |
764 | static_cast<value_type>(1-p), | |
765 | forwarding_policy()); | |
766 | return policies::checked_narrowing_cast<RealType, forwarding_policy>( | |
767 | result, | |
768 | function); | |
769 | } | |
770 | template <class A, class B, class C> | |
771 | static RealType find_degrees_of_freedom(const complemented3_type<A,B,C>& c) | |
772 | { | |
773 | const char* function = "non_central_t<%1%>::find_degrees_of_freedom"; | |
774 | typedef typename policies::evaluation<RealType, Policy>::type value_type; | |
775 | typedef typename policies::normalise< | |
776 | Policy, | |
777 | policies::promote_float<false>, | |
778 | policies::promote_double<false>, | |
779 | policies::discrete_quantile<>, | |
780 | policies::assert_undefined<> >::type forwarding_policy; | |
781 | value_type result = detail::find_t_degrees_of_freedom( | |
782 | static_cast<value_type>(c.dist), | |
783 | static_cast<value_type>(c.param1), | |
784 | static_cast<value_type>(1-c.param2), | |
785 | static_cast<value_type>(c.param2), | |
786 | forwarding_policy()); | |
787 | return policies::checked_narrowing_cast<RealType, forwarding_policy>( | |
788 | result, | |
789 | function); | |
790 | } | |
791 | static RealType find_non_centrality(RealType v, RealType x, RealType p) | |
792 | { | |
793 | const char* function = "non_central_t<%1%>::find_t_non_centrality"; | |
794 | typedef typename policies::evaluation<RealType, Policy>::type value_type; | |
795 | typedef typename policies::normalise< | |
796 | Policy, | |
797 | policies::promote_float<false>, | |
798 | policies::promote_double<false>, | |
799 | policies::discrete_quantile<>, | |
800 | policies::assert_undefined<> >::type forwarding_policy; | |
801 | value_type result = detail::find_t_non_centrality( | |
802 | static_cast<value_type>(v), | |
803 | static_cast<value_type>(x), | |
804 | static_cast<value_type>(p), | |
805 | static_cast<value_type>(1-p), | |
806 | forwarding_policy()); | |
807 | return policies::checked_narrowing_cast<RealType, forwarding_policy>( | |
808 | result, | |
809 | function); | |
810 | } | |
811 | template <class A, class B, class C> | |
812 | static RealType find_non_centrality(const complemented3_type<A,B,C>& c) | |
813 | { | |
814 | const char* function = "non_central_t<%1%>::find_t_non_centrality"; | |
815 | typedef typename policies::evaluation<RealType, Policy>::type value_type; | |
816 | typedef typename policies::normalise< | |
817 | Policy, | |
818 | policies::promote_float<false>, | |
819 | policies::promote_double<false>, | |
820 | policies::discrete_quantile<>, | |
821 | policies::assert_undefined<> >::type forwarding_policy; | |
822 | value_type result = detail::find_t_non_centrality( | |
823 | static_cast<value_type>(c.dist), | |
824 | static_cast<value_type>(c.param1), | |
825 | static_cast<value_type>(1-c.param2), | |
826 | static_cast<value_type>(c.param2), | |
827 | forwarding_policy()); | |
828 | return policies::checked_narrowing_cast<RealType, forwarding_policy>( | |
829 | result, | |
830 | function); | |
831 | } | |
832 | #endif | |
833 | private: | |
834 | // Data member, initialized by constructor. | |
835 | RealType v; // degrees of freedom | |
836 | RealType ncp; // non-centrality parameter | |
837 | }; // template <class RealType, class Policy> class non_central_t_distribution | |
838 | ||
839 | typedef non_central_t_distribution<double> non_central_t; // Reserved name of type double. | |
840 | ||
841 | // Non-member functions to give properties of the distribution. | |
842 | ||
843 | template <class RealType, class Policy> | |
844 | inline const std::pair<RealType, RealType> range(const non_central_t_distribution<RealType, Policy>& /* dist */) | |
845 | { // Range of permissible values for random variable k. | |
846 | using boost::math::tools::max_value; | |
847 | return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>()); | |
848 | } | |
849 | ||
850 | template <class RealType, class Policy> | |
851 | inline const std::pair<RealType, RealType> support(const non_central_t_distribution<RealType, Policy>& /* dist */) | |
852 | { // Range of supported values for random variable k. | |
853 | // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. | |
854 | using boost::math::tools::max_value; | |
855 | return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>()); | |
856 | } | |
857 | ||
858 | template <class RealType, class Policy> | |
859 | inline RealType mode(const non_central_t_distribution<RealType, Policy>& dist) | |
860 | { // mode. | |
861 | static const char* function = "mode(non_central_t_distribution<%1%> const&)"; | |
862 | RealType v = dist.degrees_of_freedom(); | |
863 | RealType l = dist.non_centrality(); | |
864 | RealType r; | |
865 | if(!detail::check_df_gt0_to_inf( | |
866 | function, | |
867 | v, &r, Policy()) | |
868 | || | |
869 | !detail::check_finite( | |
870 | function, | |
871 | l, | |
872 | &r, | |
873 | Policy())) | |
874 | return (RealType)r; | |
875 | ||
876 | BOOST_MATH_STD_USING | |
877 | ||
878 | RealType m = v < 3 ? 0 : detail::mean(v, l, Policy()); | |
879 | RealType var = v < 4 ? 1 : detail::variance(v, l, Policy()); | |
880 | ||
881 | return detail::generic_find_mode( | |
882 | dist, | |
883 | m, | |
884 | function, | |
885 | sqrt(var)); | |
886 | } | |
887 | ||
888 | template <class RealType, class Policy> | |
889 | inline RealType mean(const non_central_t_distribution<RealType, Policy>& dist) | |
890 | { | |
891 | BOOST_MATH_STD_USING | |
892 | const char* function = "mean(const non_central_t_distribution<%1%>&)"; | |
893 | typedef typename policies::evaluation<RealType, Policy>::type value_type; | |
894 | typedef typename policies::normalise< | |
895 | Policy, | |
896 | policies::promote_float<false>, | |
897 | policies::promote_double<false>, | |
898 | policies::discrete_quantile<>, | |
899 | policies::assert_undefined<> >::type forwarding_policy; | |
900 | RealType v = dist.degrees_of_freedom(); | |
901 | RealType l = dist.non_centrality(); | |
902 | RealType r; | |
903 | if(!detail::check_df_gt0_to_inf( | |
904 | function, | |
905 | v, &r, Policy()) | |
906 | || | |
907 | !detail::check_finite( | |
908 | function, | |
909 | l, | |
910 | &r, | |
911 | Policy())) | |
912 | return (RealType)r; | |
913 | if(v <= 1) | |
914 | return policies::raise_domain_error<RealType>( | |
915 | function, | |
916 | "The non-central t distribution has no defined mean for degrees of freedom <= 1: got v=%1%.", v, Policy()); | |
917 | // return l * sqrt(v / 2) * tgamma_delta_ratio((v - 1) * 0.5f, RealType(0.5f)); | |
918 | return policies::checked_narrowing_cast<RealType, forwarding_policy>( | |
919 | detail::mean(static_cast<value_type>(v), static_cast<value_type>(l), forwarding_policy()), function); | |
920 | ||
921 | } // mean | |
922 | ||
923 | template <class RealType, class Policy> | |
924 | inline RealType variance(const non_central_t_distribution<RealType, Policy>& dist) | |
925 | { // variance. | |
926 | const char* function = "variance(const non_central_t_distribution<%1%>&)"; | |
927 | typedef typename policies::evaluation<RealType, Policy>::type value_type; | |
928 | typedef typename policies::normalise< | |
929 | Policy, | |
930 | policies::promote_float<false>, | |
931 | policies::promote_double<false>, | |
932 | policies::discrete_quantile<>, | |
933 | policies::assert_undefined<> >::type forwarding_policy; | |
934 | BOOST_MATH_STD_USING | |
935 | RealType v = dist.degrees_of_freedom(); | |
936 | RealType l = dist.non_centrality(); | |
937 | RealType r; | |
938 | if(!detail::check_df_gt0_to_inf( | |
939 | function, | |
940 | v, &r, Policy()) | |
941 | || | |
942 | !detail::check_finite( | |
943 | function, | |
944 | l, | |
945 | &r, | |
946 | Policy())) | |
947 | return (RealType)r; | |
948 | if(v <= 2) | |
949 | return policies::raise_domain_error<RealType>( | |
950 | function, | |
951 | "The non-central t distribution has no defined variance for degrees of freedom <= 2: got v=%1%.", v, Policy()); | |
952 | return policies::checked_narrowing_cast<RealType, forwarding_policy>( | |
953 | detail::variance(static_cast<value_type>(v), static_cast<value_type>(l), forwarding_policy()), function); | |
954 | } | |
955 | ||
956 | // RealType standard_deviation(const non_central_t_distribution<RealType, Policy>& dist) | |
957 | // standard_deviation provided by derived accessors. | |
958 | ||
959 | template <class RealType, class Policy> | |
960 | inline RealType skewness(const non_central_t_distribution<RealType, Policy>& dist) | |
961 | { // skewness = sqrt(l). | |
962 | const char* function = "skewness(const non_central_t_distribution<%1%>&)"; | |
963 | typedef typename policies::evaluation<RealType, Policy>::type value_type; | |
964 | typedef typename policies::normalise< | |
965 | Policy, | |
966 | policies::promote_float<false>, | |
967 | policies::promote_double<false>, | |
968 | policies::discrete_quantile<>, | |
969 | policies::assert_undefined<> >::type forwarding_policy; | |
970 | RealType v = dist.degrees_of_freedom(); | |
971 | RealType l = dist.non_centrality(); | |
972 | RealType r; | |
973 | if(!detail::check_df_gt0_to_inf( | |
974 | function, | |
975 | v, &r, Policy()) | |
976 | || | |
977 | !detail::check_finite( | |
978 | function, | |
979 | l, | |
980 | &r, | |
981 | Policy())) | |
982 | return (RealType)r; | |
983 | if(v <= 3) | |
984 | return policies::raise_domain_error<RealType>( | |
985 | function, | |
986 | "The non-central t distribution has no defined skewness for degrees of freedom <= 3: got v=%1%.", v, Policy());; | |
987 | return policies::checked_narrowing_cast<RealType, forwarding_policy>( | |
988 | detail::skewness(static_cast<value_type>(v), static_cast<value_type>(l), forwarding_policy()), function); | |
989 | } | |
990 | ||
991 | template <class RealType, class Policy> | |
992 | inline RealType kurtosis_excess(const non_central_t_distribution<RealType, Policy>& dist) | |
993 | { | |
994 | const char* function = "kurtosis_excess(const non_central_t_distribution<%1%>&)"; | |
995 | typedef typename policies::evaluation<RealType, Policy>::type value_type; | |
996 | typedef typename policies::normalise< | |
997 | Policy, | |
998 | policies::promote_float<false>, | |
999 | policies::promote_double<false>, | |
1000 | policies::discrete_quantile<>, | |
1001 | policies::assert_undefined<> >::type forwarding_policy; | |
1002 | RealType v = dist.degrees_of_freedom(); | |
1003 | RealType l = dist.non_centrality(); | |
1004 | RealType r; | |
1005 | if(!detail::check_df_gt0_to_inf( | |
1006 | function, | |
1007 | v, &r, Policy()) | |
1008 | || | |
1009 | !detail::check_finite( | |
1010 | function, | |
1011 | l, | |
1012 | &r, | |
1013 | Policy())) | |
1014 | return (RealType)r; | |
1015 | if(v <= 4) | |
1016 | return policies::raise_domain_error<RealType>( | |
1017 | function, | |
1018 | "The non-central t distribution has no defined kurtosis for degrees of freedom <= 4: got v=%1%.", v, Policy());; | |
1019 | return policies::checked_narrowing_cast<RealType, forwarding_policy>( | |
1020 | detail::kurtosis_excess(static_cast<value_type>(v), static_cast<value_type>(l), forwarding_policy()), function); | |
1021 | } // kurtosis_excess | |
1022 | ||
1023 | template <class RealType, class Policy> | |
1024 | inline RealType kurtosis(const non_central_t_distribution<RealType, Policy>& dist) | |
1025 | { | |
1026 | return kurtosis_excess(dist) + 3; | |
1027 | } | |
1028 | ||
1029 | template <class RealType, class Policy> | |
1030 | inline RealType pdf(const non_central_t_distribution<RealType, Policy>& dist, const RealType& t) | |
1031 | { // Probability Density/Mass Function. | |
1032 | const char* function = "pdf(non_central_t_distribution<%1%>, %1%)"; | |
1033 | typedef typename policies::evaluation<RealType, Policy>::type value_type; | |
1034 | typedef typename policies::normalise< | |
1035 | Policy, | |
1036 | policies::promote_float<false>, | |
1037 | policies::promote_double<false>, | |
1038 | policies::discrete_quantile<>, | |
1039 | policies::assert_undefined<> >::type forwarding_policy; | |
1040 | ||
1041 | RealType v = dist.degrees_of_freedom(); | |
1042 | RealType l = dist.non_centrality(); | |
1043 | RealType r; | |
1044 | if(!detail::check_df_gt0_to_inf( | |
1045 | function, | |
1046 | v, &r, Policy()) | |
1047 | || | |
1048 | !detail::check_finite( | |
1049 | function, | |
1050 | l, | |
1051 | &r, | |
1052 | Policy()) | |
1053 | || | |
1054 | !detail::check_x( | |
1055 | function, | |
1056 | t, | |
1057 | &r, | |
1058 | Policy())) | |
1059 | return (RealType)r; | |
1060 | return policies::checked_narrowing_cast<RealType, forwarding_policy>( | |
1061 | detail::non_central_t_pdf(static_cast<value_type>(v), | |
1062 | static_cast<value_type>(l), | |
1063 | static_cast<value_type>(t), | |
1064 | Policy()), | |
1065 | function); | |
1066 | ||
1067 | ||
1068 | template <class RealType, class Policy> | |
1069 | RealType cdf(const non_central_t_distribution<RealType, Policy>& dist, const RealType& x) | |
1070 | { | |
1071 | const char* function = "boost::math::cdf(non_central_t_distribution<%1%>&, %1%)"; | |
1072 | // was const char* function = "boost::math::non_central_t_distribution<%1%>::cdf(%1%)"; | |
1073 | typedef typename policies::evaluation<RealType, Policy>::type value_type; | |
1074 | typedef typename policies::normalise< | |
1075 | Policy, | |
1076 | policies::promote_float<false>, | |
1077 | policies::promote_double<false>, | |
1078 | policies::discrete_quantile<>, | |
1079 | policies::assert_undefined<> >::type forwarding_policy; | |
1080 | ||
1081 | RealType v = dist.degrees_of_freedom(); | |
1082 | RealType l = dist.non_centrality(); | |
1083 | RealType r; | |
1084 | if(!detail::check_df_gt0_to_inf( | |
1085 | function, | |
1086 | v, &r, Policy()) | |
1087 | || | |
1088 | !detail::check_finite( | |
1089 | function, | |
1090 | l, | |
1091 | &r, | |
1092 | Policy()) | |
1093 | || | |
1094 | !detail::check_x( | |
1095 | function, | |
1096 | x, | |
1097 | &r, | |
1098 | Policy())) | |
1099 | return (RealType)r; | |
1100 | if ((boost::math::isinf)(v)) | |
1101 | { // Infinite degrees of freedom, so use normal distribution located at delta. | |
1102 | normal_distribution<RealType, Policy> n(l, 1); | |
1103 | cdf(n, x); | |
1104 | //return cdf(normal_distribution<RealType, Policy>(l, 1), x); | |
1105 | } | |
1106 | ||
1107 | if(l == 0) | |
1108 | { // NO non-centrality, so use Student's t instead. | |
1109 | return cdf(students_t_distribution<RealType, Policy>(v), x); | |
1110 | } | |
1111 | return policies::checked_narrowing_cast<RealType, forwarding_policy>( | |
1112 | detail::non_central_t_cdf( | |
1113 | static_cast<value_type>(v), | |
1114 | static_cast<value_type>(l), | |
1115 | static_cast<value_type>(x), | |
1116 | false, Policy()), | |
1117 | function); | |
1118 | } // cdf | |
1119 | ||
1120 | template <class RealType, class Policy> | |
1121 | RealType cdf(const complemented2_type<non_central_t_distribution<RealType, Policy>, RealType>& c) | |
1122 | { // Complemented Cumulative Distribution Function | |
1123 | // was const char* function = "boost::math::non_central_t_distribution<%1%>::cdf(%1%)"; | |
1124 | const char* function = "boost::math::cdf(const complement(non_central_t_distribution<%1%>&), %1%)"; | |
1125 | typedef typename policies::evaluation<RealType, Policy>::type value_type; | |
1126 | typedef typename policies::normalise< | |
1127 | Policy, | |
1128 | policies::promote_float<false>, | |
1129 | policies::promote_double<false>, | |
1130 | policies::discrete_quantile<>, | |
1131 | policies::assert_undefined<> >::type forwarding_policy; | |
1132 | ||
1133 | non_central_t_distribution<RealType, Policy> const& dist = c.dist; | |
1134 | RealType x = c.param; | |
1135 | RealType v = dist.degrees_of_freedom(); | |
1136 | RealType l = dist.non_centrality(); // aka delta | |
1137 | RealType r; | |
1138 | if(!detail::check_df_gt0_to_inf( | |
1139 | function, | |
1140 | v, &r, Policy()) | |
1141 | || | |
1142 | !detail::check_finite( | |
1143 | function, | |
1144 | l, | |
1145 | &r, | |
1146 | Policy()) | |
1147 | || | |
1148 | !detail::check_x( | |
1149 | function, | |
1150 | x, | |
1151 | &r, | |
1152 | Policy())) | |
1153 | return (RealType)r; | |
1154 | ||
1155 | if ((boost::math::isinf)(v)) | |
1156 | { // Infinite degrees of freedom, so use normal distribution located at delta. | |
1157 | normal_distribution<RealType, Policy> n(l, 1); | |
1158 | return cdf(complement(n, x)); | |
1159 | } | |
1160 | if(l == 0) | |
1161 | { // zero non-centrality so use Student's t distribution. | |
1162 | return cdf(complement(students_t_distribution<RealType, Policy>(v), x)); | |
1163 | } | |
1164 | return policies::checked_narrowing_cast<RealType, forwarding_policy>( | |
1165 | detail::non_central_t_cdf( | |
1166 | static_cast<value_type>(v), | |
1167 | static_cast<value_type>(l), | |
1168 | static_cast<value_type>(x), | |
1169 | true, Policy()), | |
1170 | function); | |
1171 | } // ccdf | |
1172 | ||
1173 | template <class RealType, class Policy> | |
1174 | inline RealType quantile(const non_central_t_distribution<RealType, Policy>& dist, const RealType& p) | |
1175 | { // Quantile (or Percent Point) function. | |
1176 | static const char* function = "quantile(const non_central_t_distribution<%1%>, %1%)"; | |
1177 | RealType v = dist.degrees_of_freedom(); | |
1178 | RealType l = dist.non_centrality(); | |
1179 | return detail::non_central_t_quantile(function, v, l, p, RealType(1-p), Policy()); | |
1180 | } // quantile | |
1181 | ||
1182 | template <class RealType, class Policy> | |
1183 | inline RealType quantile(const complemented2_type<non_central_t_distribution<RealType, Policy>, RealType>& c) | |
1184 | { // Quantile (or Percent Point) function. | |
1185 | static const char* function = "quantile(const complement(non_central_t_distribution<%1%>, %1%))"; | |
1186 | non_central_t_distribution<RealType, Policy> const& dist = c.dist; | |
1187 | RealType q = c.param; | |
1188 | RealType v = dist.degrees_of_freedom(); | |
1189 | RealType l = dist.non_centrality(); | |
1190 | return detail::non_central_t_quantile(function, v, l, RealType(1-q), q, Policy()); | |
1191 | } // quantile complement. | |
1192 | ||
1193 | } // namespace math | |
1194 | } // namespace boost | |
1195 | ||
1196 | // This include must be at the end, *after* the accessors | |
1197 | // for this distribution have been defined, in order to | |
1198 | // keep compilers that support two-phase lookup happy. | |
1199 | #include <boost/math/distributions/detail/derived_accessors.hpp> | |
1200 | ||
1201 | #endif // BOOST_MATH_SPECIAL_NON_CENTRAL_T_HPP | |
1202 |