Tag Archives: normalité

C’est normal ! (partie 2) et si la normalité n’existait pas?

Avec un peu de retard, je vais tenter de poursuivre le billet sur la normalité, mis en ligne en avril dernier. Mon point de départ reprenait une discussion autour de la « moyenne » comme outil central pour définir la norme, ainsi que des idées que l’on peut retrouver dans le livre de Todd Rose, The End of Averages. Pour compléter mon précédant billet, je pourrais reprendre un exemple mentionné dans ce livre, revenant sur des problèmes rencontrés par l’Armée Américaine dans les années 1950. Lors de la conception des postes de pilotage des avions de chasse, les ingénieurs avaient utilisé les dimensions de plus de 4000 pilotes pour positionner de manière optimale le siège par rapport aux pédales, le manche à balai, la hauteur du pare-brise, mais aussi la forme du siège, du casque, etc. Ces mesures ont permis de calculer les mensurations du pilote « médian » dans une dizaine de dimensions. Par exemple la taille moyenne des pilotes était de 179cm, ce qui a permis de définir la taille d’un pilote moyen entre 175 et 185cm. Si une majorité des pilotes est de taille moyenne, parmi les 4000 pilotes, aucun (parmi les 4000 !) n’était « moyen » dans toutes les dimensions. Comme l’affirmait Daniels (1952), « concevoir un poste de pilotage pour le pilote moyen revenait en fait à n’en concevoir pour aucun ». Je faisais ensuite le lien avec le concept de malédiction de la dimension. Mais je voulais revenir aujourd’hui sur le sens (plus philosophique) de la « normalité ».

De la normalité statistique à la norme, de la régularité à la règle

D’un point de vue empirique, descriptif, être dans la norme, ne signifie rien d’autre qu’être dans la moyenne, ne pas se détacher trop de cette moyenne. On aura alors tendance à définir la norme comme la fréquence de ce qui se produit le plus souvent, comme l’attitude la plus fréquemment rencontrée ou de la préférence la plus régulièrement manifestée. Mais cette normalité n’est pas la normativité, et « être dans la norme », être exemplaire, relève alors d’une dimension différente, qui se rattache cette fois non plus à une description du réel mais à une identification de ce vers quoi il doit tendre. On passe donc du registre de l’être à celui du devoir-être, du « is » au « ought » pour reprendre la terminologie de Hume (1739). Il est en effet difficile d’envisager le modèle (ou la normalité) sans glisser vers ce second sens qu’on peut trouver au concept de norme, et qui déploie quant à lui une dimension à proprement parler normative. Cette vision mène à une confusion entre normes et lois, même si toute la normativité n’est pas épuisée par les lois. Hume constate ainsi que, dans tous les systèmes de morale, les auteurs passent de constat de faits, c’est-à-dire énonciatifs de type « il y a », à des propositions qui comportent une expression normative, comme « il faut », « on doit ». Ce que Hume conteste, c’est le passage d’un type d’affirmation à un autre : pour lui, ce sont là deux types d’énoncés qui n’ont rien à voir les uns avec les autres, et qu’on ne peut donc pas enchaîner logiquement les uns avec les autres, en particulier d’une norme empirique vers une règle normative. Pour Hume, une affirmation qui ne serait pas normatives ne peut pas donner lieu à une conclusion de type normatif. Cette affirmation de Hume a suscité nombre de commentaires et d’interprétations, en particulier parce qu’énoncé telle quel, elle semble être un obstacle à toute tentative de naturalisation de la morale (comme le détaillent Mc Intyre (1959) ou Rescher (1990)). En ce sens, on trouve la distinction forte entre la norme dans la régularité (normalité) et la règle (normativité).

Les lois statistiques, du micro au macro

