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