Mercurial > hg > octave-lojdl
changeset 14816:0a868d90436b
New function: betaincinv (bug #34364)
author | Axel Mathéi <axel.mathei@gmail.com> |
---|---|
date | Fri, 22 Jun 2012 16:51:31 +0200 |
parents | 95b93a728603 |
children | 67897baaa05f |
files | NEWS doc/interpreter/arith.txi liboctave/lo-specfun.cc liboctave/lo-specfun.h scripts/help/unimplemented.m src/mappers.cc |
diffstat | 6 files changed, 765 insertions(+), 6 deletions(-) [+] |
line wrap: on
line diff
--- a/NEWS +++ b/NEWS @@ -74,9 +74,10 @@ ** Other new functions added in 3.8.0: - colorcube lines splinefit - erfcinv rgbplot tetramesh - findfigs shrinkfaces + betaincinv lines tetramesh + colorcube rgbplot + erfcinv shrinkfaces + findfigs splinefit ** Deprecated functions.
--- a/doc/interpreter/arith.txi +++ b/doc/interpreter/arith.txi @@ -255,6 +255,8 @@ @DOCSTRING(betainc) +@DOCSTRING(betaincinv) + @DOCSTRING(betaln) @DOCSTRING(bincoeff)
--- a/liboctave/lo-specfun.cc +++ b/liboctave/lo-specfun.cc @@ -615,7 +615,7 @@ { FloatComplex retval; - float r = x.real (), i = x.imag(); + float r = x.real (), i = x.imag (); if (fabs (r) < 0.5 && fabs (i) < 0.5) { @@ -2158,6 +2158,28 @@ d1_str.c_str (), d2_str.c_str (), d3_str.c_str ()); } +static void +gripe_betaincinv_nonconformant (octave_idx_type r1, octave_idx_type c1, octave_idx_type r2, octave_idx_type c2, octave_idx_type r3, + octave_idx_type c3) +{ + (*current_liboctave_error_handler) + ("betaincinv: nonconformant arguments (x is %dx%d, a is %dx%d, b is %dx%d)", + r1, c1, r2, c2, r3, c3); +} + +static void +gripe_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) { @@ -2858,7 +2880,7 @@ (*current_liboctave_error_handler) ("gammainc: nonconformant arguments (arg 1 is %s, arg 2 is %s)", - x_str.c_str (), a_str. c_str ()); + x_str.c_str (), a_str.c_str ()); } done: @@ -3124,3 +3146,476 @@ { return erfcx_impl (x); } + +// +// 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 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 = (int) (qq + cx * psq); + + // Use the Soper reduction formula. + + rx = xx / cx; + temp = qq - ai; + if (ns == 0) + { + rx = xx; + } + + for ( ; ; ) + { + term = term * temp * rx / (pp + ai); + value = value + term; + temp = fabs (term); + + if (temp <= acu && temp <= acu * value) + { + value = value * exp (pp * log (xx) + + (qq - 1.0) * log (cx) - beta) / pp; + + if (indx) + { + value = 1.0 - value; + } + break; + } + + ai = ai + 1.0; + ns = ns - 1; + + if (0 <= ns) + { + temp = qq - ai; + if (ns == 0) + { + rx = xx; + } + } + else + { + temp = psq; + 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 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 = 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 = sqrt (- 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 * 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 * pow (1.0 - t + ycur * sqrt (t), 3); + + if (t <= 0.0) + { + value = 1.0 - exp ((log ((1.0 - a) * qq) + beta) / qq); + } + else + { + t = (4.0 * pp + r - 2.0) / t; + + if (t <= 1.0) + { + value = exp ((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 / pow (a, 0.2) - 13.0, sae); + + acu = pow (10.0, iex); + + for ( ; ; ) + { + ycur = betain (value, pp, qq, beta, err); + + if (err) + { + return value; + } + + xin = value; + ycur = (ycur - a) * exp (beta + r * log (xin) + t * 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 = 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 = g / 3.0; + } + + if (tx == value) + { + break; + } + + value = tx; + yprev = ycur; + } + + if (indx) + { + value = 1.0 - value; + } + + return value; +} + +Matrix +betaincinv (double x, double a, const Matrix& b) +{ + octave_idx_type nr = b.rows (); + octave_idx_type nc = b.cols (); + + Matrix retval (nr, nc); + + for (octave_idx_type j = 0; j < nc; j++) + for (octave_idx_type i = 0; i < nr; i++) + retval(i,j) = betaincinv (x, a, b(i,j)); + + return retval; +} + +Matrix +betaincinv (double x, const Matrix& a, double b) +{ + octave_idx_type nr = a.rows (); + octave_idx_type nc = a.cols (); + + Matrix retval (nr, nc); + + for (octave_idx_type j = 0; j < nc; j++) + for (octave_idx_type i = 0; i < nr; i++) + retval(i,j) = betaincinv (x, a(i,j), b); + + return retval; +} + +Matrix +betaincinv (double x, const Matrix& a, const Matrix& b) +{ + Matrix retval; + + octave_idx_type a_nr = a.rows (); + octave_idx_type a_nc = a.cols (); + + octave_idx_type b_nr = b.rows (); + octave_idx_type b_nc = b.cols (); + + if (a_nr == b_nr && a_nc == b_nc) + { + retval.resize (a_nr, a_nc); + + for (octave_idx_type j = 0; j < a_nc; j++) + for (octave_idx_type i = 0; i < a_nr; i++) + retval(i,j) = betaincinv (x, a(i,j), b(i,j)); + } + else + gripe_betaincinv_nonconformant (1, 1, a_nr, a_nc, b_nr, b_nc); + + return retval; +} + +Matrix +betaincinv (const Matrix& x, double a, double b) +{ + octave_idx_type nr = x.rows (); + octave_idx_type nc = x.cols (); + + Matrix retval (nr, nc); + + for (octave_idx_type j = 0; j < nc; j++) + for (octave_idx_type i = 0; i < nr; i++) + retval(i,j) = betaincinv (x(i,j), a, b); + + return retval; +} + +Matrix +betaincinv (const Matrix& x, double a, const Matrix& b) +{ + Matrix retval; + + octave_idx_type nr = x.rows (); + octave_idx_type nc = x.cols (); + + octave_idx_type b_nr = b.rows (); + octave_idx_type b_nc = b.cols (); + + if (nr == b_nr && nc == b_nc) + { + retval.resize (nr, nc); + + for (octave_idx_type j = 0; j < nc; j++) + for (octave_idx_type i = 0; i < nr; i++) + retval(i,j) = betaincinv (x(i,j), a, b(i,j)); + } + else + gripe_betaincinv_nonconformant (nr, nc, 1, 1, b_nr, b_nc); + + return retval; +} + +Matrix +betaincinv (const Matrix& x, const Matrix& a, double b) +{ + Matrix retval; + + octave_idx_type nr = x.rows (); + octave_idx_type nc = x.cols (); + + octave_idx_type a_nr = a.rows (); + octave_idx_type a_nc = a.cols (); + + if (nr == a_nr && nc == a_nc) + { + retval.resize (nr, nc); + + for (octave_idx_type j = 0; j < nc; j++) + for (octave_idx_type i = 0; i < nr; i++) + retval(i,j) = betaincinv (x(i,j), a(i,j), b); + } + else + gripe_betaincinv_nonconformant (nr, nc, a_nr, a_nc, 1, 1); + + return retval; +} + +Matrix +betaincinv (const Matrix& x, const Matrix& a, const Matrix& b) +{ + Matrix retval; + + octave_idx_type nr = x.rows (); + octave_idx_type nc = x.cols (); + + octave_idx_type a_nr = a.rows (); + octave_idx_type a_nc = a.cols (); + + octave_idx_type b_nr = b.rows (); + octave_idx_type b_nc = b.cols (); + + if (nr == a_nr && nr == b_nr && nc == a_nc && nc == b_nc) + { + retval.resize (nr, nc); + + for (octave_idx_type j = 0; j < nc; j++) + for (octave_idx_type i = 0; i < nr; i++) + retval(i,j) = betaincinv (x(i,j), a(i,j), b(i,j)); + } + else + gripe_betaincinv_nonconformant (nr, nc, a_nr, a_nc, b_nr, b_nc); + + return retval; +}
--- a/liboctave/lo-specfun.h +++ b/liboctave/lo-specfun.h @@ -581,4 +581,16 @@ extern OCTAVE_API double erfcx (double x); extern OCTAVE_API float erfcx (float x); +extern OCTAVE_API double betaincinv (double y, double p, double q); +extern OCTAVE_API Matrix betaincinv (double x, double a, const Matrix& b); +extern OCTAVE_API Matrix betaincinv (double x, const Matrix& a, double b); +extern OCTAVE_API Matrix betaincinv (double x, const Matrix& a, const Matrix& b); + +extern OCTAVE_API Matrix betaincinv (const Matrix& x, double a, double b); +extern OCTAVE_API Matrix betaincinv (const Matrix& x, double a, const Matrix& b); +extern OCTAVE_API Matrix betaincinv (const Matrix& x, const Matrix& a, double b); +extern OCTAVE_API Matrix betaincinv (const Matrix& x, const Matrix& a, const Matrix& b); + +extern OCTAVE_API double betain (double x, double p, double q, double beta, bool& err); + #endif
--- a/scripts/help/unimplemented.m +++ b/scripts/help/unimplemented.m @@ -98,7 +98,6 @@ "bar3", "bar3h", "bench", - "betaincinv", "bicgstabl", "brush", "builddocsearchdb",
--- a/src/mappers.cc +++ b/src/mappers.cc @@ -1945,6 +1945,256 @@ %!error tanh (1, 2) */ +DEFUN (betaincinv, args, , + "-*- texinfo -*-\n\ +@deftypefn {Mapping Function} {} betaincinv (@var{y}, @var{z}, @var{w})\n\ +Compute the inverse of the incomplete Beta function, i.e., @var{x} such that\n\ +\n\ +@example\n\ +@var{y} == betainc (@var{x}, @var{z}, @var{w}) \n\ +@end example\n\ +@seealso{beta, betainc, betaln}\n\ +@end deftypefn") +{ + octave_value retval; + + int nargin = args.length (); + + if (nargin == 3) + { + 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 (! error_state) + { + if (b_arg.is_scalar_type ()) + { + float b = b_arg.float_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + else + { + FloatNDArray b = b_arg.float_array_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + } + } + else + { + FloatNDArray a = a_arg.float_array_value (); + + if (! error_state) + { + if (b_arg.is_scalar_type ()) + { + float b = b_arg.float_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + else + { + FloatNDArray b = b_arg.float_array_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + } + } + } + else + { + FloatNDArray x = x_arg.float_array_value (); + + if (a_arg.is_scalar_type ()) + { + float a = a_arg.float_value (); + + if (! error_state) + { + if (b_arg.is_scalar_type ()) + { + float b = b_arg.float_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + else + { + FloatNDArray b = b_arg.float_array_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + } + } + else + { + FloatNDArray a = a_arg.float_array_value (); + + if (! error_state) + { + if (b_arg.is_scalar_type ()) + { + float b = b_arg.float_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + else + { + FloatNDArray b = b_arg.float_array_value (); + + if (! error_state) + retval = betaincinv (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 (! error_state) + { + if (b_arg.is_scalar_type ()) + { + double b = b_arg.double_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + else + { + NDArray b = b_arg.array_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + } + } + else + { + NDArray a = a_arg.array_value (); + + if (! error_state) + { + if (b_arg.is_scalar_type ()) + { + double b = b_arg.double_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + else + { + NDArray b = b_arg.array_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + } + } + } + else + { + NDArray x = x_arg.array_value (); + + if (a_arg.is_scalar_type ()) + { + double a = a_arg.double_value (); + + if (! error_state) + { + if (b_arg.is_scalar_type ()) + { + double b = b_arg.double_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + else + { + NDArray b = b_arg.array_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + } + } + else + { + NDArray a = a_arg.array_value (); + + if (! error_state) + { + if (b_arg.is_scalar_type ()) + { + double b = b_arg.double_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + else + { + NDArray b = b_arg.array_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + } + } + } + } + } + else + print_usage (); + + 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)); + +%!assert (betaincinv (0, 42, 42), 0, sqrt (eps)); +%!assert (betaincinv (1, 42, 42), 1, sqrt (eps)); + +%!error betaincinv() +%!error betaincinv(1, 2) +*/ + DEFUNX ("toascii", Ftoascii, args, , "-*- texinfo -*-\n\ @deftypefn {Mapping Function} {} toascii (@var{s})\n\