[Octave-bug-tracker] [bug #39309] poisscdf fails for large arguments due to the algorithm in D9LGIT

[Dan Sebald]

Follow-up Comment #1, bug #39309 (project octave):

Just as an initial investigation, I increased the allowable number of
iterations of the core loop in the following way:


--- a/liboctave/cruft/slatec-fn/d9lgit.f
+++ b/liboctave/cruft/slatec-fn/d9lgit.f
@@ -46,7 +46,7 @@
       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)


and then found that the largest value for similar integer inputs is 4241,
i.e.:


>> poisscdf (4241,4241)
ans =  0.50408
>> poisscdf (4242,4242)
 ***MESSAGE FROM ROUTINE D9LGIT IN LIBRARY SLATEC.
 ***FATAL ERROR, PROG ABORTED, TRACEBACK REQUESTED
 *  NO CONVERGENCE IN 200 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 20         3         2  4

error: gammainc: exception encountered in Fortran subroutine xgammainc_
error: called from:
error:  
/usr/local/src/octave/octave/build-gui-11/../octave/scripts/statistics/distributions/poisscdf.m
at line 61, column 12
>>


Convergence obviously slows down with increasing magnitude of input, but
running "poisscdf (4241,4241)" doesn't seem particularly longer.  I'd say
there might be an artificial limit on the number of iterations.

We could increase the allowable number of iterations.  Typically, whenever in
numerical routines I've experienced non-convergence, the algorithm runs for a
few seconds.  So, I see no problems with extending the limit to be much higher
if the algorithm doesn't take but a fraction of a second.  That is sort of the
simplest approach without having to do much analysis.  I just increased the
limit to 5000 and can get up to:


>> poisscdf (100000,100000)
ans =  0.50084


[We can see this is converging to 0.5, probably converging to some other
distribution as argument tends to infinity.]

On the other hand, a more analytic approach would be to recognize some
characteristics for the input in which the computation should be translated
into a different problem which has faster convergence.

BUT, before doing any of that, I'm wondering if slatec:

http://www.netlib.org/slatec/index.html

is actively maintained.  It doesn't appear to be.  If a patch is made to the
routine, it should probably be renamed to avoid confusion, or maybe even
written is C.

One last thing is that the same problem exists in d9lgic.f, which apparently
is ultimately called if the second input is larger than the first.  (The "c"
means complimentary.):


>> poisscdf (1e10,1e10+1e5)
 ***MESSAGE FROM ROUTINE D9LGIC IN LIBRARY SLATEC.
 ***FATAL ERROR, PROG ABORTED, TRACEBACK REQUESTED
 *  NO CONVERGENCE IN 300 TERMS OF CONTINUED FRACTION
 *  ERROR NUMBER = 1
 *
 ***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 20         3         2 31
 SLATEC     XGMAINC    RESULT LT HALF PRECI         1         1  7
 SLATEC     D9LGIC     NO CONVERGENCE IN 30         1         2  1

error: gammainc: exception encountered in Fortran subroutine xgammainc_
error: called from:
error:  
/usr/local/src/octave/octave/build-gui-11/../octave/scripts/statistics/distributions/poisscdf.m
at line 61, column 12


I will think this over a bit.  I'm staring at the integral equation for the
incomplete gamma function at the moment and don't see an easy way to translate
for very large numbers.  I can see why the equation becomes so sensitive.
t^(largenumber - 1) grows very fast.  There is the normal gamma recursion
translation, but that still doesn't get one very far.


    _______________________________________________________

Reply to this item at:

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

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