Mercurial > hg > octave-nkf
changeset 20817: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 | f9c991dc5c1a |
children | 9d2023d1a63c |
files | scripts/statistics/distributions/betainv.m scripts/statistics/distributions/gaminv.m |
diffstat | 2 files changed, 19 insertions(+), 14 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/statistics/distributions/betainv.m +++ b/scripts/statistics/distributions/betainv.m @@ -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
--- a/scripts/statistics/distributions/gaminv.m +++ b/scripts/statistics/distributions/gaminv.m @@ -79,8 +79,10 @@ y(l) = sqrt (myeps) * ones (length (l), 1); endif - y_old = y; - for i = 1 : 100 + y_new = y; + loopcnt = 0; + do + y_old = y_new; h = (gamcdf (y_old, a, b) - x) ./ gampdf (y_old, a, b); y_new = y_old - h; ind = find (y_new <= myeps); @@ -88,13 +90,14 @@ y_new(ind) = y_old(ind) / 10; h = y_old - y_new; endif - 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 ("gaminv: calculation failed to converge for some values"); + endif inv(k) = y_new; + endif endfunction