Ejercicio 1. Simualción bootstrap

Sea \(X_1,\dots,X_{43}\) i.i.d. tal que \(X_i\sim N(3,1), i=1,\dots,43.\) Considere \(Z=\overline{X}^3=(\frac{\sum_{i=1}^n X_i}{n})^3,\) y obtenga la varianza muestral de \(Z\):

  1. A partir de los conocimientos estudiados en IE1 e IEII ¿Es posible dar una expresión cerrada para calcular \(Var(Z)\) de forma exacta? ¿Es posible calcular \(Var(Z)\) de forma aproximada?

  2. Por simulación \(Var(Z).\)

# La simualcion permite aproximar cualquier caracteristica de la distribucion
n<-43
veces<-1000
y<-matrix(rnorm(n*veces,3,1),n,veces)
z<-apply(y,2,mean)
var(z^3)
## [1] 16.81317
par(mfrow=c(1,1))
hist(z^3)

  1. Si únicamente se dispone de la muestra \(X_1,\dots,X_{43}\) y no se tiene información sobre su distribución.
# Bootstrap no parametrico (original)
set.seed(345)
x<-rnorm(n,3,1)
B<-1000

bootP<-matrix(0,B,length(x))
zStar<-c()
for(b in 1:B){
  bootP[b,]<-sample(x,length(x),replace=TRUE)
  zStar<-c(zStar,mean(bootP[b,])^3)
}
mean((zStar-mean(zStar))^2)
## [1] 15.84146
((B-1)/B)*var(zStar)# R trabaja con la quasi varianza, la teoria esta para la varianza
## [1] 15.84146
  1. Si se dispone de la muestra \(X_1,\dots,X_{43}\) y además se sabe que su distribución es normal.
# Bootstrap  parametrico (original)
#x<-rnorm(n,3,1)
B<-1000
muHat<-mean(x)
varHat<-var(x)

bootNP<-matrix(0,B,length(x))
zStarNP<-c()
for(b in 1:B){
  bootNP[b,]<-rnorm(n,muHat,sqrt(varHat))
  zStarNP<-c(zStarNP,mean(bootNP[b,])^3)
}
mean((zStarNP-mean(zStarNP))^2)
## [1] 17.63455
((B-1)/B)*var(zStarNP)
## [1] 17.63455

e.- Considere \(W=(Me(X_1,\dots,X_n))^3\), obtener el estimador bootstrap para \(Sesgo(W)\). ¿Cómo es posible corroborar que el valor que obtenemos es correcto?

# Bootstrap no parametrico (original)
set.seed(9321258)
x<-rnorm(n,3,1)
B<-1000

bootNP<-matrix(0,B,length(x))
wStar<-c()
for(b in 1:B){
  bootNP[b,]<-sample(x,length(x),replace=TRUE)
  wStar<-c(wStar,median(bootNP[b,])^3)
}
mean(wStar)-median(x)^3
## [1] -0.2830232
mean(apply(y,2,median)^3)-3^3
## [1] 0.2517645

Ejercicio 2. Bootstrap. Modelo exponencial

Considerar una muestra de tamaño \(n=35\) de una distribución exponencial negativa de parámetro \(\lambda=0.5\) con función de densidad \(f(x,\lambda)=\lambda\exp^{-\lambda x}, x>0.\)

  1. Obtener el estimador bootstrap el sesgo y de la varianza del EMV de \(\lambda\) mediante bootstrap no paramétrico
# Verdadero valor del parametro
lambda<-0.5 
n<-35
x<-rexp(n,lambda)# datos
emv<-1/mean(x)#estimador

# Bootstrap no parametrico
B<-1000
mBootNP<-matrix(0,B,n)  
emvStar<-rep(0,B)       #las B versiones bootstrap del estimador
for(b in 1:B){
    mBootNP[b,]<-sample(x,n,replace=T)      
    emvStar[b]<-1/mean(mBootNP[b,])
}

# Estimador bootstrap de la varianza del EMV
mean((emvStar-mean(emvStar))^2)
## [1] 0.005845289
vboot<-((B-1)/B)*var(emvStar)       
vboot
## [1] 0.005845289
# Estimador bootstrap del sesgo del EMV
sesgoboot<-mean(emvStar)-emv    
sesgoboot
## [1] 0.0104515
  1. Obtener el estimador bootstrap el sesgo y de la varianza del EMV mediante bootstrap paramétrico
emv<-1/mean(x)#estimador

# Bootstrap no parametrico
B<-1000
mBootP<-matrix(0,B,n)   
emvStarP<-rep(0,B)      #las B versiones bootstrap del estimador
for(b in 1:B){
    mBootP[b,]<-rexp(n,emv)     
    emvStarP[b]<-1/mean(mBootP[b,])
}

