Mercurial > octave

annotate scripts/specfun/betainc.m @ 27919:1891570abac8

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