Mercurial > octave-nkf

comparison scripts/statistics/distributions/betainv.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
79 l = find (y > 1 - myeps); 79 l = find (y > 1 - myeps);
80 if (any (l)) 80 if (any (l))
81 y(l) = 1 - sqrt (myeps) * ones (length (l), 1); 81 y(l) = 1 - sqrt (myeps) * ones (length (l), 1);
82 endif 82 endif
83 83
84 y_old = y; 84 y_new = y;
85 for i = 1 : 10000 85 loopcnt = 0;
86 do
87 y_old = y_new;
86 h = (betacdf (y_old, a, b) - x) ./ betapdf (y_old, a, b); 88 h = (betacdf (y_old, a, b) - x) ./ betapdf (y_old, a, b);
87 y_new = y_old - h; 89 y_new = y_old - h;
88 ind = find (y_new <= myeps); 90 ind = find (y_new <= myeps);
89 if (any (ind)) 91 if (any (ind))
90 y_new(ind) = y_old(ind) / 10; 92 y_new(ind) = y_old(ind) / 10;
92 ind = find (y_new >= 1 - myeps); 94 ind = find (y_new >= 1 - myeps);
93 if (any (ind)) 95 if (any (ind))
94 y_new(ind) = 1 - (1 - y_old(ind)) / 10; 96 y_new(ind) = 1 - (1 - y_old(ind)) / 10;
95 endif 97 endif
96 h = y_old - y_new; 98 h = y_old - y_new;
97 if (max (abs (h)) < sqrt (myeps)) 99 until (max (abs (h)) < sqrt (myeps) || ++loopcnt == 40)
98 break; 100
99 endif 101 if (loopcnt == 40)
100 y_old = y_new; 102 warning ("betainv: calculation failed to converge for some values");
101 endfor 103 endif
102 104
103 inv(k) = y_new; 105 inv(k) = y_new;
104 endif 106 endif
105 107
106 endfunction 108 endfunction