octave: scripts/specfun/gammaincinv.m annotate

annotate scripts/specfun/gammaincinv.m @ 27984:b09432b20a84

maint: Remove special cases of old version control keywords in code base. * common_size.m, cumtrapz.m, interp1.m, polyarea.m, trapz.m, inpolygon.m, summer.m, javamem.m, duplication_matrix.m, subspace.m, starting_stepsize.m, text.m, errorbar.m, loglogerr.m, peaks.m, __errplot__.m, semilogxerr.m, semilogyerr.m, stemleaf.m, streamtube.m, __pltopt__.m, printd.m, pchip.m, arch_fit.m, arch_rnd.m, arch_test.m, arma_rnd.m, autoreg_matrix.m, bartlett.m, blackman.m, diffpara.m, durbinlevinson.m, fractdiff.m, hamming.m, hanning.m, hurst.m, ifftshift.m, periodogram.m, rectangle_lw.m, rectangle_sw.m, triangle_lw.m, triangle_sw.m, sinetone.m, sinewave.m, spectral_adf.m, spectral_xdf.m, spencer.m, stft.m, synthesis.m, yulewalker.m, bicg.m, bicgstab.m, cgs.m, gmres.m, pcr.m, sprand.m, tfqmr.m, betaincinv.m, betaln.m, gammaincinv.m, hadamard.m, center.m, cov.m, discrete_inv.m, discrete_pdf.m, discrete_rnd.m, empirical_cdf.m, empirical_inv.m, empirical_pdf.m, empirical_rnd.m, iqr.m, kendall.m, mean.m, meansq.m, moment.m, prctile.m, quantile.m, range.m, run_count.m, spearman.m, statistics.m, var.m, zscore.m, datenum.m, datevec.m: Remove special cases of old version control keywords in code base.

author	Rik <rik@octave.org>
date	Tue, 21 Jan 2020 14:04:31 -0800
parents	a4268efb7334
children	28de41192f3c 0a5b15007766

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 ##
27919 1891570abac8 update Octave Project Developers copyright for the new year John W. Eaton <jwe@octave.org> parents: 27918 diff changeset	3 ## Copyright (C) 2017-2020 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
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	12 ## 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	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
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	90 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	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
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 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);
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	259 l = numel (y);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	260 res = -F (y, a, x0) ./ JF (a, x0);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	261 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	262 x = x0;
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	263 it = 0;
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	264 while (any (todo) && (it++ < maxit))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	265 x(todo) += res(todo);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	266 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	267 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	268 endwhile
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	269 x += res;
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	270 endfunction
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	271
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 = [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	275 %! 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	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, -3e-14);
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 = [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	282 %! 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	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, -3e-14);
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	286
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	287 %!test
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	288 %! 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	289 %! 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	290 %! [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	291 %! 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	292 %! 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	293
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	294 %!test
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	295 %! 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	296 %! 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	297 %! [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	298 %! 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	299 %! 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	300
27411 9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453). Michele Ginesi <michele.ginesi@gmail.com> parents: 26376 diff changeset	301 %!test <*56453>
9a946e34553c gammaincinv.m: Improve accuracy for gammaincinv (bug #56453). Michele Ginesi <michele.ginesi@gmail.com> parents: 26376 diff changeset	302 %! 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	303 %! 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	304
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	305 ## 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	306 %!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	307 %!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	308 %!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	309 %!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	310 %!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	311 %!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	312 %!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	313
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	314 ## Test input validation
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	315 %!error gammaincinv ()
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	316 %!error gammaincinv (1)
662faf9de127 Added the inverse of the incomplete gamma function (see bug #48036) Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	317 %!error gammaincinv (1, 2, 3, 4)
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	318 %!error <must be of common size or scalars>
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	319 %! gammaincinv (ones (2,2), ones (1,2), 1);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	320 %!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	321 %!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	322 %!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	323 %!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	324 %!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	325 %! y = ones (1, 1, 2);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	326 %! y(1,1,2) = -1;
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	327 %! gammaincinv (y,1);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	328 %!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	329 %!error <A must be strictly positive>
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	330 %! a = ones (1, 1, 2);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	331 %! a(1,1,2) = 0;
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	332 %! gammaincinv (1,a,1);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24905 diff changeset	333 %!error <invalid value for TAIL> gammaincinv (1,2, "foobar")