function y = betainc(x,a,b)
%BETAINC Incomplete beta function.
%   Y = BETAINC(X,Z,W) computes the incomplete beta function for
%   corresponding elements of X, Z, and W.  The elements of X must be
%   in the closed interval [0,1].  The arguments X, Z and W must all
%   be the same size (or any of them can be scalar).
%
%   The incomplete beta function is defined as
%
%     I_x(z,w) = 1./BETA(z,w) .* 
%                 integral from 0 to x of t.^(z-1) .* (1-t).^(w-1) dt
%
%   See also BETA, BETALN.

%   Ref: Abramowitz & Stegun, Handbook of Mathematical Functions, sec. 26.5.

%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 5.9 $  $Date: 1998/12/18 16:29:04 $

if any(x(:) < 0 | x(:) > 1) | ~isreal(x),
   error('X must be in the interval [0,1].')
end

if nargin < 3
   error('Requires three input arguments.')
end

try
  % Preallocate y (using the size rules for plus)
  y = x + a + b; 
  y(:) = 0;
catch
  error('X, A and B must all the same size (or any of them can be scalar).')
end

if ~isempty(y)
  bt = exp(gammaln(a+b)-gammaln(a)-gammaln(b) + ...
       a.*log(x+(x==0)) + b.*log(1-x+(x==1)));

  k = find(x < (a+1) ./ (a+b+2));
  if length(k) > 0
     if length(x) == 1, xk = x; else, xk = x(k); end
     if length(a) == 1, ak = a; else, ak = a(k); end
     if length(b) == 1, bk = b; else, bk = b(k); end
     y(k) = bt(k) .* betacore(xk,ak,bk) ./ ak;
  end

  k = find(x >= (a+1) ./ (a+b+2));
  if length(k) > 0
     if length(x) == 1, xk = x; else, xk = x(k); end
     if length(a) == 1, ak = a; else, ak = a(k); end
     if length(b) == 1, bk = b; else, bk = b(k); end
     y(k) = 1 - bt(k) .* betacore(1-xk,bk,ak) ./ bk;
  end

  k = find(x == 0);
  if length(k) > 0
     if length(x) == 1, y(:) = 0; else, y(k) = 0; end
  end

  k = find(x == 1);
  if length(k) > 0
     if length(x) == 1, y(:) = 1; else, y(k) = 1; end
  end
end