*
* .kumac pour comparer les distributions binomiales,
*  poissonienne et gaussienne pour une petit probabilite' p
*  et un nombre croissante des essaies N
*

macro b N Ngra
*						Distribution Binomiale
*						avec parametres N,p
*
 v/del *
 p=0.05
 si den  = [p]*array([N]+1,[N]+1#1)
 si nom  = array([N]+1,0#[N])
 vec/copy nom k
*						0! = 1 !!
 v/in den(1) 1
 v/in nom(1) 1
*						Divise chaque terme
*						pour eviter overflow
 vdiv den nom rat 
 si rrat = prod(rat)
 si bino = rrat*(1-[p])**([N]-k)
*					 	Nombre reel sur 16 byte !
 si bino = max(bino,1.7015E-38)
*v/pr bino
 1d 1000 'Distribution binomiale' [N]+1 -0.5 [N]+0.5
 put/cont 1000 bino
 h/pl 1000(:[Ngra])
return
*
macro p N Ngra
*						Distribution Poissonienne
*						avec parametres mu=Np
*
 v/del *
 p=0.05
 mu=[p]*[N]
 si den  = [mu]*array([N]+1)
 si nom  = array([N]+1,0#[N])
 vec/copy nom k
 v/in den(1) 1
 v/in nom(1) 1
 vdiv den nom rat 
 si rrat = prod(rat)
 si pois = exp(-[mu])*rrat
 si pois = max(pois,1.7015E-38)
*v/pr pois
 1d 2000 'Distribution poissonienne' [N]+1 -0.5 [N]+0.5
 put/cont 2000 pois
 h/pl 2000(:[Ngra])
 
return
*
macro g N Ngra
*						Distribution Gaussienne
*						avec parametres mu=Np
*						et s=sqrt(mu(1-p))
*
 v/del *
 p=0.05
 mu=[p]*[N]
 s=sqrt([mu]*(1-[p]))
 si k = array([N]+1,0#[N])
 si gaus = exp(-(k-[mu])**2/2/[s]/[s])/2/sqrt(acos(0.))/[s]
*v/pr gaus
 1d 3000 'Distribution gaussienne' [N]+1 -0.5 [N]+0.5
 put/cont 3000 gaus
 h/pl 3000(:[Ngra])
return
*
macro do N Ngra
 zone 2 2
 h/del 0
 exec b [N] [Ngra] 
 exec p [N] [Ngra] 
 exec g [N] [Ngra] 
 null 0 1 0 1 AB
 itx .2 .8 'p=0.05'
return
*
macro go
 set stat 1100
 title 'Comparaison entre distributions'
 exec do 10 5
 itx .2 .6 'n=10'
 wait
 exec do 20 8
 itx .2 .6 'n=18'
 wait
 exec do 100 20
 itx .2 .6 'n=100'
 wait
 exec do 400 50
 itx .2 .6 'n=400'
 wait
 exec do 1500 150
 itx .2 .6 'n=1500'
 wait
return