QUELQUES ELEMENTS DE CORRECTION DU TP1 SUR L'ESTIMATION

1.1 Simulations: l'EMV

 1. > curve(dexp(x), xlim=c(0,10), col='red', lwd=1, main='Densité de la loi exponentielle')
 2. > X <- rexp(20,2) ; hist(X,probability=T, main='Loi exponentielle'); u<- seq(0,5,by=0.1); lines(u,dexp(u,2),col='red',lwd=3)
 3. > V <- function(theta) {prod(dexp(X,theta)) } # la vraisemblance ou
    > LV <- function(theta) {sum(log(dexp(X,theta))) } # la log-vraisemblance
 5. > T<-seq(0.5,4,by=0.01) ; MV<-T; for(i in 1 :length(T)){ MV[i]<-V(T[i])}; plot(T,MV,main="Vraisemblance de l’échantillon")
    > LV <- function(theta){-sum(log(dexp(X,theta)))}; optim(2,LV)
 6. on utilise LV plutôt que V pr des histoires de pb numériques !!
 7. > EMV<-matrix(0,200,1); for(k in 1 :200){X <- rexp(20,2); LV <- function(theta) {-sum(log(dexp(X,theta)))} ; EMV[k]<-optim(2,LV)$par } ; plot(density(EMV))

1.2 Médiane

 1. > X<-rexp(100,.5)
 2. > M<-median(X); m=qexp(.5)
 3. > n=100 ; B=150 ; Mv<-matrix(0,B,1) Z <- Mv, for(i in 1 :B){Mv[i]<-median(rexp(n,.5)) ; Z[i]<-sqrt(n)*(Mv[i]-m)}; var(Mv)
    > hist(Z, probability=T, main=’Loi de l’estimateur de la médiane’) ; u<-seq(-10,10,by=0.1) ; lines(u,dnorm(u), col=’red’, lwd=1, lty=2) ; lines(u,dnorm(u,0,sqrt(var(Z))), col=’red’, lwd=2, lty=1)
 4. > n<-100 ; B<-150 ; Mv=Z=matrix(0,B,1) ; for(i in 1 :B){Mv[i]<-median(rexp(n,.5))}; T<-sqrt(n)*(Mv-mean(Mv))/sd(Mv)

2 Données réelles

 1. > library(MASS) ; data(quine) ; ?quine
    > base<-data(quine); obs<-base$Days
 2. > hist(obs) # alors loi exponentielle ? De Poisson (avantage : c'est du discret)
 3. > V <- function(theta) {prod(dpois(obs,theta)) } # ou LV avec l'analogue de l'exo 1.1
    > optim(-LV,20)
 4. > w<-hist(obs,probability=T) ; for(i in 1 :(length(w$breaks)-1){h<-(ppois(w$breaks[i+1],m)-ppois(w$breaks[i],m))/(w$breaks[i+1]-w$breaks[i]); segments(w$breaks[i],h,w$breaks[i+1],h,col=’red’,lwd=2)} # option 1 avec une loi de Poisson; on peut recommencer la même chose avec la loi exponentielle par exemple, ça marche bien mieux...
    > u<-seq(0.25,0.975,by=0.25); plot(quantile(obs,u),qexp(u,m)); segments(0,0,100,100,col=’blue’) #

3 Estimateurs pour la loi uniforme Unif([0;a])

 > theta0=2 ; X=runif(20,min=0,max=theta0)
 > m1=2*mean(X); m2=n/(n-1)*max(X) # puis comparer les variances empiriques de ces estimateurs. Plus généralement leur distribution.

4.1 Intervalles de confiance exacts (donnéees gaussiennes)

 > n<-100 ;B<-2000; X<-matrix( rnorm(n*B,2,1), n,B); Xbarre<-matrix(0,1,B) ; VX<-matrix(0,1,B)
 > for(i in 1 :B) {Xbarre[i]>-mean(X[,i]) ; VX[i]=var(X[,i])}
 > IC<-qt(0.95,df=n)*sqrt(VX/n); Ibinf<-Xbarre-IC ; Ibsup<-Xbarre+IC
 > m<-2 ; I<-(m>Ibinf)*(m<Ibsup) ; alpha<-sum(I)
 > plot(Xbarre,1 :B, xlim=c(1,3)) ; for(i in 1 :B){segments(Ibinf[i],i, Ibsup[i],i)} ; segments(m,1,m,B,col=“red”)

 > X<-rnorm(100,2,1) ; Xb<-mean(X) ; VX<-var(X) ; plot(matrix(Xb,10,1),1:10) ; for(i in 1:10){a<-qnorm(1-i/40); segments(Xb-a*sqrt(VX/100), i,Xb+a*sqrt(VX/100),i,col="blue",lwd=2)} ; segments(2,1,2,10,col=“red”)

4.2 IC approchés (données non gaussiennes)

 petite aide: si Xi~Exp(t) pr tt i, alors E[Xi]=1/t et Var(Xi)=1/t^2.
 Si on utilise le Th Central Limit pr avoir une loi asymptotique sur la moyenne, on en déduit un intervalle de la forme [Xbar +/- \phi^{-1}(0.95).Sbar/sqrt(n)] où s est l'écart type des Xi

4.3 Influence de n

 on trace les IC (asymptotiques) à 95% en coupant la pop en 2. Visiblement les données ne sont pas indépendantes (on pourra faire un test dégalité de moyenne sur ces 2 groupes qd on aura vu les tests). Une idée pourrait être de rééchantilloner. Les données deviendraient :
 > RANG<-rank(runif(length(obs))) ; obs=obs[RANG]