[ Retour aux Travaux Dirigés de SAS ]
/* Statistiques avec SAS, Corrigé du TD 6 */
/* Jean-Sebastien Roy (js@jeannot.org), 2000 */
/* Estimation non-paramétrique d'une densité */
/* observations tirées de deux lois normales */
data deuxn;
do i=1 to 200;
xi=normal(0); output;
xi=normal(0)*1.5+4; output;
end;
keep xi;
run;
/* Estimation de la densité avec Insight */
proc insight deuxn ;
dist xi;
run;
/* La grille des points sur lesquels
on va chercher à estimer la densité */
data grid;
xmin=-4;
xmax=9;
nb=100;
step=(xmax-xmin)/(nb-1);
do x=xmin to xmax by step;
output;
end;
keep x;
run;
/* Calcul de la valeur exacte de la densité */
data exact;
set grid;
PI=3.14159265359;
y1=1/(sqrt(2*PI)*1)*exp(-0.5*((x-0)/1)**2);
y2=1/(sqrt(2*PI)*1.5)*exp(-0.5*((x-4)/1.5)**2);
y=(y1+y2)/2;
name="EXACT";
keep name x y;
run;
/* Quelques statistiques utiles sur
la distribution de notre loi */
proc univariate deuxn;
var xi;
output stats
n=n /* Nombre d'observations */
std=std /* Ecart type */
qrange=qrange; /* Intervale inter-quartiles */
run;
/* Fenêtres par rapport auquelles
on va réaliser l'estimation */
data h;
set stats;
/* Une fenêtre pour donner
l'ordre de grandeur */
snr=std*(4/(3*n))**0.2;
fmin=log(0.05);
fmax=log(3);
nb=10;
step=(fmax-fmin)/(nb-1);
do f=fmin to fmax by step;
h=exp(f)*snr;
output;
end;
keep h;
run;
/* Pour l'estimation, si la proc SQL de
SAS était très efficace on écrirait :
proc sql;
create table kesti as
select h, x, 1/(sqrt(2*3.14159265359)*h)*mean(exp(-0.5*((x-xi)/h)**2)) as y
from grid, deuxn, h
group by h,x;
quit;
Malheureusement elle est assez lente,
alors on procède en deux temps.
*/
/* 1: Creation d'une vue contenant toutes
les valeurs de K((x-xi)/h)
On utilise une vue pour éviter que SAS
crée physiquement la table sur disque,
cette table pouvant etre très grosse
(nb obs*nb pts de grille*nb de bandes passsantes observations)
*/
proc sql;
create view cross as
select h, x, 1/(sqrt(2*3.14159265359)*h)*exp(-0.5*((x-xi)/h)**2) as y
from grid, deuxn, h;
quit;
/* 2: Calcul des moyennes de ces valeurs */
proc summary cross ;
class h x;
var y;
output kesti(_type_ _freq_) ;
run;
/* Tracé des estimations */
proc gplot kesti;
symbol join; /* Des lignes */
plot y*x=h; /* Une courbe par valeur de h */
run;
quit;
/* Calcul de l'erreur par rapport a la
densité exacte */
proc sql;
create table erreurs as
select h,
sqrt(mean((kesti.y-exact.y)**2)) as rmse
from kesti, exact
where kesti.x=exact.x
group by h;
quit;
proc gplot erreurs;
plot rmse*h;
run;
quit;
/* Quelques heuristiques de choix de la
fenêtre
On associe a chaque valeur un nom
pour pouvoir la reperer dans le graphique
*/
data h;
set stats;
length name $ 10;
h=std*(4/(3*n))**0.2; name="SNR="||putn(h,'BEST6.'); output;
h=0.9*min(std,qrange/1.34)*n**-0.2; name="SROT="||putn(h,'BEST6.'); output;
h=3*std*(1/(70*sqrt(3.14159265359*n)))**0.2; name="OS="||putn(h,'BEST6.'); output;
keep h name;
run;
/* On refait l'estimation.
Pas besoin de recreer la vue */
proc summary cross ;
class h x;
var y;
output kesti(_type_ _freq_) ;
run;
/* On récupère le nom associé à h
pour le graphique */
proc sort h; by h; run;
data kesti;
merge kesti h;
by h;
run;
proc gplot kesti;
plot y*x=name;
run;
quit;
/* Sous SAS v8, estimation de la densité
par la meme méthode que précédement
(infiniment plus rapide)
et optimisation de la fenêtre.
*/
proc kde deuxn kde((xi=x density=y) xi density);
var xi;
run;
data kde;
set kde;
name="KDE";
run;
/* On trace l'ensemble */
data ;
set kesti exact kde;
run;
proc gplot ;
plot y*x=name;
run;
quit;
/* Estimation fonctionnelle
non paramétrique */
/* Des observation d'une fonction
non linéaire : sur [0,1]
1/2-2*x si x<1/5
0.1+x-1/5 sinon
auquel on ajoute une perturbation
normale (sigma = 1/10)
*/
data points;
do i=1 to 100;
xi=uniform(0);
if xi<0.2 then yi=0.5-xi*2;
else yi=0.1+(xi-0.2);
yi=yi+normal(0)/10;
output;
end;
run;
/* Memes statistiques que précédement */
proc univariate points;
var xi;
output stats n=n std=std qrange=qrange;
run;
/* Memes fenetres que précédement */
data h;
set stats;
srn=std*(4/(3*n))**0.2;
fmin=log(0.05);
fmax=log(3);
nb=10;
step=(fmax-fmin)/(nb-1);
do f=fmin to fmax by step;
h=exp(f)*srn;
output;
end;
keep h;
run;
/* La aussi, une grille de points
ou evaluer la fonction */
data grid;
xmin=0;
xmax=1;
nb=100;
step=(xmax-xmin)/(nb-1);
do x=xmin to xmax by step;
output;
end;
run;
/* Meme probleme que précedement avec SAS.
On procède en deux étapes */
/* 1: Creation des numérateurs et dénominateurs
y=yi et w=K((xi-x)/h) */
proc sql;
create view cross as
select h, x, yi as y,
exp(-0.5*((x-xi)/h)**2) as w
from grid, points, h;
quit;
/* 2: Moyenne pondérée par w des y */
proc summary cross ;
class h x;
var y;
weight w;
output festi(_type_ _freq_) ;
run;
/* On rajoute les points observées
pour le graphique */
data ;
set festi points((xi=x yi=y));
output;
/* Si c'est un point observé (h=.) */
if h=. then do;
/* On fait succéder à ce point un
point manquant pour interrompre
le tracé */
y=.;
output;
end;
run;
proc sort ;
by h x y;
run;
proc gplot ;
symbol join dot 0.1 2;
plot y*x=h / ;
run;
quit;
/* Estimation fonctionnelle à l'aide
de SAS/Insight */
proc insight points ;
fit yi=xi;
run;
/* Dans SAS v8, deux procédures permettent
de faire des esimations fonctionnelles
non paramétriques :
LOESS et TPSPLINE
*/
/* Calcul de la validation croisée
Meme principe :
Calcul du toutes les valeurs de
w=K((xi-xj)/h) avec i!=j
et de cv=yi-yj
(Estimation en xj en enlevant le point j)
*/
proc sql ;
create view cross as
select h, t1.i as i,
exp(-0.5*((t1.xi-t2.xi)/h)**2) as w,
t2.yi-t1.yi as cv
from points as t1, points as t2, h
where t1.i ne t2.i;
quit;
/* Moyenne pondérée des cv pour tous les i */
proc summary cross ;
class h i;
var cv;
weight w;
output temp(_type_ _freq_ i) ;
run;
/* Et finalement, moyenne des carrés des cv */
data temp;
set temp;
if cv ne . then cv=cv**2;
run;
proc summary temp ;
class h;
var cv;
output cv(_type_ _freq_) ;
run;
/* Visualisation de la validation croisée
en fonction de h */
proc gplot cv;
symbol join 1 none;
plot cv*h;
run;
quit;