Groupes an1/an2

Lecture du fichier

Télécharger le fichier "www.douillet.info/~douillet/maotp/maotp_ds10/dat_maotp_ds10.txt"

Déclarer le nom de la copie locale du fichier

> fich:= "F:/docs/Ensait/maotp/maotp_ds10/dat_maotp_ds10a.txt";

fich :=

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':

hdl := 0

n := 1000

On en fait une liste

> li2:= [seq(lire[k], k=1..n)]:

Paramètres de dispersion

> macro(moy= stats[describe, mean], var= stats[describe, variance],
ran= stats[describe, range], nmb= stats[describe, count] );

moy, var, ran, nmb

> m:= moy(li2);
s, s2:= sqrt(var(li2)), var(li2);
vals:= ran(li2);
ASSERT(n=nmb(li2));

m := 1711.508388

s, s2 := 1484.902860, 2204936.504

vals := 138.3919284 .. 10563.07935

Visualisation

> macro(histo = xhisto, talto=stats[transform, tallyinto] );

moy, var, ran, nmb, histo, talto

Détermination du nombre de barres

> histo(li2, area=1); pl1a:= %:

[Maple Plot]

> histo(li2, area=1, numbars=100); pl1b:= %:

[Maple Plot]

> nbx:= 40: histo(li2, area=1, numbars=nbx); pl1:= %:

[Maple Plot]

> xima:= proc(pl); op(1,pl): convert(%,list): map(op,%): map2(op,2,%); (max@op)(%); end:

> xima(pl1a), xima(pl1b), xima(pl1); yx:= %[3];

.4408765285e-3, .6331125082e-3, .5486975080e-3

yx := .5486975080e-3

> 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]

Modèle lognormal

Paramètres loi lognormale

> M:='M': K:='K': m=M*sqrt(K), (s/m)^2 = K-1; solve({%}); assign(%);

> M, K;

1711.508388 = M*sqrt(K), .7527279176 = K-1

{K = 1.752727918, M = 1292.771532}

1292.771532, 1.752727918

> pl5:= plot(norlaw(ln(M), sqrt(ln(K)), ln(x))/x, x=0..rhs(vals), color=blue, linestyle=16):

> top:= exp(ln(M)+2.5*sqrt(ln(K)));

top := 8411.300144

> cut:= proc(li, num);
stats[transform, split[num]](stats[transform, statsort](li)):
map(z-> Weight(ran(z), nmb(z)), %);
end;

cut := proc (li, num) stats[transform,split[num]](s...
cut := proc (li, num) stats[transform,split[num]](s...
cut := proc (li, num) stats[transform,split[num]](s...
cut := proc (li, num) stats[transform,split[num]](s...

> tal2:= cut(li2,20);

tal2 := [Weight(138.3919284 .. 369.4617675,50), Wei...
tal2 := [Weight(138.3919284 .. 369.4617675,50), Wei...
tal2 := [Weight(138.3919284 .. 369.4617675,50), Wei...
tal2 := [Weight(138.3919284 .. 369.4617675,50), Wei...
tal2 := [Weight(138.3919284 .. 369.4617675,50), Wei...
tal2 := [Weight(138.3919284 .. 369.4617675,50), Wei...
tal2 := [Weight(138.3919284 .. 369.4617675,50), Wei...
tal2 := [Weight(138.3919284 .. 369.4617675,50), Wei...
tal2 := [Weight(138.3919284 .. 369.4617675,50), Wei...
tal2 := [Weight(138.3919284 .. 369.4617675,50), Wei...

> pl3:= histo(tal2, area=1, color=yellow):

> display(pl1,pl2,pl5, view=[0..top, DEFAULT]); display(pl3,pl2,pl5, view=[0..top, DEFAULT]);

[Maple Plot]

[Maple Plot]

>

Loi "en masse"

> newcode:= unapply(Weight(z,z/m), z);

newcode := proc (z) options operator, arrow; Weight...

> li3:= map(newcode, sort(li2)):
nmb(li3);
m_:= moy(li3); v_:= var(li3): s_:= sqrt(%);

999.9999927

m_ := 2999.808537

s_ := 2229.433270

> pl21:=histo(li3, area=1, numbars=nbx): yx:= xima(%):

> pl22:= plot([[m_,y,y=0..yx], [m_-s_,y,y=0..yx], [m_+s_,y,y=0..yx]], color=red, linestyle=[24,8,8]):

> tal3:= cut(li3,20);

tal3 := [Weight(138.3919284 .. 656.5789293,49.99999...
tal3 := [Weight(138.3919284 .. 656.5789293,49.99999...
tal3 := [Weight(138.3919284 .. 656.5789293,49.99999...
tal3 := [Weight(138.3919284 .. 656.5789293,49.99999...
tal3 := [Weight(138.3919284 .. 656.5789293,49.99999...
tal3 := [Weight(138.3919284 .. 656.5789293,49.99999...
tal3 := [Weight(138.3919284 .. 656.5789293,49.99999...
tal3 := [Weight(138.3919284 .. 656.5789293,49.99999...
tal3 := [Weight(138.3919284 .. 656.5789293,49.99999...
tal3 := [Weight(138.3919284 .. 656.5789293,49.99999...

> pl23:= histo(tal3, area=1, color=yellow):

> jx:=nbx/2: [0, seq(exp(INorlaw(ln(K*M), sqrt(ln(K)), j/jx)), j=1..jx-1), rhs(vals)*1.1]:
les_interv:= seq(%[j]..%[j+1], j=1..nops(%)-1);
tal3th:= talto(li3, [les_interv]):

les_interv := 0 .. 660.8536550, 660.8536550 .. 867....
les_interv := 0 .. 660.8536550, 660.8536550 .. 867....
les_interv := 0 .. 660.8536550, 660.8536550 .. 867....
les_interv := 0 .. 660.8536550, 660.8536550 .. 867....
les_interv := 0 .. 660.8536550, 660.8536550 .. 867....

> pl24:= histo(tal3th, area=1, color=cyan):

> pl25:= plot(norlaw(ln(M*K), sqrt(ln(K)), ln(x))/x, x=vals, color=blue, linestyle=16):

> display( pl21, pl22,pl25); display( pl23, pl22,pl25);display( pl24, pl22,pl25);

[Maple Plot]

[Maple Plot]

[Maple Plot]

> display(pl25,pl21,pl23,pl24);

[Maple Plot]

>

> m_; (s/m)^2+1; m*%;

2999.808537

1.752727918

2999.808534

> 1+(s/m)^2; 1+(s_/m_)^2;

1.752727918

1.552334133

>