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 39There are four [@http://mathworld.wolfram.com/IncompleteGammaFunction.html incomplete gamma function] 40inverses which either compute 41/x/ given /a/ and /p/ or /q/, 42or else compute /a/ given /x/ and either /p/ or /q/. 43 44The return type of these functions is computed using the __arg_promotion_rules 45when 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 50gamma function, they are talking about inverting on parameter /x/. 51These are implemented here as `gamma_p_inv` and `gamma_q_inv`, and are by 52far the most efficient of the inverses presented here. 53 54The inverse on the /a/ parameter finds use in some statistical 55applications but has to be computed by rather brute force numerical 56techniques and is consequently several times slower. 57These 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 66Returns a value x such that: `q = gamma_q(a, x);` 67 68Requires: /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 76Returns a value x such that: `p = gamma_p(a, x);` 77 78Requires: /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 86Returns a value a such that: `q = gamma_q(a, x);` 87 88Requires: /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 96Returns a value a such that: `p = gamma_p(a, x);` 97 98Requires: /x > 0/ and /1 >= p,q >= 0/. 99 100[h4 Accuracy] 101 102The accuracy of these functions doesn't vary much by platform or by 103the type T. Given that these functions are computed by iterative methods, 104they are deliberately "detuned" so as not to be too accurate: it is in 105any case impossible for these function to be more accurate than the 106regular forward incomplete gamma functions. In practice, the accuracy 107of these functions is very similar to that of __gamma_p and __gamma_q 108functions: 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 120There 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 124generous tolerances: in general both the incomplete gamma, and its 125inverses, change so rapidly that round tripping to more than a couple 126of 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 129back where you started. 130* Accuracy checks using high precision test values. These measure 131the accuracy of the result, given exact input values. 132 133[h4 Implementation] 134 135The functions `gamma_p_inv` and [@http://functions.wolfram.com/GammaBetaErf/InverseGammaRegularized/ `gamma_q_inv`] 136share a common implementation. 137 138First an initial approximation is computed using the methodology described 139in: 140 141[@http://portal.acm.org/citation.cfm?id=23109&coll=portal&dl=ACM 142A. R. Didonato and A. H. Morris, Computation of the Incomplete Gamma 143Function Ratios and their Inverse, ACM Trans. Math. Software 12 (1986), 377-393.] 144 145Finally, the last few bits are cleaned up using Halley iteration, the iteration 146limit is set to 2/3 of the number of bits in T, which by experiment is 147sufficient to ensure that the inverses are at least as accurate as the normal 148incomplete gamma functions. In testing, no more than 3 iterations are required 149to produce a result as accurate as the forward incomplete gamma function, and 150in many cases only one iteration is required. 151 152The functions `gamma_p_inva` and `gamma_q_inva` also share a common implementation 153but are handled separately from `gamma_p_inv` and `gamma_q_inv`. 154 155An 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 157for a generic derivative-free root finding algorithm. As a consequence, 158these 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 160in 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] 170