| [section:ibeta_inv_function The Incomplete Beta Function Inverses] |
| |
| `` |
| #include <boost/math/special_functions/beta.hpp> |
| `` |
| |
| namespace boost{ namespace math{ |
| |
| template <class T1, class T2, class T3> |
| ``__sf_result`` ibeta_inv(T1 a, T2 b, T3 p); |
| |
| template <class T1, class T2, class T3, class ``__Policy``> |
| ``__sf_result`` ibeta_inv(T1 a, T2 b, T3 p, const ``__Policy``&); |
| |
| template <class T1, class T2, class T3, class T4> |
| ``__sf_result`` ibeta_inv(T1 a, T2 b, T3 p, T4* py); |
| |
| template <class T1, class T2, class T3, class T4, class ``__Policy``> |
| ``__sf_result`` ibeta_inv(T1 a, T2 b, T3 p, T4* py, const ``__Policy``&); |
| |
| template <class T1, class T2, class T3> |
| ``__sf_result`` ibetac_inv(T1 a, T2 b, T3 q); |
| |
| template <class T1, class T2, class T3, class ``__Policy``> |
| ``__sf_result`` ibetac_inv(T1 a, T2 b, T3 q, const ``__Policy``&); |
| |
| template <class T1, class T2, class T3, class T4> |
| ``__sf_result`` ibetac_inv(T1 a, T2 b, T3 q, T4* py); |
| |
| template <class T1, class T2, class T3, class T4, class ``__Policy``> |
| ``__sf_result`` ibetac_inv(T1 a, T2 b, T3 q, T4* py, const ``__Policy``&); |
| |
| template <class T1, class T2, class T3> |
| ``__sf_result`` ibeta_inva(T1 b, T2 x, T3 p); |
| |
| template <class T1, class T2, class T3, class ``__Policy``> |
| ``__sf_result`` ibeta_inva(T1 b, T2 x, T3 p, const ``__Policy``&); |
| |
| template <class T1, class T2, class T3> |
| ``__sf_result`` ibetac_inva(T1 b, T2 x, T3 q); |
| |
| template <class T1, class T2, class T3, class ``__Policy``> |
| ``__sf_result`` ibetac_inva(T1 b, T2 x, T3 q, const ``__Policy``&); |
| |
| template <class T1, class T2, class T3> |
| ``__sf_result`` ibeta_invb(T1 a, T2 x, T3 p); |
| |
| template <class T1, class T2, class T3, class ``__Policy``> |
| ``__sf_result`` ibeta_invb(T1 a, T2 x, T3 p, const ``__Policy``&); |
| |
| template <class T1, class T2, class T3> |
| ``__sf_result`` ibetac_invb(T1 a, T2 x, T3 q); |
| |
| template <class T1, class T2, class T3, class ``__Policy``> |
| ``__sf_result`` ibetac_invb(T1 a, T2 x, T3 q, const ``__Policy``&); |
| |
| }} // namespaces |
| |
| [h4 Description] |
| |
| |
| There are six [@http://functions.wolfram.com/GammaBetaErf/ incomplete beta function inverses] |
| which allow you solve for |
| any of the three parameters to the incomplete beta, starting from either |
| the result of the incomplete beta (p) or its complement (q). |
| |
| [optional_policy] |
| |
| [tip When people normally talk about the inverse of the incomplete |
| beta function, they are talking about inverting on parameter /x/. |
| These are implemented here as ibeta_inv and ibetac_inv, and are by |
| far the most efficient of the inverses presented here. |
| |
| The inverses on the /a/ and /b/ parameters find use in some statistical |
| applications, but have to be computed by rather brute force numerical |
| techniques and are consequently several times slower. |
| These are implemented here as ibeta_inva and ibeta_invb, |
| and complement versions ibetac_inva and ibetac_invb.] |
| |
| The return type of these functions is computed using the __arg_pomotion_rules |
| when called with arguments T1...TN of different types. |
| |
| template <class T1, class T2, class T3> |
| ``__sf_result`` ibeta_inv(T1 a, T2 b, T3 p); |
| |
| template <class T1, class T2, class T3, class ``__Policy``> |
| ``__sf_result`` ibeta_inv(T1 a, T2 b, T3 p, const ``__Policy``&); |
| |
| template <class T1, class T2, class T3, class T4> |
| ``__sf_result`` ibeta_inv(T1 a, T2 b, T3 p, T4* py); |
| |
| template <class T1, class T2, class T3, class T4, class ``__Policy``> |
| ``__sf_result`` ibeta_inv(T1 a, T2 b, T3 p, T4* py, const ``__Policy``&); |
| |
| Returns a value /x/ such that: `p = ibeta(a, b, x);` |
| and sets `*py = 1 - x` when the `py` parameter is provided and is non-null. |
| Note that internally this function computes whichever is the smaller of |
| `x` and `1-x`, and therefore the value assigned to `*py` is free from |
| cancellation errors. That means that even if the function returns `1`, the |
| value stored in `*py` may be non-zero, albeit very small. |
| |
| Requires: /a,b > 0/ and /0 <= p <= 1/. |
| |
| [optional_policy] |
| |
| template <class T1, class T2, class T3> |
| ``__sf_result`` ibetac_inv(T1 a, T2 b, T3 q); |
| |
| template <class T1, class T2, class T3, class ``__Policy``> |
| ``__sf_result`` ibetac_inv(T1 a, T2 b, T3 q, const ``__Policy``&); |
| |
| template <class T1, class T2, class T3, class T4> |
| ``__sf_result`` ibetac_inv(T1 a, T2 b, T3 q, T4* py); |
| |
| template <class T1, class T2, class T3, class T4, class ``__Policy``> |
| ``__sf_result`` ibetac_inv(T1 a, T2 b, T3 q, T4* py, const ``__Policy``&); |
| |
| Returns a value /x/ such that: `q = ibetac(a, b, x);` |
| and sets `*py = 1 - x` when the `py` parameter is provided and is non-null. |
| Note that internally this function computes whichever is the smaller of |
| `x` and `1-x`, and therefore the value assigned to `*py` is free from |
| cancellation errors. That means that even if the function returns `1`, the |
| value stored in `*py` may be non-zero, albeit very small. |
| |
| Requires: /a,b > 0/ and /0 <= q <= 1/. |
| |
| [optional_policy] |
| |
| template <class T1, class T2, class T3> |
| ``__sf_result`` ibeta_inva(T1 b, T2 x, T3 p); |
| |
| template <class T1, class T2, class T3, class ``__Policy``> |
| ``__sf_result`` ibeta_inva(T1 b, T2 x, T3 p, const ``__Policy``&); |
| |
| Returns a value /a/ such that: `p = ibeta(a, b, x);` |
| |
| Requires: /b > 0/, /0 < x < 1/ and /0 <= p <= 1/. |
| |
| [optional_policy] |
| |
| template <class T1, class T2, class T3> |
| ``__sf_result`` ibetac_inva(T1 b, T2 x, T3 p); |
| |
| template <class T1, class T2, class T3, class ``__Policy``> |
| ``__sf_result`` ibetac_inva(T1 b, T2 x, T3 p, const ``__Policy``&); |
| |
| Returns a value /a/ such that: `q = ibetac(a, b, x);` |
| |
| Requires: /b > 0/, /0 < x < 1/ and /0 <= q <= 1/. |
| |
| [optional_policy] |
| |
| template <class T1, class T2, class T3> |
| ``__sf_result`` ibeta_invb(T1 b, T2 x, T3 p); |
| |
| template <class T1, class T2, class T3, class ``__Policy``> |
| ``__sf_result`` ibeta_invb(T1 b, T2 x, T3 p, const ``__Policy``&); |
| |
| Returns a value /b/ such that: `p = ibeta(a, b, x);` |
| |
| Requires: /a > 0/, /0 < x < 1/ and /0 <= p <= 1/. |
| |
| [optional_policy] |
| |
| template <class T1, class T2, class T3> |
| ``__sf_result`` ibetac_invb(T1 b, T2 x, T3 p); |
| |
| template <class T1, class T2, class T3, class ``__Policy``> |
| ``__sf_result`` ibetac_invb(T1 b, T2 x, T3 p, const ``__Policy``&); |
| |
| Returns a value /b/ such that: `q = ibetac(a, b, x);` |
| |
| Requires: /a > 0/, /0 < x < 1/ and /0 <= q <= 1/. |
| |
| [optional_policy] |
| |
| [h4 Accuracy] |
| |
| The accuracy of these functions should closely follow that |
| of the regular forward incomplete beta functions. However, |
| note that in some parts of their domain, these functions can |
| be extremely sensitive to changes in input, particularly when |
| the argument /p/ (or it's complement /q/) is very close to `0` or `1`. |
| |
| [h4 Testing] |
| |
| There are two sets of tests: |
| |
| * Basic sanity checks attempt to "round-trip" from |
| /a, b/ and /x/ to /p/ or /q/ and back again. These tests have quite |
| generous tolerances: in general both the incomplete beta 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 of ibeta_inv and ibetac_inv] |
| |
| These two functions share a common implementation. |
| |
| First an initial approximation to x is computed then the |
| last few bits are cleaned up using |
| [@http://en.wikipedia.org/wiki/Simple_rational_approximation Halley iteration]. |
| The iteration limit is set to 1/2 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 beta functions. Up to 5 iterations may be |
| required in extreme cases, although normally only one or two are required. |
| Further, the number of iterations required decreases with increasing /a/ and |
| /b/ (which generally form the more important use cases). |
| |
| The initial guesses used for iteration are obtained as follows: |
| |
| Firstly recall that: |
| |
| [equation ibeta_inv5] |
| |
| We may wish to start from either p or q, and to calculate either x or y. |
| In addition at |
| any stage we can exchange a for b, p for q, and x for y if it results in a |
| more manageable problem. |
| |
| For `a+b >= 5` the initial guess is computed using the methods described in: |
| |
| Asymptotic Inversion of the Incomplete Beta Function, |
| by N. M. [@http://homepages.cwi.nl/~nicot/ Temme]. |
| Journal of Computational and Applied Mathematics 41 (1992) 145-157. |
| |
| The nearly symmetrical case (section 2 of the paper) is used for |
| |
| [equation ibeta_inv2] |
| |
| and involves solving the inverse error function first. The method is accurate |
| to at least 2 decimal digits when [^a = 5] rising to at least 8 digits when |
| [^a = 10[super 5]]. |
| |
| The general error function case (section 3 of the paper) is used for |
| |
| [equation ibeta_inv3] |
| |
| and again expresses the inverse incomplete beta in terms of the |
| inverse of the error function. The method is accurate to at least |
| 2 decimal digits when [^a+b = 5] rising to 11 digits when [^a+b = 10[super 5]]. |
| However, when the result is expected to be very small, and when a+b is |
| also small, then its accuracy tails off, in this case when p[super 1/a] < 0.0025 |
| then it is better to use the following as an initial estimate: |
| |
| [equation ibeta_inv4] |
| |
| Finally the for all other cases where `a+b > 5` the method of section |
| 4 of the paper is used. This expresses the inverse incomplete beta in terms |
| of the inverse of the incomplete gamma function, and is therefore significantly |
| more expensive to compute than the other cases. However the method is accurate |
| to at least 3 decimal digits when [^a = 5] rising to at least 10 digits when |
| [^a = 10[super 5]]. This method is limited to a > b, and therefore we need to perform |
| an exchange a for b, p for q and x for y when this is not the case. In addition |
| when p is close to 1 the method is inaccurate should we actually want y rather |
| than x as output. Therefore when q is small ([^q[super 1/p] < 10[super -3]]) we use: |
| |
| [equation ibeta_inv6] |
| |
| which is both cheaper to compute than the full method, and a more accurate |
| estimate on q. |
| |
| When a and b are both small there is a distinct lack of information in the |
| literature on how to proceed. I am extremely grateful to Prof Nico Temme |
| who provided the following information with a great deal of patience and |
| explanation on his part. Any errors that follow are entirely my own, and not |
| Prof Temme's. |
| |
| When a and b are both less than 1, then there is a point of inflection in |
| the incomplete beta at point `xs = (1 - a) / (2 - a - b)`. Therefore if |
| [^p > I[sub x](a,b)] we swap a for b, p for q and x for y, so that now we always |
| look for a point x below the point of inflection `xs`, and on a convex curve. |
| An initial estimate for x is made with: |
| |
| [equation ibeta_inv7] |
| |
| which is provably below the true value for x: |
| [@http://en.wikipedia.org/wiki/Newton%27s_method Newton iteration] will |
| therefore smoothly converge on x without problems caused by overshooting etc. |
| |
| When a and b are both greater than 1, but a+b is too small to use the other |
| methods mentioned above, we proceed as follows. Observe that there is a point |
| of inflection in the incomplete beta at `xs = (1 - a) / (2 - a - b)`. Therefore |
| if [^p > I[sub x](a,b)] we swap a for b, p for q and x for y, so that now we always |
| look for a point x below the point of inflection `xs`, and on a concave curve. |
| An initial estimate for x is made with: |
| |
| [equation ibeta_inv4] |
| |
| which can be improved somewhat to: |
| |
| [equation ibeta_inv1] |
| |
| when b and x are both small (I've used b < a and x < 0.2). This |
| actually under-estimates x, which drops us on the wrong side of x for Newton |
| iteration to converge monotonically. However, use of higher derivatives |
| and Halley iteration keeps everything under control. |
| |
| The final case to be considered if when one of a and b is less than or equal |
| to 1, and the other greater that 1. Here, if b < a we swap a for b, p for q |
| and x for y. Now the curve of the incomplete beta is convex with no points |
| of inflection in [0,1]. For small p, x can be estimated using |
| |
| [equation ibeta_inv4] |
| |
| which under-estimates x, and drops us on the right side of the true value |
| for Newton iteration to converge monotonically. However, when p is large |
| this can quite badly underestimate x. This is especially an issue when we |
| really want to find y, in which case this method can be an arbitrary number |
| of order of magnitudes out, leading to very poor convergence during iteration. |
| |
| Things can be improved by considering the incomplete beta as a distorted |
| quarter circle, and estimating y from: |
| |
| [equation ibeta_inv8] |
| |
| This doesn't guarantee that we will drop in on the right side of x for |
| monotonic convergence, but it does get us close enough that Halley iteration |
| rapidly converges on the true value. |
| |
| [h4 Implementation of inverses on the a and b parameters] |
| |
| These four functions share a common implementation. |
| |
| First an initial approximation is computed for /a/ or /b/: |
| where possible this uses a Cornish-Fisher expansion for the |
| negative binomial distribution to get within around 1 of the |
| result. However, when /a/ or /b/ are very small the Cornish Fisher |
| expansion is not usable, in this case the initial approximation |
| is chosen so that |
| I[sub x](a, b) is near the middle of the range [0,1]. |
| |
| This initial guess is then |
| used as a starting value for a generic root finding algorithm. The |
| algorithm converges rapidly on the root once it has been |
| bracketed, but bracketing the root may take several iterations. |
| A better initial approximation for /a/ or /b/ would improve these |
| functions quite substantially: currently 10-20 incomplete beta |
| function invocations are required to find the root. |
| |
| [endsect][/section:ibeta_inv_function The Incomplete Beta 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). |
| ] |