octave: scripts/specfun/betainc.m annotate

annotate scripts/specfun/betainc.m @ 27923:bd51beb6205e

update formatting of copyright notices * Use <https://octave.org/copyright/> instead of <https://octave.org/COPYRIGHT.html/>. * For consistency with other comments in the Octave sources, use C++-style comments for copyright blocks in C and C++ files. * Use delimiters above and below copyright blocks that are appropriate for the language used in the file. * Eliminate extra spacing inside copyright blocks. * lex.ll (looks_like_copyright): Also allow newlines and carriage returns before the word "Copyright". * scripts/mk-doc.pl (gethelp): Also skip empty comment lines. * bp-table.cc, type.m: Adjust tests.

author	John W. Eaton <jwe@octave.org>
date	Wed, 08 Jan 2020 11:59:41 -0500
parents	1891570abac8
children	28de41192f3c 0a5b15007766

rev	line source
27923 bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	1 ########################################################################
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	2 ##
27919 1891570abac8 update Octave Project Developers copyright for the new year John W. Eaton <jwe@octave.org> parents: 27918 diff changeset	3 ## Copyright (C) 2018-2020 The Octave Project Developers
27918 b442ec6dda5c use centralized file for copyright info for individual contributors John W. Eaton <jwe@octave.org> parents: 26946 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: 24923 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
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	12 ## 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	13 ## (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	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
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	17 ## 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	18 ## 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	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: 24923 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 --
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	27 ## @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	28 ## @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	29 ## 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	30 ##
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	31 ## 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	32 ## @tex
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	33 ## $$
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	34 ## 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	35 ## $$
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	36 ## @end tex
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	37 ## @ifnottex
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	38 ##
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	39 ## @example
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	40 ## @group
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	41 ## x
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	42 ## /
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	43 ## \|
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	44 ## 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	45 ## \|
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	46 ## /
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	47 ## 0
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	48 ## @end group
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	49 ## @end example
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	50 ##
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	51 ## @end ifnottex
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	52 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	53 ## 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	54 ## 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	55 ## 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	56 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	57 ## 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	58 ## 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	59 ## 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	60 ## 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	61 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	62 ## 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	63 ## 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	64 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	65 ## @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	66 ## 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	67 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	68 ## 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	69 ## 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	70 ## ch.@: 18.
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	71 ##
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	72 ## @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	73 ## @end deftypefn
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 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	76
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	77 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	78 print_usage ();
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	79 endif
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	80
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	81 [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	82 if (err > 0)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	83 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	84 endif
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	85
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	86 if (iscomplex (x) \|\| iscomplex (a) \|\| iscomplex (b))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	87 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	88 endif
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	89
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	90 ## 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	91 orig_sz = size (x);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	92 x = x(:);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	93 a = a(:);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	94 b = b(:);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	95
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	96 if (any ((x < 0) \| (x > 1)))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	97 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	98 endif
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	99
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	100 if (any (a <= 0))
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	101 error ("betainc: A must be strictly positive");
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	102 endif
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	103
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	104 if (any (b <= 0))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	105 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	106 endif
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	107
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	108 ## 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	109 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	110 \|\| strcmp (class (b), "single"))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	111 a = single (a);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	112 b = single (b);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	113 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	114 endif
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	115
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	116 ## 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	117 if (isinteger (x))
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	118 y = double (x);
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 if (isinteger (a))
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	121 a = double (a);
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	122 endif
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	123 if (isinteger (b))
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	124 b = double (b);
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	125 endif
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	126
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	127 ## Initialize output array
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	128 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	129
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	130 ## 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	131 ## 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	132 ## 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	133 ## 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	134
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	135 if (strcmpi (tail, "lower"))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	136 fflag = (x > a./(a+b));
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	137 x(fflag) = 1 - x(fflag);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	138 [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	139 elseif (strcmpi (tail, "upper"))
24923 40ab8be7d7ec Fixed style in specfun scripts Michele Ginesi <michele.ginesi@gmail.com> parents: 24921 diff changeset	140 fflag = (x < (a ./ (a + b)));
40ab8be7d7ec Fixed style in specfun scripts Michele Ginesi <michele.ginesi@gmail.com> parents: 24921 diff changeset	141 x(! fflag) = 1 - x(! fflag);
40ab8be7d7ec Fixed style in specfun scripts Michele Ginesi <michele.ginesi@gmail.com> parents: 24921 diff changeset	142 [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	143 else
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	144 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	145 endif
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	146
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	147 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	148
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	149 ## 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	150 ## 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	151 ## 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	152
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	153 f = __betainc__ (x, a, b);
24919 ed6f6bbed604 betainc: vectorized the Lentz's algorithm Michele Ginesi <michele.ginesi@gmail.com> parents: 24907 diff changeset	154
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	155 ## 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	156 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	157 + (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	158 y = real (exp (y));
24923 40ab8be7d7ec Fixed style in specfun scripts Michele Ginesi <michele.ginesi@gmail.com> parents: 24921 diff changeset	159 y(fflag) = 1 - y(fflag);
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	160
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	161 ## Restore original shape
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	162 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	163
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	164 endfunction
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	165
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	166
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	167 ## Double precision
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	168 %!test
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	169 %! 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	170 %! b = [4, 3, 2, 1];
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	171 %! 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	172 %! 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	173 %! 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	174 %! v3 = betainc (x, a, b);
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	175 %! 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	176 %! 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	177 %! 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	178
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	179 ## Single precision
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	180 %!test
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	181 %! 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	182 %! b = single ([4, 3, 2, 1]);
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	183 %! 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	184 %! 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	185 %! 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	186 %! v3 = betainc (x, a, b);
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	187 %! 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	188 %! 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	189 %! 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	190
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	191 ## 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	192 %!test
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	193 %! 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	194 %! 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	195 %! 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	196 %! 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	197 %! 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	198 %! 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	199 %! 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	200 %! 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	201 %! 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	202
26946 04e5cb5e2cb3 update bug status in tests John W. Eaton <jwe@octave.org> parents: 26514 diff changeset	203 %!test <*51157>
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	204 %! 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	205 %! 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	206 %! 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	207
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	208 %!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	209 %!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	210 %!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	211 %!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	212 %!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	213 %!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	214
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	215 ## 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	216 %!test
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	217 %! [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	218 %! assert (betainc (0, a, b), zeros (20));
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	219 %! 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	220
26946 04e5cb5e2cb3 update bug status in tests John W. Eaton <jwe@octave.org> parents: 26514 diff changeset	221 %!test <*34405>
25015 baa7e37453b1 Added more tests for betainc and expint. Michele Ginesi <michele.ginesi@gmail.com> parents: 24996 diff changeset	222 %! assert (betainc (NaN, 1, 2), NaN);
baa7e37453b1 Added more tests for betainc and expint. Michele Ginesi <michele.ginesi@gmail.com> parents: 24996 diff changeset	223 %! 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	224
24907 bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	225 ## Test input validation
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	226 %!error betainc ()
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	227 %!error betainc (1)
bd89440407aa Incomplete beta function moved to a .m file, fixing accuracy and Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	228 %!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	229 %!error betainc (1,2,3,4,5)
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	230 %!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	231 %!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	232 %!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	233 %!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	234 %!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	235 %!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	236 %!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	237 %! x = ones (1, 1, 2);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	238 %! x(1,1,2) = -1;
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	239 %! betainc (x,1,1);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	240 %!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	241 %!error <A must be strictly positive>
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	242 %! a = ones (1, 1, 2);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	243 %! a(1,1,2) = 0;
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	244 %! betainc (1,a,1);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	245 %!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	246 %!error <B must be strictly positive>
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	247 %! b = ones (1, 1, 2);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	248 %! b(1,1,2) = 0;
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	249 %! betainc (1,1,b);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24923 diff changeset	250 %!error <invalid value for TAIL> betainc (1,2,3, "foobar")

Mercurial > octave

annotate scripts/specfun/betainc.m @ 27923:bd51beb6205e