PP plot, QQ plot et estimation à noyau

Reprenons ici le code pour récupérer la base de la taille (et du poids) d’un ensemble de 200 personnes.

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)

La première chose que l’on avait vu la semaine passée est le PP-plot,

https://latex.codecogs.com/gif.latex?\left\{\frac{i}{n},F_\star(X_{i:n})\right\}On va ici ‘tester’ une loi normale (de même moyenne et de même variance que notre échantillon… faisons simple)

PP=data.frame(empirique=(1:n)/n,
theorique=pnorm(sort(X),mean(X),sd(X)))

Comme dans tous les posts, je mets de lignes de code basées sur
le package ggplot2 (alors qu’en cours, j’utilise les fonctions de
bases de graphiques). Le code le plus simple pour tracer le PP plot est

library(ggplot2)
qplot(empirique,theorique,data=PP)

ou, pour un code qui pourra plus facilement être modifié (pour changer les couleurs, les titres, ou rajouter des courbes)

ggplot(PP,aes(x=empirique,y=theorique)) +
geom_abline(intercept=0,slope=1,colour="red") + 
geom_point()

La seconde chose que l’on avait vu la semaine passée est le QQ-plot,

https://latex.codecogs.com/gif.latex?\left\{F_\star^{-1}\left(\frac{i}{n}\right),X_{i:n}\right\}

Le point de départ est la table, avec les coordonnées, abscisse et ordonnnée,

QQ=data.frame(empirique=sort(X),
theorique=qnorm((1:n)/n,mean(X),sd(X)))

puis on fait un graphique

ggplot(QQ,aes(x=empirique,y=theorique)) +
geom_abline(intercept=0,slope=1,colour="red") +
geom_point()

Ici, on utilise le quantile de la loi normale, de même moyenne et même variance que notre échantillon. Mais comme on l’avait vu en cours, on peut en fait mettre une loi normale quelconque (disons centrée et réduite), et on a un test de normalité,

QQ=data.frame(empirique=sort(X),
theorique=qnorm((1:n)/n))

que l’on visualise avec

ggplot(QQ,aes(x=theorique,y=empirique)) +
geom_abline(intercept = mean(X), slope = sd(X),colour="red") +
geom_point()

Comme je l’avais dit en cours, on a ici une méthode simple d’ajustement de loi, pour tester

https://latex.codecogs.com/gif.latex?H_0:F=F_\star contre https://latex.codecogs.com/gif.latex?H_1:F\neq%20F_\star

essentiellement visuelle.

Si on veut aller un peu plus loin, en particulier pour avoir un
intervalle de confiance, on peut faire des simulations. Par exemple,
si la loi que l’on cherche est une loi https://latex.codecogs.com/gif.latex?\mathcal{N}(0,1).

PPs=NULL
for(s in 1:100){PPs=rbind(PPs,data.frame("empirique"=(1:n)/n,"theorique"=pnorm(sort(rnorm(n))),sample=s))}
ggplot(PPs,aes(x=empirique,y=theorique,group=sample))+
geom_line(colour="light blue")+geom_abline(intercept = 0, slope = 1,colour="red")

pour le PP plot, et pour le QQ plot

QQs=NULL
for(s in 1:100){QQs=rbind(QQs,data.frame("empirique"=sort(rnorm(n)),"theorique"=qnorm((1:n)/n),sample=s))}
ggplot(QQs,aes(x=empirique,y=theorique,group=sample))+
geom_line(colour="light blue")+
geom_abline(intercept = 0, slope = 1,colour="red")

Dans la second partie du cours, on avait parlé d’estimation à noyau de la densité.Classiquement, tout le monde connaît l’histogramme.

dX=data.frame(X=X)
H=ggplot(dX,aes(X))
H+geom_histogram(aes(y=..density..),fill="light blue")

On peut considérer plusieurs formes de noyaux. Un noyau rectangulaire

https://latex.codecogs.com/gif.latex?K(u)%20=%20\frac{1}{2}\%20\boldsymbol1_{(|u|\leq1)}

qui donnera un histogramme glissant. Le code est ici

H+geom_density(kernel="rectangular",fill="light blue",colour=NA)


On peut considérer un noyau triangulaire

https://latex.codecogs.com/gif.latex?K(u)%20=%20(1-|u|)\%20\boldsymbol1_{(|u|\leq1)}

H+geom_density(kernel="triangular",fill="light blue",colour=NA)

un noyau d’Epanechnikov

https://latex.codecogs.com/gif.latex?K(u)%20=%20\frac{3}{4}(1-u^2)\%20\boldsymbol1_{(|u|\leq1)}

H+geom_density(kernel="epanechnikov",fill="light blue",colour=NA)

ou un noyau Gaussien

https://latex.codecogs.com/gif.latex?K(u)%20=%20\frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}u^2}

H+geom_density(kernel="gaussian",fill="light blue",colour=NA)

 

Ces deux derniers donnant des visualisation très très proches. On peut ensuite jouer sur la fenêtre. Pour l’instant, on avait

H+geom_density(kernel="gaussian",adjust=1,fill="light blue",colour=NA)

mais on peut prendre une fenêtre plus petite

H+geom_density(kernel="gaussian",adjust=.25,fill="light blue",colour=NA)

Ce qui va augmenter la variance, mais baisser le biais (on va en reparler cette semaine, promis).


Leave a Reply

Your email address will not be published. Required fields are marked *