Mercurial > octave
view scripts/specfun/gammainc.m @ 24904:740df3fce3fb stable
New strategy to compute the incomplete gamma function (see bug #47800).
--
added libinterp/corefcn/__gammainc_lentz__.cc
added scripts/specfun/gammainc.m
changed libinterp/corefcn/module.mk
changed liboctave/external/slatec-fn/module.mk
changed liboctave/numeric/lo-specfun.cc
changed scripts/specfun/module.mk
removed liboctave/external/slatec-fn/dgami.f
removed liboctave/external/slatec-fn/dgamit.f
removed liboctave/external/slatec-fn/gami.f
removed liboctave/external/slatec-fn/gamit.f
removed liboctave/external/slatec-fn/xdgami.f
removed liboctave/external/slatec-fn/xdgamit.f
removed liboctave/external/slatec-fn/xgmainc.f
removed liboctave/external/slatec-fn/xsgmainc.f
author | Michele Ginesi <michele.ginesi@gmail.com> |
---|---|
date | Thu, 07 Sep 2017 15:47:33 +0200 |
parents | |
children | 662faf9de127 |
line wrap: on
line source
## Copyright (C) 2016 Marco Caliari ## Copyright (C) 2016 Nir Krakauer ## Copyright (C) 2017 Michele Ginesi ## ## 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 ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {} {} gammainc (@var{x}, @var{a}) ## @deftypefnx {} {} 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, 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 nonnegative. ## ## By default or if @var{tail} is @qcode{"lower"} 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 divided by the same quantity. ## ## References: ## ## @nospell{M. Abramowitz and I. 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, 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 >= 4 || nargin <= 1) print_usage (); endif if ((! isscalar (x)) || (! isscalar (a))) [err, x, a] = common_size (x, a); if (err > 0) error ("gammainc: x, a must be of common size or scalars"); endif endif if (any (a < 0) || any (imag (a) != 0)) error ("gammainc: a must be real and non negative"); endif if (any (imag (x) != 0)) error ("gammainc: x must be real"); endif ## Initialize output array. if (isinteger (x)) x = double (x); endif if (strcmpi (class (a), "single")) x = single (x); endif y = zeros (size (x), class (x)); ## Different x, a combinations are handled by different subfunctions. i_done = false (size (x)); # Track which elements have been calculated. ## Case 0: x == Inf, a == Inf ii = ((x == Inf) & (a == Inf)); if (any (ii(:))) y(ii) = NaN; i_done(ii) = true; endif ## Case 1: x == 0, a == 0. ii = ((x == 0) & (a == 0)); if (any (ii(:))) y(ii) = gammainc_00 (tail); i_done(ii) = true; endif ## Case 2: x == 0. ii = ((! i_done) & (x == 0)); if (any (ii(:))) y(ii) = gammainc_x0 (tail); i_done(ii) = true; endif ## Case 3: x = Inf ii = ((! i_done) & (x == Inf)); if (any (ii(:))) y(ii) = gammainc_x_inf (tail); i_done(ii) = true; endif ## Case 4: a = Inf ii = ((! i_done) & (a == Inf)); if (any (ii(:))) y(ii) = gammainc_a_inf (tail); i_done(ii) = true; endif ## Case 5: a == 0. ii = ((! i_done) & (a == 0)); if (any (ii(:))) y(ii) = gammainc_a0 (x(ii), tail); i_done(ii) = true; endif ## Case 6: a == 1. ii = ((! i_done) & (a == 1)); if (any (ii(:))) y(ii) = gammainc_a1 (x(ii), tail); i_done(ii) = true; endif flag_an = ((a > 1) & (a == fix (a)) & (x <= 36) & (abs (x) >= 1e-01) & ... (a <= 18)); ## Case 7: positive integer a; exp (x) and a! both under 1/eps. ii = ((! i_done) & flag_an); if (any (ii(:))) y(ii) = gammainc_an (x(ii), a(ii), tail); i_done(ii) = true; endif ## For a < 2, x < 0, we increment a by 2 and use a recurrence formula after ## the computations. flag_a_small = ((abs (a) < 2) & (abs(a) > 0) & (! i_done) & (x < 0)); a(flag_a_small) += 2; flag_s = (((x + 0.25 < a | x < 0 | a < 5) & (x > -20)) | (abs (x) < 1)); ## Case 8: x, a relatively small. ii = ((! i_done) & flag_s); if (any (ii(:))) y(ii) = gammainc_s (x(ii), a(ii), tail); i_done(ii) = true; endif ## Case 9: x positive and large relative to a. ii = (! i_done); if (any (ii(:))) y(ii) = gammainc_l (x(ii), a(ii), tail); i_done(ii) = true; endif if (any (flag_a_small)) if (strcmpi (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 (strcmpi (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 (strcmpi (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 (strcmpi (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 (strcmpi (tail, "upper") || strcmpi (tail, "scaledupper")) y = 0; else y = 1; endif endfunction ## x == 0. function y = gammainc_x0 (tail) if (strcmpi (tail, "upper") || strcmpi (tail, "scaledlower")) y = 1; elseif (strcmpi (tail, "lower")) y = 0; else y = Inf; endif endfunction ## x == Inf. function y = gammainc_x_inf (tail) if (strcmpi (tail, "lower")) y = 1; elseif (strcmpi (tail, "upper") || strcmpi (tail, "scaledupper")) y = 0; else y = Inf; endif endfunction ## a == Inf. function y = gammainc_a_inf (tail) if (strcmpi (tail, "lower")) y = 0; elseif (strcmpi (tail, "upper") || strcmpi (tail, "scaledlower")) y = 1; else y = Inf; endif endfunction ## a == 0. function y = gammainc_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 = gammainc_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 = 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 (strcmpi (tail, "upper")) y .*= exp (-x); elseif (strcmpi (tail, "lower")) y = 1 - exp (-x) .* y; elseif (strcmpi (tail, "scaledupper")) y .*= exp (-x) ./ D(x, a); elseif (strcmpi (tail, "scaledlower")) y = (1 - exp (-x) .* y) ./ D(x, a); endif endfunction ## x + 0.25 < a | x < 0 | x not real | abs(x) < 1 | a < 5. ## Numerical Recipes in Fortran 77 (6.2.5) ## series function y = gammainc_s (x, a, tail) if (strcmpi (tail, "scaledlower") || strcmpi (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 (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 ## __gammainc_lentz__ in libinterp/corefcn/__gammainc_lentz__.cc function y = gammainc_l (x, a, tail) n = numel (x); y = zeros (size (x), class (x)); for i = 1:n y(i) = __gammainc_lentz__ (x(i), a(i)); endfor 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), class (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: case 1,2,5 %!test %! assert (gammainc ([0, 0, 1], [0, 1, 0]), [1, 0, 1]); %!test %! assert (gammainc ([0, 0, 1], [0, 1, 0], "upper"), [0, 1, 0]); %!test %! assert (gammainc ([0, 0, 1], [0, 1, 0], "scaledlower"), [1, 1, exp(1)]); %!test %! assert (gammainc ([0, 0, 1], [0, 1, 0], "scaledupper"), [0, Inf, 0]); ## Test: case 3,4 %!test %! assert (gammainc ([2, Inf], [Inf, 2]), [0, 1]); %!test %! assert (gammainc ([2, Inf], [Inf, 2], "upper"), [1, 0]); %!test %! assert (gammainc ([2, Inf], [Inf, 2], "scaledlower"), [1, Inf]); %!test %! assert (gammainc ([2, Inf], [Inf, 2], "scaledupper"), [Inf, 0]); ## Test: case 5 %!test %! % Here Matlab fails %! assert (gammainc (-100,1,"upper"),exp (100),-eps); ## Test: case 6 %!test %! assert (gammainc ([1, 2, 3], 1), 1 - exp (-[1, 2, 3])); %!test %! assert (gammainc ([1, 2, 3], 1, "upper"), exp (- [1, 2, 3])); %!test %! 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 %!test %! assert (gammainc (2, 2, "lower"), 0.593994150290162, -2e-15); %!test %! assert (gammainc (2, 2, "upper"), 0.406005849709838, -2e-15); %!test %! assert (gammainc (2, 2, "scaledlower"), 2.194528049465325, -2e-15); %!test %! assert (gammainc (2, 2, "scaledupper"), 1.500000000000000, -2e-15); %!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: case 8 %!test %! assert (gammainc (-10, 10), 3.112658265341493126871617e7, -2 * eps); %!test %! % Here Matlab fails %! assert (isreal (gammainc (-10, 10)), true); %!test %! assert (gammainc (-10, 10.1, "upper"),... %! -2.9582761911890713293e7-1i * 9.612022339061679758e6,-30 * eps); %!test %! assert (gammainc (-10, 10, "upper"), -3.112658165341493126871616e7, ... %! -2 * eps); %!test %! assert (gammainc (-10, 10, "scaledlower"), ... %! 0.5128019364747265,-1e-14); %!test %! assert (gammainc (-10, 10, "scaledupper"), ... %! -0.5128019200000000, -1e-14); %!test %! assert (gammainc (200, 201, "upper"),... %! 0.518794309678684497, -2 * 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); ## Test: case 9 %!test <47800> %! assert (gammainc (60, 6, "upper"), ... %! 6.18022358081160257327264261e-20, -10 * eps); %!test %! % Here Matlab is better %! assert (gammainc (751, 750, "upper"),... %! 0.4805914320558831327179457887, -12 * 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 (201, 200, "upper"), 0.46249244908276709524913736667,... %! -7 * eps); ## Test small argument %!test %! assert (gammainc ([1e-05, 1e-07,1e-10,1e-14], 0.1), ... %! [0.33239840504050, 0.20972940370977, 0.10511370061022, ... %! 0.041846517936723], 1e-13); %!test %! 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); %!xtest %! assert (gammainc (-20, 1.1, "upper"),... %! 6.50986687074979e8 + 2.11518396291149e8*i, -1e-13); ## Test on the 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 gammainc () %!error gammainc (1) %!error gammainc (1,2,3,4) %!error gammainc (1, [0, 1i,1]) %!error gammainc (1, [0, -1, 1]) %!error gammainc ([0 0],[0; 0]) %!error gammainc ([1 2 3], [1 2]) %!error gammainc (1i,1) %!error gammainc (1,1i)