La loi statistique relève de ce qui est « est » parce qu’on l’a observé ainsi (par exemple « les hommes sont plus grands que les chiens »). La loi humaine (divine, ou judiciaire) relève de ce qui est « est » parce qu’on l’a décrété, et donc « doit être » (« les hommes sont libres et égaux » ou « l’homme est bon »). Enfin, la loi physique relève de ce qui « est » parce qu’on peut le montrer (« les planètes sont attirées entre eux »), souvent dans le cadre d’hypothèses. On voit que les trois concepts peuvent être liés. Par exemple la « loi de Kepler » a été établie historiquement à l’aide d’observations (et rentrait historiquement dans la première catégorie), avant d’être démontrée dans le cadre du modèle Copernicien (et de passer alors dans la troisième). On peut d’ailleurs associer un concept d’équilibre à cette loi, cette « norme ». Toutefois, comme le souligne Hilpinen (1971), les lois probabilistes posent toutefois bon nombre de questions, il suffit de penser aux lancers de dés, ou aux attentes : que signifie « il est normal d’attendre 5 minutes le bus à l’arrêt », ou plus dérangeant d’un point de vue éthique, « il est normal qu’une personne placée en détention provisoire soit emprisonnée 18 mois » ?

On peut voir la « norme » comme une régularité de cas, observée à l’aide de fréquences (ou de moyennes), par exemple, sur la taille des individus, la durée du sommeil, autrement dit des données qui constituent la description d’individus. Les données anthropométriques ont ainsi permis de définir par exemple une taille moyenne des individus dans une population donnée, en fonction de leur âge ; par rapport à cette taille moyenne, un écart de 20% en plus ou en moins détermine le gigantisme ou le nanisme. Si l’on pense aux accidents de la route, il peut être considéré comme « anormal » d’avoir un accident de la route une année donné, à un niveau individuel (micro), car la majorité des conducteurs n’ont pas d’accident. Néanmoins, du point de vue de l’assureur (macro), la norme est que 10 % des conducteurs aient un accident. Il serait donc anormal que personne n’ait d’accident. C’est l’argument que l’on retrouve dans Durkheim (1897). De l’acte singulier qu’est le suicide, s’il est considéré du point de vue de l’individu qui le commet, Durkheim tente de le voir comme un acte social, relevant alors d’une réelle régularité, au sein d’une société donnée. Dès lors, le suicide devient, selon Durkheim, un phénomène « normal ». Les statistiques permettent alors de quantifier la tendance au suicide dans une société donnée, dès lors que l’on observe non plus l’irrégularité qui apparaît dans la singularité d’une histoire individuelle, mais une « normalité sociale » du suicide.

Norme, convention et aspects éthiques

Si on prend un point de vue évolutionniste, ce qui est « normal » est ce qui est le plus capable de s’adapter, de répondre aux besoins, de fournir un modèle pour la résolution des situations (la nature faisant disparaître l’anormalité), et la normalité tend vers la normativité, et il devient difficile de faire la distinction entre les deux aspects. En fait, David Hume aborde ce point dans l’exemple bien connu des rameurs, qui montent dans le même bateau pour traverser un fleuve et qui rament en cadence (cet exemple est longuement discuté dans Mackie (1980)). Les deux rameurs ajustent peu à peu leurs coups de rame, l’un par rapport à l’autre, et il n’est pas nécessaire d’obtenir un accord explicite (qui formulerait la « norme ») qu’ils respecteraient. La loi, qui consiste à imposer une norme peut être utile en cas de conflit (si un des rameurs refuse de ramer, ou deux rameurs de capacités physiques très différentes), mais bien souvent, il n’est pas nécessaire de formuler explicitement cette norme inhérente à leur conduite. L’observateur extérieur observera une régularité (lorsque le rythme de croisière sera atteinte) qu’il pourra modéliser, mais ce rythme « normal » observé n’est pas forcément imposé par une loi. Dans le cas de rameurs, on retrouve la notion d’équilibre évoquée précédemment. Construire un modèle, c’est extraire le signal du bruit (pour reprendre la distinction de Silver (2015)), c’est chercher une norme, au sens statistique. Mais cela va plus loin s’il l’on construit un modèle prédictif, la réalité devant alors se conformer au modèle, comme l’espèrent souvent les économètres (et autre modélisateurs).

