math/distributions/non_central_chi_squared.hpp

0001 // boost\math\distributions\non_central_chi_squared.hpp
0002
0003 // Copyright John Maddock 2008.
0004
0005 // Use, modification and distribution are subject to the
0006 // Boost Software License, Version 1.0.
0007 // (See accompanying file LICENSE_1_0.txt
0008 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0009
0010 #ifndef BOOST_MATH_SPECIAL_NON_CENTRAL_CHI_SQUARE_HPP
0011 #define BOOST_MATH_SPECIAL_NON_CENTRAL_CHI_SQUARE_HPP
0012
0013 #include <boost/math/distributions/fwd.hpp>
0014 #include <boost/math/special_functions/gamma.hpp> // for incomplete gamma. gamma_q
0015 #include <boost/math/special_functions/bessel.hpp> // for cyl_bessel_i
0016 #include <boost/math/special_functions/round.hpp> // for llround
0017 #include <boost/math/distributions/complement.hpp> // complements
0018 #include <boost/math/distributions/chi_squared.hpp> // central distribution
0019 #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks
0020 #include <boost/math/special_functions/fpclassify.hpp> // isnan.
0021 #include <boost/math/tools/roots.hpp> // for root finding.
0022 #include <boost/math/distributions/detail/generic_mode.hpp>
0023 #include <boost/math/distributions/detail/generic_quantile.hpp>
0024
0025 namespace boost
0026 {
0027    namespace math
0028    {
0029
0030       template <class RealType, class Policy>
0031       class non_central_chi_squared_distribution;
0032
0033       namespace detail{
0034
0035          template <class T, class Policy>
0036          T non_central_chi_square_q(T x, T f, T theta, const Policy& pol, T init_sum = 0)
0037          {
0038             //
0039             // Computes the complement of the Non-Central Chi-Square
0040             // Distribution CDF by summing a weighted sum of complements
0041             // of the central-distributions.  The weighting factor is
0042             // a Poisson Distribution.
0043             //
0044             // This is an application of the technique described in:
0045             //
0046             // Computing discrete mixtures of continuous
0047             // distributions: noncentral chisquare, noncentral t
0048             // and the distribution of the square of the sample
0049             // multiple correlation coefficient.
0050             // D. Benton, K. Krishnamoorthy.
0051             // Computational Statistics & Data Analysis 43 (2003) 249 - 267
0052             //
0053             BOOST_MATH_STD_USING
0054
0055             // Special case:
0056             if(x == 0)
0057                return 1;
0058
0059             //
0060             // Initialize the variables we'll be using:
0061             //
0062             T lambda = theta / 2;
0063             T del = f / 2;
0064             T y = x / 2;
0065             std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0066             T errtol = boost::math::policies::get_epsilon<T, Policy>();
0067             T sum = init_sum;
0068             //
0069             // k is the starting location for iteration, we'll
0070             // move both forwards and backwards from this point.
0071             // k is chosen as the peek of the Poisson weights, which
0072             // will occur *before* the largest term.
0073             //
0074             long long k = llround(lambda, pol);
0075             // Forwards and backwards Poisson weights:
0076             T poisf = boost::math::gamma_p_derivative(static_cast<T>(1 + k), lambda, pol);
0077             T poisb = poisf * k / lambda;
0078             // Initial forwards central chi squared term:
0079             T gamf = boost::math::gamma_q(del + k, y, pol);
0080             // Forwards and backwards recursion terms on the central chi squared:
0081             T xtermf = boost::math::gamma_p_derivative(del + 1 + k, y, pol);
0082             T xtermb = xtermf * (del + k) / y;
0083             // Initial backwards central chi squared term:
0084             T gamb = gamf - xtermb;
0085
0086             //
0087             // Forwards iteration first, this is the
0088             // stable direction for the gamma function
0089             // recurrences:
0090             //
0091             long long i;
0092             for(i = k; static_cast<std::uintmax_t>(i-k) < max_iter; ++i)
0093             {
0094                T term = poisf * gamf;
0095                sum += term;
0096                poisf *= lambda / (i + 1);
0097                gamf += xtermf;
0098                xtermf *= y / (del + i + 1);
0099                if(((sum == 0) || (fabs(term / sum) < errtol)) && (term >= poisf * gamf))
0100                   break;
0101             }
0102             //Error check:
0103             if(static_cast<std::uintmax_t>(i-k) >= max_iter)
0104                return policies::raise_evaluation_error(
0105                   "cdf(non_central_chi_squared_distribution<%1%>, %1%)",
0106                   "Series did not converge, closest value was %1%", sum, pol);
0107             //
0108             // Now backwards iteration: the gamma
0109             // function recurrences are unstable in this
0110             // direction, we rely on the terms diminishing in size
0111             // faster than we introduce cancellation errors.
0112             // For this reason it's very important that we start
0113             // *before* the largest term so that backwards iteration
0114             // is strictly converging.
0115             //
0116             for(i = k - 1; i >= 0; --i)
0117             {
0118                T term = poisb * gamb;
0119                sum += term;
0120                poisb *= i / lambda;
0121                xtermb *= (del + i) / y;
0122                gamb -= xtermb;
0123                if((sum == 0) || (fabs(term / sum) < errtol))
0124                   break;
0125             }
0126
0127             return sum;
0128          }
0129
0130          template <class T, class Policy>
0131          T non_central_chi_square_p_ding(T x, T f, T theta, const Policy& pol, T init_sum = 0)
0132          {
0133             //
0134             // This is an implementation of:
0135             //
0136             // Algorithm AS 275:
0137             // Computing the Non-Central #2 Distribution Function
0138             // Cherng G. Ding
0139             // Applied Statistics, Vol. 41, No. 2. (1992), pp. 478-482.
0140             //
0141             // This uses a stable forward iteration to sum the
0142             // CDF, unfortunately this can not be used for large
0143             // values of the non-centrality parameter because:
0144             // * The first term may underflow to zero.
0145             // * We may need an extra-ordinary number of terms
0146             //   before we reach the first *significant* term.
0147             //
0148             BOOST_MATH_STD_USING
0149             // Special case:
0150             if(x == 0)
0151                return 0;
0152             T tk = boost::math::gamma_p_derivative(f/2 + 1, x/2, pol);
0153             T lambda = theta / 2;
0154             T vk = exp(-lambda);
0155             T uk = vk;
0156             T sum = init_sum + tk * vk;
0157             if(sum == 0)
0158                return sum;
0159
0160             std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0161             T errtol = boost::math::policies::get_epsilon<T, Policy>();
0162
0163             int i;
0164             T lterm(0), term(0);
0165             for(i = 1; static_cast<std::uintmax_t>(i) < max_iter; ++i)
0166             {
0167                tk = tk * x / (f + 2 * i);
0168                uk = uk * lambda / i;
0169                vk = vk + uk;
0170                lterm = term;
0171                term = vk * tk;
0172                sum += term;
0173                if((fabs(term / sum) < errtol) && (term <= lterm))
0174                   break;
0175             }
0176             //Error check:
0177             if(static_cast<std::uintmax_t>(i) >= max_iter)
0178                return policies::raise_evaluation_error(
0179                   "cdf(non_central_chi_squared_distribution<%1%>, %1%)",
0180                   "Series did not converge, closest value was %1%", sum, pol);
0181             return sum;
0182          }
0183
0184
0185          template <class T, class Policy>
0186          T non_central_chi_square_p(T y, T n, T lambda, const Policy& pol, T init_sum)
0187          {
0188             //
0189             // This is taken more or less directly from:
0190             //
0191             // Computing discrete mixtures of continuous
0192             // distributions: noncentral chisquare, noncentral t
0193             // and the distribution of the square of the sample
0194             // multiple correlation coefficient.
0195             // D. Benton, K. Krishnamoorthy.
0196             // Computational Statistics & Data Analysis 43 (2003) 249 - 267
0197             //
0198             // We're summing a Poisson weighting term multiplied by
0199             // a central chi squared distribution.
0200             //
0201             BOOST_MATH_STD_USING
0202             // Special case:
0203             if(y == 0)
0204                return 0;
0205             std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0206             T errtol = boost::math::policies::get_epsilon<T, Policy>();
0207             T errorf(0), errorb(0);
0208
0209             T x = y / 2;
0210             T del = lambda / 2;
0211             //
0212             // Starting location for the iteration, we'll iterate
0213             // both forwards and backwards from this point.  The
0214             // location chosen is the maximum of the Poisson weight
0215             // function, which ocurrs *after* the largest term in the
0216             // sum.
0217             //
0218             long long k = llround(del, pol);
0219             T a = n / 2 + k;
0220             // Central chi squared term for forward iteration:
0221             T gamkf = boost::math::gamma_p(a, x, pol);
0222
0223             if(lambda == 0)
0224                return gamkf;
0225             // Central chi squared term for backward iteration:
0226             T gamkb = gamkf;
0227             // Forwards Poisson weight:
0228             T poiskf = gamma_p_derivative(static_cast<T>(k+1), del, pol);
0229             // Backwards Poisson weight:
0230             T poiskb = poiskf;
0231             // Forwards gamma function recursion term:
0232             T xtermf = boost::math::gamma_p_derivative(a, x, pol);
0233             // Backwards gamma function recursion term:
0234             T xtermb = xtermf * x / a;
0235             T sum = init_sum + poiskf * gamkf;
0236             if(sum == 0)
0237                return sum;
0238             int i = 1;
0239             //
0240             // Backwards recursion first, this is the stable
0241             // direction for gamma function recurrences:
0242             //
0243             while(i <= k)
0244             {
0245                xtermb *= (a - i + 1) / x;
0246                gamkb += xtermb;
0247                poiskb = poiskb * (k - i + 1) / del;
0248                errorf = errorb;
0249                errorb = gamkb * poiskb;
0250                sum += errorb;
0251                if((fabs(errorb / sum) < errtol) && (errorb <= errorf))
0252                   break;
0253                ++i;
0254             }
0255             i = 1;
0256             //
0257             // Now forwards recursion, the gamma function
0258             // recurrence relation is unstable in this direction,
0259             // so we rely on the magnitude of successive terms
0260             // decreasing faster than we introduce cancellation error.
0261             // For this reason it's vital that k is chosen to be *after*
0262             // the largest term, so that successive forward iterations
0263             // are strictly (and rapidly) converging.
0264             //
0265             do
0266             {
0267                xtermf = xtermf * x / (a + i - 1);
0268                gamkf = gamkf - xtermf;
0269                poiskf = poiskf * del / (k + i);
0270                errorf = poiskf * gamkf;
0271                sum += errorf;
0272                ++i;
0273             }while((fabs(errorf / sum) > errtol) && (static_cast<std::uintmax_t>(i) < max_iter));
0274
0275             //Error check:
0276             if(static_cast<std::uintmax_t>(i) >= max_iter)
0277                return policies::raise_evaluation_error(
0278                   "cdf(non_central_chi_squared_distribution<%1%>, %1%)",
0279                   "Series did not converge, closest value was %1%", sum, pol);
0280
0281             return sum;
0282          }
0283
0284          template <class T, class Policy>
0285          T non_central_chi_square_pdf(T x, T n, T lambda, const Policy& pol)
0286          {
0287             //
0288             // As above but for the PDF:
0289             //
0290             BOOST_MATH_STD_USING
0291             std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0292             T errtol = boost::math::policies::get_epsilon<T, Policy>();
0293             T x2 = x / 2;
0294             T n2 = n / 2;
0295             T l2 = lambda / 2;
0296             T sum = 0;
0297             long long k = lltrunc(l2);
0298             T pois = gamma_p_derivative(static_cast<T>(k + 1), l2, pol) * gamma_p_derivative(static_cast<T>(n2 + k), x2);
0299             if(pois == 0)
0300                return 0;
0301             T poisb = pois;
0302             for(long long i = k; ; ++i)
0303             {
0304                sum += pois;
0305                if(pois / sum < errtol)
0306                   break;
0307                if(static_cast<std::uintmax_t>(i - k) >= max_iter)
0308                   return policies::raise_evaluation_error(
0309                      "pdf(non_central_chi_squared_distribution<%1%>, %1%)",
0310                      "Series did not converge, closest value was %1%", sum, pol);
0311                pois *= l2 * x2 / ((i + 1) * (n2 + i));
0312             }
0313             for(long long i = k - 1; i >= 0; --i)
0314             {
0315                poisb *= (i + 1) * (n2 + i) / (l2 * x2);
0316                sum += poisb;
0317                if(poisb / sum < errtol)
0318                   break;
0319             }
0320             return sum / 2;
0321          }
0322
0323          template <class RealType, class Policy>
0324          inline RealType non_central_chi_squared_cdf(RealType x, RealType k, RealType l, bool invert, const Policy&)
0325          {
0326             typedef typename policies::evaluation<RealType, Policy>::type value_type;
0327             typedef typename policies::normalise<
0328                Policy,
0329                policies::promote_float<false>,
0330                policies::promote_double<false>,
0331                policies::discrete_quantile<>,
0332                policies::assert_undefined<> >::type forwarding_policy;
0333
0334             BOOST_MATH_STD_USING
0335             value_type result;
0336             if(l == 0)
0337               return invert == false ? cdf(boost::math::chi_squared_distribution<RealType, Policy>(k), x) : cdf(complement(boost::math::chi_squared_distribution<RealType, Policy>(k), x));
0338             else if(x > k + l)
0339             {
0340                // Complement is the smaller of the two:
0341                result = detail::non_central_chi_square_q(
0342                   static_cast<value_type>(x),
0343                   static_cast<value_type>(k),
0344                   static_cast<value_type>(l),
0345                   forwarding_policy(),
0346                   static_cast<value_type>(invert ? 0 : -1));
0347                invert = !invert;
0348             }
0349             else if(l < 200)
0350             {
0351                // For small values of the non-centrality parameter
0352                // we can use Ding's method:
0353                result = detail::non_central_chi_square_p_ding(
0354                   static_cast<value_type>(x),
0355                   static_cast<value_type>(k),
0356                   static_cast<value_type>(l),
0357                   forwarding_policy(),
0358                   static_cast<value_type>(invert ? -1 : 0));
0359             }
0360             else
0361             {
0362                // For largers values of the non-centrality
0363                // parameter Ding's method will consume an
0364                // extra-ordinary number of terms, and worse
0365                // may return zero when the result is in fact
0366                // finite, use Krishnamoorthy's method instead:
0367                result = detail::non_central_chi_square_p(
0368                   static_cast<value_type>(x),
0369                   static_cast<value_type>(k),
0370                   static_cast<value_type>(l),
0371                   forwarding_policy(),
0372                   static_cast<value_type>(invert ? -1 : 0));
0373             }
0374             if(invert)
0375                result = -result;
0376             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
0377                result,
0378                "boost::math::non_central_chi_squared_cdf<%1%>(%1%, %1%, %1%)");
0379          }
0380
0381          template <class T, class Policy>
0382          struct nccs_quantile_functor
0383          {
0384             nccs_quantile_functor(const non_central_chi_squared_distribution<T,Policy>& d, T t, bool c)
0385                : dist(d), target(t), comp(c) {}
0386
0387             T operator()(const T& x)
0388             {
0389                return comp ?
0390                   target - cdf(complement(dist, x))
0391                   : cdf(dist, x) - target;
0392             }
0393
0394          private:
0395             non_central_chi_squared_distribution<T,Policy> dist;
0396             T target;
0397             bool comp;
0398          };
0399
0400          template <class RealType, class Policy>
0401          RealType nccs_quantile(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& p, bool comp)
0402          {
0403             BOOST_MATH_STD_USING
0404             static const char* function = "quantile(non_central_chi_squared_distribution<%1%>, %1%)";
0405             typedef typename policies::evaluation<RealType, Policy>::type value_type;
0406             typedef typename policies::normalise<
0407                Policy,
0408                policies::promote_float<false>,
0409                policies::promote_double<false>,
0410                policies::discrete_quantile<>,
0411                policies::assert_undefined<> >::type forwarding_policy;
0412
0413             value_type k = dist.degrees_of_freedom();
0414             value_type l = dist.non_centrality();
0415             value_type r;
0416             if(!detail::check_df(
0417                function,
0418                k, &r, Policy())
0419                ||
0420             !detail::check_non_centrality(
0421                function,
0422                l,
0423                &r,
0424                Policy())
0425                ||
0426             !detail::check_probability(
0427                function,
0428                static_cast<value_type>(p),
0429                &r,
0430                Policy()))
0431                   return static_cast<RealType>(r);
0432             //
0433             // Special cases get short-circuited first:
0434             //
0435             if(p == 0)
0436                return comp ? policies::raise_overflow_error<RealType>(function, 0, Policy()) : 0;
0437             if(p == 1)
0438                return comp ? 0 : policies::raise_overflow_error<RealType>(function, 0, Policy());
0439             //
0440             // This is Pearson's approximation to the quantile, see
0441             // Pearson, E. S. (1959) "Note on an approximation to the distribution of
0442             // noncentral chi squared", Biometrika 46: 364.
0443             // See also:
0444             // "A comparison of approximations to percentiles of the noncentral chi2-distribution",
0445             // Hardeo Sahai and Mario Miguel Ojeda, Revista de Matematica: Teoria y Aplicaciones 2003 10(1-2) : 57-76.
0446             // Note that the latter reference refers to an approximation of the CDF, when they really mean the quantile.
0447             //
0448             value_type b = -(l * l) / (k + 3 * l);
0449             value_type c = (k + 3 * l) / (k + 2 * l);
0450             value_type ff = (k + 2 * l) / (c * c);
0451             value_type guess;
0452             if(comp)
0453             {
0454                guess = b + c * quantile(complement(chi_squared_distribution<value_type, forwarding_policy>(ff), p));
0455             }
0456             else
0457             {
0458                guess = b + c * quantile(chi_squared_distribution<value_type, forwarding_policy>(ff), p);
0459             }
0460             //
0461             // Sometimes guess goes very small or negative, in that case we have
0462             // to do something else for the initial guess, this approximation
0463             // was provided in a private communication from Thomas Luu, PhD candidate,
0464             // University College London.  It's an asymptotic expansion for the
0465             // quantile which usually gets us within an order of magnitude of the
0466             // correct answer.
0467             // Fast and accurate parallel computation of quantile functions for random number generation,
0468             // Thomas LuuDoctorial Thesis 2016
0469             // http://discovery.ucl.ac.uk/1482128/
0470             //
0471             if(guess < 0.005)
0472             {
0473                value_type pp = comp ? 1 - p : p;
0474                //guess = pow(pow(value_type(2), (k / 2 - 1)) * exp(l / 2) * pp * k, 2 / k);
0475                guess = pow(pow(value_type(2), (k / 2 - 1)) * exp(l / 2) * pp * k * boost::math::tgamma(k / 2, forwarding_policy()), (2 / k));
0476                if(guess == 0)
0477                   guess = tools::min_value<value_type>();
0478             }
0479             value_type result = detail::generic_quantile(
0480                non_central_chi_squared_distribution<value_type, forwarding_policy>(k, l),
0481                p,
0482                guess,
0483                comp,
0484                function);
0485
0486             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
0487                result,
0488                function);
0489          }
0490
0491          template <class RealType, class Policy>
0492          RealType nccs_pdf(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
0493          {
0494             BOOST_MATH_STD_USING
0495             static const char* function = "pdf(non_central_chi_squared_distribution<%1%>, %1%)";
0496             typedef typename policies::evaluation<RealType, Policy>::type value_type;
0497             typedef typename policies::normalise<
0498                Policy,
0499                policies::promote_float<false>,
0500                policies::promote_double<false>,
0501                policies::discrete_quantile<>,
0502                policies::assert_undefined<> >::type forwarding_policy;
0503
0504             value_type k = dist.degrees_of_freedom();
0505             value_type l = dist.non_centrality();
0506             value_type r;
0507             if(!detail::check_df(
0508                function,
0509                k, &r, Policy())
0510                ||
0511             !detail::check_non_centrality(
0512                function,
0513                l,
0514                &r,
0515                Policy())
0516                ||
0517             !detail::check_positive_x(
0518                function,
0519                (value_type)x,
0520                &r,
0521                Policy()))
0522                   return static_cast<RealType>(r);
0523
0524          if(l == 0)
0525             return pdf(boost::math::chi_squared_distribution<RealType, forwarding_policy>(dist.degrees_of_freedom()), x);
0526
0527          // Special case:
0528          if(x == 0)
0529             return 0;
0530          if(l > 50)
0531          {
0532             r = non_central_chi_square_pdf(static_cast<value_type>(x), k, l, forwarding_policy());
0533          }
0534          else
0535          {
0536             r = log(x / l) * (k / 4 - 0.5f) - (x + l) / 2;
0537             if(fabs(r) >= tools::log_max_value<RealType>() / 4)
0538             {
0539                r = non_central_chi_square_pdf(static_cast<value_type>(x), k, l, forwarding_policy());
0540             }
0541             else
0542             {
0543                r = exp(r);
0544                r = 0.5f * r
0545                   * boost::math::cyl_bessel_i(k/2 - 1, sqrt(l * x), forwarding_policy());
0546             }
0547          }
0548          return policies::checked_narrowing_cast<RealType, forwarding_policy>(
0549                r,
0550                function);
0551          }
0552
0553          template <class RealType, class Policy>
0554          struct degrees_of_freedom_finder
0555          {
0556             degrees_of_freedom_finder(
0557                RealType lam_, RealType x_, RealType p_, bool c)
0558                : lam(lam_), x(x_), p(p_), comp(c) {}
0559
0560             RealType operator()(const RealType& v)
0561             {
0562                non_central_chi_squared_distribution<RealType, Policy> d(v, lam);
0563                return comp ?
0564                   RealType(p - cdf(complement(d, x)))
0565                   : RealType(cdf(d, x) - p);
0566             }
0567          private:
0568             RealType lam;
0569             RealType x;
0570             RealType p;
0571             bool comp;
0572          };
0573
0574          template <class RealType, class Policy>
0575          inline RealType find_degrees_of_freedom(
0576             RealType lam, RealType x, RealType p, RealType q, const Policy& pol)
0577          {
0578             const char* function = "non_central_chi_squared<%1%>::find_degrees_of_freedom";
0579             if((p == 0) || (q == 0))
0580             {
0581                //
0582                // Can't a thing if one of p and q is zero:
0583                //
0584                return policies::raise_evaluation_error<RealType>(function,
0585                   "Can't find degrees of freedom when the probability is 0 or 1, only possible answer is %1%",
0586                   RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy());
0587             }
0588             degrees_of_freedom_finder<RealType, Policy> f(lam, x, p < q ? p : q, p < q ? false : true);
0589             tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
0590             std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
0591             //
0592             // Pick an initial guess that we know will give us a probability
0593             // right around 0.5.
0594             //
0595             RealType guess = x - lam;
0596             if(guess < 1)
0597                guess = 1;
0598             std::pair<RealType, RealType> ir = tools::bracket_and_solve_root(
0599                f, guess, RealType(2), false, tol, max_iter, pol);
0600             RealType result = ir.first + (ir.second - ir.first) / 2;
0601             if(max_iter >= policies::get_max_root_iterations<Policy>())
0602             {
0603                return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:"
0604                   " or there is no answer to problem.  Current best guess is %1%", result, Policy());
0605             }
0606             return result;
0607          }
0608
0609          template <class RealType, class Policy>
0610          struct non_centrality_finder
0611          {
0612             non_centrality_finder(
0613                RealType v_, RealType x_, RealType p_, bool c)
0614                : v(v_), x(x_), p(p_), comp(c) {}
0615
0616             RealType operator()(const RealType& lam)
0617             {
0618                non_central_chi_squared_distribution<RealType, Policy> d(v, lam);
0619                return comp ?
0620                   RealType(p - cdf(complement(d, x)))
0621                   : RealType(cdf(d, x) - p);
0622             }
0623          private:
0624             RealType v;
0625             RealType x;
0626             RealType p;
0627             bool comp;
0628          };
0629
0630          template <class RealType, class Policy>
0631          inline RealType find_non_centrality(
0632             RealType v, RealType x, RealType p, RealType q, const Policy& pol)
0633          {
0634             const char* function = "non_central_chi_squared<%1%>::find_non_centrality";
0635             if((p == 0) || (q == 0))
0636             {
0637                //
0638                // Can't do a thing if one of p and q is zero:
0639                //
0640                return policies::raise_evaluation_error<RealType>(function,
0641                   "Can't find non centrality parameter when the probability is 0 or 1, only possible answer is %1%",
0642                   RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy());
0643             }
0644             non_centrality_finder<RealType, Policy> f(v, x, p < q ? p : q, p < q ? false : true);
0645             tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
0646             std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
0647             //
0648             // Pick an initial guess that we know will give us a probability
0649             // right around 0.5.
0650             //
0651             RealType guess = x - v;
0652             if(guess < 1)
0653                guess = 1;
0654             std::pair<RealType, RealType> ir = tools::bracket_and_solve_root(
0655                f, guess, RealType(2), false, tol, max_iter, pol);
0656             RealType result = ir.first + (ir.second - ir.first) / 2;
0657             if(max_iter >= policies::get_max_root_iterations<Policy>())
0658             {
0659                return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:"
0660                   " or there is no answer to problem.  Current best guess is %1%", result, Policy());
0661             }
0662             return result;
0663          }
0664
0665       }
0666
0667       template <class RealType = double, class Policy = policies::policy<> >
0668       class non_central_chi_squared_distribution
0669       {
0670       public:
0671          typedef RealType value_type;
0672          typedef Policy policy_type;
0673
0674          non_central_chi_squared_distribution(RealType df_, RealType lambda) : df(df_), ncp(lambda)
0675          {
0676             const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::non_central_chi_squared_distribution(%1%,%1%)";
0677             RealType r;
0678             detail::check_df(
0679                function,
0680                df, &r, Policy());
0681             detail::check_non_centrality(
0682                function,
0683                ncp,
0684                &r,
0685                Policy());
0686          } // non_central_chi_squared_distribution constructor.
0687
0688          RealType degrees_of_freedom() const
0689          { // Private data getter function.
0690             return df;
0691          }
0692          RealType non_centrality() const
0693          { // Private data getter function.
0694             return ncp;
0695          }
0696          static RealType find_degrees_of_freedom(RealType lam, RealType x, RealType p)
0697          {
0698             const char* function = "non_central_chi_squared<%1%>::find_degrees_of_freedom";
0699             typedef typename policies::evaluation<RealType, Policy>::type eval_type;
0700             typedef typename policies::normalise<
0701                Policy,
0702                policies::promote_float<false>,
0703                policies::promote_double<false>,
0704                policies::discrete_quantile<>,
0705                policies::assert_undefined<> >::type forwarding_policy;
0706             eval_type result = detail::find_degrees_of_freedom(
0707                static_cast<eval_type>(lam),
0708                static_cast<eval_type>(x),
0709                static_cast<eval_type>(p),
0710                static_cast<eval_type>(1-p),
0711                forwarding_policy());
0712             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
0713                result,
0714                function);
0715          }
0716          template <class A, class B, class C>
0717          static RealType find_degrees_of_freedom(const complemented3_type<A,B,C>& c)
0718          {
0719             const char* function = "non_central_chi_squared<%1%>::find_degrees_of_freedom";
0720             typedef typename policies::evaluation<RealType, Policy>::type eval_type;
0721             typedef typename policies::normalise<
0722                Policy,
0723                policies::promote_float<false>,
0724                policies::promote_double<false>,
0725                policies::discrete_quantile<>,
0726                policies::assert_undefined<> >::type forwarding_policy;
0727             eval_type result = detail::find_degrees_of_freedom(
0728                static_cast<eval_type>(c.dist),
0729                static_cast<eval_type>(c.param1),
0730                static_cast<eval_type>(1-c.param2),
0731                static_cast<eval_type>(c.param2),
0732                forwarding_policy());
0733             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
0734                result,
0735                function);
0736          }
0737          static RealType find_non_centrality(RealType v, RealType x, RealType p)
0738          {
0739             const char* function = "non_central_chi_squared<%1%>::find_non_centrality";
0740             typedef typename policies::evaluation<RealType, Policy>::type eval_type;
0741             typedef typename policies::normalise<
0742                Policy,
0743                policies::promote_float<false>,
0744                policies::promote_double<false>,
0745                policies::discrete_quantile<>,
0746                policies::assert_undefined<> >::type forwarding_policy;
0747             eval_type result = detail::find_non_centrality(
0748                static_cast<eval_type>(v),
0749                static_cast<eval_type>(x),
0750                static_cast<eval_type>(p),
0751                static_cast<eval_type>(1-p),
0752                forwarding_policy());
0753             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
0754                result,
0755                function);
0756          }
0757          template <class A, class B, class C>
0758          static RealType find_non_centrality(const complemented3_type<A,B,C>& c)
0759          {
0760             const char* function = "non_central_chi_squared<%1%>::find_non_centrality";
0761             typedef typename policies::evaluation<RealType, Policy>::type eval_type;
0762             typedef typename policies::normalise<
0763                Policy,
0764                policies::promote_float<false>,
0765                policies::promote_double<false>,
0766                policies::discrete_quantile<>,
0767                policies::assert_undefined<> >::type forwarding_policy;
0768             eval_type result = detail::find_non_centrality(
0769                static_cast<eval_type>(c.dist),
0770                static_cast<eval_type>(c.param1),
0771                static_cast<eval_type>(1-c.param2),
0772                static_cast<eval_type>(c.param2),
0773                forwarding_policy());
0774             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
0775                result,
0776                function);
0777          }
0778       private:
0779          // Data member, initialized by constructor.
0780          RealType df; // degrees of freedom.
0781          RealType ncp; // non-centrality parameter
0782       }; // template <class RealType, class Policy> class non_central_chi_squared_distribution
0783
0784       typedef non_central_chi_squared_distribution<double> non_central_chi_squared; // Reserved name of type double.
0785
0786       #ifdef __cpp_deduction_guides
0787       template <class RealType>
0788       non_central_chi_squared_distribution(RealType,RealType)->non_central_chi_squared_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0789       #endif
0790
0791       // Non-member functions to give properties of the distribution.
0792
0793       template <class RealType, class Policy>
0794       inline const std::pair<RealType, RealType> range(const non_central_chi_squared_distribution<RealType, Policy>& /* dist */)
0795       { // Range of permissible values for random variable k.
0796          using boost::math::tools::max_value;
0797          return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // Max integer?
0798       }
0799
0800       template <class RealType, class Policy>
0801       inline const std::pair<RealType, RealType> support(const non_central_chi_squared_distribution<RealType, Policy>& /* dist */)
0802       { // Range of supported values for random variable k.
0803          // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
0804          using boost::math::tools::max_value;
0805          return std::pair<RealType, RealType>(static_cast<RealType>(0),  max_value<RealType>());
0806       }
0807
0808       template <class RealType, class Policy>
0809       inline RealType mean(const non_central_chi_squared_distribution<RealType, Policy>& dist)
0810       { // Mean of poisson distribution = lambda.
0811          const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::mean()";
0812          RealType k = dist.degrees_of_freedom();
0813          RealType l = dist.non_centrality();
0814          RealType r;
0815          if(!detail::check_df(
0816             function,
0817             k, &r, Policy())
0818             ||
0819          !detail::check_non_centrality(
0820             function,
0821             l,
0822             &r,
0823             Policy()))
0824                return static_cast<RealType>(r);
0825          return k + l;
0826       } // mean
0827
0828       template <class RealType, class Policy>
0829       inline RealType mode(const non_central_chi_squared_distribution<RealType, Policy>& dist)
0830       { // mode.
0831          static const char* function = "mode(non_central_chi_squared_distribution<%1%> const&)";
0832
0833          RealType k = dist.degrees_of_freedom();
0834          RealType l = dist.non_centrality();
0835          RealType r;
0836          if(!detail::check_df(
0837             function,
0838             k, &r, Policy())
0839             ||
0840          !detail::check_non_centrality(
0841             function,
0842             l,
0843             &r,
0844             Policy()))
0845                return static_cast<RealType>(r);
0846          bool asymptotic_mode = k < l/4;
0847          RealType starting_point = asymptotic_mode ? k + l - RealType(3) : RealType(1) + k;
0848          return detail::generic_find_mode(dist, starting_point, function);
0849       }
0850
0851       template <class RealType, class Policy>
0852       inline RealType variance(const non_central_chi_squared_distribution<RealType, Policy>& dist)
0853       { // variance.
0854          const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::variance()";
0855          RealType k = dist.degrees_of_freedom();
0856          RealType l = dist.non_centrality();
0857          RealType r;
0858          if(!detail::check_df(
0859             function,
0860             k, &r, Policy())
0861             ||
0862          !detail::check_non_centrality(
0863             function,
0864             l,
0865             &r,
0866             Policy()))
0867                return static_cast<RealType>(r);
0868          return 2 * (2 * l + k);
0869       }
0870
0871       // RealType standard_deviation(const non_central_chi_squared_distribution<RealType, Policy>& dist)
0872       // standard_deviation provided by derived accessors.
0873
0874       template <class RealType, class Policy>
0875       inline RealType skewness(const non_central_chi_squared_distribution<RealType, Policy>& dist)
0876       { // skewness = sqrt(l).
0877          const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::skewness()";
0878          RealType k = dist.degrees_of_freedom();
0879          RealType l = dist.non_centrality();
0880          RealType r;
0881          if(!detail::check_df(
0882             function,
0883             k, &r, Policy())
0884             ||
0885          !detail::check_non_centrality(
0886             function,
0887             l,
0888             &r,
0889             Policy()))
0890                return static_cast<RealType>(r);
0891          BOOST_MATH_STD_USING
0892             return pow(2 / (k + 2 * l), RealType(3)/2) * (k + 3 * l);
0893       }
0894
0895       template <class RealType, class Policy>
0896       inline RealType kurtosis_excess(const non_central_chi_squared_distribution<RealType, Policy>& dist)
0897       {
0898          const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::kurtosis_excess()";
0899          RealType k = dist.degrees_of_freedom();
0900          RealType l = dist.non_centrality();
0901          RealType r;
0902          if(!detail::check_df(
0903             function,
0904             k, &r, Policy())
0905             ||
0906          !detail::check_non_centrality(
0907             function,
0908             l,
0909             &r,
0910             Policy()))
0911                return static_cast<RealType>(r);
0912          return 12 * (k + 4 * l) / ((k + 2 * l) * (k + 2 * l));
0913       } // kurtosis_excess
0914
0915       template <class RealType, class Policy>
0916       inline RealType kurtosis(const non_central_chi_squared_distribution<RealType, Policy>& dist)
0917       {
0918          return kurtosis_excess(dist) + 3;
0919       }
0920
0921       template <class RealType, class Policy>
0922       inline RealType pdf(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
0923       { // Probability Density/Mass Function.
0924          return detail::nccs_pdf(dist, x);
0925       } // pdf
0926
0927       template <class RealType, class Policy>
0928       RealType cdf(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
0929       {
0930          const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::cdf(%1%)";
0931          RealType k = dist.degrees_of_freedom();
0932          RealType l = dist.non_centrality();
0933          RealType r;
0934          if(!detail::check_df(
0935             function,
0936             k, &r, Policy())
0937             ||
0938          !detail::check_non_centrality(
0939             function,
0940             l,
0941             &r,
0942             Policy())
0943             ||
0944          !detail::check_positive_x(
0945             function,
0946             x,
0947             &r,
0948             Policy()))
0949                return static_cast<RealType>(r);
0950
0951          return detail::non_central_chi_squared_cdf(x, k, l, false, Policy());
0952       } // cdf
0953
0954       template <class RealType, class Policy>
0955       RealType cdf(const complemented2_type<non_central_chi_squared_distribution<RealType, Policy>, RealType>& c)
0956       { // Complemented Cumulative Distribution Function
0957          const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::cdf(%1%)";
0958          non_central_chi_squared_distribution<RealType, Policy> const& dist = c.dist;
0959          RealType x = c.param;
0960          RealType k = dist.degrees_of_freedom();
0961          RealType l = dist.non_centrality();
0962          RealType r;
0963          if(!detail::check_df(
0964             function,
0965             k, &r, Policy())
0966             ||
0967          !detail::check_non_centrality(
0968             function,
0969             l,
0970             &r,
0971             Policy())
0972             ||
0973          !detail::check_positive_x(
0974             function,
0975             x,
0976             &r,
0977             Policy()))
0978                return static_cast<RealType>(r);
0979
0980          return detail::non_central_chi_squared_cdf(x, k, l, true, Policy());
0981       } // ccdf
0982
0983       template <class RealType, class Policy>
0984       inline RealType quantile(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& p)
0985       { // Quantile (or Percent Point) function.
0986          return detail::nccs_quantile(dist, p, false);
0987       } // quantile
0988
0989       template <class RealType, class Policy>
0990       inline RealType quantile(const complemented2_type<non_central_chi_squared_distribution<RealType, Policy>, RealType>& c)
0991       { // Quantile (or Percent Point) function.
0992          return detail::nccs_quantile(c.dist, c.param, true);
0993       } // quantile complement.
0994
0995    } // namespace math
0996 } // namespace boost
0997
0998 // This include must be at the end, *after* the accessors
0999 // for this distribution have been defined, in order to
1000 // keep compilers that support two-phase lookup happy.
1001 #include <boost/math/distributions/detail/derived_accessors.hpp>
1002
1003 #endif // BOOST_MATH_SPECIAL_NON_CENTRAL_CHI_SQUARE_HPP