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/math/tools/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          T ap2 = a + 2;
 279          if(w < 0.15f * ap1)
 280          {
 281             // DiDonato and Morris Eq 35:
 282             T v = log(p) + boost::math::lgamma(ap1, pol);
 283             z = exp((v + w) / a);
 284             s = boost::math::log1p(z / ap1 * (1 + z / ap2));
 285             z = exp((v + z - s) / a);
 286             s = boost::math::log1p(z / ap1 * (1 + z / ap2));
 287             z = exp((v + z - s) / a);
 288             s = boost::math::log1p(z / ap1 * (1 + z / ap2 * (1 + z / (a + 3))));
 289             z = exp((v + z - s) / a);
 290             BOOST_MATH_INSTRUMENT_VARIABLE(z);
 291          }
 292
 293          if((z <= 0.01 * ap1) || (z > 0.7 * ap1))
 294          {
 295             result = z;
 296             if(z <= 0.002 * ap1)
 297                *p_has_10_digits = true;
 298             BOOST_MATH_INSTRUMENT_VARIABLE(result);
 299          }
 300          else
 301          {
 302             // DiDonato and Morris Eq 36:
 303             T ls = log(didonato_SN(a, z, 100, T(1e-4)));
 304             T v = log(p) + boost::math::lgamma(ap1, pol);
 305             z = exp((v + z - ls) / a);
 306             result = z * (1 - (a * log(z) - z - v + ls) / (a - z));
 307
 308             BOOST_MATH_INSTRUMENT_VARIABLE(result);
 309          }
 310       }
 311    }
 312    return result;
 313 }
 314
 315 template <class T, class Policy>
 316 struct gamma_p_inverse_func
 317 {
 318    gamma_p_inverse_func(T a_, T p_, bool inv) : a(a_), p(p_), invert(inv)
 319    {
 320       //
 321       // If p is too near 1 then P(x) - p suffers from cancellation
 322       // errors causing our root-finding algorithms to "thrash", better
 323       // to invert in this case and calculate Q(x) - (1-p) instead.
 324       //
 325       // Of course if p is *very* close to 1, then the answer we get will
 326       // be inaccurate anyway (because there's not enough information in p)
 327       // but at least we will converge on the (inaccurate) answer quickly.
 328       //
 329       if(p > 0.9)
 330       {
 331          p = 1 - p;
 332          invert = !invert;
 333       }
 334    }
 335
 336    boost::math::tuple<T, T, T> operator()(const T& x)const
 337    {
 338       BOOST_FPU_EXCEPTION_GUARD
 339       //
 340       // Calculate P(x) - p and the first two derivates, or if the invert
 341       // flag is set, then Q(x) - q and it's derivatives.
 342       //
 343       typedef typename policies::evaluation<T, Policy>::type value_type;
 344       // typedef typename lanczos::lanczos<T, Policy>::type evaluation_type;
 345       typedef typename policies::normalise<
 346          Policy,
 347          policies::promote_float<false>,
 348          policies::promote_double<false>,
 349          policies::discrete_quantile<>,
 350          policies::assert_undefined<> >::type forwarding_policy;
 351
 352       BOOST_MATH_STD_USING  // For ADL of std functions.
 353
 354       T f, f1;
 355       value_type ft;
 356       f = static_cast<T>(boost::math::detail::gamma_incomplete_imp(
 357                static_cast<value_type>(a),
 358                static_cast<value_type>(x),
 359                true, invert,
 360                forwarding_policy(), &ft));
 361       f1 = static_cast<T>(ft);
 362       T f2;
 363       T div = (a - x - 1) / x;
 364       f2 = f1;
 365       if((fabs(div) > 1) && (tools::max_value<T>() / fabs(div) < f2))
 366       {
 367          // overflow:
 368          f2 = -tools::max_value<T>() / 2;
 369       }
 370       else
 371       {
 372          f2 *= div;
 373       }
 374
 375       if(invert)
 376       {
 377          f1 = -f1;
 378          f2 = -f2;
 379       }
 380
 381       return boost::math::make_tuple(static_cast<T>(f - p), f1, f2);
 382    }
 383 private:
 384    T a, p;
 385    bool invert;
 386 };
 387
 388 template <class T, class Policy>
 389 T gamma_p_inv_imp(T a, T p, const Policy& pol)
 390 {
 391    BOOST_MATH_STD_USING  // ADL of std functions.
 392
 393    static const char* function = "boost::math::gamma_p_inv<%1%>(%1%, %1%)";
 394
 395    BOOST_MATH_INSTRUMENT_VARIABLE(a);
 396    BOOST_MATH_INSTRUMENT_VARIABLE(p);
 397
 398    if(a <= 0)
 399       policies::raise_domain_error<T>(function, "Argument a in the incomplete gamma function inverse must be >= 0 (got a=%1%).", a, pol);
 400    if((p < 0) || (p > 1))
 401       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);
 402    if(p == 1)
 403       return tools::max_value<T>();
 404    if(p == 0)
 405       return 0;
 406    bool has_10_digits;
 407    T guess = detail::find_inverse_gamma<T>(a, p, 1 - p, pol, &has_10_digits);
 408    if((policies::digits<T, Policy>() <= 36) && has_10_digits)
 409       return guess;
 410    T lower = tools::min_value<T>();
 411    if(guess <= lower)
 412       guess = tools::min_value<T>();
 413    BOOST_MATH_INSTRUMENT_VARIABLE(guess);
 414    //
 415    // Work out how many digits to converge to, normally this is
 416    // 2/3 of the digits in T, but if the first derivative is very
 417    // large convergence is slow, so we'll bump it up to full
 418    // precision to prevent premature termination of the root-finding routine.
 419    //
 420    unsigned digits = policies::digits<T, Policy>();
 421    if(digits < 30)
 422    {
 423       digits *= 2;
 424       digits /= 3;
 425    }
 426    else
 427    {
 428       digits /= 2;
 429       digits -= 1;
 430    }
 431    if((a < 0.125) && (fabs(gamma_p_derivative(a, guess, pol)) > 1 / sqrt(tools::epsilon<T>())))
 432       digits = policies::digits<T, Policy>() - 2;
 433    //
 434    // Go ahead and iterate:
 435    //
 436    boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
 437    guess = tools::halley_iterate(
 438       detail::gamma_p_inverse_func<T, Policy>(a, p, false),
 439       guess,
 440       lower,
 441       tools::max_value<T>(),
 442       digits,
 443       max_iter);
 444    policies::check_root_iterations<T>(function, max_iter, pol);
 445    BOOST_MATH_INSTRUMENT_VARIABLE(guess);
 446    if(guess == lower)
 447       guess = policies::raise_underflow_error<T>(function, "Expected result known to be non-zero, but is smaller than the smallest available number.", pol);
 448    return guess;
 449 }
 450
 451 template <class T, class Policy>
 452 T gamma_q_inv_imp(T a, T q, const Policy& pol)
 453 {
 454    BOOST_MATH_STD_USING  // ADL of std functions.
 455
 456    static const char* function = "boost::math::gamma_q_inv<%1%>(%1%, %1%)";
 457
 458    if(a <= 0)
 459       policies::raise_domain_error<T>(function, "Argument a in the incomplete gamma function inverse must be >= 0 (got a=%1%).", a, pol);
 460    if((q < 0) || (q > 1))
 461       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);
 462    if(q == 0)
 463       return tools::max_value<T>();
 464    if(q == 1)
 465       return 0;
 466    bool has_10_digits;
 467    T guess = detail::find_inverse_gamma<T>(a, 1 - q, q, pol, &has_10_digits);
 468    if((policies::digits<T, Policy>() <= 36) && has_10_digits)
 469       return guess;
 470    T lower = tools::min_value<T>();
 471    if(guess <= lower)
 472       guess = tools::min_value<T>();
 473    //
 474    // Work out how many digits to converge to, normally this is
 475    // 2/3 of the digits in T, but if the first derivative is very
 476    // large convergence is slow, so we'll bump it up to full
 477    // precision to prevent premature termination of the root-finding routine.
 478    //
 479    unsigned digits = policies::digits<T, Policy>();
 480    if(digits < 30)
 481    {
 482       digits *= 2;
 483       digits /= 3;
 484    }
 485    else
 486    {
 487       digits /= 2;
 488       digits -= 1;
 489    }
 490    if((a < 0.125) && (fabs(gamma_p_derivative(a, guess, pol)) > 1 / sqrt(tools::epsilon<T>())))
 491       digits = policies::digits<T, Policy>();
 492    //
 493    // Go ahead and iterate:
 494    //
 495    boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
 496    guess = tools::halley_iterate(
 497       detail::gamma_p_inverse_func<T, Policy>(a, q, true),
 498       guess,
 499       lower,
 500       tools::max_value<T>(),
 501       digits,
 502       max_iter);
 503    policies::check_root_iterations<T>(function, max_iter, pol);
 504    if(guess == lower)
 505       guess = policies::raise_underflow_error<T>(function, "Expected result known to be non-zero, but is smaller than the smallest available number.", pol);
 506    return guess;
 507 }
 508
 509 } // namespace detail
 510
 511 template <class T1, class T2, class Policy>
 512 inline typename tools::promote_args<T1, T2>::type
 513    gamma_p_inv(T1 a, T2 p, const Policy& pol)
 514 {
 515    typedef typename tools::promote_args<T1, T2>::type result_type;
 516    return detail::gamma_p_inv_imp(
 517       static_cast<result_type>(a),
 518       static_cast<result_type>(p), pol);
 519 }
 520
 521 template <class T1, class T2, class Policy>
 522 inline typename tools::promote_args<T1, T2>::type
 523    gamma_q_inv(T1 a, T2 p, const Policy& pol)
 524 {
 525    typedef typename tools::promote_args<T1, T2>::type result_type;
 526    return detail::gamma_q_inv_imp(
 527       static_cast<result_type>(a),
 528       static_cast<result_type>(p), pol);
 529 }
 530
 531 template <class T1, class T2>
 532 inline typename tools::promote_args<T1, T2>::type
 533    gamma_p_inv(T1 a, T2 p)
 534 {
 535    return gamma_p_inv(a, p, policies::policy<>());
 536 }
 537
 538 template <class T1, class T2>
 539 inline typename tools::promote_args<T1, T2>::type
 540    gamma_q_inv(T1 a, T2 p)
 541 {
 542    return gamma_q_inv(a, p, policies::policy<>());
 543 }
 544
 545 } // namespace math
 546 } // namespace boost
 547
 548 #endif // BOOST_MATH_SPECIAL_FUNCTIONS_IGAMMA_INVERSE_HPP
 549
 550
 551