[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.
```