bugGNU Octave - Bugs: bug #39309, poisscdf fails for large arguments...


bug #39309: poisscdf fails for large arguments due to the algorithm in D9LGIT

Submitted by:  Thomas Weber <tweber>
Submitted on:  Fri 21 Jun 2013 11:28:33 AM UTC  
Category: LibrariesSeverity: 3 - Normal
Priority: 3 - LowItem Group: Inaccurate Result
Status: ConfirmedAssigned to: None
Originator Name: Open/Closed: Open
Release: devOperating System: GNU/Linux

Add a New Comment(Rich Markup)

You are not logged in

Please log in, so followups can be emailed to you.


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>
Project 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>
Project Administrator
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:

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:

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

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:

[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:


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.):

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:

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.

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.


Thomas Weber <tweber>


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

Attach File(s):

No files currently attached


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -unavailable- added by mtmiller (Posted a comment)
  • -unavailable- added by sebald (Posted a comment)
  • -unavailable- added by rik5 (Updated the item)
  • -unavailable- added by tweber (Submitted the item)

    Do you think this task is very important?
    If so, you can click here to add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.


    Please enter the title of George Orwell's famous dystopian book (it's a date):



    Follow 3 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Wed 16 Nov 2016 08:57:59 PM UTCmtmillerPriority5 - Normal=>3 - Low
    Fri 21 Jun 2013 11:57:06 PM UTCrik5StatusNone=>Confirmed

    Back to the top

    Powered by Savane 3.1-cleanup1