[Scipy-tickets] [SciPy] #1893: Binomial PMF returns incorrect answer
SciPy Trac
scipy-tickets@scipy....
Wed Apr 17 20:22:47 CDT 2013
#1893: Binomial PMF returns incorrect answer
-------------------------+--------------------------------------------------
Reporter: uscjeremy | Owner: rgommers
Type: defect | Status: new
Priority: normal | Milestone: Unscheduled
Component: scipy.stats | Version: 0.12.0
Keywords: |
-------------------------+--------------------------------------------------
Comment(by uscjeremy):
Digging in a little further, it appears that this bug is because when k=n
and p=1, one term of the PMF formula reduces to 0**0, causing an
exception. This is a case where the math doesn't comport with the logical
meaning of the binomial distribution, though; if the probability of each
trial succeeding is 1, then the probability of n trials yielding n
successes is 1 and the probability of yielding any other number of
successes is 0.
Mathematica and R both give 1 as the answer for binom.pmf(1,1,1) where
scipy gives nan. It looks like R's implementation has a bunch of special
code to deal with the edge cases where p=0 or 1 and k=0 or n, which I
reproduce below (R v3.0.0, src/nmath/dbinom.c). Note they use 'x' where
scipy uses 'k', and q=1-p.
{{{
if (p == 0) return((x == 0) ? R_D__1 : R_D__0);
if (q == 0) return((x == n) ? R_D__1 : R_D__0);
if (x == 0) {
if(n == 0) return R_D__1;
lc = (p < 0.1) ? -bd0(n,n*q) - n*p : n*log(q);
return( R_D_exp(lc) );
}
if (x == n) {
lc = (q < 0.1) ? -bd0(n,n*p) - n*q : n*log(p);
return( R_D_exp(lc) );
}
if (x < 0 || x > n) return( R_D__0 );
}}}
--
Ticket URL: <http://projects.scipy.org/scipy/ticket/1893#comment:1>
SciPy <http://www.scipy.org>
SciPy is open-source software for mathematics, science, and engineering.
More information about the Scipy-tickets
mailing list