Mercurial > octave-nkf

changeset 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 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	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
--- a/scripts/statistics/distributions/gaminv.m	Sun Oct 11 18:53:05 2015 -0400
+++ b/scripts/statistics/distributions/gaminv.m	Sun Oct 11 16:55:17 2015 -0700
@@ -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