1[section:ibeta_inv_function The Incomplete Beta Function Inverses] 2 3`` 4#include <boost/math/special_functions/beta.hpp> 5`` 6 7 namespace boost{ namespace math{ 8 9 template <class T1, class T2, class T3> 10 ``__sf_result`` ibeta_inv(T1 a, T2 b, T3 p); 11 12 template <class T1, class T2, class T3, class ``__Policy``> 13 ``__sf_result`` ibeta_inv(T1 a, T2 b, T3 p, const ``__Policy``&); 14 15 template <class T1, class T2, class T3, class T4> 16 ``__sf_result`` ibeta_inv(T1 a, T2 b, T3 p, T4* py); 17 18 template <class T1, class T2, class T3, class T4, class ``__Policy``> 19 ``__sf_result`` ibeta_inv(T1 a, T2 b, T3 p, T4* py, const ``__Policy``&); 20 21 template <class T1, class T2, class T3> 22 ``__sf_result`` ibetac_inv(T1 a, T2 b, T3 q); 23 24 template <class T1, class T2, class T3, class ``__Policy``> 25 ``__sf_result`` ibetac_inv(T1 a, T2 b, T3 q, const ``__Policy``&); 26 27 template <class T1, class T2, class T3, class T4> 28 ``__sf_result`` ibetac_inv(T1 a, T2 b, T3 q, T4* py); 29 30 template <class T1, class T2, class T3, class T4, class ``__Policy``> 31 ``__sf_result`` ibetac_inv(T1 a, T2 b, T3 q, T4* py, const ``__Policy``&); 32 33 template <class T1, class T2, class T3> 34 ``__sf_result`` ibeta_inva(T1 b, T2 x, T3 p); 35 36 template <class T1, class T2, class T3, class ``__Policy``> 37 ``__sf_result`` ibeta_inva(T1 b, T2 x, T3 p, const ``__Policy``&); 38 39 template <class T1, class T2, class T3> 40 ``__sf_result`` ibetac_inva(T1 b, T2 x, T3 q); 41 42 template <class T1, class T2, class T3, class ``__Policy``> 43 ``__sf_result`` ibetac_inva(T1 b, T2 x, T3 q, const ``__Policy``&); 44 45 template <class T1, class T2, class T3> 46 ``__sf_result`` ibeta_invb(T1 a, T2 x, T3 p); 47 48 template <class T1, class T2, class T3, class ``__Policy``> 49 ``__sf_result`` ibeta_invb(T1 a, T2 x, T3 p, const ``__Policy``&); 50 51 template <class T1, class T2, class T3> 52 ``__sf_result`` ibetac_invb(T1 a, T2 x, T3 q); 53 54 template <class T1, class T2, class T3, class ``__Policy``> 55 ``__sf_result`` ibetac_invb(T1 a, T2 x, T3 q, const ``__Policy``&); 56 57 }} // namespaces 58 59[h4 Description] 60 61 62There are six [@http://functions.wolfram.com/GammaBetaErf/ incomplete beta function inverses] 63which allow you solve for 64any of the three parameters to the incomplete beta, starting from either 65the result of the incomplete beta (p) or its complement (q). 66 67[optional_policy] 68 69[tip When people normally talk about the inverse of the incomplete 70beta function, they are talking about inverting on parameter /x/. 71These are implemented here as `ibeta_inv` and `ibetac_inv`, and are by 72far the most efficient of the inverses presented here. 73 74The inverses on the /a/ and /b/ parameters find use in some statistical 75applications, but have to be computed by rather brute force numerical 76techniques and are consequently several times slower. 77These are implemented here as `ibeta_inva` and `ibeta_invb`, 78and complement versions `ibetac_inva` and `ibetac_invb`.] 79 80The return type of these functions is computed using the __arg_promotion_rules 81when called with arguments T1...TN of different types. 82 83 template <class T1, class T2, class T3> 84 ``__sf_result`` ibeta_inv(T1 a, T2 b, T3 p); 85 86 template <class T1, class T2, class T3, class ``__Policy``> 87 ``__sf_result`` ibeta_inv(T1 a, T2 b, T3 p, const ``__Policy``&); 88 89 template <class T1, class T2, class T3, class T4> 90 ``__sf_result`` ibeta_inv(T1 a, T2 b, T3 p, T4* py); 91 92 template <class T1, class T2, class T3, class T4, class ``__Policy``> 93 ``__sf_result`` ibeta_inv(T1 a, T2 b, T3 p, T4* py, const ``__Policy``&); 94 95Returns a value /x/ such that: `p = ibeta(a, b, x);` 96and sets `*py = 1 - x` when the `py` parameter is provided and is non-null. 97Note that internally this function computes whichever is the smaller of 98`x` and `1-x`, and therefore the value assigned to `*py` is free from 99cancellation errors. That means that even if the function returns `1`, the 100value stored in `*py` may be non-zero, albeit very small. 101 102Requires: /a,b > 0/ and /0 <= p <= 1/. 103 104[optional_policy] 105 106 template <class T1, class T2, class T3> 107 ``__sf_result`` ibetac_inv(T1 a, T2 b, T3 q); 108 109 template <class T1, class T2, class T3, class ``__Policy``> 110 ``__sf_result`` ibetac_inv(T1 a, T2 b, T3 q, const ``__Policy``&); 111 112 template <class T1, class T2, class T3, class T4> 113 ``__sf_result`` ibetac_inv(T1 a, T2 b, T3 q, T4* py); 114 115 template <class T1, class T2, class T3, class T4, class ``__Policy``> 116 ``__sf_result`` ibetac_inv(T1 a, T2 b, T3 q, T4* py, const ``__Policy``&); 117 118Returns a value /x/ such that: `q = ibetac(a, b, x);` 119and sets `*py = 1 - x` when the `py` parameter is provided and is non-null. 120Note that internally this function computes whichever is the smaller of 121`x` and `1-x`, and therefore the value assigned to `*py` is free from 122cancellation errors. That means that even if the function returns `1`, the 123value stored in `*py` may be non-zero, albeit very small. 124 125Requires: /a,b > 0/ and /0 <= q <= 1/. 126 127[optional_policy] 128 129 template <class T1, class T2, class T3> 130 ``__sf_result`` ibeta_inva(T1 b, T2 x, T3 p); 131 132 template <class T1, class T2, class T3, class ``__Policy``> 133 ``__sf_result`` ibeta_inva(T1 b, T2 x, T3 p, const ``__Policy``&); 134 135Returns a value /a/ such that: `p = ibeta(a, b, x);` 136 137Requires: /b > 0/, /0 < x < 1/ and /0 <= p <= 1/. 138 139[optional_policy] 140 141 template <class T1, class T2, class T3> 142 ``__sf_result`` ibetac_inva(T1 b, T2 x, T3 p); 143 144 template <class T1, class T2, class T3, class ``__Policy``> 145 ``__sf_result`` ibetac_inva(T1 b, T2 x, T3 p, const ``__Policy``&); 146 147Returns a value /a/ such that: `q = ibetac(a, b, x);` 148 149Requires: /b > 0/, /0 < x < 1/ and /0 <= q <= 1/. 150 151[optional_policy] 152 153 template <class T1, class T2, class T3> 154 ``__sf_result`` ibeta_invb(T1 b, T2 x, T3 p); 155 156 template <class T1, class T2, class T3, class ``__Policy``> 157 ``__sf_result`` ibeta_invb(T1 b, T2 x, T3 p, const ``__Policy``&); 158 159Returns a value /b/ such that: `p = ibeta(a, b, x);` 160 161Requires: /a > 0/, /0 < x < 1/ and /0 <= p <= 1/. 162 163[optional_policy] 164 165 template <class T1, class T2, class T3> 166 ``__sf_result`` ibetac_invb(T1 b, T2 x, T3 p); 167 168 template <class T1, class T2, class T3, class ``__Policy``> 169 ``__sf_result`` ibetac_invb(T1 b, T2 x, T3 p, const ``__Policy``&); 170 171Returns a value /b/ such that: `q = ibetac(a, b, x);` 172 173Requires: /a > 0/, /0 < x < 1/ and /0 <= q <= 1/. 174 175[optional_policy] 176 177[h4 Accuracy] 178 179The accuracy of these functions should closely follow that 180of the regular forward incomplete beta functions. However, 181note that in some parts of their domain, these functions can 182be extremely sensitive to changes in input, particularly when 183the argument /p/ (or it's complement /q/) is very close to `0` or `1`. 184 185Comparisons to other libraries are shown below, note that our test data 186exercises some rather extreme cases in the incomplete beta function 187which many other libraries fail to handle: 188 189[table_ibeta_inv] 190 191[table_ibetac_inv] 192 193[table_ibeta_inva] 194 195[table_ibetac_inva] 196 197[table_ibeta_invb] 198 199[table_ibetac_invb] 200 201[h4 Testing] 202 203There are two sets of tests: 204 205* Basic sanity checks attempt to "round-trip" from 206/a, b/ and /x/ to /p/ or /q/ and back again. These tests have quite 207generous tolerances: in general both the incomplete beta and its 208inverses change so rapidly, that round tripping to more than a couple 209of significant digits isn't possible. This is especially true when 210/p/ or /q/ is very near one: in this case there isn't enough 211"information content" in the input to the inverse function to get 212back where you started. 213* Accuracy checks using high precision test values. These measure 214the accuracy of the result, given exact input values. 215 216[h4 Implementation of ibeta_inv and ibetac_inv] 217 218These two functions share a common implementation. 219 220First an initial approximation to x is computed then the 221last few bits are cleaned up using 222[@http://en.wikipedia.org/wiki/Simple_rational_approximation Halley iteration]. 223The iteration limit is set to 1/2 of the number of bits in T, which by experiment is 224sufficient to ensure that the inverses are at least as accurate as the normal 225incomplete beta functions. Up to 5 iterations may be 226required in extreme cases, although normally only one or two are required. 227Further, the number of iterations required decreases with increasing /a/ and 228/b/ (which generally form the more important use cases). 229 230The initial guesses used for iteration are obtained as follows: 231 232Firstly recall that: 233 234[equation ibeta_inv5] 235 236We may wish to start from either p or q, and to calculate either x or y. 237In addition at 238any stage we can exchange a for b, p for q, and x for y if it results in a 239more manageable problem. 240 241For `a+b >= 5` the initial guess is computed using the methods described in: 242 243Asymptotic Inversion of the Incomplete Beta Function, 244by N. M. [@http://homepages.cwi.nl/~nicot/ Temme]. 245Journal of Computational and Applied Mathematics 41 (1992) 145-157. 246 247The nearly symmetrical case (section 2 of the paper) is used for 248 249[equation ibeta_inv2] 250 251and involves solving the inverse error function first. The method is accurate 252to at least 2 decimal digits when [^a = 5] rising to at least 8 digits when 253[^a = 10[super 5]]. 254 255The general error function case (section 3 of the paper) is used for 256 257[equation ibeta_inv3] 258 259and again expresses the inverse incomplete beta in terms of the 260inverse of the error function. The method is accurate to at least 2612 decimal digits when [^a+b = 5] rising to 11 digits when [^a+b = 10[super 5]]. 262However, when the result is expected to be very small, and when a+b is 263also small, then its accuracy tails off, in this case when p[super 1/a] < 0.0025 264then it is better to use the following as an initial estimate: 265 266[equation ibeta_inv4] 267 268Finally the for all other cases where `a+b > 5` the method of section 2694 of the paper is used. This expresses the inverse incomplete beta in terms 270of the inverse of the incomplete gamma function, and is therefore significantly 271more expensive to compute than the other cases. However the method is accurate 272to at least 3 decimal digits when [^a = 5] rising to at least 10 digits when 273[^a = 10[super 5]]. This method is limited to a > b, and therefore we need to perform 274an exchange a for b, p for q and x for y when this is not the case. In addition 275when p is close to 1 the method is inaccurate should we actually want y rather 276than x as output. Therefore when q is small ([^q[super 1/p] < 10[super -3]]) we use: 277 278[equation ibeta_inv6] 279 280which is both cheaper to compute than the full method, and a more accurate 281estimate on q. 282 283When a and b are both small there is a distinct lack of information in the 284literature on how to proceed. I am extremely grateful to Prof Nico Temme 285who provided the following information with a great deal of patience and 286explanation on his part. Any errors that follow are entirely my own, and not 287Prof Temme's. 288 289When a and b are both less than 1, then there is a point of inflection in 290the incomplete beta at point `xs = (1 - a) / (2 - a - b)`. Therefore if 291[^p > I[sub x](a,b)] we swap a for b, p for q and x for y, so that now we always 292look for a point x below the point of inflection `xs`, and on a convex curve. 293An initial estimate for x is made with: 294 295[equation ibeta_inv7] 296 297which is provably below the true value for x: 298[@http://en.wikipedia.org/wiki/Newton%27s_method Newton iteration] will 299therefore smoothly converge on x without problems caused by overshooting etc. 300 301When a and b are both greater than 1, but a+b is too small to use the other 302methods mentioned above, we proceed as follows. Observe that there is a point 303of inflection in the incomplete beta at `xs = (1 - a) / (2 - a - b)`. Therefore 304if [^p > I[sub x](a,b)] we swap a for b, p for q and x for y, so that now we always 305look for a point x below the point of inflection `xs`, and on a concave curve. 306An initial estimate for x is made with: 307 308[equation ibeta_inv4] 309 310which can be improved somewhat to: 311 312[equation ibeta_inv1] 313 314when b and x are both small (I've used b < a and x < 0.2). This 315actually under-estimates x, which drops us on the wrong side of x for Newton 316iteration to converge monotonically. However, use of higher derivatives 317and Halley iteration keeps everything under control. 318 319The final case to be considered if when one of a and b is less than or equal 320to 1, and the other greater that 1. Here, if b < a we swap a for b, p for q 321and x for y. Now the curve of the incomplete beta is convex with no points 322of inflection in [0,1]. For small p, x can be estimated using 323 324[equation ibeta_inv4] 325 326which under-estimates x, and drops us on the right side of the true value 327for Newton iteration to converge monotonically. However, when p is large 328this can quite badly underestimate x. This is especially an issue when we 329really want to find y, in which case this method can be an arbitrary number 330of order of magnitudes out, leading to very poor convergence during iteration. 331 332Things can be improved by considering the incomplete beta as a distorted 333quarter circle, and estimating y from: 334 335[equation ibeta_inv8] 336 337This doesn't guarantee that we will drop in on the right side of x for 338monotonic convergence, but it does get us close enough that Halley iteration 339rapidly converges on the true value. 340 341[h4 Implementation of inverses on the a and b parameters] 342 343These four functions share a common implementation. 344 345First an initial approximation is computed for /a/ or /b/: 346where possible this uses a Cornish-Fisher expansion for the 347negative binomial distribution to get within around 1 of the 348result. However, when /a/ or /b/ are very small the Cornish Fisher 349expansion is not usable, in this case the initial approximation 350is chosen so that 351I[sub x](a, b) is near the middle of the range [0,1]. 352 353This initial guess is then 354used as a starting value for a generic root finding algorithm. The 355algorithm converges rapidly on the root once it has been 356bracketed, but bracketing the root may take several iterations. 357A better initial approximation for /a/ or /b/ would improve these 358functions quite substantially: currently 10-20 incomplete beta 359function invocations are required to find the root. 360 361[endsect][/section:ibeta_inv_function The Incomplete Beta Function Inverses] 362 363[/ 364 Copyright 2006 John Maddock and Paul A. Bristow. 365 Distributed under the Boost Software License, Version 1.0. 366 (See accompanying file LICENSE_1_0.txt or copy at 367 http://www.boost.org/LICENSE_1_0.txt). 368] 369