boost/math/special_functions/detail/gamma_inva.hpp

   1 //  (C) Copyright John Maddock 2006.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 //
   7 // This is not a complete header file, it is included by gamma.hpp
   8 // after it has defined it's definitions.  This inverts the incomplete
   9 // gamma functions P and Q on the first parameter "a" using a generic
  10 // root finding algorithm (TOMS Algorithm 748).
  11 //
  12
  13 #ifndef BOOST_MATH_SP_DETAIL_GAMMA_INVA
  14 #define BOOST_MATH_SP_DETAIL_GAMMA_INVA
  15
  16 #ifdef _MSC_VER
  17 #pragma once
  18 #endif
  19
  20 #include <boost/math/tools/toms748_solve.hpp>
  21 #include <boost/cstdint.hpp>
  22
  23 namespace boost{ namespace math{ namespace detail{
  24
  25 template <class T, class Policy>
  26 struct gamma_inva_t
  27 {
  28    gamma_inva_t(T z_, T p_, bool invert_) : z(z_), p(p_), invert(invert_) {}
  29    T operator()(T a)
  30    {
  31       return invert ? p - boost::math::gamma_q(a, z, Policy()) : boost::math::gamma_p(a, z, Policy()) - p;
  32    }
  33 private:
  34    T z, p;
  35    bool invert;
  36 };
  37
  38 template <class T, class Policy>
  39 T inverse_poisson_cornish_fisher(T lambda, T p, T q, const Policy& pol)
  40 {
  41    BOOST_MATH_STD_USING
  42    // mean:
  43    T m = lambda;
  44    // standard deviation:
  45    T sigma = sqrt(lambda);
  46    // skewness
  47    T sk = 1 / sigma;
  48    // kurtosis:
  49    // T k = 1/lambda;
  50    // Get the inverse of a std normal distribution:
  51    T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p, pol) * constants::root_two<T>();
  52    // Set the sign:
  53    if(p < 0.5)
  54       x = -x;
  55    T x2 = x * x;
  56    // w is correction term due to skewness
  57    T w = x + sk * (x2 - 1) / 6;
  58    /*
  59    // Add on correction due to kurtosis.
  60    // Disabled for now, seems to make things worse?
  61    //
  62    if(lambda >= 10)
  63       w += k * x * (x2 - 3) / 24 + sk * sk * x * (2 * x2 - 5) / -36;
  64    */
  65    w = m + sigma * w;
  66    return w > tools::min_value<T>() ? w : tools::min_value<T>();
  67 }
  68
  69 template <class T, class Policy>
  70 T gamma_inva_imp(const T& z, const T& p, const T& q, const Policy& pol)
  71 {
  72    BOOST_MATH_STD_USING  // for ADL of std lib math functions
  73    //
  74    // Special cases first:
  75    //
  76    if(p == 0)
  77    {
  78       return tools::max_value<T>();
  79    }
  80    if(q == 0)
  81    {
  82       return tools::min_value<T>();
  83    }
  84    //
  85    // Function object, this is the functor whose root
  86    // we have to solve:
  87    //
  88    gamma_inva_t<T, Policy> f(z, (p < q) ? p : q, (p < q) ? false : true);
  89    //
  90    // Tolerance: full precision.
  91    //
  92    tools::eps_tolerance<T> tol(policies::digits<T, Policy>());
  93    //
  94    // Now figure out a starting guess for what a may be,
  95    // we'll start out with a value that'll put p or q
  96    // right bang in the middle of their range, the functions
  97    // are quite sensitive so we should need too many steps
  98    // to bracket the root from there:
  99    //
 100    T guess;
 101    T factor = 8;
 102    if(z >= 1)
 103    {
 104       //
 105       // We can use the relationship between the incomplete
 106       // gamma function and the poisson distribution to
 107       // calculate an approximate inverse, for large z
 108       // this is actually pretty accurate, but it fails badly
 109       // when z is very small.  Also set our step-factor according
 110       // to how accurate we think the result is likely to be:
 111       //
 112       guess = 1 + inverse_poisson_cornish_fisher(z, q, p, pol);
 113       if(z > 5)
 114       {
 115          if(z > 1000)
 116             factor = 1.01f;
 117          else if(z > 50)
 118             factor = 1.1f;
 119          else if(guess > 10)
 120             factor = 1.25f;
 121          else
 122             factor = 2;
 123          if(guess < 1.1)
 124             factor = 8;
 125       }
 126    }
 127    else if(z > 0.5)
 128    {
 129       guess = z * 1.2f;
 130    }
 131    else
 132    {
 133       guess = -0.4f / log(z);
 134    }
 135    //
 136    // Max iterations permitted:
 137    //
 138    boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
 139    //
 140    // Use our generic derivative-free root finding procedure.
 141    // We could use Newton steps here, taking the PDF of the
 142    // Poisson distribution as our derivative, but that's
 143    // even worse performance-wise than the generic method :-(
 144    //
 145    std::pair<T, T> r = bracket_and_solve_root(f, guess, factor, false, tol, max_iter, pol);
 146    if(max_iter >= policies::get_max_root_iterations<Policy>())
 147       policies::raise_evaluation_error<T>("boost::math::gamma_p_inva<%1%>(%1%, %1%)", "Unable to locate the root within a reasonable number of iterations, closest approximation so far was %1%", r.first, pol);
 148    return (r.first + r.second) / 2;
 149 }
 150
 151 } // namespace detail
 152
 153 template <class T1, class T2, class Policy>
 154 inline typename tools::promote_args<T1, T2>::type
 155    gamma_p_inva(T1 x, T2 p, const Policy& pol)
 156 {
 157    typedef typename tools::promote_args<T1, T2>::type result_type;
 158    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 159    typedef typename policies::normalise<
 160       Policy,
 161       policies::promote_float<false>,
 162       policies::promote_double<false>,
 163       policies::discrete_quantile<>,
 164       policies::assert_undefined<> >::type forwarding_policy;
 165
 166    if(p == 0)
 167    {
 168       return tools::max_value<result_type>();
 169    }
 170    if(p == 1)
 171    {
 172       return tools::min_value<result_type>();
 173    }
 174
 175    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
 176       detail::gamma_inva_imp(
 177          static_cast<value_type>(x),
 178          static_cast<value_type>(p),
 179          static_cast<value_type>(1 - static_cast<value_type>(p)),
 180          pol), "boost::math::gamma_p_inva<%1%>(%1%, %1%)");
 181 }
 182
 183 template <class T1, class T2, class Policy>
 184 inline typename tools::promote_args<T1, T2>::type
 185    gamma_q_inva(T1 x, T2 q, const Policy& pol)
 186 {
 187    typedef typename tools::promote_args<T1, T2>::type result_type;
 188    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 189    typedef typename policies::normalise<
 190       Policy,
 191       policies::promote_float<false>,
 192       policies::promote_double<false>,
 193       policies::discrete_quantile<>,
 194       policies::assert_undefined<> >::type forwarding_policy;
 195
 196    if(q == 1)
 197    {
 198       return tools::max_value<result_type>();
 199    }
 200    if(q == 0)
 201    {
 202       return tools::min_value<result_type>();
 203    }
 204
 205    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
 206       detail::gamma_inva_imp(
 207          static_cast<value_type>(x),
 208          static_cast<value_type>(1 - static_cast<value_type>(q)),
 209          static_cast<value_type>(q),
 210          pol), "boost::math::gamma_q_inva<%1%>(%1%, %1%)");
 211 }
 212
 213 template <class T1, class T2>
 214 inline typename tools::promote_args<T1, T2>::type
 215    gamma_p_inva(T1 x, T2 p)
 216 {
 217    return boost::math::gamma_p_inva(x, p, policies::policy<>());
 218 }
 219
 220 template <class T1, class T2>
 221 inline typename tools::promote_args<T1, T2>::type
 222    gamma_q_inva(T1 x, T2 q)
 223 {
 224    return boost::math::gamma_q_inva(x, q, policies::policy<>());
 225 }
 226
 227 } // namespace math
 228 } // namespace boost
 229
 230 #endif // BOOST_MATH_SP_DETAIL_GAMMA_INVA
 231
 232
 233