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