]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 | ] |