Illustration de la moyenne empirique:
x1=runif(100)
x2=runif(100)
x3=runif(100)
s1=mean(x1) #moyenne empirique
s2=mean(x2)
s3=mean(x3)
s1
## [1] 0.5377859
s2
## [1] 0.5200323
s3
## [1] 0.4884676
On peut améliorer les choses en choisissant d'augmenter la taille de l'échantillon:
x1=runif(10000)
x2=runif(10000)
x3=runif(10000)
mean(x1)
## [1] 0.497331
mean(x2)
## [1] 0.5058757
mean(x3)
## [1] 0.4995492
On peut aussi combiner ces échantillons (attention à ce qu'ils soient indépendants!):
mean(c(x1,x2,x3))
## [1] 0.5009186
Dans cette partie on suppose, pour faire simple, que les X suivent une loi uniforme sur [0,1]. (On peut généraliser à d'autres lois en remplaçant runif par les générateurs correspondants (éventuellement “faits maison”).)
Pour \( N \) fixé, \( S_N \) est la moyenne empirique de N tirages de X. On cherche à observer la densité de \( S_N \), donc on va faire \( k \) tirages de \( S_N \) (composés chacun de \( N \) tirages de \( X \)).
On fixe N=1. La moyenne empirique est donc égale à l'unique observation. On doit donc observer la densité de X. Pour X uniforme, c'est donc un rectangle [0,1]x[0,1].
n=500000#pour observer une jolie densité, il faut bcp de tirages.
S1=runif(n) #N=1 donc chaque tirage de S1 est égal à un unique tirage de X.
plot(density(S1))
Pour \( N=2 \), \( S_N \) est la moyenne de 2 tirages \( X_1 \) et \( X_2 \). Si les \( X \) suivent une loi uniforme, alors \( S_2 \) suit une loi triangulaire :-)
S2=runif(n)+runif(n) #N=2 donc chaque tirage de S2 est égal la somme de 2 tirages de X (divisée par 2 pour normaliser).
S2=S2/2
plot(density(S2))
Pour \( N \) quelconque, \( S_N \) est la moyenne de \( N \) tirages \( X_1 \) à \( X_N \). Cette fois il faut faire une boucle, paramétrée par \( N \)
tcl<-function(n=30){
S=rep(0,n) #création d'un vecteur rempli de zéros (élément neutre pour la somme).
for (i in 1:n){
S=S+runif(n) #n tirages de X additionnés pour cha
}
S=S/n
plot(density(S))
}
On peut maintenant tracer les densités pour \( N \) quelconque.
Afin de conserver ces données et d'en faire une comparaison, on peut construire un dataframe avec toutes ces échantillons de \( S_N \) pour plusieurs valeurs de \( N \): (On peut consulter une version plus générale de cette fonction dans les sources des slides de A. Legrand, avec une boucle sur plusieurs lois)
gen_unif <- function(k=50000,N=c(1,2,5,10,30,60,100,120)) { #attention k et N ont un sens très différent. k sert à observer la densité de SN.
df=data.frame();
for(n in N) {
S=rep.int(0,k);
for(i in 1:n) {
S = S + runif(k)
}
S = S/n;
df=rbind(df,data.frame(N=n,SN=S));
}
df;
}
Génération du dataframe
fig=gen_unif(k=1000) #génération du dataframe avec la fonction gen_unif
Ensuite on trace les densité obtenues dans le dataframe:
#manipulation du dataframe:
#fig[fig$N==1,] #données pour N=1
#fig[fig$N==1,SN] #colonne des k valeurs de SN pour N=1
#observation des densités
par(mfrow=c(1,4))#figure divisée en 4 courbes
plot(density(fig[fig$N==1,"SN"]),ylim=c(0,5)) #trace densité observée de SN pour N=1
plot(density(fig[fig$N==2,"SN"]),ylim=c(0,5)) #trace densité observée de SN pour N=2
plot(density(fig[fig$N==5,"SN"]),ylim=c(0,5)) #trace densité observée de SN pour N=5
plot(density(fig[fig$N==10,"SN"]),ylim=c(0,5)) #trace densité observée de SN pour N=10
On continue avec de plus grandes valeurs de N
par(mfrow=c(1,4)) #figure divisée en 4 courbes
plot(density(fig[fig$N==1,"SN"]),ylim=c(0,15)) #trace densité observée de SN pour N=1
plot(density(fig[fig$N==10,"SN"]),ylim=c(0,15)) #trace densité observée de SN pour N=10
plot(density(fig[fig$N==30,"SN"]),ylim=c(0,15)) #trace densité observée de SN pour N=30
plot(density(fig[fig$N==100,"SN"]),ylim=c(0,15)) #trace densité observée de SN pour N=100
Observation de l'évolution de la variance
par(mfrow=c(1,2))
#plot(density(fig[fig$N==30,"SN"]),ylim=c(0,15),xlim=c(0.2,0.8)) #trace densité observée de SN pour N=30
plot(density(fig[fig$N==60,"SN"]),ylim=c(0,15),xlim=c(0.2,0.8)) #trace densité observée de SN pour N=60: on double
#par(new=TRUE)
plot(density(fig[fig$N==120,"SN"]),ylim=c(0,15),xlim=c(0.2,0.8)) #trace densité observée de SN pour N=120 on double à nouveau
lci<-function(x){#x est la colonne de valeurs sur laquelle calculer l'IC
n=length(x) #taille d'échantillon
s=qnorm(0.975) #le quantile de 0.975 - par symétrie de l'intervalle. Vaut environ 1.96
m=mean(x)
sd=sd(x) #ecart-type empirique : attention DANGER ! Même pour n grand, cela peut sous-estimer l'écart-type théorique, donc sous-estimer la largeur de l'IC. Quand c'est possible on préfèrera l'écart-type théorique:
sd=1/sqrt(12) #ecart-type théorique pour la loi uniforme sur [0,1]
largeur_CI=s*sd/sqrt(n)
c(m,largeur_CI,m-largeur_CI,m+largeur_CI) #cette fonction renvoie la moyenne, la largeur d'intervalle, puis les bornes basse et haute
}
Application pour calculer les intervalles de confiance observables avec un échantillon de taille n:
x=unique(fig$N) #ensemble des valeurs de N testées; ce sera aussi le vecteur d'abcisses du graphique.
#on calcule les IC autour de la première valeur mesurée; on restreint la taille de l'échantillon considéré aux N premières mesures (IC théorique que nous offre une campagne de N mesures)
#maintenant les intervalles de confiance pour chaque SN:
ci_sup=c() #largeur d'IC
ci_inf=c() #moyenne observée sur N valeurs
sn=c()#vecteur des moyennes empiriques observées pour différents N
k=1
for (i in x){
y=fig[fig$N==i,"SN"] #ensemble des mesures faites pour N=i
z=y[1:i] #mesures considérées (limitées à N=i valeurs)
u=lci(z)
ci_sup[k]=u[4]
ci_inf[k]=u[3]
sn[k]=u[1]
k=k+1
}
sn
## [1] 0.5548040 0.2232114 0.3998293 0.5571188 0.5013308 0.4979751 0.5066747
## [8] 0.4973085
ci_sup
## [1] 1.1205968 0.6232874 0.6528596 0.7360382 0.6046299 0.5710187 0.5632539
## [8] 0.5489581
ci_inf
## [1] -0.0109889 -0.1768646 0.1467990 0.3781994 0.3980316 0.4249316
## [7] 0.4500954 0.4456589
plot(x,ci_inf,type="p",pch='+',col='red',ylim=c(-0.5,1.5),main='',xlab='',ylab='')
par(new=TRUE)
plot(x,ci_sup,type="p",pch='+',col='red',ylim=c(-0.5,1.5),main="95% Confidence intervals for the sample mean",xlab='Nb of samples',ylab='95% confidence interval (with sample mean in blue)')
par(new=TRUE)
plot(x,sn,type="p",pch='+',col='blue',ylim=c(-0.5,1.5),main='',xlab='',ylab='',cex=1.5)
#lines(ci_inf,ci_sup)
for (i in 1:length(x)){
lines(c(x[i],x[i]),c(ci_sup[i],ci_inf[i]),lwd=2,lty=1,col='red')
}
axis(1,pos=0)
On observe que pour de très faibles valeurs de N (1 ou 2 ici) l'intervalle de confiance peut même dépasser les valeurs possibles (ici l'intervalle [0,1]) !