[SciPy-User] gamma.pdf numerically unstable?

josef.pktd@gmai... josef.pktd@gmai...
Mon Aug 6 09:04:12 CDT 2012


On Mon, Aug 6, 2012 at 9:29 AM, Tom Furmston <tfurmston@googlemail.com>wrote:

> Hi,
>
> I'm trying to use the gamma distribution in the SciPy library, but it
> appears to be numerically unstable.
>
> First, I've imported the library with
>
> from scipy.stats import gamma
>
> I want to evaluate the gamma distribution are x = 1.6, where the gamma
> distribution has the parameters shape = 142, scale = 0.0098. When I
> type
>
> gamma.pdf(1.6, 142, loc=0, scale=0.0098)
>
> I get a value of inf.
>
> According to matlab the value should be around 0.7030.
>
> Am I using this distribution correctly? Is this distribution known to
> suffer from numerical issues?
>

looks like numerical instability in the pdf, I guess it should use
exp(logpdf) instead

>>> stats.gamma.pdf(1.6, 142, loc=0, scale=0.0098)
inf
>>> np.exp(stats.gamma.logpdf(1.6, 142, loc=0, scale=0.0098))
0.69705038205771896

would be worth a check with a multiprecision library.
>>> x, a = 1.6, 142; loc=0; scale=0.0098
>>> mpmath.power(x/scale, (a-1))*mpmath.exp(-x/scale)/mpmath.gamma(a) /
scale
mpf('0.6970503820577183')

a bit lower than matlab

>>> np.exp(stats.gamma.logpdf(1.6, 142, loc=0, scale=0.098))
4.5494987629369171e-79
>>> stats.gamma.pdf(1.6, 142, loc=0, scale=0.098)
4.5494987629366552e-79

It's not a known numerical problem, some numerical problems only surface
when someone tries unusual parameters.

Josef


>
> TAI
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20120806/73ee9a6b/attachment.html 


More information about the SciPy-User mailing list