non_central_t.hpp

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

// boost\math\distributions\non_central_t.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_T_HPP
#define BOOST_MATH_SPECIAL_NON_CENTRAL_T_HPP

#include <boost/math/distributions/fwd.hpp>
#include <boost/math/distributions/non_central_beta.hpp> // for nc beta
#include <boost/math/distributions/normal.hpp> // for normal CDF and quantile
#include <boost/math/distributions/students_t.hpp>
#include <boost/math/distributions/detail/generic_quantile.hpp> // quantile

namespace boost
{
   namespace math
   {

      template <class RealType, class Policy>
      class non_central_t_distribution;

      namespace detail{

         template <class T, class Policy>
         T non_central_t2_p(T n, T delta, T x, T y, const Policy& pol, T init_val)
         {
            BOOST_MATH_STD_USING
            //
            // Variables come first:
            //
            boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
            T errtol = ldexp(1.0f, -boost::math::policies::digits<T, Policy>());
            T d2 = delta * delta / 2;
            //
            // k is the starting point for iteration, and is the
            // maximum of the poisson weighting term:
            //
            int k = boost::math::itrunc(d2);
            // Starting Poisson weight:
            T pois = gamma_p_derivative(T(k+1), d2, pol) 
               * tgamma_delta_ratio(T(k + 1), T(0.5f))
               * delta / constants::root_two<T>();
            if(pois == 0)
               return init_val;
            // Starting beta term:
            T beta = x < y
               ? ibeta(T(k + 1), n / 2, x, pol)
               : ibetac(n / 2, T(k + 1), y, pol);
            // Recurance term:
            T xterm = x < y
               ? ibeta_derivative(T(k + 1), n / 2, x, pol)
               : ibeta_derivative(n / 2, T(k + 1), y, pol);
            xterm *= y / (n / 2 + k);
            T poisf(pois), betaf(beta), xtermf(xterm);
            T sum = init_val;
            if((xterm == 0) && (beta == 0))
               return init_val;

            //
            // Backwards recursion first, this is the stable
            // direction for recursion:
            //
            boost::uintmax_t count = 0;
            for(int i = k; i >= 0; --i)
            {
               T term = beta * pois;
               sum += term;
               if(fabs(term/sum) < errtol)
                  break;
               pois *= (i + 0.5f) / d2;
               beta += xterm;
               xterm *= (i) / (x * (n / 2 + i - 1));
               ++count;
            }
            for(int i = k + 1; ; ++i)
            {
               poisf *= d2 / (i + 0.5f);
               xtermf *= (x * (n / 2 + i - 1)) / (i);
               betaf -= xtermf;
               T term = poisf * betaf;
               sum += term;
               if(fabs(term/sum) < errtol)
                  break;
               ++count;
               if(count > max_iter)
               {
                  return policies::raise_evaluation_error(
                     "cdf(non_central_t_distribution<%1%>, %1%)", 
                     "Series did not converge, closest value was %1%", sum, pol);
               }
            }
            return sum;
         }

         template <class T, class Policy>
         T non_central_t2_q(T n, T delta, T x, T y, const Policy& pol, T init_val)
         {
            BOOST_MATH_STD_USING
            //
            // Variables come first:
            //
            boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
            T errtol = ldexp(1.0f, -boost::math::policies::digits<T, Policy>());
            T d2 = delta * delta / 2;
            //
            // k is the starting point for iteration, and is the
            // maximum of the poisson weighting term:
            //
            int k = boost::math::itrunc(d2);
            // Starting Poisson weight:
            T pois = gamma_p_derivative(T(k+1), d2, pol) 
               * tgamma_delta_ratio(T(k + 1), T(0.5f))
               * delta / constants::root_two<T>();
            if(pois == 0)
               return init_val;
            // Starting beta term:
            T beta = x < y 
               ? ibetac(T(k + 1), n / 2, x, pol)
               : ibeta(n / 2, T(k + 1), y, pol);
            // Recurance term:
            T xterm = x < y
               ? ibeta_derivative(T(k + 1), n / 2, x, pol)
               : ibeta_derivative(n / 2, T(k + 1), y, pol);
            xterm *= y / (n / 2 + k);
            T poisf(pois), betaf(beta), xtermf(xterm);
            T sum = init_val;
            if((xterm == 0) && (beta == 0))
               return init_val;

            //
            // Forward recursion first, this is the stable direction:
            //
            boost::uintmax_t count = 0;
            for(int i = k + 1; ; ++i)
            {
               poisf *= d2 / (i + 0.5f);
               xtermf *= (x * (n / 2 + i - 1)) / (i);
               betaf += xtermf;

               T term = poisf * betaf;
               sum += term;
               if(fabs(term/sum) < errtol)
                  break;
               if(count > max_iter)
               {
                  return policies::raise_evaluation_error(
                     "cdf(non_central_t_distribution<%1%>, %1%)", 
                     "Series did not converge, closest value was %1%", sum, pol);
               }
               ++count;
            }
            //
            // Backwards recursion:
            //
            for(int i = k; i >= 0; --i)
            {
               T term = beta * pois;
               sum += term;
               if(fabs(term/sum) < errtol)
                  break;
               pois *= (i + 0.5f) / d2;
               beta -= xterm;
               xterm *= (i) / (x * (n / 2 + i - 1));
               ++count;
               if(count > max_iter)
               {
                  return policies::raise_evaluation_error(
                     "cdf(non_central_t_distribution<%1%>, %1%)", 
                     "Series did not converge, closest value was %1%", sum, pol);
               }
            }
            return sum;
         }

