function y = betainc(x,a,b)
%BETAINC Incomplete beta function.
%	y = betainc(x,a,b).
%	The elements of x must be in the interval [0,1].
%	a and b must be scalars.
%
%	See also BETA.

%	Ref: Abramowitz & Stegun, Handbook of Mathemtical Functions, sec. 26.5.
%	Copyright (c) 1984-93 by The MathWorks, Inc.

if any(any(x < 0 | x > 1))
        error('In betainc(X,a,b), X must be in the interval [0,1].')
end
if max(size(a)) > 1 | max(size(b)) > 1
        error('In betainc(X,a,b), a and b must be scalars.');
end
if nargin < 3
        error('In betainc(X,a,b), all three input arguments must be given.');
end

[m,n] = size(x);
y = zeros(m,n);
bt = y;
k = find(x > 0 & x < 1);
bt(k) = exp(gammaln(a+b)-gammaln(a)-gammaln(b) + a*log(x(k)) + b*log(1-x(k)));
k = find(x < (a+1)/(a+b+2));
y(k) = bt(k) .* betacore(x(k),a,b) / a;
k = find(x >= (a+1)/(a+b+2));
y(k) = 1 - bt(k) .* betacore(1-x(k),b,a) / b;
