123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169 |
- [section:igamma_inv Incomplete Gamma Function Inverses]
- [h4 Synopsis]
- ``
- #include <boost/math/special_functions/gamma.hpp>
- ``
- namespace boost{ namespace math{
-
- template <class T1, class T2>
- ``__sf_result`` gamma_q_inv(T1 a, T2 q);
-
- template <class T1, class T2, class ``__Policy``>
- ``__sf_result`` gamma_q_inv(T1 a, T2 q, const ``__Policy``&);
-
- template <class T1, class T2>
- ``__sf_result`` gamma_p_inv(T1 a, T2 p);
-
- template <class T1, class T2, class ``__Policy``>
- ``__sf_result`` gamma_p_inv(T1 a, T2 p, const ``__Policy``&);
-
- template <class T1, class T2>
- ``__sf_result`` gamma_q_inva(T1 x, T2 q);
-
- template <class T1, class T2, class ``__Policy``>
- ``__sf_result`` gamma_q_inva(T1 x, T2 q, const ``__Policy``&);
-
- template <class T1, class T2>
- ``__sf_result`` gamma_p_inva(T1 x, T2 p);
-
- template <class T1, class T2, class ``__Policy``>
- ``__sf_result`` gamma_p_inva(T1 x, T2 p, const ``__Policy``&);
-
- }} // namespaces
-
- [h4 Description]
- There are four [@http://mathworld.wolfram.com/IncompleteGammaFunction.html incomplete gamma function]
- inverses which either compute
- /x/ given /a/ and /p/ or /q/,
- or else compute /a/ given /x/ and either /p/ or /q/.
- The return type of these functions is computed using the __arg_promotion_rules
- when T1 and T2 are different types, otherwise the return type is simply T1.
- [optional_policy]
- [tip When people normally talk about the inverse of the incomplete
- gamma function, they are talking about inverting on parameter /x/.
- These are implemented here as `gamma_p_inv` and `gamma_q_inv`, and are by
- far the most efficient of the inverses presented here.
- The inverse on the /a/ parameter finds use in some statistical
- applications but has to be computed by rather brute force numerical
- techniques and is consequently several times slower.
- These are implemented here as `gamma_p_inva` and `gamma_q_inva`.]
- template <class T1, class T2>
- ``__sf_result`` gamma_q_inv(T1 a, T2 q);
- template <class T1, class T2, class ``__Policy``>
- ``__sf_result`` gamma_q_inv(T1 a, T2 q, const ``__Policy``&);
- Returns a value x such that: `q = gamma_q(a, x);`
- Requires: /a > 0/ and /1 >= p,q >= 0/.
- template <class T1, class T2>
- ``__sf_result`` gamma_p_inv(T1 a, T2 p);
-
- template <class T1, class T2, class ``__Policy``>
- ``__sf_result`` gamma_p_inv(T1 a, T2 p, const ``__Policy``&);
-
- Returns a value x such that: `p = gamma_p(a, x);`
- Requires: /a > 0/ and /1 >= p,q >= 0/.
- template <class T1, class T2>
- ``__sf_result`` gamma_q_inva(T1 x, T2 q);
- template <class T1, class T2, class ``__Policy``>
- ``__sf_result`` gamma_q_inva(T1 x, T2 q, const ``__Policy``&);
- Returns a value a such that: `q = gamma_q(a, x);`
- Requires: /x > 0/ and /1 >= p,q >= 0/.
- template <class T1, class T2>
- ``__sf_result`` gamma_p_inva(T1 x, T2 p);
-
- template <class T1, class T2, class ``__Policy``>
- ``__sf_result`` gamma_p_inva(T1 x, T2 p, const ``__Policy``&);
-
- Returns a value a such that: `p = gamma_p(a, x);`
- Requires: /x > 0/ and /1 >= p,q >= 0/.
- [h4 Accuracy]
- The accuracy of these functions doesn't vary much by platform or by
- the type T. Given that these functions are computed by iterative methods,
- they are deliberately "detuned" so as not to be too accurate: it is in
- any case impossible for these function to be more accurate than the
- regular forward incomplete gamma functions. In practice, the accuracy
- of these functions is very similar to that of __gamma_p and __gamma_q
- functions:
- [table_gamma_p_inv]
- [table_gamma_q_inv]
- [table_gamma_p_inva]
- [table_gamma_q_inva]
- [h4 Testing]
- There are two sets of tests:
- * Basic sanity checks attempt to "round-trip" from
- /a/ and /x/ to /p/ or /q/ and back again. These tests have quite
- generous tolerances: in general both the incomplete gamma, and its
- inverses, change so rapidly that round tripping to more than a couple
- of significant digits isn't possible. This is especially true when
- /p/ or /q/ is very near one: in this case there isn't enough
- "information content" in the input to the inverse function to get
- back where you started.
- * Accuracy checks using high precision test values. These measure
- the accuracy of the result, given exact input values.
- [h4 Implementation]
- The functions `gamma_p_inv` and [@http://functions.wolfram.com/GammaBetaErf/InverseGammaRegularized/ `gamma_q_inv`]
- share a common implementation.
- First an initial approximation is computed using the methodology described
- in:
- [@http://portal.acm.org/citation.cfm?id=23109&coll=portal&dl=ACM
- A. R. Didonato and A. H. Morris, Computation of the Incomplete Gamma
- Function Ratios and their Inverse, ACM Trans. Math. Software 12 (1986), 377-393.]
- Finally, the last few bits are cleaned up using Halley iteration, the iteration
- limit is set to 2/3 of the number of bits in T, which by experiment is
- sufficient to ensure that the inverses are at least as accurate as the normal
- incomplete gamma functions. In testing, no more than 3 iterations are required
- to produce a result as accurate as the forward incomplete gamma function, and
- in many cases only one iteration is required.
- The functions `gamma_p_inva` and `gamma_q_inva` also share a common implementation
- but are handled separately from `gamma_p_inv` and `gamma_q_inv`.
- An initial approximation for /a/ is computed very crudely so that
- /gamma_p(a, x) ~ 0.5/, this value is then used as a starting point
- for a generic derivative-free root finding algorithm. As a consequence,
- these two functions are rather more expensive to compute than the
- `gamma_p_inv` or `gamma_q_inv` functions. Even so, the root is usually found
- in fewer than 10 iterations.
- [endsect] [/section The Incomplete Gamma Function Inverses]
- [/
- Copyright 2006 John Maddock and Paul A. Bristow.
- Distributed under 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).
- ]
|