## 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 = gammainc (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 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 (x,a) .* exp (-x); elseif strcmpi (tail,"lower") y = D (x,a) .* (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 (x,a); 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 (x,a) - 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 (x,a); elseif strcmpi (tail,"lower") y = 1 - y .* D (x,a); elseif strcmpi (tail,"scaledlower") y = 1 ./ D (x,a) - y; endif endfunction function y = D(x,a) ## Compute exp(-x)*x^a/Gamma(a+1) in a stable way for x and a large ## ## L. Knusel, Computation of the Chi-square and Poisson distribution, ## SIAM J. Sci. Stat. Comput., 7(3), 1986 ## which quotes Section 5, Abramowitz&Stegun 6.1.40, 6.1.41 athresh = 10; ## FIXME: can be better tuned? y = zeros(size(x)); i_done = false (size (x)); i_done(x == 0) = true; ii = ~i_done & (x > 0) & (a > athresh) & (a >= x); if any (ii(:)) lnGa = log (2 * pi * a(ii)) / 2 + 1 ./ (12 * a(ii)) - ... 1 ./ (360 * a(ii) .^ 3) + 1 ./ (1260 * a(ii) .^ 5) - ... 1 ./ (1680 * a(ii) .^ 7) + 1 ./ (1188 * a(ii) .^ 9)- ... 691 ./ (87360 * a(ii) .^ 11) + 1 ./ (156 * a(ii) .^ 13) - ... 3617 ./ (122400 * a(ii) .^ 15) + ... 43867 ./ (244188 * a(ii) .^ 17) - 174611 ./ (125400 * a(ii) .^ 19); lns = log1p ((a(ii) - x(ii)) ./ x(ii)); y(ii) = exp ((a(ii) - x(ii)) - a(ii) .* lns - lnGa); i_done(ii) = true; endif ii = ~i_done & (x > 0) & (a > athresh) & (a < x); if any (ii(:)) lnGa = log (2 * pi * a(ii)) / 2 + 1 ./ (12 * a(ii)) - ... 1 ./ (360 * a(ii) .^ 3) + 1 ./ (1260 * a(ii) .^ 5) - ... 1 ./ (1680 * a(ii) .^ 7) + 1 ./ (1188 * a(ii) .^ 9)- ... 691 ./ (87360 * a(ii) .^ 11) + 1 ./ (156 * a(ii) .^ 13) - ... 3617 ./ (122400 * a(ii) .^ 15) + ... 43867 ./ (244188 * a(ii) .^ 17) - 174611 ./ (125400 * a(ii) .^ 19); lns = -log1p ((x(ii) - a(ii)) ./ a(ii)); y(ii) = exp ((a(ii) - x(ii)) - a(ii) .* lns - lnGa); i_done(ii) = true; endif ii = ~i_done & ((x <= 0) | (a <= athresh)); if any (ii(:)) ## standard formula for a not so large y(ii) = exp (a(ii) .* log (x(ii)) - x(ii) - gammaln (a(ii) + 1)); i_done(ii) = true; endif ii = (x < 0) & (a == fix (a)); if any(ii(:)) ## remove spurious imaginary part y(ii) = real (y(ii)); endif endfunction %!test %! assert (gammainc([0 0],0),[1 1]); %!test %! assert (gammainc (0,[0; 0],"upper"),[0; 0]); %!test %!error (gammainc([0 0],[0; 0])); %!test %! assert (gammainc ([1 10],[1000 0]),[0 1]); %%!test %%! assert (gammainc ([3 2 36],[2 3 18], "upper"),[4/exp(3) 5*exp(-2) (4369755579265807723 / 2977975)/exp(36)]); %!test %! assert (gammainc (10,10),1 - (5719087 / 567) * exp (-10),-eps); %!test %! assert (gammainc (10,10, "upper"),(5719087 / 567) * exp (-10),-eps); %!test %! assert (gammainc (-10,10),3.112658265341493126871617e7,-2 * eps); %!test %! % Here Matlab fails %! assert (isreal (gammainc (-10,10)),true); %!test %! assert (gammainc (60,6,"upper"),6.18022358081160257327264261e-20,-10 * eps); %!test %! assert (gammainc (-10,10.1,"upper"),... %! -2.9582761911890713293e7-1i * 9.612022339061679758e6,-30 * eps); %!test %! % Here Matlab is better %! assert (gammainc (751,750,"upper"),0.4805914320558831327179457887,-12 * eps); %!test %! % Here Matlab fails %! assert (gammainc (-100,1,"upper"),exp (100),-eps); %!test %! assert (gammainc (200,200,"upper"),0.49059658199276367497217454,-4 * eps); %!test %! assert (gammainc (200,200),0.509403418007236325027825459574527043,-3 * eps); %!test %! assert (gammainc (200,200,"scaledupper"),17.3984438553791505135122900,-eps); %!test %! assert (gammainc (200,200,"scaledlower"),18.065406676779221643065,-6 * eps); %!test %! assert (gammainc (200,201,"upper"),0.5187943096786844967708766027,-2 * eps); %!test %! assert (gammainc(201,200,"upper"),0.46249244908276709524913736667,-7 * eps); %!test %! assert (gammainc (200,201,"scaledupper"),18.4904360746560462660798514,-eps); %!test %! % here we are very good (no D (x,a)) involved %! assert (gammainc(1000,1000.5,"scaledlower"),39.48467539583672271,-2* eps); %!test %! assert (gammainc (709,1000,"upper"), 0.99999999999999999999999954358, -eps);