[SciPy-dev] A couple problems of stats.poisson.pmf(k,mu)

[Patrick Marks]

Hi,
I've come across a couple problems in stats.poisson.pmf:

In [1]: from scipy.stats import *
In [2]: from scipy import special
In [3]: from numpy import asarray, exp

In [4]: poisson.pmf(10,10)
Out[4]: array(-0.026867175591182967)  # Negative PMF value!

In [5]: 10**10 * exp(-10) / asarray(special.gamma(11))
Out[5]: 0.1251100357211333

In [6]: poisson.pmf(10,10.0)
Out[6]: array(0.1251100357211333) # OK when mu is not an integer

In [7]: poisson.pmf(150,150.0)
Out[7]: array(1.#INF)   # Some sort of arithmetic overflow

The overflow could be avoided by computing the PMF in log space:
>poisson_gen._pmf(self, k, mu):
>   logPk = k * log(mu) - mu - arr(special.gammaln(k+1))
>    return exp(logPk)
which works to much higher numbers, and has a fractional deviation of
<10e-13 and typically == 0.0  from the current implementation in my
tests.

I'm not sure why the pmf is wrong when mu is an integer, although i
guess it's an easy fix.
Should i make a ticket in trac for this?

Patrick Marks
University of Illinois