boost/math/special_functions/detail/igamma_inverse.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 #ifndef BOOST_MATH_SPECIAL_FUNCTIONS_IGAMMA_INVERSE_HPP
   7 #define BOOST_MATH_SPECIAL_FUNCTIONS_IGAMMA_INVERSE_HPP
   8
   9 #ifdef _MSC_VER
  10 #pragma once
  11 #endif
  12
  13 #include <boost/tr1/tuple.hpp>
  14 #include <boost/math/special_functions/gamma.hpp>
  15 #include <boost/math/special_functions/sign.hpp>
  16 #include <boost/math/tools/roots.hpp>
  17 #include <boost/math/policies/error_handling.hpp>
  18
  19 namespace boost{ namespace math{
  20
  21 namespace detail{
  22
  23 template <class T>
  24 T find_inverse_s(T p, T q)
  25 {
  26    //
  27    // Computation of the Incomplete Gamma Function Ratios and their Inverse
  28    // ARMIDO R. DIDONATO and ALFRED H. MORRIS, JR.
  29    // ACM Transactions on Mathematical Software, Vol. 12, No. 4,
  30    // December 1986, Pages 377-393.
  31    //
  32    // See equation 32.
  33    //
  34    BOOST_MATH_STD_USING
  35    T t;
  36    if(p < 0.5)
  37    {
  38       t = sqrt(-2 * log(p));
  39    }
  40    else
  41    {
  42       t = sqrt(-2 * log(q));
  43    }
  44    static const double a[4] = { 3.31125922108741, 11.6616720288968, 4.28342155967104, 0.213623493715853 };
  45    static const double b[5] = { 1, 6.61053765625462, 6.40691597760039, 1.27364489782223, 0.3611708101884203e-1 };
  46    T s = t - tools::evaluate_polynomial(a, t) / tools::evaluate_polynomial(b, t);
  47    if(p < 0.5)
  48       s = -s;
  49    return s;
  50 }
  51
  52 template <class T>
  53 T didonato_SN(T a, T x, unsigned N, T tolerance = 0)
  54 {
  55    //
  56    // Computation of the Incomplete Gamma Function Ratios and their Inverse
  57    // ARMIDO R. DIDONATO and ALFRED H. MORRIS, JR.
  58    // ACM Transactions on Mathematical Software, Vol. 12, No. 4,
  59    // December 1986, Pages 377-393.
  60    //
  61    // See equation 34.
  62    //
  63    T sum = 1;
  64    if(N >= 1)
  65    {
  66       T partial = x / (a + 1);
  67       sum += partial;
  68       for(unsigned i = 2; i <= N; ++i)
  69       {
  70          partial *= x / (a + i);
  71          sum += partial;
  72          if(partial < tolerance)
  73             break;
  74       }
  75    }
  76    return sum;
  77 }
  78
  79 template <class T, class Policy>
  80 inline T didonato_FN(T p, T a, T x, unsigned N, T tolerance, const Policy& pol)
  81 {
  82    //
  83    // Computation of the Incomplete Gamma Function Ratios and their Inverse
  84    // ARMIDO R. DIDONATO and ALFRED H. MORRIS, JR.
  85    // ACM Transactions on Mathematical Software, Vol. 12, No. 4,
  86    // December 1986, Pages 377-393.
  87    //
  88    // See equation 34.
  89    //
  90    BOOST_MATH_STD_USING
  91    T u = log(p) + boost::math::lgamma(a + 1, pol);
  92    return exp((u + x - log(didonato_SN(a, x, N, tolerance))) / a);
  93 }
  94
  95 template <class T, class Policy>
  96 T find_inverse_gamma(T a, T p, T q, const Policy& pol, bool* p_has_10_digits)
  97 {
  98    //
  99    // In order to understand what's going on here, you will
 100    // need to refer to:
 101    //
 102    // Computation of the Incomplete Gamma Function Ratios and their Inverse
 103    // ARMIDO R. DIDONATO and ALFRED H. MORRIS, JR.
 104    // ACM Transactions on Mathematical Software, Vol. 12, No. 4,
 105    // December 1986, Pages 377-393.
 106    //
 107    BOOST_MATH_STD_USING
 108
 109    T result;
 110    *p_has_10_digits = false;
 111
 112    if(a == 1)
 113    {
 114       result = -log(q);
 115       BOOST_MATH_INSTRUMENT_VARIABLE(result);
 116    }
 117    else if(a < 1)
 118    {
 119       T g = boost::math::tgamma(a, pol);
 120       T b = q * g;
 121       BOOST_MATH_INSTRUMENT_VARIABLE(g);
 122       BOOST_MATH_INSTRUMENT_VARIABLE(b);
 123       if((b > 0.6) || ((b >= 0.45) && (a >= 0.3)))
 124       {
 125          // DiDonato & Morris Eq 21:
 126          //
 127          // There is a slight variation from DiDonato and Morris here:
 128          // the first form given here is unstable when p is close to 1,
 129          // making it impossible to compute the inverse of Q(a,x) for small
 130          // q.  Fortunately the second form works perfectly well in this case.
 131          //
 132          T u;
 133          if((b * q > 1e-8) && (q > 1e-5))
 134          {
 135             u = pow(p * g * a, 1 / a);
 136             BOOST_MATH_INSTRUMENT_VARIABLE(u);
 137          }
 138          else
 139          {
 140             u = exp((-q / a) - constants::euler<T>());
 141             BOOST_MATH_INSTRUMENT_VARIABLE(u);
 142          }
 143          result = u / (1 - (u / (a + 1)));
 144          BOOST_MATH_INSTRUMENT_VARIABLE(result);
 145       }
 146       else if((a < 0.3) && (b >= 0.35))
 147       {
 148          // DiDonato & Morris Eq 22:
 149          T t = exp(-constants::euler<T>() - b);
 150          T u = t * exp(t);
 151          result = t * exp(u);
 152          BOOST_MATH_INSTRUMENT_VARIABLE(result);
 153       }
 154       else if((b > 0.15) || (a >= 0.3))
 155       {
 156          // DiDonato & Morris Eq 23:
 157          T y = -log(b);
 158          T u = y - (1 - a) * log(y);
 159          result = y - (1 - a) * log(u) - log(1 + (1 - a) / (1 + u));
 160          BOOST_MATH_INSTRUMENT_VARIABLE(result);
 161       }
 162       else if (b > 0.1)
 163       {
 164          // DiDonato & Morris Eq 24:
 165          T y = -log(b);
 166          T u = y - (1 - a) * log(y);
 167          result = y - (1 - a) * log(u) - log((u * u + 2 * (3 - a) * u + (2 - a) * (3 - a)) / (u * u + (5 - a) * u + 2));
 168          BOOST_MATH_INSTRUMENT_VARIABLE(result);
 169       }
 170       else
 171       {
 172          // DiDonato & Morris Eq 25:
 173          T y = -log(b);
 174          T c1 = (a - 1) * log(y);
 175          T c1_2 = c1 * c1;
 176          T c1_3 = c1_2 * c1;
 177          T c1_4 = c1_2 * c1_2;
 178          T a_2 = a * a;
 179          T a_3 = a_2 * a;
 180
 181          T c2 = (a - 1) * (1 + c1);
 182          T c3 = (a - 1) * (-(c1_2 / 2) + (a - 2) * c1 + (3 * a - 5) / 2);
 183          T c4 = (a - 1) * ((c1_3 / 3) - (3 * a - 5) * c1_2 / 2 + (a_2 - 6 * a + 7) * c1 + (11 * a_2 - 46 * a + 47) / 6);
 184          T c5 = (a - 1) * (-(c1_4 / 4)
 185                            + (11 * a - 17) * c1_3 / 6
 186                            + (-3 * a_2 + 13 * a -13) * c1_2
 187                            + (2 * a_3 - 25 * a_2 + 72 * a - 61) * c1 / 2
 188                            + (25 * a_3 - 195 * a_2 + 477 * a - 379) / 12);
 189
 190          T y_2 = y * y;
 191          T y_3 = y_2 * y;
 192          T y_4 = y_2 * y_2;
 193          result = y + c1 + (c2 / y) + (c3 / y_2) + (c4 / y_3) + (c5 / y_4);
 194          BOOST_MATH_INSTRUMENT_VARIABLE(result);
 195          if(b < 1e-28f)
 196             *p_has_10_digits = true;
 197       }
 198    }
 199    else
 200    {
 201       // DiDonato and Morris Eq 31:
 202       T s = find_inverse_s(p, q);
 203
 204       BOOST_MATH_INSTRUMENT_VARIABLE(s);
 205
 206       T s_2 = s * s;
 207       T s_3 = s_2 * s;
 208       T s_4 = s_2 * s_2;
 209       T s_5 = s_4 * s;
 210       T ra = sqrt(a);
 211
 212       BOOST_MATH_INSTRUMENT_VARIABLE(ra);
 213
 214       T w = a + s * ra + (s * s -1) / 3;
 215       w += (s_3 - 7 * s) / (36 * ra);
 216       w -= (3 * s_4 + 7 * s_2 - 16) / (810 * a);
 217       w += (9 * s_5 + 256 * s_3 - 433 * s) / (38880 * a * ra);
 218
 219       BOOST_MATH_INSTRUMENT_VARIABLE(w);
 220
 221       if((a >= 500) && (fabs(1 - w / a) < 1e-6))
 222       {
 223          result = w;
 224          *p_has_10_digits = true;
 225          BOOST_MATH_INSTRUMENT_VARIABLE(result);
 226       }
 227       else if (p > 0.5)
 228       {
 229          if(w < 3 * a)
 230          {
 231             result = w;
 232             BOOST_MATH_INSTRUMENT_VARIABLE(result);
 233          }
 234          else
 235          {
 236             T D = (std::max)(T(2), T(a * (a - 1)));
 237             T lg = boost::math::lgamma(a, pol);
 238             T lb = log(q) + lg;
 239             if(lb < -D * 2.3)
 240             {
 241                // DiDonato and Morris Eq 25:
 242                T y = -lb;
 243                T c1 = (a - 1) * log(y);
 244                T c1_2 = c1 * c1;
 245                T c1_3 = c1_2 * c1;
 246                T c1_4 = c1_2 * c1_2;
 247                T a_2 = a * a;
 248                T a_3 = a_2 * a;
 249
 250                T c2 = (a - 1) * (1 + c1);
 251                T c3 = (a - 1) * (-(c1_2 / 2) + (a - 2) * c1 + (3 * a - 5) / 2);
 252                T c4 = (a - 1) * ((c1_3 / 3) - (3 * a - 5) * c1_2 / 2 + (a_2 - 6 * a + 7) * c1 + (11 * a_2 - 46 * a + 47) / 6);
 253                T c5 = (a - 1) * (-(c1_4 / 4)
 254                                  + (11 * a - 17) * c1_3 / 6
 255                                  + (-3 * a_2 + 13 * a -13) * c1_2
 256                                  + (2 * a_3 - 25 * a_2 + 72 * a - 61) * c1 / 2
 257                                  + (25 * a_3 - 195 * a_2 + 477 * a - 379) / 12);
 258
 259                T y_2 = y * y;
 260                T y_3 = y_2 * y;
 261                T y_4 = y_2 * y_2;
 262                result = y + c1 + (c2 / y) + (c3 / y_2) + (c4 / y_3) + (c5 / y_4);
 263                BOOST_MATH_INSTRUMENT_VARIABLE(result);
 264             }
 265             else
 266             {
 267                // DiDonato and Morris Eq 33:
 268                T u = -lb + (a - 1) * log(w) - log(1 + (1 - a) / (1 + w));
 269                result = -lb + (a - 1) * log(u) - log(1 + (1 - a) / (1 + u));
 270                BOOST_MATH_INSTRUMENT_VARIABLE(result);
 271             }
 272          }
 273       }
 274       else
 275       {
 276          T z = w;
 277          T ap1 = a + 1;
 278          if(w < 0.15f * ap1)
 279          {
 280             // DiDonato and Morris Eq 35:
 281             T v = log(p) + boost::math::lgamma(ap1, pol);
 282             T s = 1;
 283             z = exp((v + w) / a);
 284             s = boost::math::log1p(z / ap1 * (1 + z / (a + 2)));
 285             z = exp((v + z - s) / a);
 286             z = exp((v + z - s) / a);
 287             s = boost::math::log1p(z / ap1 * (1 + z / (a + 2) * (1 + z / (a + 3))));
 288             z = exp((v + z - s) / a);
 289             BOOST_MATH_INSTRUMENT_VARIABLE(z);
 290          }
 291
 292          if((z <= 0.01 * ap1) || (z > 0.7 * ap1))
 293          {
 294             result = z;
 295             if(z <= 0.002 * ap1)
 296                *p_has_10_digits = true;
 297             BOOST_MATH_INSTRUMENT_VARIABLE(result);
 298          }
 299          else
 300          {
 301             // DiDonato and Morris Eq 36:
 302             T ls = log(didonato_SN(a, z, 100, T(1e-4)));
 303             T v = log(p) + boost::math::lgamma(ap1, pol);
 304             z = exp((v + z - ls) / a);
 305             result = z * (1 - (a * log(z) - z - v + ls) / (a - z));
 306
 307             BOOST_MATH_INSTRUMENT_VARIABLE(result);
 308          }
 309       }
 310    }
 311    return result;
 312 }
 313
 314 template <class T, class Policy>
 315 struct gamma_p_inverse_func
 316 {
 317    gamma_p_inverse_func(T a_, T p_, bool inv) : a(a_), p(p_), invert(inv)
 318    {
 319       //
 320       // If p is too near 1 then P(x) - p suffers from cancellation
 321       // errors causing our root-finding algorithms to "thrash", better
 322       // to invert in this case and calculate Q(x) - (1-p) instead.
 323       //
 324       // Of course if p is *very* close to 1, then the answer we get will
 325       // be inaccurate anyway (because there's not enough information in p)
 326       // but at least we will converge on the (inaccurate) answer quickly.
 327       //
 328       if(p > 0.9)
 329       {
 330          p = 1 - p;
 331          invert = !invert;
 332       }
 333    }
 334
 335    std::tr1::tuple<T, T, T> operator()(const T& x)const
 336    {
 337       BOOST_FPU_EXCEPTION_GUARD
 338       //
 339       // Calculate P(x) - p and the first two derivates, or if the invert
 340       // flag is set, then Q(x) - q and it's derivatives.
 341       //
 342       typedef typename policies::evaluation<T, Policy>::type value_type;
 343       typedef typename lanczos::lanczos<T, Policy>::type evaluation_type;
 344       typedef typename policies::normalise<
 345          Policy,
 346          policies::promote_float<false>,
 347          policies::promote_double<false>,
 348          policies::discrete_quantile<>,
 349          policies::assert_undefined<> >::type forwarding_policy;
 350
 351       BOOST_MATH_STD_USING  // For ADL of std functions.
 352
 353       T f, f1;
 354       value_type ft;
 355       f = static_cast<T>(boost::math::detail::gamma_incomplete_imp(
 356                static_cast<value_type>(a),
 357                static_cast<value_type>(x),
 358                true, invert,
 359                forwarding_policy(), &ft));
 360       f1 = static_cast<T>(ft);
 361       T f2;
 362       T div = (a - x - 1) / x;
 363       f2 = f1;
 364       if((fabs(div) > 1) && (tools::max_value<T>() / fabs(div) < f2))
 365       {
 366          // overflow:
 367          f2 = -tools::max_value<T>() / 2;
 368       }
 369       else
 370       {
 371          f2 *= div;
 372       }
 373
 374       if(invert)
 375       {
 376          f1 = -f1;
 377          f2 = -f2;
 378       }
 379
 380       return std::tr1::make_tuple(f - p, f1, f2);
 381    }
 382 private:
 383    T a, p;
 384    bool invert;
 385 };
 386
 387 template <class T, class Policy>
 388 T gamma_p_inv_imp(T a, T p, const Policy& pol)
 389 {
 390    BOOST_MATH_STD_USING  // ADL of std functions.
 391
 392    static const char* function = "boost::math::gamma_p_inv<%1%>(%1%, %1%)";
 393
 394    BOOST_MATH_INSTRUMENT_VARIABLE(a);
 395    BOOST_MATH_INSTRUMENT_VARIABLE(p);
 396
 397    if(a <= 0)
 398       policies::raise_domain_error<T>(function, "Argument a in the incomplete gamma function inverse must be >= 0 (got a=%1%).", a, pol);
 399    if((p < 0) || (p > 1))
 400       policies::raise_domain_error<T>(function, "Probabilty must be in the range [0,1] in the incomplete gamma function inverse (got p=%1%).", p, pol);
 401    if(p == 1)
 402       return tools::max_value<T>();
 403    if(p == 0)
 404       return 0;
 405    bool has_10_digits;
 406    T guess = detail::find_inverse_gamma<T>(a, p, 1 - p, pol, &has_10_digits);
 407    if((policies::digits<T, Policy>() <= 36) && has_10_digits)
 408       return guess;
 409    T lower = tools::min_value<T>();
 410    if(guess <= lower)
 411       guess = tools::min_value<T>();
 412    BOOST_MATH_INSTRUMENT_VARIABLE(guess);
 413    //
 414    // Work out how many digits to converge to, normally this is
 415    // 2/3 of the digits in T, but if the first derivative is very
 416    // large convergence is slow, so we'll bump it up to full
 417    // precision to prevent premature termination of the root-finding routine.
 418    //
 419    unsigned digits = policies::digits<T, Policy>();
 420    if(digits < 30)
 421    {
 422       digits *= 2;
 423       digits /= 3;
 424    }
 425    else
 426    {
 427       digits /= 2;
 428       digits -= 1;
 429    }
 430    if((a < 0.125) && (fabs(gamma_p_derivative(a, guess, pol)) > 1 / sqrt(tools::epsilon<T>())))
 431       digits = policies::digits<T, Policy>() - 2;
 432    //
 433    // Go ahead and iterate:
 434    //
 435    boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
 436    guess = tools::halley_iterate(
 437       detail::gamma_p_inverse_func<T, Policy>(a, p, false),
 438       guess,
 439       lower,
 440       tools::max_value<T>(),
 441       digits,
 442       max_iter);
 443    policies::check_root_iterations(function, max_iter, pol);
 444    BOOST_MATH_INSTRUMENT_VARIABLE(guess);
 445    if(guess == lower)
 446       guess = policies::raise_underflow_error<T>(function, "Expected result known to be non-zero, but is smaller than the smallest available number.", pol);
 447    return guess;
 448 }
 449
 450 template <class T, class Policy>
 451 T gamma_q_inv_imp(T a, T q, const Policy& pol)
 452 {
 453    BOOST_MATH_STD_USING  // ADL of std functions.
 454
 455    static const char* function = "boost::math::gamma_q_inv<%1%>(%1%, %1%)";
 456
 457    if(a <= 0)
 458       policies::raise_domain_error<T>(function, "Argument a in the incomplete gamma function inverse must be >= 0 (got a=%1%).", a, pol);
 459    if((q < 0) || (q > 1))
 460       policies::raise_domain_error<T>(function, "Probabilty must be in the range [0,1] in the incomplete gamma function inverse (got q=%1%).", q, pol);
 461    if(q == 0)
 462       return tools::max_value<T>();
 463    if(q == 1)
 464       return 0;
 465    bool has_10_digits;
 466    T guess = detail::find_inverse_gamma<T>(a, 1 - q, q, pol, &has_10_digits);
 467    if((policies::digits<T, Policy>() <= 36) && has_10_digits)
 468       return guess;
 469    T lower = tools::min_value<T>();
 470    if(guess <= lower)
 471       guess = tools::min_value<T>();
 472    //
 473    // Work out how many digits to converge to, normally this is
 474    // 2/3 of the digits in T, but if the first derivative is very
 475    // large convergence is slow, so we'll bump it up to full
 476    // precision to prevent premature termination of the root-finding routine.
 477    //
 478    unsigned digits = policies::digits<T, Policy>();
 479    if(digits < 30)
 480    {
 481       digits *= 2;
 482       digits /= 3;
 483    }
 484    else
 485    {
 486       digits /= 2;
 487       digits -= 1;
 488    }
 489    if((a < 0.125) && (fabs(gamma_p_derivative(a, guess, pol)) > 1 / sqrt(tools::epsilon<T>())))
 490       digits = policies::digits<T, Policy>();
 491    //
 492    // Go ahead and iterate:
 493    //
 494    boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
 495    guess = tools::halley_iterate(
 496       detail::gamma_p_inverse_func<T, Policy>(a, q, true),
 497       guess,
 498       lower,
 499       tools::max_value<T>(),
 500       digits,
 501       max_iter);
 502    policies::check_root_iterations(function, max_iter, pol);
 503    if(guess == lower)
 504       guess = policies::raise_underflow_error<T>(function, "Expected result known to be non-zero, but is smaller than the smallest available number.", pol);
 505    return guess;
 506 }
 507
 508 } // namespace detail
 509
 510 template <class T1, class T2, class Policy>
 511 inline typename tools::promote_args<T1, T2>::type
 512    gamma_p_inv(T1 a, T2 p, const Policy& pol)
 513 {
 514    typedef typename tools::promote_args<T1, T2>::type result_type;
 515    return detail::gamma_p_inv_imp(
 516       static_cast<result_type>(a),
 517       static_cast<result_type>(p), pol);
 518 }
 519
 520 template <class T1, class T2, class Policy>
 521 inline typename tools::promote_args<T1, T2>::type
 522    gamma_q_inv(T1 a, T2 p, const Policy& pol)
 523 {
 524    typedef typename tools::promote_args<T1, T2>::type result_type;
 525    return detail::gamma_q_inv_imp(
 526       static_cast<result_type>(a),
 527       static_cast<result_type>(p), pol);
 528 }
 529
 530 template <class T1, class T2>
 531 inline typename tools::promote_args<T1, T2>::type
 532    gamma_p_inv(T1 a, T2 p)
 533 {
 534    return gamma_p_inv(a, p, policies::policy<>());
 535 }
 536
 537 template <class T1, class T2>
 538 inline typename tools::promote_args<T1, T2>::type
 539    gamma_q_inv(T1 a, T2 p)
 540 {
 541    return gamma_q_inv(a, p, policies::policy<>());
 542 }
 543
 544 } // namespace math
 545 } // namespace boost
 546
 547 #endif // BOOST_MATH_SPECIAL_FUNCTIONS_IGAMMA_INVERSE_HPP
 548
 549
 550