1######################################################################## 2## 3## Copyright (C) 2017-2021 The Octave Project Developers 4## 5## See the file COPYRIGHT.md in the top-level directory of this 6## distribution or <https://octave.org/copyright/>. 7## 8## This file is part of Octave. 9## 10## Octave is free software: you can redistribute it and/or modify it 11## under the terms of the GNU General Public License as published by 12## the Free Software Foundation; either version 3 of the License, or 13## (at your option) any later version. 14## 15## Octave is distributed in the hope that it will be useful, but 16## WITHOUT ANY WARRANTY; without even the implied warranty of 17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18## GNU General Public License for more details. 19## 20## You should have received a copy of the GNU General Public License 21## along with Octave; see the file COPYING. If not, see 22## <https://www.gnu.org/licenses/>. 23## 24######################################################################## 25 26## -*- texinfo -*- 27## @deftypefn {} {} gammaincinv (@var{y}, @var{a}) 28## @deftypefnx {} {} gammaincinv (@var{y}, @var{a}, @var{tail}) 29## Compute the inverse of the normalized incomplete gamma function. 30## 31## The normalized incomplete gamma function is defined as 32## @tex 33## $$ 34## \gamma (x, a) = {1 \over {\Gamma (a)}}\displaystyle{\int_0^x t^{a-1} e^{-t} dt} 35## $$ 36## @end tex 37## @ifnottex 38## 39## @example 40## @group 41## x 42## 1 / 43## gammainc (x, a) = --------- | exp (-t) t^(a-1) dt 44## gamma (a) / 45## t=0 46## @end group 47## @end example 48## 49## @end ifnottex 50## 51## and @code{gammaincinv (gammainc (@var{x}, @var{a}), @var{a}) = @var{x}} 52## for each non-negative value of @var{x}. If @var{a} is scalar then 53## @code{gammaincinv (@var{y}, @var{a})} is returned for each element of 54## @var{y} and vice versa. 55## 56## If neither @var{y} nor @var{a} is scalar then the sizes of @var{y} and 57## @var{a} must agree, and @code{gammaincinv} is applied element-by-element. 58## The variable @var{y} must be in the interval @math{[0,1]} while @var{a} must 59## be real and positive. 60## 61## By default, @var{tail} is @qcode{"lower"} and the inverse of the incomplete 62## gamma function integrated from 0 to @var{x} is computed. If @var{tail} is 63## @qcode{"upper"}, then the complementary function integrated from @var{x} to 64## infinity is inverted. 65## 66## The function is computed with Newton's method by solving 67## @tex 68## $$ 69## y - \gamma (x, a) = 0 70## $$ 71## @end tex 72## @ifnottex 73## 74## @example 75## @var{y} - gammainc (@var{x}, @var{a}) = 0 76## @end example 77## 78## @end ifnottex 79## 80## Reference: @nospell{A. Gil, J. Segura, and N. M. Temme}, @cite{Efficient and 81## accurate algorithms for the computation and inversion of the incomplete 82## gamma function ratios}, @nospell{SIAM J. Sci.@: Computing}, pp.@: 83## A2965--A2981, Vol 34, 2012. 84## 85## @seealso{gammainc, gamma, gammaln} 86## @end deftypefn 87 88function x = gammaincinv (y, a, tail = "lower") 89 90 if (nargin < 2 || nargin > 3) 91 print_usage (); 92 endif 93 94 [err, y, a] = common_size (y, a); 95 if (err > 0) 96 error ("gammaincinv: Y and A must be of common size or scalars"); 97 endif 98 99 if (iscomplex (y) || iscomplex (a)) 100 error ("gammaincinv: all inputs must be real"); 101 endif 102 103 ## Remember original shape of data, but convert to column vector for calcs. 104 orig_sz = size (y); 105 y = y(:); 106 a = a(:); 107 108 if (any ((y < 0) | (y > 1))) 109 error ("gammaincinv: Y must be in the range [0, 1]"); 110 endif 111 112 if (any (a <= 0)) 113 error ("gammaincinv: A must be strictly positive"); 114 endif 115 116 ## If any of the arguments is single then the output should be as well. 117 if (strcmp (class (y), "single") || strcmp (class (a), "single")) 118 y = single (y); 119 a = single (a); 120 endif 121 122 ## Convert to floating point if necessary 123 if (isinteger (y)) 124 y = double (y); 125 endif 126 if (isinteger (a)) 127 a = double (a); 128 endif 129 130 ## Initialize output array 131 x = zeros (size (y), class (y)); 132 133 maxit = 20; 134 tol = eps (class (y)); 135 136 ## Special cases, a = 1 or y = 0, 1. 137 138 if (strcmpi (tail, "lower")) 139 x(a == 1) = - log1p (- y(a == 1)); 140 x(y == 0) = 0; 141 x(y == 1) = Inf; 142 p = y; 143 q = 1 - p; 144 elseif (strcmpi (tail, "upper")) 145 x(a == 1) = - log (y(a == 1)); 146 x(y == 0) = Inf; 147 x(y == 1) = 0; 148 q = y; 149 p = 1 - q; 150 else 151 error ("gammaincinv: invalid value for TAIL") 152 endif 153 154 todo = (a != 1) & (y != 0) & (y != 1); 155 156 ## Case 1: p small. 157 158 i_flag_1 = todo & (p < ((0.2 * (1 + a)) .^ a) ./ gamma (1 + a)); 159 160 if (any (i_flag_1)) 161 aa = a(i_flag_1); 162 pp = p(i_flag_1); 163 164 ## Initial guess. 165 166 r = (pp .* gamma (1 + aa)) .^ (1 ./ aa); 167 168 c2 = 1 ./ (aa + 1); 169 c3 = (3 * aa + 5) ./ (2 * (aa + 1) .^2 .* (aa + 2)); 170 c4 = (8 * aa .^ 2 + 33 * aa + 31) ./ (3 * (aa + 1) .^ 3 .* (aa + 2) .* ... 171 (aa + 3)); 172 c5 = (125 * aa .^ 4 + 1179 * aa .^ 3 + 3971 * aa.^2 + 5661 * aa + 2888) ... 173 ./ (24 * (1 + aa) .^4 .* (aa + 2) .^ 2 .* (aa + 3) .* (aa + 4)); 174 175 ## FIXME: Would polyval() be better here for more accuracy? 176 x0 = r + c2 .* r .^ 2 + c3 .* r .^ 3 + c4 .* r .^4 + c5 .* r .^ 5; 177 178 ## For this case we invert the lower version. 179 180 F = @(p, a, x) p - gammainc (x, a, "lower"); 181 JF = @(a, x) - exp (- gammaln (a) - x + (a - 1) .* log (x)); 182 x(i_flag_1) = newton_method (F, JF, pp, aa, x0, tol, maxit); 183 endif 184 185 todo(i_flag_1) = false; 186 187 ## Case 2: q small. 188 189 i_flag_2 = (q < exp (- 0.5 * a) ./ gamma (1 + a)) & (a > 0) & (a < 10); 190 i_flag_2 &= todo; 191 192 if (any (i_flag_2)) 193 aa = a(i_flag_2); 194 qq = q(i_flag_2); 195 196 ## Initial guess. 197 198 x0 = (-log (qq) - gammaln (aa)); 199 200 ## For this case, we invert the upper version. 201 202 F = @(q, a, x) q - gammainc (x, a, "upper"); 203 JF = @(a, x) exp (- gammaln (a) - x) .* x .^ (a - 1); 204 x(i_flag_2) = newton_method (F, JF, qq, aa, x0, tol, maxit); 205 endif 206 207 todo(i_flag_2) = false; 208 209 ## Case 3: a small. 210 211 i_flag_3 = todo & ((a > 0) & (a < 1)); 212 213 if (any (i_flag_3)) 214 aa = a(i_flag_3); 215 pp = p(i_flag_3); 216 217 ## Initial guess 218 219 xl = (pp .* gamma (aa + 1)) .^ (1 ./ aa); 220 x0 = xl; 221 222 ## For this case, we invert the lower version. 223 224 F = @(p, a, x) p - gammainc (x, a, "lower"); 225 JF = @(a, x) - exp (-gammaln (a) - x) .* x .^ (a - 1); 226 x(i_flag_3) = newton_method (F, JF, pp, aa, x0, tol, maxit); 227 endif 228 229 todo(i_flag_3) = false; 230 231 ## Case 4: a large. 232 233 i_flag_4 = todo; 234 235 if (any (i_flag_4)) 236 aa = a(i_flag_4); 237 qq = q(i_flag_4); 238 239 ## Initial guess 240 241 d = 1 ./ (9 * aa); 242 t = 1 - d + sqrt (2) * erfcinv (2 * qq) .* sqrt (d); 243 x0 = aa .* (t .^ 3); 244 245 ## For this case, we invert the upper version. 246 247 F = @(q, a, x) q - gammainc (x, a, "upper"); 248 JF = @(a, x) exp (- gammaln (a) - x + (a - 1) .* log (x)); 249 x(i_flag_4) = newton_method (F, JF, qq, aa, x0, tol, maxit); 250 endif 251 252 ## Restore original shape 253 x = reshape (x, orig_sz); 254 255endfunction 256 257## subfunction: Newton's Method 258function x = newton_method (F, JF, y, a, x0, tol, maxit); 259 l = numel (y); 260 res = -F (y, a, x0) ./ JF (a, x0); 261 todo = (abs (res) >= tol * abs (x0)); 262 x = x0; 263 it = 0; 264 while (any (todo) && (it++ < maxit)) 265 x(todo) += res(todo); 266 res(todo) = -F (y(todo), a(todo), x(todo)) ./ JF (a(todo), x(todo)); 267 todo = (abs (res) >= tol * abs (x)); 268 endwhile 269 x += res; 270endfunction 271 272 273%!test 274%! x = [1e-10, 1e-09, 1e-08, 1e-07]; 275%! a = [2, 3, 4]; 276%! [x, a] = ndgrid (x, a); 277%! xx = gammainc (gammaincinv (x, a), a); 278%! assert (xx, x, -3e-14); 279 280%!test 281%! x = [1e-10, 1e-09, 1e-08, 1e-07]; 282%! a = [2, 3, 4]; 283%! [x, a] = ndgrid (x, a); 284%! xx = gammainc (gammaincinv (x, a, "upper"), a, "upper"); 285%! assert (xx, x, -3e-14); 286 287%!test 288%! x = linspace (0, 1)'; 289%! a = [linspace(0.1, 1, 10), 2:5]; 290%! [x, a] = ndgrid (x, a); 291%! xx = gammainc (gammaincinv (x, a), a); 292%! assert (xx, x, -1e-13); 293 294%!test 295%! x = linspace (0, 1)'; 296%! a = [linspace(0.1, 1, 10), 2:5]; 297%! [x, a] = ndgrid (x, a); 298%! xx = gammainc (gammaincinv (x, a, "upper"), a, "upper"); 299%! assert (xx, x, -1e-13); 300 301%!test <*56453> 302%! assert (gammaincinv (1e-15, 1) * 2, 2e-15, -1e-15); 303%! assert (gammaincinv (1e-16, 1) * 2, 2e-16, -1e-15); 304 305## Test the conservation of the input class 306%!assert (class (gammaincinv (0.5, 1)), "double") 307%!assert (class (gammaincinv (single (0.5), 1)), "single") 308%!assert (class (gammaincinv (0.5, single (1))), "single") 309%!assert (class (gammaincinv (int8 (0), 1)), "double") 310%!assert (class (gammaincinv (0.5, int8 (1))), "double") 311%!assert (class (gammaincinv (int8 (0), single (1))), "single") 312%!assert (class (gammaincinv (single (0.5), int8 (1))), "single") 313 314## Test input validation 315%!error gammaincinv () 316%!error gammaincinv (1) 317%!error gammaincinv (1, 2, 3, 4) 318%!error <must be of common size or scalars> 319%! gammaincinv (ones (2,2), ones (1,2), 1); 320%!error <all inputs must be real> gammaincinv (0.5i, 1) 321%!error <all inputs must be real> gammaincinv (0, 1i) 322%!error <Y must be in the range \[0, 1\]> gammaincinv (-0.1,1) 323%!error <Y must be in the range \[0, 1\]> gammaincinv (1.1,1) 324%!error <Y must be in the range \[0, 1\]> 325%! y = ones (1, 1, 2); 326%! y(1,1,2) = -1; 327%! gammaincinv (y,1); 328%!error <A must be strictly positive> gammaincinv (0.5, 0) 329%!error <A must be strictly positive> 330%! a = ones (1, 1, 2); 331%! a(1,1,2) = 0; 332%! gammaincinv (1,a,1); 333%!error <invalid value for TAIL> gammaincinv (1,2, "foobar") 334