...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
#include <boost/math/special_functions/beta.hpp>
namespace boost{ namespace math{ template <class T1, class T2, class T3> calculatedresulttype ibeta(T1 a, T2 b, T3 x); template <class T1, class T2, class T3, class Policy> calculatedresulttype ibeta(T1 a, T2 b, T3 x, const Policy&); template <class T1, class T2, class T3> calculatedresulttype ibetac(T1 a, T2 b, T3 x); template <class T1, class T2, class T3, class Policy> calculatedresulttype ibetac(T1 a, T2 b, T3 x, const Policy&); template <class T1, class T2, class T3> calculatedresulttype beta(T1 a, T2 b, T3 x); template <class T1, class T2, class T3, class Policy> calculatedresulttype beta(T1 a, T2 b, T3 x, const Policy&); template <class T1, class T2, class T3> calculatedresulttype betac(T1 a, T2 b, T3 x); template <class T1, class T2, class T3, class Policy> calculatedresulttype betac(T1 a, T2 b, T3 x, const Policy&); }} // namespaces
There are four incomplete beta functions : two are normalised versions (also known as regularized beta functions) that return values in the range [0, 1], and two are nonnormalised and return values in the range [0, beta(a, b)]. Users interested in statistical applications should use the normalised (or regularized ) versions (ibeta and ibetac).
All of these functions require 0 <= x <= 1.
The normalized functions ibeta and ibetac require a,b >= 0, and in addition that not both a and b are zero.
The functions beta and betac require a,b > 0.
The return type of these functions is computed using the result type calculation rules when T1, T2 and T3 are different types.
The final Policy argument is optional and can be used to control the behaviour of the function: how it handles errors, what level of precision to use etc. Refer to the policy documentation for more details.
template <class T1, class T2, class T3> calculatedresulttype ibeta(T1 a, T2 b, T3 x); template <class T1, class T2, class T3, class Policy> calculatedresulttype ibeta(T1 a, T2 b, T3 x, const Policy&);
Returns the normalised incomplete beta function of a, b and x:
template <class T1, class T2, class T3> calculatedresulttype ibetac(T1 a, T2 b, T3 x); template <class T1, class T2, class T3, class Policy> calculatedresulttype ibetac(T1 a, T2 b, T3 x, const Policy&);
Returns the normalised complement of the incomplete beta function of a, b and x:
template <class T1, class T2, class T3> calculatedresulttype beta(T1 a, T2 b, T3 x); template <class T1, class T2, class T3, class Policy> calculatedresulttype beta(T1 a, T2 b, T3 x, const Policy&);
Returns the full (nonnormalised) incomplete beta function of a, b and x:
template <class T1, class T2, class T3> calculatedresulttype betac(T1 a, T2 b, T3 x); template <class T1, class T2, class T3, class Policy> calculatedresulttype betac(T1 a, T2 b, T3 x, const Policy&);
Returns the full (nonnormalised) complement of the incomplete beta function of a, b and x:
The following tables give peak and mean relative errors in over various domains of a, b and x, along with comparisons to the GSL1.9 and Cephes libraries. Note that only results for the widest floatingpoint type on the system are given as narrower types have effectively zero error.
Note that the results for 80 and 128bit long doubles are noticeably higher than for doubles: this is because the wider exponent range of these types allow more extreme test cases to be tested. For example expected results that are zero at double precision, may be finite but exceptionally small with the wider exponent range of the long double types.
Table 6.18. Error rates for ibeta
GNU C++ version 7.1.0 
GNU C++ version 7.1.0 
Sun compiler version 0x5150 
Microsoft Visual C++ version 14.1 


Incomplete Beta Function: Small Values 
Max = 0ε (Mean = 0ε) 
Max = 8.97ε (Mean = 2.09ε) 
Max = 21.3ε (Mean = 2.75ε) 
Max = 8.4ε (Mean = 1.93ε) 
Incomplete Beta Function: Medium Values 
Max = 0ε (Mean = 0ε) 
Max = 50ε (Mean = 12.1ε) 
Max = 124ε (Mean = 18.4ε) 
Max = 106ε (Mean = 16.3ε) 
Incomplete Beta Function: Large and Diverse Values 
Max = 1.26ε (Mean = 0.063ε) 
Max = 1.96e+04ε (Mean = 997ε) 
Max = 4.98e+04ε (Mean = 2.07e+03ε) 
Max = 1.32e+03ε (Mean = 68.5ε) 
Incomplete Beta Function: Small Integer Values 
Max = 0ε (Mean = 0ε) 
Max = 4.45ε (Mean = 0.814ε) 
Max = 44.5ε (Mean = 10.1ε) 
Max = 3.85ε (Mean = 0.791ε) 
Table 6.19. Error rates for ibetac
GNU C++ version 7.1.0 
GNU C++ version 7.1.0 
Sun compiler version 0x5150 
Microsoft Visual C++ version 14.1 


