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