Incomplete beta function in MATLAB

Question

Mathematica is capable of evaluating the incomplete beta function

$$ \mathrm{Beta}[z,a,b] = \int_0^z u^{a-1}\left(1-u\right)^{b-1}\,du $$

even when the argument $z$ is negative. MATLAB's function betainc(z,a,b), however, only allows $z\in[0,1]$. Do there exist any work-arounds for this (analytical tricks, or better MATLAB functions)?

Answer 1

Another possibility is to evaluate the integral that defines the incomplete beta function using quadrature. In Matlab you can use the integral function (use quadgk if you have an older version of Matlab)

y = integral(@(t)t.^(a-1).*(1-t).^(b-1),0,z)./beta(a,b)

where the output is divided by beta(a,b) to obtain the regularized incomplete beta function matching betainc. This can be vectorized across z with

f = @(t)t.^(a-1).*(1-t).^(b-1);
for i = numel(z):-1:1
    y(i) = integral(f,0,z(i));
end
y = y./beta(a,b);

The result is about the same speed as using hypergeom on my computer.

Also, note the complete definition of the incomplete beta function and the regularized beta function given at the functions.Wolfram.com site. For certain special combinations of parameters, both the hypergeometric form, and the equivalent integral, need to be evaluated differently. Also, Matlab's beta function will produce errors in these situations as well, as it's based on gammaln under the hood in order to be numerically careful. The symbolic beta function can be used as a substitute: double(beta(sym(a),sym(b))) or mfun('beta',a,b).

Answer 2

You can use the equivalence with Gauss' hypergeometric function to write, instead of betainc(z,a,b),

z.^a/a.*hypergeom([a,1-b],a+1,z) / beta(a,b)

This requires scalar a and b, and is slower than betainc (for the values for which the latter is defined, that is).