libs/math/doc/sf/ibeta.qbk

   1 [section:ibeta_function Incomplete Beta Functions]
   2
   3 [h4 Synopsis]
   4
   5 ``
   6 #include <boost/math/special_functions/beta.hpp>
   7 ``
   8
   9    namespace boost{ namespace math{
  10
  11    template <class T1, class T2, class T3>
  12    ``__sf_result`` ibeta(T1 a, T2 b, T3 x);
  13
  14    template <class T1, class T2, class T3, class ``__Policy``>
  15    ``__sf_result`` ibeta(T1 a, T2 b, T3 x, const ``__Policy``&);
  16
  17    template <class T1, class T2, class T3>
  18    ``__sf_result`` ibetac(T1 a, T2 b, T3 x);
  19
  20    template <class T1, class T2, class T3, class ``__Policy``>
  21    ``__sf_result`` ibetac(T1 a, T2 b, T3 x, const ``__Policy``&);
  22
  23    template <class T1, class T2, class T3>
  24    ``__sf_result`` beta(T1 a, T2 b, T3 x);
  25
  26    template <class T1, class T2, class T3, class ``__Policy``>
  27    ``__sf_result`` beta(T1 a, T2 b, T3 x, const ``__Policy``&);
  28
  29    template <class T1, class T2, class T3>
  30    ``__sf_result`` betac(T1 a, T2 b, T3 x);
  31
  32    template <class T1, class T2, class T3, class ``__Policy``>
  33    ``__sf_result`` betac(T1 a, T2 b, T3 x, const ``__Policy``&);
  34
  35    }} // namespaces
  36
  37 [h4 Description]
  38
  39 There are four [@http://en.wikipedia.org/wiki/Incomplete_beta_function incomplete beta functions]
  40 : two are normalised versions (also known as ['regularized] beta functions)
  41 that return values in the range [0, 1], and two are non-normalised and
  42 return values in the range [0, __beta(a, b)].  Users interested in statistical
  43 applications should use the normalised (or
  44 [@http://mathworld.wolfram.com/RegularizedBetaFunction.html regularized]
  45 ) versions (ibeta and ibetac).
  46
  47 All of these functions require /0 <= x <= 1/.
  48
  49 The normalized functions __ibeta and __ibetac require /a,b >= 0/, and in addition that
  50 not both /a/ and /b/ are zero.
  51
  52 The functions __beta and __betac require /a,b > 0/.
  53
  54 The return type of these functions is computed using the __arg_pomotion_rules
  55 when T1, T2 and T3 are different types.
  56
  57 [optional_policy]
  58
  59    template <class T1, class T2, class T3>
  60    ``__sf_result`` ibeta(T1 a, T2 b, T3 x);
  61
  62    template <class T1, class T2, class T3, class ``__Policy``>
  63    ``__sf_result`` ibeta(T1 a, T2 b, T3 x, const ``__Policy``&);
  64
  65 Returns the normalised incomplete beta function of a, b and x:
  66
  67 [equation ibeta3]
  68
  69 [graph ibeta]
  70
  71    template <class T1, class T2, class T3>
  72    ``__sf_result`` ibetac(T1 a, T2 b, T3 x);
  73
  74    template <class T1, class T2, class T3, class ``__Policy``>
  75    ``__sf_result`` ibetac(T1 a, T2 b, T3 x, const ``__Policy``&);
  76
  77 Returns the normalised complement of the incomplete beta function of a, b and x:
  78
  79 [equation ibeta4]
  80
  81    template <class T1, class T2, class T3>
  82    ``__sf_result`` beta(T1 a, T2 b, T3 x);
  83
  84    template <class T1, class T2, class T3, class ``__Policy``>
  85    ``__sf_result`` beta(T1 a, T2 b, T3 x, const ``__Policy``&);
  86
  87 Returns the full (non-normalised) incomplete beta function of a, b and x:
  88
  89 [equation ibeta1]
  90
  91    template <class T1, class T2, class T3>
  92    ``__sf_result`` betac(T1 a, T2 b, T3 x);
  93
  94    template <class T1, class T2, class T3, class ``__Policy``>
  95    ``__sf_result`` betac(T1 a, T2 b, T3 x, const ``__Policy``&);
  96
  97 Returns the full (non-normalised) complement of the incomplete beta function of a, b and x:
  98
  99 [equation ibeta2]
 100
 101 [h4 Accuracy]
 102
 103 The following tables give peak and mean relative errors in over various domains of
 104 a, b and x, along with comparisons to the __gsl and __cephes libraries.
 105 Note that only results for the widest floating-point type on the system are given as
 106 narrower types have __zero_error.
 107
 108 Note that the results for 80 and 128-bit long doubles are noticeably higher than
 109 for doubles: this is because the wider exponent range of these types allow
 110 more extreme test cases to be tested.  For example expected results that
 111 are zero at double precision, may be finite but exceptionally small with
 112 the wider exponent range of the long double types.
 113
 114 [table Errors In the Function ibeta(a,b,x)
 115 [[Significand Size] [Platform and Compiler] [0 < a,b < 10
 116
 117 and
 118
 119 0 < x < 1]  [0 < a,b < 100
 120
 121 and
 122
 123 0 < x < 1][1x10[super -5] < a,b < 1x10[super 5]
 124
 125 and
 126
 127 0 < x < 1]]
 128 [[53] [Win32, Visual C++ 8]
 129       [Peak=42.3 Mean=2.9
 130
 131 (GSL Peak=682 Mean=32.5)
 132
 133 (__cephes Peak=42.7 Mean=7.0)]
 134       [Peak=108 Mean=16.6
 135
 136 (GSL Peak=690 Mean=151)
 137
 138 (__cephes Peak=1545 Mean=218)]
 139       [Peak=4x10[super 3][space] Mean=203
 140
 141 (GSL Peak~3x10[super 5][space] Mean~2x10[super 4][space])
 142
 143 (__cephes Peak~5x10[super 5][space] Mean~2x10[super 4][space])]]
 144 [[64] [Redhat Linux IA32, gcc-3.4.4] [Peak=21.9 Mean=3.1]
 145       [Peak=270.7 Mean=26.8] [Peak~5x10[super 4][space] Mean=3x10[super 3][space] ]]
 146 [[64] [Redhat Linux IA64, gcc-3.4.4] [Peak=15.4 Mean=3.0] [Peak=112.9 Mean=14.3]
 147       [Peak~5x10[super 4][space] Mean=3x10[super 3][space]]]
 148 [[113] [HPUX IA64, aCC A.06.06] [Peak=20.9 Mean=2.6] [Peak=88.1 Mean=14.3]
 149       [Peak~2x10[super 4][space] Mean=1x10[super 3][space] ]]
 150 ]
 151
 152 [table Errors In the Function ibetac(a,b,x)
 153 [[Significand Size] [Platform and Compiler] [0 < a,b < 10
 154
 155 and
 156
 157 0 < x < 1]  [0 < a,b < 100
 158
 159 and
 160
 161 0 < x < 1][1x10[super -5] < a,b < 1x10[super 5]
 162
 163 and
 164
 165 0 < x < 1]]
 166 [[53] [Win32, Visual C++ 8] [Peak=13.9 Mean=2.0]
 167       [Peak=56.2 Mean=14] [Peak=3x10[super 3][space] Mean=159]]
 168 [[64] [Redhat Linux IA32, gcc-3.4.4] [Peak=21.1 Mean=3.6]
 169       [Peak=221.7 Mean=25.8]
 170       [Peak~9x10[super 4][space] Mean=3x10[super 3][space] ]]
 171 [[64] [Redhat Linux IA64, gcc-3.4.4] [Peak=10.6 Mean=2.2]
 172       [Peak=73.9 Mean=11.9]
 173       [Peak~9x10[super 4][space] Mean=3x10[super 3][space] ]]
 174 [[113] [HPUX IA64, aCC A.06.06] [Peak=9.9 Mean=2.6]
 175       [Peak=117.7 Mean=15.1]
 176       [Peak~3x10[super 4][space] Mean=1x10[super 3][space] ]]
 177 ]
 178
 179 [table Errors In the Function beta(a, b, x)
 180 [[Significand Size] [Platform and Compiler] [0 < a,b < 10
 181
 182 and
 183
 184 0 < x < 1]  [0 < a,b < 100
 185
 186 and
 187
 188 0 < x < 1][1x10[super -5] < a,b < 1x10[super 5]
 189
 190 and
 191
 192 0 < x < 1]]
 193 [[53] [Win32, Visual C++ 8] [Peak=39 Mean=2.9] [Peak=91 Mean=12.7] [Peak=635 Mean=25]]
 194 [[64] [Redhat Linux IA32, gcc-3.4.4] [Peak=26 Mean=3.6] [Peak=180.7 Mean=30.1] [Peak~7x10[super 4][space] Mean=3x10[super 3][space] ]]
 195 [[64] [Redhat Linux IA64, gcc-3.4.4] [Peak=13 Mean=2.4] [Peak=67.1 Mean=13.4] [Peak~7x10[super 4][space] Mean=3x10[super 3][space] ]]
 196 [[113] [HPUX IA64, aCC A.06.06] [Peak=27.3 Mean=3.6] [Peak=49.8 Mean=9.1] [Peak~6x10[super 4][space] Mean=3x10[super 3][space] ]]
 197 ]
 198
 199 [table Errors In the Function betac(a,b,x)
 200 [[Significand Size] [Platform and Compiler] [0 < a,b < 10
 201
 202 and
 203
 204 0 < x < 1]  [0 < a,b < 100
 205
 206 and
 207
 208 0 < x < 1][1x10[super -5] < a,b < 1x10[super 5]
 209
 210 and
 211
 212 0 < x < 1]]
 213 [[53] [Win32, Visual C++ 8] [Peak=12.0 Mean=2.4] [Peak=91 Mean=15] [Peak=4x10[super 3][space] Mean=113]]
 214 [[64] [Redhat Linux IA32, gcc-3.4.4] [Peak=19.8 Mean=3.8] [Peak=295.1 Mean=33.9] [Peak~1x10[super 5][space] Mean=5x10[super 3][space] ]]
 215 [[64] [Redhat Linux IA64, gcc-3.4.4] [Peak=11.2 Mean=2.4] [Peak=63.5 Mean=13.6] [Peak~1x10[super 5][space] Mean=5x10[super 3][space]  ]]
 216 [[113] [HPUX IA64, aCC A.06.06] [Peak=15.6 Mean=3.5] [Peak=39.8 Mean=8.9] [Peak~9x10[super 4][space] Mean=5x10[super 3][space]  ]]
 217 ]
 218
 219 [h4 Testing]
 220
 221 There are two sets of tests: spot tests compare values taken from
 222 [@http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp?name=BetaRegularized Mathworld's online function evaluator]
 223 with this implementation: they provide a basic "sanity check"
 224 for the implementation, with one spot-test in each implementation-domain
 225 (see implementation notes below).
 226
 227 Accuracy tests use data generated at very high precision
 228 (with [@http://shoup.net/ntl/doc/RR.txt NTL RR class] set at 1000-bit precision),
 229 using the "textbook" continued fraction representation (refer to the first continued
 230 fraction in the implementation discussion below).
 231 Note that this continued fraction is /not/ used in the implementation,
 232 and therefore we have test data that is fully independent of the code.
 233
 234 [h4 Implementation]
 235
 236 This implementation is closely based upon
 237 [@http://portal.acm.org/citation.cfm?doid=131766.131776 "Algorithm 708; Significant digit computation of the incomplete beta function ratios", DiDonato and Morris, ACM, 1992.]
 238
 239 All four of these functions share a common implementation: this is passed both
 240 x and y, and can return either p or q where these are related by:
 241
 242 [equation ibeta_inv5]
 243
 244 so at any point we can swap a for b, x for y and p for q if this results in
 245 a more favourable position.  Generally such swaps are performed so that we always
 246 compute a value less than 0.9: when required this can then be subtracted from 1
 247 without undue cancellation error.
 248
 249 The following continued fraction representation is found in many textbooks
 250 but is not used in this implementation - it's both slower and less accurate than
 251 the alternatives - however it is used to generate test data:
 252
 253 [equation ibeta5]
 254
 255 The following continued fraction is due to [@http://portal.acm.org/citation.cfm?doid=131766.131776 Didonato and Morris],
 256 and is used in this implementation when a and b are both greater than 1:
 257
 258 [equation ibeta6]
 259
 260 For smallish b and x then a series representation can be used:
 261
 262 [equation ibeta7]
 263
 264 When b << a then the transition from 0 to 1 occurs very close to x = 1
 265 and some care has to be taken over the method of computation, in that case
 266 the following series representation is used:
 267
 268 [equation ibeta8]
 269 [/[equation ibeta9]]
 270
 271 Where Q(a,x) is an [@http://functions.wolfram.com/GammaBetaErf/Gamma2/ incomplete gamma function].
 272 Note that this method relies
 273 on keeping a table of all the p[sub n ] previously computed, which does limit
 274 the precision of the method, depending upon the size of the table used.
 275
 276 When /a/ and /b/ are both small integers, then we can relate the incomplete
 277 beta to the binomial distribution and use the following finite sum:
 278
 279 [equation ibeta12]
 280
 281 Finally we can sidestep difficult areas, or move to an area with a more
 282 efficient means of computation, by using the duplication formulae:
 283
 284 [equation ibeta10]
 285
 286 [equation ibeta11]
 287
 288 The domains of a, b and x for which the various methods are used are identical
 289 to those described in the
 290 [@http://portal.acm.org/citation.cfm?doid=131766.131776 Didonato and Morris TOMS 708 paper].
 291
 292 [endsect][/section:ibeta_function The Incomplete Beta Function]
 293
 294 [/
 295   Copyright 2006 John Maddock and Paul A. Bristow.
 296   Distributed under the Boost Software License, Version 1.0.
 297   (See accompanying file LICENSE_1_0.txt or copy at
 298   http://www.boost.org/LICENSE_1_0.txt).
 299 ]