import matplotlib.pyplot as plt
%matplotlib inline
from collections import Counter
import numpy as np
from scipy.misc import factorial
rand = np.random.rand
def Poisson(m,z=1): return np.exp(-z)*z**m/factorial(m)
This simulates throwing nballs
balls into nbins
bins, with average z = nballs/nbins
:
nbins=1000
nballs=1000
z=nballs/float(nbins)
balls_in_bin=Counter(map(int,nbins*rand(nballs)))
c=Counter(balls_in_bin.values())
c[0]= nbins-sum(c.values()) #bins with no balls
m=np.arange(max(c.keys())+1) #an array of values [0,1,2,...]
p=Poisson(m,z) #so also an array of values
yerr=np.sqrt(nbins*p*(1-p)) #and also an array of values
plt.bar(m,[c[m0] for m0 in m],color='g',align='center',alpha=.5)
plt.errorbar(m,nbins*p,yerr,fmt='o')
plt.xlim(-.4,m[-1]+.5)
plt.ylim(0,)
plt.xlabel('m')
plt.ylabel('#bins with m balls')
plt.title('{} balls in {} bins, mean={}'.format(nballs,nbins,z));
We can also divide everything by nbins to get the probabilities (note that the stderr now has nbins in the denominator under the squareroot):
nbins=1000
nballs=2000
z=nballs/float(nbins)
balls_in_bin=Counter(map(int,nbins*rand(nballs)))
c=Counter(balls_in_bin.values())
c[0]=nbins-sum(c.values()) #bins with no balls
m=np.arange(max(c.keys())+1) #an array of values [0,1,2,...]
p=Poisson(m,z) #so also an array of values
yerr=np.sqrt(p*(1-p)/nbins) #and also an array of values
plt.bar(m,[c[m0]/float(nbins) for m0 in m],color='g',align='center',alpha=.5)
plt.errorbar(m,p,yerr,fmt='o')
plt.xlim(-.4,m[-1]+.5)
plt.ylim(0,)
plt.xlabel('m')
plt.ylabel('probability for bin to have m balls')
plt.title('{} balls in {} bins, mean={}'.format(nballs,nbins,z));
Note that the above Poisson()
function defined above is also available as a standard library function poisson.pmf()
from scipy.stats import poisson
(see http://docs.scipy.org/doc/scipy-0.17.0/reference/generated/scipy.stats.poisson.html
for various other poisson()
methods such as .rvs() for drawing random samples, .cdf() for cumulative distribution function, and so on)