[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:

 http://groups.google.ca/group/comp.lang.python/browse_thread/thread/839009574397dc37

 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
 http://groups.google.ca/group/comp.lang.python/browse_thread/thread/839009574397dc37

 '''
 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.


More information about the Scipy-tickets mailing list