ceph/src/boost/libs/math/doc/sf/ibeta_inv.qbk

   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
  62 There are six [@http://functions.wolfram.com/GammaBetaErf/ incomplete beta function inverses]
  63 which allow you solve for
  64 any of the three parameters to the incomplete beta, starting from either
  65 the 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
  70 beta function, they are talking about inverting on parameter /x/.
  71 These are implemented here as ibeta_inv and ibetac_inv, and are by
  72 far the most efficient of the inverses presented here.
  73
  74 The inverses on the /a/ and /b/ parameters find use in some statistical
  75 applications, but have to be computed by rather brute force numerical
  76 techniques and are consequently several times slower.
  77 These are implemented here as ibeta_inva and ibeta_invb,
  78 and complement versions ibetac_inva and ibetac_invb.]
  79
  80 The return type of these functions is computed using the __arg_promotion_rules
  81 when 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
  95 Returns a value /x/ such that: `p = ibeta(a, b, x);`
  96 and sets `*py = 1 - x` when the `py` parameter is provided and is non-null.
  97 Note 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
  99 cancellation errors.  That means that even if the function returns `1`, the
 100 value stored in `*py` may be non-zero, albeit very small.
 101
 102 Requires:  /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
 118 Returns a value /x/ such that: `q = ibetac(a, b, x);`
 119 and sets `*py = 1 - x` when the `py` parameter is provided and is non-null.
 120 Note 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
 122 cancellation errors.  That means that even if the function returns `1`, the
 123 value stored in `*py` may be non-zero, albeit very small.
 124
 125 Requires:  /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
 135 Returns a value /a/ such that: `p = ibeta(a, b, x);`
 136
 137 Requires:  /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
 147 Returns a value /a/ such that: `q = ibetac(a, b, x);`
 148
 149 Requires:  /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
 159 Returns a value /b/ such that: `p = ibeta(a, b, x);`
 160
 161 Requires:  /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
 171 Returns a value /b/ such that: `q = ibetac(a, b, x);`
 172
 173 Requires:  /a > 0/, /0 < x < 1/ and /0 <= q <= 1/.
 174
 175 [optional_policy]
 176
 177 [h4 Accuracy]
 178
 179 The accuracy of these functions should closely follow that
 180 of the regular forward incomplete beta functions.  However,
 181 note that in some parts of their domain, these functions can
 182 be extremely sensitive to changes in input, particularly when
 183 the argument /p/ (or it's complement /q/) is very close to `0` or `1`.
 184
 185 Comparisons to other libraries are shown below, note that our test data
 186 exercises some rather extreme cases in the incomplete beta function
 187 which 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
 203 There 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
 207 generous tolerances: in general both the incomplete beta and its
 208 inverses change so rapidly, that round tripping to more than a couple
 209 of 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
 212 back where you started.
 213 * Accuracy checks using high precision test values.  These measure
 214 the accuracy of the result, given exact input values.
 215
 216 [h4 Implementation of ibeta_inv and ibetac_inv]
 217
 218 These two functions share a common implementation.
 219
 220 First an initial approximation to x is computed then the
 221 last few bits are cleaned up using
 222 [@http://en.wikipedia.org/wiki/Simple_rational_approximation Halley iteration].
 223 The iteration limit is set to 1/2 of the number of bits in T, which by experiment is
 224 sufficient to ensure that the inverses are at least as accurate as the normal
 225 incomplete beta functions.  Up to 5 iterations may be
 226 required in extreme cases, although normally only one or two are required.
 227 Further, the number of iterations required decreases with increasing /a/ and
 228 /b/ (which generally form the more important use cases).
 229
 230 The initial guesses used for iteration are obtained as follows:
 231
 232 Firstly recall that:
 233
 234 [equation ibeta_inv5]
 235
 236 We may wish to start from either p or q, and to calculate either x or y.
 237 In addition at
 238 any stage we can exchange a for b, p for q, and x for y if it results in a
 239 more manageable problem.
 240
 241 For `a+b >= 5` the initial guess is computed using the methods described in:
 242
 243 Asymptotic Inversion of the Incomplete Beta Function,
 244 by N. M. [@http://homepages.cwi.nl/~nicot/ Temme].
 245 Journal of Computational and Applied Mathematics 41 (1992) 145-157.
 246
 247 The nearly symmetrical case (section 2 of the paper) is used for
 248
 249 [equation ibeta_inv2]
 250
 251 and involves solving the inverse error function first.  The method is accurate
 252 to at least 2 decimal digits when [^a = 5] rising to at least 8 digits when
 253 [^a = 10[super 5]].
 254
 255 The general error function case (section 3 of the paper) is used for
 256
 257 [equation ibeta_inv3]
 258
 259 and again expresses the inverse incomplete beta in terms of the
 260 inverse of the error function.  The method is accurate to at least
 261 2 decimal digits when [^a+b = 5] rising to 11 digits when [^a+b = 10[super 5]].
 262 However, when the result is expected to be very small, and when a+b is
 263 also small, then its accuracy tails off, in this case when p[super 1/a] < 0.0025
 264 then it is better to use the following as an initial estimate:
 265
 266 [equation ibeta_inv4]
 267
 268 Finally the for all other cases where `a+b > 5` the method of section
 269 4 of the paper is used.  This expresses the inverse incomplete beta in terms
 270 of the inverse of the incomplete gamma function, and is therefore significantly
 271 more expensive to compute than the other cases.  However the method is accurate
 272 to 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
 274 an exchange a for b, p for q and x for y when this is not the case.  In addition
 275 when p is close to 1 the method is inaccurate should we actually want y rather
 276 than x as output.  Therefore when q is small ([^q[super 1/p] < 10[super -3]]) we use:
 277
 278 [equation ibeta_inv6]
 279
 280 which is both cheaper to compute than the full method, and a more accurate
 281 estimate on q.
 282
 283 When a and b are both small there is a distinct lack of information in the
 284 literature on how to proceed.  I am extremely grateful to Prof Nico Temme
 285 who provided the following information with a great deal of patience and
 286 explanation on his part.  Any errors that follow are entirely my own, and not
 287 Prof Temme's.
 288
 289 When a and b are both less than 1, then there is a point of inflection in
 290 the 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
 292 look for a point x below the point of inflection `xs`, and on a convex curve.
 293 An initial estimate for x is made with:
 294
 295 [equation ibeta_inv7]
 296
 297 which is provably below the true value for x:
 298 [@http://en.wikipedia.org/wiki/Newton%27s_method Newton iteration] will
 299 therefore smoothly converge on x without problems caused by overshooting etc.
 300
 301 When a and b are both greater than 1, but a+b is too small to use the other
 302 methods mentioned above, we proceed as follows.  Observe that there is a point
 303 of inflection in the incomplete beta at `xs = (1 - a) / (2 - a - b)`.  Therefore
 304 if [^p > I[sub x](a,b)] we swap a for b, p for q and x for y, so that now we always
 305 look for a point x below the point of inflection `xs`, and on a concave curve.
 306 An initial estimate for x is made with:
 307
 308 [equation ibeta_inv4]
 309
 310 which can be improved somewhat to:
 311
 312 [equation ibeta_inv1]
 313
 314 when b and x are both small (I've used b < a and x < 0.2).  This
 315 actually under-estimates x, which drops us on the wrong side of x for Newton
 316 iteration to converge monotonically.  However, use of higher derivatives
 317 and Halley iteration keeps everything under control.
 318
 319 The final case to be considered if when one of a and b is less than or equal
 320 to 1, and the other greater that 1.  Here, if b < a we swap a for b, p for q
 321 and x for y.  Now the curve of the incomplete beta is convex with no points
 322 of inflection in [0,1].  For small p, x can be estimated using
 323
 324 [equation ibeta_inv4]
 325
 326 which under-estimates x, and drops us on the right side of the true value
 327 for Newton iteration to converge monotonically.  However, when p is large
 328 this can quite badly underestimate x.  This is especially an issue when we
 329 really want to find y, in which case this method can be an arbitrary number
 330 of order of magnitudes out, leading to very poor convergence during iteration.
 331
 332 Things can be improved by considering the incomplete beta as a distorted
 333 quarter circle, and estimating y from:
 334
 335 [equation ibeta_inv8]
 336
 337 This doesn't guarantee that we will drop in on the right side of x for
 338 monotonic convergence, but it does get us close enough that Halley iteration
 339 rapidly converges on the true value.
 340
 341 [h4 Implementation of inverses on the a and b parameters]
 342
 343 These four functions share a common implementation.
 344
 345 First an initial approximation is computed for /a/ or /b/:
 346 where possible this uses a Cornish-Fisher expansion for the
 347 negative binomial distribution to get within around 1 of the
 348 result.  However, when /a/ or /b/ are very small the Cornish Fisher
 349 expansion is not usable, in this case the initial approximation
 350 is chosen so that
 351 I[sub x](a, b) is near the middle of the range [0,1].
 352
 353 This initial guess is then
 354 used as a starting value for a generic root finding algorithm. The
 355 algorithm converges rapidly on the root once it has been
 356 bracketed, but bracketing the root may take several iterations.
 357 A better initial approximation for /a/ or /b/ would improve these
 358 functions quite substantially: currently 10-20 incomplete beta
 359 function 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 ]