math/distributions/non_central_t.hpp

0001 // boost\math\distributions\non_central_t.hpp
0002
0003 // Copyright John Maddock 2008.
0004
0005 // Use, modification and distribution are subject to the
0006 // Boost Software License, Version 1.0.
0007 // (See accompanying file LICENSE_1_0.txt
0008 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0009
0010 #ifndef BOOST_MATH_SPECIAL_NON_CENTRAL_T_HPP
0011 #define BOOST_MATH_SPECIAL_NON_CENTRAL_T_HPP
0012
0013 #include <boost/math/distributions/fwd.hpp>
0014 #include <boost/math/distributions/non_central_beta.hpp> // for nc beta
0015 #include <boost/math/distributions/normal.hpp> // for normal CDF and quantile
0016 #include <boost/math/distributions/students_t.hpp>
0017 #include <boost/math/distributions/detail/generic_quantile.hpp> // quantile
0018 #include <boost/math/special_functions/trunc.hpp>
0019 #include <boost/math/special_functions/detail/hypergeometric_series.hpp>
0020 #include <boost/math/quadrature/exp_sinh.hpp>
0021
0022 namespace boost
0023 {
0024    namespace math
0025    {
0026
0027       template <class RealType, class Policy>
0028       class non_central_t_distribution;
0029
0030       namespace detail{
0031
0032          template <class T, class Policy>
0033          T non_central_t2_p(T v, T delta, T x, T y, const Policy& pol, T init_val)
0034          {
0035             BOOST_MATH_STD_USING
0036             //
0037             // Variables come first:
0038             //
0039             std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0040             T errtol = policies::get_epsilon<T, Policy>();
0041             T d2 = delta * delta / 2;
0042             //
0043             // k is the starting point for iteration, and is the
0044             // maximum of the poisson weighting term, we don't
0045             // ever allow k == 0 as this can lead to catastrophic
0046             // cancellation errors later (test case is v = 1621286869049072.3
0047             // delta = 0.16212868690490723, x = 0.86987415482475994).
0048             //
0049             long long k = lltrunc(d2);
0050             T pois;
0051             if(k == 0) k = 1;
0052             // Starting Poisson weight:
0053             pois = gamma_p_derivative(T(k+1), d2, pol)
0054                * tgamma_delta_ratio(T(k + 1), T(0.5f))
0055                * delta / constants::root_two<T>();
0056             if(pois == 0)
0057                return init_val;
0058             T xterm, beta;
0059             // Recurrence & starting beta terms:
0060             beta = x < y
0061                ? detail::ibeta_imp(T(k + 1), T(v / 2), x, pol, false, true, &xterm)
0062                : detail::ibeta_imp(T(v / 2), T(k + 1), y, pol, true, true, &xterm);
0063
0064             while (fabs(beta * pois) < tools::min_value<T>())
0065             {
0066                if ((k == 0) || (pois == 0))
0067                   return init_val;
0068                k /= 2;
0069                pois = gamma_p_derivative(T(k + 1), d2, pol)
0070                   * tgamma_delta_ratio(T(k + 1), T(0.5f))
0071                   * delta / constants::root_two<T>();
0072                beta = x < y
0073                   ? detail::ibeta_imp(T(k + 1), T(v / 2), x, pol, false, true, &xterm)
0074                   : detail::ibeta_imp(T(v / 2), T(k + 1), y, pol, true, true, &xterm);
0075             }
0076
0077             xterm *= y / (v / 2 + k);
0078             T poisf(pois), betaf(beta), xtermf(xterm);
0079             T sum = init_val;
0080             if((xterm == 0) && (beta == 0))
0081                return init_val;
0082
0083             //
0084             // Backwards recursion first, this is the stable
0085             // direction for recursion:
0086             //
0087             std::uintmax_t count = 0;
0088             T last_term = 0;
0089             for(auto i = k; i >= 0; --i)
0090             {
0091                T term = beta * pois;
0092                sum += term;
0093                // Don't terminate on first term in case we "fixed" k above:
0094                if(((fabs(last_term) > fabs(term)) && fabs(term/sum) < errtol) || (v == 2 && i == 0))
0095                   break;
0096                last_term = term;
0097                pois *= (i + 0.5f) / d2;
0098                beta += xterm;
0099                xterm *= (i) / (x * (v / 2 + i - 1));
0100                ++count;
0101             }
0102             last_term = 0;
0103             for(auto i = k + 1; ; ++i)
0104             {
0105                poisf *= d2 / (i + 0.5f);
0106                xtermf *= (x * (v / 2 + i - 1)) / (i);
0107                betaf -= xtermf;
0108                T term = poisf * betaf;
0109                sum += term;
0110                if((fabs(last_term) >= fabs(term)) && (fabs(term/sum) < errtol))
0111                   break;
0112                last_term = term;
0113                ++count;
0114                if(count > max_iter)
0115                {
0116                   return policies::raise_evaluation_error("cdf(non_central_t_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE
0117                }
0118             }
0119             return sum;
0120          }
0121
0122          template <class T, class Policy>
0123          T non_central_t2_q(T v, T delta, T x, T y, const Policy& pol, T init_val)
0124          {
0125             BOOST_MATH_STD_USING
0126             //
0127             // Variables come first:
0128             //
0129             std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0130             T errtol = boost::math::policies::get_epsilon<T, Policy>();
0131             T d2 = delta * delta / 2;
0132             //
0133             // k is the starting point for iteration, and is the
0134             // maximum of the poisson weighting term, we don't allow
0135             // k == 0 as this can cause catastrophic cancellation errors
0136             // (test case is v = 561908036470413.25, delta = 0.056190803647041321,
0137             // x = 1.6155232703966216):
0138             //
0139             long long k = lltrunc(d2);
0140             if(k == 0) k = 1;
0141             // Starting Poisson weight:
0142             T pois;
0143             if((k < static_cast<long long>(max_factorial<T>::value)) && (d2 < tools::log_max_value<T>()) && (log(d2) * k < tools::log_max_value<T>()))
0144             {
0145                //
0146                // For small k we can optimise this calculation by using
0147                // a simpler reduced formula:
0148                //
0149                pois = exp(-d2);
0150                pois *= pow(d2, static_cast<T>(k));
0151                pois /= boost::math::tgamma(T(k + 1 + 0.5), pol);
0152                pois *= delta / constants::root_two<T>();
0153             }
0154             else
0155             {
0156                pois = gamma_p_derivative(T(k+1), d2, pol)
0157                   * tgamma_delta_ratio(T(k + 1), T(0.5f))
0158                   * delta / constants::root_two<T>();
0159             }
0160             if(pois == 0)
0161                return init_val;
0162             // Recurance term:
0163             T xterm;
0164             T beta;
0165             // Starting beta term:
0166             if(k != 0)
0167             {
0168                beta = x < y
0169                   ? detail::ibeta_imp(T(k + 1), T(v / 2), x, pol, true, true, &xterm)
0170                   : detail::ibeta_imp(T(v / 2), T(k + 1), y, pol, false, true, &xterm);
0171
0172                xterm *= y / (v / 2 + k);
0173             }
0174             else
0175             {
0176                beta = pow(y, v / 2);
0177                xterm = beta;
0178             }
0179             T poisf(pois), betaf(beta), xtermf(xterm);
0180             T sum = init_val;
0181             if((xterm == 0) && (beta == 0))
0182                return init_val;
0183
0184             //
0185             // Fused forward and backwards recursion:
0186             //
0187             std::uintmax_t count = 0;
0188             T last_term = 0;
0189             for(auto i = k + 1, j = k; ; ++i, --j)
0190             {
0191                poisf *= d2 / (i + 0.5f);
0192                xtermf *= (x * (v / 2 + i - 1)) / (i);
0193                betaf += xtermf;
0194                T term = poisf * betaf;
0195
0196                if(j >= 0)
0197                {
0198                   term += beta * pois;
0199                   pois *= (j + 0.5f) / d2;
0200                   beta -= xterm;
0201                   if(!(v == 2 && j == 0))
0202                      xterm *= (j) / (x * (v / 2 + j - 1));
0203                }
0204
0205                sum += term;
0206                // Don't terminate on first term in case we "fixed" the value of k above:
0207                if((fabs(last_term) > fabs(term)) && fabs(term/sum) < errtol)
0208                   break;
0209                last_term = term;
0210                if(count > max_iter)
0211                {
0212                   return policies::raise_evaluation_error("cdf(non_central_t_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE
0213                }
0214                ++count;
0215             }
0216             return sum;
0217          }
0218
0219          template <class T, class Policy>
0220          T non_central_t_cdf(T v, T delta, T t, bool invert, const Policy& pol)
0221          {
0222             BOOST_MATH_STD_USING
0223             if ((boost::math::isinf)(v))
0224             { // Infinite degrees of freedom, so use normal distribution located at delta.
0225                normal_distribution<T, Policy> n(delta, 1);
0226                return cdf(n, t);
0227             }
0228             //
0229             // Otherwise, for t < 0 we have to use the reflection formula:
0230             if(t < 0)
0231             {
0232                t = -t;
0233                delta = -delta;
0234                invert = !invert;
0235             }
0236             if(fabs(delta / (4 * v)) < policies::get_epsilon<T, Policy>())
0237             {
0238                // Approximate with a Student's T centred on delta,
0239                // the crossover point is based on eq 2.6 from
0240                // "A Comparison of Approximations To Percentiles of the
0241                // Noncentral t-Distribution".  H. Sahai and M. M. Ojeda,
0242                // Revista Investigacion Operacional Vol 21, No 2, 2000.
0243                // Original sources referenced in the above are:
0244                // "Some Approximations to the Percentage Points of the Noncentral
0245                // t-Distribution". C. van Eeden. International Statistical Review, 29, 4-31.
0246                // "Continuous Univariate Distributions".  N.L. Johnson, S. Kotz and
0247                // N. Balkrishnan. 1995. John Wiley and Sons New York.
0248                T result = cdf(students_t_distribution<T, Policy>(v), t - delta);
0249                return invert ? 1 - result : result;
0250             }
0251             //
0252             // x and y are the corresponding random
0253             // variables for the noncentral beta distribution,
0254             // with y = 1 - x:
0255             //
0256             T x = t * t / (v + t * t);
0257             T y = v / (v + t * t);
0258             T d2 = delta * delta;
0259             T a = 0.5f;
0260             T b = v / 2;
0261             T c = a + b + d2 / 2;
0262             //
0263             // Crossover point for calculating p or q is the same
0264             // as for the noncentral beta:
0265             //
0266             T cross = 1 - (b / c) * (1 + d2 / (2 * c * c));
0267             T result;
0268             if(x < cross)
0269             {
0270                //
0271                // Calculate p:
0272                //
0273                if(x != 0)
0274                {
0275                   result = non_central_beta_p(a, b, d2, x, y, pol);
0276                   result = non_central_t2_p(v, delta, x, y, pol, result);
0277                   result /= 2;
0278                }
0279                else
0280                   result = 0;
0281                if (invert)
0282                {
0283                   result = cdf(complement(boost::math::normal_distribution<T, Policy>(), -delta)) - result;
0284                   invert = false;
0285                }
0286                else
0287                   result += cdf(boost::math::normal_distribution<T, Policy>(), -delta);
0288             }
0289             else
0290             {
0291                //
0292                // Calculate q:
0293                //
0294                invert = !invert;
0295                if(x != 0)
0296                {
0297                   result = non_central_beta_q(a, b, d2, x, y, pol);
0298                   result = non_central_t2_q(v, delta, x, y, pol, result);
0299                   result /= 2;
0300                }
0301                else // x == 0
0302                   result = cdf(complement(boost::math::normal_distribution<T, Policy>(), -delta));
0303             }
0304             if(invert)
0305                result = 1 - result;
0306             return result;
0307          }
0308
0309          template <class T, class Policy>
0310          T non_central_t_quantile(const char* function, T v, T delta, T p, T q, const Policy&)
0311          {
0312             BOOST_MATH_STD_USING
0313      //       static const char* function = "quantile(non_central_t_distribution<%1%>, %1%)";
0314      // now passed as function
0315             typedef typename policies::evaluation<T, Policy>::type value_type;
0316             typedef typename policies::normalise<
0317                Policy,
0318                policies::promote_float<false>,
0319                policies::promote_double<false>,
0320                policies::discrete_quantile<>,
0321                policies::assert_undefined<> >::type forwarding_policy;
0322
0323                T r;
0324                if(!detail::check_df_gt0_to_inf(
0325                   function,
0326                   v, &r, Policy())
0327                   ||
0328                !detail::check_non_centrality(
0329                   function,
0330                   static_cast<T>(delta * delta),
0331                   &r,
0332                   Policy())
0333                   ||
0334                !detail::check_probability(
0335                   function,
0336                   p,
0337                   &r,
0338                   Policy()))
0339                      return r;
0340
0341
0342             value_type guess = 0;
0343             if ( ((boost::math::isinf)(v)) || (v > 1 / boost::math::tools::epsilon<T>()) )
0344             { // Infinite or very large degrees of freedom, so use normal distribution located at delta.
0345                normal_distribution<T, Policy> n(delta, 1);
0346                if (p < q)
0347                {
0348                  return quantile(n, p);
0349                }
0350                else
0351                {
0352                  return quantile(complement(n, q));
0353                }
0354             }
0355             else if(v > 3)
0356             { // Use normal distribution to calculate guess.
0357                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));
0358                value_type var = (v > 1 / policies::get_epsilon<T, Policy>()) ? value_type(1) : (((delta * delta + 1) * v) / (v - 2) - mean * mean);
0359                if(p < q)
0360                   guess = quantile(normal_distribution<value_type, forwarding_policy>(mean, var), p);
0361                else
0362                   guess = quantile(complement(normal_distribution<value_type, forwarding_policy>(mean, var), q));
0363             }
0364             //
0365             // We *must* get the sign of the initial guess correct,
0366             // or our root-finder will fail, so double check it now:
0367             //
0368             value_type pzero = non_central_t_cdf(
0369                static_cast<value_type>(v),
0370                static_cast<value_type>(delta),
0371                static_cast<value_type>(0),
0372                !(p < q),
0373                forwarding_policy());
0374             int s;
0375             if(p < q)
0376                s = boost::math::sign(p - pzero);
0377             else
0378                s = boost::math::sign(pzero - q);
0379             if(s != boost::math::sign(guess))
0380             {
0381                guess = static_cast<T>(s);
0382             }
0383
0384             value_type result = detail::generic_quantile(
0385                non_central_t_distribution<value_type, forwarding_policy>(v, delta),
0386                (p < q ? p : q),
0387                guess,
0388                (p >= q),
0389                function);
0390             return policies::checked_narrowing_cast<T, forwarding_policy>(
0391                result,
0392                function);
0393          }
0394
0395          template <class T, class Policy>
0396          T non_central_t_pdf_integral(T x, T v, T mu, const Policy& pol)
0397          {
0398             BOOST_MATH_STD_USING
0399             boost::math::quadrature::exp_sinh<T, Policy> integrator;
0400             T integral = pow(v, v / 2) * exp(-v * mu * mu / (2 * (x * x + v)));
0401             if (integral != 0)
0402             {
0403                integral *= integrator.integrate([&x, v, mu](T y)
0404                   {
0405                      T p;
0406                      if (v * log(y) < tools::log_max_value<T>())
0407                         p = pow(y, v) * exp(boost::math::pow<2>((y - mu * x / sqrt(x * x + v))) / -2);
0408                      else
0409                         p = exp(log(y) * v + boost::math::pow<2>((y - mu * x / sqrt(x * x + v))) / -2);
0410                      return p;
0411                   });
0412             }
0413             integral /= boost::math::constants::root_pi<T>() * boost::math::tgamma(v / 2, pol) * pow(T(2), (v - 1) / 2) * pow(x * x + v, (v + 1) / 2);
0414             return integral;
0415          }
0416
0417          template <class T, class Policy>
0418          T non_central_t_pdf_hypergeometric(T x, T v, T mu, const Policy& pol)
0419          {
0420             BOOST_MATH_STD_USING
0421             long long scale = 0;
0422             const char* function = "non central T PDF";
0423             //
0424             // We only call this routine when we know that the series form of 1F1 is cheap to evaluate,
0425             // so no need to call the whole 1F1 function, just the series will do:
0426             //
0427             T Av = hypergeometric_1F1_generic_series(static_cast<T>((v + 1) / 2), boost::math::constants::half<T>(), static_cast<T>(mu * mu * x * x / (2 * (x * x + v))), pol, scale, function);
0428             Av = ldexp(Av, static_cast<int>(scale));
0429             scale = 0;
0430             T Bv = hypergeometric_1F1_generic_series(static_cast<T>(v / 2 + T(1)), static_cast<T>(T(3) / 2), static_cast<T>(mu * mu * x * x / (2 * (x * x + v))), pol, scale, function);
0431             Bv = ldexp(Bv, static_cast<int>(scale));
0432             Bv *= boost::math::tgamma_delta_ratio(v / 2 + T(1), -constants::half<T>(), pol);
0433             Bv *= boost::math::constants::root_two<T>() * mu * x / sqrt(x * x + v);
0434
0435             T tolerance = tools::root_epsilon<T>() * Av * 4;
0436             Av += Bv;
0437
0438             if (Av < tolerance)
0439             {
0440                // More than half the digits have cancelled, fall back to the integral method:
0441                return non_central_t_pdf_integral(x, v, mu, pol);
0442             }
0443
0444             Av *= exp(-mu * mu / 2) * pow(1 + x * x / v, -(v + 1) / 2) * boost::math::tgamma_delta_ratio(v / 2 + constants::half<T>(), -constants::half<T>(), pol);
0445             Av /= sqrt(v) * boost::math::constants::root_pi<T>();
0446             return Av;
0447          }
0448
0449          template <class T, class Policy>
0450          T non_central_t2_pdf(T n, T delta, T x, T y, const Policy& pol, T init_val)
0451          {
0452             BOOST_MATH_STD_USING
0453             //
0454             // Variables come first:
0455             //
0456             std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0457             T errtol = boost::math::policies::get_epsilon<T, Policy>();
0458             T d2 = delta * delta / 2;
0459             //
0460             // k is the starting point for iteration, and is the
0461             // maximum of the poisson weighting term:
0462             //
0463             long long k = lltrunc(d2);
0464             T pois, xterm;
0465             if(k == 0)
0466                k = 1;
0467             // Starting Poisson weight:
0468             pois = gamma_p_derivative(T(k+1), d2, pol)
0469                * tgamma_delta_ratio(T(k + 1), T(0.5f))
0470                * delta / constants::root_two<T>();
0471             BOOST_MATH_INSTRUMENT_VARIABLE(pois);
0472             // Starting beta term:
0473             xterm = x < y
0474                ? ibeta_derivative(T(k + 1), n / 2, x, pol)
0475                : ibeta_derivative(n / 2, T(k + 1), y, pol);
0476             BOOST_MATH_INSTRUMENT_VARIABLE(xterm);
0477
0478             while (fabs(xterm * pois) < tools::min_value<T>())
0479             {
0480                if (k == 0)
0481                   return init_val;
0482                k /= 2;
0483                pois = gamma_p_derivative(T(k + 1), d2, pol)
0484                   * tgamma_delta_ratio(T(k + 1), T(0.5f))
0485                   * delta / constants::root_two<T>();
0486                BOOST_MATH_INSTRUMENT_VARIABLE(pois);
0487                // Starting beta term:
0488                xterm = x < y
0489                   ? ibeta_derivative(T(k + 1), n / 2, x, pol)
0490                   : ibeta_derivative(n / 2, T(k + 1), y, pol);
0491                BOOST_MATH_INSTRUMENT_VARIABLE(xterm);
0492             }
0493
0494             T poisf(pois), xtermf(xterm);
0495             T sum = init_val;
0496
0497             //
0498             // Backwards recursion first, this is the stable
0499             // direction for recursion:
0500             //
0501             std::uintmax_t count = 0;
0502             T old_ratio = 1;  // arbitrary "large" value
0503             for(auto i = k; i >= 0; --i)
0504             {
0505                T term = xterm * pois;
0506                sum += term;
0507                BOOST_MATH_INSTRUMENT_VARIABLE(sum);
0508                T ratio = fabs(term / sum);
0509                if(((ratio < errtol) && (i != k) && (ratio < old_ratio)) || (term == 0))
0510                   break;
0511                old_ratio = ratio;
0512                pois *= (i + 0.5f) / d2;
0513                xterm *= (i) / (x * (n / 2 + i));
0514                ++count;
0515                if(count > max_iter)
0516                {
0517                   return policies::raise_evaluation_error("pdf(non_central_t_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE
0518                }
0519             }
0520             BOOST_MATH_INSTRUMENT_VARIABLE(sum);
0521             old_ratio = 0;
0522             for(auto i = k + 1; ; ++i)
0523             {
0524                poisf *= d2 / (i + 0.5f);
0525                xtermf *= (x * (n / 2 + i)) / (i);
0526                T term = poisf * xtermf;
0527                sum += term;
0528                T ratio = fabs(term / sum);
0529                if(((ratio < errtol) && (ratio < old_ratio)) || (term == 0))
0530                   break;
0531                ++count;
0532                old_ratio = ratio;
0533                if(count > max_iter)
0534                {
0535                   return policies::raise_evaluation_error("pdf(non_central_t_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE
0536                }
0537             }
0538             BOOST_MATH_INSTRUMENT_VARIABLE(sum);
0539             return sum;
0540          }
0541
0542          template <class T, class Policy>
0543          T non_central_t_pdf(T n, T delta, T t, const Policy& pol)
0544          {
0545             BOOST_MATH_STD_USING
0546             if ((boost::math::isinf)(n))
0547             { // Infinite degrees of freedom, so use normal distribution located at delta.
0548                normal_distribution<T, Policy> norm(delta, 1);
0549                return pdf(norm, t);
0550             }
0551             if(t * t < tools::epsilon<T>())
0552             {
0553                //
0554                // Handle this as a special case, using the formula
0555                // from Weisstein, Eric W.
0556                // "Noncentral Student's t-Distribution."
0557                // From MathWorld--A Wolfram Web Resource.
0558                // http://mathworld.wolfram.com/NoncentralStudentst-Distribution.html
0559                //
0560                // The formula is simplified thanks to the relation
0561                // 1F1(a,b,0) = 1.
0562                //
0563                return tgamma_delta_ratio(n / 2 + 0.5f, T(0.5f))
0564                   * sqrt(n / constants::pi<T>())
0565                   * exp(-delta * delta / 2) / 2;
0566             }
0567             if(fabs(delta / (4 * n)) < policies::get_epsilon<T, Policy>())
0568             {
0569                // Approximate with a Student's T centred on delta,
0570                // the crossover point is based on eq 2.6 from
0571                // "A Comparison of Approximations To Percentiles of the
0572                // Noncentral t-Distribution".  H. Sahai and M. M. Ojeda,
0573                // Revista Investigacion Operacional Vol 21, No 2, 2000.
0574                // Original sources referenced in the above are:
0575                // "Some Approximations to the Percentage Points of the Noncentral
0576                // t-Distribution". C. van Eeden. International Statistical Review, 29, 4-31.
0577                // "Continuous Univariate Distributions".  N.L. Johnson, S. Kotz and
0578                // N. Balkrishnan. 1995. John Wiley and Sons New York.
0579                return pdf(students_t_distribution<T, Policy>(n), t - delta);
0580             }
0581             //
0582             // Figure out if the hypergeometric formula will be efficient or not,
0583             // based on where the summit of the series.
0584             //
0585             T a = (n + 1) / 2;
0586             T x = delta * delta * t * t / (2 * (t * t + n));
0587             T summit = (sqrt(x * (4 * a + x)) + x) / 2;
0588             if (summit < 40)
0589                return non_central_t_pdf_hypergeometric(t, n, delta, pol);
0590             //
0591             // Otherwise, for t < 0 we have to use the reflection formula:
0592             //
0593             if (t < 0)
0594             {
0595                t = -t;
0596                delta = -delta;
0597             }
0598             //
0599             // x and y are the corresponding random
0600             // variables for the noncentral beta distribution,
0601             // with y = 1 - x:
0602             //
0603             x = t * t / (n + t * t);
0604             T y = n / (n + t * t);
0605             a = 0.5f;
0606             T b = n / 2;
0607             T d2 = delta * delta;
0608             //
0609             // Calculate pdf:
0610             //
0611             T dt = n * t / (n * n + 2 * n * t * t + t * t * t * t);
0612             BOOST_MATH_INSTRUMENT_VARIABLE(dt);
0613             T result = non_central_beta_pdf(a, b, d2, x, y, pol);
0614             BOOST_MATH_INSTRUMENT_VARIABLE(result);
0615             T tol = tools::root_epsilon<T>() * result;
0616             result = non_central_t2_pdf(n, delta, x, y, pol, result);
0617             BOOST_MATH_INSTRUMENT_VARIABLE(result);
0618             result *= dt;
0619             if (result <= tol)
0620             {
0621                // More than half the digits in the result have cancelled,
0622                // Try direct integration... super slow but reliable as far as we can tell...
0623                if (delta < 0)
0624                {
0625                   // reflect back:
0626                   delta = -delta;
0627                   t = -t;
0628                }
0629                result = non_central_t_pdf_integral(t, n, delta, pol);
0630             }
0631             return result;
0632          }
0633
0634          template <class T, class Policy>
0635          T mean(T v, T delta, const Policy& pol)
0636          {
0637             if ((boost::math::isinf)(v))
0638             {
0639                return delta;
0640             }
0641             BOOST_MATH_STD_USING
0642             if (v > 1 / boost::math::tools::epsilon<T>() )
0643             {
0644               //normal_distribution<T, Policy> n(delta, 1);
0645               //return boost::math::mean(n);
0646               return delta;
0647             }
0648             else
0649             {
0650              return delta * sqrt(v / 2) * tgamma_delta_ratio((v - 1) * 0.5f, T(0.5f), pol);
0651             }
0652             // Other moments use mean so using normal distribution is propagated.
0653          }
0654
0655          template <class T, class Policy>
0656          T variance(T v, T delta, const Policy& pol)
0657          {
0658             if ((boost::math::isinf)(v))
0659             {
0660                return 1;
0661             }
0662             if (delta == 0)
0663             {  // == Student's t
0664               return v / (v - 2);
0665             }
0666             T result = ((delta * delta + 1) * v) / (v - 2);
0667             T m = mean(v, delta, pol);
0668             result -= m * m;
0669             return result;
0670          }
0671
0672          template <class T, class Policy>
0673          T skewness(T v, T delta, const Policy& pol)
0674          {
0675             BOOST_MATH_STD_USING
0676             if ((boost::math::isinf)(v))
0677             {
0678                return 0;
0679             }
0680             if(delta == 0)
0681             { // == Student's t
0682               return 0;
0683             }
0684             T mean = boost::math::detail::mean(v, delta, pol);
0685             T l2 = delta * delta;
0686             T var = ((l2 + 1) * v) / (v - 2) - mean * mean;
0687             T result = -2 * var;
0688             result += v * (l2 + 2 * v - 3) / ((v - 3) * (v - 2));
0689             result *= mean;
0690             result /= pow(var, T(1.5f));
0691             return result;
0692          }
0693
0694          template <class T, class Policy>
0695          T kurtosis_excess(T v, T delta, const Policy& pol)
0696          {
0697             BOOST_MATH_STD_USING
0698             if ((boost::math::isinf)(v))
0699             {
0700                return 1;
0701             }
0702             if (delta == 0)
0703             { // == Student's t
0704               return 1;
0705             }
0706             T mean = boost::math::detail::mean(v, delta, pol);
0707             T l2 = delta * delta;
0708             T var = ((l2 + 1) * v) / (v - 2) - mean * mean;
0709             T result = -3 * var;
0710             result += v * (l2 * (v + 1) + 3 * (3 * v - 5)) / ((v - 3) * (v - 2));
0711             result *= -mean * mean;
0712             result += v * v * (l2 * l2 + 6 * l2 + 3) / ((v - 4) * (v - 2));
0713             result /= var * var;
0714             result -= static_cast<T>(3);
0715             return result;
0716          }
0717
0718 #if 0
0719          //
0720          // This code is disabled, since there can be multiple answers to the
0721          // question, and it's not clear how to find the "right" one.
0722          //
0723          template <class RealType, class Policy>
0724          struct t_degrees_of_freedom_finder
0725          {
0726             t_degrees_of_freedom_finder(
0727                RealType delta_, RealType x_, RealType p_, bool c)
0728                : delta(delta_), x(x_), p(p_), comp(c) {}
0729
0730             RealType operator()(const RealType& v)
0731             {
0732                non_central_t_distribution<RealType, Policy> d(v, delta);
0733                return comp ?
0734                   p - cdf(complement(d, x))
0735                   : cdf(d, x) - p;
0736             }
0737          private:
0738             RealType delta;
0739             RealType x;
0740             RealType p;
0741             bool comp;
0742          };
0743
0744          template <class RealType, class Policy>
0745          inline RealType find_t_degrees_of_freedom(
0746             RealType delta, RealType x, RealType p, RealType q, const Policy& pol)
0747          {
0748             const char* function = "non_central_t<%1%>::find_degrees_of_freedom";
0749             if((p == 0) || (q == 0))
0750             {
0751                //
0752                // Can't a thing if one of p and q is zero:
0753                //
0754                return policies::raise_evaluation_error<RealType>(function, "Can't find degrees of freedom when the probability is 0 or 1, only possible answer is %1%", // LCOV_EXCL_LINE
0755                   RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy()); // LCOV_EXCL_LINE
0756             }
0757             t_degrees_of_freedom_finder<RealType, Policy> f(delta, x, p < q ? p : q, p < q ? false : true);
0758             tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
0759             std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
0760             //
0761             // Pick an initial guess:
0762             //
0763             RealType guess = 200;
0764             std::pair<RealType, RealType> ir = tools::bracket_and_solve_root(
0765                f, guess, RealType(2), false, tol, max_iter, pol);
0766             RealType result = ir.first + (ir.second - ir.first) / 2;
0767             if(max_iter >= policies::get_max_root_iterations<Policy>())
0768             {
0769                return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE
0770                   " or there is no answer to problem.  Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE
0771             }
0772             return result;
0773          }
0774
0775          template <class RealType, class Policy>
0776          struct t_non_centrality_finder
0777          {
0778             t_non_centrality_finder(
0779                RealType v_, RealType x_, RealType p_, bool c)
0780                : v(v_), x(x_), p(p_), comp(c) {}
0781
0782             RealType operator()(const RealType& delta)
0783             {
0784                non_central_t_distribution<RealType, Policy> d(v, delta);
0785                return comp ?
0786                   p - cdf(complement(d, x))
0787                   : cdf(d, x) - p;
0788             }
0789          private:
0790             RealType v;
0791             RealType x;
0792             RealType p;
0793             bool comp;
0794          };
0795
0796          template <class RealType, class Policy>
0797          inline RealType find_t_non_centrality(
0798             RealType v, RealType x, RealType p, RealType q, const Policy& pol)
0799          {
0800             const char* function = "non_central_t<%1%>::find_t_non_centrality";
0801             if((p == 0) || (q == 0))
0802             {
0803                //
0804                // Can't do a thing if one of p and q is zero:
0805                //
0806                return policies::raise_evaluation_error<RealType>(function, "Can't find non-centrality parameter when the probability is 0 or 1, only possible answer is %1%", // LCOV_EXCL_LINE
0807                   RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy()); // LCOV_EXCL_LINE
0808             }
0809             t_non_centrality_finder<RealType, Policy> f(v, x, p < q ? p : q, p < q ? false : true);
0810             tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
0811             std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
0812             //
0813             // Pick an initial guess that we know is the right side of
0814             // zero:
0815             //
0816             RealType guess;
0817             if(f(0) < 0)
0818                guess = 1;
0819             else
0820                guess = -1;
0821             std::pair<RealType, RealType> ir = tools::bracket_and_solve_root(
0822                f, guess, RealType(2), false, tol, max_iter, pol);
0823             RealType result = ir.first + (ir.second - ir.first) / 2;
0824             if(max_iter >= policies::get_max_root_iterations<Policy>())
0825             {
0826                return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE
0827                   " or there is no answer to problem.  Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE
0828             }
0829             return result;
0830          }
0831 #endif
0832       } // namespace detail ======================================================================
0833
0834       template <class RealType = double, class Policy = policies::policy<> >
0835       class non_central_t_distribution
0836       {
0837       public:
0838          typedef RealType value_type;
0839          typedef Policy policy_type;
0840
0841          non_central_t_distribution(RealType v_, RealType lambda) : v(v_), ncp(lambda)
0842          {
0843             const char* function = "boost::math::non_central_t_distribution<%1%>::non_central_t_distribution(%1%,%1%)";
0844             RealType r;
0845             detail::check_df_gt0_to_inf(
0846                function,
0847                v, &r, Policy());
0848             detail::check_non_centrality(
0849                function,
0850                static_cast<value_type>(lambda * lambda),
0851                &r,
0852                Policy());
0853          } // non_central_t_distribution constructor.
0854
0855          RealType degrees_of_freedom() const
0856          { // Private data getter function.
0857             return v;
0858          }
0859          RealType non_centrality() const
0860          { // Private data getter function.
0861             return ncp;
0862          }
0863 #if 0
0864          //
0865          // This code is disabled, since there can be multiple answers to the
0866          // question, and it's not clear how to find the "right" one.
0867          //
0868          static RealType find_degrees_of_freedom(RealType delta, RealType x, RealType p)
0869          {
0870             const char* function = "non_central_t<%1%>::find_degrees_of_freedom";
0871             typedef typename policies::evaluation<RealType, Policy>::type value_type;
0872             typedef typename policies::normalise<
0873                Policy,
0874                policies::promote_float<false>,
0875                policies::promote_double<false>,
0876                policies::discrete_quantile<>,
0877                policies::assert_undefined<> >::type forwarding_policy;
0878             value_type result = detail::find_t_degrees_of_freedom(
0879                static_cast<value_type>(delta),
0880                static_cast<value_type>(x),
0881                static_cast<value_type>(p),
0882                static_cast<value_type>(1-p),
0883                forwarding_policy());
0884             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
0885                result,
0886                function);
0887          }
0888          template <class A, class B, class C>
0889          static RealType find_degrees_of_freedom(const complemented3_type<A,B,C>& c)
0890          {
0891             const char* function = "non_central_t<%1%>::find_degrees_of_freedom";
0892             typedef typename policies::evaluation<RealType, Policy>::type value_type;
0893             typedef typename policies::normalise<
0894                Policy,
0895                policies::promote_float<false>,
0896                policies::promote_double<false>,
0897                policies::discrete_quantile<>,
0898                policies::assert_undefined<> >::type forwarding_policy;
0899             value_type result = detail::find_t_degrees_of_freedom(
0900                static_cast<value_type>(c.dist),
0901                static_cast<value_type>(c.param1),
0902                static_cast<value_type>(1-c.param2),
0903                static_cast<value_type>(c.param2),
0904                forwarding_policy());
0905             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
0906                result,
0907                function);
0908          }
0909          static RealType find_non_centrality(RealType v, RealType x, RealType p)
0910          {
0911             const char* function = "non_central_t<%1%>::find_t_non_centrality";
0912             typedef typename policies::evaluation<RealType, Policy>::type value_type;
0913             typedef typename policies::normalise<
0914                Policy,
0915                policies::promote_float<false>,
0916                policies::promote_double<false>,
0917                policies::discrete_quantile<>,
0918                policies::assert_undefined<> >::type forwarding_policy;
0919             value_type result = detail::find_t_non_centrality(
0920                static_cast<value_type>(v),
0921                static_cast<value_type>(x),
0922                static_cast<value_type>(p),
0923                static_cast<value_type>(1-p),
0924                forwarding_policy());
0925             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
0926                result,
0927                function);
0928          }
0929          template <class A, class B, class C>
0930          static RealType find_non_centrality(const complemented3_type<A,B,C>& c)
0931          {
0932             const char* function = "non_central_t<%1%>::find_t_non_centrality";
0933             typedef typename policies::evaluation<RealType, Policy>::type value_type;
0934             typedef typename policies::normalise<
0935                Policy,
0936                policies::promote_float<false>,
0937                policies::promote_double<false>,
0938                policies::discrete_quantile<>,
0939                policies::assert_undefined<> >::type forwarding_policy;
0940             value_type result = detail::find_t_non_centrality(
0941                static_cast<value_type>(c.dist),
0942                static_cast<value_type>(c.param1),
0943                static_cast<value_type>(1-c.param2),
0944                static_cast<value_type>(c.param2),
0945                forwarding_policy());
0946             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
0947                result,
0948                function);
0949          }
0950 #endif
0951       private:
0952          // Data member, initialized by constructor.
0953          RealType v;   // degrees of freedom
0954          RealType ncp; // non-centrality parameter
0955       }; // template <class RealType, class Policy> class non_central_t_distribution
0956
0957       typedef non_central_t_distribution<double> non_central_t; // Reserved name of type double.
0958
0959       #ifdef __cpp_deduction_guides
0960       template <class RealType>
0961       non_central_t_distribution(RealType,RealType)->non_central_t_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0962       #endif
0963
0964       // Non-member functions to give properties of the distribution.
0965
0966       template <class RealType, class Policy>
0967       inline const std::pair<RealType, RealType> range(const non_central_t_distribution<RealType, Policy>& /* dist */)
0968       { // Range of permissible values for random variable k.
0969          using boost::math::tools::max_value;
0970          return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>());
0971       }
0972
0973       template <class RealType, class Policy>
0974       inline const std::pair<RealType, RealType> support(const non_central_t_distribution<RealType, Policy>& /* dist */)
0975       { // Range of supported values for random variable k.
0976          // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
0977          using boost::math::tools::max_value;
0978          return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>());
0979       }
0980
0981       template <class RealType, class Policy>
0982       inline RealType mode(const non_central_t_distribution<RealType, Policy>& dist)
0983       { // mode.
0984          static const char* function = "mode(non_central_t_distribution<%1%> const&)";
0985          RealType v = dist.degrees_of_freedom();
0986          RealType l = dist.non_centrality();
0987          RealType r;
0988          if(!detail::check_df_gt0_to_inf(
0989             function,
0990             v, &r, Policy())
0991             ||
0992          !detail::check_non_centrality(
0993             function,
0994             static_cast<RealType>(l * l),
0995             &r,
0996             Policy()))
0997                return static_cast<RealType>(r);
0998
0999          BOOST_MATH_STD_USING
1000
1001          RealType m = v < 3 ? 0 : detail::mean(v, l, Policy());
1002          RealType var = v < 4 ? 1 : detail::variance(v, l, Policy());
1003
1004          return detail::generic_find_mode(
1005             dist,
1006             m,
1007             function,
1008             sqrt(var));
1009       }
1010
1011       template <class RealType, class Policy>
1012       inline RealType mean(const non_central_t_distribution<RealType, Policy>& dist)
1013       {
1014          BOOST_MATH_STD_USING
1015          const char* function = "mean(const non_central_t_distribution<%1%>&)";
1016          typedef typename policies::evaluation<RealType, Policy>::type value_type;
1017          typedef typename policies::normalise<
1018             Policy,
1019             policies::promote_float<false>,
1020             policies::promote_double<false>,
1021             policies::discrete_quantile<>,
1022             policies::assert_undefined<> >::type forwarding_policy;
1023          RealType v = dist.degrees_of_freedom();
1024          RealType l = dist.non_centrality();
1025          RealType r;
1026          if(!detail::check_df_gt0_to_inf(
1027             function,
1028             v, &r, Policy())
1029             ||
1030          !detail::check_non_centrality(
1031             function,
1032             static_cast<RealType>(l * l),
1033             &r,
1034             Policy()))
1035                return static_cast<RealType>(r);
1036          if(v <= 1)
1037             return policies::raise_domain_error<RealType>(
1038                function,
1039                "The non-central t distribution has no defined mean for degrees of freedom <= 1: got v=%1%.", v, Policy());
1040          // return l * sqrt(v / 2) * tgamma_delta_ratio((v - 1) * 0.5f, RealType(0.5f));
1041          return policies::checked_narrowing_cast<RealType, forwarding_policy>(
1042             detail::mean(static_cast<value_type>(v), static_cast<value_type>(l), forwarding_policy()), function);
1043
1044       } // mean
1045
1046       template <class RealType, class Policy>
1047       inline RealType variance(const non_central_t_distribution<RealType, Policy>& dist)
1048       { // variance.
1049          const char* function = "variance(const non_central_t_distribution<%1%>&)";
1050          typedef typename policies::evaluation<RealType, Policy>::type value_type;
1051          typedef typename policies::normalise<
1052             Policy,
1053             policies::promote_float<false>,
1054             policies::promote_double<false>,
1055             policies::discrete_quantile<>,
1056             policies::assert_undefined<> >::type forwarding_policy;
1057          BOOST_MATH_STD_USING
1058          RealType v = dist.degrees_of_freedom();
1059          RealType l = dist.non_centrality();
1060          RealType r;
1061          if(!detail::check_df_gt0_to_inf(
1062             function,
1063             v, &r, Policy())
1064             ||
1065          !detail::check_non_centrality(
1066             function,
1067             static_cast<RealType>(l * l),
1068             &r,
1069             Policy()))
1070                return static_cast<RealType>(r);
1071          if(v <= 2)
1072             return policies::raise_domain_error<RealType>(
1073                function,
1074                "The non-central t distribution has no defined variance for degrees of freedom <= 2: got v=%1%.", v, Policy());
1075          return policies::checked_narrowing_cast<RealType, forwarding_policy>(
1076             detail::variance(static_cast<value_type>(v), static_cast<value_type>(l), forwarding_policy()), function);
1077       }
1078
1079       // RealType standard_deviation(const non_central_t_distribution<RealType, Policy>& dist)
1080       // standard_deviation provided by derived accessors.
1081
1082       template <class RealType, class Policy>
1083       inline RealType skewness(const non_central_t_distribution<RealType, Policy>& dist)
1084       { // skewness = sqrt(l).
1085          const char* function = "skewness(const non_central_t_distribution<%1%>&)";
1086          typedef typename policies::evaluation<RealType, Policy>::type value_type;
1087          typedef typename policies::normalise<
1088             Policy,
1089             policies::promote_float<false>,
1090             policies::promote_double<false>,
1091             policies::discrete_quantile<>,
1092             policies::assert_undefined<> >::type forwarding_policy;
1093          RealType v = dist.degrees_of_freedom();
1094          RealType l = dist.non_centrality();
1095          RealType r;
1096          if(!detail::check_df_gt0_to_inf(
1097             function,
1098             v, &r, Policy())
1099             ||
1100          !detail::check_non_centrality(
1101             function,
1102             static_cast<RealType>(l * l),
1103             &r,
1104             Policy()))
1105                return static_cast<RealType>(r);
1106          if(v <= 3)
1107             return policies::raise_domain_error<RealType>(
1108                function,
1109                "The non-central t distribution has no defined skewness for degrees of freedom <= 3: got v=%1%.", v, Policy());;
1110          return policies::checked_narrowing_cast<RealType, forwarding_policy>(
1111             detail::skewness(static_cast<value_type>(v), static_cast<value_type>(l), forwarding_policy()), function);
1112       }
1113
1114       template <class RealType, class Policy>
1115       inline RealType kurtosis_excess(const non_central_t_distribution<RealType, Policy>& dist)
1116       {
1117          const char* function = "kurtosis_excess(const non_central_t_distribution<%1%>&)";
1118          typedef typename policies::evaluation<RealType, Policy>::type value_type;
1119          typedef typename policies::normalise<
1120             Policy,
1121             policies::promote_float<false>,
1122             policies::promote_double<false>,
1123             policies::discrete_quantile<>,
1124             policies::assert_undefined<> >::type forwarding_policy;
1125          RealType v = dist.degrees_of_freedom();
1126          RealType l = dist.non_centrality();
1127          RealType r;
1128          if(!detail::check_df_gt0_to_inf(
1129             function,
1130             v, &r, Policy())
1131             ||
1132          !detail::check_non_centrality(
1133             function,
1134             static_cast<RealType>(l * l),
1135             &r,
1136             Policy()))
1137                return static_cast<RealType>(r);
1138          if(v <= 4)
1139             return policies::raise_domain_error<RealType>(
1140                function,
1141                "The non-central t distribution has no defined kurtosis for degrees of freedom <= 4: got v=%1%.", v, Policy());;
1142          return policies::checked_narrowing_cast<RealType, forwarding_policy>(
1143             detail::kurtosis_excess(static_cast<value_type>(v), static_cast<value_type>(l), forwarding_policy()), function);
1144       } // kurtosis_excess
1145
1146       template <class RealType, class Policy>
1147       inline RealType kurtosis(const non_central_t_distribution<RealType, Policy>& dist)
1148       {
1149          return kurtosis_excess(dist) + 3;
1150       }
1151
1152       template <class RealType, class Policy>
1153       inline RealType pdf(const non_central_t_distribution<RealType, Policy>& dist, const RealType& t)
1154       { // Probability Density/Mass Function.
1155          const char* function = "pdf(non_central_t_distribution<%1%>, %1%)";
1156          typedef typename policies::evaluation<RealType, Policy>::type value_type;
1157          typedef typename policies::normalise<
1158             Policy,
1159             policies::promote_float<false>,
1160             policies::promote_double<false>,
1161             policies::discrete_quantile<>,
1162             policies::assert_undefined<> >::type forwarding_policy;
1163
1164          RealType v = dist.degrees_of_freedom();
1165          RealType l = dist.non_centrality();
1166          RealType r;
1167          if(!detail::check_df_gt0_to_inf(
1168             function,
1169             v, &r, Policy())
1170             ||
1171          !detail::check_non_centrality(
1172             function,
1173             static_cast<RealType>(l * l), // we need l^2 to be countable.
1174             &r,
1175             Policy())
1176             ||
1177          !detail::check_x(
1178             function,
1179             t,
1180             &r,
1181             Policy()))
1182                return static_cast<RealType>(r);
1183          return policies::checked_narrowing_cast<RealType, forwarding_policy>(
1184             detail::non_central_t_pdf(static_cast<value_type>(v),
1185                static_cast<value_type>(l),
1186                static_cast<value_type>(t),
1187                Policy()),
1188             function);
1189       } // pdf
1190
1191       template <class RealType, class Policy>
1192       RealType cdf(const non_central_t_distribution<RealType, Policy>& dist, const RealType& x)
1193       {
1194          const char* function = "boost::math::cdf(non_central_t_distribution<%1%>&, %1%)";
1195 //   was const char* function = "boost::math::non_central_t_distribution<%1%>::cdf(%1%)";
1196          typedef typename policies::evaluation<RealType, Policy>::type value_type;
1197          typedef typename policies::normalise<
1198             Policy,
1199             policies::promote_float<false>,
1200             policies::promote_double<false>,
1201             policies::discrete_quantile<>,
1202             policies::assert_undefined<> >::type forwarding_policy;
1203
1204          RealType v = dist.degrees_of_freedom();
1205          RealType l = dist.non_centrality();
1206          RealType r;
1207          if(!detail::check_df_gt0_to_inf(
1208             function,
1209             v, &r, Policy())
1210             ||
1211          !detail::check_non_centrality(
1212             function,
1213             static_cast<RealType>(l * l),
1214             &r,
1215             Policy())
1216             ||
1217          !detail::check_x(
1218             function,
1219             x,
1220             &r,
1221             Policy()))
1222                return static_cast<RealType>(r);
1223          if ((boost::math::isinf)(v))
1224           { // Infinite degrees of freedom, so use normal distribution located at delta.
1225              normal_distribution<RealType, Policy> n(l, 1);
1226              cdf(n, x);
1227               //return cdf(normal_distribution<RealType, Policy>(l, 1), x);
1228           }
1229
1230          if(l == 0)
1231          { // NO non-centrality, so use Student's t instead.
1232             return cdf(students_t_distribution<RealType, Policy>(v), x);
1233          }
1234          return policies::checked_narrowing_cast<RealType, forwarding_policy>(
1235             detail::non_central_t_cdf(
1236                static_cast<value_type>(v),
1237                static_cast<value_type>(l),
1238                static_cast<value_type>(x),
1239                false, Policy()),
1240             function);
1241       } // cdf
1242
1243       template <class RealType, class Policy>
1244       RealType cdf(const complemented2_type<non_central_t_distribution<RealType, Policy>, RealType>& c)
1245       { // Complemented Cumulative Distribution Function
1246   // was       const char* function = "boost::math::non_central_t_distribution<%1%>::cdf(%1%)";
1247          const char* function = "boost::math::cdf(const complement(non_central_t_distribution<%1%>&), %1%)";
1248          typedef typename policies::evaluation<RealType, Policy>::type value_type;
1249          typedef typename policies::normalise<
1250             Policy,
1251             policies::promote_float<false>,
1252             policies::promote_double<false>,
1253             policies::discrete_quantile<>,
1254             policies::assert_undefined<> >::type forwarding_policy;
1255
1256          non_central_t_distribution<RealType, Policy> const& dist = c.dist;
1257          RealType x = c.param;
1258          RealType v = dist.degrees_of_freedom();
1259          RealType l = dist.non_centrality(); // aka delta
1260          RealType r;
1261          if(!detail::check_df_gt0_to_inf(
1262             function,
1263             v, &r, Policy())
1264             ||
1265          !detail::check_non_centrality(
1266             function,
1267             static_cast<RealType>(l * l),
1268             &r,
1269             Policy())
1270             ||
1271          !detail::check_x(
1272             function,
1273             x,
1274             &r,
1275             Policy()))
1276                return static_cast<RealType>(r);
1277
1278          if ((boost::math::isinf)(v))
1279          { // Infinite degrees of freedom, so use normal distribution located at delta.
1280              normal_distribution<RealType, Policy> n(l, 1);
1281              return cdf(complement(n, x));
1282          }
1283          if(l == 0)
1284          { // zero non-centrality so use Student's t distribution.
1285             return cdf(complement(students_t_distribution<RealType, Policy>(v), x));
1286          }
1287          return policies::checked_narrowing_cast<RealType, forwarding_policy>(
1288             detail::non_central_t_cdf(
1289                static_cast<value_type>(v),
1290                static_cast<value_type>(l),
1291                static_cast<value_type>(x),
1292                true, Policy()),
1293             function);
1294       } // ccdf
1295
1296       template <class RealType, class Policy>
1297       inline RealType quantile(const non_central_t_distribution<RealType, Policy>& dist, const RealType& p)
1298       { // Quantile (or Percent Point) function.
1299          static const char* function = "quantile(const non_central_t_distribution<%1%>, %1%)";
1300          RealType v = dist.degrees_of_freedom();
1301          RealType l = dist.non_centrality();
1302          return detail::non_central_t_quantile(function, v, l, p, RealType(1-p), Policy());
1303       } // quantile
1304
1305       template <class RealType, class Policy>
1306       inline RealType quantile(const complemented2_type<non_central_t_distribution<RealType, Policy>, RealType>& c)
1307       { // Quantile (or Percent Point) function.
1308          static const char* function = "quantile(const complement(non_central_t_distribution<%1%>, %1%))";
1309          non_central_t_distribution<RealType, Policy> const& dist = c.dist;
1310          RealType q = c.param;
1311          RealType v = dist.degrees_of_freedom();
1312          RealType l = dist.non_centrality();
1313          return detail::non_central_t_quantile(function, v, l, RealType(1-q), q, Policy());
1314       } // quantile complement.
1315
1316    } // namespace math
1317 } // namespace boost
1318
1319 // This include must be at the end, *after* the accessors
1320 // for this distribution have been defined, in order to
1321 // keep compilers that support two-phase lookup happy.
1322 #include <boost/math/distributions/detail/derived_accessors.hpp>
1323
1324 #endif // BOOST_MATH_SPECIAL_NON_CENTRAL_T_HPP
1325