Mercurial > octave-nkf
diff scripts/statistics/distributions/betainv.m @ 20641:c3c052b9192a
Improve performance and error reporting of betainv, gaminv (bug #34363).
* betainv.m: Replace for loop with do-until loop. Shorten max loop cycles to
40, rather than 10,000. Issue warning if algorithm fails.
* gaminv.m: Replace for loop with do-until loop. Shorten max loop cycles to
40, rather than 100. Issue warning if algorithm fails.
author | Rik <rik@octave.org> |
---|---|
date | Sun, 11 Oct 2015 16:55:17 -0700 |
parents | d9341b422488 |
children | 4e307c55a2b5 |
line wrap: on
line diff
--- a/scripts/statistics/distributions/betainv.m Sun Oct 11 18:53:05 2015 -0400 +++ b/scripts/statistics/distributions/betainv.m Sun Oct 11 16:55:17 2015 -0700 @@ -81,8 +81,10 @@ y(l) = 1 - sqrt (myeps) * ones (length (l), 1); endif - y_old = y; - for i = 1 : 10000 + y_new = y; + loopcnt = 0; + do + y_old = y_new; h = (betacdf (y_old, a, b) - x) ./ betapdf (y_old, a, b); y_new = y_old - h; ind = find (y_new <= myeps); @@ -94,11 +96,11 @@ y_new(ind) = 1 - (1 - y_old(ind)) / 10; endif h = y_old - y_new; - if (max (abs (h)) < sqrt (myeps)) - break; - endif - y_old = y_new; - endfor + until (max (abs (h)) < sqrt (myeps) || ++loopcnt == 40) + + if (loopcnt == 40) + warning ("betainv: calculation failed to converge for some values"); + endif inv(k) = y_new; endif