libs/math/doc/sf/igamma_inv.qbk

   1 [section:igamma_inv Incomplete Gamma Function Inverses]
   2
   3 [h4 Synopsis]
   4
   5 ``
   6 #include <boost/math/special_functions/gamma.hpp>
   7 ``
   8
   9    namespace boost{ namespace math{
  10
  11    template <class T1, class T2>
  12    ``__sf_result`` gamma_q_inv(T1 a, T2 q);
  13
  14    template <class T1, class T2, class ``__Policy``>
  15    ``__sf_result`` gamma_q_inv(T1 a, T2 q, const ``__Policy``&);
  16
  17    template <class T1, class T2>
  18    ``__sf_result`` gamma_p_inv(T1 a, T2 p);
  19
  20    template <class T1, class T2, class ``__Policy``>
  21    ``__sf_result`` gamma_p_inv(T1 a, T2 p, const ``__Policy``&);
  22
  23    template <class T1, class T2>
  24    ``__sf_result`` gamma_q_inva(T1 x, T2 q);
  25
  26    template <class T1, class T2, class ``__Policy``>
  27    ``__sf_result`` gamma_q_inva(T1 x, T2 q, const ``__Policy``&);
  28
  29    template <class T1, class T2>
  30    ``__sf_result`` gamma_p_inva(T1 x, T2 p);
  31
  32    template <class T1, class T2, class ``__Policy``>
  33    ``__sf_result`` gamma_p_inva(T1 x, T2 p, const ``__Policy``&);
  34
  35    }} // namespaces
  36
  37 [h4 Description]
  38
  39 There are four [@http://mathworld.wolfram.com/IncompleteGammaFunction.html incomplete gamma function]
  40 inverses which either compute
  41 /x/ given /a/ and /p/ or /q/,
  42 or else compute /a/ given /x/ and either /p/ or /q/.
  43
  44 The return type of these functions is computed using the __arg_promotion_rules
  45 when T1 and T2 are different types, otherwise the return type is simply T1.
  46
  47 [optional_policy]
  48
  49 [tip When people normally talk about the inverse of the incomplete
  50 gamma function, they are talking about inverting on parameter /x/.
  51 These are implemented here as `gamma_p_inv` and `gamma_q_inv`, and are by
  52 far the most efficient of the inverses presented here.
  53
  54 The inverse on the /a/ parameter finds use in some statistical
  55 applications but has to be computed by rather brute force numerical
  56 techniques and is consequently several times slower.
  57 These are implemented here as `gamma_p_inva` and `gamma_q_inva`.]
  58
  59
  60    template <class T1, class T2>
  61    ``__sf_result`` gamma_q_inv(T1 a, T2 q);
  62
  63    template <class T1, class T2, class ``__Policy``>
  64    ``__sf_result`` gamma_q_inv(T1 a, T2 q, const ``__Policy``&);
  65
  66 Returns a value x such that: `q = gamma_q(a, x);`
  67
  68 Requires: /a > 0/ and /1 >= p,q >= 0/.
  69
  70    template <class T1, class T2>
  71    ``__sf_result`` gamma_p_inv(T1 a, T2 p);
  72
  73    template <class T1, class T2, class ``__Policy``>
  74    ``__sf_result`` gamma_p_inv(T1 a, T2 p, const ``__Policy``&);
  75
  76 Returns a value x such that: `p = gamma_p(a, x);`
  77
  78 Requires: /a > 0/ and /1 >= p,q >= 0/.
  79
  80    template <class T1, class T2>
  81    ``__sf_result`` gamma_q_inva(T1 x, T2 q);
  82
  83    template <class T1, class T2, class ``__Policy``>
  84    ``__sf_result`` gamma_q_inva(T1 x, T2 q, const ``__Policy``&);
  85
  86 Returns a value a such that: `q = gamma_q(a, x);`
  87
  88 Requires: /x > 0/ and /1 >= p,q >= 0/.
  89
  90    template <class T1, class T2>
  91    ``__sf_result`` gamma_p_inva(T1 x, T2 p);
  92
  93    template <class T1, class T2, class ``__Policy``>
  94    ``__sf_result`` gamma_p_inva(T1 x, T2 p, const ``__Policy``&);
  95
  96 Returns a value a such that: `p = gamma_p(a, x);`
  97
  98 Requires: /x > 0/ and /1 >= p,q >= 0/.
  99
 100 [h4 Accuracy]
 101
 102 The accuracy of these functions doesn't vary much by platform or by
 103 the type T.  Given that these functions are computed by iterative methods,
 104 they are deliberately "detuned" so as not to be too accurate: it is in
 105 any case impossible for these function to be more accurate than the
 106 regular forward incomplete gamma functions.  In practice, the accuracy
 107 of these functions is very similar to that of __gamma_p and __gamma_q
 108 functions:
 109
 110 [table_gamma_p_inv]
 111
 112 [table_gamma_q_inv]
 113
 114 [table_gamma_p_inva]
 115
 116 [table_gamma_q_inva]
 117
 118 [h4 Testing]
 119
 120 There are two sets of tests:
 121
 122 * Basic sanity checks attempt to "round-trip" from
 123 /a/ and /x/ to /p/ or /q/ and back again.  These tests have quite
 124 generous tolerances: in general both the incomplete gamma, and its
 125 inverses, change so rapidly that round tripping to more than a couple
 126 of significant digits isn't possible.  This is especially true when
 127 /p/ or /q/ is very near one: in this case there isn't enough
 128 "information content" in the input to the inverse function to get
 129 back where you started.
 130 * Accuracy checks using high precision test values.  These measure
 131 the accuracy of the result, given exact input values.
 132
 133 [h4 Implementation]
 134
 135 The functions `gamma_p_inv` and [@http://functions.wolfram.com/GammaBetaErf/InverseGammaRegularized/ `gamma_q_inv`]
 136 share a common implementation.
 137
 138 First an initial approximation is computed using the methodology described
 139 in:
 140
 141 [@http://portal.acm.org/citation.cfm?id=23109&coll=portal&dl=ACM
 142 A. R. Didonato and A. H. Morris, Computation of the Incomplete Gamma
 143 Function Ratios and their Inverse, ACM Trans. Math. Software 12 (1986), 377-393.]
 144
 145 Finally, the last few bits are cleaned up using Halley iteration, the iteration
 146 limit is set to 2/3 of the number of bits in T, which by experiment is
 147 sufficient to ensure that the inverses are at least as accurate as the normal
 148 incomplete gamma functions.  In testing, no more than 3 iterations are required
 149 to produce a result as accurate as the forward incomplete gamma function, and
 150 in many cases only one iteration is required.
 151
 152 The functions `gamma_p_inva` and `gamma_q_inva` also share a common implementation
 153 but are handled separately from `gamma_p_inv` and `gamma_q_inv`.
 154
 155 An initial approximation for /a/ is computed very crudely so that
 156 /gamma_p(a, x) ~ 0.5/, this value is then used as a starting point
 157 for a generic derivative-free root finding algorithm.  As a consequence,
 158 these two functions are rather more expensive to compute than the
 159 `gamma_p_inv` or `gamma_q_inv` functions.  Even so, the root is usually found
 160 in fewer than 10 iterations.
 161
 162 [endsect] [/section The Incomplete Gamma Function Inverses]
 163
 164 [/
 165   Copyright 2006 John Maddock and Paul A. Bristow.
 166   Distributed under the Boost Software License, Version 1.0.
 167   (See accompanying file LICENSE_1_0.txt or copy at
 168   http://www.boost.org/LICENSE_1_0.txt).
 169 ]