*
* KUMAC pour montrer l'estimation par maximum de vraisemblance
* Distribution exponentielle avec fonction de resolution normale
*
his/cre/1d 100 'Distribution en distance vraie        ' 60 -5. 10.
his/cre/1d 200 'Distribution en distance observee r=.1' 60 -5. 10.
his/cre/1d 201 'Distribution theorique                ' 60 -5. 10.
his/cre/1d 300 'Distribution en distance observee r=.5' 60 -5. 10.
his/cre/1d 301 'Distribution theorique                ' 60 -5. 10.
his/cre/1d 400 'Distribution en distance observee r=1.' 60 -5. 10.
his/cre/1d 401 'Distribution theorique                ' 60 -5. 10.
*                             Generer echantillon exponentiel
ne = 10000
sigma lam=1.
sigma x=rndm(array([ne]))
sigma x=-lam*log(x)
vec/hfill x  100
*                             Resolution gaussienne
sigma r=.2
sigma r1=rndm(array([ne])) ; sigma r2=rndm(array([ne]))
sigma xp=sin(2.*pi*r1)*sqrt(-2.*log(r2))*r+x
sigma t=array(60,-4.875#9.875)
sigma y=[ne]*0.125*exp(0.5*(r/lam)**2-t/lam)*erfc((r/lam-t/r)/sqrt(2.))/lam
vec/hfill xp 200
his/put/cont 201 y
*
sigma r=.5
sigma xp=sin(2.*pi*r1)*sqrt(-2.*log(r2))*r+x
sigma t=array(60,-4.875#9.875)
sigma y=[ne]*0.125*exp(0.5*(r/lam)**2-t/lam)*erfc((r/lam-t/r)/sqrt(2.))/lam
vec/hfill xp 300
his/put/cont 301 y
*
sigma r=1.
sigma xp=sin(2.*pi*r1)*sqrt(-2.*log(r2))*r+x
sigma t=array(60,-4.875#9.875)
sigma y=[ne]*0.125*exp(0.5*(r/lam)**2-t/lam)*erfc((r/lam-t/r)/sqrt(2.))/lam
vec/hfill xp 400
his/put/cont 401 y
*                             Presentation
opt stat
opt logy
set stat 1111
zone 2 2
his/plot 100
his/plot 200
his/plot 201 'cs'
his/plot 300
his/plot 301 'cs'
his/plot 400
his/plot 401 'cs'