* Nom du fichier: tvie.kumac
*
* KUMAC pour montrer l'estimation par maximum de vraisemblance
*                                      MP 3/95
his/cre/1d 100 'Distribution en distance x (n=1000)' 50  0. 5.
his/cre/1d 101 'Distribution avec l max L          ' 50  0. 5.
his/cre/1d 200 'Vraisemblance vs Lambda    (n=1000)' 50  .9 1.1
his/cre/1d 110 'Distribution en distance x (n=10000)' 50  0. 5.
his/cre/1d 111 'Distribution avec l max L           ' 50  0. 5.
his/cre/1d 210 'Vraisemblance vs Lambda    (n=10000)' 50  .9 1.1
his/cre/1d 120 'Distribution en distance x (n=50000)' 50  0. 5.
his/cre/1d 121 'Distribution avec l max L           ' 50  0. 5.
his/cre/1d 220 'Vraisemblance vs Lambda    (n=50000)' 50  .9 1.1
*                             Generer echantillon exponentiel
ne = 50000
sigma lam=1.
sigma x3=rndm(array([ne]))
sigma x3=-lam*log(x3)
vec/copy x3(1:1000) x1
vec/copy x3(1:10000) x2
vec/hfill x1 100
vec/hfill x2 110
vec/hfill x3 120
*                             Estimations max L  (detecteur infini)
sigma xb1=vsum(x1)/1000.
sigma eb1=xb1/sqrt(1000.)
sigma sb1=$hinfo(100,'rms')/sqrt(1000.)
sigma ex1=1000.*.1*exp(-array(50,0.#5.)/xb1)/xb1
his/put/cont 101 ex1
sigma xb2=vsum(x2)/10000.
sigma eb2=xb2/sqrt(10000.)
sigma sb2=$hinfo(110,'rms')/sqrt(10000.)
sigma ex2=10000.*.1*exp(-array(50,0.#5.)/xb2)/xb2
his/put/cont 111 ex2
sigma xb3=vsum(x3)/50000.
sigma eb3=xb3/sqrt(50000.)
sigma sb3=$hinfo(120,'rms')/sqrt(50000.)
sigma ex3=50000.*.1*exp(-array(50,0.#5.)/xb3)/xb3
his/put/cont 121 ex3
*                             L(lambda)
sigma l=array(50,.9#1.1)
sigma lnl=-1000.*log(l)-1000.*xb1/l
sigma lma=-1000.*log(xb1)-1000.
sigma lnl=exp(lnl-lma)
his/put/cont 200 lnl
sigma lnl=-10000.*log(l)-10000.*xb2/l
sigma lma=-10000.*log(xb2)-10000.
sigma lnl=exp(lnl-lma)
his/put/cont 210 lnl
sigma lnl=-50000.*log(l)-50000.*xb3/l
sigma lma=-50000.*log(xb3)-50000.
sigma lnl=exp(lnl-lma)
his/put/cont 220 lnl
*                             Presentation
opt stat
set stat 1111
zone 2 2
his/plot 100
his/plot 101 'sc'
itx 1.50 80.00 'Lambda       = '//$format(xb1,f6.3)
itx 1.50 75.00 'Erreur dist   = '//$format(eb1,f6.3)
itx 1.50 70.00 'Erreur echa = '//$format(sb1,f6.3)
*
his/plot 110
his/plot 111 'sc'
itx 1.50 800.00 'Lambda       = '//$format(xb2,f6.3)
itx 1.50 750.00 'Erreur dist   = '//$format(eb2,f6.3)
itx 1.50 700.00 'Erreur echa = '//$format(sb2,f6.3)
*
his/plot 120
his/plot 121 'sc'
itx 1.50 4000.00 'Lambda       = '//$format(xb3,f6.3)
itx 1.50 3750.00 'Erreur dist   = '//$format(eb3,f6.3)
itx 1.50 3500.00 'Erreur echa = '//$format(sb3,f6.3)
*
opt nsta
his/plot 220 'c'
his/plot 210 'cs'
his/plot 200 'cs'