## Copyright (C) 2016 Marco Caliari ## Copyright (C) 2016 Nir Krakauer ## ## This file is part of Octave. ## ## Octave is free software; you can redistribute it and/or modify it ## under the terms of the GNU General Public License as published by ## the Free Software Foundation; either version 3 of the License, or (at ## your option) any later version. ## ## Octave is distributed in the hope that it will be useful, but ## WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ## General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {} {} gammainc (@var{x}, @var{a}) ## @deftypefnx {} {} gammainc (@var{x}, @var{a}, \"lower\") ## @deftypefnx {} {} gammainc (@var{x}, @var{a}, \"upper\") ## @deftypefnx {} {} gammainc (@var{x}, @var{a}, \"scaledlower\") ## @deftypefnx {} {} gammainc (@var{x}, @var{a}, \"scaledupper\") ## Compute the normalized incomplete gamma function. ## ## This is defined as ## @tex ## $$ ## \\gamma (x, a) = {1 \\over {\\Gamma (a)}}\\displaystyle{\\int_0^x t^{a-1} e^{-t} dt} ## $$ ## @end tex ## @ifnottex ## ## @example ## @group ## x ## 1 / ## gammainc (x, a) = --------- | exp (-t) t^(a-1) dt ## gamma (a) / ## t=0 ## @end group ## @end example ## ## @end ifnottex ## with the limiting value of 1 as @var{x} approaches infinity. ## The standard notation is @math{P(a,x)}, e.g., @nospell{Abramowitz} and ## @nospell{Stegun} (6.5.1). ## ## If @var{a} is scalar, then @code{gammainc (@var{x}, @var{a})} is returned ## for each element of @var{x} and vice versa. ## ## If neither @var{x} nor @var{a} is scalar, the sizes of @var{x} and ## @var{a} must agree, and @code{gammainc} is applied element-by-element. ## ## By default the incomplete gamma function integrated from 0 to @var{x} is ## computed. If @qcode{\"upper\"} is given then the complementary function ## integrated from @var{x} to infinity is calculated. ## ## Algorithms used include the following. ## 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. ## ## @seealso{gamma, gammaln} ## @end deftypefn ## P(a,x) = gamma(a,x)/Gamma(a), upper ## 1-P(a,x)=Q(a,x)=Gamma(a,x)/Gamma(a), lower function y = ngammainc (x,a,tail) if (nargin == 2) tail = "lower"; endif if (! isscalar (x) || ! isscalar (a)) [err, x, a] = common_size (x, a); if (err > 0) error ("ngammainc: x, a must be of common size or scalars"); endif endif ## Initialize output array. y = nan (size (x)); ## Different x, a combinations are handled by different subfunctions. i_done = false (size (x)); # Track which elements have been calculated. ## Case 1: x == 0 ii = (x == 0); if any (ii(:)) y(ii) = ngammainc_x0 (tail); i_done(ii) = true; endif ## Case 2: a == 0 ii = ~i_done & (a == 0); if any (ii(:)) y(ii) = ngammainc_a0 (x(ii), tail); i_done(ii) = true; endif ## Case 3: a == 1 ii = ~i_done & (a == 1); if any (ii(:)) y(ii) = ngammainc_a1 (x(ii), tail); i_done(ii) = true; endif ## Case 4: positive integer a; exp (x) and a! both under 1/eps ii = ~i_done & (a > 1) & (a == fix (a)) & (x <= 36) & (a <= 18); if any (ii(:)) y(ii) = ngammainc_an (x(ii), a(ii), tail); i_done(ii) = true; endif ## Case 5: x relatively small ii = ~i_done & (x + 0.25 < a | x < 0); if any (ii(:)) y(ii) = ngammainc_s (x(ii), a(ii), tail); i_done(ii) = true; endif ## Case 6: x positive and large relative to a ii = ~i_done; if any (ii(:)) y(ii) = ngammainc_l (x(ii), a(ii), tail); i_done(ii) = true; endif endfunction ##subfunctions to handle each case: ## x == 0 function y = ngammainc_x0 (tail) y = 1; if (strcmpi (tail,"upper") || strcmpi (tail,"scaledupper")) y = 0; endif endfunction ## a == 0 function y = ngammainc_a0 (x, tail) if strcmpi (tail,"lower") y = 1; elseif strcmpi (tail,"scaledlower") y = exp (x); else y = 0; endif endfunction ## a == 1 function y = ngammainc_a1 (x, tail) if strcmpi (tail,"lower") y = 1 - exp (-x); elseif strcmpi (tail,"scaledlower") if (abs (x) < 1/2) y = expm1 (x) / x; else y = (exp (x) - 1) / x; endif elseif strcmpi (tail,"upper") y = exp (-x); else strcmpi (tail,"scaledupper") y = 1 / x; endif endfunction ## positive integer a; exp (x) and a! both under 1/eps ## uses closed-form expressions for nonnegative integer a ## -- http://mathworld.wolfram.com/IncompleteGammaFunction.html function y = ngammainc_an (x, a, tail) y = t = ones (size (x)); i = 1; while any (a(:) > i) jj = (a > i); t(jj) .*= (x(jj) / i); y(jj) += t(jj); i++; endwhile if strcmpi (tail,"upper") y .*= exp (-x); elseif strcmpi (tail,"lower") y = 1 - exp (-x) .* y; elseif strcmpi (tail,"scaledupper") y .*= D (a,x) .* exp (-x); elseif strcmpi (tail,"lower") y = D (a,x) .* (1 - exp (-x) .* y); endif endfunction ## x + 0.25 < a | x < 0 ## Numerical Recipes in Fortran 77 (6.2.5) ## series function y = ngammainc_s (x, a, tail) if (strcmpi (tail,"scaledlower") || strcmpi (tail,"scaledupper")) y = ones (size (x)); term = x ./ (a + 1); else ## Of course it is possible to scale at the end, but some tests fail. ## And try gammainc (1,1000), it take 0 iterations if you scale now. y = D (a,x); term = y .* x ./ (a + 1); endif n = 1; while any (abs (term(:)) > abs (y(:)) * eps) ## y can be zero from the beginning (gammainc (1,1000)) jj = abs (term) > abs (y) * eps; n += 1; y(jj) += term(jj); term(jj) .*= x(jj) ./ (a(jj) + n); endwhile if strcmpi (tail,"upper") y = 1 - y; elseif strcmpi (tail,"scaledupper") y = 1 ./ D (a,x) - y; endif endfunction ## x positive and large relative to a ## NRF77 (6.2.7) ## Gamma (a,x)/Gamma (a) ## Lentz's algorithm function y = ngammainc_l (x, a, tail) tiny = 2^(-100); y = tiny; C0 = y; D0 = 0; bj = x + 1 - 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 any (abs ((Deltaj(:) - 1) ./ y(:)) > eps) Dj = bj + aj .* D0; Cj = bj + aj ./ C0; Dj = 1 ./ Dj; Deltaj = Cj .* Dj; y .*= Deltaj; D0 = Dj; C0 = Cj; bj += 2; aj = -j * (j - a); j += 1; endwhile if strcmpi (tail,"upper") y .*= D (a,x); elseif strcmpi (tail,"lower") y = 1 - y .* D (a,x); elseif strcmpi (tail,"scaledlower") y = 1 ./ D (a,x) - y; endif endfunction function y = D (a,x) y = zeros (size (x)); j = x < 0; if any (j(:)) ## to prevent imaginary parts when x negative and a integer y(j) = x(j).^a(j) .* exp (-x(j) - gammaln (a(j) + 1)); endif j = x > 0; if any (j(:)) ## 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(j) = exp (a(j) .* log (x(j)) - x(j) - gammaln (a(j) + 1)); endif endfunction %!test %! assert (ngammainc([0 0],0),[1 1]); %!test %! assert (ngammainc (0,[0; 0],"upper"),[0; 0]); %!test %!error (ngammainc([0 0],[0; 0])); %!test %! assert (ngammainc ([1 10],[1000 0]),[0 1]); %!test %! assert (ngammainc ([3 2 36],[2 3 18], "upper"),[4/exp(3) 5*exp(-2) (4369755579265807723 / 2977975)/exp(36)]); %!test %! assert (ngammainc (10,10),1 - (5719087 / 567) * exp (-10),-eps); %!test %! assert (ngammainc (10,10, "upper"),(5719087 / 567) * exp (-10),-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 (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 (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); %!test %! assert (ngammainc (709,1000,"upper"), 0.999999999999999999999999543586098796378631194583297890832, -eps);