5

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)?

Doubt
  • 1,767
  • 1
    Did you look what Mathematica's incomplete beta function's codes are? or Mathematica does not allow people to get a sneakpeak of the codes for built-in functions. – Shuhao Cao Jun 04 '13 at 22:13
  • Unfortunately no. I do not know if Mathematica discloses its algorithms for built-in functions. – Doubt Jun 05 '13 at 16:58
  • 2
    Note that Matlab's betainc is actually the regularized incomplete beta function. The Mathematica equivalent is actually BetaRegularized. Also, betainc is a purely numeric function. – horchler Sep 06 '13 at 17:14

2 Answers2

5

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).

horchler
  • 3,258
  • 2
  • 27
  • 41
4

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).

Luis Mendo
  • 1,894