ibeta.qbk

Boost provides free peer-reviewed portable C++ source libraries. We emphasize libraries that work

[section:ibeta_function Incomplete Beta Functions]

[h4 Synopsis]

``
#include <boost/math/special_functions/gamma.hpp>
``

   namespace boost{ namespace math{
   
   template <class T1, class T2, class T3>
   ``__sf_result`` ibeta(T1 a, T2 b, T3 x);
   
   template <class T1, class T2, class T3, class ``__Policy``>
   ``__sf_result`` ibeta(T1 a, T2 b, T3 x, const ``__Policy``&);
   
   template <class T1, class T2, class T3>
   ``__sf_result`` ibetac(T1 a, T2 b, T3 x);
   
   template <class T1, class T2, class T3, class ``__Policy``>
   ``__sf_result`` ibetac(T1 a, T2 b, T3 x, const ``__Policy``&);
   
   template <class T1, class T2, class T3>
   ``__sf_result`` beta(T1 a, T2 b, T3 x);
   
   template <class T1, class T2, class T3, class ``__Policy``>
   ``__sf_result`` beta(T1 a, T2 b, T3 x, const ``__Policy``&);
   
   template <class T1, class T2, class T3>
   ``__sf_result`` betac(T1 a, T2 b, T3 x);
   
   template <class T1, class T2, class T3, class ``__Policy``>
   ``__sf_result`` betac(T1 a, T2 b, T3 x, const ``__Policy``&);
   
   }} // namespaces
   
[h4 Description]

There are four [@http://en.wikipedia.org/wiki/Incomplete_beta_function incomplete beta functions]
: two are normalised versions (also known as ['regularized] beta functions)
that return values in the range [0, 1], and two are non-normalised and
return values in the range [0, __beta(a, b)].  Users interested in statistical
applications should use the normalised (or
[@http://mathworld.wolfram.com/RegularizedBetaFunction.html regularized]
) versions (ibeta and ibetac).

All of these functions require /a > 0/, /b > 0/ and /0 <= x <= 1/.

The return type of these functions is computed using the __arg_pomotion_rules
when T1, T2 and T3 are different types.

[optional_policy]

   template <class T1, class T2, class T3>
   ``__sf_result`` ibeta(T1 a, T2 b, T3 x);

   template <class T1, class T2, class T3, class ``__Policy``>
   ``__sf_result`` ibeta(T1 a, T2 b, T3 x, const ``__Policy``&);

Returns the normalised incomplete beta function of a, b and x:

[equation ibeta3]

[graph ibeta]

   template <class T1, class T2, class T3>
   ``__sf_result`` ibetac(T1 a, T2 b, T3 x);
   
   template <class T1, class T2, class T3, class ``__Policy``>
   ``__sf_result`` ibetac(T1 a, T2 b, T3 x, const ``__Policy``&);
   
Returns the normalised complement of the incomplete beta function of a, b and x:

[equation ibeta4]

   template <class T1, class T2, class T3>
   ``__sf_result`` beta(T1 a, T2 b, T3 x);

   template <class T1, class T2, class T3, class ``__Policy``>
   ``__sf_result`` beta(T1 a, T2 b, T3 x, const ``__Policy``&);

Returns the full (non-normalised) incomplete beta function of a, b and x:

[equation ibeta1]

   template <class T1, class T2, class T3>
   ``__sf_result`` betac(T1 a, T2 b, T3 x);

   template <class T1, class T2, class T3, class ``__Policy``>
   ``__sf_result`` betac(T1 a, T2 b, T3 x, const ``__Policy``&);

Returns the full (non-normalised) complement of the incomplete beta function of a, b and x:

[equation ibeta2]

[h4 Accuracy]

The following tables give peak and mean relative errors in over various domains of
a, b and x, along with comparisons to the __gsl and __cephes libraries.  
Note that only results for the widest floating-point type on the system are given as
narrower types have __zero_error.

Note that the results for 80 and 128-bit 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 Errors In the Function ibeta(a,b,x)
[[Significand Size] [Platform and Compiler] [0 < a,b < 10

and

0 < x < 1]  [0 < a,b < 100

and

0 < x < 1][1x10[super -5] < a,b < 1x10[super 5]

and

0 < x < 1]]
[[53] [Win32, Visual C++ 8] 
      [Peak=42.3 Mean=2.9

(GSL Peak=682 Mean=32.5)

(__cephes Peak=42.7 Mean=7.0)]  
      [Peak=108 Mean=16.6
      
(GSL Peak=690 Mean=151)

(__cephes Peak=1545 Mean=218)] 
      [Peak=4x10[super 3][space] Mean=203
      
(GSL Peak~3x10[super 5][space] Mean~2x10[super 4][space])

(__cephes Peak~5x10[super 5][space] Mean~2x10[super 4][space])]]
[[64] [Redhat Linux IA32, gcc-3.4.4] [Peak=21.9 Mean=3.1] 
      [Peak=270.7 Mean=26.8] [Peak~5x10[super 4][space] Mean=3x10[super 3][space] ]]
[[64] [Redhat Linux IA64, gcc-3.4.4] [Peak=15.4 Mean=3.0] [Peak=112.9 Mean=14.3] 
      [Peak~5x10[super 4][space] Mean=3x10[super 3][space]]]
[[113] [HPUX IA64, aCC A.06.06] [Peak=20.9 Mean=2.6] [Peak=88.1 Mean=14.3] 
      [Peak~2x10[super 4][space] Mean=1x10[super 3][space] ]]
]

[table Errors In the Function ibetac(a,b,x)
[[Significand Size] [Platform and Compiler] [0 < a,b < 10

and

0 < x < 1]  [0 < a,b < 100

and

0 < x < 1][1x10[super -5] < a,b < 1x10[super 5]

and

0 < x < 1]]
[[53] [Win32, Visual C++ 8] [Peak=13.9 Mean=2.0]  
      [Peak=56.2 Mean=14] [Peak=3x10[super 3][space] Mean=159]]
[[64] [Redhat Linux IA32, gcc-3.4.4] [Peak=21.1 Mean=3.6] 
      [Peak=221.7 Mean=25.8] 
      [Peak~9x10[super 4][space] Mean=3x10[super 3][space] ]]
[[64] [Redhat Linux IA64, gcc-3.4.4] [Peak=10.6 Mean=2.2] 
      [Peak=73.9 Mean=11.9] 
      [Peak~9x10[super 4][space] Mean=3x10[super 3][space] ]]
[[113] [HPUX IA64, aCC A.06.06] [Peak=9.9 Mean=2.6] 
      [Peak=117.7 Mean=15.1] 
      [Peak~3x10[super 4][space] Mean=1x10[super 3][space] ]]
]

[table Errors In the Function beta(a, b, x)
[[Significand Size] [Platform and Compiler] [0 < a,b < 10

and

0 < x < 1]  [0 < a,b < 100

and

0 < x < 1][1x10[super -5] < a,b < 1x10[super 5]

and

0 < x < 1]]
[[53] [Win32, Visual C++ 8] [Peak=39 Mean=2.9] [Peak=91 Mean=12.7] [Peak=635 Mean=25]]
[[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] ]]
[[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] ]]
[[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] ]]
]

