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 Mathemtical Functions, sec. 26.5.

%   Copyright (c) 1984-97 by The MathWorks, Inc.
%   $Revision: 5.5 $  $Date: 1997/04/08 06:41:38 $

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

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

siz = max([size(x);size(a);size(b)]);

if (any(siz~=size(x)) & prod(size(x))~=1) | ...
   (any(siz~=size(a)) & prod(size(a))~=1) | ...
   (any(siz~=size(b)) & prod(size(b))~=1)
  error("X, A and B must all the same size (or any of them can be scalar).")
end

y = zeros(siz);
bt = exp(gammaln(a+b)-gammaln(a)-gammaln(b) + a.*log(x) + b.*log(1-x));

k = find(x < (a+1) ./ (a+b+2));
if prod(size(k)) > 0
   if prod(size(x)) == 1
      xk = x;
   else
      xk = x(k);
   end
   if prod(size(a)) == 1
      ak = a;
   else
      ak = a(k);
   end
   if prod(size(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 prod(size(k)) > 0
   if prod(size(x)) == 1
      xk = x;
   else
      xk = x(k);
   end
   if prod(size(a)) == 1
      ak = a;
   else
      ak = a(k);
   end
   if prod(size(b)) == 1; bk = b; else; bk = b(k); end
   y(k) = 1 - bt(k) .* betacore(1-xk,bk,ak) ./ bk;
end

end
