[SciPy-User] incomplete beta function B(x; k+1, 0) ?

[Pauli Virtanen]

Sat, 26 Jun 2010 17:59:09 -0400, josef.pktd wrote:
> Is there a incomplete beta function with a zero argument  B(x; k+1, 0)
>  available in scipy?
> 
>>>> special.betainc( np.arange(5)+1, 0, 0.5)
> array([ NaN,  NaN,  NaN,  NaN,  NaN])

>>> print special.betainc.__doc__
betainc(x1, x2, x3[, out])
y=betainc(a,b,x) returns the incomplete beta integral of the
arguments, evaluated from zero to x: gamma(a+b) / (gamma(a)*gamma(b))
* integral(t**(a-1) (1-t)**(b-1), t=0..x).

So the function in Scipy has an additional normalization prefactor. The 
prefactor however seems to be zero for b=0, so this function in Scipy 
probably can't easily be used for finding the value of the integral at 
b=0.

But if you just set b=1e-99 (or anything smaller than the machine 
epsilon) and divide the prefactor away, you should end up with the result 
you are looking for.

At b=0 the integral is logarithmically divergent towards x->1, and a 
finite b > 0 cuts this off around x ~ 1 - exp(-1/b) --- so it shouldn't 
matter. The rest of the integrand should also saturate at b ~ machine 
epsilon.

But apparently, there's something fishy in the betainc algorithm for x 
close to 1 and b close to 0, for example this discontinuity:

	>>> sc.betainc(1, 1.000002e-6, 1 - 1e-6)
	1.3815442755027441e-05
	>>> sc.betainc(1, 1.000001e-6, 1 - 1e-6)
	1.2397031262149499e-05

So while things in principle are OK, you shouldn't probably trust things 
beyond x > 1 - 1e-3.

-- 
Pauli Virtanen