function y = newgammainc(x,a,tail) % incomplete gamma function % % P(a,x) = 1/Gamma(a) * int_x^\infty exp(-t)*t^(a-1)dt % % a>=0, x real % % tail is either upper or lower % Numerical Recipes in Fortran77, \S 6.2 %Gautschi, W. (1979). A computational procedure for incomplete gamma functions. ACM Transactions on Mathematical Software (TOMS), 5(4), 466-481. % P(a,x) = gamma(a,x)/Gamma(a) % 1-P(a,x)=Q(a,x)=Gamma(a,x)/Gamma(a) %parse the input flag if (nargin == 2) tail = 'lower'; end if (numel(tail) >= 6) && strcmp(tail(1:6),'scaled') scaled = true; %return P or Q times scale = gamma(a+1)*exp(x)*x^-a tail = tail(7:end); else scaled = false; end if (x > 1.5) && (x > 20) & (x < 300) && ((a - x) <= 0.25) %Gautschi, Sec. 4.3 [Legendre continued fraction for y = (x - a + 1) * x^-a * exp(x) * Gamma(a,x) -- works for large positive inputs, but unfortunately incurs a good bit of roundoff error] p = 0; q = (x - a - 1) * (x - a + 1); r = 4 * (x - a + 1); s = 1 - a; rho = 0; t = 1; y = t; while abs(t / y) > eps p += s; q += r; r += 8; s += 2; tau = p * (1 + rho); rho = tau / (q - tau); t *= rho; y += t; end if scaled y *= a / (x - a + 1); %scale*Q if strcmp(tail, 'lower') y = exp(x - a*log(x) + gammaln(a+1)) - y; %scale*P end else y *= exp(-x + a*log(x) - gammaln(a)) / (x - a + 1); %Q if strcmp(tail, 'lower') y = 1 - y; %P end end return end if (x >= 0) && ((a - x) > 0.25) %Gautschi, Sec. 5 [power series for y = gamma(a+1) * x^-a * exp(x) * P = scale*P] t = 1; y = t; n = 0; while abs(t / y) > eps n++; t *= x / (a + n); y += t; end if scaled if strcmp(tail, 'upper') y = exp(x - a*log(x) + gammaln(a+1)) - y; %scale*Q end else y *= exp(-x + a*log(x) - gammaln(a+1)); %P if strcmp(tail, 'upper') y = 1 - y; %Q end end return end if (abs(x) <= a+1 || x <= 0) % some corner cases if (x == 0) y = 0; elseif (a == 0) y = 1; elseif (a == 1) y = 1-exp(-x); elseif (abs(x) >= 300) % to be fixed % (6.2.5) % series gamma(a,x)/Gamma(a) for large numbers % take exp(log(each term of the series)) % not very accurate, Kahan summation does not help term = exp(a*log(x)-x-gammaln(a+1)); y = 0; n = 0; while (abs(term) > eps*abs(y)) y = y+term; n = n+1; term = exp((a+n)*log(x)-x-gammaln(a+1+n)); end else % (6.2.5) % series if ~scaled gammaa = gamma(a); gamma1n = a*gammaa; xn = exp(-x)*x^a; % formula xn = exp(-x+a*log(x)) does not work with x<0, e.g. x=-10,a=10 % but for positive x it may help y = xn/gamma1n; gamma1n = (a+1)*gamma1n; xn = xn*x; n = 1; while (abs(xn/gamma1n) > eps*abs(y)) y = y+xn/gamma1n; gamma1n = (a+n+1)*gamma1n; xn = xn*x; n = n+1; end else % a scaled version was required % series, but compute a*gamma(a,x)/(exp(-x)*x^a) and simplify % Gamma(a)/Gamma(a+1+n) by 1/(a+n)! xn = 1; y = xn; rgammainv = (a+1); xn = xn*x; n = 1; while (abs(xn/rgammainv) > eps*abs(y)) y = y+xn/rgammainv; rgammainv = (a+n+1)*rgammainv; xn = xn*x; n = n+1; end end end if strcmp(tail,'upper') if ~scaled y = 1-y; else y = exp(-a*log(x)+x+gammaln(a+1))-y; end end else % (6.2.7) % Gamma(a,x)/Gamma(a) % Lentz's algorithm % here x is positive % scaled versions not implemented yet % Could it be that continued fraction (6.2.6), although slower, is % more accurate? tiny = 2^(-100); b0 = 0; f0 = tiny; C0 = f0; D0 = 0; bj = x+1-a; aj = 1; Deltaj = 0; j = 1; while (abs(Deltaj-1) > eps) Dj = bj+aj*D0; Cj = bj+aj/C0; Dj = 1/Dj; Deltaj = Cj*Dj; fj = f0*Deltaj; D0 = Dj; C0 = Cj; f0 = fj; bj = x+2*j+1-a; aj = -j*(j-a); j = j+1; end % try to change next line to % y = exp(-x-gammaln(a)+a*log(x))*fj; % one more test will fail! y = exp(a*log(x)-x-gammaln(a))*fj; if strcmp(tail,'lower') y = 1-y; end end %!test %! assert(newgammainc(0,0),0) %!test %! assert(newgammainc(0,0,'upper'),1) %!test %! assert(newgammainc(10,10),0.5420702855281477916858351,-3*eps) %!test %! % Here Matlab fails %! assert(newgammainc(-10,10),3.112658265341493126871617e7,-2*eps) %!test %! assert(newgammainc(60,6,'upper'),6.18022358081160257327264261e-20,-13*eps) %!test %! assert(newgammainc(-10,10.1,'upper'),... %! -2.9582761911890713293e7-1i*9.612022339061679758e6,-27*eps) %!test %! % Here Matlab is better %! assert(newgammainc(751,750,'upper'),0.4805914320558831327179457887,-1800*eps) %!test %! % Here Matlab fails %! assert(newgammainc(-100,1,'upper'),exp(100),-eps) %!test %! assert(newgammainc(10,0),1,-eps) %!test %! assert(newgammainc(200,200,'upper'),0.49059658199276367497217454,-300*eps) %!test %! assert(newgammainc(200,201,'upper'),0.518794309678684496770876602764835968845235402065399405783,-1000*eps) %!test %! assert(newgammainc(201,200,'upper'),0.462492449082767095249137366673156681609762469721457129052,-1000*eps)