         template <class T, class Policy>
         T non_central_t_cdf(T n, T delta, T t, bool invert, const Policy& pol)
         {
            //
            // For t < 0 we have to use reflect:
            //
            if(t < 0)
            {
               t = -t;
               delta = -delta;
               invert = !invert;
            }
            //
            // x and y are the corresponding random
            // variables for the noncentral beta distribution,
            // with y = 1 - x:
            //
            T x = t * t / (n + t * t);
            T y = n / (n + t * t);
            T d2 = delta * delta;
            T a = 0.5f;
            T b = n / 2;
            T c = a + b + d2 / 2;
            //
            // Crossover point for calculating p or q is the same
            // as for the noncentral beta:
            //
            T cross = 1 - (b / c) * (1 + d2 / (2 * c * c));
            T result;
            if(x < cross)
            {
               //
               // Calculate p:
               //
               if(x != 0)
               {
                  result = non_central_beta_p(a, b, d2, x, y, pol);
                  result = non_central_t2_p(n, delta, x, y, pol, result);
                  result /= 2;
               }
               else
                  result = 0;
               result += cdf(boost::math::normal_distribution<T, Policy>(), -delta);
            }
            else
            {
               //
               // Calculate q:
               //
               invert = !invert;
               if(x != 0)
               {
                  result = non_central_beta_q(a, b, d2, x, y, pol);
                  result = non_central_t2_q(n, delta, x, y, pol, result);
                  result /= 2;
               }
               else
                  result = cdf(complement(boost::math::normal_distribution<T, Policy>(), -delta));
            }
            if(invert)
               result = 1 - result;
            return result;
         }

         template <class T, class Policy>
         T non_central_t_quantile(T v, T delta, T p, T q, const Policy&)
         {
            BOOST_MATH_STD_USING
            static const char* function = "quantile(non_central_t_distribution<%1%>, %1%)";
            typedef typename policies::evaluation<T, 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;

               T r;
               if(!detail::check_df(
                  function,
                  v, &r, Policy())
                  ||
               !detail::check_finite(
                  function,
                  delta,
                  &r,
                  Policy())
                  ||
               !detail::check_probability(
                  function,
                  p,
                  &r,
                  Policy()))
                     return r;

            value_type guess = 0;
            if(v > 3)
            {
               value_type mean = delta * sqrt(v / 2) * tgamma_delta_ratio((v - 1) * 0.5f, T(0.5f));
               value_type var = ((delta * delta + 1) * v) / (v - 2) - mean * mean;
               if(p < q)
                  guess = quantile(normal_distribution<value_type, forwarding_policy>(mean, var), p);
               else
                  guess = quantile(complement(normal_distribution<value_type, forwarding_policy>(mean, var), q));
            }
            //
            // We *must* get the sign of the initial guess correct, 
            // or our root-finder will fail, so double check it now:
            //
            value_type pzero = non_central_t_cdf(
               static_cast<value_type>(v), 
               static_cast<value_type>(delta), 
               static_cast<value_type>(0), 
               !(p < q), 
               forwarding_policy());
            int s;
            if(p < q)
               s = boost::math::sign(p - pzero);
            else
               s = boost::math::sign(pzero - q);
            if(s != boost::math::sign(guess))
            {
               guess = s;
            }

            value_type result = detail::generic_quantile(
               non_central_t_distribution<value_type, forwarding_policy>(v, delta), 
               (p < q ? p : q), 
               guess, 
               (p >= q), 
               function);
            return policies::checked_narrowing_cast<T, forwarding_policy>(
               result, 
               function);
         }

         template <class T, class Policy>
         T non_central_t2_pdf(T n, T delta, T x, T y, const Policy& pol, T init_val)
         {
            BOOST_MATH_STD_USING
            //
            // Variables come first:
            //
            boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
            T errtol = ldexp(1.0f, -boost::math::policies::digits<T, Policy>());
            T d2 = delta * delta / 2;
            //
            // k is the starting point for iteration, and is the
            // maximum of the poisson weighting term:
            //
            int k = boost::math::itrunc(d2);
            // Starting Poisson weight:
            T pois = gamma_p_derivative(T(k+1), d2, pol) 
               * tgamma_delta_ratio(T(k + 1), T(0.5f))
               * delta / constants::root_two<T>();
            // Starting beta term:
            T xterm = x < y
               ? ibeta_derivative(T(k + 1), n / 2, x, pol)
               : ibeta_derivative(n / 2, T(k + 1), y, pol);
            T poisf(pois), xtermf(xterm);
            T sum = init_val;
            if((pois == 0) || (xterm == 0))
               return init_val;

            //
            // Backwards recursion first, this is the stable
            // direction for recursion:
            //
            boost::uintmax_t count = 0;
            for(int i = k; i >= 0; --i)
            {
               T term = xterm * pois;
               sum += term;
               if((fabs(term/sum) < errtol) || (term == 0))
                  break;
               pois *= (i + 0.5f) / d2;
               xterm *= (i) / (x * (n / 2 + i));