> restart:
E2 - Tests d'hypothèses
Evaluation du 06/12/2004
On charge quelques bibliothèques
Vérifier les numéros de version des bibliothèmes disponibles
>
with(pldx); # il faut la version 24
ever_loaded;
>
with(simul); # il faut la version 16
ever_loaded;
Warning, the name ever_loaded has been redefined
> kernelopts(ASSERT=true):
Création du fichier
M, K:= 1200., 1.9;
ra:= stats[random, normald[ln(M), sqrt(ln(K))]];
n:= 1200; li:= map(exp, [ ra( n )]):
fich:= "F:/docs/Ensait/decis/decis_ds10/dat_ds_10.txt";
hdl:= fopen(fich, WRITE);
for i to n do fprintf(hdl,"%a\n", li[i]); od:
fclose(fich):
Commandes communes
>
macro(moy= stats[describe, mean], var= stats[describe, variance],
ran= stats[describe, range], nmb= stats[describe, count] );
>
macro(histo = xhisto );
xima:=proc(pl); op(1,pl) : convert(%,list): map(op,%): map2(op,2,%): max(op(%)): end:
> macro(talto=stats[transform, tallyinto]);
>
recode:= proc(item) global m; local a,b, c;
a,b,c:= op(op(1,item)), op(2, item);
Weight(a..b, c*(a+b)/2/m)
end :
>
>
Groupe A
A0-Lecture du fichier
Télécharger le fichier "www.douillet.info/~douillet/decis_ds10/dat_ds_10.txt"
Déclarer le nom de la copie locale du fichier
> fich:= "F:/docs/Ensait/decis/decis_ds10/dat_ds_10a.txt";
On lit tous les nombres
>
lire:= table(): tmp:= 'tmp':
hdl:= fopen(fich, READ);
for i to 5000 while tmp <> [ ] do
tmp:= fscanf(hdl,"%a");
lire[i]:= op (tmp);
od:
fclose(fich); n:= i-2 ; i:='i':
On en fait une liste
> li:= [ seq(lire[k], k=1..n)]:
A1-Paramètres de dispersion
> m, s, vals, nn:= moy(li), sqrt(var(li)), ran(li), nmb(li); ASSERT(n=nn);
Histogramme automatique
> kx:= 40: pl1:= histo(li, area=1, numbars=kx):
>
yx:=xima(pl1) :
pl2:= plot([[m,y,y=0..yx], [m-s,y,y=0..yx], [m+s,y,y=0..yx]], color=red, linestyle=[24,8,8]):
> display(pl1, pl2);
![[Maple Plot]](images_maple/maple18.gif)
A2-Modèle lognormal
>
M:='M': K:='K': m=M*sqrt(K), (s/m)^2 = K-1; solve({%});
assign(%);
> pl5:= plot(norlaw(ln(M), sqrt(ln(K)), ln(x))/x, x=vals, color=blue, linestyle=16):
> display(pl1,pl2,pl5);
![[Maple Plot]](images_maple/maple21.gif)
> kx:=40;
>
kz:= 13: deb, fin:= op(vals); pas:= (fin-deb)/kx: mid1:= deb+pas*15: mid2:= deb+pas*19:
les_interv_d:= seq(deb+pas*(k-1)..deb+pas*(k), k=1..kz):
les_interv_f := deb+pas*kz..mid1, mid1..mid2, mid2..fin+0.001:
les_interv := [les_interv_d, les_interv_f]:
> li2:= talto(li, les_interv);
> map2(op,2,li2); [seq(`if`(%[k]>=10, NULL,k), k=1..nops(%))]: ASSERT(%=[ ], %=[ ]);
> display(pl1, histo(li2, color=yellow), view= [DEFAULT, 0.. 0.0001]);
![[Maple Plot]](images_maple/maple31.gif)
> Law1:= unapply(Norlaw(ln(M), sqrt(ln(K)), ln(x)), x);
>
proc(z) local a,b: a,b:= op(z); Weight(a..b, n*(Law1(b)-Law1(a))) end:
th2:= map(%, les_interv);
> display( histo(li2, color=yellow), histo(th2));
![[Maple Plot]](images_maple/maple41.gif)
> nu:= nops(li2)-3; add((op(2, th2[k])-op(2, li2[k]))^2/op(2, th2[k]), k=1..nops(li2));
A3-Loi "en masse"
> li3:= map(recode, li2);
> nmb(li3);
> m_, s_:= moy(li3), sqrt(var(li3));
>
pl23:=histo(li3, area=1):
yx:=xima(%) :
pl24:= plot([[m_,y,y=0..yx], [m_-s_,y,y=0..yx], [m_+s_,y,y=0..yx]], color=red, linestyle=[24,8,8]):
Modèle issu du modèle "en nombre"
>
Law3:= unapply(Norlaw(ln(M*K), sqrt(ln(K)), ln(x)), x):
law3:= unapply(norlaw(ln(M*K), sqrt(ln(K)), ln(x))/x, x):
pl25:= plot(law3(x), x=vals, color=blue, linestyle=16):
> display( pl23, pl24,pl25);
![[Maple Plot]](images_maple/maple54.gif)
> m_; (s/m)^2+1; m*%;
>
mk_th3:= proc(z) local a,b: a,b:= op(z); Weight(a..b, n*(Law3(b)-Law3(a))) end:
th30:= map(mk_th3, les_interv): n - nmb(th30);
subsop(2= op(2, th30[-1]) +%, th30[-1]): th3:= subsop(-1=%, th30):
ASSERT(nmb(th3)=n);
Test du Chi2 sans correction
>
nu:= nops(li3)-1; add((op(2, th30[k])-op(2, li3[k]))^2/op(2, th30[k]), k=1..nops(li3));
chi_red=(%-nu)/sqrt(2.*nu);
>
Test du Chi2
>
nu:= nops(li3)-1; add((op(2, th3[k])-op(2, li3[k]))^2/op(2, th3[k]), k=1..nops(li3));
chi_red=(%-nu)/sqrt(2.*nu);
Test du Chi2 (sans le dernier intervalle)
>
nu:= nops(li3)-2; add((op(2, th3[k])-op(2, li3[k]))^2/op(2, th3[k]), k=1..nops(li3)-1);
chi_red=(%-nu)/sqrt(2.*nu);
>
>
>
Tentative de détermination d'un modèle à partir de m_ et s_
> M_:='M_': K_:='K_': m_=M_*sqrt(K_), (s_/m_)^2 = K_-1; solve({%}); assign(%); M_/K_;
> pl27:= plot(norlaw(ln(M_), sqrt(ln(K_)), ln(x))/x, x=vals, color=red, linestyle=16):
> display( pl23, histo(th3, color=yellow), pl25, pl27, view=[0..fin, DEFAULT]);
![[Maple Plot]](images_maple/maple71.gif)
> m,s; m_, s_; m*K-m_; %/m_/sqrt(K-1); %*sqrt(n*1.);
> evalf(%/(s_/sqrt(n)));
>
>
>
>
Groupe B
B0-Lecture du fichier
Télécharger le fichier "www.douillet.info/~douillet/decis_ds10/dat_ds_10.txt"
Déclarer le nom de la copie locale du fichier
> fich:= "F:/docs/Ensait/decis/decis_ds10/dat_ds_10b.txt";
On lit tous les nombres
>
lire:= table(): tmp:= 'tmp':
hdl:= fopen(fich, READ);
for i to 5000 while tmp <> [ ] do
tmp:= fscanf(hdl,"%a");
lire[i]:= op (tmp);
od:
fclose(fich); n:= i-2 ; i:='i':
On en fait une liste
> li:= [ seq(lire[k], k=1..n)]:
B1-Paramètres de dispersion
> m, s, vals, nn:= moy(li), sqrt(var(li)), ran(li), nmb(li); ASSERT(n=nn);
Histogramme automatique
> kx:= 40: pl1:= histo(li, area=1, numbars=kx):
>
yx:=xima(pl1) :
pl2:= plot([[m,y,y=0..yx], [m-s,y,y=0..yx], [m+s,y,y=0..yx]], color=red, linestyle=[24,8,8]):
> display(pl1, pl2);
![[Maple Plot]](images_maple/maple82.gif)
B2-Modèle lognormal
>
M:='M': K:='K': m=M*sqrt(K), (s/m)^2 = K-1; solve({%});
assign(%);
> pl5:= plot(norlaw(ln(M), sqrt(ln(K)), ln(x))/x, x=vals, color=blue, linestyle=16):
> display(pl1,pl2,pl5);
![[Maple Plot]](images_maple/maple85.gif)
> kx:=45;
>
kz:= 12: deb, fin:= op(vals); pas:= (fin-deb)/kx: mid1:= deb+pas*(kz+2) : mid2:= deb+pas*(kz+5): mid3:= deb+pas*(kz+8):
les_interv_d:= seq(deb+pas*(k-1)..deb+pas*(k), k=1..kz):
les_interv_f := deb+pas*kz..mid1, mid1..mid2, mid2..mid3, mid3..fin+0.001:
les_interv := [les_interv_d, les_interv_f]:
> li2:= talto(li, les_interv);
> map2(op,2,li2); [seq(`if`(%[k]>=10, NULL,k), k=1..nops(%))]: ASSERT(%=[ ], %=[ ]);
> display(pl1, histo(li2, color=yellow), view= [DEFAULT, 0.. 0.0001]);
![[Maple Plot]](images_maple/maple95.gif)
> Law1:= unapply(Norlaw(ln(M), sqrt(ln(K)), ln(x)), x);
>
proc(z) local a,b: a,b:= op(z); Weight(a..b, n*(Law1(b)-Law1(a))) end:
th2:= map(%, les_interv);
> display( histo(li2, color=yellow), histo(th2));
![[Maple Plot]](images_maple/maple105.gif)
> nu:= nops(li2)-3; add((op(2, th2[k])-op(2, li2[k]))^2/op(2, th2[k]), k=1..nops(li2));
B3-Loi "en masse"
> li3:= map(recode, li2);
> nmb(li3);
> m_, s_:= moy(li3), sqrt(var(li3));
>
pl23:=histo(li3, area=1):
yx:=xima(%) :
pl24:= plot([[m_,y,y=0..yx], [m_-s_,y,y=0..yx], [m_+s_,y,y=0..yx]], color=red, linestyle=[24,8,8]):
Modèle issu du modèle "en nombre"
>
Law3:= unapply(Norlaw(ln(M*K), sqrt(ln(K)), ln(x)), x):
law3:= unapply(norlaw(ln(M*K), sqrt(ln(K)), ln(x))/x, x):
pl25:= plot(law3(x), x=vals, color=blue, linestyle=16):
> display( pl23, pl24,pl25);
![[Maple Plot]](images_maple/maple118.gif)
> m_; (s/m)^2+1; m*%;
>
mk_th3:= proc(z) local a,b: a,b:= op(z); Weight(a..b, n*(Law3(b)-Law3(a))) end:
th30:= map(mk_th3, les_interv): n - nmb(th30);
subsop(2= op(2, th30[-1]) +%, th30[-1]): th3:= subsop(-1=%, th30):
ASSERT(nmb(th3)=n);
Test du Chi2
>
nu:= nops(li3)-1; add((op(2, th3[k])-op(2, li3[k]))^2/op(2, th3[k]), k=1..nops(li3));
chi_red=(%-nu)/sqrt(2.*nu);
Test du Chi2 sans correction
>
nu:= nops(li3)-1; add((op(2, th30[k])-op(2, li3[k]))^2/op(2, th30[k]), k=1..nops(li3));
chi_red=(%-nu)/sqrt(2.*nu);
>
>
Test du Chi2 (sans le dernier intervalle)
>
nu:= nops(li3)-2; add((op(2, th3[k])-op(2, li3[k]))^2/op(2, th3[k]), k=1..nops(li3)-1);
chi_red=(%-nu)/sqrt(2.*nu);
>
>
>
Tentative de détermination d'un modèle à partir de m_ et s_
> M_:='M_': K_:='K_': m_=M_*sqrt(K_), (s_/m_)^2 = K_-1; solve({%}); assign(%); M_/K_;
> pl27:= plot(norlaw(ln(M_), sqrt(ln(K_)), ln(x))/x, x=vals, color=red, linestyle=16):
> display( pl23, histo(th3, color=yellow), pl25, pl27, view=[0..fin, DEFAULT]);
![[Maple Plot]](images_maple/maple135.gif)
>
>
File d'attente
Procédure simul_wait (un seul serveur)
> _seed:= 11112003; lambda:= 5.5; mu:=8; N:=5000;
>
with(stats): with(simul): macro(histo = xhisto, talto=stats[transform, tallyinto] ):
macro(moy= stats[describe, mean], var= stats[describe, variance], ran= stats[describe, range], nmb= stats[describe, count] );
Warning, these names have been redefined: anova, describe, fit, importdata, random, statevalf, statplots, transform
>
>
gx:= proc() global _seed; _seed := irem(427419669081*_seed, 999999999989);
-1/alpha*ln(evalf(_seed/999999999989)) end proc:
ga:= subs(alpha=lambda, eval(gx)): gb:= subs(alpha=mu, eval(gx)):
>
simul_wait:= proc(N) global mem, ga, gb:
local load,k,ech,i,laps,event,n,gg,compar;
load:=0: k:=0: ech:=[[0,0,0]]: mem:= table():
compar:= (A,B)-> evalb(A[1]<B[1]);
for i to 2*N do
laps, event, n := op(ech[1]);
ech:= map(A->[A[1]-laps,A[2],A[3]], subsop(1=NULL, ech));
load:= `if`(load>0,load-laps,0);
if event=0 # arrivée
then n:= n+1; k:=k+1; gg:= gb();
load:= load+gg; mem[n]:= [load,gg,k];
ech:= sort([ [ga(),0,n],[load,1,n], op(ech) ], compar);
else # départ
k:=k-1;
fi;
od: end:
>
compar_histo:= proc(ex, th) global mihis; local nn, mi; nn:= nops(th);
if nops(ex)<>nn then error `: les deux structures n'ont pas la même longueur` fi;
mihis:= [seq(Weight( op(1,th[j]), min(op(2,th[j]),op(2,ex[j]))), j=1..nn)];
display(histo(mihis, color=grey, area=count), histo(ex, color=cyan, area=count)
, histo(th, color=yellow, area=count));
end:
> stamp:= time(): simul_wait(N): print(time()-stamp);
Exploitation
> Matrix([seq(mem[j], j=1..10)]);
![_rtable[33054252]](images_maple/maple143.gif)
Séjour
>
tmp:= sort([seq(mem[j][1],j=1..N)]):
law:= unapply((mu-lambda)*exp(-(mu-lambda)*t),t);
display(histo(tmp, area=1), plot(law(t), t=0..4));
![[Maple Plot]](images_maple/maple145.gif)
> vals:= ran(tmp);
> Law:= unapply(1- exp(t*(-mu+lambda)),t);
>
nu:=9; seq(solve(k/(nu+1)=Law(t)), k=0..nu): boun:= %, 1.3*max(%[-1]*2, rhs(vals));
les_interv:= [seq(boun[j]..boun[j+1], j=1..nu+1)]: nops(%);
> exper:= talto(tmp, les_interv);
>
proc(z) local a,b: a,b:= op(z); Weight(a..b, N*(Law(b)-Law(a))) end:
theor:= map(%, les_interv);
> display(compar_histo(exper, theor), plot(N*law(t), t=0..4), view=[0..2, DEFAULT]);
![[Maple Plot]](images_maple/maple162.gif)
> nops(theor); nu;
>
add( (op(2,theor[j])-op(2,exper[j]))^2/op(2,theor[j]), j=1..nops(theor)-1);
(%-nu)/sqrt(2.*nu); (%%-nu-1)/sqrt(2.*(nu-1));
>
>
Hors sujet
Service
>
tmp:= [seq(mem[j][2],j=1..N)]: histo(tmp, area=1, numbars=20):
display(%, plot(mu*exp(-mu*t), t=0..4), view=[0..2, DEFAULT]);
![[Maple Plot]](images_maple/maple168.gif)
Nb clients
> tmp:= [seq(mem[j][3],j=1..N)]: stats[transform, tally](%): jx:= (max@op@map2)(op,1,%):
>
exper:= stats[transform, tallyinto](tmp, [seq(j-1/2..j+1/2, j=1..jx)]):
pl_nb:= xhisto(%, area=1): evalf([moy(tmp), var(tmp)]);
![[3.435000000, 7.339775000]](images_maple/maple169.gif){kind=small}
> rho:= lambda/mu*1.; theor:= [seq(Weight(k+1/2..k+3/2, N*(1-rho)*rho^k), k=0..jx-1)]:
> compar_histo(exper, theor);
![[Maple Plot]](images_maple/maple171.gif)
>
>