import numpy as np
import matplotlib.pyplot as plt
import math as m

hist=np.zeros(10)
poisson=np.zeros(1000)
s=0
for l in range(1000):
 r=np.random.rand(100)
 for k in range(10):
   for i in range(100):
     if (k/10)<r[i]<((k+1)/10): s=s+1
   hist[k]=s
   s=0
 poisson[l]=hist[9]
plt.hist(poisson,bins=30,range=(-0.5,29.5))	# opraveno J. Matousek

n=1000
N=100
p=0.1
def P_binomial(k,N,p):
  	P = 1.0
  	for i in range(1,N+1):
  		if i <= k and i <= N-k:
  			P = P / i
  		elif i > k and i > N-k:
  			P = P * i
  	return P * pow(p,k) * pow(1-p,N-k)

def P_poisson(k,nu):
	P = m.exp(-nu)
	for i in range(1,k+1):
		P *= nu / i
	return P

binroz=np.empty(N)
nExpPois=np.empty(N)
for k in range(N):
  	binroz[k]   = P_binomial(k,N,p) * n
for k in range(N):
    nExpPois[k] = P_poisson(k,N*p) * n
plt.plot(range(0,N),binroz,color='red',marker='.',label='Binomial')
plt.plot(range(0,N),nExpPois,color='black',marker='.',label='Poisson')
plt.xlim(0,30)
plt.legend()
plt.show()