Skip to content
This repository

BUG: fix precision issue with nbinom.pmf. Closes ticket 1779. #371

Merged
merged 1 commit into from over 1 year ago

3 participants

Ralf Gommers Josef Perktold Pauli Virtanen
Ralf Gommers
Owner

No description provided.

Josef Perktold
Collaborator

Thanks Ralf and bioinformed, Looks good to me

@pv a "special" question since you looked into similar issues recently
Do you know if it makes a difference for some relevant parameters?
betaln could replace 3 calls to gammaln, but I haven't seen any case yet in scipy.stats and in statsmodels where gammaln wasn't doing a good job
http://comments.gmane.org/gmane.comp.python.scientific.user/32968

Pauli Virtanen
Owner
pv commented December 02, 2012

@josef-pkt: if you want the binomial coefficient, you can use scipy.special.binom (will be in Scipy 0.12.0).

betaln is implemented using gammaln

Josef Perktold
Collaborator

Thanks Pauli, then I'd rather stick with gammaln, or whichever is closer to the formulas in the references.

Pauli Virtanen
Owner
pv commented December 02, 2012

@josef-pkt: the gammaln formula gives wrong results if |n| >> |x| or |x| >> |n|. PR #370 fixes this in the beta function.

Pauli Virtanen
Owner
pv commented December 02, 2012

But this PR looks good to me, too.

Ralf Gommers
Owner

Thanks for having a look. Merging.

Ralf Gommers rgommers merged commit 15e9ac0 into from December 02, 2012
Ralf Gommers rgommers closed this December 02, 2012
Clemens ClemensFMN referenced this pull request from a commit August 06, 2013
Commit has since been removed from the repository and is no longer available.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Showing 1 unique commit by 1 author.

Dec 02, 2012
Ralf Gommers BUG: fix precision issue with nbinom.pmf. Closes ticket 1779. 9fc353a
This page is out of date. Refresh to see the latest.
3  scipy/stats/distributions.py
@@ -6529,8 +6529,7 @@ def _rvs(self, n, p):
6529 6529
     def _argcheck(self, n, p):
6530 6530
         return (n >= 0) & (p >= 0) & (p <= 1)
6531 6531
     def _pmf(self, x, n, p):
6532  
-        coeff = exp(gamln(n+x) - gamln(x+1) - gamln(n))
6533  
-        return coeff * power(p,n) * power(1-p,x)
  6532
+        return exp(self._logpmf(x, n, p))
6534 6533
     def _logpmf(self, x, n, p):
6535 6534
         coeff = gamln(n+x) - gamln(x+1) - gamln(n)
6536 6535
         return coeff + n*log(p) + x*log(1-p)
6  scipy/stats/tests/test_distributions.py
@@ -160,6 +160,12 @@ def test_rvs(self):
160 160
         assert_(isinstance(val, numpy.ndarray))
161 161
         assert_(val.dtype.char in typecodes['AllInteger'])
162 162
 
  163
+    def test_pmf(self):
  164
+        # regression test for ticket 1779
  165
+        assert_allclose(np.exp(stats.nbinom.logpmf(700, 721, 0.52)),
  166
+                        stats.nbinom.pmf(700, 721, 0.52))
  167
+
  168
+
163 169
 class TestGeom(TestCase):
164 170
     def test_rvs(self):
165 171
         vals = stats.geom.rvs(0.75, size=(2, 50))
Commit_comment_tip

Tip: You can add notes to lines in a file. Hover to the left of a line to make a note

Something went wrong with that request. Please try again.