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

[Dan Sebald]

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

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.

    _______________________________________________________

Reply to this item at:

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

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