# Estimador bootstrap de la varianza del EMV
mean((emvStarP-mean(emvStarP))^2)
## [1] 0.003932946
vbootP<-((B-1)/B)*var(emvStarP)     
vbootP
## [1] 0.003932946
# Estimador bootstrap del sesgo del EMV
sesgobootP<-mean(emvStarP)-emv  
sesgobootP 
## [1] 0.01121677
  1. Calcular mediante simulación, el sesgo y la varianza del estimador propuesto basado en una muestra de tamaño \(n=35\).
# Estos cálculos se pueden hacer exactos utilizando la distribución gamma 
# Exponencial(lambda~gama(1,lambda)

# Se obtiene: sesgo(emv)=lambda/(n-1) =0.0147  (para n=35, y lambda=0.5), desigualdad Jensen
# var(emv)=(lambda^2)*n^2/((n-1)^2)*(n-2)=0.0080    (para n=35, y lambda=0.5)


lambda<-0.5 
n<-35
x<-rexp(n,lambda)# datos
emv<-1/mean(x)#estimador


lambda<-0.5
veces<-1000
y<-matrix(rexp(veces*n,lambda),veces,n) #generar 1000 muestra de tamaño n=35 de la exponencial negativa 0.5.
lSim<-apply(y,1,mean)
lSim<-1/lSim            
varestimador<-((veces-1)/veces)*var(lSim)
sesgoestimador<-mean(lSim)-lambda


varestimador # (lambda^2)*n^2/(((n-1)^2)*(n-2))=0.008
## [1] 0.008057114
sesgoestimador # lambda/(n-1) =0.0147
## [1] 0.01386419
# Varianza
(lambda^2)*n^2/(((n-1)^2)*(n-2))
## [1] 0.008027944
# Sesgo
lambda/(n-1)
## [1] 0.01470588

Ejercicio 3. IC y CH Boostrap. Modelo normal

Considerar una muestra de tamaño \(n=30\) de una distribución normal estándar.

  1. Calcular IC (percentil) bootstrap paramétrico de nivel \(\alpha=0.05\) para la media poblacional. Compara con el IC estandar que utilizarías en este caso.
set.seed(12345)
x<-rnorm(30)
n<-30
B<-2000

mu<-mean(x)
ds<-sd(x)
Z<-matrix(rnorm(n*B,mu,ds),B,n)

muStar<-apply(Z,1,mean)


ext.inferior1<-quantile(muStar,0.025)
ext.superior1<-quantile(muStar,0.975)
ext.inferior1
##       2.5% 
## -0.2600999
ext.superior1
##    97.5% 
## 0.405254
#IC (t-Student) para la media es:
mu+c(-1,1)*qt(0.975,n-1)*ds/sqrt(n-1)
## [1] -0.2775214  0.4351354
# Distribucion (histograma)
par(mfrow=c(1,2))
hist(muStar)
abline(v=c(ext.inferior1,ext.superior1),lty=2)#dashed
abline(v=c(mu+c(-1,1)*qt(0.975,n-1)*ds/sqrt(n-1)),lty=3)#dotted

  1. Calcular IC (percentil) bootstrap no paramétrico de nivel \(\alpha=0.05\) para la media poblacional. Compara con el IC estandar que utilizarías en este caso.
Z<-matrix(sample(x,n*B,replace=TRUE),B,n)
T<-apply(Z,1,mean)

ext.inferior2<-quantile(T,0.025)
ext.superior2<-quantile(T,0.975)
ext.inferior2
##       2.5% 
## -0.2418123
ext.superior2
##     97.5% 
## 0.3981169
#IC (t-Student) para la media es:
mu+c(-1,1)*qt(0.975,n-1)*ds/sqrt(n-1)
## [1] -0.2775214  0.4351354
# Distribucion (histograma)
par(mfrow=c(1,2))
hist(T)
abline(v=c(ext.inferior2,ext.superior2),lty=2)
abline(v=c(mu+c(-1,1)*qt(0.975,n-1)*ds/sqrt(n-1)),lty=3)

  1. Considere \(\varphi=\mu^3+\mu^5\). Calcular IC (percentil) bootstrap no paramétrico de nivel \(\alpha=0.05\) para \(\varphi\). Compara con el IC estandar que utilizarías en este caso. ¿Por qué hay diferencias? ¿Cuál es mejor?
Z<-c()#marix(0,B,length(x))
for(j in 1:B){
  Z<-c(Z,(mean(sample(x,length(x),replace=TRUE)))^3+(mean(sample(x,length(x),replace=TRUE)))^5)
}

