# [Scipy-tickets] [SciPy] #620: scipy.stats.distributions.binom.pmf returns incorrect values

SciPy scipy-tickets@scipy....
Fri Dec 5 11:10:02 CST 2008

```#620: scipy.stats.distributions.binom.pmf returns incorrect values
----------------------------+-----------------------------------------------
Reporter:  robertbjornson  |        Owner:  josefpktd
Type:  defect          |       Status:  new
Priority:  normal          |    Milestone:  0.8
Component:  scipy.stats     |      Version:
Severity:  normal          |   Resolution:
Keywords:                  |
----------------------------+-----------------------------------------------
Changes (by josefpktd):

* owner:  somebody => josefpktd
* milestone:  0.7.0 => 0.8

Comment:

I found a related discussion about the hypergeometric distribution:

Robert Kern proposed to use gammaln, which seems to work also for the
binomial very well.

Accuracy looks pretty good, comparable to using scipy.comb for small
values (run script below)

{{{
'''
implementing binomial pdf with gammaln

see: Robert Kerns discussion in

'''
import numpy as np
from scipy import special, comb

def bpmf(k,n, p):
#this does not work for large n
return comb(n, k)*(p**k)*((1-p)**(n-k))

# proposed version using gammaln
def combinln(n,k):
return (special.gammaln(n+1) - (special.gammaln(k+1) +
special.gammaln(n-k+1)))

def bpmfln(k,n, p):
return np.exp(combinln(n, k) + k*np.log(p) + (n-k)*np.log(1-p))

n=10;
p=1.0e-5;

print "using gammaln"
print bpmfln(np.arange(11),10, p)
print "using comb"
print bpmf(np.arange(11),10, p)
print "difference"
print bpmfln(np.arange(11),10, p)- bpmf(np.arange(11),10, p)

pmfnln = bpmfln(np.arange(5001),5000, 0.99)
print 'n = 5000'
print 'nans', np.sum(np.isnan(pmfnln))
print 'sum', np.sum(pmfnln)
print 'sum (repr)', repr(np.sum(pmfnln))

}}}

--
Ticket URL: <http://scipy.org/scipy/scipy/ticket/620#comment:5>
SciPy <http://www.scipy.org/>
SciPy is open-source software for mathematics, science, and engineering.
```