Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP

Loading…

Fixing a special.gamma overrun in gamma.pdf #5

Closed
wants to merge 1 commit into from

3 participants

@pv
Owner
pv commented

Copied from scipy/scipy-svn#5


wesm opened this pull request February 21, 2011
Fixing a special.gamma overrun in gamma.pdf
There are more to be found in distributions.py-- I'll fix them as I find them.

@josef-pkt
Collaborator

(a-1)*log(x) - x - gamln(a) breaks if a=1 and x=0, corner case that is not ruled out.

the change in this pull request is good, however that it also propagates the 0log0 error from _logpdf to _pdf.

overall it looks like the change to _logpdf introduced several cases of possible 0*log(0) errors

@wesm

Likely no choice here but to check if x != 0? I'm not sure what's the best way to add a commit to this pull request

@josef-pkt
Collaborator

I would just pull this into scipy, and create separate pull requests or commits for the 0*log0 fixes

(Of all possible solution that we use for 0log0 I haven't found one I'm really happy with.)

@wesm

Why don't we do that and close this pull request, we should separately make a pass through and fix the 0*log(0) issues as you mention. @pv is that OK with you?

@pv
Owner
pv commented

Pulled.

Personally, I'd just make a separate ylogx function for y*log(x) with appropriately defined y -> 0 limit and use it everywhere where needed.

@pv pv closed this
@josef-pkt
Collaborator

Pauli,

having a correct xlogy function available would solve the problem and would be useful in many cases in scipy and statsmodels. I proposed it once on the numpy mailing list (Jan 14) but without reponse.

@Nikratio Nikratio referenced this pull request
Closed

Segfault in bisplrep #2446

@jnothman jnothman referenced this pull request from a commit
Commit has since been removed from the repository and is no longer available.
@jnothman jnothman referenced this pull request from a commit
Commit has since been removed from the repository and is no longer available.
@rgommers rgommers referenced this pull request from a commit
@rgommers rgommers BLD: fix MinGW Windows build issue by splitting linalg.interpolate so…
…urces.

Type of errors that this fixes::

    Argument #5 (named `ifac') of `zfftb1' is one type at (2) but is some
    other type at (1) [info -f g77 M GLOBALS]

Warning shows up with gfortran as::

    Warning: Type mismatch in argument 'ifac' at (1); passed REAL(8) to INTEGER(4)

Closes gh-2709.
72dc1d6
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Commits on Mar 10, 2011
  1. @wesm @pv

    BUG: fix overrun in gamma pdf

    wesm authored pv committed
This page is out of date. Refresh to see the latest.
View
2  scipy/stats/distributions.py
@@ -3061,7 +3061,7 @@ class gamma_gen(rv_continuous):
def _rvs(self, a):
return mtrand.standard_gamma(a, self._size)
def _pdf(self, x, a):
- return x**(a-1)*exp(-x)/special.gamma(a)
+ return exp(self._logpdf(x, a))
def _logpdf(self, x, a):
return (a-1)*log(x) - x - gamln(a)
def _cdf(self, x, a):
View
11 scipy/stats/tests/test_distributions.py
@@ -332,6 +332,17 @@ def test_cdf(self):
assert_almost_equal(stats.skellam.cdf(k, mu1, mu2), skcdfR, decimal=5)
+class TestGamma(TestCase):
+
+ def test_pdf(self):
+ # a few test cases to compare with R
+ pdf = stats.gamma.pdf(90, 394, scale=1./5)
+ assert_almost_equal(pdf, 0.002312341)
+
+ pdf = stats.gamma.pdf(3, 10, scale=1./5)
+ assert_almost_equal(pdf, 0.1620358)
+
+
class TestHypergeom(TestCase):
def test_precision(self):
# comparison number from mpmath
Something went wrong with that request. Please try again.