Sondages et prévisions de séries temporelles

Ce soir, @imparibus proposait sur son blog un billet passionnant sur l’élection présidentielle en France, avec des graphiques superbes, basé sur un lissage temporel des résultats des sondages pour le premier tour de l’élection présidentielle en France (qu’il a fait sous excel, on saluera la performance !). En reprenant les données du site http://www.lemonde.fr/, on peut obtenir la même chose assez simplement (j’ai repris – ou presque – les couleurs utilisées sur le site, François Hollande en rose, Nicolas Sarkozy en bleu, François Bayrou en orange, et Marine Le Pen en noir…). Pour commencer, le code en R pour importer les données et les manipuler est le suivant

sondage=read.table("http://freakonometrics.blog.free.fr/
public/data/sondage1ertour2012.csv",
header=TRUE,sep=";",dec=",")
sondage[36,1]="11/03/12"
sondeur=unique(sondage$SONDEUR)
sondage$date=as.Date(as.character(sondage$DATE),
"%d/%m/%y")

(je fais manuellement une correction de date car je me suis trompé en saisissant les chiffres à la main). A partir de là, on peut reprendre l’idée de faire une régression locale pour trouver une tendance,

CL=brewer.pal(6, "RdBu")
couleur=c(CL[1],CL[6],"orange","grey")
datefin=as.Date("2012/04/12")
k=1
plot(sondage$date,sondage[,k+2],col=couleur[k],
pch=19,cex=.7,xlab="",ylab="",ylim=c(0,40),
xlim=c(min(sondage$date),datefin))
rl=lowess(sondage$date,sondage[,k+2])
lines(rl,lwd=2,col=couleur[k])
for(k in 2:4){
points(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7)
rl=lowess(sondage$date,sondage[,k+2])
lines(rl,lwd=2,col=couleur[k])
}

Le code a l’air long mais il faut définir les couleurs, générer une fenêtre graphique, mettre les régressions locales dedans, etc. En tant que telle, la commande qui permet de faire le lissage est juste

rl=lowess(sondage$date,sondage[,k+2])

Bon, on peut d’ailleurs faire toutes sortes de lissages, avec des splines par exemple (ce que j’ai davantage tendance à utiliser),

