function retval = gamm_inv(p,a,b)
% Usage: gamminvp = gamm_inv(p,a,b)
% Description: This program accepts matrix argument p of
% probabilities, scalars a and b and returns the inverse
% c.d.f. for the gamma distribution evaluated at p with
% parameters:
%   a = alpha
%   b = beta
% in the standard description so that the distribution has
% mean alpha*beta and variance alpha*beta^2.

% Contact: T. Shores at tshores1@math.unl.edu

if (nargin < 3)
  error('gamm_inv: need three arguments');
end
if ((a < 0) | (b < 0))
  error('gamm_inv: parameters a, b must be positive');
end
if (sum(sum((p > 1) | (p < 0))))
  error('gamm_inv: p must have values between zero and one');
end
tol = 512*eps;
maxiter = 32;
ctr = 0;
retval = p;
ndx0 = find(retval==0);
ndx1 = find(retval==1);
ndx = find(0<retval & retval<1);
if (length(ndx)>0)
  x = retval(ndx);
  pndx = x;
  dx = ones(size(pndx));
  while (sum(abs(dx)>tol*pndx) > 0 & ctr < maxiter)
    dx = (gamm_cdf(x,a,b)-pndx)./gamm_pdf(x,a,b);
    dx = dx./(1+(abs(dx)>4*x).*abs(dx));
    x = x - dx;
    ctr = ctr+1;
  end
retval(ndx) = x;
end;
retval(ndx1) = Inf;
retval(ndx0) = 0;
