octave: scripts/specfun/betaincinv.m annotate

annotate scripts/specfun/betaincinv.m @ 32069:217a35f2a589 stable

betaincinv.m: Slightly increase test tolerance (for macOS 13). * scripts/specfun/betaincinv.m: Slightly increase test tolerance to pass on macOS 13.

author	Markus Mützel <markus.muetzel@gmx.de>
date	Tue, 02 May 2023 11:32:58 +0200
parents	597f3ee61a48
children	2e484f9f1f18

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 ##
31706 597f3ee61a48 update Octave Project Developers copyright for the new year John W. Eaton <jwe@octave.org> parents: 31548 diff changeset	3 ## Copyright (C) 2017-2023 The Octave Project Developers
27918 b442ec6dda5c use centralized file for copyright info for individual contributors John W. Eaton <jwe@octave.org> parents: 26376 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/>.
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	7 ##
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	8 ## This file is part of Octave.
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	9 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	10 ## Octave is free software: you can redistribute it and/or modify it
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	11 ## under the terms of the GNU General Public License as published by
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	12 ## the Free Software Foundation, either version 3 of the License, or
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	13 ## (at your option) any later version.
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	14 ##
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	15 ## Octave is distributed in the hope that it will be useful, but
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	16 ## WITHOUT ANY WARRANTY; without even the implied warranty of
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	17 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	18 ## GNU General Public License for more details.
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	19 ##
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	20 ## You should have received a copy of the GNU General Public License
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and 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: 24907 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 ########################################################################
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	25
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	26 ## -- texinfo --
30875 5d3faba0342e doc: Ensure documentation lists output argument when it exists for all m-files. Rik <rik@octave.org> parents: 30564 diff changeset	27 ## @deftypefn {} {@var{x} =} betaincinv (@var{y}, @var{a}, @var{b})
5d3faba0342e doc: Ensure documentation lists output argument when it exists for all m-files. Rik <rik@octave.org> parents: 30564 diff changeset	28 ## @deftypefnx {} {@var{x} =} betaincinv (@var{y}, @var{a}, @var{b}, "lower")
5d3faba0342e doc: Ensure documentation lists output argument when it exists for all m-files. Rik <rik@octave.org> parents: 30564 diff changeset	29 ## @deftypefnx {} {@var{x} =} betaincinv (@var{y}, @var{a}, @var{b}, "upper")
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	30 ## Compute the inverse of the normalized incomplete beta function.
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	31 ##
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	32 ## The normalized incomplete beta function is defined as
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	33 ## @tex
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	34 ## $$
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	35 ## I_x (a, b) = {1 \over {B(a,b)}} \displaystyle{\int_0^x t^{a-1} (1-t)^{b-1} dt}
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	36 ## $$
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	37 ## @end tex
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	38 ## @ifnottex
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	39 ##
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	40 ## @example
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	41 ## @group
29664 2a1f57067fbf betainc.m, betaincinv.m: Correct non-Tex definition of beta incomplete integral. Rik <rik@octave.org> parents: 29656 diff changeset	42 ## x
2a1f57067fbf betainc.m, betaincinv.m: Correct non-Tex definition of beta incomplete integral. Rik <rik@octave.org> parents: 29656 diff changeset	43 ## /
2a1f57067fbf betainc.m, betaincinv.m: Correct non-Tex definition of beta incomplete integral. Rik <rik@octave.org> parents: 29656 diff changeset	44 ## 1 \|
2a1f57067fbf betainc.m, betaincinv.m: Correct non-Tex definition of beta incomplete integral. Rik <rik@octave.org> parents: 29656 diff changeset	45 ## I_x (a, b) = ---------- \| t^(a-1) (1-t)^(b-1) dt
2a1f57067fbf betainc.m, betaincinv.m: Correct non-Tex definition of beta incomplete integral. Rik <rik@octave.org> parents: 29656 diff changeset	46 ## beta (a,b) \|
2a1f57067fbf betainc.m, betaincinv.m: Correct non-Tex definition of beta incomplete integral. Rik <rik@octave.org> parents: 29656 diff changeset	47 ## /
2a1f57067fbf betainc.m, betaincinv.m: Correct non-Tex definition of beta incomplete integral. Rik <rik@octave.org> parents: 29656 diff changeset	48 ## 0
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	49 ## @end group
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	50 ## @end example
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	51 ##
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	52 ## @end ifnottex
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	53 ##
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	54 ## If two inputs are scalar, then @code{betaincinv (@var{y}, @var{a}, @var{b})}
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	55 ## is returned for each of the other inputs.
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	56 ##
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	57 ## If two or more inputs are not scalar, the sizes of them must agree, and
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	58 ## @code{betaincinv} is applied element-by-element.
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	59 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	60 ## The variable @var{y} must be in the interval [0,1], while @var{a} and
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	61 ## @var{b} must be real and strictly positive.
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	62 ##
25579 07c2c42f457e doc: Miscellaneous documentation fixes all over the manual (bug #54288). Rik <rik@octave.org> parents: 25054 diff changeset	63 ## By default, @var{tail} is @qcode{"lower"} and the inverse of the incomplete
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	64 ## beta 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: 24907 diff changeset	65 ## @qcode{"upper"} then the complementary function integrated from @var{x} to 1
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	66 ## is inverted.
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	67 ##
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	68 ## The function is computed by standard Newton's method, by solving
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	69 ## @tex
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	70 ## $$
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	71 ## y - I_x (a, b) = 0
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	72 ## $$
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	73 ## @end tex
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	74 ## @ifnottex
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	75 ##
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	76 ## @example
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	77 ## @var{y} - betainc (@var{x}, @var{a}, @var{b}) = 0
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	78 ## @end example
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	79 ##
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	80 ## @end ifnottex
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	81 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	82 ## @seealso{betainc, beta, betaln}
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	83 ## @end deftypefn
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	84
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	85 function x = betaincinv (y, a, b, tail = "lower")
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	86
29656 7f5bd197fea6 maint: use std::size_t in more instances (bug #60471) Rik <rik@octave.org> parents: 29616 diff changeset	87 if (nargin < 3)
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	88 print_usage ();
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	89 endif
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	90
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	91 [err, y, a, b] = common_size (y, a, b);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	92 if (err > 0)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	93 error ("betaincinv: Y, A, and B must be of common size or scalars");
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	94 endif
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	95
29656 7f5bd197fea6 maint: use std::size_t in more instances (bug #60471) Rik <rik@octave.org> parents: 29616 diff changeset	96 if (! (isfloat (y) && isfloat (a) && isfloat (b)
7f5bd197fea6 maint: use std::size_t in more instances (bug #60471) Rik <rik@octave.org> parents: 29616 diff changeset	97 && isreal (y) && isreal (a) && isreal (b)))
7f5bd197fea6 maint: use std::size_t in more instances (bug #60471) Rik <rik@octave.org> parents: 29616 diff changeset	98 error ("betaincinv: Y, A, and B must be real, floating point values");
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	99 endif
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	100
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	101 ## 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: 24907 diff changeset	102 orig_sz = size (y);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	103 y = y(:);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	104 a = a(:);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	105 b = b(:);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	106
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	107 if (any ((y < 0) \| (y > 1)))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	108 error ("betaincinv: Y must be in the range [0, 1]");
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	109 endif
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	110
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	111 if (any (a <= 0))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	112 error ("betaincinv: A must be strictly positive");
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	113 endif
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	114
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	115 if (any (b <= 0))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	116 error ("betaincinv: B must be strictly positive");
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	117 endif
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	118
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	119 ## If any of the arguments is single then the output should be as well.
29656 7f5bd197fea6 maint: use std::size_t in more instances (bug #60471) Rik <rik@octave.org> parents: 29616 diff changeset	120 if (isa (y, "single") \|\| isa (a, "single") \|\| isa (b, "single"))
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	121 y = single (y);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	122 a = single (a);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	123 b = single (b);
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	124 endif
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	125
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	126 if (strcmpi (tail, "lower"))
30017 26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	127 ys = y;
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	128 elseif (strcmpi (tail, "upper"))
30017 26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	129 ys = 1 - y; # only for computation of initial points, no loss of accuracy
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	130 else
28928 ae7ce8358953 maint: Add semicolon to end of all warning() and error() invocations. Rik <rik@octave.org> parents: 28912 diff changeset	131 error ("betaincinv: invalid value for TAIL");
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	132 endif
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	133
30017 26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	134 ## Choose starting point for Newton's Method to guarantee convergence.
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	135 ## If (a-1)*(b-1) > 0, F has a point of inflection at x = (a-1)/(a+b-2).
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	136 ## In this case, it is convex on (0,x) and concave on (x,1) if a>1; otherwise
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	137 ## it is the other way round. If (a-1)*(b-1) <= 0, there is no point of
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	138 ## inflection, and it is everywhere convex for a>1 and concave otherwise.
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	139 ## We thus choose our starting x for the Newton iterations so that we stay
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	140 ## within a region of constant sign of curvature and on the correct side of
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	141 ## the eventual solution, guaranteeing convergence. Curvatures above are to
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	142 ## be understood under the condition tail=="lower".
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	143
30017 26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	144 ## Initialize output array
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	145 x = x_i = y_i = zeros (size (y), class (y));
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	146
30017 26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	147 ## Have point of inflection
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	148 idx = find ((a - 1) .* (b - 1) > 0);
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	149 if (! isempty (idx))
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	150 x_i(idx) = (a(idx) - 1) ./ (a(idx) + b(idx) - 2);
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	151 y_i(idx) = betainc (x_i(idx), a(idx), b(idx));
30379 363fb10055df maint: Style check m-files ahead of 7.1 release. Rik <rik@octave.org> parents: 30029 diff changeset	152 endif
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	153
30017 26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	154 ## Converge outwards
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	155 tmpidx = find (a(idx) > 1);
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	156 if (! isempty (tmpidx))
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	157 x(idx(tmpidx)) = x_i(idx(tmpidx));
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	158 endif
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	159 ## Converge inwards
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	160 ## To the left of inflection point
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	161 tmpidx = idx(find ((a(idx) <= 1) & (y_i(idx) >= ys(idx))));
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	162 if (! isempty (tmpidx))
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	163 x(tmpidx) = (ys(tmpidx) ./ y_i(tmpidx)).^(1 ./ a(tmpidx)) .* x_i(tmpidx);
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	164 endif
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	165 ## To the right of inflection point
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	166 tmpidx = idx(find ((a(idx) <= 1) & (y_i(idx) < ys(idx))));
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	167 if (! isempty (tmpidx))
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	168 x(tmpidx) = 1 - ...
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	169 ((1 - ys(tmpidx)) ./ (1 - y_i(tmpidx))).^(1 ./ b(tmpidx)) ...
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	170 .* (1 - x_i(tmpidx));
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	171 endif
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	172
30017 26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	173 ## Have no point of inflection
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	174 idx = find ((a - 1) .* (b - 1) <= 0);
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	175
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	176 ## Negative curvature
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	177 tmpidx = idx(find (a(idx) < 1));
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	178 if (! isempty (tmpidx))
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	179 x(tmpidx) = (ys(tmpidx) .* beta (a(tmpidx), b(tmpidx)) .* a(tmpidx)) ...
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	180 .^ (1 ./ a(tmpidx));
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	181 endif
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	182 ## Positive curvature
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	183 tmpidx = idx(find (a(idx) >= 1));
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	184 if (! isempty (tmpidx))
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	185 x(tmpidx) = 1 - ...
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	186 ((1 - ys(tmpidx)) .* beta (a(tmpidx), b(tmpidx)) .* b(tmpidx)) ...
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	187 .^ (1 ./ b(tmpidx));
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	188 endif
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	189
30017 26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	190 ## Cleanup memory before continuing
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	191 clear ys x_i y_i idx tmpidx
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	192
30017 26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	193 if (strcmpi (tail, "lower"))
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	194 x(y == 0) = 0;
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	195 x(y == 1) = 1;
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	196 F = @(x, a, b, y) y - betainc (x, a, b);
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	197 JF = @(x, a, b, Bln) -exp ((a-1) .* log (x) + (b-1) .* log1p (-x) - Bln);
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	198 else
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	199 x(y == 0) = 1;
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	200 x(y == 1) = 0;
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	201 F = @(x, a, b, y) y - betainc (x, a, b, "upper");
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	202 JF = @(x, a, b, Bln) exp ((a-1) .* log (x) + (b-1) .* log1p (-x) - Bln);
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	203 endif
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	204
30017 26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	205 x = newton_method (F, JF, x, a, b, y);
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	206
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	207 ## Restore original shape
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	208 x = reshape (x, orig_sz);
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	209
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	210 endfunction
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	211
31548 c8ad083a5802 maint: Clean up m-files before Octave 8.1 release. Rik <rik@octave.org> parents: 30875 diff changeset	212 function x = newton_method (F, JF, x, a, b, y)
28945 6e460773bdda maint: Use newlines after "function" and before "endfunction" for clarity. Rik <rik@octave.org> parents: 28928 diff changeset	213
30017 26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	214 Bln = betaln (a, b);
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	215 ## Exclude special values that have been already computed.
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	216 todo = find ((y != 0) & (y != 1));
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	217 step = -F (x(todo), a(todo), b(todo), y(todo)) ./ ...
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	218 JF (x(todo), a(todo), b(todo), Bln(todo));
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	219 x_old = x(todo);
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	220 x(todo) += step;
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	221 dx = x(todo) - x_old;
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	222 idx = (dx != 0);
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	223 todo = todo(idx);
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	224 dx_old = dx(idx);
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	225 while (! isempty (todo))
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	226 step = -F (x(todo), a(todo), b(todo), y(todo)) ./ ...
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	227 JF (x(todo), a(todo), b(todo), Bln(todo));
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	228 x_old = x(todo);
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	229 x(todo) += step;
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	230 dx = x(todo) - x_old;
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	231 idx = (abs (dx) < abs (dx_old)); # Converging if dx is getting smaller
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	232 todo = todo(idx);
26bb2cbf6da2 betaincinv.m: Overhaul function for performance and accuracy (bug #60539). Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> parents: 29665 diff changeset	233 dx_old = dx(idx);
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	234 endwhile
28945 6e460773bdda maint: Use newlines after "function" and before "endfunction" for clarity. Rik <rik@octave.org> parents: 28928 diff changeset	235
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	236 endfunction
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	237
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	238
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	239 %!test
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	240 %! x = linspace (0.1, 0.9, 11);
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	241 %! a = [2, 3, 4];
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	242 %! [x,a,b] = ndgrid (x,a,a);
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	243 %! xx = betaincinv (betainc (x, a, b), a, b);
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	244 %! assert (xx, x, 3e-15);
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	245
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	246 %!test
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	247 %! x = linspace (0.1, 0.9, 11);
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	248 %! a = [2, 3, 4];
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	249 %! [x,a,b] = ndgrid (x,a,a);
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	250 %! xx = betaincinv (betainc (x, a, b, "upper"), a, b, "upper");
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	251 %! assert (xx, x, 3e-15);
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	252
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	253 %!test
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	254 %! x = linspace (0.1, 0.9, 11);
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	255 %! a = [0.1:0.1:1];
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	256 %! [x,a,b] = ndgrid (x,a,a);
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	257 %! xx = betaincinv (betainc (x, a, b), a, b);
24996 f80e68529bbf test: relax tolerances on specfun tests to pass on i386 (bug #53437) Mike Miller <mtmiller@octave.org> parents: 24927 diff changeset	258 %! assert (xx, x, 5e-15);
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	259
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	260 %!test
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	261 %! x = linspace (0.1, 0.9, 11);
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	262 %! a = [0.1:0.1:1];
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	263 %! [x,a,b] = ndgrid (x,a,a);
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	264 %! xx = betaincinv (betainc (x, a, b, "upper"), a, b, "upper");
32069 217a35f2a589 betaincinv.m: Slightly increase test tolerance (for macOS 13). Markus Mützel <markus.muetzel@gmx.de> parents: 31706 diff changeset	265 %! assert (xx, x, 24*eps);
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	266
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	267 ## Test the conservation of the input class
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	268 %!assert (class (betaincinv (0.5, 1, 1)), "double")
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	269 %!assert (class (betaincinv (single (0.5), 1, 1)), "single")
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	270 %!assert (class (betaincinv (0.5, single (1), 1)), "single")
29656 7f5bd197fea6 maint: use std::size_t in more instances (bug #60471) Rik <rik@octave.org> parents: 29616 diff changeset	271 %!assert (class (betaincinv (0.5, 1, single (1))), "single")
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	272
30023 374e68721324 betaincinv.m: Add BIST tests for extreme values of inputs. Rik <rik@octave.org> parents: 30017 diff changeset	273 ## Extreme values for y, a, b that really test the algorithm
374e68721324 betaincinv.m: Add BIST tests for extreme values of inputs. Rik <rik@octave.org> parents: 30017 diff changeset	274 %!assert (betaincinv ([0, 1], 1, 3), [0, 1])
30029 40157299fdeb betaincinv.m: Update reference values in BISTs (bug #60528). Markus Mützel <markus.muetzel@gmx.de> parents: 30023 diff changeset	275 %!assert <60528> (betaincinv (1e-6, 1, 3), 3.3333344444450617e-7, 2eps)
40157299fdeb betaincinv.m: Update reference values in BISTs (bug #60528). Markus Mützel <markus.muetzel@gmx.de> parents: 30023 diff changeset	276 %!assert <60528> (betaincinv (1-1e-6, 3, 1), 0.9999996666665555, 2eps)
30023 374e68721324 betaincinv.m: Add BIST tests for extreme values of inputs. Rik <rik@octave.org> parents: 30017 diff changeset	277 %!assert (betainc (betaincinv (0.9, 1e-3, 1), 1e-3, 1), 0.9, 2*eps)
374e68721324 betaincinv.m: Add BIST tests for extreme values of inputs. Rik <rik@octave.org> parents: 30017 diff changeset	278 %!assert (betainc (betaincinv (.01, 1, 1e-3), 1, 1e-3), .01, 6*eps)
374e68721324 betaincinv.m: Add BIST tests for extreme values of inputs. Rik <rik@octave.org> parents: 30017 diff changeset	279 %!assert (betainc (betaincinv (0.5, 100, 1), 100, 1), 0.5, 8*eps)
374e68721324 betaincinv.m: Add BIST tests for extreme values of inputs. Rik <rik@octave.org> parents: 30017 diff changeset	280 %!assert (betainc (betaincinv (0.5, 1, 100), 1, 100), 0.5, 22*eps)
374e68721324 betaincinv.m: Add BIST tests for extreme values of inputs. Rik <rik@octave.org> parents: 30017 diff changeset	281 %!assert (betaincinv ([0, 1], 1, 3, "upper"), [1, 0])
374e68721324 betaincinv.m: Add BIST tests for extreme values of inputs. Rik <rik@octave.org> parents: 30017 diff changeset	282 %!assert <60528> (betaincinv (1e-6, 1, 3, "upper"), 0.99, 2eps)
30029 40157299fdeb betaincinv.m: Update reference values in BISTs (bug #60528). Markus Mützel <markus.muetzel@gmx.de> parents: 30023 diff changeset	283 %!assert <60528> (betaincinv (1-1e-6, 3, 1,"upper"), .01, 250eps)
30023 374e68721324 betaincinv.m: Add BIST tests for extreme values of inputs. Rik <rik@octave.org> parents: 30017 diff changeset	284 %!assert (betainc (betaincinv (0.1, 1e-3, 1, "upper"), 1e-3, 1, "upper"),
374e68721324 betaincinv.m: Add BIST tests for extreme values of inputs. Rik <rik@octave.org> parents: 30017 diff changeset	285 %! 0.1, 2*eps)
374e68721324 betaincinv.m: Add BIST tests for extreme values of inputs. Rik <rik@octave.org> parents: 30017 diff changeset	286 %!assert (betainc (betaincinv (.99, 1, 1e-3, "upper"), 1, 1e-3, "upper"),
374e68721324 betaincinv.m: Add BIST tests for extreme values of inputs. Rik <rik@octave.org> parents: 30017 diff changeset	287 %! .99, 6*eps)
374e68721324 betaincinv.m: Add BIST tests for extreme values of inputs. Rik <rik@octave.org> parents: 30017 diff changeset	288 %!assert (betainc (betaincinv (0.5, 100, 1, "upper"), 100, 1, "upper"),
374e68721324 betaincinv.m: Add BIST tests for extreme values of inputs. Rik <rik@octave.org> parents: 30017 diff changeset	289 %! 0.5, 8*eps)
374e68721324 betaincinv.m: Add BIST tests for extreme values of inputs. Rik <rik@octave.org> parents: 30017 diff changeset	290 %!assert (betainc (betaincinv (0.5, 1, 100, "upper"), 1, 100, "upper"),
374e68721324 betaincinv.m: Add BIST tests for extreme values of inputs. Rik <rik@octave.org> parents: 30017 diff changeset	291 %! 0.5, 22*eps)
29614 26ba91f0eea7 betaincinv.m: Correctly handle small inputs (bug #60528) Rik <rik@octave.org> parents: 29358 diff changeset	292
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	293 ## Test input validation
28896 90fea9cc9caa test: Add expected error message <Invalid call> to BIST tests for nargin. Rik <rik@octave.org> parents: 28789 diff changeset	294 %!error <Invalid call> betaincinv ()
90fea9cc9caa test: Add expected error message <Invalid call> to BIST tests for nargin. Rik <rik@octave.org> parents: 28789 diff changeset	295 %!error <Invalid call> betaincinv (1)
90fea9cc9caa test: Add expected error message <Invalid call> to BIST tests for nargin. Rik <rik@octave.org> parents: 28789 diff changeset	296 %!error <Invalid call> betaincinv (1,2)
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	297 %!error <must be of common size or scalars>
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	298 %! betaincinv (ones (2,2), ones (1,2), 1);
29656 7f5bd197fea6 maint: use std::size_t in more instances (bug #60471) Rik <rik@octave.org> parents: 29616 diff changeset	299 %!error <must be .* floating point> betaincinv ('a', 1, 2)
7f5bd197fea6 maint: use std::size_t in more instances (bug #60471) Rik <rik@octave.org> parents: 29616 diff changeset	300 %!error <must be .* floating point> betaincinv (0, int8 (1), 1)
7f5bd197fea6 maint: use std::size_t in more instances (bug #60471) Rik <rik@octave.org> parents: 29616 diff changeset	301 %!error <must be .* floating point> betaincinv (0, 1, true)
7f5bd197fea6 maint: use std::size_t in more instances (bug #60471) Rik <rik@octave.org> parents: 29616 diff changeset	302 %!error <must be real> betaincinv (0.5i, 1, 2)
7f5bd197fea6 maint: use std::size_t in more instances (bug #60471) Rik <rik@octave.org> parents: 29616 diff changeset	303 %!error <must be real> betaincinv (0, 1i, 1)
7f5bd197fea6 maint: use std::size_t in more instances (bug #60471) Rik <rik@octave.org> parents: 29616 diff changeset	304 %!error <must be real> betaincinv (0, 1, 1i)
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	305 %!error <Y must be in the range \[0, 1\]> betaincinv (-0.1,1,1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	306 %!error <Y must be in the range \[0, 1\]> betaincinv (1.1,1,1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	307 %!error <Y must be in the range \[0, 1\]>
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	308 %! y = ones (1, 1, 2);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	309 %! y(1,1,2) = -1;
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	310 %! betaincinv (y,1,1);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	311 %!error <A must be strictly positive> betaincinv (0.5,0,1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	312 %!error <A must be strictly positive>
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	313 %! a = ones (1, 1, 2);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	314 %! a(1,1,2) = 0;
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	315 %! betaincinv (1,a,1);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	316 %!error <B must be strictly positive> betaincinv (0.5,1,0)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	317 %!error <B must be strictly positive>
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	318 %! b = ones (1, 1, 2);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	319 %! b(1,1,2) = 0;
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	320 %! betaincinv (1,1,b);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24907 diff changeset	321 %!error <invalid value for TAIL> betaincinv (1,2,3, "foobar")

Mercurial > octave

annotate scripts/specfun/betaincinv.m @ 32069:217a35f2a589 stable