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