D=data.frame(date=min(sondage$date)+0:148)
library(splines)
k=1
plot(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
prl=predict(rs,newdata=D)
prlse=sqrt(prl/100*(1-prl/100)/1000)*100
polygon(c(D$date,rev(D$date)),c(prl+2*prlse,
rev(prl-2*prlse)),col=CL[3],border=NA)
lines(D$date,prl,lwd=2,col=CL[1])
points(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7,)
for(k in 2:4){
points(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7)
}

avec, là encore, la commande suivante pour lisser

rs=lm(sondage[,k+2]~bs(date),data=sondage)

Cette fois le code est un peu plus long, parce que j’ai tracé un intervalle de confiance à 95% (intervalle de confiance classique, en supposant que 1000 personnes ont été interrogées, comme évoqué dans d’anciens billets),

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic10.gif

On a ici la courbe suivante

avec l’intervalle de confiance pour François Hollande, mais on peut faire la même chose pour Nicolas Sarkozy,

k=2
plot(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
prl=predict(rs,newdata=D)
prlse=sqrt(prl/100*(1-prl/100)/1000)*100
polygon(c(D$date,rev(D$date)),c(prl+2*prlse,
rev(prl-2*prlse)),col=CL[4],border=NA)
lines(D$date,prl,lwd=2,col=CL[6])
points(sondage$date,sondage[,k+2],col=
couleur[k],pch=19,cex=.7,)
for(k in c(1,3:4)){
points(sondage$date,sondage[,k+2],col=
couleur[k],pch=19,cex=.7)
}

Mais comme auparavant, on peut aussi visualiser les “tendances” des quatre candidats en tête,

Voilà pour le début. En fait, idéalement, on voudrait faire un peu de prévision… Le hic est que les code usuel pour faire de la prévision (avec des ARIMA, du lissage exponentiel ou tout autre modèle classique) nécessitent de travailler avec des séries temporelles, telles qu’elles sont classiquement définies, c’est à dire “régulièrement espacées dans le temps“.

Je ne vais pas commencer à réclamer plus de sondages, loin de moins cette idée, mais pour mon modèle, j’avoue que j’aurais préféré avoir plus de points. Beaucoup plus de points. Un sondage par jour aurait été idéal en fait…

Qu’à cela ne tienne, on va simuler des sondages. L’idée est simple: on a une tendance qui nous donne une probabilité jour par jour (c’est exactement ce qui a été calculé pour tracer les courbes), via

D=data.frame(date=min(sondage$date)+0:148)
prl=predict(rs,newdata=D)

On va ensuite utiliser une hypothèse de normalité multivariée sur nos 4 candidats, auquel j’en ai en rajouté un cinquième fictif (mais ça ne servait à rien) en utilisant l’analogue multivariée de l’expression suivante (évoqué dans le premier billet de l’année sur les élections)

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/ic-sondage01.gif

On peut alors générer, jour après jour, des sondages, en simulant des vecteurs Gaussiens. Je vais les générer indépendants les uns des autres, parce que c’est plus simple, et que cela ne me semble pas être une trop grosse hypothèse. Pour générer un ensemble de sondages, on utilise la fonction suivante

library(mnormt)
simulation=function(S){
proba=c(S,100-sum(S))/100
variance= -proba%*%t(proba)
diag(variance)=proba*(1-proba)
variance=variance/1000
return(rmnorm(1,proba[1:4],variance[1:4,1:4]))}
simulsondages=function(M){
Z=rep(NA,ncol(M))
for(i in 1:nrow(M)){
Z=rbind(Z,simulation(M[i,]))
}
return(Z[-1,])}
prediction4=matrix(NA,nrow(D),4)
for(k in 1:4){
rs=lm(sondage[,k+2]~bs(date),data=sondage)
prediction4[,k]=predict(rs,newdata=D)
}

Par exemple, si on la fait tourner une fois, on obtient le graphique suivant

set.seed(1)
S4=simulsondages(prediction4)
S100=100*S4
k=1
plot(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
for(k in 2:4){
points(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7)
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
}

(les courbes lissées sont celles obtenues sur les vrais sondages). Là on peut être heureux, parce qu’on a des vraies séries temporelles. On peut alors faire de la prévision, par exemple en faisant du lissage exponentiel (optimisé),

library(forecast)
k=1
X=S100[,k]
ETS=ets(X)
F=forecast(ETS,h=60)
fdate=max(D$date)+1:60
k=1
plot(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
for(k in 2:4){
points(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7)
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),lwd=2,col=couleur[k])
}
polygon(c(fdate,rev(fdate)),c(as.numeric(F$lower[,2]),
rev(as.numeric(F$upper[,2]))),col=CL[3],border=NA)
polygon(c(fdate,rev(fdate)),
c(as.numeric(F$lower[,1]),rev(as.numeric(F$upper[,1]))),
col=CL[2],border=NA)
lines(fdate,as.numeric(F$mean),lwd=2,col=CL[1])

Là encore le code peut paraître long, mais c’est surtout la partie associée au graphique qui prend de la place (par soucis purement esthétique, on passe du temps sur les codes graphiques depuis le début). Le code qui modélise la série, et qui la projette, est ici donné par les deux lignes suivantes

ETS=ets(X)
F=forecast(ETS,h=60)

Pour François Hollande, sur la simulation des sondages passés que l’on vient d’effectuer, on obtient la prévision suivante

On peut bien entendu faire une projection similaire pour Nicolas Sarkozy,

k=2
X=S100[,k]
ETS = ets(X)
F=forecast(ETS,h=60)
k=1
plot(D$date,S100[,k],col=couleur[k],pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
for(k in 2:4){
points(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7)
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
}
polygon(c(fdate,rev(fdate)),c(as.numeric(F$lower[,2]),
rev(as.numeric(F$upper[,2]))),col=CL[4],border=NA)
polygon(c(fdate,rev(fdate)),
c(as.numeric(F$lower[,1]),rev(as.numeric(F$upper[,1]))),
col=CL[5],border=NA)
lines(fdate,as.numeric(F$mean),lwd=2,col=CL[6])

On notera que je ne me suis pas trop fatigué sur la projection: autant sur la génération de sondages passés on a tenu compte de corrélation (i.e. si un candidat a un score élevé, ça se fera au détriment des autres, d’où la corrélation négative utilisée dans la simulation), autant ici, les projections sont faites de manière complétement indépendantes. En se fatiguant un peu plus, on pourrait se lancer dans un modèle vectoriel gaussien. Ou mieux, on pourrait travailler sur des processus de Dirichlet (surtout que R a un package dédié à ce genre de modèle) parce qu’on travaille depuis le début sur des taux, comme évoqué auparavant. Mais commençons par faire simple pour un premier billet rapide sur ce sujet.

On peut ensuite s’amuser à générer plusieurs jeux de sondages, et de lancer des prévisions dessus,

set.seed(1)
for(sim in 1:20){
Ssim=simulsondages(prediction4)
F1=forecast(ets(100*Ssim[,1]),h=60)
F2=forecast(ets(100*Ssim[,2]),h=60)
lines(fdate,as.numeric(F1$mean),col=CL[1])
lines(fdate,as.numeric(F2$mean),col=CL[6])}

Une fois qu’on a fait tout ça, on a presque fini. On peut regarder au 22 avril qui est en tête (voire, par scénario, calculer la probabilité qu’un des deux candidats soit en tête), c’est à dire au soir du 1er tour

set.seed(1)
VICTOIRE=rep(NA,1000)
for(sim in 1:1000){
Ssim=simulsondages(prediction4)
F1=forecast(ets(100*Ssim[,1]),h=60)
F2=forecast(ets(100*Ssim[,2]),h=60)
VICTOIRE[sim]=(F1$mean[30]>F2$mean[30])}
mean(VICTOIRE)

Et voilà. Ah oui, je n’ai pas laissé le résultat. Tout d’abord parce que je suis un statisticien, pas unprédicateur. Mais aussi parce que si ce genre de pronostic amuse des gens… ils n’ont qu’à se mettre à R pour faire tourner les codes !

La normalité ? avant tout une histoire de taille

Un billet rapide pour préciser le fond de ma pensée lorsque j’explique aux étudiants (qui travaillent sur leur devoir) que sur les bases de données avec beaucoup d’observations, il est très dur d’accepter l’hypothèse de normalité.

Pour ça, simulons des jeux de données: on a une seule variable explicative, un bruit gaussien, homoscédastique. Mais on va mal spécifier le modèle. En particulier, le vrai modèle ne sera pas linéaire. Voilà ce qu’on simule, et que l’on va chercher à ajuster.

n=50
X1=runif(n)
epsilon=rnorm(n)
Y=-2+8*X1-4*(X1<.2)+epsilon
base=data.frame(Y,X1)
reg=lm(Y~X1,data=base)

Autrement dit, les résidus (du modèle linéaire) ne devraient pas être gaussiens, car on n’ajuste pas le bon modèle. Maintenant, mon point est le suivant: avec peu de données, la plupart des tests de normalité ne vont pas rejeter l’hypothèse de normalité. Alors qu’avec beaucoup de données, l’hypothèse de normalité sera presque toujours rejetée.
Prenons par exemple le test Lillie, évoqué dans un vieux billet, qui est une adaptation du test de Kolmogorov-Smirnov quand on compare la distribution empirique avec la loi normale dont les paramètres ont été obtenues par maximum de vraisemblance. Sur le graphique ci-dessous, à gauche, on a en abscisse la taille de l’échantillon (qui va de 10 à 1000), et en ordonnée la p-value du test Lillie. Ou plus précisément, la courbe bleue est la région où s’est trouvée la p-value dans 90% des cas (la courbe supérieure étant le quantile à 95% de la p-value sur 10,000 échantillons simulés). La courbe rouge est la valeur moyenne de la p-value. A droite, on a – toujours en fonction de la taille de l’échantillon – la proportion des échantillons simulés pour lesquels la p-value était supérieure à 5%.

Autrement dit, plus la taille de l’échantillon est grande, plus on a de chance de rejeter (fort justement d’ailleurs) l’hypothèse de normalité. Et ce, uniquement sur un effet taille. Ou pour être plus rigoureux, sur les petits échantillons, on accepte beaucoup trop souvent (à tort cette fois) l’hypothèse de normalité.
On peut aussi regarder le test de Shapiro-Wilk

ou d’Anderson-Darling,

ou encore de Cramér-von Mises

Le code R pour faire ces tests est relativement simple,

X1=runif(n)
Y=-2+8*X1-4*(X1<.2)+rnorm(n)
base=data.frame(Y,X1)
reg=lm(Y~X1,data=base)
AD[n,i]=ad.test(residuals(reg))$p.value
SW[n,i]=shapiro.test(residuals(reg))$p.value
CV[n,i]=cvm.test(residuals(reg))$p.value
LL[n,i]=lillie.test(residuals(reg))$p.value

Pour ceux qui veulent aller plus loin, se cache derrière la notion de puissance du test…

Normalité des estimateurs dans une régression

Cette semaine en cours, on a évoqué la normalité des estimateurs dans une régression linéaire. En supposant la normalité des résidus, on a vu en cours que http://freakonometrics.blog.free.fr/public/perso5/ols-002.gif était un estimateur Gaussien. En particulier, chacun des estimateurs est alors Gaussien, au sens où http://freakonometrics.blog.free.fr/public/perso5/ols--0003.gif, ce qui peut se visualiser sur le graphique suivant (la constante est en abscisse, et la pente en ordonnée), avec un intervalle de confiance à 95%,

En fait, la variance étant inconnue (mais pouvant être estimée), si on remplace la variance des résidus par leur estimateur, la loi de l’estimateur est une loi de Student. De même pour l’estimateur de la pente de la droite de régression, http://freakonometrics.blog.free.fr/public/perso5/ols---ooo4.gif

Le fait que le couple soit également Gaussien fait que l’on peut construire non plus un intervalle de confiance (on n’est plus en dimension 1) mais une ellipse de confiance. On retient une forme elliptique car c’est la plus petite région dans laquelle on se trouvera avec une probabilité de 95% (comme discuté dans un vieux billet).

Mais pour mieux comprendre cette notion d’ellipse de confiance, le plus simple est d’avoir recours à du rééchantillonnage dans la base de données. En effet, comme on ne dispose que d’un jeu de données, on a un estimateur. Et la discussion sur la loi de notre estimateur est purement théorique. Pour visualiser la loi de notre estimateur, il faudrait des centaines, voire des milliers de bases similaires. Que l’on n’a pas. La solution est alors de tirer au hasard des points de notre échantillon (avec remise), i.e. de faire du bootstrap,

set.seed(1)
COEF=matrix(NA,10000,2)
for(s in 1:nrow(COEF)){
I=sample(1:nrow(cars),nrow(cars),replace=TRUE)
COEF[s,]=lm(dist~speed,data=cars[I,])$coefficients
}

Si on regarde ce qui se passe, tirage après tirage, on génère un paquet d’échantillon, et pour chaque échantillon, on ajuste une droite de régression.

http://freakonometrics.blog.free.fr/public/perso5/BOOOT.gif

Les distributions – sur tous ces échantillons – de la constante et de la pente semblent effectivement Gaussiennes,

hist(COEF[,1],col="light blue",prob=TRUE)
u=seq(min(COEF[,1]),max(COEF[,1]),length=500)
v=dnorm(u,mean(COEF[,1]),sd(COEF[,1]))
lines(u,v,lwd=3,col="red")
hist(COEF[,2],col="light blue",prob=TRUE)
u=seq(min(COEF[,2]),max(COEF[,2]),length=500)
v=dnorm(u,mean(COEF[,2]),sd(COEF[,2]))
lines(u,v,lwd=3,col="red")

pour la distribution de l’estimateur de la pente, alors que pour l’estimateur de la constante, on a la distribution suivante,

Si on regarde maintenant la loi jointe,

on retrouve une forme elliptique pour le nuage des points (et une forte corrélation négative entre les deux estimateurs). Mais si on creuse un peu

library(ellipse)
reg=lm(dist~speed,data=cars)
e=ellipse(reg)
plot(e,type="l",lwd=2)
polygon(e,col="light blue")
points(COEF,cex=.5)

l’ellipse ne coïncide pas (tout à fait) avec celle obtenu théoriquement. Ce qui laisse à penser que l’hypothèse de normalité des résidus est peut-être à revoir..,

TCL et approximation

Pour répondre (rapidement) à une autre question, sur l’intérêt du théorème central limite, il y tout d’abord l’approximation gaussienne pour les sommes. Pour rappel, le théorème central limite nous permet d’écrire, pour un échantillon i.i.d. de variance finie

http://freakonometrics.blog.free.fr/public/maths/TCL-approx-2.png

lorsque http://freakonometrics.blog.free.fr/public/maths/TCL-approx-8.png, où http://freakonometrics.blog.free.fr/public/maths/TCL-approx-10.png et http://freakonometrics.blog.free.fr/public/maths/TCL-approx-11.png. On peut aussi l’écrire non pas sur la moyenne empirique, mais sur la somme (c’est cette forme que l’on a démontrée),

http://freakonometrics.blog.free.fr/public/maths/TCL-approx-1.png

Pour parler d’approximations, rappelons que si http://freakonometrics.blog.free.fr/public/maths/TCL-approx-3.png, alors http://freakonometrics.blog.free.fr/public/maths/TCL-approx-4.png. Et de plus, si http://freakonometrics.blog.free.fr/public/maths/TCL-approx-6.png, alors http://freakonometrics.blog.free.fr/public/maths/TCL-approx-5.png. Bref, en bricolant un peu, on peut en déduire une approximation de la forme suivante

http://freakonometrics.blog.free.fr/public/maths/TCL-approx-7.png

pour http://freakonometrics.blog.free.fr/public/maths/TCL-approx-9.png suffisamment grand. Empiriquement, si on s’amuse à sommer des variables de Bernoulli indépendantes (car on connaît très exactement la loi de la somme) on obtient, pour la fonction de probabilité

http://freakonometrics.free.fr/blog/TCL-cv-1.gif

et pour la fonction de répartition

http://freakonometrics.free.fr/blog/TCL-cv-1.gif

Bref, on a une belle convergence vers la loi normale….

Mais au delà de ces approximations numériques (très utiles à une époque où les ordinateurs n’étaient pas aussi avancés qu’aujourd’hui), il y a les propriétés de normalité asymptotique de plusieurs estimateurs, avec des résultats de la forme

http://freakonometrics.blog.free.fr/public/maths/TCL-approx-12.png

Ceci permettra de faire des tests, ou de construire des intervalles de confiance.Mais nous verrons ça (très) bientôt…