Mercurial > octave

comparison scripts/specfun/betainc.m @ 30947:fe898ae23e0e stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
betainc.m: Use sophisticated technique for calculating exponents to avoid innacuracies (bug #62329) * betainc.m: Use expm1() and log1p() to replace exponentiation operator ".^" in special case code when a or b equals 1. Add BIST tests for bug #62329.
author Nir Krakauer <nkrakauer@ccny.cuny.edu>
date Tue, 19 Apr 2022 15:09:27 -0700
parents 796f54d4ddbf
children 74f8854958d2
comparison
equal deleted inserted replaced
30938:9c7822f4fc9c 30947:fe898ae23e0e
137 a_one &= b_not_one; 137 a_one &= b_not_one;
138 b_one &= a_not_one; 138 b_one &= a_not_one;
139 139
140 if (strcmpi (tail, "lower")) 140 if (strcmpi (tail, "lower"))
141 y(a_b_one) = x(a_b_one); 141 y(a_b_one) = x(a_b_one);
142 y(a_one) = 1 - (1 - x(a_one)) .^ b(a_one); 142 ## See bug #62329.
143 ## equivalent to "1 - (1 - x(a_one)) .^ b(a_one)", but less roundoff error
144 y(a_one) = - expm1 (log1p (- x(a_one)) .* b(a_one));
143 y(b_one) = x(b_one) .^ a(b_one); 145 y(b_one) = x(b_one) .^ a(b_one);
144 elseif (strcmpi (tail, "upper")) 146 elseif (strcmpi (tail, "upper"))
145 y(a_b_one) = 1 - x(a_b_one); 147 y(a_b_one) = 1 - x(a_b_one);
146 y(a_one) = (1 - x(a_one)) .^ b(a_one); 148 ## equivalent to "(1 - x(a_one)) .^ b(a_one)", but less roundoff error
147 y(b_one) = 1 - x(b_one) .^ a(b_one); 149 y(a_one) = exp (log1p (- x(a_one)) .* b(a_one));
150 ## equivalent to "1 - x(b_one) .^ a(b_one)", but less roundoff error
151 y(b_one) = - expm1 (log (x(b_one)) .* a(b_one));
148 endif 152 endif
149 153
150 ## Non-Trivial cases 154 ## Non-Trivial cases
151 ## In the following, we use the fact that the continued fraction Octave uses 155 ## In the following, we use the fact that the continued fraction Octave uses
152 ## is more efficient when x <= a / (a + b). Moreover, to compute the upper 156 ## is more efficient when x <= a / (a + b). Moreover, to compute the upper
243 ## Test trivial values 247 ## Test trivial values
244 %!test 248 %!test
245 %! [a,b] = ndgrid (linspace (1e-4, 100, 20), linspace (1e-4, 100, 20)); 249 %! [a,b] = ndgrid (linspace (1e-4, 100, 20), linspace (1e-4, 100, 20));
246 %! assert (betainc (0, a, b), zeros (20)); 250 %! assert (betainc (0, a, b), zeros (20));
247 %! assert (betainc (1, a, b), ones (20)); 251 %! assert (betainc (1, a, b), ones (20));
252 %! assert (betainc (0, a, b, "upper"), ones (20));
253 %! assert (betainc (1, a, b, "upper"), zeros (20));
248 254
249 %!test <*34405> 255 %!test <*34405>
250 %! assert (betainc (NaN, 1, 2), NaN); 256 %! assert (betainc (NaN, 1, 2), NaN);
251 %! assert (betainc (0.5, 1, Inf), 1); 257 %! assert (betainc (0.5, 1, Inf), 1);
258
259 %!test <*62329>
260 %! assert (betainc (2e-20, 1, 0.5), 1e-20, -1e-15)
261 %! assert (betainc (2e-5, 1, 0.5), 2e-5 / (1 + sqrt (1 - 2e-5)), -1e-15)
262 %! assert (betainc (0.99, 1, 0.5, "upper"), 0.1, -1e-15)
263 %! assert (betainc (0.99, 0.5, 1, "upper"), - expm1 (log (0.99)/2), -1e-15)
252 264
253 ## Test input validation 265 ## Test input validation
254 %!error <Invalid call> betainc () 266 %!error <Invalid call> betainc ()
255 %!error <Invalid call> betainc (1) 267 %!error <Invalid call> betainc (1)
256 %!error <Invalid call> betainc (1,2) 268 %!error <Invalid call> betainc (1,2)