[SciPy-dev] A couple problems of stats.poisson.pmf(k,mu)
Patrick Marks
pmarks at gmail.com
Mon Jul 3 14:24:26 CDT 2006
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
More information about the Scipy-dev
mailing list