Actividad 2 Métodos y Simulación Estadística - Simulación
Johan Berrio y Jorge Agredo
2024-03-18
Problema 1
Estimación de π con simulación
A partir de un cuadrardo de 1x1 y un circulo inscrito en su interior, se requiere estimar el numero π.
Inicialmente se generan n datos aleatorios con distribución uniforme para las variables x y y entre los valores 0 y 1, que conforman los puntos (x,y).
Posteriormente, teniendo en cuenta que el centro de las figuras se encuentra en el punto (0.5 , 0.5) y que si la distancia entre el punto generado y el centro es menor o igual a 0.5 el punto se encuentra al interior del circulo, se calcula la distancia d respecto al centro.
Una vez establecida la distacia se procede a contar los puntos que se encuentran al interior del circulo i.
## [1] 785709Es importante identificar que el cuadrado tiene area de 1 (1x1), mientras que el circulo en el interior tiene area igual a π*r^2, conociendo que r=0.5, entonces el área es igual 0.25π ó π/4, entonces estas figuras tienen una relación de sus áreas 1:π/4. de lo anterior tenemos que:
## [1] 3.142836Para esta simulación la diferencia es de:
## [1] 0.001243346A medida que n crece se acerca más a π =3.1415927.
Tabla con Resultados
| n | π estimado |
|---|---|
| 1000 | 3.128 |
| 10000 | 3.1428 |
| 100000 | 3.13952 |
| 1000000 | 3.142808 |
| 10000000 | 3.1425516 |
Problema 2
Propiedades de los estimadores
La simulación ayuda a entender y validar las propiedades de los estimadores estadísticos como son: insesgadez, eficiencia y la consistencia principalmente. El siguiente problema permite evidenciar las principales características de un grupo de estimadores propuestos para la estimación de un parámetro asociado a un modelo de probabilidad.
estimacion = function(n,lambda){
x1 = rexp(n,lambda)
x2 = rexp(n,lambda)
x3 = rexp(n,lambda)
x4 = rexp(n,lambda)
df = data.frame(x1,x2,x3,x4)
dmin = apply(df,1,min)
dmax = apply(df,1,max)
estim = data.frame(
t1 = (x1+x2)/6 + (x3+x4)/3,
t2 = (x1+2*x2+3*x3+4*x4)/5,
t3 = (x1+x2+x3+x4)/4,
t4 = (dmin+dmax)/2
)
boxplot(estim, main = bquote(paste("Tamaño de muestra: ", .(n))))
abline(h=1/lambda)
print(summary(estim))
}
list_muestra <- list(20, 50, 100, 1000) # Tamaño de la muestra
lambda <- 5
for (n in list_muestra){
estimacion(n, lambda)
}
## t1 t2 t3 t4
## Min. :0.02872 Min. :0.04848 Min. :0.02768 Min. :0.03079
## 1st Qu.:0.17015 1st Qu.:0.31855 1st Qu.:0.16894 1st Qu.:0.18047
## Median :0.18765 Median :0.40165 Median :0.18893 Median :0.21385
## Mean :0.20757 Mean :0.41715 Mean :0.20765 Mean :0.24194
## 3rd Qu.:0.22959 3rd Qu.:0.49438 3rd Qu.:0.23890 3rd Qu.:0.29901
## Max. :0.39599 Max. :0.74194 Max. :0.38194 Max. :0.49467
## t1 t2 t3 t4
## Min. :0.03265 Min. :0.0658 Min. :0.03998 Min. :0.03598
## 1st Qu.:0.12671 1st Qu.:0.2652 1st Qu.:0.12781 1st Qu.:0.16136
## Median :0.17397 Median :0.3443 Median :0.17224 Median :0.20784
## Mean :0.18884 Mean :0.3794 Mean :0.19304 Mean :0.23326
## 3rd Qu.:0.21806 3rd Qu.:0.4627 3rd Qu.:0.23086 3rd Qu.:0.28650
## Max. :0.47315 Max. :0.9568 Max. :0.52325 Max. :0.78327
## t1 t2 t3 t4
## Min. :0.02377 Min. :0.04105 Min. :0.02694 Min. :0.0329
## 1st Qu.:0.11323 1st Qu.:0.21572 1st Qu.:0.11516 1st Qu.:0.1274
## Median :0.17544 Median :0.36337 Median :0.18459 Median :0.2022
## Mean :0.19853 Mean :0.39545 Mean :0.19833 Mean :0.2220
## 3rd Qu.:0.25888 3rd Qu.:0.50290 3rd Qu.:0.24847 3rd Qu.:0.2696
## Max. :0.65746 Max. :1.27960 Max. :0.52667 Max. :0.7638
## t1 t2 t3 t4
## Min. :0.01746 Min. :0.02927 Min. :0.02049 Min. :0.02669
## 1st Qu.:0.12029 1st Qu.:0.23461 1st Qu.:0.12564 1st Qu.:0.14273
## Median :0.18218 Median :0.35753 Median :0.18342 Median :0.20787
## Mean :0.19609 Mean :0.39128 Mean :0.19669 Mean :0.23357
## 3rd Qu.:0.24902 3rd Qu.:0.50218 3rd Qu.:0.24986 3rd Qu.:0.29697
## Max. :0.78741 Max. :1.33214 Max. :0.71222 Max. :0.93206Se observa que, para cada una de las longitudes de muestra evaluadas con un lambda de 5 para la distribución exponencial, el estimador t2 es el que en definitiva menos se acerca a la validación de sus propiedades, estas se explican a continuación:
- Insesgadez: t2 es quien menos se acerca al valor esperado, en este caso a 0.2.
- Consistencia: A medida que se aumenta el tamaño de muestra. t2 tiende a disminuir muy poco el error respecto al valor esperado si se compara con los demás estimadores.
- Eficiencia: t2 es quien presenta mayor dispersión de los datos y por ende mayor varianza
Teniendo en cuenta el análisis anterior, los demás estimadores cumplen mejor las propiedades y puntualmente t3 es quien más se acerca al estimador ideal.
Problema 3
Teorema del Límite Central
Para abordar la verificación que hay una aproximación buena apartir de n>30 comenzamos generando una población de 1000 plantas con una proporción del 50% enfermas, es decir p=0.5, a partir de la cual se seleccionan 500 muestras aleatorias de tamaño m=500, a las cuales se les calcula la media y la desviación estandar de proporción muestral estimada phat.
n = 1000
poblacion <- c(rep(0, n/2), rep(1, n/2))
muestra <- function(poblacion, m) {
muestra <- sample(poblacion, size = m, replace = FALSE)
proporcion_muestral <- mean(muestra)
return(proporcion_muestral)
}La media de las proporciones estimadas phat en 500 simulaciones es 0.500076 y la desviación estándar 0.0160366, siendo un estimador con poca varabilidad y una alta asertividad.
hist(resultados,
main = "Histograma de proporciones muestrales",
xlab = "Proporción muestral",
ylab = "Frecuencia",
breaks = 20,
col = "skyblue",
border = "black")
Respecto a la simetria en la anterior gráfica se puede observar la forma de campana caracteristica de la distribucion normal. Se estima el estadistico de Shapiro-Wilks el cual muestra un p valor alto por no existiría suficiente evidencia para rezhazar la hipotesis nula que provienen de una distribuión normal.
##
## Shapiro-Wilk normality test
##
## data: resultados
## W = 0.99661, p-value = 0.3771Diferentes tamaños de muestra
A continuacíon realizamos el mismo ejecicio para diferentes tamaños de muestras m, en particular para 5, 10, 15, 20, 30, 50, 60, 100, 200 y 500.
tamanos_muestra <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)
resultados <- matrix(NA, nrow = 500, ncol = length(tamanos_muestra))
shapiro_test_pvalues <- numeric(length(tamanos_muestra))
for (i in 1:length(tamanos_muestra)) {
resultados[, i] <- replicate(500, muestra(poblacion, tamanos_muestra[i]))
shapiro_test_pvalues[i] <- shapiro.test(resultados[, i])$p.value
}
medias <- apply(resultados, 2, mean)
desviaciones <- apply(resultados, 2, sd)
diferencia <- medias - 0.5
resultados_df <- data.frame(Tamaño_de_muestra = tamanos_muestra,
Media = round(medias, 3),
Desviación_estándar = round(desviaciones, 3),
Diferencia_con_0.5 = round(diferencia, 3),
Shapiro_Wilk_P_Value = round(shapiro_test_pvalues,3)
)Del ejercicio anterior se obtienen los siguientes resultados:
Mostrando que el p valor de la prueba de Shapiro Wilk comienza a aumentar despues del tamaño de muestra 30, al igual que la desviación estandar comienza a ser cada vez más pequeña. Vale la pena mencionar, que la diferencia con la p de la población es pequeña en todos los niveles, por lo que podria decirse que las 500 simulación vuelven el parametro estimado insesgado.
Graficamente puede observar como los puntos se ubican cerca a la diagonal siendo cada más parecida en la medida que el tamaño de muestra aumenta.
par(mfrow=c(2,3))
for (i in 1:length(tamanos_muestra)) {
qqnorm(resultados[,i], main = paste("QQ Plot - n =", tamanos_muestra[i]))
}
Diferentes proporciones
Realizamos el mismo ejercicio cambiando la proporción de la población a 10% de plantas enfermas.
n = 1000
p = 0.1
poblacion <- c(rep(0, n*(1-p)), rep(1,n*p))
muestra <- function(poblacion, m) {
muestra <- sample(poblacion, size = m, replace = FALSE)
proporcion_muestral <- mean(muestra)
return(proporcion_muestral)
}
tamanos_muestra <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)
resultados <- matrix(NA, nrow = 500, ncol = length(tamanos_muestra))
shapiro_test_pvalues <- numeric(length(tamanos_muestra))
for (i in 1:length(tamanos_muestra)) {
resultados[, i] <- replicate(500, muestra(poblacion, tamanos_muestra[i]))
shapiro_test_pvalues[i] <- shapiro.test(resultados[, i])$p.value
}
medias <- apply(resultados, 2, mean)
desviaciones <- apply(resultados, 2, sd)
diferencia <- medias - p
resultados_df <- data.frame(Tamaño_de_muestra = tamanos_muestra,
Media = round(medias, 3),
Desviación_estándar = round(desviaciones, 3),
Diferencia_con_p = round(diferencia, 3),
Shapiro_Wilk_P_Value = round(shapiro_test_pvalues,3)
)
resultados_dfAl igual que el ejercicio inicial, la desviación estandar del parametro es más pequeña a medida que aumenta el m, por otro lado la diferencia del parametro estimado contra el p de la población sólo no es significativa con m= 5. El estádistico Shapiro Wilks para este ejercicio muestra p_valores cercanos a cero, rezachando la hipotesis nula que provienen de una muestra con distribución normal.
Graficamente, se puede evidenciar el menor ajuste a la diagonal en el gráfico comparado con el ejercicio inicial.
par(mfrow=c(2,3))
for (i in 1:length(tamanos_muestra)) {
qqnorm(resultados[,i], main = paste("QQ Plot - n =", tamanos_muestra[i]))
}
Por ultimo, realizamos el ejercicio con la proporción de la población a 90% de plantas enfermas.
n = 1000
p = 0.9
poblacion <- c(rep(0, n*(1-p)), rep(1,n*p))
muestra <- function(poblacion, m) {
muestra <- sample(poblacion, size = m, replace = FALSE)
proporcion_muestral <- mean(muestra)
return(proporcion_muestral)
}
tamanos_muestra <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)
resultados <- matrix(NA, nrow = 500, ncol = length(tamanos_muestra))
shapiro_test_pvalues <- numeric(length(tamanos_muestra))
for (i in 1:length(tamanos_muestra)) {
resultados[, i] <- replicate(500, muestra(poblacion, tamanos_muestra[i]))
shapiro_test_pvalues[i] <- shapiro.test(resultados[, i])$p.value
}
medias <- apply(resultados, 2, mean)
desviaciones <- apply(resultados, 2, sd)
diferencia <- medias - p
resultados_df <- data.frame(Tamaño_de_muestra = tamanos_muestra,
Media = round(medias, 3),
Desviación_estándar = round(desviaciones, 3),
Diferencia_con_p = round(diferencia, 3),
Shapiro_Wilk_P_Value = round(shapiro_test_pvalues,3)
)
resultados_dfAl igual que el ejercicio inicial, la desviación estandar del parametro es más pequeña a medida que aumenta el m, por otro lado la diferencia contra lel parametro estimado con la p de la población no es significativa con ningun tamaño de muestra. El estádistico Shapiro Wilks para este ejercicio muestra p_valores cercanos a cero, rezachando la hipotesis nula que provienen de una muestra con distribución normal excepto por las simulaciones con m= 500 donde el p_valor es inferior a 0.05.
Graficamente, se puede evidenciar el menor ajuste a la diagonal en el gráfico comparado con el ejercicio inicial.
par(mfrow=c(2,3))
for (i in 1:length(tamanos_muestra)) {
qqnorm(resultados[,i], main = paste("QQ Plot - n =", tamanos_muestra[i]))
}
Problema 4
Estimación boostrap
Cuando se extrae una muestra de una población que no es normal y se requiere estimar un intervalo de confianza se pueden utilizar los métodos de estimación bootstrap. Esta metodología supone que se puede reconstruir la población objeto de estudio mediante un muestreo con reemplazo de la muestra que se tiene. Existen varias versiones del método. Una presentación básica del método se describe a continuación:
El artículo de In-use Emissions from Heavy Duty Dissel Vehicles (J.Yanowitz, 2001) presenta las mediciones de eficiencia de combustible en millas/galón de una muestra de siete camiones. Los datos obtenidos son los siguientes: 7.69, 4.97, 4.56, 6.49, 4.34, 6.24 y 4.45. Se supone que es una muestra aleatoria de camiones y que se desea construir un intervalo de confianza del 95 % para la media de la eficiencia de combustible de esta población. No se tiene información de la distribución de los datos. El método bootstrap permite construir intervalos de confianza del 95 %.
x = c( 7.69, 4.97, 4.56, 6.49, 4.34, 6.24, 4.45) # datos muestra
boot = sample(x,7000,replace=TRUE) # se extraen n x m muestras
b = matrix(boot,nrow=1000,ncol=7) # se construye matriz de n x m
mx = apply(b,1,mean) # se calculan m medias por fila
ic1=quantile(mx, probs=c(0.025, 0.975)) # se calcula IC método 1
cat("Intervalo de confianza por el metodo 1: \n")## Intervalo de confianza por el metodo 1:## 2.5% 97.5%
## 4.705429 6.385714ic2=c(2*mean(mx)-ic1[2], 2*mean(mx)-ic1[1]) # se calcula IC método 2
cat("\n Intervalo de confianza por el metodo 2: \n")##
## Intervalo de confianza por el metodo 2:## 97.5% 2.5%
## 4.673917 6.354203hist(mx, las=1, main="Intervalos de confianza con boostrap", ylab = " ", xlab = " ", col="#034A94")
abline(v=ic1, col="#FF7F00",lwd=2)
abline(v=ic2, col="#0EB0C6",lwd=2)
A partir del remuestreo de la población de mediciones de eficiencia de combustible en los 7 camiones, el método boostrap permite estimar intervalos de confianza más precisos al usar más datos. Con 95% de confianza, el promedio de rendimiento de la gasolina medida en millas/galón se encuentra según el primero método de estimación no paramétrica entre 4.7 y 6.4 y, según el segundo entre 4.6 y 6.3.
Al indagar sobre el consumo real se observó que el promedio es de 6.5 millas galón, quiere decir que el método 1 se acercó más a este valor y, por otro lado, que el consumo promedio de los camiones muestreados tiende a estar por debajo del general.