Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP

Loading…

Fix for Ticket #1842: cornercase: scipy.stats.binom.pfm(100,100,1) #461

Closed
wants to merge 6 commits into from

3 participants

Per A. Brodtkorb Pauli Virtanen Josef Perktold
Per A. Brodtkorb

This PR fixes the bug reported in Ticket #1842 ( http://projects.scipy.org/scipy/ticket/1842 ):

wrong value returned by scipy.stats.binom probability mass function when true probability is 1.

Also added a small test for this

scipy/stats/distributions.py
@@ -6599,7 +6604,12 @@ def _logpmf(self, x, n, p):
k = floor(x)
combiln = (gamln(n+1) - (gamln(k+1) +
gamln(n-k+1)))
- return combiln + k*np.log(p) + (n-k)*np.log(1-p)
+ logp = where((p==0) & (k==0), 1, log(p))
+ log1mp = where((p==1) & (k==n), 1, log1p(-p))
+ klogp_nmklog1mp = where(k==0, n * log1mp,
+ where(k==n, k * logp,
+ k*logp + (n-k)*log1mp))
Josef Perktold Collaborator

wouldn't it be easier to create the klogp_nmklog1mp first and then conditionally assign to it, instead of the many and nested where

Per A. Brodtkorb
pbrod added a note

You are right the nested "where" statements are not necessary and can be replaced with klogp + (n-k)log1mp directly. I will fix that.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Pauli Virtanen
Owner
pv commented

Thanks, merged in 17edfda --- but without the log1p redefinition

Pauli Virtanen pv closed this
Clemens ClemensFMN referenced this pull request from a commit
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
Commits on Mar 9, 2013
  1. Per A. Brodtkorb

    Fix for Ticket #1842 : wrong value returned by scipy.stats.binom prob…

    pbrod authored
    …ability mass function when true probability is 1
    
    + a test for it.
  2. Per A. Brodtkorb

    Removed try except

    pbrod authored
  3. Per A. Brodtkorb
Commits on Mar 10, 2013
  1. Per A. Brodtkorb
Commits on Mar 16, 2013
  1. Per A. Brodtkorb
  2. Per A. Brodtkorb
This page is out of date. Refresh to see the latest.
17 scipy/stats/distributions.py
View
@@ -25,14 +25,16 @@
from numpy import atleast_1d, polyval, ceil, place, extract, \
any, argsort, argmax, vectorize, r_, asarray, nan, inf, pi, isinf, \
NINF, empty
+_log1p = log1p
import numpy
import numpy as np
import numpy.random as mtrand
from numpy import flatnonzero as nonzero
+
from . import vonmises_cython
from ._tukeylambda_stats import tukeylambda_variance as _tlvar, \
- tukeylambda_kurtosis as _tlkurt
-
+ tukeylambda_kurtosis as _tlkurt
+
__all__ = [
'rv_continuous',
'ksone', 'kstwobign', 'norm', 'alpha', 'anglit', 'arcsine',
@@ -71,7 +73,11 @@
# Python 3
def instancemethod(func, obj, cls):
return types.MethodType(func, obj)
-
+
+def log1p(x):
+ '''avoids warnings for x==-1'''
+ mx = where(x==-1, 0, x)
+ return where(x==-1, -inf, _log1p(mx))
# These are the docstring parts used for substitution in specific
# distribution docstrings.
@@ -6599,7 +6605,7 @@ def _logpmf(self, x, n, p):
k = floor(x)
combiln = (gamln(n+1) - (gamln(k+1) +
gamln(n-k+1)))
- return combiln + k*np.log(p) + (n-k)*np.log(1-p)
+ return combiln + special.xlogy(k,p) + special.xlog1py(n-k, -p)
def _pmf(self, x, n, p):
return exp(self._logpmf(x, n, p))
def _cdf(self, x, n, p):
@@ -7282,3 +7288,6 @@ def _stats(self, mu1, mu2):
return mean, var, g1, g2
skellam = skellam_gen(a=-np.inf, name="skellam", longname='A Skellam',
shapes="mu1,mu2")
+
+
+
8 scipy/stats/tests/test_distributions.py
View
@@ -136,7 +136,13 @@ def test_rvs(self):
val = stats.binom(10, 0.75).rvs(3)
assert_(isinstance(val, numpy.ndarray))
assert_(val.dtype.char in typecodes['AllInteger'])
-
+
+ def test_pmf(self):
+ # regression test for Ticket #1842
+ vals1 = stats.binom.pmf(100, 100,1)
+ vals2 = stats.binom.pmf(0, 100,0)
+ assert_allclose(vals1, 1.0, rtol=1e-15, atol=0)
+ assert_allclose(vals2, 1.0, rtol=1e-15, atol=0)
class TestBernoulli(TestCase):
def test_rvs(self):
Something went wrong with that request. Please try again.