[Octave-bug-tracker] [bug #52569] unexptected error occured when using dynare with octave

[Dan Sebald]

Follow-up Comment #2, bug #52569 (project octave):

Thanks for the links.

I've looked into the code a bit, just to get a handle on what the issue is. 
SLATEC D9LGIT is a very short routine, so not difficult to comprehend.  SLATEC
is a collection of lesser libraries, of which d9lgit belongs to a library
called The Fullerton Function Library FN, which has been updated and ported to
C/C++ in some way but I didn't look at much of that:

https://people.sc.fsu.edu/~jburkardt/f77_src/fn/fn.html
https://people.sc.fsu.edu/~jburkardt/cpp_src/fn/fn.html
http://www.netlib.org/fn/d9lgit.f

The last reference is the actual, more up-to-date code.  However, it isn't
much different, aside from using the double versions of some of the routines. 
I'm attaching a diff file that makes Octave's d9lgit.f in line with the one in
the link above.  (Don't apply this to the repository, though.  I don't think
we want to start freely modifying third-party-sourced code.)  I don't see any
difference in convergence behavior when changing the ABS to DABS.

The issue is that poisscdf (A,A+E) doesn't converge well when A grows large
and E shrinks, e.g., poisscdf(1000,1001) versus poisscdf(1000,1002).  I don't
have time to think of the mathematical ramifications of this property (e.g.,
does this algorithm outright lose accuracy for the mathematical function as A
becomes large?).

That being the case, it appears that simply bumping up the allowable number of
iterations at least avoids the FORTRAN error.  That is, I made this change:


diff --git a/liboctave/external/slatec-fn/d9lgit.f
b/liboctave/external/slatec-f
--- a/liboctave/external/slatec-fn/d9lgit.f
+++ b/liboctave/external/slatec-fn/d9lgit.f
@@ -46,7 +46,7 @@ C
       R = 0.D0
       P = 1.D0
       S = P
-      DO 20 K=1,200
+      DO 20 K=1,500
         FK = K
         T = (A+FK)*X*(1.D0+R)
         R = T/((AX+FK)*(A1X+FK)-T)
@@ -55,7 +55,7 @@ C
         IF (ABS(P).LT.EPS*S) GO TO 30
  20   CONTINUE
       CALL XERMSG ('SLATEC', 'D9LGIT',
-     +   'NO CONVERGENCE IN 200 TERMS OF CONTINUED FRACTION', 3, 2)
+     +   'NO CONVERGENCE IN 500 TERMS OF CONTINUED FRACTION', 3, 2)
 C
  30   HSTAR = 1.0D0 - X*S/A1X
       IF (HSTAR .LT. SQEPS) CALL XERMSG ('SLATEC', 'D9LGIT',
(END)


and can then get results for larger A.  For example:


octave:5> poisscdf (4000,4000)
ans =  0.50421
octave:6> poisscdf (5000,5000)
 ***MESSAGE FROM ROUTINE D9LGIT IN LIBRARY SLATEC.
 ***FATAL ERROR, PROG ABORTED, TRACEBACK REQUESTED
 *  NO CONVERGENCE IN 500 TERMS OF CONTINUED FRACTION
 *  ERROR NUMBER = 3
 *   
 ***END OF MESSAGE
 
 ***JOB ABORT DUE TO FATAL ERROR.
0          ERROR MESSAGE SUMMARY
 LIBRARY    SUBROUTINE MESSAGE START             NERR     LEVEL     COUNT
 SLATEC     D9LGIT     NO CONVERGENCE IN 50         3         2         2

error: gammainc: exception encountered in Fortran subroutine xgammainc_
error: called from
    poisscdf at line 60 column 12
octave:6>


So, obviously this algorithm is going to fall apart at some point.  Perhaps
the original submitter could give us some insight as to reasonable
expectations with regard to applications.

This "error: Due to a bug in Octave" isn't quite a bug, it would seem.  It's
more like a limitation on numerical methods for approximating the function
which perhaps we could stretch out but never totally eliminate.

(file #42546)
    _______________________________________________________

Additional Item Attachment:

File name: octave-d9lgit_convergence-djs2017dec01.diff Size:1 KB


    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?52569>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/