octave: scripts/specfun/gammaincinv.m annotate

annotate scripts/specfun/gammaincinv.m @ 30564:796f54d4ddbf stable

update Octave Project Developers copyright for the new year In files that have the "Octave Project Developers" copyright notice, update for 2021. In all .txi and .texi files except gpl.txi and gpl.texi in the doc/liboctave and doc/interpreter directories, change the copyright to "Octave Project Developers", the same as used for other source files. Update copyright notices for 2022 (not done since 2019). For gpl.txi and gpl.texi, change the copyright notice to be "Free Software Foundation, Inc." and leave the date at 2007 only because this file only contains the text of the GPL, not anything created by the Octave Project Developers. Add Paul Thomas to contributors.in.

author	John W. Eaton <jwe@octave.org>
date	Tue, 28 Dec 2021 18:22:40 -0500
parents	7854d5752dd2
children	5d3faba0342e

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 --
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	27 ## @deftypefn {} {} gammaincinv (@var{y}, @var{a})
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	28 ## @deftypefnx {} {} gammaincinv (@var{y}, @var{a}, @var{tail})
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
24905 662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: 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")