Incomplete Beta Function: Small Values 
Max = 0ε (Mean = 0ε) 
Max = 10.6ε (Mean = 2.22ε) 
Max = 13.8ε (Mean = 2.68ε) 
Max = 6.94ε (Mean = 1.71ε) 
Incomplete Beta Function: Medium Values 
Max = 0ε (Mean = 0ε) 
Max = 73.9ε (Mean = 11.2ε) 
Max = 132ε (Mean = 19.8ε) 
Max = 56.7ε (Mean = 14.3ε) 
Incomplete Beta Function: Large and Diverse Values 
Max = 0.981ε (Mean = 0.0573ε) 
Max = 3.45e+04ε (Mean = 1.32e+03ε) 
Max = 6.31e+04ε (Mean = 2.04e+03ε) 
Max = 1.88e+03ε (Mean = 82.7ε) 
Incomplete Beta Function: Small Integer Values 
Max = 0ε (Mean = 0ε) 
Max = 5.34ε (Mean = 1.11ε) 
Max = 107ε (Mean = 17.1ε) 
Max = 6.37ε (Mean = 1.03ε) 
Table 6.20. Error rates for beta (incomplete)
GNU C++ version 7.1.0 
GNU C++ version 7.1.0 
Sun compiler version 0x5150 
Microsoft Visual C++ version 14.1 


Incomplete Beta Function: Small Values 
Max = 0ε (Mean = 0ε) 
Max = 11.1ε (Mean = 2.32ε) 
Max = 18.7ε (Mean = 3.19ε) 
Max = 9.94ε (Mean = 2.17ε) 
Incomplete Beta Function: Medium Values 
Max = 0.568ε (Mean = 0.0254ε) 
Max = 69.2ε (Mean = 13.4ε) 
Max = 174ε (Mean = 25ε) 
Max = 90ε (Mean = 12.7ε) 
Incomplete Beta Function: Large and Diverse Values 
Max = 0.999ε (Mean = 0.0325ε) 
Max = 6.84e+04ε (Mean = 2.76e+03ε) 
Max = 6.86e+04ε (Mean = 2.79e+03ε) 
Max = 633ε (Mean = 29.7ε) 
Incomplete Beta Function: Small Integer Values 
Max = 0.786ε (Mean = 0.0323ε) 
Max = 11.6ε (Mean = 3.6ε) 
Max = 51.8ε (Mean = 11ε) 
Max = 26ε (Mean = 6.28ε) 
Table 6.21. Error rates for betac
GNU C++ version 7.1.0 
GNU C++ version 7.1.0 
Sun compiler version 0x5150 
Microsoft Visual C++ version 14.1 


Incomplete Beta Function: Small Values 
Max = 0.676ε (Mean = 0.0302ε) 
Max = 9.92ε (Mean = 2.3ε) 
Max = 11.2ε (Mean = 2.94ε) 
Max = 8.94ε (Mean = 2.06ε) 
Incomplete Beta Function: Medium Values 
Max = 0.949ε (Mean = 0.098ε) 
Max = 63.5ε (Mean = 13.5ε) 
Max = 97.6ε (Mean = 24.3ε) 
Max = 90.6ε (Mean = 14.8ε) 
Incomplete Beta Function: Large and Diverse Values 
Max = 1.12ε (Mean = 0.0458ε) 
Max = 1.05e+05ε (Mean = 5.45e+03ε) 
Max = 1.04e+05ε (Mean = 5.46e+03ε) 
Max = 3.72e+03ε (Mean = 113ε) 
Incomplete Beta Function: Small Integer Values 
Max = 0.586ε (Mean = 0.0314ε) 
Max = 11.1ε (Mean = 3.65ε) 
Max = 103ε (Mean = 17.4ε) 
Max = 26.2ε (Mean = 6.36ε) 
There are two sets of tests: spot tests compare values taken from Mathworld's online function evaluator with this implementation: they provide a basic "sanity check" for the implementation, with one spottest in each implementationdomain (see implementation notes below).
Accuracy tests use data generated at very high precision (with NTL RR class set at 1000bit precision), using the "textbook" continued fraction representation (refer to the first continued fraction in the implementation discussion below). Note that this continued fraction is not used in the implementation, and therefore we have test data that is fully independent of the code.
This implementation is closely based upon "Algorithm 708; Significant digit computation of the incomplete beta function ratios", DiDonato and Morris, ACM, 1992.
All four of these functions share a common implementation: this is passed both x and y, and can return either p or q where these are related by:
so at any point we can swap a for b, x for y and p for q if this results in a more favourable position. Generally such swaps are performed so that we always compute a value less than 0.9: when required this can then be subtracted from 1 without undue cancellation error.
The following continued fraction representation is found in many textbooks but is not used in this implementation  it's both slower and less accurate than the alternatives  however it is used to generate test data:
The following continued fraction is due to Didonato and Morris, and is used in this implementation when a and b are both greater than 1:
For smallish b and x then a series representation can be used:
When b << a then the transition from 0 to 1 occurs very close to x = 1 and some care has to be taken over the method of computation, in that case the following series representation is used:
Where Q(a,x) is an incomplete gamma function. Note that this method relies on keeping a table of all the p_{n } previously computed, which does limit the precision of the method, depending upon the size of the table used.
When a and b are both small integers, then we can relate the incomplete beta to the binomial distribution and use the following finite sum:
Finally we can sidestep difficult areas, or move to an area with a more efficient means of computation, by using the duplication formulae:
The domains of a, b and x for which the various methods are used are identical to those described in the Didonato and Morris TOMS 708 paper.