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

   1 // boost\math\distributions\non_central_chi_squared.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_CHI_SQUARE_HPP
  11 #define BOOST_MATH_SPECIAL_NON_CENTRAL_CHI_SQUARE_HPP
  12
  13 #include <boost/math/distributions/fwd.hpp>
  14 #include <boost/math/special_functions/gamma.hpp> // for incomplete gamma. gamma_q
  15 #include <boost/math/special_functions/bessel.hpp> // for cyl_bessel_i
  16 #include <boost/math/special_functions/round.hpp> // for iround
  17 #include <boost/math/distributions/complement.hpp> // complements
  18 #include <boost/math/distributions/chi_squared.hpp> // central distribution
  19 #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks
  20 #include <boost/math/special_functions/fpclassify.hpp> // isnan.
  21 #include <boost/math/tools/roots.hpp> // for root finding.
  22 #include <boost/math/distributions/detail/generic_mode.hpp>
  23 #include <boost/math/distributions/detail/generic_quantile.hpp>
  24
  25 namespace boost
  26 {
  27    namespace math
  28    {
  29
  30       template <class RealType, class Policy>
  31       class non_central_chi_squared_distribution;
  32
  33       namespace detail{
  34
  35          template <class T, class Policy>
  36          T non_central_chi_square_q(T x, T f, T theta, const Policy& pol, T init_sum = 0)
  37          {
  38             //
  39             // Computes the complement of the Non-Central Chi-Square
  40             // Distribution CDF by summing a weighted sum of complements
  41             // of the central-distributions.  The weighting factor is
  42             // a Poisson Distribution.
  43             //
  44             // This is an application of the technique described in:
  45             //
  46             // Computing discrete mixtures of continuous
  47             // distributions: noncentral chisquare, noncentral t
  48             // and the distribution of the square of the sample
  49             // multiple correlation coeficient.
  50             // D. Benton, K. Krishnamoorthy.
  51             // Computational Statistics & Data Analysis 43 (2003) 249 - 267
  52             //
  53             BOOST_MATH_STD_USING
  54
  55             // Special case:
  56             if(x == 0)
  57                return 1;
  58
  59             //
  60             // Initialize the variables we'll be using:
  61             //
  62             T lambda = theta / 2;
  63             T del = f / 2;
  64             T y = x / 2;
  65             boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
  66             T errtol = boost::math::policies::get_epsilon<T, Policy>();
  67             T sum = init_sum;
  68             //
  69             // k is the starting location for iteration, we'll
  70             // move both forwards and backwards from this point.
  71             // k is chosen as the peek of the Poisson weights, which
  72             // will occur *before* the largest term.
  73             //
  74             int k = iround(lambda, pol);
  75             // Forwards and backwards Poisson weights:
  76             T poisf = boost::math::gamma_p_derivative(static_cast<T>(1 + k), lambda, pol);
  77             T poisb = poisf * k / lambda;
  78             // Initial forwards central chi squared term:
  79             T gamf = boost::math::gamma_q(del + k, y, pol);
  80             // Forwards and backwards recursion terms on the central chi squared:
  81             T xtermf = boost::math::gamma_p_derivative(del + 1 + k, y, pol);
  82             T xtermb = xtermf * (del + k) / y;
  83             // Initial backwards central chi squared term:
  84             T gamb = gamf - xtermb;
  85
  86             //
  87             // Forwards iteration first, this is the
  88             // stable direction for the gamma function
  89             // recurrences:
  90             //
  91             int i;
  92             for(i = k; static_cast<boost::uintmax_t>(i-k) < max_iter; ++i)
  93             {
  94                T term = poisf * gamf;
  95                sum += term;
  96                poisf *= lambda / (i + 1);
  97                gamf += xtermf;
  98                xtermf *= y / (del + i + 1);
  99                if(((sum == 0) || (fabs(term / sum) < errtol)) && (term >= poisf * gamf))
 100                   break;
 101             }
 102             //Error check:
 103             if(static_cast<boost::uintmax_t>(i-k) >= max_iter)
 104                return policies::raise_evaluation_error(
 105                   "cdf(non_central_chi_squared_distribution<%1%>, %1%)",
 106                   "Series did not converge, closest value was %1%", sum, pol);
 107             //
 108             // Now backwards iteration: the gamma
 109             // function recurrences are unstable in this
 110             // direction, we rely on the terms deminishing in size
 111             // faster than we introduce cancellation errors.
 112             // For this reason it's very important that we start
 113             // *before* the largest term so that backwards iteration
 114             // is strictly converging.
 115             //
 116             for(i = k - 1; i >= 0; --i)
 117             {
 118                T term = poisb * gamb;
 119                sum += term;
 120                poisb *= i / lambda;
 121                xtermb *= (del + i) / y;
 122                gamb -= xtermb;
 123                if((sum == 0) || (fabs(term / sum) < errtol))
 124                   break;
 125             }
 126
 127             return sum;
 128          }
 129
 130          template <class T, class Policy>
 131          T non_central_chi_square_p_ding(T x, T f, T theta, const Policy& pol, T init_sum = 0)
 132          {
 133             //
 134             // This is an implementation of:
 135             //
 136             // Algorithm AS 275:
 137             // Computing the Non-Central #2 Distribution Function
 138             // Cherng G. Ding
 139             // Applied Statistics, Vol. 41, No. 2. (1992), pp. 478-482.
 140             //
 141             // This uses a stable forward iteration to sum the
 142             // CDF, unfortunately this can not be used for large
 143             // values of the non-centrality parameter because:
 144             // * The first term may underfow to zero.
 145             // * We may need an extra-ordinary number of terms
 146             //   before we reach the first *significant* term.
 147             //
 148             BOOST_MATH_STD_USING
 149             // Special case:
 150             if(x == 0)
 151                return 0;
 152             T tk = boost::math::gamma_p_derivative(f/2 + 1, x/2, pol);
 153             T lambda = theta / 2;
 154             T vk = exp(-lambda);
 155             T uk = vk;
 156             T sum = init_sum + tk * vk;
 157             if(sum == 0)
 158                return sum;
 159
 160             boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
 161             T errtol = boost::math::policies::get_epsilon<T, Policy>();
 162
 163             int i;
 164             T lterm(0), term(0);
 165             for(i = 1; static_cast<boost::uintmax_t>(i) < max_iter; ++i)
 166             {
 167                tk = tk * x / (f + 2 * i);
 168                uk = uk * lambda / i;
 169                vk = vk + uk;
 170                lterm = term;
 171                term = vk * tk;
 172                sum += term;
 173                if((fabs(term / sum) < errtol) && (term <= lterm))
 174                   break;
 175             }
 176             //Error check:
 177             if(static_cast<boost::uintmax_t>(i) >= max_iter)
 178                return policies::raise_evaluation_error(
 179                   "cdf(non_central_chi_squared_distribution<%1%>, %1%)",
 180                   "Series did not converge, closest value was %1%", sum, pol);
 181             return sum;
 182          }
 183
 184
 185          template <class T, class Policy>
 186          T non_central_chi_square_p(T y, T n, T lambda, const Policy& pol, T init_sum)
 187          {
 188             //
 189             // This is taken more or less directly from:
 190             //
 191             // Computing discrete mixtures of continuous
 192             // distributions: noncentral chisquare, noncentral t
 193             // and the distribution of the square of the sample
 194             // multiple correlation coeficient.
 195             // D. Benton, K. Krishnamoorthy.
 196             // Computational Statistics & Data Analysis 43 (2003) 249 - 267
 197             //
 198             // We're summing a Poisson weighting term multiplied by
 199             // a central chi squared distribution.
 200             //
 201             BOOST_MATH_STD_USING
 202             // Special case:
 203             if(y == 0)
 204                return 0;
 205             boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
 206             T errtol = boost::math::policies::get_epsilon<T, Policy>();
 207             T errorf(0), errorb(0);
 208
 209             T x = y / 2;
 210             T del = lambda / 2;
 211             //
 212             // Starting location for the iteration, we'll iterate
 213             // both forwards and backwards from this point.  The
 214             // location chosen is the maximum of the Poisson weight
 215             // function, which ocurrs *after* the largest term in the
 216             // sum.
 217             //
 218             int k = iround(del, pol);
 219             T a = n / 2 + k;
 220             // Central chi squared term for forward iteration:
 221             T gamkf = boost::math::gamma_p(a, x, pol);
 222
 223             if(lambda == 0)
 224                return gamkf;
 225             // Central chi squared term for backward iteration:
 226             T gamkb = gamkf;
 227             // Forwards Poisson weight:
 228             T poiskf = gamma_p_derivative(static_cast<T>(k+1), del, pol);
 229             // Backwards Poisson weight:
 230             T poiskb = poiskf;
 231             // Forwards gamma function recursion term:
 232             T xtermf = boost::math::gamma_p_derivative(a, x, pol);
 233             // Backwards gamma function recursion term:
 234             T xtermb = xtermf * x / a;
 235             T sum = init_sum + poiskf * gamkf;
 236             if(sum == 0)
 237                return sum;
 238             int i = 1;
 239             //
 240             // Backwards recursion first, this is the stable
 241             // direction for gamma function recurrences:
 242             //
 243             while(i <= k)
 244             {
 245                xtermb *= (a - i + 1) / x;
 246                gamkb += xtermb;
 247                poiskb = poiskb * (k - i + 1) / del;
 248                errorf = errorb;
 249                errorb = gamkb * poiskb;
 250                sum += errorb;
 251                if((fabs(errorb / sum) < errtol) && (errorb <= errorf))
 252                   break;
 253                ++i;
 254             }
 255             i = 1;
 256             //
 257             // Now forwards recursion, the gamma function
 258             // recurrence relation is unstable in this direction,
 259             // so we rely on the magnitude of successive terms
 260             // decreasing faster than we introduce cancellation error.
 261             // For this reason it's vital that k is chosen to be *after*
 262             // the largest term, so that successive forward iterations
 263             // are strictly (and rapidly) converging.
 264             //
 265             do
 266             {
 267                xtermf = xtermf * x / (a + i - 1);
 268                gamkf = gamkf - xtermf;
 269                poiskf = poiskf * del / (k + i);
 270                errorf = poiskf * gamkf;
 271                sum += errorf;
 272                ++i;
 273             }while((fabs(errorf / sum) > errtol) && (static_cast<boost::uintmax_t>(i) < max_iter));
 274
 275             //Error check:
 276             if(static_cast<boost::uintmax_t>(i) >= max_iter)
 277                return policies::raise_evaluation_error(
 278                   "cdf(non_central_chi_squared_distribution<%1%>, %1%)",
 279                   "Series did not converge, closest value was %1%", sum, pol);
 280
 281             return sum;
 282          }
 283
 284          template <class T, class Policy>
 285          T non_central_chi_square_pdf(T x, T n, T lambda, const Policy& pol)
 286          {
 287             //
 288             // As above but for the PDF:
 289             //
 290             BOOST_MATH_STD_USING
 291             boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
 292             T errtol = boost::math::policies::get_epsilon<T, Policy>();
 293             T x2 = x / 2;
 294             T n2 = n / 2;
 295             T l2 = lambda / 2;
 296             T sum = 0;
 297             int k = itrunc(l2);
 298             T pois = gamma_p_derivative(static_cast<T>(k + 1), l2, pol) * gamma_p_derivative(static_cast<T>(n2 + k), x2);
 299             if(pois == 0)
 300                return 0;
 301             T poisb = pois;
 302             for(int i = k; ; ++i)
 303             {
 304                sum += pois;
 305                if(pois / sum < errtol)
 306                   break;
 307                if(static_cast<boost::uintmax_t>(i - k) >= max_iter)
 308                   return policies::raise_evaluation_error(
 309                      "pdf(non_central_chi_squared_distribution<%1%>, %1%)",
 310                      "Series did not converge, closest value was %1%", sum, pol);
 311                pois *= l2 * x2 / ((i + 1) * (n2 + i));
 312             }
 313             for(int i = k - 1; i >= 0; --i)
 314             {
 315                poisb *= (i + 1) * (n2 + i) / (l2 * x2);
 316                sum += poisb;
 317                if(poisb / sum < errtol)
 318                   break;
 319             }
 320             return sum / 2;
 321          }
 322
 323          template <class RealType, class Policy>
 324          inline RealType non_central_chi_squared_cdf(RealType x, RealType k, RealType l, bool invert, const Policy&)
 325          {
 326             typedef typename policies::evaluation<RealType, Policy>::type value_type;
 327             typedef typename policies::normalise<
 328                Policy,
 329                policies::promote_float<false>,
 330                policies::promote_double<false>,
 331                policies::discrete_quantile<>,
 332                policies::assert_undefined<> >::type forwarding_policy;
 333
 334             BOOST_MATH_STD_USING
 335             value_type result;
 336             if(l == 0)
 337               return invert == false ? cdf(boost::math::chi_squared_distribution<RealType, Policy>(k), x) : cdf(complement(boost::math::chi_squared_distribution<RealType, Policy>(k), x));
 338             else if(x > k + l)
 339             {
 340                // Complement is the smaller of the two:
 341                result = detail::non_central_chi_square_q(
 342                   static_cast<value_type>(x),
 343                   static_cast<value_type>(k),
 344                   static_cast<value_type>(l),
 345                   forwarding_policy(),
 346                   static_cast<value_type>(invert ? 0 : -1));
 347                invert = !invert;
 348             }
 349             else if(l < 200)
 350             {
 351                // For small values of the non-centrality parameter
 352                // we can use Ding's method:
 353                result = detail::non_central_chi_square_p_ding(
 354                   static_cast<value_type>(x),
 355                   static_cast<value_type>(k),
 356                   static_cast<value_type>(l),
 357                   forwarding_policy(),
 358                   static_cast<value_type>(invert ? -1 : 0));
 359             }
 360             else
 361             {
 362                // For largers values of the non-centrality
 363                // parameter Ding's method will consume an
 364                // extra-ordinary number of terms, and worse
 365                // may return zero when the result is in fact
 366                // finite, use Krishnamoorthy's method instead:
 367                result = detail::non_central_chi_square_p(
 368                   static_cast<value_type>(x),
 369                   static_cast<value_type>(k),
 370                   static_cast<value_type>(l),
 371                   forwarding_policy(),
 372                   static_cast<value_type>(invert ? -1 : 0));
 373             }
 374             if(invert)
 375                result = -result;
 376             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
 377                result,
 378                "boost::math::non_central_chi_squared_cdf<%1%>(%1%, %1%, %1%)");
 379          }
 380
 381          template <class T, class Policy>
 382          struct nccs_quantile_functor
 383          {
 384             nccs_quantile_functor(const non_central_chi_squared_distribution<T,Policy>& d, T t, bool c)
 385                : dist(d), target(t), comp(c) {}
 386
 387             T operator()(const T& x)
 388             {
 389                return comp ?
 390                   target - cdf(complement(dist, x))
 391                   : cdf(dist, x) - target;
 392             }
 393
 394          private:
 395             non_central_chi_squared_distribution<T,Policy> dist;
 396             T target;
 397             bool comp;
 398          };
 399
 400          template <class RealType, class Policy>
 401          RealType nccs_quantile(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& p, bool comp)
 402          {
 403             BOOST_MATH_STD_USING
 404             static const char* function = "quantile(non_central_chi_squared_distribution<%1%>, %1%)";
 405             typedef typename policies::evaluation<RealType, Policy>::type value_type;
 406             typedef typename policies::normalise<
 407                Policy,
 408                policies::promote_float<false>,
 409                policies::promote_double<false>,
 410                policies::discrete_quantile<>,
 411                policies::assert_undefined<> >::type forwarding_policy;
 412
 413             value_type k = dist.degrees_of_freedom();
 414             value_type l = dist.non_centrality();
 415             value_type r;
 416             if(!detail::check_df(
 417                function,
 418                k, &r, Policy())
 419                ||
 420             !detail::check_non_centrality(
 421                function,
 422                l,
 423                &r,
 424                Policy())
 425                ||
 426             !detail::check_probability(
 427                function,
 428                static_cast<value_type>(p),
 429                &r,
 430                Policy()))
 431                   return (RealType)r;
 432             //
 433             // Special cases get short-circuited first:
 434             //
 435             if(p == 0)
 436                return comp ? policies::raise_overflow_error<RealType>(function, 0, Policy()) : 0;
 437             if(p == 1)
 438                return comp ? 0 : policies::raise_overflow_error<RealType>(function, 0, Policy());
 439             //
 440             // This is Pearson's approximation to the quantile, see
 441             // Pearson, E. S. (1959) "Note on an approximation to the distribution of
 442             // noncentral chi squared", Biometrika 46: 364.
 443             // See also:
 444             // "A comparison of approximations to percentiles of the noncentral chi2-distribution",
 445             // Hardeo Sahai and Mario Miguel Ojeda, Revista de Matematica: Teoria y Aplicaciones 2003 10(1-2) : 57-76.
 446             // Note that the latter reference refers to an approximation of the CDF, when they really mean the quantile.
 447             //
 448             value_type b = -(l * l) / (k + 3 * l);
 449             value_type c = (k + 3 * l) / (k + 2 * l);
 450             value_type ff = (k + 2 * l) / (c * c);
 451             value_type guess;
 452             if(comp)
 453             {
 454                guess = b + c * quantile(complement(chi_squared_distribution<value_type, forwarding_policy>(ff), p));
 455             }
 456             else
 457             {
 458                guess = b + c * quantile(chi_squared_distribution<value_type, forwarding_policy>(ff), p);
 459             }
 460             //
 461             // Sometimes guess goes very small or negative, in that case we have
 462             // to do something else for the initial guess, this approximation
 463             // was provided in a private communication from Thomas Luu, PhD candidate,
 464             // University College London.  It's an asymptotic expansion for the
 465             // quantile which usually gets us within an order of magnitude of the
 466             // correct answer.
 467             // Fast and accurate parallel computation of quantile functions for random number generation,
 468             // Thomas LuuDoctorial Thesis 2016
 469             // http://discovery.ucl.ac.uk/1482128/
 470             //
 471             if(guess < 0.005)
 472             {
 473                value_type pp = comp ? 1 - p : p;
 474                //guess = pow(pow(value_type(2), (k / 2 - 1)) * exp(l / 2) * pp * k, 2 / k);
 475                guess = pow(pow(value_type(2), (k / 2 - 1)) * exp(l / 2) * pp * k * boost::math::tgamma(k / 2, forwarding_policy()), (2 / k));
 476                if(guess == 0)
 477                   guess = tools::min_value<value_type>();
 478             }
 479             value_type result = detail::generic_quantile(
 480                non_central_chi_squared_distribution<value_type, forwarding_policy>(k, l),
 481                p,
 482                guess,
 483                comp,
 484                function);
 485
 486             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
 487                result,
 488                function);
 489          }
 490
 491          template <class RealType, class Policy>
 492          RealType nccs_pdf(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
 493          {
 494             BOOST_MATH_STD_USING
 495             static const char* function = "pdf(non_central_chi_squared_distribution<%1%>, %1%)";
 496             typedef typename policies::evaluation<RealType, Policy>::type value_type;
 497             typedef typename policies::normalise<
 498                Policy,
 499                policies::promote_float<false>,
 500                policies::promote_double<false>,
 501                policies::discrete_quantile<>,
 502                policies::assert_undefined<> >::type forwarding_policy;
 503
 504             value_type k = dist.degrees_of_freedom();
 505             value_type l = dist.non_centrality();
 506             value_type r;
 507             if(!detail::check_df(
 508                function,
 509                k, &r, Policy())
 510                ||
 511             !detail::check_non_centrality(
 512                function,
 513                l,
 514                &r,
 515                Policy())
 516                ||
 517             !detail::check_positive_x(
 518                function,
 519                (value_type)x,
 520                &r,
 521                Policy()))
 522                   return (RealType)r;
 523
 524          if(l == 0)
 525             return pdf(boost::math::chi_squared_distribution<RealType, forwarding_policy>(dist.degrees_of_freedom()), x);
 526
 527          // Special case:
 528          if(x == 0)
 529             return 0;
 530          if(l > 50)
 531          {
 532             r = non_central_chi_square_pdf(static_cast<value_type>(x), k, l, forwarding_policy());
 533          }
 534          else
 535          {
 536             r = log(x / l) * (k / 4 - 0.5f) - (x + l) / 2;
 537             if(fabs(r) >= tools::log_max_value<RealType>() / 4)
 538             {
 539                r = non_central_chi_square_pdf(static_cast<value_type>(x), k, l, forwarding_policy());
 540             }
 541             else
 542             {
 543                r = exp(r);
 544                r = 0.5f * r
 545                   * boost::math::cyl_bessel_i(k/2 - 1, sqrt(l * x), forwarding_policy());
 546             }
 547          }
 548          return policies::checked_narrowing_cast<RealType, forwarding_policy>(
 549                r,
 550                function);
 551          }
 552
 553          template <class RealType, class Policy>
 554          struct degrees_of_freedom_finder
 555          {
 556             degrees_of_freedom_finder(
 557                RealType lam_, RealType x_, RealType p_, bool c)
 558                : lam(lam_), x(x_), p(p_), comp(c) {}
 559
 560             RealType operator()(const RealType& v)
 561             {
 562                non_central_chi_squared_distribution<RealType, Policy> d(v, lam);
 563                return comp ?
 564                   RealType(p - cdf(complement(d, x)))
 565                   : RealType(cdf(d, x) - p);
 566             }
 567          private:
 568             RealType lam;
 569             RealType x;
 570             RealType p;
 571             bool comp;
 572          };
 573
 574          template <class RealType, class Policy>
 575          inline RealType find_degrees_of_freedom(
 576             RealType lam, RealType x, RealType p, RealType q, const Policy& pol)
 577          {
 578             const char* function = "non_central_chi_squared<%1%>::find_degrees_of_freedom";
 579             if((p == 0) || (q == 0))
 580             {
 581                //
 582                // Can't a thing if one of p and q is zero:
 583                //
 584                return policies::raise_evaluation_error<RealType>(function,
 585                   "Can't find degrees of freedom when the probability is 0 or 1, only possible answer is %1%",
 586                   RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy());
 587             }
 588             degrees_of_freedom_finder<RealType, Policy> f(lam, x, p < q ? p : q, p < q ? false : true);
 589             tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
 590             boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
 591             //
 592             // Pick an initial guess that we know will give us a probability
 593             // right around 0.5.
 594             //
 595             RealType guess = x - lam;
 596             if(guess < 1)
 597                guess = 1;
 598             std::pair<RealType, RealType> ir = tools::bracket_and_solve_root(
 599                f, guess, RealType(2), false, tol, max_iter, pol);
 600             RealType result = ir.first + (ir.second - ir.first) / 2;
 601             if(max_iter >= policies::get_max_root_iterations<Policy>())
 602             {
 603                return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:"
 604                   " or there is no answer to problem.  Current best guess is %1%", result, Policy());
 605             }
 606             return result;
 607          }
 608
 609          template <class RealType, class Policy>
 610          struct non_centrality_finder
 611          {
 612             non_centrality_finder(
 613                RealType v_, RealType x_, RealType p_, bool c)
 614                : v(v_), x(x_), p(p_), comp(c) {}
 615
 616             RealType operator()(const RealType& lam)
 617             {
 618                non_central_chi_squared_distribution<RealType, Policy> d(v, lam);
 619                return comp ?
 620                   RealType(p - cdf(complement(d, x)))
 621                   : RealType(cdf(d, x) - p);
 622             }
 623          private:
 624             RealType v;
 625             RealType x;
 626             RealType p;
 627             bool comp;
 628          };
 629
 630          template <class RealType, class Policy>
 631          inline RealType find_non_centrality(
 632             RealType v, RealType x, RealType p, RealType q, const Policy& pol)
 633          {
 634             const char* function = "non_central_chi_squared<%1%>::find_non_centrality";
 635             if((p == 0) || (q == 0))
 636             {
 637                //
 638                // Can't do a thing if one of p and q is zero:
 639                //
 640                return policies::raise_evaluation_error<RealType>(function,
 641                   "Can't find non centrality parameter when the probability is 0 or 1, only possible answer is %1%",
 642                   RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy());
 643             }
 644             non_centrality_finder<RealType, Policy> f(v, x, p < q ? p : q, p < q ? false : true);
 645             tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
 646             boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
 647             //
 648             // Pick an initial guess that we know will give us a probability
 649             // right around 0.5.
 650             //
 651             RealType guess = x - v;
 652             if(guess < 1)
 653                guess = 1;
 654             std::pair<RealType, RealType> ir = tools::bracket_and_solve_root(
 655                f, guess, RealType(2), false, tol, max_iter, pol);
 656             RealType result = ir.first + (ir.second - ir.first) / 2;
 657             if(max_iter >= policies::get_max_root_iterations<Policy>())
 658             {
 659                return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:"
 660                   " or there is no answer to problem.  Current best guess is %1%", result, Policy());
 661             }
 662             return result;
 663          }
 664
 665       }
 666
 667       template <class RealType = double, class Policy = policies::policy<> >
 668       class non_central_chi_squared_distribution
 669       {
 670       public:
 671          typedef RealType value_type;
 672          typedef Policy policy_type;
 673
 674          non_central_chi_squared_distribution(RealType df_, RealType lambda) : df(df_), ncp(lambda)
 675          {
 676             const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::non_central_chi_squared_distribution(%1%,%1%)";
 677             RealType r;
 678             detail::check_df(
 679                function,
 680                df, &r, Policy());
 681             detail::check_non_centrality(
 682                function,
 683                ncp,
 684                &r,
 685                Policy());
 686          } // non_central_chi_squared_distribution constructor.
 687
 688          RealType degrees_of_freedom() const
 689          { // Private data getter function.
 690             return df;
 691          }
 692          RealType non_centrality() const
 693          { // Private data getter function.
 694             return ncp;
 695          }
 696          static RealType find_degrees_of_freedom(RealType lam, RealType x, RealType p)
 697          {
 698             const char* function = "non_central_chi_squared<%1%>::find_degrees_of_freedom";
 699             typedef typename policies::evaluation<RealType, Policy>::type eval_type;
 700             typedef typename policies::normalise<
 701                Policy,
 702                policies::promote_float<false>,
 703                policies::promote_double<false>,
 704                policies::discrete_quantile<>,
 705                policies::assert_undefined<> >::type forwarding_policy;
 706             eval_type result = detail::find_degrees_of_freedom(
 707                static_cast<eval_type>(lam),
 708                static_cast<eval_type>(x),
 709                static_cast<eval_type>(p),
 710                static_cast<eval_type>(1-p),
 711                forwarding_policy());
 712             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
 713                result,
 714                function);
 715          }
 716          template <class A, class B, class C>
 717          static RealType find_degrees_of_freedom(const complemented3_type<A,B,C>& c)
 718          {
 719             const char* function = "non_central_chi_squared<%1%>::find_degrees_of_freedom";
 720             typedef typename policies::evaluation<RealType, Policy>::type eval_type;
 721             typedef typename policies::normalise<
 722                Policy,
 723                policies::promote_float<false>,
 724                policies::promote_double<false>,
 725                policies::discrete_quantile<>,
 726                policies::assert_undefined<> >::type forwarding_policy;
 727             eval_type result = detail::find_degrees_of_freedom(
 728                static_cast<eval_type>(c.dist),
 729                static_cast<eval_type>(c.param1),
 730                static_cast<eval_type>(1-c.param2),
 731                static_cast<eval_type>(c.param2),
 732                forwarding_policy());
 733             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
 734                result,
 735                function);
 736          }
 737          static RealType find_non_centrality(RealType v, RealType x, RealType p)
 738          {
 739             const char* function = "non_central_chi_squared<%1%>::find_non_centrality";
 740             typedef typename policies::evaluation<RealType, Policy>::type eval_type;
 741             typedef typename policies::normalise<
 742                Policy,
 743                policies::promote_float<false>,
 744                policies::promote_double<false>,
 745                policies::discrete_quantile<>,
 746                policies::assert_undefined<> >::type forwarding_policy;
 747             eval_type result = detail::find_non_centrality(
 748                static_cast<eval_type>(v),
 749                static_cast<eval_type>(x),
 750                static_cast<eval_type>(p),
 751                static_cast<eval_type>(1-p),
 752                forwarding_policy());
 753             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
 754                result,
 755                function);
 756          }
 757          template <class A, class B, class C>
 758          static RealType find_non_centrality(const complemented3_type<A,B,C>& c)
 759          {
 760             const char* function = "non_central_chi_squared<%1%>::find_non_centrality";
 761             typedef typename policies::evaluation<RealType, Policy>::type eval_type;
 762             typedef typename policies::normalise<
 763                Policy,
 764                policies::promote_float<false>,
 765                policies::promote_double<false>,
 766                policies::discrete_quantile<>,
 767                policies::assert_undefined<> >::type forwarding_policy;
 768             eval_type result = detail::find_non_centrality(
 769                static_cast<eval_type>(c.dist),
 770                static_cast<eval_type>(c.param1),
 771                static_cast<eval_type>(1-c.param2),
 772                static_cast<eval_type>(c.param2),
 773                forwarding_policy());
 774             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
 775                result,
 776                function);
 777          }
 778       private:
 779          // Data member, initialized by constructor.
 780          RealType df; // degrees of freedom.
 781          RealType ncp; // non-centrality parameter
 782       }; // template <class RealType, class Policy> class non_central_chi_squared_distribution
 783
 784       typedef non_central_chi_squared_distribution<double> non_central_chi_squared; // Reserved name of type double.
 785
 786       // Non-member functions to give properties of the distribution.
 787
 788       template <class RealType, class Policy>
 789       inline const std::pair<RealType, RealType> range(const non_central_chi_squared_distribution<RealType, Policy>& /* dist */)
 790       { // Range of permissible values for random variable k.
 791          using boost::math::tools::max_value;
 792          return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // Max integer?
 793       }
 794
 795       template <class RealType, class Policy>
 796       inline const std::pair<RealType, RealType> support(const non_central_chi_squared_distribution<RealType, Policy>& /* dist */)
 797       { // Range of supported values for random variable k.
 798          // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
 799          using boost::math::tools::max_value;
 800          return std::pair<RealType, RealType>(static_cast<RealType>(0),  max_value<RealType>());
 801       }
 802
 803       template <class RealType, class Policy>
 804       inline RealType mean(const non_central_chi_squared_distribution<RealType, Policy>& dist)
 805       { // Mean of poisson distribution = lambda.
 806          const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::mean()";
 807          RealType k = dist.degrees_of_freedom();
 808          RealType l = dist.non_centrality();
 809          RealType r;
 810          if(!detail::check_df(
 811             function,
 812             k, &r, Policy())
 813             ||
 814          !detail::check_non_centrality(
 815             function,
 816             l,
 817             &r,
 818             Policy()))
 819                return r;
 820          return k + l;
 821       } // mean
 822
 823       template <class RealType, class Policy>
 824       inline RealType mode(const non_central_chi_squared_distribution<RealType, Policy>& dist)
 825       { // mode.
 826          static const char* function = "mode(non_central_chi_squared_distribution<%1%> const&)";
 827
 828          RealType k = dist.degrees_of_freedom();
 829          RealType l = dist.non_centrality();
 830          RealType r;
 831          if(!detail::check_df(
 832             function,
 833             k, &r, Policy())
 834             ||
 835          !detail::check_non_centrality(
 836             function,
 837             l,
 838             &r,
 839             Policy()))
 840                return (RealType)r;
 841          return detail::generic_find_mode(dist, 1 + k, function);
 842       }
 843
 844       template <class RealType, class Policy>
 845       inline RealType variance(const non_central_chi_squared_distribution<RealType, Policy>& dist)
 846       { // variance.
 847          const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::variance()";
 848          RealType k = dist.degrees_of_freedom();
 849          RealType l = dist.non_centrality();
 850          RealType r;
 851          if(!detail::check_df(
 852             function,
 853             k, &r, Policy())
 854             ||
 855          !detail::check_non_centrality(
 856             function,
 857             l,
 858             &r,
 859             Policy()))
 860                return r;
 861          return 2 * (2 * l + k);
 862       }
 863
 864       // RealType standard_deviation(const non_central_chi_squared_distribution<RealType, Policy>& dist)
 865       // standard_deviation provided by derived accessors.
 866
 867       template <class RealType, class Policy>
 868       inline RealType skewness(const non_central_chi_squared_distribution<RealType, Policy>& dist)
 869       { // skewness = sqrt(l).
 870          const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::skewness()";
 871          RealType k = dist.degrees_of_freedom();
 872          RealType l = dist.non_centrality();
 873          RealType r;
 874          if(!detail::check_df(
 875             function,
 876             k, &r, Policy())
 877             ||
 878          !detail::check_non_centrality(
 879             function,
 880             l,
 881             &r,
 882             Policy()))
 883                return r;
 884          BOOST_MATH_STD_USING
 885             return pow(2 / (k + 2 * l), RealType(3)/2) * (k + 3 * l);
 886       }
 887
 888       template <class RealType, class Policy>
 889       inline RealType kurtosis_excess(const non_central_chi_squared_distribution<RealType, Policy>& dist)
 890       {
 891          const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::kurtosis_excess()";
 892          RealType k = dist.degrees_of_freedom();
 893          RealType l = dist.non_centrality();
 894          RealType r;
 895          if(!detail::check_df(
 896             function,
 897             k, &r, Policy())
 898             ||
 899          !detail::check_non_centrality(
 900             function,
 901             l,
 902             &r,
 903             Policy()))
 904                return r;
 905          return 12 * (k + 4 * l) / ((k + 2 * l) * (k + 2 * l));
 906       } // kurtosis_excess
 907
 908       template <class RealType, class Policy>
 909       inline RealType kurtosis(const non_central_chi_squared_distribution<RealType, Policy>& dist)
 910       {
 911          return kurtosis_excess(dist) + 3;
 912       }
 913
 914       template <class RealType, class Policy>
 915       inline RealType pdf(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
 916       { // Probability Density/Mass Function.
 917          return detail::nccs_pdf(dist, x);
 918       } // pdf
 919
 920       template <class RealType, class Policy>
 921       RealType cdf(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
 922       {
 923          const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::cdf(%1%)";
 924          RealType k = dist.degrees_of_freedom();
 925          RealType l = dist.non_centrality();
 926          RealType r;
 927          if(!detail::check_df(
 928             function,
 929             k, &r, Policy())
 930             ||
 931          !detail::check_non_centrality(
 932             function,
 933             l,
 934             &r,
 935             Policy())
 936             ||
 937          !detail::check_positive_x(
 938             function,
 939             x,
 940             &r,
 941             Policy()))
 942                return r;
 943
 944          return detail::non_central_chi_squared_cdf(x, k, l, false, Policy());
 945       } // cdf
 946
 947       template <class RealType, class Policy>
 948       RealType cdf(const complemented2_type<non_central_chi_squared_distribution<RealType, Policy>, RealType>& c)
 949       { // Complemented Cumulative Distribution Function
 950          const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::cdf(%1%)";
 951          non_central_chi_squared_distribution<RealType, Policy> const& dist = c.dist;
 952          RealType x = c.param;
 953          RealType k = dist.degrees_of_freedom();
 954          RealType l = dist.non_centrality();
 955          RealType r;
 956          if(!detail::check_df(
 957             function,
 958             k, &r, Policy())
 959             ||
 960          !detail::check_non_centrality(
 961             function,
 962             l,
 963             &r,
 964             Policy())
 965             ||
 966          !detail::check_positive_x(
 967             function,
 968             x,
 969             &r,
 970             Policy()))
 971                return r;
 972
 973          return detail::non_central_chi_squared_cdf(x, k, l, true, Policy());
 974       } // ccdf
 975
 976       template <class RealType, class Policy>
 977       inline RealType quantile(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& p)
 978       { // Quantile (or Percent Point) function.
 979          return detail::nccs_quantile(dist, p, false);
 980       } // quantile
 981
 982       template <class RealType, class Policy>
 983       inline RealType quantile(const complemented2_type<non_central_chi_squared_distribution<RealType, Policy>, RealType>& c)
 984       { // Quantile (or Percent Point) function.
 985          return detail::nccs_quantile(c.dist, c.param, true);
 986       } // quantile complement.
 987
 988    } // namespace math
 989 } // namespace boost
 990
 991 // This include must be at the end, *after* the accessors
 992 // for this distribution have been defined, in order to
 993 // keep compilers that support two-phase lookup happy.
 994 #include <boost/math/distributions/detail/derived_accessors.hpp>
 995
 996 #endif // BOOST_MATH_SPECIAL_NON_CENTRAL_CHI_SQUARE_HPP
 997
 998
 999