Mercurial > jwe > octave
view scripts/specfun/gammainc.m @ 31170:06c76eca5f95 default tip @
maint: merge stable to default
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Fri, 29 Jul 2022 12:02:42 -0400 |
parents | 5d3faba0342e |
children |
line wrap: on
line source
######################################################################## ## ## Copyright (C) 2016-2022 The Octave Project Developers ## ## See the file COPYRIGHT.md in the top-level directory of this ## distribution or <https://octave.org/copyright/>. ## ## 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 ## <https://www.gnu.org/licenses/>. ## ######################################################################## ## -*- texinfo -*- ## @deftypefn {} {@var{y} =} gammainc (@var{x}, @var{a}) ## @deftypefnx {} {@var{y} =} gammainc (@var{x}, @var{a}, @var{tail}) ## 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 then the sizes of @var{x} and ## @var{a} must agree, and @code{gammainc} is applied element-by-element. ## The elements of @var{a} must be non-negative. ## ## By default, @var{tail} is @qcode{"lower"} and the incomplete gamma function ## integrated from 0 to @var{x} is computed. If @var{tail} is @qcode{"upper"} ## then the complementary function integrated from @var{x} to infinity is ## calculated. ## ## If @var{tail} is @qcode{"scaledlower"}, then the lower incomplete gamma ## function is multiplied by ## @tex ## $\Gamma(a+1)\exp(x)x^{-a}$. ## @end tex ## @ifnottex ## @math{gamma(a+1)*exp(x)/(x^a)}. ## @end ifnottex ## If @var{tail} is @qcode{"scaledupper"}, then the upper incomplete gamma ## function is multiplied by the same quantity. ## ## References: ## ## @nospell{M. Abramowitz and I.A. Stegun}, ## @cite{Handbook of mathematical functions}, ## @nospell{Dover publications, Inc.}, 1972. ## ## @nospell{W. Gautschi}, ## @cite{A computational procedure for incomplete gamma functions}, ## @nospell{ACM Trans.@: Math Software}, pp.@: 466--481, Vol 5, No.@: 4, 2012. ## ## @nospell{W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery}, ## @cite{Numerical Recipes in Fortran 77}, ch.@: 6.2, Vol 1, 1992. ## ## @seealso{gamma, gammaincinv, 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 = "lower") if (nargin < 2) print_usage (); endif [err, x, a] = common_size (x, a); if (err > 0) error ("gammainc: X and A must be of common size or scalars"); endif if (iscomplex (x) || iscomplex (a)) error ("gammainc: all inputs must be real"); endif ## Remember original shape of data, but convert to column vector for calcs. x_sz = size (x); x = x(:); a = a(:); if (any (a < 0)) error ("gammainc: A must be non-negative"); endif if (nargin == 3 && ! any (strcmpi (tail, {"lower","upper","scaledlower","scaledupper"}))) error ("gammainc: invalid value for TAIL"); endif tail = tolower (tail); ## If any of the arguments is single then the output should be as well. if (strcmp (class (x), "single") || strcmp (class (a), "single")) x = single (x); a = single (a); endif ## Convert to floating point if necessary if (isinteger (x)) x = double (x); endif if (isinteger (a)) a = double (a); endif ## Initialize output array y = zeros (x_sz, class (x)); ## Different x, a combinations are handled by different subfunctions. todo = true (size (x)); # Track which elements need to be calculated. ## Case 0: x == Inf, a == Inf idx = (x == Inf) & (a == Inf); if (any (idx)) y(idx) = NaN; todo(idx) = false; endif ## Case 1: x == 0, a == 0. idx = (x == 0) & (a == 0); if (any (idx)) y(idx) = gammainc_00 (tail); todo(idx) = false; endif ## Case 2: x == 0. idx = todo & (x == 0); if (any (idx)) y(idx) = gammainc_x0 (tail); todo(idx) = false; endif ## Case 3: x = Inf idx = todo & (x == Inf); if (any (idx)) y(idx) = gammainc_x_inf (tail); todo(idx) = false; endif ## Case 4: a = Inf idx = todo & (a == Inf); if (any (idx)) y(idx) = gammainc_a_inf (tail); todo(idx) = false; endif ## Case 5: a == 0. idx = todo & (a == 0); if (any (idx)) y(idx) = gammainc_a0 (x(idx), tail); todo(idx) = false; endif ## Case 6: a == 1. idx = todo & (a == 1); if (any (idx)) y(idx) = gammainc_a1 (x(idx), tail); todo(idx) = false; endif ## Case 7: positive integer a; exp (x) and a! both under 1/eps. idx = (todo & (a == fix (a)) & (a > 1) & (a <= 18) & (x <= 36) & (abs (x) >= .1)); if (any (idx)) y(idx) = gammainc_an (x(idx), a(idx), tail); todo(idx) = false; endif ## For a < 2, x < 0, we increment a by 2 and use a recurrence formula after ## the computations. flag_a_small = todo & (abs (a) > 0) & (abs (a) < 2) & (x < 0); a(flag_a_small) += 2; flag_s = (((x + 0.25 < a) | (x < 0)) & (x > -20)) | (abs (x) < 1); ## Case 8: x, a relatively small. idx = todo & flag_s; if (any (idx)) y(idx) = gammainc_s (x(idx), a(idx), tail); todo(idx) = false; endif ## Case 9: x positive and large relative to a. idx = todo; if (any (idx)) y(idx) = gammainc_l (x(idx), a(idx), tail); todo(idx) = false; endif if (any (flag_a_small)) if (strcmp (tail, "lower")) y(flag_a_small) += D (x(flag_a_small), a(flag_a_small) - 1) + ... D (x(flag_a_small), a(flag_a_small) - 2); elseif (strcmp (tail, "upper")) y(flag_a_small) -= D (x(flag_a_small), a(flag_a_small) - 1) + ... D (x(flag_a_small), a(flag_a_small) - 2); elseif (strcmp (tail, "scaledlower")) y(flag_a_small) = y(flag_a_small) .* (x(flag_a_small) .^ 2) ./ ... (a(flag_a_small) .* (a(flag_a_small) - 1)) + (x(flag_a_small) ./ ... (a(flag_a_small) - 1)) + 1; elseif (strcmp (tail, "scaledupper")) y(flag_a_small) = y(flag_a_small) .* (x(flag_a_small) .^ 2) ./ ... (a(flag_a_small) .* (a(flag_a_small) - 1)) - (x(flag_a_small) ./ ... (a(flag_a_small) - 1)) - 1; endif endif endfunction ## Subfunctions to handle each case: ## x == 0, a == 0. function y = gammainc_00 (tail) if (strcmp (tail, "upper") || strcmp (tail, "scaledupper")) y = 0; else y = 1; endif endfunction ## x == 0. function y = gammainc_x0 (tail) if (strcmp (tail, "lower")) y = 0; elseif (strcmp (tail, "upper") || strcmp (tail, "scaledlower")) y = 1; else y = Inf; endif endfunction ## x == Inf. function y = gammainc_x_inf (tail) if (strcmp (tail, "lower")) y = 1; elseif (strcmp (tail, "upper") || strcmp (tail, "scaledupper")) y = 0; else y = Inf; endif endfunction ## a == Inf. function y = gammainc_a_inf (tail) if (strcmp (tail, "lower")) y = 0; elseif (strcmp (tail, "upper") || strcmp (tail, "scaledlower")) y = 1; else y = Inf; endif endfunction ## a == 0. function y = gammainc_a0 (x, tail) if (strcmp (tail, "lower")) y = 1; elseif (strcmp (tail, "scaledlower")) y = exp (x); else y = 0; endif endfunction ## a == 1. function y = gammainc_a1 (x, tail) if (strcmp (tail, "lower")) if (abs (x) < 1/2) y = - expm1 (-x); else y = 1 - exp (-x); endif elseif (strcmp (tail, "upper")) y = exp (-x); elseif (strcmp (tail, "scaledlower")) if (abs (x) < 1/2) y = expm1 (x) ./ x; else y = (exp (x) - 1) ./ x; endif 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 = gammainc_an (x, a, tail) y = t = ones (size (x), class (x)); i = 1; while (any (a(:) > i)) jj = (a > i); t(jj) .*= (x(jj) / i); y(jj) += t(jj); i++; endwhile if (strcmp (tail, "lower")) y = 1 - exp (-x) .* y; elseif (strcmp (tail, "upper")) y .*= exp (-x); elseif (strcmp (tail, "scaledlower")) y = (1 - exp (-x) .* y) ./ D(x, a); elseif (strcmp (tail, "scaledupper")) y .*= exp (-x) ./ D(x, a); endif endfunction ## x + 0.25 < a | x < 0 | abs(x) < 1. ## Numerical Recipes in Fortran 77 (6.2.5) ## series function y = gammainc_s (x, a, tail) if (strcmp (tail, "scaledlower") || strcmp (tail, "scaledupper")) y = ones (size (x), class (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 (strcmp (tail, "upper")) y = 1 - y; elseif (strcmp (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 ## __gammainc__ in libinterp/corefcn/__gammainc__.cc function y = gammainc_l (x, a, tail) y = __gammainc__ (x, a); if (strcmp (tail, "lower")) y = 1 - y .* D (x, a); elseif (strcmp (tail, "upper")) y .*= D (x, a); elseif (strcmp (tail, "scaledlower")) y = 1 ./ D (x, a) - y; endif endfunction ## 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. function y = D (x, a) athresh = 10; # FIXME: can this be better tuned? y = zeros (size (x), class (x)); todo = true (size (x)); todo(x == 0) = false; ii = todo & (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); todo(ii) = false; endif ii = todo & (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); todo(ii) = false; endif ii = todo & ((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)); todo(ii) = false; endif ii = (x < 0) & (a == fix (a)); if (any (ii)) # remove spurious imaginary part. y(ii) = real (y(ii)); endif endfunction ## Test: case 1,2,5 %!assert (gammainc ([0, 0, 1], [0, 1, 0]), [1, 0, 1]) %!assert (gammainc ([0, 0, 1], [0, 1, 0], "upper"), [0, 1, 0]) %!assert (gammainc ([0, 0, 1], [0, 1, 0], "scaledlower"), [1, 1, exp(1)]) %!assert (gammainc ([0, 0, 1], [0, 1, 0], "scaledupper"), [0, Inf, 0]) ## Test: case 3,4 %!assert (gammainc ([2, Inf], [Inf, 2]), [0, 1]) %!assert (gammainc ([2, Inf], [Inf, 2], "upper"), [1, 0]) %!assert (gammainc ([2, Inf], [Inf, 2], "scaledlower"), [1, Inf]) %!assert (gammainc ([2, Inf], [Inf, 2], "scaledupper"), [Inf, 0]) ## Test: case 5 ## Matlab fails for this test %!assert (gammainc (-100,1,"upper"), exp (100), -eps) ## Test: case 6 %!assert (gammainc ([1, 2, 3], 1), 1 - exp (-[1, 2, 3])) %!assert (gammainc ([1, 2, 3], 1, "upper"), exp (- [1, 2, 3])) %!assert (gammainc ([1, 2, 3], 1, "scaledlower"), ... %! (exp ([1, 2, 3]) - 1) ./ [1, 2, 3]) %!assert (gammainc ([1, 2, 3], 1, "scaledupper"), 1 ./ [1, 2, 3]) ## Test: case 7 %!assert (gammainc (2, 2, "lower"), 0.593994150290162, -2e-15) %!assert (gammainc (2, 2, "upper"), 0.406005849709838, -2e-15) %!assert (gammainc (2, 2, "scaledlower"), 2.194528049465325, -2e-15) %!assert (gammainc (2, 2, "scaledupper"), 1.500000000000000, -2e-15) %!assert (gammainc ([3 2 36],[2 3 18], "upper"), ... %! [4/exp(3) 5*exp(-2) (4369755579265807723 / 2977975)/exp(36)], -eps) %!assert (gammainc (10, 10), 1 - (5719087 / 567) * exp (-10), -eps) %!assert (gammainc (10, 10, "upper"), (5719087 / 567) * exp (-10), -eps) ## Test: case 8 %!assert (gammainc (-10, 10), 3.112658265341493126871617e7, -2*eps) ## Matlab fails this next one%! %! %!assert (isreal (gammainc (-10, 10)), true) %!assert (gammainc (-10, 10.1, "upper"), ... %! -2.9582761911890713293e7-1i * 9.612022339061679758e6, -30*eps) %!assert (gammainc (-10, 10, "upper"), -3.112658165341493126871616e7, ... %! -2*eps) %!assert (gammainc (-10, 10, "scaledlower"), 0.5128019364747265, -1e-14) %!assert (gammainc (-10, 10, "scaledupper"), -0.5128019200000000, -1e-14) %!assert (gammainc (200, 201, "upper"), 0.518794309678684497, -2 * eps) %!assert (gammainc (200, 201, "scaledupper"), %! 18.4904360746560462660798514, -eps) ## Here we are very good (no D (x,a)) involved %!assert (gammainc (1000, 1000.5, "scaledlower"), 39.48467539583672271, -2*eps) %!assert (gammainc (709, 1000, "upper"), 0.99999999999999999999999954358, -eps) ## Test: case 9 %!test <*47800> %! assert (gammainc (60, 6, "upper"), 6.18022358081160257327264261e-20, %! -10*eps); ## Matlab is better here than Octave %!assert (gammainc (751, 750, "upper"), 0.4805914320558831327179457887, -12*eps) %!assert (gammainc (200, 200, "upper"), 0.49059658199276367497217454, -6*eps) %!assert (gammainc (200, 200), 0.509403418007236325027825459574527043, -5*eps) %!assert (gammainc (200, 200, "scaledupper"), 17.3984438553791505135122900, %! -3*eps) %!assert (gammainc (200, 200, "scaledlower"), 18.065406676779221643065, -8*eps) %!assert (gammainc (201, 200, "upper"), 0.46249244908276709524913736667, %! -7*eps) %!assert <*54550> (gammainc (77, 2), 1) %!assert (gammainc (77, 2, "upper"), 0, -eps) %!assert (gammainc (1000, 3.1), 1) %!assert (gammainc (1000, 3.1, "upper"), 0) ## Test small argument %!assert (gammainc ([1e-05, 1e-07,1e-10,1e-14], 0.1), ... %! [0.33239840504050, 0.20972940370977, 0.10511370061022, ... %! 0.041846517936723], 1e-13); %!assert (gammainc ([1e-05, 1e-07,1e-10,1e-14], 0.2), ... %! [0.10891226058559, 0.043358823442178, 0.010891244210402, ... %! 0.0017261458806785], 1e-13); %!test %!assert (gammainc ([1e-02, 1e-03, 1e-5, 1e-9, 1e-14], 0.9), ... %! [0.016401189184068, 0.0020735998660840, 0.000032879756964708, ... %! 8.2590606569241e-9, 2.6117443021738e-13], -1e-12); %!test %!assert (gammainc ([1e-02, 1e-03, 1e-5, 1e-9, 1e-14], 2), ... %! [0.0000496679133402659, 4.99666791633340e-7, 4.99996666679167e-11, ... %! 4.99999999666667e-19, 4.99999999999997e-29], -1e-12); %!test <*53543> %! y_exp = 9.995001666250085e-04; %! assert (gammainc (1/1000, 1), y_exp, -eps); %!test <53612> %! assert (gammainc (-20, 1.1, "upper"), ... %! 6.50986687074979e8 + 2.11518396291149e8*i, -1e-13); ## Test conservation of the class (five tests for each subroutine). %!assert (class (gammainc (0, 1)) == "double") %!assert (class (gammainc (single (0), 1)) == "single") %!assert (class (gammainc (int8 (0), 1)) == "double") %!assert (class (gammainc (0, single (1))) == "single") %!assert (class (gammainc (0, int8 (1))) == "double") %!assert (class (gammainc (1, 0)) == "double") %!assert (class (gammainc (single (1), 0)) == "single") %!assert (class (gammainc (int8 (1), 0)) == "double") %!assert (class (gammainc (1, single (0))) == "single") %!assert (class (gammainc (1, int8 (0))) == "double") %!assert (class (gammainc (1, 1)) == "double") %!assert (class (gammainc (single (1), 1)) == "single") %!assert (class (gammainc (int8 (1), 1)) == "double") %!assert (class (gammainc (1, single (1))) == "single") %!assert (class (gammainc (1, int8 (1))) == "double") %!assert (class (gammainc (1, 2)) == "double") %!assert (class (gammainc (single (1), 2)) == "single") %!assert (class (gammainc (int8 (1), 2)) == "double") %!assert (class (gammainc (1, single (2))) == "single") %!assert (class (gammainc (1, int8 (2))) == "double") %!assert (class (gammainc (-1, 0.5)) == "double") %!assert (class (gammainc (single (-1), 0.5)) == "single") %!assert (class (gammainc (int8 (-1), 0.5)) == "double") %!assert (class (gammainc (-1, single (0.5))) == "single") %!assert (class (gammainc (-1, int8 (0.5))) == "double") %!assert (class (gammainc (1, 0.5)) == "double") %!assert (class (gammainc (single (1), 0.5)) == "single") %!assert (class (gammainc (int8 (1), 0.5)) == "double") %!assert (class (gammainc (1, single (0.5))) == "single") %!assert (class (gammainc (1, int8 (0.5))) == "double") ## Test input validation %!error <Invalid call> gammainc () %!error <Invalid call> gammainc (1) %!error <must be of common size or scalars> gammainc ([0, 0],[0; 0]) %!error <must be of common size or scalars> gammainc ([1 2 3], [1 2]) %!error <all inputs must be real> gammainc (2+i, 1) %!error <all inputs must be real> gammainc (1, 2+i) %!error <A must be non-negative> gammainc (1, [0, -1, 1]) %!error <A must be non-negative> %! a = ones (2,2,2); %! a(1,1,2) = -1; %! gammainc (1, a); %!error <invalid value for TAIL> gammainc (1,2, "foobar")