Mercurial > octave

annotate scripts/specfun/gammaincinv.m @ 31548:c8ad083a5802 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Clean up m-files before Octave 8.1 release. * external.txi, oop.txi, Table.h, documentation.cc, gui-preferences-ed.h, lo-specfun.cc, range.tst : Eliminate triple newlines. * Map.m, MemoizedFunction.m, delaunayn.m, inputParser.m, __publish_latex_output__.m, publish.m, unpack.m, fminbnd.m, __add_default_menu__.m, gammainc.m, gallery.m, hadamard.m, weboptions.m: Add newline after keyword "function" or before keyword "endfunction" for readability. * getaudiodata.m, pkg.m : Add semicolon to end of line for error() statement. * movegui.m: Combine mutliple calls to set() into one for performance. * __unimplemented__.m (missing_functions): Remove missing functions that have been implemented. * __vectorize__.m, check_default_input.m, betaincinv.m, gammaincinv.m: Remove semicolon at end of line with "function" declaration. * weboptions.m: Remove semicolon at end of line with "if" keyword. * integrate_adaptive.m, factor.m: Use keyword "endif" rather than bare "end".
author Rik <rik@octave.org>
date Fri, 25 Nov 2022 21:23:54 -0800
parents 5d3faba0342e
children 597f3ee61a48
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
27923
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
1 ########################################################################
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
2 ##
30564
796f54d4ddbf update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents: 29359
diff changeset
3 ## Copyright (C) 2017-2022 The Octave Project Developers
27918
b442ec6dda5c use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents: 27800
diff changeset
4 ##
27923
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
5 ## See the file COPYRIGHT.md in the top-level directory of this
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
6 ## distribution or <https://octave.org/copyright/>.
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
7 ##
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
8 ## This file is part of Octave.
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
9 ##
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
10 ## Octave is free software: you can redistribute it and/or modify it
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
11 ## under the terms of the GNU General Public License as published by
28944
d39b09b4c5db maint: Use identical Copyright text in scripts/.
Rik <rik@octave.org>
parents: 28928
diff changeset
12 ## the Free Software Foundation, either version 3 of the License, or
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
13 ## (at your option) any later version.
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
14 ##
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
15 ## Octave is distributed in the hope that it will be useful, but
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
16 ## WITHOUT ANY WARRANTY; without even the implied warranty of
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
17 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
18 ## GNU General Public License for more details.
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
19 ##
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
20 ## You should have received a copy of the GNU General Public License
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
21 ## along with Octave; see the file COPYING. If not, see
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
22 ## <https://www.gnu.org/licenses/>.
27923
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
23 ##
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
24 ########################################################################
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
25
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
26 ## -*- texinfo -*-
30875
5d3faba0342e doc: Ensure documentation lists output argument when it exists for all m-files.
Rik <rik@octave.org>
parents: 30564
diff changeset
27 ## @deftypefn {} {@var{x} =} gammaincinv (@var{y}, @var{a})
5d3faba0342e doc: Ensure documentation lists output argument when it exists for all m-files.
Rik <rik@octave.org>
parents: 30564
diff changeset
28 ## @deftypefnx {} {@var{x} =} gammaincinv (@var{y}, @var{a}, @var{tail})
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
29 ## Compute the inverse of the normalized incomplete gamma function.
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
30 ##
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
31 ## The normalized incomplete gamma function is defined as
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
32 ## @tex
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
33 ## $$
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
34 ## \gamma (x, a) = {1 \over {\Gamma (a)}}\displaystyle{\int_0^x t^{a-1} e^{-t} dt}
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
35 ## $$
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
36 ## @end tex
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
37 ## @ifnottex
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
38 ##
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
39 ## @example
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
40 ## @group
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
41 ## x
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
42 ## 1 /
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
43 ## gammainc (x, a) = --------- | exp (-t) t^(a-1) dt
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
44 ## gamma (a) /
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
45 ## t=0
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
46 ## @end group
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
47 ## @end example
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
48 ##
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
49 ## @end ifnottex
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
50 ##
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
51 ## and @code{gammaincinv (gammainc (@var{x}, @var{a}), @var{a}) = @var{x}}
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
52 ## for each non-negative value of @var{x}. If @var{a} is scalar then
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
53 ## @code{gammaincinv (@var{y}, @var{a})} is returned for each element of
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
54 ## @var{y} and vice versa.
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
55 ##
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
56 ## If neither @var{y} nor @var{a} is scalar then the sizes of @var{y} and
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
57 ## @var{a} must agree, and @code{gammaincinv} is applied element-by-element.
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
58 ## The variable @var{y} must be in the interval @math{[0,1]} while @var{a} must
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
59 ## be real and positive.
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
60 ##
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
61 ## By default, @var{tail} is @qcode{"lower"} and the inverse of the incomplete
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
62 ## gamma function integrated from 0 to @var{x} is computed. If @var{tail} is
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
63 ## @qcode{"upper"}, then the complementary function integrated from @var{x} to
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
64 ## infinity is inverted.
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
65 ##
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
66 ## The function is computed with Newton's method by solving
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
67 ## @tex
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
68 ## $$
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
69 ## y - \gamma (x, a) = 0
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
70 ## $$
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
71 ## @end tex
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
72 ## @ifnottex
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
73 ##
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
74 ## @example
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
75 ## @var{y} - gammainc (@var{x}, @var{a}) = 0
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
76 ## @end example
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
77 ##
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
78 ## @end ifnottex
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
79 ##
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
80 ## Reference: @nospell{A. Gil, J. Segura, and N. M. Temme}, @cite{Efficient and
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
81 ## accurate algorithms for the computation and inversion of the incomplete
27800
5a6a19a4e3da doc: Use Texinfo non-sentence ending periods in citations.
Rik <rik@octave.org>
parents: 27411
diff changeset
82 ## gamma function ratios}, @nospell{SIAM J. Sci.@: Computing}, pp.@:
5a6a19a4e3da doc: Use Texinfo non-sentence ending periods in citations.
Rik <rik@octave.org>
parents: 27411
diff changeset
83 ## A2965--A2981, Vol 34, 2012.
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
84 ##
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
85 ## @seealso{gammainc, gamma, gammaln}
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
86 ## @end deftypefn
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
87
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
88 function x = gammaincinv (y, a, tail = "lower")
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
89
28789
28de41192f3c Eliminate unneeded verification of nargin, nargout in m-files.
Rik <rik@octave.org>
parents: 27984
diff changeset
90 if (nargin < 2)
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
91 print_usage ();
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
92 endif
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
93
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
94 [err, y, a] = common_size (y, a);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
95 if (err > 0)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
96 error ("gammaincinv: Y and A must be of common size or scalars");
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
97 endif
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
98
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
99 if (iscomplex (y) || iscomplex (a))
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
100 error ("gammaincinv: all inputs must be real");
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
101 endif
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
102
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
103 ## Remember original shape of data, but convert to column vector for calcs.
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
104 orig_sz = size (y);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
105 y = y(:);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
106 a = a(:);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
107
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
108 if (any ((y < 0) | (y > 1)))
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
109 error ("gammaincinv: Y must be in the range [0, 1]");
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
110 endif
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
111
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
112 if (any (a <= 0))
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
113 error ("gammaincinv: A must be strictly positive");
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
114 endif
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
115
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
116 ## If any of the arguments is single then the output should be as well.
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
117 if (strcmp (class (y), "single") || strcmp (class (a), "single"))
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
118 y = single (y);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
119 a = single (a);
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
120 endif
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
121
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
122 ## Convert to floating point if necessary
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
123 if (isinteger (y))
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
124 y = double (y);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
125 endif
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
126 if (isinteger (a))
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
127 a = double (a);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
128 endif
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
129
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
130 ## Initialize output array
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
131 x = zeros (size (y), class (y));
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
132
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
133 maxit = 20;
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
134 tol = eps (class (y));
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
135
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
136 ## Special cases, a = 1 or y = 0, 1.
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
137
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
138 if (strcmpi (tail, "lower"))
27411
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
139 x(a == 1) = - log1p (- y(a == 1));
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
140 x(y == 0) = 0;
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
141 x(y == 1) = Inf;
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
142 p = y;
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
143 q = 1 - p;
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
144 elseif (strcmpi (tail, "upper"))
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
145 x(a == 1) = - log (y(a == 1));
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
146 x(y == 0) = Inf;
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
147 x(y == 1) = 0;
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
148 q = y;
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
149 p = 1 - q;
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
150 else
28928
ae7ce8358953 maint: Add semicolon to end of all warning() and error() invocations.
Rik <rik@octave.org>
parents: 28896
diff changeset
151 error ("gammaincinv: invalid value for TAIL");
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
152 endif
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
153
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
154 todo = (a != 1) & (y != 0) & (y != 1);
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
155
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
156 ## Case 1: p small.
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
157
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
158 i_flag_1 = todo & (p < ((0.2 * (1 + a)) .^ a) ./ gamma (1 + a));
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
159
27411
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
160 if (any (i_flag_1))
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
161 aa = a(i_flag_1);
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
162 pp = p(i_flag_1);
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
163
27411
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
164 ## Initial guess.
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
165
27411
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
166 r = (pp .* gamma (1 + aa)) .^ (1 ./ aa);
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
167
27411
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
168 c2 = 1 ./ (aa + 1);
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
169 c3 = (3 * aa + 5) ./ (2 * (aa + 1) .^2 .* (aa + 2));
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
170 c4 = (8 * aa .^ 2 + 33 * aa + 31) ./ (3 * (aa + 1) .^ 3 .* (aa + 2) .* ...
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
171 (aa + 3));
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
172 c5 = (125 * aa .^ 4 + 1179 * aa .^ 3 + 3971 * aa.^2 + 5661 * aa + 2888) ...
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
173 ./ (24 * (1 + aa) .^4 .* (aa + 2) .^ 2 .* (aa + 3) .* (aa + 4));
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
174
27411
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
175 ## FIXME: Would polyval() be better here for more accuracy?
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
176 x0 = r + c2 .* r .^ 2 + c3 .* r .^ 3 + c4 .* r .^4 + c5 .* r .^ 5;
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
177
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
178 ## For this case we invert the lower version.
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
179
27411
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
180 F = @(p, a, x) p - gammainc (x, a, "lower");
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
181 JF = @(a, x) - exp (- gammaln (a) - x + (a - 1) .* log (x));
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
182 x(i_flag_1) = newton_method (F, JF, pp, aa, x0, tol, maxit);
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
183 endif
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
184
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
185 todo(i_flag_1) = false;
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
186
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
187 ## Case 2: q small.
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
188
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
189 i_flag_2 = (q < exp (- 0.5 * a) ./ gamma (1 + a)) & (a > 0) & (a < 10);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
190 i_flag_2 &= todo;
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
191
27411
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
192 if (any (i_flag_2))
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
193 aa = a(i_flag_2);
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
194 qq = q(i_flag_2);
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
195
27411
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
196 ## Initial guess.
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
197
27411
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
198 x0 = (-log (qq) - gammaln (aa));
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
199
27411
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
200 ## For this case, we invert the upper version.
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
201
27411
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
202 F = @(q, a, x) q - gammainc (x, a, "upper");
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
203 JF = @(a, x) exp (- gammaln (a) - x) .* x .^ (a - 1);
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
204 x(i_flag_2) = newton_method (F, JF, qq, aa, x0, tol, maxit);
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
205 endif
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
206
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
207 todo(i_flag_2) = false;
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
208
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
209 ## Case 3: a small.
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
210
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
211 i_flag_3 = todo & ((a > 0) & (a < 1));
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
212
27411
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
213 if (any (i_flag_3))
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
214 aa = a(i_flag_3);
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
215 pp = p(i_flag_3);
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
216
27411
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
217 ## Initial guess
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
218
27411
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
219 xl = (pp .* gamma (aa + 1)) .^ (1 ./ aa);
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
220 x0 = xl;
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
221
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
222 ## For this case, we invert the lower version.
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
223
27411
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
224 F = @(p, a, x) p - gammainc (x, a, "lower");
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
225 JF = @(a, x) - exp (-gammaln (a) - x) .* x .^ (a - 1);
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
226 x(i_flag_3) = newton_method (F, JF, pp, aa, x0, tol, maxit);
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
227 endif
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
228
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
229 todo(i_flag_3) = false;
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
230
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
231 ## Case 4: a large.
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
232
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
233 i_flag_4 = todo;
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
234
27411
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
235 if (any (i_flag_4))
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
236 aa = a(i_flag_4);
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
237 qq = q(i_flag_4);
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
238
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
239 ## Initial guess
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
240
27411
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
241 d = 1 ./ (9 * aa);
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
242 t = 1 - d + sqrt (2) * erfcinv (2 * qq) .* sqrt (d);
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
243 x0 = aa .* (t .^ 3);
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
244
27411
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
245 ## For this case, we invert the upper version.
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
246
27411
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
247 F = @(q, a, x) q - gammainc (x, a, "upper");
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
248 JF = @(a, x) exp (- gammaln (a) - x + (a - 1) .* log (x));
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
249 x(i_flag_4) = newton_method (F, JF, qq, aa, x0, tol, maxit);
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
250 endif
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
251
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
252 ## Restore original shape
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
253 x = reshape (x, orig_sz);
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
254
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
255 endfunction
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
256
27984
b09432b20a84 maint: Remove special cases of old version control keywords in code base.
Rik <rik@octave.org>
parents: 27978
diff changeset
257 ## subfunction: Newton's Method
31548
c8ad083a5802 maint: Clean up m-files before Octave 8.1 release.
Rik <rik@octave.org>
parents: 30875
diff changeset
258 function x = newton_method (F, JF, y, a, x0, tol, maxit)
28945
6e460773bdda maint: Use newlines after "function" and before "endfunction" for clarity.
Rik <rik@octave.org>
parents: 28944
diff changeset
259
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
260 l = numel (y);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
261 res = -F (y, a, x0) ./ JF (a, x0);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
262 todo = (abs (res) >= tol * abs (x0));
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
263 x = x0;
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
264 it = 0;
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
265 while (any (todo) && (it++ < maxit))
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
266 x(todo) += res(todo);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
267 res(todo) = -F (y(todo), a(todo), x(todo)) ./ JF (a(todo), x(todo));
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
268 todo = (abs (res) >= tol * abs (x));
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
269 endwhile
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
270 x += res;
28945
6e460773bdda maint: Use newlines after "function" and before "endfunction" for clarity.
Rik <rik@octave.org>
parents: 28944
diff changeset
271
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
272 endfunction
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
273
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
274
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
275 %!test
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
276 %! x = [1e-10, 1e-09, 1e-08, 1e-07];
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
277 %! a = [2, 3, 4];
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
278 %! [x, a] = ndgrid (x, a);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
279 %! xx = gammainc (gammaincinv (x, a), a);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
280 %! assert (xx, x, -3e-14);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
281
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
282 %!test
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
283 %! x = [1e-10, 1e-09, 1e-08, 1e-07];
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
284 %! a = [2, 3, 4];
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
285 %! [x, a] = ndgrid (x, a);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
286 %! xx = gammainc (gammaincinv (x, a, "upper"), a, "upper");
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
287 %! assert (xx, x, -3e-14);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
288
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
289 %!test
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
290 %! x = linspace (0, 1)';
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
291 %! a = [linspace(0.1, 1, 10), 2:5];
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
292 %! [x, a] = ndgrid (x, a);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
293 %! xx = gammainc (gammaincinv (x, a), a);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
294 %! assert (xx, x, -1e-13);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
295
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
296 %!test
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
297 %! x = linspace (0, 1)';
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
298 %! a = [linspace(0.1, 1, 10), 2:5];
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
299 %! [x, a] = ndgrid (x, a);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
300 %! xx = gammainc (gammaincinv (x, a, "upper"), a, "upper");
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
301 %! assert (xx, x, -1e-13);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
302
27411
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
303 %!test <*56453>
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
304 %! assert (gammaincinv (1e-15, 1) * 2, 2e-15, -1e-15);
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
305 %! assert (gammaincinv (1e-16, 1) * 2, 2e-16, -1e-15);
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 26376
diff changeset
306
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
307 ## Test the conservation of the input class
24905
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
308 %!assert (class (gammaincinv (0.5, 1)), "double")
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
309 %!assert (class (gammaincinv (single (0.5), 1)), "single")
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
310 %!assert (class (gammaincinv (0.5, single (1))), "single")
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
311 %!assert (class (gammaincinv (int8 (0), 1)), "double")
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
312 %!assert (class (gammaincinv (0.5, int8 (1))), "double")
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
313 %!assert (class (gammaincinv (int8 (0), single (1))), "single")
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
314 %!assert (class (gammaincinv (single (0.5), int8 (1))), "single")
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
315
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036)
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff changeset
316 ## Test input validation
28896
90fea9cc9caa test: Add expected error message <Invalid call> to BIST tests for nargin.
Rik <rik@octave.org>
parents: 28789
diff changeset
317 %!error <Invalid call> gammaincinv ()
90fea9cc9caa test: Add expected error message <Invalid call> to BIST tests for nargin.
Rik <rik@octave.org>
parents: 28789
diff changeset
318 %!error <Invalid call> gammaincinv (1)
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
319 %!error <must be of common size or scalars>
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
320 %! gammaincinv (ones (2,2), ones (1,2), 1);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
321 %!error <all inputs must be real> gammaincinv (0.5i, 1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
322 %!error <all inputs must be real> gammaincinv (0, 1i)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
323 %!error <Y must be in the range \[0, 1\]> gammaincinv (-0.1,1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
324 %!error <Y must be in the range \[0, 1\]> gammaincinv (1.1,1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
325 %!error <Y must be in the range \[0, 1\]>
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
326 %! y = ones (1, 1, 2);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
327 %! y(1,1,2) = -1;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
328 %! gammaincinv (y,1);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
329 %!error <A must be strictly positive> gammaincinv (0.5, 0)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
330 %!error <A must be strictly positive>
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
331 %! a = ones (1, 1, 2);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
332 %! a(1,1,2) = 0;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
333 %! gammaincinv (1,a,1);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24905
diff changeset
334 %!error <invalid value for TAIL> gammaincinv (1,2, "foobar")