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 39There are four [@http://en.wikipedia.org/wiki/Incomplete_beta_function incomplete beta functions] 40: two are normalised versions (also known as ['regularized] beta functions) 41that return values in the range [0, 1], and two are non-normalised and 42return values in the range [0, __beta(a, b)]. Users interested in statistical 43applications should use the normalised (or 44[@http://mathworld.wolfram.com/RegularizedBetaFunction.html regularized] 45) versions (ibeta and ibetac). 46 47All of these functions require /0 <= x <= 1/. 48 49The normalized functions __ibeta and __ibetac require /a,b >= 0/, and in addition that 50not both /a/ and /b/ are zero. 51 52The functions __beta and __betac require /a,b > 0/. 53 54The return type of these functions is computed using the __arg_promotion_rules 55when 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 65Returns 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 77Returns 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 87Returns 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 97Returns the full (non-normalised) complement of the incomplete beta function of a, b and x: 98 99[equation ibeta2] 100 101[h4 Accuracy] 102 103The following tables give peak and mean relative errors in over various domains of 104a, b and x, along with comparisons to the __gsl and __cephes libraries. 105Note that only results for the widest floating-point type on the system are given as 106narrower types have __zero_error. 107 108Note that the results for 80 and 128-bit long doubles are noticeably higher than 109for doubles: this is because the wider exponent range of these types allow 110more extreme test cases to be tested. For example expected results that 111are zero at double precision, may be finite but exceptionally small with 112the 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 124There 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] 126with this implementation: they provide a basic "sanity check" 127for the implementation, with one spot-test in each implementation-domain 128(see implementation notes below). 129 130Accuracy 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), 132using the "textbook" continued fraction representation (refer to the first continued 133fraction in the implementation discussion below). 134Note that this continued fraction is /not/ used in the implementation, 135and therefore we have test data that is fully independent of the code. 136 137[h4 Implementation] 138 139This 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 142All four of these functions share a common implementation: this is passed both 143x and y, and can return either p or q where these are related by: 144 145[equation ibeta_inv5] 146 147so at any point we can swap a for b, x for y and p for q if this results in 148a more favourable position. Generally such swaps are performed so that we always 149compute a value less than 0.9: when required this can then be subtracted from 1 150without undue cancellation error. 151 152The following continued fraction representation is found in many textbooks 153but is not used in this implementation - it's both slower and less accurate than 154the alternatives - however it is used to generate test data: 155 156[equation ibeta5] 157 158The following continued fraction is due to [@http://portal.acm.org/citation.cfm?doid=131766.131776 Didonato and Morris], 159and is used in this implementation when a and b are both greater than 1: 160 161[equation ibeta6] 162 163For smallish b and x then a series representation can be used: 164 165[equation ibeta7] 166 167When b << a then the transition from 0 to 1 occurs very close to x = 1 168and some care has to be taken over the method of computation, in that case 169the following series representation is used: 170 171[equation ibeta8] 172[/[equation ibeta9]] 173 174Where Q(a,x) is an [@http://functions.wolfram.com/GammaBetaErf/Gamma2/ incomplete gamma function]. 175Note that this method relies 176on keeping a table of all the p[sub n ] previously computed, which does limit 177the precision of the method, depending upon the size of the table used. 178 179When /a/ and /b/ are both small integers, then we can relate the incomplete 180beta to the binomial distribution and use the following finite sum: 181 182[equation ibeta12] 183 184Finally we can sidestep difficult areas, or move to an area with a more 185efficient means of computation, by using the duplication formulae: 186 187[equation ibeta10] 188 189[equation ibeta11] 190 191The domains of a, b and x for which the various methods are used are identical 192to 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] 203