ceph/src/boost/libs/math/include/boost/math/distributions/non_central_t.hpp

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