ceph/src/boost/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_promotion_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_ibeta]
 115
 116 [table_ibetac]
 117
 118 [table_beta_incomplete_]
 119
 120 [table_betac]
 121
 122 [h4 Testing]
 123
 124 There are two sets of tests: spot tests compare values taken from
 125 [@http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp?name=BetaRegularized Mathworld's online function evaluator]
 126 with this implementation: they provide a basic "sanity check"
 127 for the implementation, with one spot-test in each implementation-domain
 128 (see implementation notes below).
 129
 130 Accuracy tests use data generated at very high precision
 131 (with [@http://shoup.net/ntl/doc/RR.txt NTL RR class] set at 1000-bit precision),
 132 using the "textbook" continued fraction representation (refer to the first continued
 133 fraction in the implementation discussion below).
 134 Note that this continued fraction is /not/ used in the implementation,
 135 and therefore we have test data that is fully independent of the code.
 136
 137 [h4 Implementation]
 138
 139 This implementation is closely based upon
 140 [@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.]
 141
 142 All four of these functions share a common implementation: this is passed both
 143 x and y, and can return either p or q where these are related by:
 144
 145 [equation ibeta_inv5]
 146
 147 so at any point we can swap a for b, x for y and p for q if this results in
 148 a more favourable position.  Generally such swaps are performed so that we always
 149 compute a value less than 0.9: when required this can then be subtracted from 1
 150 without undue cancellation error.
 151
 152 The following continued fraction representation is found in many textbooks
 153 but is not used in this implementation - it's both slower and less accurate than
 154 the alternatives - however it is used to generate test data:
 155
 156 [equation ibeta5]
 157
 158 The following continued fraction is due to [@http://portal.acm.org/citation.cfm?doid=131766.131776 Didonato and Morris],
 159 and is used in this implementation when a and b are both greater than 1:
 160
 161 [equation ibeta6]
 162
 163 For smallish b and x then a series representation can be used:
 164
 165 [equation ibeta7]
 166
 167 When b << a then the transition from 0 to 1 occurs very close to x = 1
 168 and some care has to be taken over the method of computation, in that case
 169 the following series representation is used:
 170
 171 [equation ibeta8]
 172 [/[equation ibeta9]]
 173
 174 Where Q(a,x) is an [@http://functions.wolfram.com/GammaBetaErf/Gamma2/ incomplete gamma function].
 175 Note that this method relies
 176 on keeping a table of all the p[sub n ] previously computed, which does limit
 177 the precision of the method, depending upon the size of the table used.
 178
 179 When /a/ and /b/ are both small integers, then we can relate the incomplete
 180 beta to the binomial distribution and use the following finite sum:
 181
 182 [equation ibeta12]
 183
 184 Finally we can sidestep difficult areas, or move to an area with a more
 185 efficient means of computation, by using the duplication formulae:
 186
 187 [equation ibeta10]
 188
 189 [equation ibeta11]
 190
 191 The domains of a, b and x for which the various methods are used are identical
 192 to those described in the
 193 [@http://portal.acm.org/citation.cfm?doid=131766.131776 Didonato and Morris TOMS 708 paper].
 194
 195 [endsect][/section:ibeta_function The Incomplete Beta Function]
 196
 197 [/
 198   Copyright 2006 John Maddock and Paul A. Bristow.
 199   Distributed under the Boost Software License, Version 1.0.
 200   (See accompanying file LICENSE_1_0.txt or copy at
 201   http://www.boost.org/LICENSE_1_0.txt).
 202 ]