Jeudi et vendredi, on finira les rappels de cours, avec quelques rappels de statistique. Les slides sont en ligne,

Sinon, histoire d’illustrer un peu les différents aspects, j’avais dit que j’aimais bien utiliser des simulations. En particulier, il sera intéressant de générer plusieurs (beaucoup) échantillons de taille donnée. Le code R est le suivant (une ligne correspond à un échantillon, de 20 observations)

> n=20
> m=1000
> M=matrix(rnorm(n*m),m,n)

On utilise ensuite la statistique que l’on souhaite étudier, qui peut être la moyenne empirique,

> Xbarn=apply(M,1,mean)

ou la variance

> Sn=apply(M,1,var)

On peut ensuite visualiser la distribution de la statistique, sur tous les échantillons. Par exemple, le fait que l’estimateur de la moyenne est gaussien, on peut faire un histogramme

> hist(Xbarn,probability=TRUE,
+ breaks=seq(-1,1,by=.1),col="light blue")

et rajouter la densité de la loi normale (car la théorie nous dit que la moyenne suit une loi normale),

> u=seq(-2,2,by=.01)
> lines(u,dnorm(u,sd=1/sqrt(n)),col="red",lwd=2)

Cela dit, comme Ewen Gallic fera essentiellement du ggplot2 dans son cours de R, je vais aussi mettre les codes en ggplot2 (et essayer de m’y convertir… ce n’est pas l’envie qui manque !)

> ggplot(data=data.frame(Xbarn = Xbarn), aes(x = Xbarn)) +
+   geom_histogram(aes(y = ..density..),
+   binwidth = .1, fill = "light blue", col = "black") +
+   stat_function(fun=function(x) dnorm(x,sd=1/sqrt(n)),
+ col = "red", size = 1) + ylab("Densite") +
+   ggtitle("Histogramme de Xbarn")

Cela dit,  au lieu de représenter la loi normale, on peut aussi visualiser un estimateur non paramétrique de la densité (on en reparlera dans quelques semaines), comme l’estimateur à noyau

> ggplot(data=data.frame(Xbarn = Xbarn), aes(x = Xbarn)) +
+   geom_histogram(aes(y = ..density..),
+   binwidth = .1, fill = "light blue", col = "black") +
+   geom_line(stat="density", col = "red", size = 1) +
+   ylab("Densite") +
+   ggtitle("Histogramme de Xbarn")

On peut aussi voir que la variance empirique suit une loi du chi-deux,

> hist((n-1)*Sn,probability=TRUE,
+ breaks=seq(0,56,by=2),
+ col="light blue")
> u=seq(0,56,by=.1)
> lines(u,dchisq(u,df=n-1),col="red",lwd=2)

avec là encore une version en ggplot2

> ggplot(data = data.frame(x = (n-1)*Sn), aes(x = x)) +
+   geom_histogram(aes(y = ..density..),
+   binwidth = 2, fill = "light blue", col = "black") +
+   geom_line(stat="density", col = "red", size = 1) +
+   xlab(expression((n-1) %*% Sn)) +
+   ylab("Densite") +
+   ggtitle(expression(paste("Histogramme de ", (n-1) %*% Sn)))

On peut aussi utiliser ces simulations pour comprendre l’idée d’intervalle de confiance à 90%,

> x1=Xbarn-1.64*1/sqrt(n)
> x2=Xbarn+1.64*1/sqrt(n)
> ind=x1*x2>0
> plot(1:m,Xbarn,ylim=c(-1,1),col="white",
+ ylab="intervalle de confiance",xlab="")
> segments(1:m,x1,1:m,x2,col=c("light blue","red")[ind+1],lwd=1.5)
> abline(h=0,lty=2)

ou pour la version ggplot2 (mais on reparlera de tout ça en cours, pas forcément cette semaine cela dit)

> ggplot(data = df) +
+   # Les segments. La couleur doit etre fonction d'un logique
+   geom_segment(aes(x = x, xend = x, y = x1, yend = x2, col = x1*x2>0)) +
+   # La ligne horizontale en y = 0
+   geom_hline(yintercept = 0, linetype = "dashed") +
+   # Changer les couleurs des segments
+   scale_colour_manual(values = c("light blue" ,"red")) +
+   # Modifier le titre de l'axe des y
+   ylab("Intervalle de confiance") +
+   # Modifier le titre de l'axe des x
+   xlab("") +
+   # Cacher la legende
+   theme(legend.position="none")

On note que dans 90% des scénarios (notion d’intervalle de confiance à 90%), la vraie valeur (connue ici car on simule nos échantillons) est effectivement dans l’intervalle de confiance.