boost/math/distributions/non_central_chi_squared.hpp
// 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 llround #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 coefficient. // 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; std::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. // long long k = llround(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: // long long i; for(i = k; static_cast<std::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<std::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); // LCOV_EXCL_LINE // // Now backwards iteration: the gamma // function recurrences are unstable in this // direction, we rely on the terms diminishing 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 underflow 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; std::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<std::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<std::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); // LCOV_EXCL_LINE 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 coefficient. // 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; std::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. // long long k = llround(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<std::uintmax_t>(i) < max_iter)); //Error check: if(static_cast<std::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); // LCOV_EXCL_LINE 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 std::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; long long k = lltrunc(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(long long i = k; ; ++i) { sum += pois; if(pois / sum < errtol) break; if(static_cast<std::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); // LCOV_EXCL_LINE pois *= l2 * x2 / ((i + 1) * (n2 + i)); } for(long long 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 static_cast<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 static_cast<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%", // LCOV_EXCL_LINE RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy()); // LCOV_EXCL_LINE } 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>()); std::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:" // LCOV_EXCL_LINE " or there is no answer to problem. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE } 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%", // LCOV_EXCL_LINE RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy()); // LCOV_EXCL_LINE } 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>()); std::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:" // LCOV_EXCL_LINE " or there is no answer to problem. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE } 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. #ifdef __cpp_deduction_guides template <class RealType> non_central_chi_squared_distribution(RealType,RealType)->non_central_chi_squared_distribution<typename boost::math::tools::promote_args<RealType>::type>; #endif // 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 static_cast<RealType>(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 static_cast<RealType>(r); bool asymptotic_mode = k < l/4; RealType starting_point = asymptotic_mode ? k + l - RealType(3) : RealType(1) + k; return detail::generic_find_mode(dist, starting_point, 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 static_cast<RealType>(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 static_cast<RealType>(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 static_cast<RealType>(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 static_cast<RealType>(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 static_cast<RealType>(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