octave: scripts/specfun/gammaincinv.m annotate

annotate scripts/specfun/gammaincinv.m @ 26376:00f796120a6d stable

maint: Update copyright dates in all source files.

author	John W. Eaton <jwe@octave.org>
date	Wed, 02 Jan 2019 16:32:43 -0500
parents	6652d3823428
children	9a946e34553c

rev	line source
26376 00f796120a6d maint: Update copyright dates in all source files. John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	1 ## Copyright (C) 2017-2019 Michele Ginesi
24905 662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	2 ##
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	3 ## 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	4 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	5 ## 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	6 ## under the terms of the GNU General Public License as published by
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	7 ## the Free Software Foundation; either version 3 of the License, or
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	8 ## (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	9 ##
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	10 ## 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	11 ## 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	12 ## 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	13 ## 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	14 ##
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	15 ## 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	16 ## 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	17 ## <https://www.gnu.org/licenses/>.
24905 662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	18
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	19 ## Author: Michele Ginesi <michele.ginesi@gmail.com>
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	20
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	21 ## -- texinfo --
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	22 ## @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	23 ## @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	24 ## 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	25 ##
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	26 ## 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	27 ## @tex
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	28 ## $$
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	29 ## \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	30 ## $$
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	31 ## @end tex
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	32 ## @ifnottex
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 ## @example
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	35 ## @group
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	36 ## x
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	37 ## 1 /
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	38 ## 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	39 ## gamma (a) /
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	40 ## t=0
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	41 ## @end group
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	42 ## @end example
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	43 ##
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	44 ## @end ifnottex
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	45 ##
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	46 ## 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	47 ## 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	48 ## @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	49 ## @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	50 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	51 ## 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	52 ## @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	53 ## 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	54 ## 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	55 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	56 ## 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	57 ## 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	58 ## @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	59 ## 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	60 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	61 ## 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	62 ## @tex
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	63 ## $$
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	64 ## 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	65 ## $$
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	66 ## @end tex
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	67 ## @ifnottex
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 ## @example
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	70 ## @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	71 ## @end example
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	72 ##
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	73 ## @end ifnottex
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	74 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	75 ## 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	76 ## accurate algorithms for the computation and inversion of the incomplete
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	77 ## gamma function ratios}, @nospell{SIAM J. Sci. Computing}, pp. A2965--A2981,
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	78 ## 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	79 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	80 ## @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	81 ## @end deftypefn
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	82
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	83 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	84
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	85 if (nargin < 2 \|\| nargin > 3)
24905 662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	86 print_usage ();
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	87 endif
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	88
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	89 [err, y, a] = common_size (y, a);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	90 if (err > 0)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	91 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	92 endif
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	93
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	94 if (iscomplex (y) \|\| iscomplex (a))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	95 error ("gammaincinv: all inputs must be real");
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	96 endif
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	97
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	98 ## 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	99 orig_sz = size (y);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	100 y = y(:);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	101 a = a(:);
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 if (any ((y < 0) \| (y > 1)))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	104 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	105 endif
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	106
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	107 if (any (a <= 0))
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	108 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	109 endif
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	110
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	111 ## 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	112 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	113 y = single (y);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	114 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	115 endif
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	116
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	117 ## 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	118 if (isinteger (y))
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	119 y = double (y);
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 if (isinteger (a))
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	122 a = double (a);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	123 endif
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	124
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	125 ## Initialize output array
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	126 x = zeros (size (y), class (y));
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	127
24905 662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	128 maxit = 20;
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	129 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	130
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	131 ## 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	132
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	133 if (strcmpi (tail, "lower"))
24905 662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	134 x(a == 1) = - log (1 - y(a == 1));
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	135 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	136 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	137 p = y;
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	138 q = 1 - p;
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	139 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	140 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	141 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	142 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	143 q = y;
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	144 p = 1 - q;
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	145 else
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	146 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	147 endif
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	148
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	149 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	150
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	151 ## 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	152
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	153 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	154
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	155 aa = a(i_flag_1);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	156 pp = p(i_flag_1);
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 ## Initial guess.
24905 662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	159
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	160 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	161
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	162 c2 = 1 ./ (aa + 1);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	163 c3 = (3 * aa + 5) ./ (2 * (aa + 1) .^2 .* (aa + 2));
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	164 c4 = (8 * aa .^ 2 + 33 * aa + 31) ./ (3 * (aa + 1) .^ 3 .* (aa + 2) .* ...
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	165 (aa + 3));
24905 662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	166 c5 = (125 * aa .^ 4 + 1179 * aa .^ 3 + 3971 * aa.^2 + 5661 * aa + 2888) ...
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	167 ./ (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	168
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	169 ## FIXME: Would polyval() be better here for more accuracy?
24905 662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	170 x0 = r + c2 .* r .^ 2 + c3 .* r .^ 3 + c4 .* r .^4 + c5 .* r .^ 5;
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	171
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	172 ## 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	173
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	174 F = @(p, a, x) p - gammainc (x, a, "lower");
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	175 JF = @(a, x) - exp (-gammaln (a) - x + (a - 1) .* log (x));
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	176 x(i_flag_1) = newton_method (F, JF, pp, aa, x0, tol, maxit);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	177
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	178 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	179
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	180 ## 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	181
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	182 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	183 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	184
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	185 aa = a(i_flag_2);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	186 qq = q(i_flag_2);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	187
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	188 ## Initial guess.
24905 662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	189
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	190 x0 = (-log (qq) - gammaln (aa));
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	191
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	192 ## 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	193
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	194 F = @(q, a, x) q - gammainc (x, a, "upper");
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	195 JF = @(a, x) exp (- gammaln (a) - x) .* x .^ (a - 1);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	196 x(i_flag_2) = newton_method (F, JF, qq, aa, x0, tol, maxit);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	197
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	198 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	199
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	200 ## 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	201
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	202 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	203
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	204 aa = a(i_flag_3);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	205 pp = p(i_flag_3);
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 ## Initial guess
24905 662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	208
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	209 xl = (pp .* gamma (aa + 1)) .^ (1 ./ aa);
24905 662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	210 x0 = xl;
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	211
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	212 ## 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	213
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	214 F = @(p, a, x) p - gammainc (x, a, "lower");
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	215 JF = @(a, x) - exp (-gammaln (a) - x) .* x .^ (a - 1);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	216 x(i_flag_3) = newton_method (F, JF, pp, aa, x0, tol, maxit);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	217
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	218 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	219
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	220 ## 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	221
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	222 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	223 aa = a(i_flag_4);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	224 qq = q(i_flag_4);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	225
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	226 ## Initial guess
24905 662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	227
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	228 d = 1 ./ (9 * aa);
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	229 t = 1 - d + sqrt (2) * erfcinv (2 * qq) .* sqrt (d);
24905 662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	230 x0 = aa .* (t .^ 3);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	231
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	232 ## 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	233
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	234 F = @(q, a, x) q - gammainc (x, a, "upper");
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	235 JF = @(a, x) exp (- gammaln (a) - x + (a - 1) .* log (x));
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	236 x(i_flag_4) = newton_method (F, JF, qq, aa, x0, tol, maxit);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	237
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	238 ## Restore original shape
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	239 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	240
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	241 endfunction
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	242
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	243 ## 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	244 function x = newton_method (F, JF, y, a, x0, tol, maxit);
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	245 l = numel (y);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	246 res = -F (y, a, x0) ./ JF (a, x0);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	247 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	248 x = x0;
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	249 it = 0;
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	250 while (any (todo) && (it++ < maxit))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	251 x(todo) += res(todo);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	252 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	253 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	254 endwhile
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	255 x += res;
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	256 endfunction
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	257
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	258
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	259 %!test
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	260 %! 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	261 %! 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	262 %! [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	263 %! 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	264 %! 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	265
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	266 %!test
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	267 %! 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	268 %! 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	269 %! [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	270 %! 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	271 %! 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	272
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	273 %!test
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	274 %! 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	275 %! 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	276 %! [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	277 %! 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	278 %! 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	279
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	280 %!test
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	281 %! 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	282 %! 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	283 %! [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	284 %! 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	285 %! 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	286
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	287 ## 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	288 %!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	289 %!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	290 %!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	291 %!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	292 %!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	293 %!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	294 %!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	295
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	296 ## Test input validation
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	297 %!error gammaincinv ()
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	298 %!error gammaincinv (1)
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	299 %!error gammaincinv (1, 2, 3, 4)
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	300 %!error <must be of common size or scalars>
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	301 %! gammaincinv (ones (2,2), ones (1,2), 1);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	302 %!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	303 %!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	304 %!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	305 %!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	306 %!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	307 %! y = ones (1, 1, 2);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	308 %! y(1,1,2) = -1;
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	309 %! gammaincinv (y,1);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	310 %!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	311 %!error <A must be strictly positive>
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	312 %! a = ones (1, 1, 2);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	313 %! a(1,1,2) = 0;
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	314 %! gammaincinv (1,a,1);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	315 %!error <invalid value for TAIL> gammaincinv (1,2, "foobar")