Mercurial > octave-nkf

comparison scripts/statistics/distributions/gaminv.m @ 20641:c3c052b9192a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
comparison
equal deleted inserted replaced
20640:f9c991dc5c1a 20641:c3c052b9192a
77 l = find (x < myeps); 77 l = find (x < myeps);
78 if (any (l)) 78 if (any (l))
79 y(l) = sqrt (myeps) * ones (length (l), 1); 79 y(l) = sqrt (myeps) * ones (length (l), 1);
80 endif 80 endif
81 81
82 y_old = y; 82 y_new = y;
83 for i = 1 : 100 83 loopcnt = 0;
84 do
85 y_old = y_new;
84 h = (gamcdf (y_old, a, b) - x) ./ gampdf (y_old, a, b); 86 h = (gamcdf (y_old, a, b) - x) ./ gampdf (y_old, a, b);
85 y_new = y_old - h; 87 y_new = y_old - h;
86 ind = find (y_new <= myeps); 88 ind = find (y_new <= myeps);
87 if (any (ind)) 89 if (any (ind))
88 y_new(ind) = y_old(ind) / 10; 90 y_new(ind) = y_old(ind) / 10;
89 h = y_old - y_new; 91 h = y_old - y_new;
90 endif 92 endif
91 if (max (abs (h)) < sqrt (myeps)) 93 until (max (abs (h)) < sqrt (myeps) || ++loopcnt == 40)
92 break; 94
93 endif 95 if (loopcnt == 40)
94 y_old = y_new; 96 warning ("gaminv: calculation failed to converge for some values");
95 endfor 97 endif
96 98
97 inv(k) = y_new; 99 inv(k) = y_new;
100
98 endif 101 endif
99 102
100 endfunction 103 endfunction
101 104
102 105