octave: scripts/specfun/betainc.m annotate

annotate scripts/specfun/betainc.m @ 26946:04e5cb5e2cb3

update bug status in tests * __ichol__.cc, file-io.cc, mappers.cc, oct-map.cc, syscalls.cc, asech.m, quadgk.m, deconv.m, poly.m, betainc.m, expint.m, gammainc.m, perms.m, native2unicode.m, bug-38236.tst, bug-47680.tst, bug-50716.tst, classdef.tst, classes.tst, for.tst, nest.tst, struct.tst: Update status of fixed bugs in tests.

author	John W. Eaton <jwe@octave.org>
date	Wed, 20 Mar 2019 02:27:30 +0000
parents	02a9e1bb314a
children	b442ec6dda5c

rev	line source
26376 00f796120a6d maint: Update copyright dates in all source files. John W. Eaton <jwe@octave.org> parents: 25673 diff changeset	1 ## Copyright (C) 2018-2019 Stefan Schlögl
00f796120a6d maint: Update copyright dates in all source files. John W. Eaton <jwe@octave.org> parents: 25673 diff changeset	2 ## Copyright (C) 2018-2019 Michele Ginesi
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	3 ##
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	4 ## 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	5 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	6 ## 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	7 ## under the terms of the GNU General Public License as published by
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	8 ## the Free Software Foundation; either version 3 of the License, or
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	9 ## (at your option) any later version.
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	10 ##
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	11 ## 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	12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	14 ## GNU General Public License for more details.
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	15 ##
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	16 ## 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	17 ## along with Octave; see the file COPYING. If not, see
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	18 ## <https://www.gnu.org/licenses/>.
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 ## -- texinfo --
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	21 ## @deftypefn {} {} betainc (@var{x}, @var{a}, @var{b})
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	22 ## @deftypefnx {} {} betainc (@var{x}, @var{a}, @var{b}, @var{tail})
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	23 ## Compute the incomplete beta function.
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	24 ##
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	25 ## This is defined as
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	26 ## @tex
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	27 ## $$
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	28 ## 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	29 ## $$
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	30 ## @end tex
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	31 ## @ifnottex
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	32 ##
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	33 ## @example
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	34 ## @group
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	35 ## x
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 ## \|
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	38 ## I_x (a, b) = \| 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	39 ## \|
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	40 ## /
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	41 ## 0
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	42 ## @end group
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	43 ## @end example
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	44 ##
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	45 ## @end ifnottex
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	46 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	47 ## with real @var{x} in the range [0,1]. The inputs @var{a} and @var{b} must
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	48 ## be real and strictly positive (> 0). If one of the inputs is not a scalar
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	49 ## then the other inputs must be scalar or of compatible dimensions.
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	50 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	51 ## By default, @var{tail} is @qcode{"lower"} and the incomplete beta function
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	52 ## integrated from 0 to @var{x} is computed. If @var{tail} is @qcode{"upper"}
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	53 ## then the complementary function integrated from @var{x} to 1 is calculated.
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	54 ## The two choices are related by
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	55 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	56 ## betainc (@var{x}, @var{a}, @var{b}, @qcode{"upper"}) =
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	57 ## 1 - betainc (@var{x}, @var{a}, @var{b}, @qcode{"lower"}).
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	58 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	59 ## @code{betainc} uses a more sophisticated algorithm than subtraction to
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	60 ## get numerically accurate results when the @qcode{"lower"} value is small.
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	61 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	62 ## Reference: @nospell{A. Cuyt, V. Brevik Petersen, B. Verdonk, H. Waadeland,
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	63 ## W.B. Jones}, @cite{Handbook of Continued Fractions for Special Functions},
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	64 ## ch.@: 18.
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	65 ##
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	66 ## @seealso{beta, betaincinv, betaln}
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	67 ## @end deftypefn
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	68
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	69 function y = betainc (x, 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	70
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	71 if (nargin < 3 \|\| nargin > 4)
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	72 print_usage ();
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	73 endif
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	74
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	75 [err, x, a, b] = common_size (x, a, b);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	76 if (err > 0)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	77 error ("betainc: X, A, and B must be of common size or scalars");
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	78 endif
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	79
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	80 if (iscomplex (x) \|\| iscomplex (a) \|\| iscomplex (b))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	81 error ("betainc: all inputs must be real");
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	82 endif
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	83
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	84 ## 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: 24923 diff changeset	85 orig_sz = size (x);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	86 x = x(:);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	87 a = a(:);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	88 b = b(:);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	89
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	90 if (any ((x < 0) \| (x > 1)))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	91 error ("betainc: X 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	92 endif
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	93
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	94 if (any (a <= 0))
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	95 error ("betainc: A must be strictly positive");
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	96 endif
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	97
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	98 if (any (b <= 0))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	99 error ("betainc: 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	100 endif
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	101
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	102 ## 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: 24923 diff changeset	103 if (strcmp (class (x), "single") \|\| strcmp (class (a), "single")
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	104 \|\| strcmp (class (b), "single"))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	105 a = single (a);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	106 b = single (b);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	107 x = single (x);
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	108 endif
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	109
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	110 ## Convert to floating point if necessary
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	111 if (isinteger (x))
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	112 y = double (x);
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	113 endif
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	114 if (isinteger (a))
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	115 a = double (a);
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	116 endif
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	117 if (isinteger (b))
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	118 b = double (b);
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	119 endif
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	120
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	121 ## Initialize output array
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	122 y = zeros (size (x), class (x));
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	123
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	124 ## In the following, we use the fact that the continued fraction Octave uses
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	125 ## is more efficient when x <= a / (a + b). Moreover, to compute the upper
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	126 ## version, which is defined as I_x(a,b,"upper") = 1 - I_x(a,b) we use the
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	127 ## property I_x(a,b) + I_(1-x) (b,a) = 1.
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	128
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	129 if (strcmpi (tail, "lower"))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	130 fflag = (x > a./(a+b));
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	131 x(fflag) = 1 - x(fflag);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	132 [a(fflag), b(fflag)] = deal (b(fflag), a(fflag));
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	133 elseif (strcmpi (tail, "upper"))
24923 40ab8be7d7ec Fixed style in specfun scripts Michele Ginesi <michele.ginesi@gmail.com> parents: 24921 diff changeset	134 fflag = (x < (a ./ (a + b)));
40ab8be7d7ec Fixed style in specfun scripts Michele Ginesi <michele.ginesi@gmail.com> parents: 24921 diff changeset	135 x(! fflag) = 1 - x(! fflag);
40ab8be7d7ec Fixed style in specfun scripts Michele Ginesi <michele.ginesi@gmail.com> parents: 24921 diff changeset	136 [a(! fflag), b(! fflag)] = deal (b(! fflag), a(! fflag));
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	137 else
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	138 error ("betainc: 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	139 endif
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	140
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	141 f = zeros (size (x), class (x));
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	142
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	143 ## Continued fractions: CPVWJ, formula 18.5.20, modified Lentz algorithm
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	144 ## implemented in a separate .cc file. This particular continued fraction
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	145 ## gives (B(a,b) * I_x(a,b)) / (x^a * (1-x)^b).
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	146
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	147 f = __betainc__ (x, a, b);
24919 ed6f6bbed604 betainc: vectorized the Lentz's algorithm Michele Ginesi <michele.ginesi@gmail.com> parents: 24907 diff changeset	148
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	149 ## Divide continued fraction by B(a,b) / (x^a * (1-x)^b) to obtain I_x(a,b).
25021 592db8745f26 betainc.m: Use parentheses to control order of addition to reduce round-off error. Rik <rik@octave.org> parents: 25015 diff changeset	150 y = a .* log (x) + b .* log1p (-x) ...
592db8745f26 betainc.m: Use parentheses to control order of addition to reduce round-off error. Rik <rik@octave.org> parents: 25015 diff changeset	151 + (gammaln (a + b) - gammaln (a) - gammaln (b)) + log (f);
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	152 y = real (exp (y));
24923 40ab8be7d7ec Fixed style in specfun scripts Michele Ginesi <michele.ginesi@gmail.com> parents: 24921 diff changeset	153 y(fflag) = 1 - y(fflag);
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	154
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	155 ## Restore original shape
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	156 y = reshape (y, orig_sz);
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	157
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	158 endfunction
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	159
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	160
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	161 ## Double precision
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	162 %!test
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	163 %! a = [1, 1.5, 2, 3];
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	164 %! b = [4, 3, 2, 1];
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	165 %! v1 = betainc (1, a, b);
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	166 %! v2 = [1,1,1,1];
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	167 %! x = [.2, .4, .6, .8];
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	168 %! v3 = betainc (x, a, b);
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	169 %! v4 = 1 - betainc (1-x, b, a);
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	170 %! assert (v1, v2, sqrt (eps));
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	171 %! assert (v3, v4, sqrt (eps));
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	172
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	173 ## Single precision
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	174 %!test
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	175 %! a = single ([1, 1.5, 2, 3]);
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	176 %! b = single ([4, 3, 2, 1]);
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	177 %! v1 = betainc (1, a, b);
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	178 %! v2 = single ([1,1,1,1]);
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	179 %! x = single ([.2, .4, .6, .8]);
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	180 %! v3 = betainc (x, a, b);
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	181 %! v4 = 1 - betainc (1-x, b, a);
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	182 %! assert (v1, v2, sqrt (eps ("single")));
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	183 %! assert (v3, v4, sqrt (eps ("single")));
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	184
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	185 ## Mixed double/single precision
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	186 %!test
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	187 %! a = single ([1, 1.5, 2, 3]);
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	188 %! b = [4, 3, 2, 1];
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	189 %! v1 = betainc (1,a,b);
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	190 %! v2 = single ([1,1,1,1]);
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	191 %! x = [.2, .4, .6, .8];
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	192 %! v3 = betainc (x, a, b);
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	193 %! v4 = 1-betainc (1.-x, b, a);
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	194 %! assert (v1, v2, sqrt (eps ("single")));
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	195 %! assert (v3, v4, sqrt (eps ("single")));
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	196
26946 04e5cb5e2cb3 update bug status in tests John W. Eaton <jwe@octave.org> parents: 26514 diff changeset	197 %!test <*51157>
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	198 %! y = betainc ([0.00780;0.00782;0.00784],250.005,49750.995);
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	199 %! y_ex = [0.999999999999989; 0.999999999999992; 0.999999999999995];
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	200 %! assert (y, y_ex, -1e-14);
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	201
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	202 %!assert (betainc (0.001, 20, 30), 2.750687665855991e-47, -3e-14);
24996 f80e68529bbf test: relax tolerances on specfun tests to pass on i386 (bug #53437) Mike Miller <mtmiller@octave.org> parents: 24927 diff changeset	203 %!assert (betainc (0.0001, 20, 30), 2.819953178893307e-67, -7e-14);
26946 04e5cb5e2cb3 update bug status in tests John W. Eaton <jwe@octave.org> parents: 26514 diff changeset	204 %!assert <*54383> (betainc (0.99, 20, 30, "upper"), 1.5671643161872703e-47, -7e-14);
24996 f80e68529bbf test: relax tolerances on specfun tests to pass on i386 (bug #53437) Mike Miller <mtmiller@octave.org> parents: 24927 diff changeset	205 %!assert (betainc (0.999, 20, 30, "upper"), 1.850806276141535e-77, -7e-14);
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	206 %!assert (betainc (0.5, 200, 300), 0.9999964565197356, -1e-15);
24996 f80e68529bbf test: relax tolerances on specfun tests to pass on i386 (bug #53437) Mike Miller <mtmiller@octave.org> parents: 24927 diff changeset	207 %!assert (betainc (0.5, 200, 300, "upper"), 3.54348026439253e-06, -3e-13);
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	208
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	209 ## Test trivial values
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	210 %!test
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	211 %! [a,b] = ndgrid (linspace (1e-4, 100, 20), linspace (1e-4, 100, 20));
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	212 %! assert (betainc (0, a, b), zeros (20));
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	213 %! assert (betainc (1, a, b), ones (20));
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	214
26946 04e5cb5e2cb3 update bug status in tests John W. Eaton <jwe@octave.org> parents: 26514 diff changeset	215 %!test <*34405>
25015 baa7e37453b1 Added more tests for betainc and expint. Michele Ginesi <michele.ginesi@gmail.com> parents: 24996 diff changeset	216 %! assert (betainc (NaN, 1, 2), NaN);
baa7e37453b1 Added more tests for betainc and expint. Michele Ginesi <michele.ginesi@gmail.com> parents: 24996 diff changeset	217 %! assert (betainc (0.5, 1, Inf), NaN);
baa7e37453b1 Added more tests for betainc and expint. Michele Ginesi <michele.ginesi@gmail.com> parents: 24996 diff changeset	218
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	219 ## Test input validation
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	220 %!error betainc ()
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	221 %!error betainc (1)
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	222 %!error betainc (1,2)
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	223 %!error betainc (1,2,3,4,5)
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	224 %!error <must be of common size or scalars> betainc (ones (2,2), ones (1,2), 1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	225 %!error <all inputs must be real> betainc (0.5i, 1, 2)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	226 %!error <all inputs must be real> betainc (0, 1i, 1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	227 %!error <all inputs must be real> betainc (0, 1, 1i)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	228 %!error <X must be in the range \[0, 1\]> betainc (-0.1,1,1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	229 %!error <X must be in the range \[0, 1\]> betainc (1.1,1,1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	230 %!error <X must be in the range \[0, 1\]>
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	231 %! x = ones (1, 1, 2);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	232 %! x(1,1,2) = -1;
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	233 %! betainc (x,1,1);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	234 %!error <A must be strictly positive> betainc (0.5,0,1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	235 %!error <A must be strictly positive>
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	236 %! a = ones (1, 1, 2);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	237 %! a(1,1,2) = 0;
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	238 %! betainc (1,a,1);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	239 %!error <B must be strictly positive> betainc (0.5,1,0)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	240 %!error <B must be strictly positive>
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	241 %! b = ones (1, 1, 2);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	242 %! b(1,1,2) = 0;
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	243 %! betainc (1,1,b);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	244 %!error <invalid value for TAIL> betainc (1,2,3, "foobar")

Mercurial > octave

annotate scripts/specfun/betainc.m @ 26946:04e5cb5e2cb3