Mercurial > octave
annotate scripts/specfun/betaincinv.m @ 24927:c280560d9c96 stable
Overhaul special functions modified by GSOC2018 project.
* NEWS: Add note about new functions added. Add note explaining the changes
done to existing functions.
* __betainc__.cc: Renamed from __betainc_lentz__.cc. Use standard GPL v3
copyright block. Add missing #include "dNDArray.h". Add one-line texinfo
documentation for internal function. Remove fourth input argument to function.
Use is_single_type() to decide whether to operate with FloatNDArray or NDArray.
Delete temporary variables x_arg_s, a_arg_s, b_arg_s used in input validation.
Rename len_x, len_a, len_b to numel_[xab] for clarity. Remove input validation
that numel match (internal function, no need). Rename function output to
retval according to Octave coding conventions. Use more spacing and newlines
for readability of code.
* __expint__.cc: Renamed from __expint_lentz__.cc. Use standard GPL v3
copyright block. Cull list of #includes to just the necessary ones. Use
Complex and FloatComplex typedefs defined by Octave. Add one-line texinfo
documentation for internal function. Remove second input argument to function.
Use is_single_type() to decide whether to operate with FloatComplexNDArray or
ComplexNDArray. Delete temporary variables x_arg_s used in input validation.
Rename len_x to numel_x for clarity. Use constructor with dim_vector and
scalar value rather than fill() after creating array. Rename function output
to retval according to Octave coding conventions. Use more spacing and
newlines for readability of code.
* __gammainc__.cc: Renamed from __gammainc_lentz__.cc. Use standard GPL v3
copyright block. Add one-line texinfo documentation for internal function.
Remove third input argument to function. Remove input validation that numel
match (internal function, no need). Use is_single_type() to decide whether to
operate with FloatNDArray or NDArray. Delete temporary variables x_arg, a_arg
used in input validation. Use constructor with dim_vector and scalar value
rather than fill() after creating array. Rename function output to retval
according to Octave coding conventions. Use more spacing and newlines for
readability of code.
* __betainc_lentz__.cc, __expint_lentz__.cc, __gammainc_lentz__.cc: Removed.
* libinterp/corefcn/module.mk: Add renamed functions to build system.
* betainc.m: Use Octave standard GPL block. Rewrite parts of docstring.
Don't use array brackets around single output of function. Remove isscalar
checks on inputs because common_size() function will already handle it.
Use capital variable names in error messages to match documentation as displayed
in terminal. Reshape all inputs in to column vectors quickly so that input
validation tests that depend on all/any will pass with N-D arrays. Add
comments to code. Check for specific error messages in input validation BIST
tests.
* betaincinv.m: Use Octave standard GPL block. Rewrite parts of docstring.
Don't use array brackets around single output of function. Remove isscalar
checks on inputs because common_size() function will already handle it.
Use capital variable names in error messages to match documentation as displayed
in terminal. Reshape all inputs in to column vectors quickly so that input
validation tests that depend on all/any will pass with N-D arrays. Put most
common case of tail ("lower") first in if/elseif trees. Call functions directly
with function handle rather than using unnecessary feval() call. Use numel
in preference to length. Rename variable i_miss to todo for clarity. Add
comments to code. Check for specific error messages in input validation BIST
tests.
* cosint.m: Use Octave standard GPL block. Rewrite parts of docstring.
Don't use array brackets around single output of function. Remove isscalar
checks on inputs because common_size() function will already handle it.
Add input validation check for isnumeric value. Convert integer classes to
double before proceeding. Rename i_miss to todo for clarity. Use isinf to
detect both -Inf and +Inf rather than separate tests. Use ++it in while
loop conditional to shorten loop blocks. Add Input validation BIST tests.
* expint.m: Remove Sylvain who was not actually an author on this file.
Rewrite parts of docstring. Rename variable sparse_x to orig_sparse.
Eliminate temporary variables res_tmp, x_s_tmp, ssum_tmp. Rename i_miss to
todo. Use Octave coding conventions throughout. Add comments to code.
* gammainc.m: Use Octave standard GPL block. Rewrite parts of docstring.
Remove isscalar checks on inputs because common_size() function will already
handle it. Use capital variable names in error messages to match documentation
as displayed in terminal. Reshape all inputs in to column vectors quickly so
that input validation tests that depend on all/any will pass with N-D arrays.
Put most common case of tail ("lower") first in if/elseif trees. Add input
validation of tail. Rename variable ii to idx for clarity. Rename variable
i_done to todo and switch polarity so that the '!' operator is not required
every time the variable is updated. Use indexing and direct assignment to
update todo rather than logical operator '&' which is slower. Use tolower on
tail variable and then switch strcmpi calls to strcmp. Reformat %!test blocks
in to %!assert blocks to be more compact. Check for specific error messages in
input validation BIST tests.
* gammaincinv.m: Use Octave standard GPL block. Rewrite parts of docstring.
Don't use array brackets around single output of function. Remove isscalar
checks on inputs because common_size() function will already handle it.
Add input validation check for iscomplex value. Use capital variable names in
error messages to match documentation as displayed in terminal. Reshape all
inputs in to column vectors quickly so that input validation tests that depend
on all/any will pass with N-D arrays. Rename i_miss to todo. Use numel
in preference to length. Call functions directly with function handle rather
than using unnecessary feval() call.Rename variable i_miss to todo for clarity.
Use it++ in while loop conditional to shorten loop blocks. Add comments to
code. Check for specific error messages in input validation BIST tests.
Add input validation BIST tests for all error messages.
* sinint.m: Use Octave standard GPL block. Rewrite parts of docstring.
Add input validation for isnumeric. Convert integers to double for
calculation. Reshape input to column vector. Rename variable sz to orig_sz
for clarity. rename i_miss to todo. Reformat BIST tests to mak them more
compact. Add input validation BIST tests.
author | Rik <rik@octave.org> |
---|---|
date | Mon, 19 Mar 2018 10:01:48 -0700 |
parents | bd89440407aa |
children | f80e68529bbf |
rev | line source |
---|---|
24907
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
1 ## Copyright (C) 2017 Michele Ginesi |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
2 ## |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
3 ## 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
|
4 ## |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
5 ## 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
|
6 ## 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
|
7 ## 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
|
8 ## (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
|
9 ## |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
10 ## 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
|
11 ## 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
|
12 ## 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
|
13 ## 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
|
14 ## |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
15 ## 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
|
16 ## 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
|
17 ## <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
|
18 |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
19 ## Author: Michele Ginesi <michele.ginesi@gmail.com> |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
20 |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
21 ## -*- texinfo -*- |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
22 ## @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
|
23 ## @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
|
24 ## @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
|
25 ## 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
|
26 ## |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
27 ## 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
|
28 ## @tex |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
29 ## $$ |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
30 ## 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
|
31 ## $$ |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
32 ## @end tex |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
33 ## @ifnottex |
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 ## @example |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
36 ## @group |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
37 ## x |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
38 ## / |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
39 ## | |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
40 ## 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
|
41 ## | |
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 ## 0 |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
44 ## @end group |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
45 ## @end example |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
46 ## |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
47 ## @end ifnottex |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
48 ## |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
49 ## 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
|
50 ## 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
|
51 ## |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
52 ## 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
|
53 ## @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
|
54 ## |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
55 ## 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
|
56 ## @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
|
57 ## |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
58 ## By default, @var{tail} is @qcode{"lower")} and the inverse of the incomplete |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
59 ## 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
|
60 ## @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
|
61 ## is inverted. |
24907
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
62 ## |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
63 ## 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
|
64 ## @tex |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
65 ## $$ |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
66 ## 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
|
67 ## $$ |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
68 ## @end tex |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
69 ## @ifnottex |
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 ## @example |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
72 ## @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
|
73 ## @end example |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
74 ## |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
75 ## @end ifnottex |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
76 ## |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
77 ## @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
|
78 ## @end deftypefn |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
79 |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
80 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
|
81 |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
82 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
|
83 print_usage (); |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
84 endif |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
85 |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
86 [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
|
87 if (err > 0) |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
88 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
|
89 endif |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
90 |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
91 if (iscomplex (y) || iscomplex (a) || iscomplex (b)) |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
92 error ("betaincinv: all inputs must be real"); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
93 endif |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
94 |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
95 ## FIXME: Should there be isnumeric checking? Right now it accepts char |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
96 ## arrays, but then produces a weird error later on. |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
97 |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
98 ## 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
|
99 orig_sz = size (y); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
100 y = y(:); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
101 a = a(:); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
102 b = b(:); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
103 |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
104 if (any ((y < 0) | (y > 1))) |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
105 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
|
106 endif |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
107 |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
108 if (any (a <= 0)) |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
109 error ("betaincinv: A must be strictly positive"); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
110 endif |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
111 |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
112 if (any (b <= 0)) |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
113 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
|
114 endif |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
115 |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
116 ## 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:
24907
diff
changeset
|
117 if (strcmp (class (y), "single") || strcmp (class (a), "single") |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
118 || strcmp (class (b), "single")) |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
119 y = single (y); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
120 a = single (a); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
121 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
|
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:
24907
diff
changeset
|
124 ## 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
|
125 if (isinteger (y)) |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
126 y = double (y); |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
127 endif |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
128 if (isinteger (a)) |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
129 a = double (a); |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
130 endif |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
131 if (isinteger (b)) |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
132 b = double (b); |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
133 endif |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
134 |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
135 ## Initialize output array |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
136 x = 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
|
137 |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
138 ## Parameters for the Newton method |
24907
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
139 maxit = 20; |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
140 tol = eps (class (y)); |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
141 |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
142 if (strcmpi (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 p = y; |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
144 q = 1 - y; |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
145 x(y == 0) = 0; |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
146 x(y == 1) = 1; |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
147 elseif (strcmpi (tail, "upper")) |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
148 p = 1 - y; |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
149 q = y; |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
150 x(y == 0) = 1; |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
151 x(y == 1) = 0; |
24907
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
152 else |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
153 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
|
154 endif |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
155 |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
156 ## Special values have been already computed. |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
157 todo = (y != 0) & (y != 1); |
24907
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
158 |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
159 ## We will invert the lower version for p < 0.5 and the upper otherwise. |
24907
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
160 i_low = (p < 0.5); |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
161 i_upp = (! i_low); |
24907
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
162 |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
163 idx = todo & i_low; |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
164 if (any (idx)); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
165 n = nnz (idx); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
166 ## Function and derivative of the lower version. |
24907
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
167 F = @(x, a, b, y) y - 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
|
168 JF = @(x, a, b) - real (exp ((a-1) .* log (x) + (b-1) .* log1p (-x) + ... |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
169 gammaln (a+b) - gammaln (a) - gammaln (b))); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
170 |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
171 ## Compute the initial guess with a bisection method of 10 steps. |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
172 x0 = bisection_method (F, zeros (n,1), ones (n,1), ... |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
173 a(i_low), b(i_low), p(i_low), 10); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
174 |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
175 ## Use Newton's method to iteratively find solution. |
24907
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
176 x(i_low) = newton_method (F, JF, x0, a(i_low), b(i_low), p(i_low), ... |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
177 tol, maxit); |
24907
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
178 endif |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
179 |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
180 idx = todo & i_upp; |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
181 if (any (idx)); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
182 n = nnz (idx); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
183 ## Function and derivative of the upper version. |
24907
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
184 F = @(x, a, b, y) y - betainc (x, a, b, "upper"); |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
185 JF = @(x, a, b) real (exp ((a-1) .* log (x) + (b-1) .* log1p (-x) + ... |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
186 gammaln (a+b) - gammaln (a) - gammaln (b))); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
187 |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
188 ## Compute the initial guess with a bisection method of 10 steps. |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
189 x0 = bisection_method (F, zeros (n,1), ones (n,1), ... |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
190 a(i_upp), b(i_upp), q(i_upp), 10); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
191 |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
192 ## Use Newton's method to iteratively find solution. |
24907
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
193 x(i_upp) = newton_method (F, JF, x0, a(i_upp), b(i_upp), q(i_upp), ... |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
194 tol, maxit); |
24907
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
195 endif |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
196 |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
197 ## Restore original shape |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
198 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
|
199 |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
200 endfunction |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
201 |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
202 |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
203 ## Subfunctions: Bisection and Newton Methods |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
204 function xc = bisection_method (F, xl, xr, a, b, y, maxit) |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
205 F_l = F (xl, a, b, y); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
206 F_r = F (xr, a, b, y); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
207 for it = 1:maxit |
24907
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
208 xc = (xl + xr) / 2; |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
209 F_c = F (xc, a, b, y); |
24907
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
210 flag_l = ((F_c .* F_r) < 0); |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
211 flag_r = ((F_c .* F_l) < 0); |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
212 flag_c = (F_c == 0); |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
213 xl(flag_l) = xc(flag_l); |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
214 xr(flag_r) = xc(flag_r); |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
215 xl(flag_c) = xr(flag_c) = xc(flag_c); |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
216 F_l(flag_l) = F_c(flag_l); |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
217 F_r(flag_r) = F_c(flag_r); |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
218 F_l(flag_c) = F_r(flag_c) = 0; |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
219 endfor |
24907
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
220 endfunction |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
221 |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
222 function x = newton_method (F, JF, x0, a, b, y, tol, maxit); |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
223 l = numel (y); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
224 res = -F (x0, a, b, y) ./ JF (x0, a, b); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
225 todo = (abs(res) >= tol * abs (x0)); |
24907
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
226 x = x0; |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
227 it = 0; |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
228 while (any (todo) && (it < maxit)) |
24907
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
229 it++; |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
230 x(todo) += res(todo); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
231 res(todo) = -F(x(todo), a(todo), b(todo), y(todo)) ... |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
232 ./ JF (x(todo), a(todo), b(todo)); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
233 todo = (abs(res) >= tol * abs (x)); |
24907
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
234 endwhile |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
235 x += res; |
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); |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
258 %! 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
|
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"); |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
265 %! 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
|
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") |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
271 %!assert (class (betaincinv (int8 (0), 1, 1)), "double") |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
272 %!assert (class (betaincinv (0.5, int8 (1), 1)), "double") |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
273 %!assert (class (betaincinv (int8 (0), single (1), 1)), "single") |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
274 %!assert (class (betaincinv (single (0.5), int8 (1), 1)), "single") |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
275 |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
276 ## Test input validation |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
277 %!error betaincinv () |
bd89440407aa
Incomplete beta function moved to a .m file, fixing accuracy and
Michele Ginesi <michele.ginesi@gmail.com>
parents:
diff
changeset
|
278 %!error betaincinv (1) |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
279 %!error betaincinv (1,2) |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
280 %!error betaincinv (1,2,3,4,5) |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
281 %!error <must be of common size or scalars> |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
282 %! betaincinv (ones (2,2), ones (1,2), 1); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
283 %!error <all inputs must be real> betaincinv (0.5i, 1, 2) |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
284 %!error <all inputs must be real> betaincinv (0, 1i, 1) |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
285 %!error <all inputs must be real> betaincinv (0, 1, 1i) |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
286 %!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
|
287 %!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
|
288 %!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
|
289 %! y = ones (1, 1, 2); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
290 %! y(1,1,2) = -1; |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
291 %! betaincinv (y,1,1); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
292 %!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
|
293 %!error <A must be strictly positive> |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
294 %! a = ones (1, 1, 2); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
295 %! a(1,1,2) = 0; |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
296 %! betaincinv (1,a,1); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
297 %!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
|
298 %!error <B must be strictly positive> |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
299 %! b = ones (1, 1, 2); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
300 %! b(1,1,2) = 0; |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
301 %! betaincinv (1,1,b); |
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24907
diff
changeset
|
302 %!error <invalid value for TAIL> betaincinv (1,2,3, "foobar") |