special: handle the edge case lambda=0 in pdtr, pdtrc and pdtrik

[WarrenWeckesser]

This PR makes a few changes so that special.pdtr, special.pdtrc and special.pdtrik give the correct limiting value when the Poisson parameter lambda is 0.

In the case of pdtrik, the function now immediately returns 0 for lambda < 0.01 and p < 0.975. This is a lot more than the simple edge case lambda = 0, but it is based on the following figure showing the results of computing pdtrik from master on a grid of values:

The gray dots are where the result is 0. The red dots are where pdtrik returns 1e+100; the "correct" value is 0 at those points.

The script to generate the plot is in a gist: https://gist.github.com/WarrenWeckesser/7995139

[coveralls]

Coverage remained the same when pulling bf715ac on WarrenWeckesser:special-poisson-edge-cases into 34ae412 on scipy:master.

[WarrenWeckesser]

I haven't been able to reproduce the failure that is occurring in travis with python 2.7. With the changes in the PR, cephes.pdtrik([[0], [0.25], [0.95]], [0, 1e-20, 1e-6]) should return all zeros, but on travis, it returns all nans. Anyone have a suggestion for how to track this down?

[rgommers]

I can reproduce this on 32-bit linux, will have a look.

[rgommers]

This should be status = 0, that should fix the failure.

[WarrenWeckesser]

Ah. I must have looked at that line dozens of times and never noticed the typo.

Adding IMPLICIT NONE would have caught it. Apparently IMPLICIT NONE is not standard Fortran 77, but we're already using it in special/specfun/specfun.f. Any objections to adding it to this subroutine?

[rgommers]

No opinion, I don't really speak Fortran.

[ev-br]

Unless implicit none breaks in some of our more esoteric fortan 77 setups, please please please add it!

[ev-br]

I personally am more wary of Fortran code without implicit none --- standard or not, as long as it actually compiles.

[ev-br]

gh-3168

[rgommers]

Still strange that the failure depends on Python version (or even the way I call it; within IPython it doesn't fail). So something else may still be wrong.

[WarrenWeckesser]

@rgommers wrote:

Still strange that the failure depends on Python version (or even the way I call it; within IPython it doesn't fail). So something else may still be wrong.

Because of the typo, the variable was never assigned in the Fortran code. The variable was originally declared on the stack in the C wrapper, so it was uninitialized. If the value happened to be negative, 3 or 4, the wrapper function returned nan. So it is not surprising that it would fail in random ways. (Of course, I may be overconfident here--the tests on Travis for the updated code haven't finished yet. :)

[WarrenWeckesser]

We already use "implicit none" quite a bit (search for "implicit none"), so apparently it is not a problem, even though it is not standard Fortran 77. I added it to the cdfpoi subroutine in the fixed version of the PR.

[coveralls]

Coverage remained the same when pulling f2d32d0 on WarrenWeckesser:special-poisson-edge-cases into 34ae412 on scipy:master.

[coveralls]

Coverage remained the same when pulling 6e61391 on WarrenWeckesser:special-poisson-edge-cases into 34ae412 on scipy:master.

[WarrenWeckesser]

All the tests pass now. Yesterday the Python 3.3 build had a timeout error.

[WarrenWeckesser]

@pv: As the resident special specialist, could you take a look at this when you get a chance?

[ev-br]

A weird comment very very late: do we actually need to have pdtrik returning a float? The only use of it I can find seems to be in https://github.com/scipy/scipy/blob/master/scipy/stats/_discrete_distns.py#L446
where it's being cast to a whole number anyway.

[josef-pkt]

do we actually need to have pdtrik returning a float?
I know that the Poisson log-likelihood is popular also for estimating non-negative continuous data.
But I don't see what a continuous extension of pdtrik should mean.

In some cases it's useful to have a continuous extension for calculations, even if we convert to integers at the end. Maybe it would be easy to generate Poisson random numbers with a continuous pdtrik return.

[pv]

LGTM. Some checking shows that pdtr is a bit buggy for large lambda, however (underflows):

In [75]: p = np.r_[0, np.logspace(-300, 0, 2000)]
In [76]: m = np.r_[0, np.logspace(-300, 5, 2000)]
In [77]: v = (pdtr(np.ceil(pdtrik(p[:,None], m)), m) < p[:,None])
In [78]: np.where(v)
(array([  1,   1,   1, ..., 927, 929, 929]),
 array([1987, 1988, 1989, ..., 1995, 1997, 1998]))
In [80]: k = 123; pdtr(np.ceil(pdtrik(p[i[k]], m[j[k]])), m[j[k]]), p[i[k]], m[j[k]], np.ceil(pdtrik(p[i[k]], m[j[k]]))
(0.0, 1.5870854134187277e-299, 17263.090976752374, 9771.0)