Chris@16: // boost\math\distributions\non_central_chi_squared.hpp
Chris@16: 
Chris@16: // Copyright John Maddock 2008.
Chris@16: 
Chris@16: // Use, modification and distribution are subject to the
Chris@16: // Boost Software License, Version 1.0.
Chris@16: // (See accompanying file LICENSE_1_0.txt
Chris@16: // or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_MATH_SPECIAL_NON_CENTRAL_CHI_SQUARE_HPP
Chris@16: #define BOOST_MATH_SPECIAL_NON_CENTRAL_CHI_SQUARE_HPP
Chris@16: 
Chris@16: #include <boost/math/distributions/fwd.hpp>
Chris@16: #include <boost/math/special_functions/gamma.hpp> // for incomplete gamma. gamma_q
Chris@16: #include <boost/math/special_functions/bessel.hpp> // for cyl_bessel_i
Chris@16: #include <boost/math/special_functions/round.hpp> // for iround
Chris@16: #include <boost/math/distributions/complement.hpp> // complements
Chris@16: #include <boost/math/distributions/chi_squared.hpp> // central distribution
Chris@16: #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks
Chris@16: #include <boost/math/special_functions/fpclassify.hpp> // isnan.
Chris@16: #include <boost/math/tools/roots.hpp> // for root finding.
Chris@16: #include <boost/math/distributions/detail/generic_mode.hpp>
Chris@16: #include <boost/math/distributions/detail/generic_quantile.hpp>
Chris@16: 
Chris@16: namespace boost
Chris@16: {
Chris@16:    namespace math
Chris@16:    {
Chris@16: 
Chris@16:       template <class RealType, class Policy>
Chris@16:       class non_central_chi_squared_distribution;
Chris@16: 
Chris@16:       namespace detail{
Chris@16: 
Chris@16:          template <class T, class Policy>
Chris@16:          T non_central_chi_square_q(T x, T f, T theta, const Policy& pol, T init_sum = 0)
Chris@16:          {
Chris@16:             //
Chris@16:             // Computes the complement of the Non-Central Chi-Square
Chris@16:             // Distribution CDF by summing a weighted sum of complements
Chris@16:             // of the central-distributions.  The weighting factor is
Chris@16:             // a Poisson Distribution.
Chris@16:             //
Chris@16:             // This is an application of the technique described in:
Chris@16:             //
Chris@16:             // Computing discrete mixtures of continuous
Chris@16:             // distributions: noncentral chisquare, noncentral t
Chris@16:             // and the distribution of the square of the sample
Chris@16:             // multiple correlation coeficient.
Chris@16:             // D. Benton, K. Krishnamoorthy.
Chris@16:             // Computational Statistics & Data Analysis 43 (2003) 249 - 267
Chris@16:             //
Chris@16:             BOOST_MATH_STD_USING
Chris@16: 
Chris@16:             // Special case:
Chris@16:             if(x == 0)
Chris@16:                return 1;
Chris@16: 
Chris@16:             //
Chris@16:             // Initialize the variables we'll be using:
Chris@16:             //
Chris@16:             T lambda = theta / 2;
Chris@16:             T del = f / 2;
Chris@16:             T y = x / 2;
Chris@16:             boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
Chris@16:             T errtol = boost::math::policies::get_epsilon<T, Policy>();
Chris@16:             T sum = init_sum;
Chris@16:             //
Chris@16:             // k is the starting location for iteration, we'll
Chris@16:             // move both forwards and backwards from this point.
Chris@16:             // k is chosen as the peek of the Poisson weights, which
Chris@16:             // will occur *before* the largest term.
Chris@16:             //
Chris@16:             int k = iround(lambda, pol);
Chris@16:             // Forwards and backwards Poisson weights:
Chris@16:             T poisf = boost::math::gamma_p_derivative(1 + k, lambda, pol);
Chris@16:             T poisb = poisf * k / lambda;
Chris@16:             // Initial forwards central chi squared term:
Chris@16:             T gamf = boost::math::gamma_q(del + k, y, pol);
Chris@16:             // Forwards and backwards recursion terms on the central chi squared:
Chris@16:             T xtermf = boost::math::gamma_p_derivative(del + 1 + k, y, pol);
Chris@16:             T xtermb = xtermf * (del + k) / y;
Chris@16:             // Initial backwards central chi squared term:
Chris@16:             T gamb = gamf - xtermb;
Chris@16: 
Chris@16:             //
Chris@16:             // Forwards iteration first, this is the
Chris@16:             // stable direction for the gamma function
Chris@16:             // recurrences:
Chris@16:             //
Chris@16:             int i;
Chris@16:             for(i = k; static_cast<boost::uintmax_t>(i-k) < max_iter; ++i)
Chris@16:             {
Chris@16:                T term = poisf * gamf;
Chris@16:                sum += term;
Chris@16:                poisf *= lambda / (i + 1);
Chris@16:                gamf += xtermf;
Chris@16:                xtermf *= y / (del + i + 1);
Chris@16:                if(((sum == 0) || (fabs(term / sum) < errtol)) && (term >= poisf * gamf))
Chris@16:                   break;
Chris@16:             }
Chris@16:             //Error check:
Chris@16:             if(static_cast<boost::uintmax_t>(i-k) >= max_iter)
Chris@101:                return policies::raise_evaluation_error(
Chris@16:                   "cdf(non_central_chi_squared_distribution<%1%>, %1%)",
Chris@16:                   "Series did not converge, closest value was %1%", sum, pol);
Chris@16:             //
Chris@16:             // Now backwards iteration: the gamma
Chris@16:             // function recurrences are unstable in this
Chris@16:             // direction, we rely on the terms deminishing in size
Chris@16:             // faster than we introduce cancellation errors.
Chris@16:             // For this reason it's very important that we start
Chris@16:             // *before* the largest term so that backwards iteration
Chris@16:             // is strictly converging.
Chris@16:             //
Chris@16:             for(i = k - 1; i >= 0; --i)
Chris@16:             {
Chris@16:                T term = poisb * gamb;
Chris@16:                sum += term;
Chris@16:                poisb *= i / lambda;
Chris@16:                xtermb *= (del + i) / y;
Chris@16:                gamb -= xtermb;
Chris@16:                if((sum == 0) || (fabs(term / sum) < errtol))
Chris@16:                   break;
Chris@16:             }
Chris@16: 
Chris@16:             return sum;
Chris@16:          }
Chris@16: 
Chris@16:          template <class T, class Policy>
Chris@16:          T non_central_chi_square_p_ding(T x, T f, T theta, const Policy& pol, T init_sum = 0)
Chris@16:          {
Chris@16:             //
Chris@16:             // This is an implementation of:
Chris@16:             //
Chris@16:             // Algorithm AS 275:
Chris@16:             // Computing the Non-Central #2 Distribution Function
Chris@16:             // Cherng G. Ding
Chris@16:             // Applied Statistics, Vol. 41, No. 2. (1992), pp. 478-482.
Chris@16:             //
Chris@16:             // This uses a stable forward iteration to sum the
Chris@16:             // CDF, unfortunately this can not be used for large
Chris@16:             // values of the non-centrality parameter because:
Chris@16:             // * The first term may underfow to zero.
Chris@16:             // * We may need an extra-ordinary number of terms
Chris@16:             //   before we reach the first *significant* term.
Chris@16:             //
Chris@16:             BOOST_MATH_STD_USING
Chris@16:             // Special case:
Chris@16:             if(x == 0)
Chris@16:                return 0;
Chris@16:             T tk = boost::math::gamma_p_derivative(f/2 + 1, x/2, pol);
Chris@16:             T lambda = theta / 2;
Chris@16:             T vk = exp(-lambda);
Chris@16:             T uk = vk;
Chris@16:             T sum = init_sum + tk * vk;
Chris@16:             if(sum == 0)
Chris@16:                return sum;
Chris@16: 
Chris@16:             boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
Chris@16:             T errtol = boost::math::policies::get_epsilon<T, Policy>();
Chris@16: 
Chris@16:             int i;
Chris@16:             T lterm(0), term(0);
Chris@16:             for(i = 1; static_cast<boost::uintmax_t>(i) < max_iter; ++i)
Chris@16:             {
Chris@16:                tk = tk * x / (f + 2 * i);
Chris@16:                uk = uk * lambda / i;
Chris@16:                vk = vk + uk;
Chris@16:                lterm = term;
Chris@16:                term = vk * tk;
Chris@16:                sum += term;
Chris@16:                if((fabs(term / sum) < errtol) && (term <= lterm))
Chris@16:                   break;
Chris@16:             }
Chris@16:             //Error check:
Chris@16:             if(static_cast<boost::uintmax_t>(i) >= max_iter)
Chris@101:                return policies::raise_evaluation_error(
Chris@16:                   "cdf(non_central_chi_squared_distribution<%1%>, %1%)",
Chris@16:                   "Series did not converge, closest value was %1%", sum, pol);
Chris@16:             return sum;
Chris@16:          }
Chris@16: 
Chris@16: 
Chris@16:          template <class T, class Policy>
Chris@16:          T non_central_chi_square_p(T y, T n, T lambda, const Policy& pol, T init_sum)
Chris@16:          {
Chris@16:             //
Chris@16:             // This is taken more or less directly from:
Chris@16:             //
Chris@16:             // Computing discrete mixtures of continuous
Chris@16:             // distributions: noncentral chisquare, noncentral t
Chris@16:             // and the distribution of the square of the sample
Chris@16:             // multiple correlation coeficient.
Chris@16:             // D. Benton, K. Krishnamoorthy.
Chris@16:             // Computational Statistics & Data Analysis 43 (2003) 249 - 267
Chris@16:             //
Chris@16:             // We're summing a Poisson weighting term multiplied by
Chris@16:             // a central chi squared distribution.
Chris@16:             //
Chris@16:             BOOST_MATH_STD_USING
Chris@16:             // Special case:
Chris@16:             if(y == 0)
Chris@16:                return 0;
Chris@16:             boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
Chris@16:             T errtol = boost::math::policies::get_epsilon<T, Policy>();
Chris@16:             T errorf(0), errorb(0);
Chris@16: 
Chris@16:             T x = y / 2;
Chris@16:             T del = lambda / 2;
Chris@16:             //
Chris@16:             // Starting location for the iteration, we'll iterate
Chris@16:             // both forwards and backwards from this point.  The
Chris@16:             // location chosen is the maximum of the Poisson weight
Chris@16:             // function, which ocurrs *after* the largest term in the
Chris@16:             // sum.
Chris@16:             //
Chris@16:             int k = iround(del, pol);
Chris@16:             T a = n / 2 + k;
Chris@16:             // Central chi squared term for forward iteration:
Chris@16:             T gamkf = boost::math::gamma_p(a, x, pol);
Chris@16: 
Chris@16:             if(lambda == 0)
Chris@16:                return gamkf;
Chris@16:             // Central chi squared term for backward iteration:
Chris@16:             T gamkb = gamkf;
Chris@16:             // Forwards Poisson weight:
Chris@16:             T poiskf = gamma_p_derivative(k+1, del, pol);
Chris@16:             // Backwards Poisson weight:
Chris@16:             T poiskb = poiskf;
Chris@16:             // Forwards gamma function recursion term:
Chris@16:             T xtermf = boost::math::gamma_p_derivative(a, x, pol);
Chris@16:             // Backwards gamma function recursion term:
Chris@16:             T xtermb = xtermf * x / a;
Chris@16:             T sum = init_sum + poiskf * gamkf;
Chris@16:             if(sum == 0)
Chris@16:                return sum;
Chris@16:             int i = 1;
Chris@16:             //
Chris@16:             // Backwards recursion first, this is the stable
Chris@16:             // direction for gamma function recurrences:
Chris@16:             //
Chris@16:             while(i <= k)
Chris@16:             {
Chris@16:                xtermb *= (a - i + 1) / x;
Chris@16:                gamkb += xtermb;
Chris@16:                poiskb = poiskb * (k - i + 1) / del;
Chris@16:                errorf = errorb;
Chris@16:                errorb = gamkb * poiskb;
Chris@16:                sum += errorb;
Chris@16:                if((fabs(errorb / sum) < errtol) && (errorb <= errorf))
Chris@16:                   break;
Chris@16:                ++i;
Chris@16:             }
Chris@16:             i = 1;
Chris@16:             //
Chris@16:             // Now forwards recursion, the gamma function
Chris@16:             // recurrence relation is unstable in this direction,
Chris@16:             // so we rely on the magnitude of successive terms
Chris@16:             // decreasing faster than we introduce cancellation error.
Chris@16:             // For this reason it's vital that k is chosen to be *after*
Chris@16:             // the largest term, so that successive forward iterations
Chris@16:             // are strictly (and rapidly) converging.
Chris@16:             //
Chris@16:             do
Chris@16:             {
Chris@16:                xtermf = xtermf * x / (a + i - 1);
Chris@16:                gamkf = gamkf - xtermf;
Chris@16:                poiskf = poiskf * del / (k + i);
Chris@16:                errorf = poiskf * gamkf;
Chris@16:                sum += errorf;
Chris@16:                ++i;
Chris@16:             }while((fabs(errorf / sum) > errtol) && (static_cast<boost::uintmax_t>(i) < max_iter));
Chris@16: 
Chris@16:             //Error check:
Chris@16:             if(static_cast<boost::uintmax_t>(i) >= max_iter)
Chris@101:                return policies::raise_evaluation_error(
Chris@16:                   "cdf(non_central_chi_squared_distribution<%1%>, %1%)",
Chris@16:                   "Series did not converge, closest value was %1%", sum, pol);
Chris@16: 
Chris@16:             return sum;
Chris@16:          }
Chris@16: 
Chris@16:          template <class T, class Policy>
Chris@16:          T non_central_chi_square_pdf(T x, T n, T lambda, const Policy& pol)
Chris@16:          {
Chris@16:             //
Chris@16:             // As above but for the PDF:
Chris@16:             //
Chris@16:             BOOST_MATH_STD_USING
Chris@16:             boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
Chris@16:             T errtol = boost::math::policies::get_epsilon<T, Policy>();
Chris@16:             T x2 = x / 2;
Chris@16:             T n2 = n / 2;
Chris@16:             T l2 = lambda / 2;
Chris@16:             T sum = 0;
Chris@16:             int k = itrunc(l2);
Chris@16:             T pois = gamma_p_derivative(k + 1, l2, pol) * gamma_p_derivative(n2 + k, x2);
Chris@16:             if(pois == 0)
Chris@16:                return 0;
Chris@16:             T poisb = pois;
Chris@16:             for(int i = k; ; ++i)
Chris@16:             {
Chris@16:                sum += pois;
Chris@16:                if(pois / sum < errtol)
Chris@16:                   break;
Chris@16:                if(static_cast<boost::uintmax_t>(i - k) >= max_iter)
Chris@16:                   return policies::raise_evaluation_error(
Chris@16:                      "pdf(non_central_chi_squared_distribution<%1%>, %1%)",
Chris@16:                      "Series did not converge, closest value was %1%", sum, pol);
Chris@16:                pois *= l2 * x2 / ((i + 1) * (n2 + i));
Chris@16:             }
Chris@16:             for(int i = k - 1; i >= 0; --i)
Chris@16:             {
Chris@16:                poisb *= (i + 1) * (n2 + i) / (l2 * x2);
Chris@16:                sum += poisb;
Chris@16:                if(poisb / sum < errtol)
Chris@16:                   break;
Chris@16:             }
Chris@16:             return sum / 2;
Chris@16:          }
Chris@16: 
Chris@16:          template <class RealType, class Policy>
Chris@16:          inline RealType non_central_chi_squared_cdf(RealType x, RealType k, RealType l, bool invert, const Policy&)
Chris@16:          {
Chris@16:             typedef typename policies::evaluation<RealType, Policy>::type value_type;
Chris@16:             typedef typename policies::normalise<
Chris@16:                Policy,
Chris@16:                policies::promote_float<false>,
Chris@16:                policies::promote_double<false>,
Chris@16:                policies::discrete_quantile<>,
Chris@16:                policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:             BOOST_MATH_STD_USING
Chris@16:             value_type result;
Chris@16:             if(l == 0)
Chris@16:                result = cdf(boost::math::chi_squared_distribution<RealType, Policy>(k), x);
Chris@16:             else if(x > k + l)
Chris@16:             {
Chris@16:                // Complement is the smaller of the two:
Chris@16:                result = detail::non_central_chi_square_q(
Chris@16:                   static_cast<value_type>(x),
Chris@16:                   static_cast<value_type>(k),
Chris@16:                   static_cast<value_type>(l),
Chris@16:                   forwarding_policy(),
Chris@16:                   static_cast<value_type>(invert ? 0 : -1));
Chris@16:                invert = !invert;
Chris@16:             }
Chris@16:             else if(l < 200)
Chris@16:             {
Chris@16:                // For small values of the non-centrality parameter
Chris@16:                // we can use Ding's method:
Chris@16:                result = detail::non_central_chi_square_p_ding(
Chris@16:                   static_cast<value_type>(x),
Chris@16:                   static_cast<value_type>(k),
Chris@16:                   static_cast<value_type>(l),
Chris@16:                   forwarding_policy(),
Chris@16:                   static_cast<value_type>(invert ? -1 : 0));
Chris@16:             }
Chris@16:             else
Chris@16:             {
Chris@16:                // For largers values of the non-centrality
Chris@16:                // parameter Ding's method will consume an
Chris@16:                // extra-ordinary number of terms, and worse
Chris@16:                // may return zero when the result is in fact
Chris@16:                // finite, use Krishnamoorthy's method instead:
Chris@16:                result = detail::non_central_chi_square_p(
Chris@16:                   static_cast<value_type>(x),
Chris@16:                   static_cast<value_type>(k),
Chris@16:                   static_cast<value_type>(l),
Chris@16:                   forwarding_policy(),
Chris@16:                   static_cast<value_type>(invert ? -1 : 0));
Chris@16:             }
Chris@16:             if(invert)
Chris@16:                result = -result;
Chris@16:             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
Chris@16:                result,
Chris@16:                "boost::math::non_central_chi_squared_cdf<%1%>(%1%, %1%, %1%)");
Chris@16:          }
Chris@16: 
Chris@16:          template <class T, class Policy>
Chris@16:          struct nccs_quantile_functor
Chris@16:          {
Chris@16:             nccs_quantile_functor(const non_central_chi_squared_distribution<T,Policy>& d, T t, bool c)
Chris@16:                : dist(d), target(t), comp(c) {}
Chris@16: 
Chris@16:             T operator()(const T& x)
Chris@16:             {
Chris@16:                return comp ?
Chris@16:                   target - cdf(complement(dist, x))
Chris@16:                   : cdf(dist, x) - target;
Chris@16:             }
Chris@16: 
Chris@16:          private:
Chris@16:             non_central_chi_squared_distribution<T,Policy> dist;
Chris@16:             T target;
Chris@16:             bool comp;
Chris@16:          };
Chris@16: 
Chris@16:          template <class RealType, class Policy>
Chris@16:          RealType nccs_quantile(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& p, bool comp)
Chris@16:          {
Chris@16:             BOOST_MATH_STD_USING
Chris@16:             static const char* function = "quantile(non_central_chi_squared_distribution<%1%>, %1%)";
Chris@16:             typedef typename policies::evaluation<RealType, Policy>::type value_type;
Chris@16:             typedef typename policies::normalise<
Chris@16:                Policy,
Chris@16:                policies::promote_float<false>,
Chris@16:                policies::promote_double<false>,
Chris@16:                policies::discrete_quantile<>,
Chris@16:                policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:             value_type k = dist.degrees_of_freedom();
Chris@16:             value_type l = dist.non_centrality();
Chris@16:             value_type r;
Chris@16:             if(!detail::check_df(
Chris@16:                function,
Chris@16:                k, &r, Policy())
Chris@16:                ||
Chris@16:             !detail::check_non_centrality(
Chris@16:                function,
Chris@16:                l,
Chris@16:                &r,
Chris@16:                Policy())
Chris@16:                ||
Chris@16:             !detail::check_probability(
Chris@16:                function,
Chris@16:                static_cast<value_type>(p),
Chris@16:                &r,
Chris@16:                Policy()))
Chris@16:                   return (RealType)r;
Chris@16:             //
Chris@16:             // Special cases get short-circuited first:
Chris@16:             //
Chris@16:             if(p == 0)
Chris@101:                return comp ? policies::raise_overflow_error<RealType>(function, 0, Policy()) : 0;
Chris@16:             if(p == 1)
Chris@101:                return comp ? 0 : policies::raise_overflow_error<RealType>(function, 0, Policy());
Chris@16:             //
Chris@16:             // This is Pearson's approximation to the quantile, see
Chris@16:             // Pearson, E. S. (1959) "Note on an approximation to the distribution of 
Chris@16:             // noncentral chi squared", Biometrika 46: 364.
Chris@16:             // See also:
Chris@16:             // "A comparison of approximations to percentiles of the noncentral chi2-distribution",
Chris@16:             // Hardeo Sahai and Mario Miguel Ojeda, Revista de Matematica: Teoria y Aplicaciones 2003 10(1–2) : 57–76.
Chris@16:             // Note that the latter reference refers to an approximation of the CDF, when they really mean the quantile.
Chris@16:             //
Chris@16:             value_type b = -(l * l) / (k + 3 * l);
Chris@16:             value_type c = (k + 3 * l) / (k + 2 * l);
Chris@16:             value_type ff = (k + 2 * l) / (c * c);
Chris@16:             value_type guess;
Chris@16:             if(comp)
Chris@16:             {
Chris@16:                guess = b + c * quantile(complement(chi_squared_distribution<value_type, forwarding_policy>(ff), p));
Chris@16:             }
Chris@16:             else
Chris@16:             {
Chris@16:                guess = b + c * quantile(chi_squared_distribution<value_type, forwarding_policy>(ff), p);
Chris@16:             }
Chris@16:             //
Chris@16:             // Sometimes guess goes very small or negative, in that case we have
Chris@16:             // to do something else for the initial guess, this approximation
Chris@16:             // was provided in a private communication from Thomas Luu, PhD candidate, 
Chris@16:             // University College London.  It's an asymptotic expansion for the
Chris@16:             // quantile which usually gets us within an order of magnitude of the
Chris@16:             // correct answer.
Chris@16:             //
Chris@16:             if(guess < 0.005)
Chris@16:             {
Chris@16:                value_type pp = comp ? 1 - p : p;
Chris@16:                //guess = pow(pow(value_type(2), (k / 2 - 1)) * exp(l / 2) * pp * k, 2 / k);
Chris@16:                guess = pow(pow(value_type(2), (k / 2 - 1)) * exp(l / 2) * pp * k * boost::math::tgamma(k / 2, forwarding_policy()), (2 / k));
Chris@16:                if(guess == 0)
Chris@16:                   guess = tools::min_value<value_type>();
Chris@16:             }
Chris@16:             value_type result = detail::generic_quantile(
Chris@16:                non_central_chi_squared_distribution<value_type, forwarding_policy>(k, l),
Chris@16:                p,
Chris@16:                guess,
Chris@16:                comp,
Chris@16:                function);
Chris@16: 
Chris@16:             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
Chris@16:                result,
Chris@16:                function);
Chris@16:          }
Chris@16: 
Chris@16:          template <class RealType, class Policy>
Chris@16:          RealType nccs_pdf(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
Chris@16:          {
Chris@16:             BOOST_MATH_STD_USING
Chris@16:             static const char* function = "pdf(non_central_chi_squared_distribution<%1%>, %1%)";
Chris@16:             typedef typename policies::evaluation<RealType, Policy>::type value_type;
Chris@16:             typedef typename policies::normalise<
Chris@16:                Policy,
Chris@16:                policies::promote_float<false>,
Chris@16:                policies::promote_double<false>,
Chris@16:                policies::discrete_quantile<>,
Chris@16:                policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:             value_type k = dist.degrees_of_freedom();
Chris@16:             value_type l = dist.non_centrality();
Chris@16:             value_type r;
Chris@16:             if(!detail::check_df(
Chris@16:                function,
Chris@16:                k, &r, Policy())
Chris@16:                ||
Chris@16:             !detail::check_non_centrality(
Chris@16:                function,
Chris@16:                l,
Chris@16:                &r,
Chris@16:                Policy())
Chris@16:                ||
Chris@16:             !detail::check_positive_x(
Chris@16:                function,
Chris@16:                (value_type)x,
Chris@16:                &r,
Chris@16:                Policy()))
Chris@16:                   return (RealType)r;
Chris@16: 
Chris@16:          if(l == 0)
Chris@16:             return pdf(boost::math::chi_squared_distribution<RealType, forwarding_policy>(dist.degrees_of_freedom()), x);
Chris@16: 
Chris@16:          // Special case:
Chris@16:          if(x == 0)
Chris@16:             return 0;
Chris@16:          if(l > 50)
Chris@16:          {
Chris@16:             r = non_central_chi_square_pdf(static_cast<value_type>(x), k, l, forwarding_policy());
Chris@16:          }
Chris@16:          else
Chris@16:          {
Chris@16:             r = log(x / l) * (k / 4 - 0.5f) - (x + l) / 2;
Chris@16:             if(fabs(r) >= tools::log_max_value<RealType>() / 4)
Chris@16:             {
Chris@16:                r = non_central_chi_square_pdf(static_cast<value_type>(x), k, l, forwarding_policy());
Chris@16:             }
Chris@16:             else
Chris@16:             {
Chris@16:                r = exp(r);
Chris@16:                r = 0.5f * r
Chris@16:                   * boost::math::cyl_bessel_i(k/2 - 1, sqrt(l * x), forwarding_policy());
Chris@16:             }
Chris@16:          }
Chris@16:          return policies::checked_narrowing_cast<RealType, forwarding_policy>(
Chris@16:                r,
Chris@16:                function);
Chris@16:          }
Chris@16: 
Chris@16:          template <class RealType, class Policy>
Chris@16:          struct degrees_of_freedom_finder
Chris@16:          {
Chris@16:             degrees_of_freedom_finder(
Chris@16:                RealType lam_, RealType x_, RealType p_, bool c)
Chris@16:                : lam(lam_), x(x_), p(p_), comp(c) {}
Chris@16: 
Chris@16:             RealType operator()(const RealType& v)
Chris@16:             {
Chris@16:                non_central_chi_squared_distribution<RealType, Policy> d(v, lam);
Chris@16:                return comp ?
Chris@16:                   RealType(p - cdf(complement(d, x)))
Chris@16:                   : RealType(cdf(d, x) - p);
Chris@16:             }
Chris@16:          private:
Chris@16:             RealType lam;
Chris@16:             RealType x;
Chris@16:             RealType p;
Chris@16:             bool comp;
Chris@16:          };
Chris@16: 
Chris@16:          template <class RealType, class Policy>
Chris@16:          inline RealType find_degrees_of_freedom(
Chris@16:             RealType lam, RealType x, RealType p, RealType q, const Policy& pol)
Chris@16:          {
Chris@16:             const char* function = "non_central_chi_squared<%1%>::find_degrees_of_freedom";
Chris@16:             if((p == 0) || (q == 0))
Chris@16:             {
Chris@16:                //
Chris@16:                // Can't a thing if one of p and q is zero:
Chris@16:                //
Chris@16:                return policies::raise_evaluation_error<RealType>(function,
Chris@16:                   "Can't find degrees of freedom when the probability is 0 or 1, only possible answer is %1%",
Chris@16:                   RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy());
Chris@16:             }
Chris@16:             degrees_of_freedom_finder<RealType, Policy> f(lam, x, p < q ? p : q, p < q ? false : true);
Chris@16:             tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
Chris@16:             boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
Chris@16:             //
Chris@16:             // Pick an initial guess that we know will give us a probability
Chris@16:             // right around 0.5.
Chris@16:             //
Chris@16:             RealType guess = x - lam;
Chris@16:             if(guess < 1)
Chris@16:                guess = 1;
Chris@16:             std::pair<RealType, RealType> ir = tools::bracket_and_solve_root(
Chris@16:                f, guess, RealType(2), false, tol, max_iter, pol);
Chris@16:             RealType result = ir.first + (ir.second - ir.first) / 2;
Chris@16:             if(max_iter >= policies::get_max_root_iterations<Policy>())
Chris@16:             {
Chris@101:                return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:"
Chris@16:                   " or there is no answer to problem.  Current best guess is %1%", result, Policy());
Chris@16:             }
Chris@16:             return result;
Chris@16:          }
Chris@16: 
Chris@16:          template <class RealType, class Policy>
Chris@16:          struct non_centrality_finder
Chris@16:          {
Chris@16:             non_centrality_finder(
Chris@16:                RealType v_, RealType x_, RealType p_, bool c)
Chris@16:                : v(v_), x(x_), p(p_), comp(c) {}
Chris@16: 
Chris@16:             RealType operator()(const RealType& lam)
Chris@16:             {
Chris@16:                non_central_chi_squared_distribution<RealType, Policy> d(v, lam);
Chris@16:                return comp ?
Chris@16:                   RealType(p - cdf(complement(d, x)))
Chris@16:                   : RealType(cdf(d, x) - p);
Chris@16:             }
Chris@16:          private:
Chris@16:             RealType v;
Chris@16:             RealType x;
Chris@16:             RealType p;
Chris@16:             bool comp;
Chris@16:          };
Chris@16: 
Chris@16:          template <class RealType, class Policy>
Chris@16:          inline RealType find_non_centrality(
Chris@16:             RealType v, RealType x, RealType p, RealType q, const Policy& pol)
Chris@16:          {
Chris@16:             const char* function = "non_central_chi_squared<%1%>::find_non_centrality";
Chris@16:             if((p == 0) || (q == 0))
Chris@16:             {
Chris@16:                //
Chris@16:                // Can't do a thing if one of p and q is zero:
Chris@16:                //
Chris@16:                return policies::raise_evaluation_error<RealType>(function,
Chris@16:                   "Can't find non centrality parameter when the probability is 0 or 1, only possible answer is %1%",
Chris@16:                   RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy());
Chris@16:             }
Chris@16:             non_centrality_finder<RealType, Policy> f(v, x, p < q ? p : q, p < q ? false : true);
Chris@16:             tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
Chris@16:             boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
Chris@16:             //
Chris@16:             // Pick an initial guess that we know will give us a probability
Chris@16:             // right around 0.5.
Chris@16:             //
Chris@16:             RealType guess = x - v;
Chris@16:             if(guess < 1)
Chris@16:                guess = 1;
Chris@16:             std::pair<RealType, RealType> ir = tools::bracket_and_solve_root(
Chris@16:                f, guess, RealType(2), false, tol, max_iter, pol);
Chris@16:             RealType result = ir.first + (ir.second - ir.first) / 2;
Chris@16:             if(max_iter >= policies::get_max_root_iterations<Policy>())
Chris@16:             {
Chris@101:                return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:"
Chris@16:                   " or there is no answer to problem.  Current best guess is %1%", result, Policy());
Chris@16:             }
Chris@16:             return result;
Chris@16:          }
Chris@16: 
Chris@16:       }
Chris@16: 
Chris@16:       template <class RealType = double, class Policy = policies::policy<> >
Chris@16:       class non_central_chi_squared_distribution
Chris@16:       {
Chris@16:       public:
Chris@16:          typedef RealType value_type;
Chris@16:          typedef Policy policy_type;
Chris@16: 
Chris@16:          non_central_chi_squared_distribution(RealType df_, RealType lambda) : df(df_), ncp(lambda)
Chris@16:          {
Chris@16:             const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::non_central_chi_squared_distribution(%1%,%1%)";
Chris@16:             RealType r;
Chris@16:             detail::check_df(
Chris@16:                function,
Chris@16:                df, &r, Policy());
Chris@16:             detail::check_non_centrality(
Chris@16:                function,
Chris@16:                ncp,
Chris@16:                &r,
Chris@16:                Policy());
Chris@16:          } // non_central_chi_squared_distribution constructor.
Chris@16: 
Chris@16:          RealType degrees_of_freedom() const
Chris@16:          { // Private data getter function.
Chris@16:             return df;
Chris@16:          }
Chris@16:          RealType non_centrality() const
Chris@16:          { // Private data getter function.
Chris@16:             return ncp;
Chris@16:          }
Chris@16:          static RealType find_degrees_of_freedom(RealType lam, RealType x, RealType p)
Chris@16:          {
Chris@16:             const char* function = "non_central_chi_squared<%1%>::find_degrees_of_freedom";
Chris@16:             typedef typename policies::evaluation<RealType, Policy>::type value_type;
Chris@16:             typedef typename policies::normalise<
Chris@16:                Policy,
Chris@16:                policies::promote_float<false>,
Chris@16:                policies::promote_double<false>,
Chris@16:                policies::discrete_quantile<>,
Chris@16:                policies::assert_undefined<> >::type forwarding_policy;
Chris@16:             value_type result = detail::find_degrees_of_freedom(
Chris@16:                static_cast<value_type>(lam),
Chris@16:                static_cast<value_type>(x),
Chris@16:                static_cast<value_type>(p),
Chris@16:                static_cast<value_type>(1-p),
Chris@16:                forwarding_policy());
Chris@16:             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
Chris@16:                result,
Chris@16:                function);
Chris@16:          }
Chris@16:          template <class A, class B, class C>
Chris@16:          static RealType find_degrees_of_freedom(const complemented3_type<A,B,C>& c)
Chris@16:          {
Chris@16:             const char* function = "non_central_chi_squared<%1%>::find_degrees_of_freedom";
Chris@16:             typedef typename policies::evaluation<RealType, Policy>::type value_type;
Chris@16:             typedef typename policies::normalise<
Chris@16:                Policy,
Chris@16:                policies::promote_float<false>,
Chris@16:                policies::promote_double<false>,
Chris@16:                policies::discrete_quantile<>,
Chris@16:                policies::assert_undefined<> >::type forwarding_policy;
Chris@16:             value_type result = detail::find_degrees_of_freedom(
Chris@16:                static_cast<value_type>(c.dist),
Chris@16:                static_cast<value_type>(c.param1),
Chris@16:                static_cast<value_type>(1-c.param2),
Chris@16:                static_cast<value_type>(c.param2),
Chris@16:                forwarding_policy());
Chris@16:             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
Chris@16:                result,
Chris@16:                function);
Chris@16:          }
Chris@16:          static RealType find_non_centrality(RealType v, RealType x, RealType p)
Chris@16:          {
Chris@16:             const char* function = "non_central_chi_squared<%1%>::find_non_centrality";
Chris@16:             typedef typename policies::evaluation<RealType, Policy>::type value_type;
Chris@16:             typedef typename policies::normalise<
Chris@16:                Policy,
Chris@16:                policies::promote_float<false>,
Chris@16:                policies::promote_double<false>,
Chris@16:                policies::discrete_quantile<>,
Chris@16:                policies::assert_undefined<> >::type forwarding_policy;
Chris@16:             value_type result = detail::find_non_centrality(
Chris@16:                static_cast<value_type>(v),
Chris@16:                static_cast<value_type>(x),
Chris@16:                static_cast<value_type>(p),
Chris@16:                static_cast<value_type>(1-p),
Chris@16:                forwarding_policy());
Chris@16:             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
Chris@16:                result,
Chris@16:                function);
Chris@16:          }
Chris@16:          template <class A, class B, class C>
Chris@16:          static RealType find_non_centrality(const complemented3_type<A,B,C>& c)
Chris@16:          {
Chris@16:             const char* function = "non_central_chi_squared<%1%>::find_non_centrality";
Chris@16:             typedef typename policies::evaluation<RealType, Policy>::type value_type;
Chris@16:             typedef typename policies::normalise<
Chris@16:                Policy,
Chris@16:                policies::promote_float<false>,
Chris@16:                policies::promote_double<false>,
Chris@16:                policies::discrete_quantile<>,
Chris@16:                policies::assert_undefined<> >::type forwarding_policy;
Chris@16:             value_type result = detail::find_non_centrality(
Chris@16:                static_cast<value_type>(c.dist),
Chris@16:                static_cast<value_type>(c.param1),
Chris@16:                static_cast<value_type>(1-c.param2),
Chris@16:                static_cast<value_type>(c.param2),
Chris@16:                forwarding_policy());
Chris@16:             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
Chris@16:                result,
Chris@16:                function);
Chris@16:          }
Chris@16:       private:
Chris@16:          // Data member, initialized by constructor.
Chris@16:          RealType df; // degrees of freedom.
Chris@16:          RealType ncp; // non-centrality parameter
Chris@16:       }; // template <class RealType, class Policy> class non_central_chi_squared_distribution
Chris@16: 
Chris@16:       typedef non_central_chi_squared_distribution<double> non_central_chi_squared; // Reserved name of type double.
Chris@16: 
Chris@16:       // Non-member functions to give properties of the distribution.
Chris@16: 
Chris@16:       template <class RealType, class Policy>
Chris@16:       inline const std::pair<RealType, RealType> range(const non_central_chi_squared_distribution<RealType, Policy>& /* dist */)
Chris@16:       { // Range of permissible values for random variable k.
Chris@16:          using boost::math::tools::max_value;
Chris@16:          return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // Max integer?
Chris@16:       }
Chris@16: 
Chris@16:       template <class RealType, class Policy>
Chris@16:       inline const std::pair<RealType, RealType> support(const non_central_chi_squared_distribution<RealType, Policy>& /* dist */)
Chris@16:       { // Range of supported values for random variable k.
Chris@16:          // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
Chris@16:          using boost::math::tools::max_value;
Chris@16:          return std::pair<RealType, RealType>(static_cast<RealType>(0),  max_value<RealType>());
Chris@16:       }
Chris@16: 
Chris@16:       template <class RealType, class Policy>
Chris@16:       inline RealType mean(const non_central_chi_squared_distribution<RealType, Policy>& dist)
Chris@16:       { // Mean of poisson distribution = lambda.
Chris@16:          const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::mean()";
Chris@16:          RealType k = dist.degrees_of_freedom();
Chris@16:          RealType l = dist.non_centrality();
Chris@16:          RealType r;
Chris@16:          if(!detail::check_df(
Chris@16:             function,
Chris@16:             k, &r, Policy())
Chris@16:             ||
Chris@16:          !detail::check_non_centrality(
Chris@16:             function,
Chris@16:             l,
Chris@16:             &r,
Chris@16:             Policy()))
Chris@16:                return r;
Chris@16:          return k + l;
Chris@16:       } // mean
Chris@16: 
Chris@16:       template <class RealType, class Policy>
Chris@16:       inline RealType mode(const non_central_chi_squared_distribution<RealType, Policy>& dist)
Chris@16:       { // mode.
Chris@16:          static const char* function = "mode(non_central_chi_squared_distribution<%1%> const&)";
Chris@16: 
Chris@16:          RealType k = dist.degrees_of_freedom();
Chris@16:          RealType l = dist.non_centrality();
Chris@16:          RealType r;
Chris@16:          if(!detail::check_df(
Chris@16:             function,
Chris@16:             k, &r, Policy())
Chris@16:             ||
Chris@16:          !detail::check_non_centrality(
Chris@16:             function,
Chris@16:             l,
Chris@16:             &r,
Chris@16:             Policy()))
Chris@16:                return (RealType)r;
Chris@16:          return detail::generic_find_mode(dist, 1 + k, function);
Chris@16:       }
Chris@16: 
Chris@16:       template <class RealType, class Policy>
Chris@16:       inline RealType variance(const non_central_chi_squared_distribution<RealType, Policy>& dist)
Chris@16:       { // variance.
Chris@16:          const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::variance()";
Chris@16:          RealType k = dist.degrees_of_freedom();
Chris@16:          RealType l = dist.non_centrality();
Chris@16:          RealType r;
Chris@16:          if(!detail::check_df(
Chris@16:             function,
Chris@16:             k, &r, Policy())
Chris@16:             ||
Chris@16:          !detail::check_non_centrality(
Chris@16:             function,
Chris@16:             l,
Chris@16:             &r,
Chris@16:             Policy()))
Chris@16:                return r;
Chris@16:          return 2 * (2 * l + k);
Chris@16:       }
Chris@16: 
Chris@16:       // RealType standard_deviation(const non_central_chi_squared_distribution<RealType, Policy>& dist)
Chris@16:       // standard_deviation provided by derived accessors.
Chris@16: 
Chris@16:       template <class RealType, class Policy>
Chris@16:       inline RealType skewness(const non_central_chi_squared_distribution<RealType, Policy>& dist)
Chris@16:       { // skewness = sqrt(l).
Chris@16:          const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::skewness()";
Chris@16:          RealType k = dist.degrees_of_freedom();
Chris@16:          RealType l = dist.non_centrality();
Chris@16:          RealType r;
Chris@16:          if(!detail::check_df(
Chris@16:             function,
Chris@16:             k, &r, Policy())
Chris@16:             ||
Chris@16:          !detail::check_non_centrality(
Chris@16:             function,
Chris@16:             l,
Chris@16:             &r,
Chris@16:             Policy()))
Chris@16:                return r;
Chris@16:          BOOST_MATH_STD_USING
Chris@16:             return pow(2 / (k + 2 * l), RealType(3)/2) * (k + 3 * l);
Chris@16:       }
Chris@16: 
Chris@16:       template <class RealType, class Policy>
Chris@16:       inline RealType kurtosis_excess(const non_central_chi_squared_distribution<RealType, Policy>& dist)
Chris@16:       {
Chris@16:          const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::kurtosis_excess()";
Chris@16:          RealType k = dist.degrees_of_freedom();
Chris@16:          RealType l = dist.non_centrality();
Chris@16:          RealType r;
Chris@16:          if(!detail::check_df(
Chris@16:             function,
Chris@16:             k, &r, Policy())
Chris@16:             ||
Chris@16:          !detail::check_non_centrality(
Chris@16:             function,
Chris@16:             l,
Chris@16:             &r,
Chris@16:             Policy()))
Chris@16:                return r;
Chris@16:          return 12 * (k + 4 * l) / ((k + 2 * l) * (k + 2 * l));
Chris@16:       } // kurtosis_excess
Chris@16: 
Chris@16:       template <class RealType, class Policy>
Chris@16:       inline RealType kurtosis(const non_central_chi_squared_distribution<RealType, Policy>& dist)
Chris@16:       {
Chris@16:          return kurtosis_excess(dist) + 3;
Chris@16:       }
Chris@16: 
Chris@16:       template <class RealType, class Policy>
Chris@16:       inline RealType pdf(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
Chris@16:       { // Probability Density/Mass Function.
Chris@16:          return detail::nccs_pdf(dist, x);
Chris@16:       } // pdf
Chris@16: 
Chris@16:       template <class RealType, class Policy>
Chris@16:       RealType cdf(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
Chris@16:       {
Chris@16:          const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::cdf(%1%)";
Chris@16:          RealType k = dist.degrees_of_freedom();
Chris@16:          RealType l = dist.non_centrality();
Chris@16:          RealType r;
Chris@16:          if(!detail::check_df(
Chris@16:             function,
Chris@16:             k, &r, Policy())
Chris@16:             ||
Chris@16:          !detail::check_non_centrality(
Chris@16:             function,
Chris@16:             l,
Chris@16:             &r,
Chris@16:             Policy())
Chris@16:             ||
Chris@16:          !detail::check_positive_x(
Chris@16:             function,
Chris@16:             x,
Chris@16:             &r,
Chris@16:             Policy()))
Chris@16:                return r;
Chris@16: 
Chris@16:          return detail::non_central_chi_squared_cdf(x, k, l, false, Policy());
Chris@16:       } // cdf
Chris@16: 
Chris@16:       template <class RealType, class Policy>
Chris@16:       RealType cdf(const complemented2_type<non_central_chi_squared_distribution<RealType, Policy>, RealType>& c)
Chris@16:       { // Complemented Cumulative Distribution Function
Chris@16:          const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::cdf(%1%)";
Chris@16:          non_central_chi_squared_distribution<RealType, Policy> const& dist = c.dist;
Chris@16:          RealType x = c.param;
Chris@16:          RealType k = dist.degrees_of_freedom();
Chris@16:          RealType l = dist.non_centrality();
Chris@16:          RealType r;
Chris@16:          if(!detail::check_df(
Chris@16:             function,
Chris@16:             k, &r, Policy())
Chris@16:             ||
Chris@16:          !detail::check_non_centrality(
Chris@16:             function,
Chris@16:             l,
Chris@16:             &r,
Chris@16:             Policy())
Chris@16:             ||
Chris@16:          !detail::check_positive_x(
Chris@16:             function,
Chris@16:             x,
Chris@16:             &r,
Chris@16:             Policy()))
Chris@16:                return r;
Chris@16: 
Chris@16:          return detail::non_central_chi_squared_cdf(x, k, l, true, Policy());
Chris@16:       } // ccdf
Chris@16: 
Chris@16:       template <class RealType, class Policy>
Chris@16:       inline RealType quantile(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& p)
Chris@16:       { // Quantile (or Percent Point) function.
Chris@16:          return detail::nccs_quantile(dist, p, false);
Chris@16:       } // quantile
Chris@16: 
Chris@16:       template <class RealType, class Policy>
Chris@16:       inline RealType quantile(const complemented2_type<non_central_chi_squared_distribution<RealType, Policy>, RealType>& c)
Chris@16:       { // Quantile (or Percent Point) function.
Chris@16:          return detail::nccs_quantile(c.dist, c.param, true);
Chris@16:       } // quantile complement.
Chris@16: 
Chris@16:    } // namespace math
Chris@16: } // namespace boost
Chris@16: 
Chris@16: // This include must be at the end, *after* the accessors
Chris@16: // for this distribution have been defined, in order to
Chris@16: // keep compilers that support two-phase lookup happy.
Chris@16: #include <boost/math/distributions/detail/derived_accessors.hpp>
Chris@16: 
Chris@16: #endif // BOOST_MATH_SPECIAL_NON_CENTRAL_CHI_SQUARE_HPP
Chris@16: 
Chris@16: 
Chris@16: