function y = ngammainc(x,a,tail) % incomplete (upper) 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 % Regions roughly as in 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), upper % 1-P(a,x)=Q(a,x)=Gamma(a,x)/Gamma(a), lower if (nargin == 2) tail = 'lower'; end % some corner cases if (x == 0) y = 1; if (strcmp(tail,'upper') || strcmp(tail,'scaledupper')) y = 0; end elseif (a == 0) if strcmp(tail,'lower') y = 1; elseif strcmp(tail,'scaledlower') y = exp(x); else y = 0; end %elseif (a == fix(a)) % a integer %http://mathworld.wolfram.com/IncompleteGammaFunction.html elseif (a == 1) if strcmp(tail,'lower') y = 1-exp(-x); elseif strcmp(tail,'scaledlower') if (abs(x) < 1/2) y = expm1(x)/x; else y = (exp(x)-1)/x; end elseif strcmp(tail,'upper') y = exp(-x); else strcmp(tail,'scaledupper') y = 1/x; end elseif (x + 0.25 < a || x < 0) % NRF77 (6.2.5) % series if (strcmp(tail,'scaledlower') || strcmp(tail,'scaledupper')) y = 1; term = x/(a+1); else % of course it is possibble to scale at the end. Try gammainc(1,1000), % it take 0 iterations if you scale now. % moreover, some test fail y = D(a,x); term = y*x/(a+1); end n = 1; while (abs(term) > abs(y)*eps) % y can be zero from the beginning (gammainc(1,1000)) n = n+1; y = y+term; term = term*(x/(a+n)); end if strcmp(tail,'upper') y = 1-y; elseif strcmp(tail,'scaledupper') y = 1/D(a,x)-y; end else % NRF77 (6.2.7) % Gamma(a,x)/Gamma(a) % Lentz's algorithm % here x is positive tiny = 2^(-100); y = tiny; C0 = y; D0 = 0; bj = x+1-a; aj = a; % if 'lower' or 'upper', it should be % % aj = a*exp(a*log(x)-x-gammaln(a+1)) % % and no scaling at the end. But it suffers cancellation % Better start with aj = a always. aj = a; Deltaj = 0; j = 1; while (abs((Deltaj-1)/y) > eps) Dj = bj+aj*D0; Cj = bj+aj/C0; Dj = 1/Dj; Deltaj = Cj*Dj; y = y*Deltaj; D0 = Dj; C0 = Cj; bj = bj+2; aj = -j*(j-a); j = j+1; end if strcmp(tail,'upper') y = y*D(a,x); elseif strcmp(tail,'lower') y = 1-y*D(a,x); elseif strcmp(tail,'scaledlower') y = 1/D(a,x)-y; end end end function y = D(a,x) if (x == 0) y = 0; elseif (x < 0) % to prevent imaginary parts when x negative and a integer y = x^a*exp(-x-gammaln(a+1)); else % Here, Temme's formula should be used % N. M. Temme, A set of algorithms for the incomplete Gamma functions, % Probab. Engrg. Inform. Sci., 8, 1994. y = exp(a*log(x)-x-gammaln(a+1)); end end %!test %! assert(ngammainc(0,0),1) %!test %! assert(ngammainc(0,0,'upper'),0) %!test %! assert(ngammainc(10,10),0.5420702855281477916858351,-7*eps) %!test %! assert(ngammainc(-10,10),3.112658265341493126871617e7,-11*eps) %!test %! % Here Matlab fails %! assert(isreal(ngammainc(-10,10)),true) %!test %! assert(ngammainc(60,6,'upper'),6.18022358081160257327264261e-20,-10*eps) %!test %! assert(ngammainc(-10,10.1,'upper'),... %! -2.9582761911890713293e7-1i*9.612022339061679758e6,-30*eps) %!test %! % Here Matlab is better %! assert(ngammainc(751,750,'upper'),0.4805914320558831327179457887,-512*eps) %!test %! % Here Matlab fails %! assert(ngammainc(-100,1,'upper'),exp(100),-eps) %!test %! assert(ngammainc(10,0),1,-eps) %!test %! assert(ngammainc(200,200,'upper'),0.49059658199276367497217454,-39*eps) %!test %! assert(ngammainc(200,200),0.5094034180072363250278254595745270431549,-38*eps) %!test %! assert(ngammainc(200,200,'scaledupper'),17.3984438553791505135122900,-eps) %!test %! assert(ngammainc(200,200,'scaledlower'),18.0654066767792216430658593,-78*eps) %!test %! assert(ngammainc(200,201,'upper'),0.5187943096786844967708766027648,-808*eps) %!test %! assert(ngammainc(201,200,'upper'),0.4624924490827670952491373666731,-698*eps) %!test %! assert(ngammainc(1,1000),0) %!test %! assert(ngammainc(200,201,'scaledupper'),18.4904360746560462660798514700,-1678*eps) %!test %! % here we are very good (no D(a,x)) involved %! assert(ngammainc(1000,1000.5,'scaledlower'),39.4846753958367227115137,-2*eps)