bugGNU Octave - Bugs: bug #39309, [octave forge] (statistics)...

 
 

bug #39309: [octave forge] (statistics) poisscdf fails for large arguments due to the algorithm in D9LGIT

Submitter:  Thomas Weber <tweber>
Submitted:  Fri 21 Jun 2013 11:28:33 AM UTC
   
 
Category:  Octave Package Severity:  3 - Normal
Priority:  3 - Low Item Group:  Inaccurate Result
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Tue 26 Jul 2022 06:35:26 PM UTC, comment #5: 

looking at the current Statistics 1.4.3 with Octave 7.1.0-

poisscdf is largely just a call to gammainc, which is in core octave, and has undergone at least a few revisions in the last decade.

the comment #0, comment #1, comment #2 codes no longer produce fortran errors:


comment #0:

>> poisscdf (604,604)
ans = 0.510819531737768
>> poisscdf (605,605)
ans = 0.510810590057861
>> poisscdf (4241,4241)
ans = 0.504083866960345
>> poisscdf (4242,4242)
ans = 0.504083385599947


comment #1:

>> poisscdf (100000,100000)
ans = 0.500841043099343
>> poisscdf (100001,100001)
ans = 0.500841038894170
>> poisscdf (1e10,1e10+1e5)
ans = 0.158657673645873


comment #2

>> poisscdf (10,10)
ans = 0.583039750192986
>> poisscdf (100,100)
ans = 0.526562198529999
>> poisscdf (1000,1000)
ans = 0.508409367168506
>> poisscdf (100000,100000)
ans = 0.500841043099343
>> poisscdf (100000000,100000000)
ans = 0.500026596152138
>> poisscdf (100000000000,100000000000)
ans = 0.500000841049293
>> poisscdf (100000000000000,100000000000000)
ans = 0.500000026770566
>> poisscdf (100000000000000000,100000000000000000)
ans = 84167.84956743612


soooo....  it seems that the original issue has largely been fixed.

the last test did go way too fast, so I assume it did hit some internal convergence error, although there were no fortran errors/warnings.  that part is a bit concerning, as I would have hoped some non-convergence warning might have been issued.  but again, that would seem to be an issue with gammainc and not poisscdf. so I'm going to close this old bug report.

Nicholas Jankowski <nrjank>
Group Member
Wed 16 Nov 2016 10:28:49 PM UTC, comment #4: 

Before the next major release (4.4.0) I want to trim down Octave's core.  For historical reasons, we have many more functions in the core than Matlab and many of these are only of interest to a narrow group of users.  These functions should be shifted to the appropriate packages on Octave-Forge so that users can load on demand these infrequent functions.

Among the functions to move out of core are the extensive number of statistical functions like poisscdf.  There are also probably too many "signal" functions which should be moved to the signal package.

Rik <rik5>
Group administrator
Wed 16 Nov 2016 08:57:59 PM UTC, comment #3: 

Just updating and re-confirming that this bug is still present in the development version (4.3.0+), no changes have been made to the underlying Fortran routines.

Mike Miller <mtmiller>
Group Member
Sat 22 Jun 2013 08:08:08 AM UTC, comment #2: 

I increased the limit on the number of iterations to a very large number, 1000000000.  That allows using larger numbers as inputs:


>> poisscdf (10,10)
ans =  0.58304
>> poisscdf (100,100)
ans =  0.52656
>> poisscdf (100000,100000)
ans =  0.50084
>> poisscdf (100000000,100000000)
ans =  0.50003
>> poisscdf (100000000000,100000000000)
ans =  0.49986
>> poisscdf (100000000000000,100000000000000)
ans =  0.39347
>> poisscdf (100000000000000000,100000000000000000)
 ***MESSAGE FROM ROUTINE D9LGIT IN LIBRARY SLATEC.
 ***POTENTIALLY RECOVERABLE ERROR, PROG ABORTED, TRACEBACK REQUESTED
 *  RESULT LESS THAN HALF PRECISION
 *  ERROR NUMBER = 1
 *
 ***END OF MESSAGE

 ***JOB ABORT DUE TO UNRECOVERED ERROR.
0          ERROR MESSAGE SUMMARY
 LIBRARY    SUBROUTINE MESSAGE START             NERR LEVEL     COUNT
 SLATEC     D9LGIT     RESULT LESS THAN HAL         1     1         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


but notice that a different type of error is beginning to surface.  Also, near that point of failure the answer is starting to move away from 0.5.

ans =  0.39347

Without doing a mathematical convergence proof, my gut feeling is that the CDF for a Poisson r.v. when k = lambda should tend to 0.5 for large lambda.  The Poisson r.v. begins to approximate a Gaussian r.v. with mean lambda and variance lambda.

So, at some point for large lambda this routine is experiencing numerical problems.  Again, it isn't surprising because of the way gamma is defined.  If one were to use the definition of CDF (accumulating PDF values), the probability would be spread so thin that I'm guessing that too would have numerical problems as well.

I think what we can do is increase the limit on that DO loop some, but at some point the user is going to have to realize there has to be reasonable bound on those input values or switch over to an approximation such as a Gaussian r.v., or scale their problem somehow.

Dan Sebald <sebald>
Sat 22 Jun 2013 04:03:42 AM UTC, comment #1: 

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.

Dan Sebald <sebald>
Fri 21 Jun 2013 11:28:33 AM UTC, original submission:  

Hi,
this is an old bug report from Debian that I forward here in case someone wants to spend some time on Fortran or GSL. Please note that John already commented in the original bug report years ago.
http://bugs.debian.org/cgi-bin/bugreport.cgi?bug=418158

So, the issue is that
poisscdf (604,604)
works, but
poisscdf (605,605)
fails with a Fortran error, as the algorithm did not converge. This is ultimately due to the algorithm in libcruft/slatec-fn/d9lgit.f.

There was the suggestion to use the gamma_inc_P function from the gsl package, but as far as I can see, this was never implemented.

Regards
Thomas

Thomas Weber <tweber>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by nrjank (Posted a comment)
  • -email is unavailable- added by sebald (Posted a comment)
  • -email is unavailable- added by tweber (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 10 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-07-26 nrjank StatusConfirmed Fixed
        Open/ClosedOpen Closed
    2019-06-04 mtmiller Summary[Octave-Forge] poisscdf fails for large arguments due to the algorithm in D9LGIT [octave forge] (statistics) poisscdf fails for large arguments due to the algorithm in D9LGIT
    2019-06-04 mtmiller Carbon-CopyRemoved 80942 -
    2019-06-04 rik5 Carbon-CopyRemoved 72865 -
    2019-06-04 rik5 CategoryLibraries Octave Package
        Summarypoisscdf fails for large arguments due to the algorithm in D9LGIT [Octave-Forge] poisscdf fails for large arguments due to the algorithm in D9LGIT
    2016-11-16 mtmiller Priority5 - Normal 3 - Low
        Release3.6.4 dev
    2013-06-21 rik5 StatusNone Confirmed

    Back to the top

    Powered by Savane 3.13-3230.
    Corresponding source code