non_central_chi_squared.hpp

来自「Boost provides free peer-reviewed portab」· HPP 代码 · 共 965 行 · 第 1/3 页

// 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 = ldexp(1.0, -boost::math::policies::digits<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(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)
               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 = ldexp(1.0, -boost::math::policies::digits<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)
               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 = ldexp(1.0, -boost::math::policies::digits<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(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)
               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 = ldexp(1.0, -boost::math::policies::digits<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(k + 1, l2, pol) * gamma_p_derivative(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;
         }