gamma.pdf numerically unstable?
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? TAI
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?
What SciPy version? [~/] [1]: from scipy.stats import gamma [~/] [2]: gamma.pdf(1.6, 142, loc=0, scale=.0098) [2]: 0.69705038205764014 [~/] [3]: np.exp(gamma.logpdf(1.6, 142, loc=0, scale=.0098)) [3]: 0.69705038205763981 [~/] [4]: import scipy [~/] [5]: scipy.__version__ [5]: '0.12.0.dev-8647010' In R
dgamma(1.6, 142, scale=.0098) [1] 0.6970504
Skipper
On Mon, Aug 6, 2012 at 9:48 AM, Skipper Seabold <jsseabold@gmail.com> wrote:
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?
What SciPy version?
import scipy scipy.__version__ '0.9.0'
already fixed then
[~/] [1]: from scipy.stats import gamma
[~/] [2]: gamma.pdf(1.6, 142, loc=0, scale=.0098) [2]: 0.69705038205764014
[~/] [3]: np.exp(gamma.logpdf(1.6, 142, loc=0, scale=.0098)) [3]: 0.69705038205763981
[~/] [4]: import scipy
[~/] [5]: scipy.__version__ [5]: '0.12.0.dev-8647010'
In R
dgamma(1.6, 142, scale=.0098) [1] 0.6970504
Skipper _______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
Thanks for all the responses. My SciPy version is 0.7.2, so it looks like it might be a version issue. Will update and so what happens. On 8/6/12, josef.pktd@gmail.com <josef.pktd@gmail.com> wrote:
On Mon, Aug 6, 2012 at 9:48 AM, Skipper Seabold <jsseabold@gmail.com> wrote:
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?
What SciPy version?
import scipy scipy.__version__ '0.9.0'
already fixed then
[~/] [1]: from scipy.stats import gamma
[~/] [2]: gamma.pdf(1.6, 142, loc=0, scale=.0098) [2]: 0.69705038205764014
[~/] [3]: np.exp(gamma.logpdf(1.6, 142, loc=0, scale=.0098)) [3]: 0.69705038205763981
[~/] [4]: import scipy
[~/] [5]: scipy.__version__ [5]: '0.12.0.dev-8647010'
In R
dgamma(1.6, 142, scale=.0098) [1] 0.6970504
Skipper _______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
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
participants (3)
-
josef.pktd＠gmail.com -
Skipper Seabold -
Tom Furmston