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 à .
> Davis=read.table(
+ "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt")
> Davis[12,c(2,3)]=Davis[12,c(3,2)]
> X=Davis$height
> 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
> 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 points
et celui pour le QQ plot serait
> 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 points
où 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 est
Nous 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
> m0=170
> s0=9
> F0=function(x) pnorm(x,m0,s0)
> Q0=function(x) qnorm(x,m0,s0)
Le graphique probabilité-probabilité est ici
> plot(PP(X,F=F0),xlab="Theoretical Probabilities",
+ ylab="Empirical Probabilities")
> abline(a=0,b=1,col=rgb(1,0,0,.4))
alors que le graphique quantile-quantile est
> plot(QQ(X,Q=Q0),xlab="Theoretical Quantiles",
+ ylab="Empirical Quantiles")
> 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
> MVPP=MVQQ=matrix(NA,500,length(X))
> for(s in 1:500) MVPP[s,]=PP(Q0(runif(length(X))),F0)$x
> for(s in 1:500) MVQQ[s,]=QQ(Q0(runif(length(X))),Q0)$y
Commençons par le graphique quantile-quantile
> plot(QQ(X,Q=Q0),xlab="Theoretical Quantiles",
+ ylab="Empirical Quantiles",col="white")
> abscisse=Q0(seq(1/n/2,1-1/n/2,by=1/n))
> for(s in 1:500) lines(abscisse,MVQQ[s,],col="light blue",lwd=.5)
> points(QQ(X,Q=Q0))
> abline(a=0,b=1,col=rgb(1,0,0,.4))
On peut calculer les quantiles empiriques, vus verticalement
> Q1=apply(MVQQ,2,function(x) quantile(x,.025))
> Q2=apply(MVQQ,2,function(x) quantile(x,.975))
> lines(abscisse,Q1,col="red")
> 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
> plot(QQ(X,Q=Q0),xlab="Theoretical Quantiles",
+ ylab="Empirical Quantiles",col="white")
> polygon(c(abscisse,rev(abscisse)),c(Q1,rev(Q2)),border=NA,col="yellow")
> id=((QQ(X,Q=Q0)$y>=Q1)&(QQ(X,Q=Q0)$y<=Q2))+1
> 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
> plot(PP(X,F=F0),xlab="Theoretical Probabilities",
+ ylab="Empirical Probabilities",col="white")
> ordonnee=seq(1/n/2,1-1/n/2,by=1/n)
> for(s in 1:500) lines(MVPP[s,],ordonnee,col="light blue",lwd=.5)
> points(PP(X,F=F0))
> 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)
> Q1=apply(MVPP,2,function(x) quantile(x,.025))
> Q2=apply(MVPP,2,function(x) quantile(x,.975))
> lines(Q1,ordonnee,col="red")
> 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)
> plot(PP(X,F=F0))
> polygon(c(Q1,rev(Q2)),c(ordonnee,rev(ordonnee)), border=NA,col="yellow")
> id=((PP(X,F=F0)$x>=Q1)&(PP(X,F=F0)$x<=Q2))+1
> 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
> 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))
+ }
> plot(Femp(X),type="s",xlim=c(150,200))
Une alternative est d’uilitser la fonction suivante
> F=Vectorize(function(x) mean((X<=x)))
> x=seq(150,200,by=.5)
> plot(x,F(x),type="s",ylab=" Probabilities")
Notons que ce graphique est très proche des PP et QQ plot
> 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,
> 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
> MVF=matrix(NA,500,length(x))
> for(s in 1:500){
+ Xs=Q0(runif(length(X)))
+ Fs=Vectorize(function(x) mean((Xs<=x)))
+ MVF[s,]=Fs(x)
+ }
On peeut superposer les courbes
> plot(x,F(x),type="s",xlab="Values",ylab="Cumulated Probabilities",col="white")
> for(s in 1:500) lines(x,MVF[s,],col="light blue",lwd=.5,type="s")
> lines(x,F(x),type="s")
> lines(x,F0(x),type="l",col=rgb(1,0,0,.5))
La distance utilisée pour le test est ici
Comme nous l’avions vu, cette distance peut être vu comme un maximum, en cherchant la valeur maximale des
> D=function(Y,F0){
+ F=Vectorize(function(x) mean((Y<=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
> D(X,F0)
[1] 0.06296857
On retrouve cette statistique à l’aide de la fonction
> ks.test(X,F0)
One-sample Kolmogorov-Smirnov test
data: X
D = 0.062969, p-value = 0.406
alternative hypothesis: two-sided
Pour comprendre d’où sort la p-value, là encore, on peut utiliser nos simulations
> VD1=rep(NA,500)
> 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)
> mean(D(X,F0)<=VD1)
[1] 0.38
ce qui est comparable avec les 40% donnés ici. On peut visualiser cette p-value ci-dessous
> dst=density(VD1)
> plot(dst,col="darkblue")
> id=which(dst$x>D(X,F0))
> 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))
> abline(v=D(X,F0),lty=2)
OpenEdition suggests that you cite this post as follows:
Arthur Charpentier (November 14, 2016). Tests non-paramétriques et simulations. Freakonometrics. Retrieved October 8, 2024 from https://doi.org/10.58079/ov5i
One thought on “Tests non-paramétriques et simulations”