Mercurial > octave
changeset 24907:bd89440407aa stable
Incomplete beta function moved to a .m file, fixing accuracy and
input validation (see bugs #51157 and #34405). It's inverse also
has been rewritten as .m file.
--
added libinterp/corefcn/__betainc_lentz__.cc
added scripts/specfun/betainc.m
added scripts/specfun/betaincinv.m
changed libinterp/corefcn/module.mk
changed liboctave/external/slatec-fn/module.mk
changed liboctave/numeric/lo-specfun.cc
changed liboctave/numeric/lo-specfun.h
changed scripts/specfun/module.mk
changed scripts/statistics/distributions/betainv.m
changed scripts/statistics/distributions/binocdf.m
removed libinterp/corefcn/betainc.cc
removed liboctave/external/slatec-fn/betai.f
removed liboctave/external/slatec-fn/dbetai.f
removed liboctave/external/slatec-fn/xbetai.f
removed liboctave/external/slatec-fn/xdbetai.f
| author | Michele Ginesi <michele.ginesi@gmail.com> |
|---|---|
| date | Thu, 07 Sep 2017 16:13:16 +0200 |
| parents | 451f4f291f46 |
| children | 6c082a43abd8 |
| files | libinterp/corefcn/__betainc_lentz__.cc libinterp/corefcn/betainc.cc libinterp/corefcn/module.mk liboctave/external/slatec-fn/betai.f liboctave/external/slatec-fn/dbetai.f liboctave/external/slatec-fn/module.mk liboctave/external/slatec-fn/xbetai.f liboctave/external/slatec-fn/xdbetai.f liboctave/numeric/lo-specfun.cc liboctave/numeric/lo-specfun.h scripts/specfun/betainc.m scripts/specfun/betaincinv.m scripts/specfun/module.mk |
| diffstat | 13 files changed, 625 insertions(+), 1546 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libinterp/corefcn/__betainc_lentz__.cc Thu Sep 07 16:13:16 2017 +0200 @@ -0,0 +1,106 @@ +// Copyright (C) 2017 Michele Ginesi +// +// This file is part of Octave. +// +// Octave is free software; you can redistribute it and/or modify it +// under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 3 of the License, or +// (at your option) any later version. +// +// Octave is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with Octave; see the file COPYING. If not, see +// <http://www.gnu.org/licenses/>. + +#if defined (HAVE_CONFIG_H) +# include "config.h" +#endif + +#include "defun.h" + +DEFUN (__betainc_lentz__, args, , "Continued fraction for incomplete gamma function") +{ + int nargin = args.length (); + + if (nargin != 3) + print_usage (); + else + { + octave_value x_arg = args(0); + octave_value a_arg = args(1); + octave_value b_arg = args(2); + if (x_arg.is_single_type () || a_arg.is_single_type () || b_arg.is_single_type ()) + { + const float x = args(0).float_value (); + const float a = args(1).float_value (); + const float b = args(2).float_value (); + static const float tiny = pow (2, -50); + static const float eps = std::numeric_limits<float>::epsilon(); + float f = tiny; + float C = f; + float D = 0; + float alpha_m = 1; + float beta_m = a - (a * (a+b)) / (a + 1) * x; + float x2 = x * x; + float Delta = 0; + int m = 1; + int maxit = 100; + while((std::abs(Delta - 1) > eps) & (m < maxit)) + { + D = beta_m + alpha_m * D; + if (D == 0) + D = tiny; + C = beta_m + alpha_m / C; + if (C == 0) + C = tiny; + D = 1 / D; + Delta = C * D; + f *= Delta; + alpha_m = ((a + m - 1) * (a + b + m - 1) * (b - m) * m) / ((a + 2 * m - 1) * (a + 2 * m - 1)) * x2; + beta_m = a + 2 * m + ((m * (b - m)) / (a + 2 * m - 1) - ((a + m) * (a + b + m)) / (a + 2 * m + 1)) * x; + m++; + } + if (! error_state) + return octave_value (f); + } + else + { + const double x = args(0).double_value (); + const double a = args(1).double_value (); + const double b = args(2).double_value (); + static const double tiny = pow (2, -100); + static const double eps = std::numeric_limits<double>::epsilon(); + double f = tiny; + double C = f; + double D = 0; + double alpha_m = 1; + double beta_m = a - (a * (a+b)) / (a + 1) * x; + double x2 = x * x; + double Delta = 0; + int m = 1; + int maxit = 200; + while((std::abs(Delta - 1) > eps) & (m < maxit)) + { + D = beta_m + alpha_m * D; + if (D == 0) + D = tiny; + C = beta_m + alpha_m / C; + if (C == 0) + C = tiny; + D = 1 / D; + Delta = C * D; + f *= Delta; + alpha_m = ((a + m - 1) * (a + b + m - 1) * (b - m) * m) / ((a + 2 * m - 1) * (a + 2 * m - 1)) * x2; + beta_m = a + 2 * m + ((m * (b - m)) / (a + 2 * m - 1) - ((a + m) * (a + b + m)) / (a + 2 * m + 1)) * x; + m++; + } + if (! error_state) + return octave_value (f); + } + } + return octave_value_list (); +}
--- a/libinterp/corefcn/betainc.cc Thu Sep 07 15:58:26 2017 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,424 +0,0 @@ -/* - -Copyright (C) 1997-2017 John W. Eaton - -This file is part of Octave. - -Octave is free software: you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by -the Free Software Foundation, either version 3 of the License, or -(at your option) any later version. - -Octave is distributed in the hope that it will be useful, but -WITHOUT ANY WARRANTY; without even the implied warranty of -MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -GNU General Public License for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, see -<https://www.gnu.org/licenses/>. - -*/ - -#if defined (HAVE_CONFIG_H) -# include "config.h" -#endif - -#include "lo-specfun.h" - -#include "defun.h" -#include "error.h" -#include "errwarn.h" -#include "ovl.h" -#include "utils.h" - - -DEFUN (betainc, args, , - doc: /* -*- texinfo -*- -@deftypefn {} {} betainc (@var{x}, @var{a}, @var{b}) -Compute the regularized incomplete Beta function. - -The regularized incomplete Beta function is defined by -@tex -$$ - I (x, a, b) = {1 \over {B (a, b)}} \int_0^x t^{(a-z)} (1-t)^{(b-1)} dt. -$$ -@end tex -@ifnottex -@c Set example in small font to prevent overfull line - -@smallexample -@group - x - 1 / -betainc (x, a, b) = ----------- | t^(a-1) (1-t)^(b-1) dt. - beta (a, b) / - t=0 -@end group -@end smallexample - -@end ifnottex - -If @var{x} has more than one component, both @var{a} and @var{b} must be -scalars. If @var{x} is a scalar, @var{a} and @var{b} must be of -compatible dimensions. -@seealso{betaincinv, beta, betaln} -@end deftypefn */) -{ - if (args.length () != 3) - print_usage (); - - octave_value retval; - - octave_value x_arg = args(0); - octave_value a_arg = args(1); - octave_value b_arg = args(2); - - // FIXME: Can we make a template version of the duplicated code below - if (x_arg.is_single_type () || a_arg.is_single_type () - || b_arg.is_single_type ()) - { - if (x_arg.is_scalar_type ()) - { - float x = x_arg.float_value (); - - if (a_arg.is_scalar_type ()) - { - float a = a_arg.float_value (); - - if (b_arg.is_scalar_type ()) - { - float b = b_arg.float_value (); - - retval = octave::math::betainc (x, a, b); - } - else - { - Array<float> b = b_arg.float_array_value (); - - retval = octave::math::betainc (x, a, b); - } - } - else - { - Array<float> a = a_arg.float_array_value (); - - if (b_arg.is_scalar_type ()) - { - float b = b_arg.float_value (); - - retval = octave::math::betainc (x, a, b); - } - else - { - Array<float> b = b_arg.float_array_value (); - - retval = octave::math::betainc (x, a, b); - } - } - } - else - { - Array<float> x = x_arg.float_array_value (); - - if (a_arg.is_scalar_type ()) - { - float a = a_arg.float_value (); - - if (b_arg.is_scalar_type ()) - { - float b = b_arg.float_value (); - - retval = octave::math::betainc (x, a, b); - } - else - { - Array<float> b = b_arg.float_array_value (); - - retval = octave::math::betainc (x, a, b); - } - } - else - { - Array<float> a = a_arg.float_array_value (); - - if (b_arg.is_scalar_type ()) - { - float b = b_arg.float_value (); - - retval = octave::math::betainc (x, a, b); - } - else - { - Array<float> b = b_arg.float_array_value (); - - retval = octave::math::betainc (x, a, b); - } - } - } - } - else - { - if (x_arg.is_scalar_type ()) - { - double x = x_arg.double_value (); - - if (a_arg.is_scalar_type ()) - { - double a = a_arg.double_value (); - - if (b_arg.is_scalar_type ()) - { - double b = b_arg.double_value (); - - retval = octave::math::betainc (x, a, b); - } - else - { - Array<double> b = b_arg.array_value (); - - retval = octave::math::betainc (x, a, b); - } - } - else - { - Array<double> a = a_arg.array_value (); - - if (b_arg.is_scalar_type ()) - { - double b = b_arg.double_value (); - - retval = octave::math::betainc (x, a, b); - } - else - { - Array<double> b = b_arg.array_value (); - - retval = octave::math::betainc (x, a, b); - } - } - } - else - { - Array<double> x = x_arg.array_value (); - - if (a_arg.is_scalar_type ()) - { - double a = a_arg.double_value (); - - if (b_arg.is_scalar_type ()) - { - double b = b_arg.double_value (); - - retval = octave::math::betainc (x, a, b); - } - else - { - Array<double> b = b_arg.array_value (); - - retval = octave::math::betainc (x, a, b); - } - } - else - { - Array<double> a = a_arg.array_value (); - - if (b_arg.is_scalar_type ()) - { - double b = b_arg.double_value (); - - retval = octave::math::betainc (x, a, b); - } - else - { - Array<double> b = b_arg.array_value (); - - retval = octave::math::betainc (x, a, b); - } - } - } - } - - return retval; -} - -/* -## Double precision -%!test -%! a = [1, 1.5, 2, 3]; -%! b = [4, 3, 2, 1]; -%! v1 = betainc (1,a,b); -%! v2 = [1,1,1,1]; -%! x = [.2, .4, .6, .8]; -%! v3 = betainc (x, a, b); -%! v4 = 1 - betainc (1.-x, b, a); -%! assert (v1, v2, sqrt (eps)); -%! assert (v3, v4, sqrt (eps)); - -## Single precision -%!test -%! a = single ([1, 1.5, 2, 3]); -%! b = single ([4, 3, 2, 1]); -%! v1 = betainc (1,a,b); -%! v2 = single ([1,1,1,1]); -%! x = single ([.2, .4, .6, .8]); -%! v3 = betainc (x, a, b); -%! v4 = 1 - betainc (1.-x, b, a); -%! assert (v1, v2, sqrt (eps ("single"))); -%! assert (v3, v4, sqrt (eps ("single"))); - -## Mixed double/single precision -%!test -%! a = single ([1, 1.5, 2, 3]); -%! b = [4, 3, 2, 1]; -%! v1 = betainc (1,a,b); -%! v2 = single ([1,1,1,1]); -%! x = [.2, .4, .6, .8]; -%! v3 = betainc (x, a, b); -%! v4 = 1-betainc (1.-x, b, a); -%! assert (v1, v2, sqrt (eps ("single"))); -%! assert (v3, v4, sqrt (eps ("single"))); - -%!error betainc () -%!error betainc (1) -%!error betainc (1,2) -%!error betainc (1,2,3,4) -*/ - -DEFUN (betaincinv, args, , - doc: /* -*- texinfo -*- -@deftypefn {} {} betaincinv (@var{y}, @var{a}, @var{b}) -Compute the inverse of the incomplete Beta function. - -The inverse is the value @var{x} such that - -@example -@var{y} == betainc (@var{x}, @var{a}, @var{b}) -@end example -@seealso{betainc, beta, betaln} -@end deftypefn */) -{ - if (args.length () != 3) - print_usage (); - - octave_value retval; - - octave_value x_arg = args(0); - octave_value a_arg = args(1); - octave_value b_arg = args(2); - - if (x_arg.is_scalar_type ()) - { - double x = x_arg.double_value (); - - if (a_arg.is_scalar_type ()) - { - double a = a_arg.double_value (); - - if (b_arg.is_scalar_type ()) - { - double b = b_arg.double_value (); - - retval = octave::math::betaincinv (x, a, b); - } - else - { - Array<double> b = b_arg.array_value (); - - retval = octave::math::betaincinv (x, a, b); - } - } - else - { - Array<double> a = a_arg.array_value (); - - if (b_arg.is_scalar_type ()) - { - double b = b_arg.double_value (); - - retval = octave::math::betaincinv (x, a, b); - } - else - { - Array<double> b = b_arg.array_value (); - - retval = octave::math::betaincinv (x, a, b); - } - } - } - else - { - Array<double> x = x_arg.array_value (); - - if (a_arg.is_scalar_type ()) - { - double a = a_arg.double_value (); - - if (b_arg.is_scalar_type ()) - { - double b = b_arg.double_value (); - - retval = octave::math::betaincinv (x, a, b); - } - else - { - Array<double> b = b_arg.array_value (); - - retval = octave::math::betaincinv (x, a, b); - } - } - else - { - Array<double> a = a_arg.array_value (); - - if (b_arg.is_scalar_type ()) - { - double b = b_arg.double_value (); - - retval = octave::math::betaincinv (x, a, b); - } - else - { - Array<double> b = b_arg.array_value (); - - retval = octave::math::betaincinv (x, a, b); - } - } - } - - // FIXME: It would be better to have an algorithm for betaincinv which - // accepted float inputs and returned float outputs. As it is, we do - // extra work to calculate betaincinv to double precision and then throw - // that precision away. - if (x_arg.is_single_type () || a_arg.is_single_type () - || b_arg.is_single_type ()) - { - retval = Array<float> (retval.array_value ()); - } - - return retval; -} - -/* -%!assert (betaincinv ([0.875 0.6875], [1 2], 3), [0.5 0.5], sqrt (eps)) -%!assert (betaincinv (0.5, 3, 3), 0.5, sqrt (eps)) -%!assert (betaincinv (0.34375, 4, 3), 0.5, sqrt (eps)) -%!assert (betaincinv (0.2265625, 5, 3), 0.5, sqrt (eps)) -%!assert (betaincinv (0.14453125, 6, 3), 0.5, sqrt (eps)) -%!assert (betaincinv (0.08984375, 7, 3), 0.5, sqrt (eps)) -%!assert (betaincinv (0.0546875, 8, 3), 0.5, sqrt (eps)) -%!assert (betaincinv (0.03271484375, 9, 3), 0.5, sqrt (eps)) -%!assert (betaincinv (0.019287109375, 10, 3), 0.5, sqrt (eps)) - -## Test class single as well -%!assert (betaincinv ([0.875 0.6875], [1 2], single (3)), [0.5 0.5], sqrt (eps ("single"))) -%!assert (betaincinv (0.5, 3, single (3)), 0.5, sqrt (eps ("single"))) -%!assert (betaincinv (0.34375, 4, single (3)), 0.5, sqrt (eps ("single"))) - -## Extreme values -%!assert (betaincinv (0, 42, 42), 0, sqrt (eps)) -%!assert (betaincinv (1, 42, 42), 1, sqrt (eps)) - -%!error betaincinv () -%!error betaincinv (1, 2) -*/
--- a/libinterp/corefcn/module.mk Thu Sep 07 15:58:26 2017 +0200 +++ b/libinterp/corefcn/module.mk Thu Sep 07 16:13:16 2017 +0200 @@ -107,6 +107,7 @@ COREFCN_SRC = \ %reldir%/Cell.cc \ + %reldir%/__betainc_lentz__.cc \ %reldir%/__contourc__.cc \ %reldir%/__dsearchn__.cc \ %reldir%/__gammainc_lentz__.cc \ @@ -119,7 +120,6 @@ %reldir%/__qp__.cc \ %reldir%/balance.cc \ %reldir%/besselj.cc \ - %reldir%/betainc.cc \ %reldir%/bitfcns.cc \ %reldir%/bsxfun.cc \ %reldir%/c-file-ptr-stream.cc \
--- a/liboctave/external/slatec-fn/betai.f Thu Sep 07 15:58:26 2017 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,118 +0,0 @@ -*DECK BETAI - REAL FUNCTION BETAI (X, PIN, QIN) -C***BEGIN PROLOGUE BETAI -C***PURPOSE Calculate the incomplete Beta function. -C***LIBRARY SLATEC (FNLIB) -C***CATEGORY C7F -C***TYPE SINGLE PRECISION (BETAI-S, DBETAI-D) -C***KEYWORDS FNLIB, INCOMPLETE BETA FUNCTION, SPECIAL FUNCTIONS -C***AUTHOR Fullerton, W., (LANL) -C***DESCRIPTION -C -C BETAI calculates the REAL incomplete beta function. -C -C The incomplete beta function ratio is the probability that a -C random variable from a beta distribution having parameters PIN and -C QIN will be less than or equal to X. -C -C -- Input Arguments -- All arguments are REAL. -C X upper limit of integration. X must be in (0,1) inclusive. -C PIN first beta distribution parameter. PIN must be .GT. 0.0. -C QIN second beta distribution parameter. QIN must be .GT. 0.0. -C -C***REFERENCES Nancy E. Bosten and E. L. Battiste, Remark on Algorithm -C 179, Communications of the ACM 17, 3 (March 1974), -C pp. 156. -C***ROUTINES CALLED ALBETA, R1MACH, XERMSG -C***REVISION HISTORY (YYMMDD) -C 770401 DATE WRITTEN -C 890531 Changed all specific intrinsics to generic. (WRB) -C 890531 REVISION DATE from Version 3.2 -C 891214 Prologue converted to Version 4.0 format. (BAB) -C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) -C 900326 Removed duplicate information from DESCRIPTION section. -C (WRB) -C 920528 DESCRIPTION and REFERENCES sections revised. (WRB) -C***END PROLOGUE BETAI - LOGICAL FIRST - SAVE EPS, ALNEPS, SML, ALNSML, FIRST - DATA FIRST /.TRUE./ -C***FIRST EXECUTABLE STATEMENT BETAI - IF (FIRST) THEN - EPS = R1MACH(3) - ALNEPS = LOG(EPS) - SML = R1MACH(1) - ALNSML = LOG(SML) - ENDIF - FIRST = .FALSE. -C - IF (X .LT. 0. .OR. X .GT. 1.0) CALL XERMSG ('SLATEC', 'BETAI', - + 'X IS NOT IN THE RANGE (0,1)', 1, 2) - IF (PIN .LE. 0. .OR. QIN .LE. 0.) CALL XERMSG ('SLATEC', 'BETAI', - + 'P AND/OR Q IS LE ZERO', 2, 2) -C - Y = X - P = PIN - Q = QIN - IF (Q.LE.P .AND. X.LT.0.8) GO TO 20 - IF (X.LT.0.2) GO TO 20 - Y = 1.0 - Y - P = QIN - Q = PIN -C - 20 IF ((P+Q)*Y/(P+1.).LT.EPS) GO TO 80 -C -C EVALUATE THE INFINITE SUM FIRST. -C TERM WILL EQUAL Y**P/BETA(PS,P) * (1.-PS)I * Y**I / FAC(I) -C - PS = Q - AINT(Q) - IF (PS.EQ.0.) PS = 1.0 - XB = P*LOG(Y) - ALBETA(PS, P) - LOG(P) - BETAI = 0.0 - IF (XB.LT.ALNSML) GO TO 40 -C - BETAI = EXP (XB) - TERM = BETAI*P - IF (PS.EQ.1.0) GO TO 40 -C - N = MAX (ALNEPS/LOG(Y), 4.0E0) - DO 30 I=1,N - TERM = TERM*(I-PS)*Y/I - BETAI = BETAI + TERM/(P+I) - 30 CONTINUE -C -C NOW EVALUATE THE FINITE SUM, MAYBE. -C - 40 IF (Q.LE.1.0) GO TO 70 -C - XB = P*LOG(Y) + Q*LOG(1.0-Y) - ALBETA(P,Q) - LOG(Q) - IB = MAX (XB/ALNSML, 0.0E0) - TERM = EXP (XB - IB*ALNSML) - C = 1.0/(1.0-Y) - P1 = Q*C/(P+Q-1.) -C - FINSUM = 0.0 - N = Q - IF (Q.EQ.REAL(N)) N = N - 1 - DO 50 I=1,N - IF (P1.LE.1.0 .AND. TERM/EPS.LE.FINSUM) GO TO 60 - TERM = (Q-I+1)*C*TERM/(P+Q-I) -C - IF (TERM.GT.1.0) IB = IB - 1 - IF (TERM.GT.1.0) TERM = TERM*SML -C - IF (IB.EQ.0) FINSUM = FINSUM + TERM - 50 CONTINUE -C - 60 BETAI = BETAI + FINSUM - 70 IF (Y.NE.X .OR. P.NE.PIN) BETAI = 1.0 - BETAI - BETAI = MAX (MIN (BETAI, 1.0), 0.0) - RETURN -C - 80 BETAI = 0.0 - XB = P*LOG(MAX(Y,SML)) - LOG(P) - ALBETA(P,Q) - IF (XB.GT.ALNSML .AND. Y.NE.0.) BETAI = EXP (XB) - IF (Y.NE.X .OR. P.NE.PIN) BETAI = 1.0 - BETAI - RETURN -C - END
--- a/liboctave/external/slatec-fn/dbetai.f Thu Sep 07 15:58:26 2017 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,121 +0,0 @@ - -*DECK DBETAI - DOUBLE PRECISION FUNCTION DBETAI (X, PIN, QIN) -C***BEGIN PROLOGUE DBETAI -C***PURPOSE Calculate the incomplete Beta function. -C***LIBRARY SLATEC (FNLIB) -C***CATEGORY C7F -C***TYPE DOUBLE PRECISION (BETAI-S, DBETAI-D) -C***KEYWORDS FNLIB, INCOMPLETE BETA FUNCTION, SPECIAL FUNCTIONS -C***AUTHOR Fullerton, W., (LANL) -C***DESCRIPTION -C -C DBETAI calculates the DOUBLE PRECISION incomplete beta function. -C -C The incomplete beta function ratio is the probability that a -C random variable from a beta distribution having parameters PIN and -C QIN will be less than or equal to X. -C -C -- Input Arguments -- All arguments are DOUBLE PRECISION. -C X upper limit of integration. X must be in (0,1) inclusive. -C PIN first beta distribution parameter. PIN must be .GT. 0.0. -C QIN second beta distribution parameter. QIN must be .GT. 0.0. -C -C***REFERENCES Nancy E. Bosten and E. L. Battiste, Remark on Algorithm -C 179, Communications of the ACM 17, 3 (March 1974), -C pp. 156. -C***ROUTINES CALLED D1MACH, DLBETA, XERMSG -C***REVISION HISTORY (YYMMDD) -C 770701 DATE WRITTEN -C 890531 Changed all specific intrinsics to generic. (WRB) -C 890911 Removed unnecessary intrinsics. (WRB) -C 890911 REVISION DATE from Version 3.2 -C 891214 Prologue converted to Version 4.0 format. (BAB) -C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) -C 920528 DESCRIPTION and REFERENCES sections revised. (WRB) -C***END PROLOGUE DBETAI - DOUBLE PRECISION X, PIN, QIN, ALNEPS, ALNSML, C, EPS, FINSUM, P, - 1 PS, Q, SML, TERM, XB, XI, Y, D1MACH, DLBETA, P1 - LOGICAL FIRST - SAVE EPS, ALNEPS, SML, ALNSML, FIRST - DATA FIRST /.TRUE./ -C***FIRST EXECUTABLE STATEMENT DBETAI - IF (FIRST) THEN - EPS = D1MACH(3) - ALNEPS = LOG (EPS) - SML = D1MACH(1) - ALNSML = LOG (SML) - ENDIF - FIRST = .FALSE. -C - IF (X .LT. 0.D0 .OR. X .GT. 1.D0) CALL XERMSG ('SLATEC', 'DBETAI', - + 'X IS NOT IN THE RANGE (0,1)', 1, 2) - IF (PIN .LE. 0.D0 .OR. QIN .LE. 0.D0) CALL XERMSG ('SLATEC', - + 'DBETAI', 'P AND/OR Q IS LE ZERO', 2, 2) -C - Y = X - P = PIN - Q = QIN - IF (Q.LE.P .AND. X.LT.0.8D0) GO TO 20 - IF (X.LT.0.2D0) GO TO 20 - Y = 1.0D0 - Y - P = QIN - Q = PIN -C - 20 IF ((P+Q)*Y/(P+1.D0).LT.EPS) GO TO 80 -C -C EVALUATE THE INFINITE SUM FIRST. TERM WILL EQUAL -C Y**P/BETA(PS,P) * (1.-PS)-SUB-I * Y**I / FAC(I) . -C - PS = Q - AINT(Q) - IF (PS.EQ.0.D0) PS = 1.0D0 - XB = P*LOG(Y) - DLBETA(PS,P) - LOG(P) - DBETAI = 0.0D0 - IF (XB.LT.ALNSML) GO TO 40 -C - DBETAI = EXP (XB) - TERM = DBETAI*P - IF (PS.EQ.1.0D0) GO TO 40 - N = MAX (ALNEPS/LOG(Y), 4.0D0) - DO 30 I=1,N - XI = I - TERM = TERM * (XI-PS)*Y/XI - DBETAI = DBETAI + TERM/(P+XI) - 30 CONTINUE -C -C NOW EVALUATE THE FINITE SUM, MAYBE. -C - 40 IF (Q.LE.1.0D0) GO TO 70 -C - XB = P*LOG(Y) + Q*LOG(1.0D0-Y) - DLBETA(P,Q) - LOG(Q) - IB = MAX (XB/ALNSML, 0.0D0) - TERM = EXP(XB - IB*ALNSML) - C = 1.0D0/(1.D0-Y) - P1 = Q*C/(P+Q-1.D0) -C - FINSUM = 0.0D0 - N = Q - IF (Q.EQ.DBLE(N)) N = N - 1 - DO 50 I=1,N - IF (P1.LE.1.0D0 .AND. TERM/EPS.LE.FINSUM) GO TO 60 - XI = I - TERM = (Q-XI+1.0D0)*C*TERM/(P+Q-XI) -C - IF (TERM.GT.1.0D0) IB = IB - 1 - IF (TERM.GT.1.0D0) TERM = TERM*SML -C - IF (IB.EQ.0) FINSUM = FINSUM + TERM - 50 CONTINUE -C - 60 DBETAI = DBETAI + FINSUM - 70 IF (Y.NE.X .OR. P.NE.PIN) DBETAI = 1.0D0 - DBETAI - DBETAI = MAX (MIN (DBETAI, 1.0D0), 0.0D0) - RETURN -C - 80 DBETAI = 0.0D0 - XB = P*LOG(MAX(Y,SML)) - LOG(P) - DLBETA(P,Q) - IF (XB.GT.ALNSML .AND. Y.NE.0.0D0) DBETAI = EXP(XB) - IF (Y.NE.X .OR. P.NE.PIN) DBETAI = 1.0D0 - DBETAI -C - RETURN - END
--- a/liboctave/external/slatec-fn/module.mk Thu Sep 07 15:58:26 2017 +0200 +++ b/liboctave/external/slatec-fn/module.mk Thu Sep 07 16:13:16 2017 +0200 @@ -3,13 +3,11 @@ %reldir%/alngam.f \ %reldir%/alnrel.f \ %reldir%/algams.f \ - %reldir%/betai.f \ %reldir%/csevl.f \ %reldir%/d9gmit.f \ %reldir%/d9lgic.f \ %reldir%/d9lgit.f \ %reldir%/d9lgmc.f \ - %reldir%/dbetai.f \ %reldir%/dcsevl.f \ %reldir%/dgamlm.f \ %reldir%/dgamma.f \ @@ -32,8 +30,6 @@ %reldir%/r9lgmc.f \ %reldir%/r9lgit.f \ %reldir%/r9gmit.f \ - %reldir%/r9lgic.f \ - %reldir%/xdbetai.f \ - %reldir%/xbetai.f + %reldir%/r9lgic.f DIRSTAMP_FILES += %reldir%/$(octave_dirstamp)
--- a/liboctave/external/slatec-fn/xbetai.f Thu Sep 07 15:58:26 2017 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,6 +0,0 @@ - subroutine xbetai (x, a, b, result) - external betai - real x, a, b, result, betai - result = betai (x, a, b) - return - end
--- a/liboctave/external/slatec-fn/xdbetai.f Thu Sep 07 15:58:26 2017 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,6 +0,0 @@ - subroutine xdbetai (x, a, b, result) - external dbetai - double precision x, a, b, result, dbetai - result = dbetai (x, a, b) - return - end
--- a/liboctave/numeric/lo-specfun.cc Thu Sep 07 15:58:26 2017 +0200 +++ b/liboctave/numeric/lo-specfun.cc Thu Sep 07 16:13:16 2017 +0200 @@ -1374,765 +1374,6 @@ #undef NN_BESSEL #undef RC_BESSEL - OCTAVE_NORETURN static void - err_betainc_nonconformant (const dim_vector& d1, const dim_vector& d2, - const dim_vector& d3) - { - std::string d1_str = d1.str (); - std::string d2_str = d2.str (); - std::string d3_str = d3.str (); - - (*current_liboctave_error_handler) - ("betainc: nonconformant arguments (x is %s, a is %s, b is %s)", - d1_str.c_str (), d2_str.c_str (), d3_str.c_str ()); - } - - OCTAVE_NORETURN static void - err_betaincinv_nonconformant (const dim_vector& d1, const dim_vector& d2, - const dim_vector& d3) - { - std::string d1_str = d1.str (); - std::string d2_str = d2.str (); - std::string d3_str = d3.str (); - - (*current_liboctave_error_handler) - ("betaincinv: nonconformant arguments (x is %s, a is %s, b is %s)", - d1_str.c_str (), d2_str.c_str (), d3_str.c_str ()); - } - - double - betainc (double x, double a, double b) - { - double retval; - F77_XFCN (xdbetai, XDBETAI, (x, a, b, retval)); - return retval; - } - - Array<double> - betainc (double x, double a, const Array<double>& b) - { - dim_vector dv = b.dims (); - octave_idx_type nel = dv.numel (); - - Array<double> retval (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x, a, b(i)); - - return retval; - } - - Array<double> - betainc (double x, const Array<double>& a, double b) - { - dim_vector dv = a.dims (); - octave_idx_type nel = dv.numel (); - - Array<double> retval (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x, a(i), b); - - return retval; - } - - Array<double> - betainc (double x, const Array<double>& a, const Array<double>& b) - { - Array<double> retval; - dim_vector dv = a.dims (); - - if (dv != b.dims ()) - err_betainc_nonconformant (dim_vector (0, 0), dv, b.dims ()); - - octave_idx_type nel = dv.numel (); - - retval.resize (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x, a(i), b(i)); - - return retval; - } - - Array<double> - betainc (const Array<double>& x, double a, double b) - { - dim_vector dv = x.dims (); - octave_idx_type nel = dv.numel (); - - Array<double> retval (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x(i), a, b); - - return retval; - } - - Array<double> - betainc (const Array<double>& x, double a, const Array<double>& b) - { - Array<double> retval; - dim_vector dv = x.dims (); - - if (dv != b.dims ()) - err_betainc_nonconformant (dv, dim_vector (0, 0), b.dims ()); - - octave_idx_type nel = dv.numel (); - - retval.resize (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x(i), a, b(i)); - - return retval; - } - - Array<double> - betainc (const Array<double>& x, const Array<double>& a, double b) - { - Array<double> retval; - dim_vector dv = x.dims (); - - if (dv != a.dims ()) - err_betainc_nonconformant (dv, a.dims (), dim_vector (0, 0)); - - octave_idx_type nel = dv.numel (); - - retval.resize (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x(i), a(i), b); - - return retval; - } - - Array<double> - betainc (const Array<double>& x, const Array<double>& a, - const Array<double>& b) - { - Array<double> retval; - dim_vector dv = x.dims (); - - if (dv != a.dims () || dv != b.dims ()) - err_betainc_nonconformant (dv, a.dims (), b.dims ()); - - octave_idx_type nel = dv.numel (); - - retval.resize (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x(i), a(i), b(i)); - - return retval; - } - - float - betainc (float x, float a, float b) - { - float retval; - F77_XFCN (xbetai, XBETAI, (x, a, b, retval)); - return retval; - } - - Array<float> - betainc (float x, float a, const Array<float>& b) - { - dim_vector dv = b.dims (); - octave_idx_type nel = dv.numel (); - - Array<float> retval (dv); - - float *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x, a, b(i)); - - return retval; - } - - Array<float> - betainc (float x, const Array<float>& a, float b) - { - dim_vector dv = a.dims (); - octave_idx_type nel = dv.numel (); - - Array<float> retval (dv); - - float *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x, a(i), b); - - return retval; - } - - Array<float> - betainc (float x, const Array<float>& a, const Array<float>& b) - { - Array<float> retval; - dim_vector dv = a.dims (); - - if (dv != b.dims ()) - err_betainc_nonconformant (dim_vector (0, 0), dv, b.dims ()); - - octave_idx_type nel = dv.numel (); - - retval.resize (dv); - - float *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x, a(i), b(i)); - - return retval; - } - - Array<float> - betainc (const Array<float>& x, float a, float b) - { - dim_vector dv = x.dims (); - octave_idx_type nel = dv.numel (); - - Array<float> retval (dv); - - float *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x(i), a, b); - - return retval; - } - - Array<float> - betainc (const Array<float>& x, float a, const Array<float>& b) - { - Array<float> retval; - dim_vector dv = x.dims (); - - if (dv != b.dims ()) - err_betainc_nonconformant (dv, dim_vector (0, 0), b.dims ()); - - octave_idx_type nel = dv.numel (); - - retval.resize (dv); - - float *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x(i), a, b(i)); - - return retval; - } - - Array<float> - betainc (const Array<float>& x, const Array<float>& a, float b) - { - Array<float> retval; - dim_vector dv = x.dims (); - - if (dv != a.dims ()) - err_betainc_nonconformant (dv, a.dims (), dim_vector (0, 0)); - - octave_idx_type nel = dv.numel (); - - retval.resize (dv); - - float *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x(i), a(i), b); - - return retval; - } - - Array<float> - betainc (const Array<float>& x, const Array<float>& a, - const Array<float>& b) - { - Array<float> retval; - dim_vector dv = x.dims (); - - if (dv != a.dims () || dv != b.dims ()) - err_betainc_nonconformant (dv, a.dims (), b.dims ()); - - octave_idx_type nel = dv.numel (); - - retval.resize (dv); - - float *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x(i), a(i), b(i)); - - return retval; - } - - // - // Incomplete Beta function ratio - // - // Algorithm based on the one by John Burkardt. - // See http://people.sc.fsu.edu/~jburkardt/cpp_src/asa109/asa109.html - // - // The original code is distributed under the GNU LGPL v3 license. - // - // Reference: - // - // KL Majumder, GP Bhattacharjee, - // Algorithm AS 63: - // The incomplete Beta Integral, - // Applied Statistics, - // Volume 22, Number 3, 1973, pages 409-411. - // - double - betain (double x, double p, double q, double beta, bool& err) - { - double acu = 0.1E-14, ai, cx; - bool indx; - int ns; - double pp, psq, qq, rx, temp, term, value, xx; - - value = x; - err = false; - - // Check the input arguments. - - if ((p <= 0.0 || q <= 0.0) || (x < 0.0 || 1.0 < x)) - { - err = true; - return value; - } - - // Special cases. - - if (x == 0.0 || x == 1.0) - { - return value; - } - - // Change tail if necessary and determine S. - - psq = p + q; - cx = 1.0 - x; - - if (p < psq * x) - { - xx = cx; - cx = x; - pp = q; - qq = p; - indx = true; - } - else - { - xx = x; - pp = p; - qq = q; - indx = false; - } - - term = 1.0; - ai = 1.0; - value = 1.0; - ns = static_cast<int> (qq + cx * psq); - - // Use the Soper reduction formula. - - rx = xx / cx; - temp = qq - ai; - if (ns == 0) - { - rx = xx; - } - - for ( ; ; ) - { - term *= temp * rx / (pp + ai); - value += term; - temp = fabs (term); - - if (temp <= acu && temp <= acu * value) - { - value *= exp (pp * std::log (xx) - + (qq - 1.0) * std::log (cx) - beta) / pp; - - if (indx) - { - value = 1.0 - value; - } - break; - } - - ai += 1.0; - ns -= 1; - - if (0 <= ns) - { - temp = qq - ai; - if (ns == 0) - { - rx = xx; - } - } - else - { - temp = psq; - psq += 1.0; - } - } - - return value; - } - - // - // Inverse of the incomplete Beta function - // - // Algorithm based on the one by John Burkardt. - // See http://people.sc.fsu.edu/~jburkardt/cpp_src/asa109/asa109.html - // - // The original code is distributed under the GNU LGPL v3 license. - // - // Reference: - // - // GW Cran, KJ Martin, GE Thomas, - // Remark AS R19 and Algorithm AS 109: - // A Remark on Algorithms AS 63: The Incomplete Beta Integral - // and AS 64: Inverse of the Incomplete Beta Integeral, - // Applied Statistics, - // Volume 26, Number 1, 1977, pages 111-114. - // - double - betaincinv (double y, double p, double q) - { - double a, acu, adj, fpu, g, h; - int iex; - bool indx; - double pp, prev, qq, r, s, sae = -37.0, sq, t, tx, value, w, xin, ycur, yprev; - - double beta = lgamma (p) + lgamma (q) - lgamma (p + q); - bool err = false; - fpu = std::pow (10.0, sae); - value = y; - - // Test for admissibility of parameters. - - if (p <= 0.0 || q <= 0.0) - (*current_liboctave_error_handler) ("betaincinv: wrong parameters"); - if (y < 0.0 || 1.0 < y) - (*current_liboctave_error_handler) ("betaincinv: wrong parameter Y"); - - if (y == 0.0 || y == 1.0) - return value; - - // Change tail if necessary. - - if (0.5 < y) - { - a = 1.0 - y; - pp = q; - qq = p; - indx = true; - } - else - { - a = y; - pp = p; - qq = q; - indx = false; - } - - // Calculate the initial approximation. - - r = std::sqrt (- std::log (a * a)); - - ycur = r - (2.30753 + 0.27061 * r) / (1.0 + (0.99229 + 0.04481 * r) * r); - - if (1.0 < pp && 1.0 < qq) - { - r = (ycur * ycur - 3.0) / 6.0; - s = 1.0 / (pp + pp - 1.0); - t = 1.0 / (qq + qq - 1.0); - h = 2.0 / (s + t); - w = ycur * std::sqrt (h + r) / h - (t - s) * (r + 5.0 / 6.0 - 2.0 / (3.0 * h)); - value = pp / (pp + qq * exp (w + w)); - } - else - { - r = qq + qq; - t = 1.0 / (9.0 * qq); - t = r * std::pow (1.0 - t + ycur * std::sqrt (t), 3); - - if (t <= 0.0) - { - value = 1.0 - exp ((std::log ((1.0 - a) * qq) + beta) / qq); - } - else - { - t = (4.0 * pp + r - 2.0) / t; - - if (t <= 1.0) - { - value = exp ((std::log (a * pp) + beta) / pp); - } - else - { - value = 1.0 - 2.0 / (t + 1.0); - } - } - } - - // Solve for X by a modified Newton-Raphson method, - // using the function BETAIN. - - r = 1.0 - pp; - t = 1.0 - qq; - yprev = 0.0; - sq = 1.0; - prev = 1.0; - - if (value < 0.0001) - { - value = 0.0001; - } - - if (0.9999 < value) - { - value = 0.9999; - } - - iex = std::max (- 5.0 / pp / pp - 1.0 / std::pow (a, 0.2) - 13.0, sae); - - acu = std::pow (10.0, iex); - - for ( ; ; ) - { - ycur = betain (value, pp, qq, beta, err); - - if (err) - { - return value; - } - - xin = value; - ycur = (ycur - a) * exp (beta + r * std::log (xin) - + t * std::log (1.0 - xin)); - - if (ycur * yprev <= 0.0) - { - prev = std::max (sq, fpu); - } - - g = 1.0; - - for ( ; ; ) - { - for ( ; ; ) - { - adj = g * ycur; - sq = adj * adj; - - if (sq < prev) - { - tx = value - adj; - - if (0.0 <= tx && tx <= 1.0) - { - break; - } - } - g /= 3.0; - } - - if (prev <= acu) - { - if (indx) - { - value = 1.0 - value; - } - return value; - } - - if (ycur * ycur <= acu) - { - if (indx) - { - value = 1.0 - value; - } - return value; - } - - if (tx != 0.0 && tx != 1.0) - { - break; - } - - g /= 3.0; - } - - if (tx == value) - { - break; - } - - value = tx; - yprev = ycur; - } - - if (indx) - { - value = 1.0 - value; - } - - return value; - } - - Array<double> - betaincinv (double x, double a, const Array<double>& b) - { - dim_vector dv = b.dims (); - octave_idx_type nel = dv.numel (); - - Array<double> retval (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betaincinv (x, a, b(i)); - - return retval; - } - - Array<double> - betaincinv (double x, const Array<double>& a, double b) - { - dim_vector dv = a.dims (); - octave_idx_type nel = dv.numel (); - - Array<double> retval (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betaincinv (x, a(i), b); - - return retval; - } - - Array<double> - betaincinv (double x, const Array<double>& a, const Array<double>& b) - { - Array<double> retval; - dim_vector dv = a.dims (); - - if (dv != b.dims ()) - err_betaincinv_nonconformant (dim_vector (0, 0), dv, b.dims ()); - - octave_idx_type nel = dv.numel (); - - retval.resize (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betaincinv (x, a(i), b(i)); - - return retval; - } - - Array<double> - betaincinv (const Array<double>& x, double a, double b) - { - dim_vector dv = x.dims (); - octave_idx_type nel = dv.numel (); - - Array<double> retval (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betaincinv (x(i), a, b); - - return retval; - } - - Array<double> - betaincinv (const Array<double>& x, double a, const Array<double>& b) - { - Array<double> retval; - dim_vector dv = x.dims (); - - if (dv != b.dims ()) - err_betaincinv_nonconformant (dv, dim_vector (0, 0), b.dims ()); - - octave_idx_type nel = dv.numel (); - - retval.resize (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betaincinv (x(i), a, b(i)); - - return retval; - } - - Array<double> - betaincinv (const Array<double>& x, const Array<double>& a, double b) - { - Array<double> retval; - dim_vector dv = x.dims (); - - if (dv != a.dims ()) - err_betaincinv_nonconformant (dv, a.dims (), dim_vector (0, 0)); - - octave_idx_type nel = dv.numel (); - - retval.resize (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betaincinv (x(i), a(i), b); - - return retval; - } - - Array<double> - betaincinv (const Array<double>& x, const Array<double>& a, - const Array<double>& b) - { - Array<double> retval; - dim_vector dv = x.dims (); - - if (dv != a.dims () && dv != b.dims ()) - err_betaincinv_nonconformant (dv, a.dims (), b.dims ()); - - octave_idx_type nel = dv.numel (); - - retval.resize (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betaincinv (x(i), a(i), b(i)); - - return retval; - } - Complex biry (const Complex& z, bool deriv, bool scaled, octave_idx_type& ierr) { @@ -3132,6 +2373,7 @@ FloatComplexNDArray airy (const FloatComplexNDArray& z, bool deriv, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::airy (z, deriv, scaled, ierr); } FloatComplexNDArray biry (const FloatComplexNDArray& z, bool deriv, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::biry (z, deriv, scaled, ierr); } +<<<<<<< dest Array<double> betainc (double x, double a, const Array<double>& b) { return octave::math::betainc (x, a, b); } Array<double> betainc (double x, const Array<double>& a, double b) { return octave::math::betainc (x, a, b); } Array<double> betainc (double x, const Array<double>& a, const Array<double>& b) { return octave::math::betainc (x, a, b); }
--- a/liboctave/numeric/lo-specfun.h Thu Sep 07 15:58:26 2017 +0200 +++ b/liboctave/numeric/lo-specfun.h Thu Sep 07 16:13:16 2017 +0200 @@ -282,55 +282,6 @@ extern OCTAVE_API FloatComplexMatrix besselh2 (const FloatRowVector& alpha, const FloatComplexColumnVector& x, bool scaled, Array<octave_idx_type>& ierr); - extern OCTAVE_API double betainc (double x, double a, double b); - extern OCTAVE_API Array<double> betainc (double x, double a, - const Array<double>& b); - extern OCTAVE_API Array<double> betainc (double x, const Array<double>& a, - double b); - extern OCTAVE_API Array<double> betainc (double x, const Array<double>& a, - const Array<double>& b); - extern OCTAVE_API Array<double> betainc (const Array<double>& x, double a, - double b); - extern OCTAVE_API Array<double> betainc (const Array<double>& x, double a, - const Array<double>& b); - extern OCTAVE_API Array<double> betainc (const Array<double>& x, - const Array<double>& a, double b); - extern OCTAVE_API Array<double> betainc (const Array<double>& x, - const Array<double>& a, const Array<double>& b); - - extern OCTAVE_API float betainc (float x, float a, float b); - extern OCTAVE_API Array<float> betainc (float x, float a, - const Array<float>& b); - extern OCTAVE_API Array<float> betainc (float x, const Array<float>& a, - float b); - extern OCTAVE_API Array<float> betainc (float x, const Array<float>& a, - const Array<float>& b); - extern OCTAVE_API Array<float> betainc (const Array<float>& x, float a, - float b); - extern OCTAVE_API Array<float> betainc (const Array<float>& x, float a, - const Array<float>& b); - extern OCTAVE_API Array<float> betainc (const Array<float>& x, - const Array<float>& a, float b); - extern OCTAVE_API Array<float> betainc (const Array<float>& x, - const Array<float>& a, const Array<float>& b); - - extern OCTAVE_API double betaincinv (double x, double a, double b); - extern OCTAVE_API Array<double> betaincinv (double x, double a, - const Array<double>& b); - extern OCTAVE_API Array<double> betaincinv (double x, const Array<double>& a, - double b); - extern OCTAVE_API Array<double> betaincinv (double x, const Array<double>& a, - const Array<double>& b); - - extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, double a, - double b); - extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, double a, - const Array<double>& b); - extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, - const Array<double>& a, double b); - extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, - const Array<double>& a, const Array<double>& b); - extern OCTAVE_API Complex biry (const Complex& z, bool deriv, bool scaled, octave_idx_type& ierr); extern OCTAVE_API ComplexMatrix biry (const ComplexMatrix& z, bool deriv, @@ -771,42 +722,6 @@ OCTAVE_DEPRECATED (4.2, "use 'octave::math::biry' instead") extern OCTAVE_API FloatComplexNDArray biry (const FloatComplexNDArray& z, bool deriv, bool scaled, Array<octave_idx_type>& ierr); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -inline double betainc (double x, double a, double b) { return octave::math::betainc (x, a, b); } -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API Array<double> betainc (double x, double a, const Array<double>& b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API Array<double> betainc (double x, const Array<double>& a, double b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API Array<double> betainc (double x, const Array<double>& a, const Array<double>& b); - -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -inline float betainc (float x, float a, float b) { return octave::math::betainc (x, a, b); } -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API Array<double> betainc (const Array<double>& x, double a, const Array<double>& b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API Array<double> betainc (const Array<double>& x, const Array<double>& a, double b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API Array<double> betainc (const Array<double>& x, const Array<double>& a, const Array<double>& b); - -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API float betainc (float x, float a, float b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API Array<float> betainc (float x, float a, const Array<float>& b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API Array<float> betainc (float x, const Array<float>& a, float b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API Array<float> betainc (float x, const Array<float>& a, const Array<float>& b); - -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API Array<float> betainc (const Array<float>& x, float a, float b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API Array<float> betainc (const Array<float>& x, float a, const Array<float>& b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API Array<float> betainc (const Array<float>& x, const Array<float>& a, float b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API Array<float> betainc (const Array<float>& x, const Array<float>& a, const Array<float>& b); - OCTAVE_DEPRECATED (4.2, "use 'octave::math::gammainc' instead") inline double gammainc (double x, double a, bool& err) { return octave::math::gammainc (x, a, err); } OCTAVE_DEPRECATED (4.2, "use 'octave::math::gammainc' instead") @@ -887,27 +802,6 @@ OCTAVE_DEPRECATED (4.2, "use 'octave::math::dawson' instead") inline FloatComplex dawson (const FloatComplex& x) { return octave::math::dawson (x); } -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betaincinv' instead") -inline double betaincinv (double x, double a, double b) { return octave::math::betaincinv (x, a, b); } - -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betaincinv' instead") -extern OCTAVE_API double betaincinv (double x, double a, double b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betaincinv' instead") -extern OCTAVE_API Array<double> betaincinv (double x, double a, const Array<double>& b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betaincinv' instead") -extern OCTAVE_API Array<double> betaincinv (double x, const Array<double>& a, double b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betaincinv' instead") -extern OCTAVE_API Array<double> betaincinv (double x, const Array<double>& a, const Array<double>& b); - -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betaincinv' instead") -extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, double a, double b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betaincinv' instead") -extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, double a, const Array<double>& b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betaincinv' instead") -extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, const Array<double>& a, double b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betaincinv' instead") -extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, const Array<double>& a, const Array<double>& b); - OCTAVE_DEPRECATED (4.2, "use 'octave::math::ellipj' instead") inline void ellipj (double u, double m, double& sn, double& cn, double& dn, double& err) { octave::math::ellipj (u, m, sn, cn, dn, err); } OCTAVE_DEPRECATED (4.2, "use 'octave::math::ellipj' instead")
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/specfun/betainc.m Thu Sep 07 16:13:16 2017 +0200 @@ -0,0 +1,231 @@ +## Copyright (C) 2017 Michele Ginesi +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or +## (at your option) any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## Authors: Michele Ginesi <michele.ginesi@gmail.com> + +## -*- texinfo -*- +## @deftypefn {} {} betainc (@var{x}, @var{a}, @var{b}) +## @deftypefnx {} {} betainc (@var{x}, @var{a}, @var{b}, @var{tail}) +## Compute the incomplete beta function ratio. +## +## This is defined as +## @tex +## $$ +## I_x (a, b) = {1 \over {B(a,b)}} \displaystyle{\int_0^x t^{a-1} (1-t)^{b-1} dt} +## $$ +## @end tex +## @ifnottex +## +## @example +## @group +## x +## / +## | +## I_x (a, b) = | t^(a-1) (1-t)^(b-1) dt +## | +## / +## 0 +## @end group +## @end example +## +## @end ifnottex +## +## with @var{x} real in [0,1], @var{a} and @var{b} real and strictly positive. +## If one of the input has more than one components, then the others must be +## scalar or of compatible dimensions. +## +## By default or if @var{tail} is @qcode{"lower"} the incomplete beta function +## ratio integrated from 0 to @var{x} is computed. If @var{tail} is +## @qcode{"upper"} then the complementary function integrated from @var{x} to +## 1 is calculated. The two choices are related as +## +## betainc (@var{x}, @var{a}, @var{b}, @qcode{"lower"}) = +## 1 - betainc (@var{x}, @var{a}, @var{b}, @qcode{"upper"}). +## +## Reference +## +## @nospell{A. Cuyt, V. Brevik Petersen, B. Verdonk, H. Waadeland, W.B. Jones} +## @cite{Handbook of Continued Fractions for Special Functions}, +## ch. 18. +## +## @seealso{beta, betaincinv, betaln} +## +## @end deftypefn + +function [y] = betainc (x, a, b, tail = "lower") + + if (nargin > 4 || nargin < 3) + print_usage (); + endif + + if (! isscalar (x) || ! isscalar (a) || ! isscalar (b)) + [err, x, a, b] = common_size (x, a, b); + if (err > 0) + error ("betainc: x, a and b must be of common size or scalars"); + endif + endif + + if (iscomplex (x) || iscomplex (a) || iscomplex (b)) + error ("betainc: inputs must be real or integer"); + endif + + if (any (a <= 0)) + error ("betainc: a must be strictly positive"); + endif + + if (any (b <= 0)) + error ("betainc: b must be strictly positive"); + endif + + if (any (x > 1 | x < 0)) + error ("betainc: x must be between 0 and 1"); + endif + + if (isinteger (x)) + y = double (x); + endif + + if (isinteger (a)) + a = double (a); + endif + + if (isinteger (b)) + b = double (b); + endif + + y = zeros (size (x)); + + # If any of the argument is single, also the output should be + + if (strcmpi (class (y), "single") || strcmpi (class (a), "single") || strcmpi (class (b), "single")) + a = single (a); + b = single (b); + x = single (x); + y = single (y); + endif + + # In the following, we use the fact that the continued fraction we will + # use is more efficient when x <= a / (a+b). + # Moreover, to compute the upper version, which is defined as + # I_x(a,b,"upper") = 1 - I_x(a,b) we used the property + # I_x(a,b) + I_(1-x) (b,a) = 1. + + if (strcmpi (tail, "upper")) + flag = (x < a./(a+b)); + x(!flag) = 1 - x(!flag); + [a(!flag), b(!flag)] = deal (b(!flag), a(!flag)); + elseif (strcmpi (tail, "lower")) + flag = (x > a./(a+b)); + x (flag) = 1 - x(flag); + [a(flag), b(flag)] = deal (b(flag), a(flag)); + else + error ("betainc: invalid value for flag") + endif + + f = zeros (size (x), class (x)); + + ## Continued fractions: CPVWJ, formula 18.5.20, modified Lentz algorithm + ## implemented in a separate .cc file. This particular continued fraction + ## gives (B(a,b) * I_x(a,b)) / (x^a * (1-x)^b). + + for i = 1 : numel (x) + f(i) = __betainc_lentz__ (x(i), a(i), b(i)); + endfor + # We divide the continued fraction by B(a,b) / (x^a * (1-x)^b) + # to obtain I_x(a,b). + y = a .* log (x) + b .* log1p (-x) + gammaln (a + b) - ... + gammaln (a) - gammaln (b) + log (f); + y = real (exp (y)); + y(flag) = 1 - y(flag); + +endfunction + +## Tests from betainc.cc + +# Double precision +%!test +%! a = [1, 1.5, 2, 3]; +%! b = [4, 3, 2, 1]; +%! v1 = betainc (1,a,b); +%! v2 = [1,1,1,1]; +%! x = [.2, .4, .6, .8]; +%! v3 = betainc (x, a, b); +%! v4 = 1 - betainc (1.-x, b, a); +%! assert (v1, v2, sqrt (eps)); +%! assert (v3, v4, sqrt (eps)); + +# Single precision +%!test +%! a = single ([1, 1.5, 2, 3]); +%! b = single ([4, 3, 2, 1]); +%! v1 = betainc (1,a,b); +%! v2 = single ([1,1,1,1]); +%! x = single ([.2, .4, .6, .8]); +%! v3 = betainc (x, a, b); +%! v4 = 1 - betainc (1.-x, b, a); +%! assert (v1, v2, sqrt (eps ("single"))); +%! assert (v3, v4, sqrt (eps ("single"))); + +# Mixed double/single precision +%!test +%! a = single ([1, 1.5, 2, 3]); +%! b = [4, 3, 2, 1]; +%! v1 = betainc (1,a,b); +%! v2 = single ([1,1,1,1]); +%! x = [.2, .4, .6, .8]; +%! v3 = betainc (x, a, b); +%! v4 = 1-betainc (1.-x, b, a); +%! assert (v1, v2, sqrt (eps ("single"))); +%! assert (v3, v4, sqrt (eps ("single"))); + +## New test + +%!test #<51157> +%! y = betainc([0.00780;0.00782;0.00784],250.005,49750.995); +%! y_ex = [0.999999999999989; 0.999999999999992; 0.999999999999995]; +%! assert (y, y_ex, -1e-14); + +%!assert (betainc (0.001, 20, 30), 2.750687665855991e-47, -3e-14); +%!assert (betainc (0.0001, 20, 30), 2.819953178893307e-67, -3e-14); +%!assert (betainc (0.99, 20, 30, "upper"), 1.5671643161872703e-47, -3e-14); +%!assert (betainc (0.999, 20, 30, "upper"), 1.850806276141535e-77, -3e-14); +%!assert (betainc (0.5, 200, 300), 0.9999964565197356, -1e-15); +%!assert (betainc (0.5, 200, 300, "upper"), 3.54348026439253e-06, -1e-13); + +# Test trivial values + +%!test +%! [a,b] = ndgrid (linspace(1e-4, 100, 20), linspace(1e-4, 100, 20)); +%! assert (betainc (0,a,b), zeros(20)); +%! assert (betainc (1,a,b), ones(20)); + +## Test input validation + +%!error betainc () +%!error betainc (1) +%!error betainc (1,2) +%!error betainc (1.1,1,1) +%!error betainc (-0.1,1,1) +%!error betainc (0.5,0,1) +%!error betainc (0.5,1,0) +%!error betainc (1,2,3,4) +%!error betainc (1,2) +%!error betainc (1,2,3,4,5) +%!error betainc (0.5i, 1, 2) +%!error betainc (0, 1i, 1) +%!error betainc (0, 1, 1i)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/specfun/betaincinv.m Thu Sep 07 16:13:16 2017 +0200 @@ -0,0 +1,283 @@ +## Copyright (C) 2017 Michele Ginesi +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## Author: Michele Ginesi <michele.ginesi@gmail.com> + +## -*- texinfo -*- +## @deftypefn {} {} betaincinv (@var{y}, @var{a}, @var{b}) +## @deftypefnx {} {} betaincinv (@var{y}, @var{a}, @var{b}, "lower") +## @deftypefnx {} {} betaincinv (@var{y}, @var{a}, @var{b}, "upper") +## Compute the inverse of the normalized incomplete beta function. +## +## The normalized incomplete beta function is defined as +## @tex +## $$ +## I_x (a, b) = {1 \over {B(a,b)}} \displaystyle{\int_0^x t^{a-1} (1-t)^{b-1} dt} +## $$ +## @end tex +## @ifnottex +## +## @example +## @group +## x +## / +## | +## I_x (a, b) = | t^(a-1) (1-t)^(b-1) dt +## | +## / +## 0 +## @end group +## @end example +## +## @end ifnottex +## +## If two inputs are scalar, then @code{betaincinv (@var{y}, @var{a}, @var{b})} +## is returned for each of the other inputs. +## +## If two or more inputs are not scalar, the sizes of them must agree, and +## @code{betaincinv} is applied element-by-element. +## +## The variable @var{y} must be in the interval [0,1], while @var{a} and @var{b} +## must be real and strictly positive. +## +## By default the inverse of the incomplete beta function integrated from 0 +## to @var{x} is computed. If @qcode{"upper"} is given then the complementary +## function integrated from @var{x} to 1 is inverted. +## +## The function is computed by standard Newton's method, by solving +## @tex +## $$ +## y - I_x (a, b) = 0 +## $$ +## @end tex +## @ifnottex +## +## @example +## @group +## +## @var{y} - betainc (@var{x}, @var{a}, @var{b}) = 0 +## +## @end group +## @end example +## +## @end ifnottex +## +## +## @seealso{beta, betainc, betaln} +## @end deftypefn + +function [x] = betaincinv (y, a, b, tail = "lower") + + if (nargin > 4 || nargin < 3) + print_usage (); + endif + + if (! isscalar (y) || ! isscalar (a) || ! isscalar (b)) + [err, y, a, b] = common_size (y, a, b); + if (err > 0) + error ("betaincinv: y, a must be of common size or scalars"); + endif + endif + + if (iscomplex (y) || iscomplex (a) || iscomplex (b)) + error ("betaincinv: inputs must be real or integer"); + endif + + if (any (a <= 0) | any (b <= 0)) + error ("betaincinv: a must be strictly positive"); + endif + + if (any (y > 1 | y < 0)) + error ("betaincinv: y must be between 0 and 1"); + endif + + if (isinteger (y)) + y = double (y); + endif + + if (isinteger (a)) + a = double (a); + endif + + if (isinteger (b)) + b = double (b); + endif + + + # Extract the size. + sz = size (y); + # Write the inputs as two column vectors. + y = y(:); + a = a(:); + b = b(:); + l = length (y); + # Initialise the output. + x = zeros (l, 1); + + # If one of the inputs is of single type, also the output should be + + if (strcmpi (class (y), "single") || strcmpi (class (a), "single") || strcmpi (class (b), "single")) + a = single (a); + b = single (b); + y = single (y); + x = single (x); + endif + + # Parameters of the Newton method + maxit = 20; + tol = eps (class (y)); + + if (strcmpi (tail, "upper")) + p = 1 - y; + q = y; + x(y == 0) = 1; + x(y == 1) = 0; + elseif (strcmpi (tail, "lower")) + p = y; + q = 1 - y; + x(y == 0) = 0; + x(y == 1) = 1; + else + error ("betaincinv: invalid value for tail") + endif + + # Trivial values have been already computed. + i_miss = ((y != 0) & (y != 1)); + + # We will invert the lower version for p < 0.5 and the upper otherwise. + i_low = (p < 0.5); + i_upp = (!i_low); + + len_low = nnz (i_miss & i_low); + len_upp = nnz (i_miss & i_upp); + + if (any (i_miss & i_low)); + # Function and derivative of the lower version. + F = @(x, a, b, y) y - betainc (x, a, b); + JF = @(x, a, b) - real (exp ((a-1) .* log (x) + (b-1) .* log1p (-x) + ... + gammaln (a+b) - gammaln (a) - gammaln (b))); + # We compute the initial guess with a bisection method of 10 steps. + x0 = bisection_method (F, zeros (len_low,1), ones (len_low,1),... + a(i_low), b(i_low), p(i_low), 10); + x(i_low) = newton_method (F, JF, x0, a(i_low), b(i_low), p(i_low), ... + tol, maxit); + endif + + if (any (i_miss & i_upp)); + # Function and derivative of the lower version. + F = @(x, a, b, y) y - betainc (x, a, b, "upper"); + JF = @(x, a, b) real (exp ((a-1) .* log (x) + (b-1) .* log1p (-x) + ... + gammaln (a+b) - gammaln (a) - gammaln (b))); + # We compute the initial guess with a bisection method of 10 steps. + x0 = bisection_method (F, zeros (len_upp,1), ones (len_upp,1),... + a(i_upp), b(i_upp), q(i_upp), 10); + x(i_upp) = newton_method (F, JF, x0, a(i_upp), b(i_upp), q(i_upp), ... + tol, maxit); + endif + + ## Re-organize the output. + + x = reshape (x, sz); + +endfunction + + +## Subfunctions: Bisection and Newton Methods + +function xc = bisection_method (F, xl, xr, a, b, y, maxit) + F_l = feval (F, xl, a, b, y); + F_r = feval (F, xr, a, b, y); + for it = 1:maxit + xc = (xl + xr) / 2; + F_c = feval (F, xc, a, b, y); + flag_l = ((F_c .* F_r) < 0); + flag_r = ((F_c .* F_l) < 0); + flag_c = (F_c == 0); + xl(flag_l) = xc(flag_l); + xr(flag_r) = xc(flag_r); + xl(flag_c) = xr(flag_c) = xc(flag_c); + F_l(flag_l) = F_c(flag_l); + F_r(flag_r) = F_c(flag_r); + F_l(flag_c) = F_r(flag_c) = 0; + endfor +endfunction + +function x = newton_method (F, JF, x0, a, b, y, tol, maxit); + l = length (y); + res = -feval (F, x0, a, b, y) ./ feval (JF, x0, a, b); + i_miss = (abs(res) >= tol * abs (x0)); + x = x0; + it = 0; + while (any (i_miss) && (it < maxit)) + it++; + x(i_miss) += res(i_miss); + res(i_miss) = - feval (F, x(i_miss), a(i_miss), b(i_miss), ... + y(i_miss)) ./ feval (JF, x(i_miss), a(i_miss), b(i_miss)); + i_miss = (abs(res) >= tol * abs (x)); + endwhile + x += res; +endfunction + + +## Test + +%!test +%! x = linspace(0.1, 0.9, 11); +%! a = [2, 3, 4]; +%! [x,a,b] = ndgrid (x,a,a); +%! xx = betaincinv (betainc (x, a, b), a, b); +%! assert (xx, x, 3e-15); + +%!test +%! x = linspace(0.1, 0.9, 11); +%! a = [2, 3, 4]; +%! [x,a,b] = ndgrid (x,a,a); +%! xx = betaincinv (betainc (x, a, b, "upper"), a, b, "upper"); +%! assert (xx, x, 3e-15); + +%!test +%! x = linspace(0.1, 0.9, 11); +%! a = [0.1:0.1:1]; +%! [x,a,b] = ndgrid (x,a,a); +%! xx = betaincinv (betainc (x, a, b), a, b); +%! assert (xx, x, 3e-15); + +%!test +%! x = linspace(0.1, 0.9, 11); +%! a = [0.1:0.1:1]; +%! [x,a,b] = ndgrid (x,a,a); +%! xx = betaincinv (betainc (x, a, b, "upper"), a, b, "upper"); +%! assert (xx, x, 3e-15); + +## Test on the conservation of the class +%!assert (class (betaincinv (0.5, 1, 1)), "double") +%!assert (class (betaincinv (single (0.5), 1, 1)), "single") +%!assert (class (betaincinv (0.5, single (1), 1)), "single") +%!assert (class (betaincinv (int8 (0), 1, 1)), "double") +%!assert (class (betaincinv (0.5, int8 (1), 1)), "double") +%!assert (class (betaincinv (int8 (0), single (1), 1)), "single") +%!assert (class (betaincinv (single (0.5), int8 (1), 1)), "single") + +## Test input validation +%!error betaincinv () +%!error betaincinv (1) +%!error betaincinv (1,2,3,4) +%!error betaincinv (1, "2") +%!error betaincinv (0.5i, 1, 2) +%!error betaincinv (0, 1i, 1) +%!error betaincinv (0, 1, 1i)