Mercurial > octave
changeset 30016:ad6a57b215e8
betainc.m: Improve accuracy for certain special integer inputs (bug #60682).
Use exact computation for trivial values of parameters a and b. Otherwise,
fall back on continued fraction expansion (default algorithm).
* NEWS: Announce better accuracy for special cases.
* betainc.m: Check input for 3 special cases: 1) a = 1, b = 1;
2) a = 1, b = anything; 3) a = anything, b = 1. Compute these
exactly. Use indexing to select remaining non-trivial values
from input and calculate these with continued fraction expansion.
| author | Michele Ginesi <michele.ginesi@gmail.com> |
|---|---|
| date | Fri, 28 May 2021 18:26:05 +0200 |
| parents | af135277dc50 |
| children | 26bb2cbf6da2 |
| files | NEWS scripts/specfun/betainc.m |
| diffstat | 2 files changed, 36 insertions(+), 6 deletions(-) [+] |
line wrap: on
line diff
--- a/NEWS Thu Aug 19 16:36:35 2021 -0400 +++ b/NEWS Fri May 28 18:26:05 2021 +0200 @@ -81,6 +81,9 @@ Matlab-compatible and yields results which more nearly minimize `norm (A*x - b)`. Previously, Octave computed a minimum-norm solution. +- The `betainc` function now calculates an exact output for the +important special cases where a or b are 1. + - The `whos` function now displays an additional attribute 's' when the variable is a sparse type.
--- a/scripts/specfun/betainc.m Thu Aug 19 16:36:35 2021 -0400 +++ b/scripts/specfun/betainc.m Fri May 28 18:26:05 2021 +0200 @@ -127,13 +127,38 @@ ## Initialize output array y = zeros (size (x), class (x)); + ## Trivial cases (long code here trades memory for speed) + a_one = (a == 1); + b_one = (b == 1); + a_b_one = a_one & b_one; + a_not_one = ! a_one; + b_not_one = ! b_one; + non_trivial = a_not_one & b_not_one; + a_one &= b_not_one; + b_one &= a_not_one; + + if (strcmpi (tail, "lower")) + y(a_b_one) = x(a_b_one); + y(a_one) = 1 - (1 - x(a_one)) .^ b(a_one); + y(b_one) = x(b_one) .^ a(b_one); + elseif (strcmpi (tail, "upper")) + y(a_b_one) = 1 - x(a_b_one); + y(a_one) = (1 - x(a_one)) .^ b(a_one); + y(b_one) = 1 - x(b_one) .^ a(b_one); + endif + + ## Non-Trivial cases ## In the following, we use the fact that the continued fraction Octave uses ## 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 use the ## property I_x(a,b) + I_(1-x) (b,a) = 1. + + x = x(non_trivial); + a = a(non_trivial); + b = b(non_trivial); if (strcmpi (tail, "lower")) - fflag = (x > a./(a+b)); + fflag = (x > a./(a + b)); x(fflag) = 1 - x(fflag); [a(fflag), b(fflag)] = deal (b(fflag), a(fflag)); elseif (strcmpi (tail, "upper")) @@ -153,10 +178,12 @@ f = __betainc__ (x, a, b); ## Divide 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(fflag) = 1 - y(fflag); + y_nt = a .* log (x) + b .* log1p (-x) ... + + (gammaln (a + b) - gammaln (a) - gammaln (b)) + log (f); + y_nt = real (exp (y_nt)); + y_nt(fflag) = 1 - y_nt(fflag); + + y(non_trivial) = y_nt; ## Restore original shape y = reshape (y, orig_sz); @@ -220,7 +247,7 @@ %!test <*34405> %! assert (betainc (NaN, 1, 2), NaN); -%! assert (betainc (0.5, 1, Inf), NaN); +%! assert (betainc (0.5, 1, Inf), 1); ## Test input validation %!error <Invalid call> betainc ()