Aide Moyenne empirique avec logiciel R

[invite4252db6a]

Bonjour,

J'ai ce projet à réaliser à l'aide du logiciel R :

Loi de la moyenne empirique I (bootstrap).

On dispose de 1000 observations (urne initiale, à générer) d'une variable aléatoire suivant une loi normale de paramètres m et s. Tirer au hasard avec remise 1000 observations dans l'urne initiale.
Calculer les valeurs empiriques des paramètres m et s. Répéter cette procédure cinquante fois pour obtenir des estimateurs m i et s i , i=1..50.
Comparer la précision des valeurs empiriques des estimateurs basées sur l'échantillon de départ à celle des nouveaux estimateurs moyens définis par m : moyenne des m, s : moyenne des si .
Que se passerait-il si on réalisait la même procédure 500 fois plutôt que 50 ?
On donnera dans la présentation des précisions sur la notion de bootstrap.



J'ai essayé d'y répondre en réalisant ce programme, mais je pense qu'on peut l'améliorer notamment la partie où il faut comparer la précision des estimateurs à l'aide d'un graphe.
D'autre part, si quelqu'un a des notions sur le bootstrap, pourrait-il me faire part de ses connaissances ?
Voici mon programme :

projet=function(m,s) {
M=0
S=0
#pour générer l'urne de 1000 observations d'une v.a suivant une loi normale de
#paramètres m et s
x=rnorm(1000,m,s)
for(i in 1:50) {
#pour tirer au hasard avec remise 1000 observations dans l'urne initiale :
y=sample(x,1000,replace=TRUE)
#pour calculer les valeurs empiriques de m et s
mi=mean(y)
si=sd(y)

if(M[1]!=0) M=cbind(M,mi)
if(M[1]==0) M=mean(mi)
if(S[1]!=0) S=cbind(S,si)
if(S[1]==0) S=mean(si)

}

D=data.frame(M,S)
T=matrix(D,nrow=50,ncol=2)
#pour afficher le vecteur des mi
print("Les 50 estimateurs mi et si sont :")
print(T)
#pour afficher la moyenne de tous les mi
print("La moyenne des estimateurs mi est :")
print(mean(M))
#pour afficher la moyenne de tous les si
print("La moyenne des estimateurs si est :")
print(mean(S))


par(mfrow=c(2,1))
barplot(M,ylim=c(m-1,m+1),sub="50 estimateurs mi",col="blue")
abline(h=mean(M),col="red")

barplot(S,ylim=c(s-1,s+1),sub="50 estimateurs si",col="yellow")
abline(h=mean(S),col="blue")
}

Je vous remercie par avance pour votre aide.
Cordialement.
-----

[invite986312212]

salut,

avec R il faut essayer d'éviter les boucles "for".

Par exemple pour avoir 50 moyennes par bootstrap, tu peux faire:

u<-sample(x, 50000, replace=T)
u<-matrix(u,ncol=50)
moyennes<-apply(u,2,mean)

que tu peux d'ailleurs regrouper en une seule instruction.

[invitec5eb4b89]

Bonjour à tou-te-s,
OK c'est beaucoup trop tard pour répondre à cette question, mais ce prblème m'amusait, donc j'ai essayé de le coder en suivant les conseils éclairés d'ambrosio :

# Fonction de génération des b échantillons bootstrap
projet <- function(m=0,s=1,N=1000,b=50){
u <-matrix(sample(rnorm(N,m,s), N*b, replace=T),ncol=b)
list("moy" = apply(u,2,mean),"dev" = apply(u,2,sd))
}
# Estimation des moyenne et écart-type
m <- 0; s <- 1
boot50 <- projet(m=m,s=s)
moy50 <- boot50$moy ; dev50 <- boot50$dev
boot500 <- projet(m=m,s=s,b=500)
moy500 <- boot500$moy ; dev500 <- boot500$dev

# Affichage
maxigraf1 <- max(c(hist(moy50,freq=FALSE)$intensities,hist(moy500,freq=FALSE)$intensities))
maxigraf2 <- max(c(hist(dev50,freq=FALSE)$intensities,hist(dev500,freq=FALSE)$intensities))
graphics.off()
layout(1:2)
graf1 <- hist(c(moy50,moy500,m),ylim=range(0,maxigraf1),freq=FALSE,col=NA,border=NA,xlab="",ylab="",main="")
hist(moy50 ,add=TRUE,border=1,freq=FALSE,breaks=graf1$breaks)
hist(moy500,add=TRUE,border=2,freq=FALSE,breaks=graf1$breaks)
abline(v=m,lty=2,col=3)
graf2 <- hist(c(dev50,dev500,s),ylim=range(0,maxigraf2),freq=FALSE,col=NA,border=NA,xlab="",ylab="",main="")
hist(dev50 ,add=TRUE,border=1,freq=FALSE,breaks=graf2$breaks)
hist(dev500,add=TRUE,border=2,freq=FALSE,breaks=graf2$breaks)
abline(v=s,lty=2,col=3)

L'affichage des histogrammes me rend le plus de place parce que je ne suis pas très doué pour afficher plusieurs choses différentes sur le même graphe...
Encore désolé pour le retard à l'allumage,
V.