ext.inferior2<-quantile(Z,0.025)
ext.superior2<-quantile(Z,0.975)
ext.inferior2
##        2.5% 
## -0.01458542
ext.superior2
##      97.5% 
## 0.07262164
#IC (normal) para la media es:
mean(x)^3+mean(x)^5+c(-1,1)*qnorm(0.975)*sd(Z)
## [1] -0.04888045  0.04986540
# Distribucion (histograma)
par(mfrow=c(1,2))
hist(Z)
abline(v=c(ext.inferior2,ext.superior2),lty=2)
abline(v=c(mean(x)^3+mean(x)^5+c(-1,1)*qnorm(0.975)*sd(Z)),lty=3)

  1. A partir de los datos dados para el modelo normal. Obtener el valor crítico y el p-valor bootstrap para contrastar \(H_0:\mu=0\) \(H_1:\mu\neq 0\) (\(\sigma\) desconocida). Consider \(\alpha=0.05\).
# El estadístico test es T=mean(x)*sqrt(n)/sd(x). Región crítica (|T|>C). 
# Usando estimador de la quasivarianza. Si uso el de la varianza (como en teoria)
# sd(.)*(n-1)/n porque R trabaja siempre con la quasivarianza
# El valor crítico exacto para nivel 0.05 es: qt(0.975,29)=2.03
# El pvalor exacto es:  2*pt(-tobs,29), siendo tobs el valor observado de |T|

# Valor observado del estadiatico
tobs<-abs(mean(x)*sqrt(n)/sd(x))
tobs
## [1] 0.4600639
# Valor critico 
qt(0.975,n-1)
## [1] 2.04523
# Pvalor
2*pt(-tobs,n-1)
## [1] 0.6489013
# d1) Bootstrap parametrico
muB<-rep(0,B)
sdB<-rep(0,B)
tobsB<-rep(0,B)
for(b in 1:B){
    muestra<-rnorm(n,mu,ds)
    muB[b]<-mean(muestra)
    sdB[b]<-sd(muestra)
    tobsB[b]<-abs(muB[b]*sqrt(n)/sdB[b])
}
vcB<-quantile(tobsB,probs=0.95) # Podeis comprobarlo subiendo n y B en un bootstrap paramétrico que este valor es 1.96
vcB
##      95% 
## 2.222528
pvalorB1<-mean(tobsB>=tobs)
pvalorB1 # ¿Se deberia rechazar?
## [1] 0.6795
# d2) Bootstrap no parametrico
tobs<-abs(mu*sqrt(n)/ds)
muB<-rep(0,B)
sdB<-rep(0,B)
tobsB<-rep(0,B)
for(b in 1:B){
    muestra<-sample(x,length(x),replace=TRUE)
    muB[b]<-mean(muestra)
    sdB[b]<-sd(muestra)
    tobsB[b]<-abs(muB[b]*sqrt(n)/sdB[b])
}
vcB<-quantile(tobsB,0.95)
vcB
##      95% 
## 2.286952
pvalorB2<-mean(tobsB>=tobs)
pvalorB2 # ¿Deberia estar de acuerdo con el bootstrap parametrico?
## [1] 0.671

Ejercicio Propuesto 1. Ejercicio de entrega del procesador (mixtura de normales)

En un procesador se está ejecutando un proceso de simulación, en el que un procedimiento estadístico es repetido hasta completar un millón de repeticiones. En alguna de estas repeticiones el proceso falla, aunque se ha programado para que el proceso no se detenga y descarte dicha repetición sustituyéndola por otra.

Tras observar un total de \(n=272\) tiempos de fallo durante dos semanas, la experiencia del investigador le hace pensar que los fallos provienen de dos causas diferentes. Concretamente, afirma que los tiempos de espera hasta que se produce fallo (Y) pueden modelarse como una mixtura de dos normales, es decir que cada \(Y_i, i=1,\dots,272\) proviene bien de una distribución \(N(\mu,\sigma)\) con probabilidad \(p\), o bien de una distribución \(N(\nu,\tau)\) con probabilidad \(1-p\).

A partir de los datos de la muestra de tiempos tiempo de fallo facilitada (tiempoFallo.RData) y los valores iniciales \((p_0,\mu_0,\sigma_0,\nu_0,\tau_0)=(0.3,50,4,82,6)\):

  1. Obtenr IC Wald con confianza aproximada del 95% para los cinco parámetros a partir de la función optim.

  2. Obtener IC 95% bootstrap (percentil) para los cinco parámetros, comparando los resultados con los obtenidos en el apartado a).

Ejercicio Propuesto 2.

Obtener el valor crítico y el p-valor bootstrap para resolver los contrastes de hipótesis propuestos en los apartados 1d) y 2g) de la Práctica 3.