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