1 // boost\math\distributions\non_central_t.hpp
3 // Copyright John Maddock 2008.
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)
10 #ifndef BOOST_MATH_SPECIAL_NON_CENTRAL_T_HPP
11 #define BOOST_MATH_SPECIAL_NON_CENTRAL_T_HPP
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
24 template <class RealType, class Policy>
25 class non_central_t_distribution;
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)
34 // Variables come first:
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;
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).
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>();
56 // Recurrance & starting beta terms:
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);
63 if((xterm == 0) && (beta == 0))
67 // Backwards recursion first, this is the stable
68 // direction for recursion:
70 boost::uintmax_t count = 0;
72 for(int i = k; i >= 0; --i)
76 // Don't terminate on first term in case we "fixed" k above:
77 if((fabs(last_term) > fabs(term)) && fabs(term/sum) < errtol)
80 pois *= (i + 0.5f) / d2;
82 xterm *= (i) / (x * (v / 2 + i - 1));
86 for(int i = k + 1; ; ++i)
88 poisf *= d2 / (i + 0.5f);
89 xtermf *= (x * (v / 2 + i - 1)) / (i);
91 T term = poisf * betaf;
93 if((fabs(last_term) >= fabs(term)) && (fabs(term/sum) < errtol))
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);
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)
112 // Variables come first:
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;
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):
126 // Starting Poisson weight:
128 if((k < (int)(max_factorial<T>::value)) && (d2 < tools::log_max_value<T>()) && (log(d2) * k < tools::log_max_value<T>()))
131 // For small k we can optimise this calculation by using
132 // a simpler reduced formula:
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>();
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>();
150 // Starting beta term:
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);
157 xterm *= y / (v / 2 + k);
161 beta = pow(y, v / 2);
164 T poisf(pois), betaf(beta), xtermf(xterm);
166 if((xterm == 0) && (beta == 0))
170 // Fused forward and backwards recursion:
172 boost::uintmax_t count = 0;
174 for(int i = k + 1, j = k; ; ++i, --j)
176 poisf *= d2 / (i + 0.5f);
177 xtermf *= (x * (v / 2 + i - 1)) / (i);
179 T term = poisf * betaf;
184 pois *= (j + 0.5f) / d2;
186 xterm *= (j) / (x * (v / 2 + j - 1));
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)
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);
205 template <class T, class Policy>
206 T non_central_t_cdf(T v, T delta, T t, bool invert, const Policy& pol)
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);
215 // Otherwise, for t < 0 we have to use the reflection formula:
222 if(fabs(delta / (4 * v)) < policies::get_epsilon<T, Policy>())
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;
238 // x and y are the corresponding random
239 // variables for the noncentral beta distribution,
242 T x = t * t / (v + t * t);
243 T y = v / (v + t * t);
244 T d2 = delta * delta;
247 T c = a + b + d2 / 2;
249 // Crossover point for calculating p or q is the same
250 // as for the noncentral beta:
252 T cross = 1 - (b / c) * (1 + d2 / (2 * c * c));
261 result = non_central_beta_p(a, b, d2, x, y, pol);
262 result = non_central_t2_p(v, delta, x, y, pol, result);
267 result += cdf(boost::math::normal_distribution<T, Policy>(), -delta);
277 result = non_central_beta_q(a, b, d2, x, y, pol);
278 result = non_central_t2_q(v, delta, x, y, pol, result);
282 result = cdf(complement(boost::math::normal_distribution<T, Policy>(), -delta));
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&)
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<
298 policies::promote_float<false>,
299 policies::promote_double<false>,
300 policies::discrete_quantile<>,
301 policies::assert_undefined<> >::type forwarding_policy;
304 if(!detail::check_df_gt0_to_inf(
308 !detail::check_finite(
314 !detail::check_probability(
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);
328 return quantile(n, p);
332 return quantile(complement(n, q));
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);
340 guess = quantile(normal_distribution<value_type, forwarding_policy>(mean, var), p);
342 guess = quantile(complement(normal_distribution<value_type, forwarding_policy>(mean, var), q));
345 // We *must* get the sign of the initial guess correct,
346 // or our root-finder will fail, so double check it now:
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),
353 forwarding_policy());
356 s = boost::math::sign(p - pzero);
358 s = boost::math::sign(pzero - q);
359 if(s != boost::math::sign(guess))
361 guess = static_cast<T>(s);
364 value_type result = detail::generic_quantile(
365 non_central_t_distribution<value_type, forwarding_policy>(v, delta),
370 return policies::checked_narrowing_cast<T, forwarding_policy>(
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)
380 // Variables come first:
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;
386 // k is the starting point for iteration, and is the
387 // maximum of the poisson weighting term:
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:
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);
403 if((pois == 0) || (xterm == 0))
407 // Backwards recursion first, this is the stable
408 // direction for recursion:
410 boost::uintmax_t count = 0;
411 for(int i = k; i >= 0; --i)
413 T term = xterm * pois;
415 if(((fabs(term/sum) < errtol) && (i != k)) || (term == 0))
417 pois *= (i + 0.5f) / d2;
418 xterm *= (i) / (x * (n / 2 + i));
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);
427 for(int i = k + 1; ; ++i)
429 poisf *= d2 / (i + 0.5f);
430 xtermf *= (x * (n / 2 + i)) / (i);
431 T term = poisf * xtermf;
433 if((fabs(term/sum) < errtol) || (term == 0))
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);
446 template <class T, class Policy>
447 T non_central_t_pdf(T n, T delta, T t, const Policy& pol)
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);
456 // Otherwise, for t < 0 we have to use the reflection formula:
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
471 // The formula is simplified thanks to the relation
474 return tgamma_delta_ratio(n / 2 + 0.5f, T(0.5f))
475 * sqrt(n / constants::pi<T>())
476 * exp(-delta * delta / 2) / 2;
478 if(fabs(delta / (4 * n)) < policies::get_epsilon<T, Policy>())
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);
493 // x and y are the corresponding random
494 // variables for the noncentral beta distribution,
497 T x = t * t / (n + t * t);
498 T y = n / (n + t * t);
501 T d2 = delta * delta;
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);
515 template <class T, class Policy>
516 T mean(T v, T delta, const Policy& pol)
518 if ((boost::math::isinf)(v))
523 if (v > 1 / boost::math::tools::epsilon<T>() )
525 //normal_distribution<T, Policy> n(delta, 1);
526 //return boost::math::mean(n);
531 return delta * sqrt(v / 2) * tgamma_delta_ratio((v - 1) * 0.5f, T(0.5f), pol);
533 // Other moments use mean so using normal distribution is propagated.
536 template <class T, class Policy>
537 T variance(T v, T delta, const Policy& pol)
539 if ((boost::math::isinf)(v))
547 T result = ((delta * delta + 1) * v) / (v - 2);
548 T m = mean(v, delta, pol);
553 template <class T, class Policy>
554 T skewness(T v, T delta, const Policy& pol)
557 if ((boost::math::isinf)(v))
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;
569 result += v * (l2 + 2 * v - 3) / ((v - 3) * (v - 2));
571 result /= pow(var, T(1.5f));
575 template <class T, class Policy>
576 T kurtosis_excess(T v, T delta, const Policy& pol)
579 if ((boost::math::isinf)(v))
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;
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));
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.
603 template <class RealType, class Policy>
604 struct t_degrees_of_freedom_finder
606 t_degrees_of_freedom_finder(
607 RealType delta_, RealType x_, RealType p_, bool c)
608 : delta(delta_), x(x_), p(p_), comp(c) {}
610 RealType operator()(const RealType& v)
612 non_central_t_distribution<RealType, Policy> d(v, delta);
614 p - cdf(complement(d, x))
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)
628 const char* function = "non_central_t<%1%>::find_degrees_of_freedom";
629 if((p == 0) || (q == 0))
632 // Can't a thing if one of p and q is zero:
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());
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>();
642 // Pick an initial guess:
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>())
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());
656 template <class RealType, class Policy>
657 struct t_non_centrality_finder
659 t_non_centrality_finder(
660 RealType v_, RealType x_, RealType p_, bool c)
661 : v(v_), x(x_), p(p_), comp(c) {}
663 RealType operator()(const RealType& delta)
665 non_central_t_distribution<RealType, Policy> d(v, delta);
667 p - cdf(complement(d, x))
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)
681 const char* function = "non_central_t<%1%>::find_t_non_centrality";
682 if((p == 0) || (q == 0))
685 // Can't do a thing if one of p and q is zero:
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());
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>();
695 // Pick an initial guess that we know is the right side of
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>())
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());
714 } // namespace detail ======================================================================
716 template <class RealType = double, class Policy = policies::policy<> >
717 class non_central_t_distribution
720 typedef RealType value_type;
721 typedef Policy policy_type;
723 non_central_t_distribution(RealType v_, RealType lambda) : v(v_), ncp(lambda)
725 const char* function = "boost::math::non_central_t_distribution<%1%>::non_central_t_distribution(%1%,%1%)";
727 detail::check_df_gt0_to_inf(
730 detail::check_finite(
735 } // non_central_t_distribution constructor.
737 RealType degrees_of_freedom() const
738 { // Private data getter function.
741 RealType non_centrality() const
742 { // Private data getter function.
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.
750 static RealType find_degrees_of_freedom(RealType delta, RealType x, RealType p)
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<
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>(
770 template <class A, class B, class C>
771 static RealType find_degrees_of_freedom(const complemented3_type<A,B,C>& c)
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<
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>(
791 static RealType find_non_centrality(RealType v, RealType x, RealType p)
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<
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>(
811 template <class A, class B, class C>
812 static RealType find_non_centrality(const complemented3_type<A,B,C>& c)
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<
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>(
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
839 typedef non_central_t_distribution<double> non_central_t; // Reserved name of type double.
841 // Non-member functions to give properties of the distribution.
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>());
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>());
858 template <class RealType, class Policy>
859 inline RealType mode(const non_central_t_distribution<RealType, Policy>& dist)
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();
865 if(!detail::check_df_gt0_to_inf(
869 !detail::check_finite(
878 RealType m = v < 3 ? 0 : detail::mean(v, l, Policy());
879 RealType var = v < 4 ? 1 : detail::variance(v, l, Policy());
881 return detail::generic_find_mode(
888 template <class RealType, class Policy>
889 inline RealType mean(const non_central_t_distribution<RealType, Policy>& dist)
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<
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();
903 if(!detail::check_df_gt0_to_inf(
907 !detail::check_finite(
914 return policies::raise_domain_error<RealType>(
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);
923 template <class RealType, class Policy>
924 inline RealType variance(const non_central_t_distribution<RealType, Policy>& dist)
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<
930 policies::promote_float<false>,
931 policies::promote_double<false>,
932 policies::discrete_quantile<>,
933 policies::assert_undefined<> >::type forwarding_policy;
935 RealType v = dist.degrees_of_freedom();
936 RealType l = dist.non_centrality();
938 if(!detail::check_df_gt0_to_inf(
942 !detail::check_finite(
949 return policies::raise_domain_error<RealType>(
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);
956 // RealType standard_deviation(const non_central_t_distribution<RealType, Policy>& dist)
957 // standard_deviation provided by derived accessors.
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<
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();
973 if(!detail::check_df_gt0_to_inf(
977 !detail::check_finite(
984 return policies::raise_domain_error<RealType>(
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);
991 template <class RealType, class Policy>
992 inline RealType kurtosis_excess(const non_central_t_distribution<RealType, Policy>& dist)
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<
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();
1005 if(!detail::check_df_gt0_to_inf(
1009 !detail::check_finite(
1016 return policies::raise_domain_error<RealType>(
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
1023 template <class RealType, class Policy>
1024 inline RealType kurtosis(const non_central_t_distribution<RealType, Policy>& dist)
1026 return kurtosis_excess(dist) + 3;
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<
1036 policies::promote_float<false>,
1037 policies::promote_double<false>,
1038 policies::discrete_quantile<>,
1039 policies::assert_undefined<> >::type forwarding_policy;
1041 RealType v = dist.degrees_of_freedom();
1042 RealType l = dist.non_centrality();
1044 if(!detail::check_df_gt0_to_inf(
1048 !detail::check_finite(
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),
1068 template <class RealType, class Policy>
1069 RealType cdf(const non_central_t_distribution<RealType, Policy>& dist, const RealType& x)
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<
1076 policies::promote_float<false>,
1077 policies::promote_double<false>,
1078 policies::discrete_quantile<>,
1079 policies::assert_undefined<> >::type forwarding_policy;
1081 RealType v = dist.degrees_of_freedom();
1082 RealType l = dist.non_centrality();
1084 if(!detail::check_df_gt0_to_inf(
1088 !detail::check_finite(
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);
1104 //return cdf(normal_distribution<RealType, Policy>(l, 1), x);
1108 { // NO non-centrality, so use Student's t instead.
1109 return cdf(students_t_distribution<RealType, Policy>(v), x);
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),
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<
1128 policies::promote_float<false>,
1129 policies::promote_double<false>,
1130 policies::discrete_quantile<>,
1131 policies::assert_undefined<> >::type forwarding_policy;
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
1138 if(!detail::check_df_gt0_to_inf(
1142 !detail::check_finite(
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));
1161 { // zero non-centrality so use Student's t distribution.
1162 return cdf(complement(students_t_distribution<RealType, Policy>(v), x));
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),
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());
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.
1194 } // namespace boost
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>
1201 #endif // BOOST_MATH_SPECIAL_NON_CENTRAL_T_HPP