{"word_count": 718, "character_count": 7315}
Lors du dernier cours de statistique, nous avons présenter les tests d'ajustment de lois. Nous avions illustré ces tests à partir de la taille d'individus (déjà utilisé pour présenter l'ajustement de lois, et l'estimation de densité), correspondant à .&#62; Davis=read.table(+ "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt")&#62; Davis[12,c(2,3)]=Davis[12,c(3,2)]&#62; X=Davis$height&#62; n=length(X)On notera  la statistique d'ordre, au sens oùParmi les outils graphiques, nous avons vu le PP plot (graphique probabilité-probabilité) et le QQ plot (graphique quantile). Le code pour créer un PP plot peut être le suivant&#62; PP=function(Y,F=pnorm){+   n=length(Y)+   x=F(sort(Y))+   y=seq(1/n/2,1-1/n/2,by=1/n)+   return(data.frame(x=x,y=y))+ }qui représente (à un détail près) le nuage de pointset celui pour le QQ plot serait&#62; QQ=function(Y,Q=qnorm){+   n=length(Y)+   x=Q(seq(1/n/2,1-1/n/2,by=1/n))+   y=sort(Y)+   return(data.frame(x=x,y=y))+ }qui représente (toujours à un détail près) le nuage de pointsoù  est la loi que l'on cherche à tester, au sens où on a  avec comme hypothèse alternative .Rappelons que la fonction de réparition empirique estNous avions vu que .En cours, nous avions utilisé une séquence  alors qu'ici je centre ces valeurs, en prenant . La première raison est que le décallage à droite n'est pas joli, et surtout  n'a peut pas de sens (si la loi n'est pas bornée).On va ici chercher à tester une loi normale , notée &#62; m0=170&#62; s0=9&#62; F0=function(x) pnorm(x,m0,s0)&#62; Q0=function(x) qnorm(x,m0,s0)Le graphique probabilité-probabilité est ici&#62; plot(PP(X,F=F0),xlab="Theoretical Probabilities",+      ylab="Empirical Probabilities")&#62; abline(a=0,b=1,col=rgb(1,0,0,.4))alors que le graphique quantile-quantile est&#62; plot(QQ(X,Q=Q0),xlab="Theoretical Quantiles",+      ylab="Empirical Quantiles")&#62; abline(a=0,b=1,col=rgb(1,0,0,.4))Si  est la vraie loi, les PP plot et QQ plot sont des nuages de points qui doivent être alignés suivant la première bissectrice. Comme on fait des statistiques, "alignés" doit se comprendre comme "presque alignés", ou "pas trop éloignés de la première bissectrice". Pour savoir si on n'est pas trop loin, on peut utiliser des simulations. Il suffit de simuler des échantillons de même taille que \boldsymbol{x} suivant la loi , et de tracer les PP et QQ plots&#62; MVPP=MVQQ=matrix(NA,500,length(X))&#62; for(s in 1:500) MVPP[s,]=PP(Q0(runif(length(X))),F0)$x&#62; for(s in 1:500) MVQQ[s,]=QQ(Q0(runif(length(X))),Q0)$yCommençons par le graphique quantile-quantile&#62; plot(QQ(X,Q=Q0),xlab="Theoretical Quantiles",+      ylab="Empirical Quantiles",col="white")&#62; abscisse=Q0(seq(1/n/2,1-1/n/2,by=1/n))&#62; for(s in 1:500) lines(abscisse,MVQQ[s,],col="light blue",lwd=.5)&#62; points(QQ(X,Q=Q0))&#62; abline(a=0,b=1,col=rgb(1,0,0,.4))On peut calculer les quantiles empiriques, vus verticalement&#62; Q1=apply(MVQQ,2,function(x) quantile(x,.025))&#62; Q2=apply(MVQQ,2,function(x) quantile(x,.975))&#62; lines(abscisse,Q1,col="red")&#62; lines(abscisse,Q2,col="red")Ce sont les courbes rouges ci-dessus. On peut aussi représenter les points qui sortes de la région de confiance, ponctuellement&#62; plot(QQ(X,Q=Q0),xlab="Theoretical Quantiles",+      ylab="Empirical Quantiles",col="white")&#62; polygon(c(abscisse,rev(abscisse)),c(Q1,rev(Q2)),border=NA,col="yellow")&#62; id=((QQ(X,Q=Q0)$y&#62;=Q1)&#38;(QQ(X,Q=Q0)$y&#60;=Q2))+1&#62; points(QQ(X,Q=Q0),col=c('red','blue')[id])Notons qu'ici 5 points, sur 200, sortent des régions de confiance.Dans le cas du graphique probabilité-probabilité, on utilise&#62; plot(PP(X,F=F0),xlab="Theoretical Probabilities",+      ylab="Empirical Probabilities",col="white")&#62; ordonnee=seq(1/n/2,1-1/n/2,by=1/n)&#62; for(s in 1:500) lines(MVPP[s,],ordonnee,col="light blue",lwd=.5)&#62; points(PP(X,F=F0))&#62; abline(a=0,b=1,col=rgb(1,0,0,.4))Cette fois, commes les invariants dans les graphiques sont sur l'axe des ordonnées, la recherche de quantile se fait horizontalement (ce qui est peu commun)&#62; Q1=apply(MVPP,2,function(x) quantile(x,.025))&#62; Q2=apply(MVPP,2,function(x) quantile(x,.975))&#62; lines(Q1,ordonnee,col="red")&#62; lines(Q2,ordonnee,col="red")On peut visualiser des histoires de points hors de la région d'acceptation sur le graphique ci-dessous (en insistant sur la lecture horizontale des quantiles)&#62; plot(PP(X,F=F0))&#62; polygon(c(Q1,rev(Q2)),c(ordonnee,rev(ordonnee)), border=NA,col="yellow")&#62; id=((PP(X,F=F0)$x&#62;=Q1)&#38;(PP(X,F=F0)$x&#60;=Q2))+1&#62; points(PP(X,F=F0),col=c('red','blue')[id],pch=c(19,1)[id])avec le fait que les points hors de la région d'acceptation se voient horizontalement ici.Un autre outils, plus formel, est le test de Kolmogorov-Smirnov. Cet dernier est basé sur la fonction de répartition empirique. On peut tracer cette dernière à l'aide de la fonction suivante&#62; Femp = function(Y){+   n=length(Y)+   x=sort(Y)+   y=seq(1/n/2,1-1/n/2,by=1/n)+   return(data.frame(x=x,y=y))+ }&#62; plot(Femp(X),type="s",xlim=c(150,200))Une alternative est d'uilitser la fonction suivante&#62; F=Vectorize(function(x) mean((X&#60;=x)))&#62; x=seq(150,200,by=.5)&#62; plot(x,F(x),type="s",ylab=" Probabilities")Notons que ce graphique est très proche des PP et QQ plot&#62; plot(Femp(X)$x,Q0(Femp(X)$y),type="s")En transformant par l'inverse de  en ordonnées, on retrouve le QQ plot, et en transformant par  en abscisse, on retrouve le PP plot,&#62; plot(F0(Femp(X)$x),Femp(X)$y,type="s")Là encore, on peut faire des simulations, pour voir où se trouverait la fonction de répartition, si on simulait des échantillons suivant &#62; MVF=matrix(NA,500,length(x))&#62; for(s in 1:500){+   Xs=Q0(runif(length(X)))+   Fs=Vectorize(function(x) mean((Xs&#60;=x)))+   MVF[s,]=Fs(x)+ }On peeut superposer les courbes&#62; plot(x,F(x),type="s",xlab="Values",ylab="Cumulated Probabilities",col="white")&#62; for(s in 1:500) lines(x,MVF[s,],col="light blue",lwd=.5,type="s")&#62; lines(x,F(x),type="s")&#62; lines(x,F0(x),type="l",col=rgb(1,0,0,.5))La distance utilisée pour le test est iciComme nous l'avions vu, cette distance peut être vu comme un maximum, en cherchant la valeur maximale des&#62; D=function(Y,F0){+   F=Vectorize(function(x) mean((Y&#60;=x)))+   n=length(Y)+   x=sort(Y)+   d1=abs(F(x+1e-6)-F0(x+1e-6))+   d2=abs(F(x-1e-6)-F0(x-1e-6))+   return(max(c(d1,d2)))+ }Ici, la distance entre la fonction de répartition et la loi normale que l'on cherche à tester est&#62; D(X,F0)[1] 0.06296857On retrouve cette statistique à l'aide de la fonction&#62; ks.test(X,F0)One-sample Kolmogorov-Smirnov testdata:  XD = 0.062969, p-value = 0.406alternative hypothesis: two-sidedPour comprendre d'où sort la p-value, là encore, on peut utiliser nos simulations&#62; VD1=rep(NA,500)&#62; for(s in 1:500){+   Xs=rnorm(length(X),m0,s0)+   VD1[s]=D(Xs,function(x) pnorm(x,m0,s0))+ }Dans 38% des scénarios, la distance dépasse la distance observée (0.0629)&#62; mean(D(X,F0)&#60;=VD1)[1] 0.38ce qui est comparable avec les 40% donnés ici. On peut visualiser cette p-value ci-dessous&#62; dst=density(VD1)&#62; plot(dst,col="darkblue")&#62; id=which(dst$x&#62;D(X,F0))&#62; polygon(c(dst$x[id],rev(dst$x[id])),c(dst$y[id],   + rep(0,length(id))),border=NA,col=rgb(0,0,1,.4))&#62; abline(v=D(X,F0),lty=2)