[table Errors In the Function betac(a,b,x)
[[Significand Size] [Platform and Compiler] [0 < a,b < 10

and

0 < x < 1]  [0 < a,b < 100

and

0 < x < 1][1x10[super -5] < a,b < 1x10[super 5]

and

0 < x < 1]]
[[53] [Win32, Visual C++ 8] [Peak=12.0 Mean=2.4] [Peak=91 Mean=15] [Peak=4x10[super 3][space] Mean=113]]
[[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] ]]
[[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]  ]]
[[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]  ]]
]

[h4 Testing]

There are two sets of tests: spot tests compare values taken from 
[@http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp?name=BetaRegularized Mathworld's online function evaluator]
with this implementation: they provide a basic "sanity check"
for the implementation, with one spot-test in each implementation-domain 
(see implementation notes below). 
 
Accuracy tests use data generated at very high precision
(with [@http://shoup.net/ntl/doc/RR.txt NTL RR class] set at 1000-bit 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.  

[h4 Implementation]

This implementation is closely based upon 
[@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.]

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:

[equation ibeta_inv5]

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:

[equation ibeta5]

The following continued fraction is due to [@http://portal.acm.org/citation.cfm?doid=131766.131776 Didonato and Morris],
and is used in this implementation when a and b are both greater than 1:

[equation ibeta6]

For smallish b and x then a series representation can be used:

[equation ibeta7]

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:

[equation ibeta8]
[/[equation ibeta9]]

Where Q(a,x) is an [@http://functions.wolfram.com/GammaBetaErf/Gamma2/ incomplete gamma function].
Note that this method relies
on keeping a table of all the p[sub 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:

[equation ibeta12]

Finally we can sidestep difficult areas, or move to an area with a more
efficient means of computation, by using the duplication formulae:

[equation ibeta10]

[equation ibeta11]

The domains of a, b and x for which the various methods are used are identical
to those described in the
[@http://portal.acm.org/citation.cfm?doid=131766.131776 Didonato and Morris TOMS 708 paper].

[endsect][/section:ibeta_function The Incomplete Beta Function]

[/ 
  Copyright 2006 John Maddock and Paul A. Bristow.
  Distributed under the Boost Software License, Version 1.0.
  (See accompanying file LICENSE_1_0.txt or copy at
  http://www.boost.org/LICENSE_1_0.txt).
]