1. Estimación del valor de π

1.1. Problema

La siguiente figura sugiere cómo estimar el valor de π con una simulación. En la figura, un circuito con un área igual a π/4, está inscrito en un cuadrado cuya área es igual a 1. Se elige de forma aleatoria n puntos dentro del cuadrado. La probabilidad de que un punto esté dentro del círculo es igual a la fracción del área del cuadrado que abarca a este, la cual es π/4. Por tanto, se puede estimar el valor de π/4 al contar el número de puntos dentro del círculo, para obtener la estimación de π/4. De este último resultado se encuentra una aproximación para el valor de π.

  • ¿Cuántos de los puntos están dentro del círculo?
  • ¿Cuál es su estimación de π?

1.2. Resultados y Discusiones

Para dar solución a la consigna se desarrolla la función denominada estimador_pi que toma como parámetro la cantidad n de puntos simulados. Esta función genera n coordenadas (x, y) distribuidas uniformemente en el intervalo [0, 1]. Posteriormente, calcula la distancia entre cada coordenada y el centro del círculo. Las coordenadas cuya distancia al centro sea menor o igual a 0.5 se consideran dentro del área del círculo, mientras que aquellas cuya distancia sea mayor a 0.5 se clasifican fuera del círculo. Utilizando la misma función se realiza una representación gráfica de las coordenadas, asignando el color amarillo a las que se encuentran dentro del circulo y el color azul a las que están fuera. Finalmente, se estima el valor de π multiplicando por 4 la cantidad de puntos dentro del círculo y dividiendo por el número total de puntos n. La fórmula imprime una representacion gráfica con el n correspondiente, los puntos dentro del círculo y el valor de π.

estimador_pi <- function(n){
  
  set.seed(111)
  
  # Generación de coordenadas a partir de distribución uniforme
  x <- runif(n,0,1)
  y <- runif(n,0,1)
  
  # Calculo de distancia de los puntos al centro
  distancia <- (x-0.5)^2+(y-0.5)^2
  
  # Conteo de puntos que estan dentro y fuera del circulo
  puntos_dentro <- which(distancia<0.25)
  puntos_fuera <- which(distancia>=0.25)
  
  # Gráfica de puntos
  plot(x , y, pch=19, las=1,xlim = c(0,1), ylim = c(0,1),asp = 1,cex = 0, xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  points(x[puntos_dentro],y[puntos_dentro], pch =19, col = "#FFEE99", cex = 1)
  points(x[puntos_fuera],y[puntos_fuera], pch =19, col = "#4CC3FF", cex = 1)
  draw.circle(0.5, 0.5, 0.5 , border = "#FFC34C", lw = 5)
  axis(side = 1, at = seq(0, 1, by = 1), labels = c("0", "1"))
  axis(side = 2, at = seq(0, 1, by = 1), labels = c("0", "1"), las = 1)
  
  # Estimación de valor de pi
  valor_pi <- 4*(length(puntos_dentro)/n)
  
  # Especificaciones de prueba
  mtext(paste("n = ", format(n, scientific = FALSE)), side = 1, line = 2, cex = 0.9)
  mtext(paste("Puntos dentro = ", format(length(puntos_dentro), scientific = FALSE)), side = 1, line = 3, cex = 0.9)
  mtext(paste("Valor de π : ", round(valor_pi,4)), side = 1, line = 4, cex = 0.9)
}

Para verificar el nivel de precisión de la estimación de π, se presentan a continuación tres ensayos utilizando la función estimación_pi, donde la variable n toma valores de 1000, 10000 y 100000, respectivamente.

# Pruebas con n valores
par(mfrow = c(1, 3), pty = "s")
par(mar = c(1, 2, 1, 1))
estimador_pi(1000)
estimador_pi(10000)
estimador_pi(100000)

Despues de varios ciclos de simulaciones pueden observarse algunos aspectos relevantes. Por ejemplo, dado que el valor real del area del círculo (π/4) es 0.785, en todas las simulaciones la proporción de puntos que ocupan el interior del círculo se aproxima al 78%. La precisión del modelo para el calculo de π aumenta a medida que aumenta el valor de n, como lo establecen las propiedades de la estadística, entre mayor es la muestra más preciso es el estimador. Así, con un n =1000, el error está al rededor del 0.090, mientras que con la simulación de n = 100.000 está al rededor de 0.002.

2. Propiedades de los estimadores

2.1. Problema

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.

Sean X1, X2, X3 y X4, una muestra aleatoria de tamaño n=4 cuya población la conforma una distribución exponencial con parámetro θ desconocido. Determine las características de cada uno de los siguientes estimadores propuestos:

\[ \hat{\theta1} = \frac{X1 + X2}6 + \frac{X3 + X4}6\] \[ \hat{\theta2} = \frac{X 1 + 2 X2 + 3 X3 + 4X4} 5\] \[ \hat{\theta3} = \frac{X 1 + X2 + X3 + X4} 4\] \[ \hat{\theta4} = \frac{min(X 1 + X2 + X3 + X4)+max(X 1 + X2 + X3 + X4)} 2\]

  • Genere una muestras de n=20, 50, 100 y 1000 para cada uno de los estimadores planteados.

  • En cada caso evalúe las propiedades de insesgadez, eficiencia y consistencia.

2.2. Resultados y Discusiones

a) Para dar solución al problema se implementó una función llamada estimadoresθ. Primero, la función genera las 4 muestras (X1, X2, X3, X4) distintas utilizando la función rexp() para cumplir la condición de que vengan de una población con distribución exponencial. En segundo lugar, se calculan los 4 estimadores (θ1, θ2, θ3, θ4) y se organizan en un Data Frame llamado Mθ. A continuación, se utilizan los datos dentro de Mθ para generar 4 diagramas de cajas y bigotes, que brindan la posibilidad de visualizar el comportamiento de cada estimador respecto a su insesgadez, eficiencia y consistencia.

estimadoresθ <- function(n,lambda){
  
  set.seed(221)
  
  # Generación de números aleatorios
  X1 <- rexp(n,lambda)
  X2 <- rexp(n,lambda)
  X3 <- rexp(n,lambda)
  X4 <- rexp(n,lambda)
  MX = data.frame(X1,X2,X3,X4)

  # Formulas de estimadores
  θ1 <- ((X1+X2)/6) + ((X3+X4)/3)
  θ2 <- (X1+2*X2+3*X3+4*X4)/5
  θ3 <- (X1+X2+X3+X4)/4
  θ4 <- (apply(MX,1,min)+apply(MX,1,max))/2
  
  # Data frame de valores de los estimadores
  Mθ <- data.frame(θ1,θ2,θ3,θ4)
  
  # Gráfica de comparación de estimadores para sesgo y eficiencia
  boxplot(Mθ, ylim = c(0,7), las=1 , col = "#FFEE99",main = paste("n = ", n)) 
  abline(h  = lambda,  col="red") 
  
  # Transformación de data frame para verificación consistencia
  Mθ_trans <- gather(Mθ, key = "θ", value = "Estimacion")
  Mθ_trans$n <- n
  Mθ_trans <- Mθ_trans %>% 
    select(n,everything())
  
  # Devolver Mθ_trans
  return(Mθ_trans)
}

b) En la siguiente cuadricula se disponen de una ronda de simulaciones generadas para los tamaños de muestra de 20, 50, 100 y 1000. En todos los casos se reconoce fácilmente que el estimador más disperso y más sesgado es θ2. También, los estimadores θ1 y θ3 presentan comportamientos similares entre sí en cuanto a insesgadez y eficiencia, aparentemente los mejores del grupo. Por último, para θ4 aunque no se cuenta con una dispersión tan grande como la de θ2, su rango intercuartílico es mayor en comparación con los estimadores de mejor rendimiento.

par(mfrow = c(2, 2), pty = "s")
par(mar = c(2, 1, 2, 1))
Mθ20   <- estimadoresθ(20,1)
Mθ50   <- estimadoresθ(50,1)
Mθ100  <- estimadoresθ(100,1)
Mθ1000 <- estimadoresθ(1000,1)

c) Con el fin de visualizar fácilmente la insesgadez y eficiencia de cada estimador se presenta la siguiente tabla, donde se denota según el tamaño de muestra, la media y la varianza de cada uno. La media, como medida de tendencia central, muestra la insesgadez de los estimadores en tanto se aproxime al valor supuesto de θ, en este caso 1. La Varianza, como medida de dispersión, muestra la eficiencia del estimador en cuanto al considerar los estimadores insesgadoz, se resalta aquel con una dispersión menor.

Teniendo en cuenta lo anterior y los valores presentados en la siguiente tabla, puede afirmarse que los estimadores con menor sesgo son θ1,θ3 y θ4 pues en todas las muestras sus valores son los que más se aproximan a 1. Sin embargo, de los 3 es claro que el estimador θ1 es el que presenta menor sesgo pues en la muestra de 100 presenta una media de 1.064 y en la de 1000 es de 1.004, en contraste, para θ3 en la muestra de 100 se presenta una media de 1.086 y en la de 1000 de 1.012. Por otro lado, para el indicador de eficiencia se examinaron los estimadores (θ1,θ3,θ4) en tanto son los más insesgados. La tabla permite evaluar que θ3 es el más eficiente ya que posee una menor varianza (0.247) en los cuatro tamaños de muestra.

Como se evidenció en el anterior inciso θ2 es el estimador con menor desempeño, siendo el más sesgado ya que la media en todas las iteraciones fue superior a 2, y tambien fue el más disperso manteniendo su varianza por encima de 1.1.

Mθ <- rbind(Mθ20,Mθ50,Mθ100,Mθ1000)
tabla <- Mθ %>% group_by(n,θ) %>% 
  summarise(Media = round(mean(Estimacion),3), Varianza = round(var(Estimacion),3))
datatable(tabla, list(pageLength=16))

d) Por medio del siguiente grafico, es evidente la propiedad de consitencia para los distintos estimadores, en donde a mayor tamaño de muestra el valor estimado tiende al parametro poblacional. Este comportamiento es mas claro para los estimadores θ1,θ3 y θ4 ya que presentan una mejor insesgadez y eficiencia para los tamaños de muestra de 1000, como se discutió en el punto anterior. Por otro lado, también se puede observar un aumento progresivo en la cantidad de valores atípicos en los diagramas de cajas y bigotes a medida que aumenta el tamaño de la muestra. Esto puede deberse a la variabilidad natural de la población o a la presencia de outliers en la muestra.

ggplot(Mθ,aes(x = as.factor(n),y = Estimacion, color = θ)) +
  geom_boxplot() +
  geom_hline(yintercept = 1, color = "red", size = 1) +
  labs(title = "Comprobacion de consistencia de esimadores θ",
       x = "Tamaño de muestra", y = "Valor estimador") +
  theme(axis.text = element_text(color = "black")) +
  theme_minimal()   # Estilo minimalista

3. Teorema del Límite Central

3.1. Problema

El Teorema del Límite Central es uno de los más importantes en la inferencia estadística y habla sobre la convergencia de los estimadores como la proporción muestral a la distribución normal. Algunos autores afirman que esta aproximación es bastante buena a partir del umbral n/> 30. A continuación se describen los siguientes pasos para su verificación:

a) Realice una simulación en la cual genere una población de n=1000 (Lote), donde el porcentaje de individuos (supongamos plantas) enfermas sea del 50%.

b) Genere una función que permita: Obtener una muestra aleatoria de la población y Calcule el estimador de la proporción muestral p para un tamaño de muestra dado n.

c) Repita el escenario anterior (b) n=500 veces y analice los resultados en cuanto al comportamiento de los 500 resultados del estimador p. ¿Qué tan simétricos o sesgados son los resultados obtenidos? y ¿qué se puede observar en cuanto a la variabilidad?. Realice en su informe un comentario sobre los resultados obtenidos.

d) Repita los puntos b y c para tamaños de muestra n=5, 10, 15, 20, 30, 50, 60, 100, 200, 500. Compare los resultados obtenidos para los diferentes tamaños de muestra en cuanto a la normalidad. Utilice pruebas de bondad y ajuste (shapiro wilks :shspiro.test()) y métodos gráficos (gráfico de normalidad: qqnorm()). Comente en su informe los resultados obtenidos

e) Repita toda la simulación (puntos a – d), pero ahora para lotes con 10% de plantas enfermas y de nuevo para lotes con un 90% de plantas enfermas. Concluya sobre los resultados del ejercicio.

3.2. Resultados y Discusiones

a) En principio, como lo pide la consigna, se genera una población de 1000 plantas, 500 de ellas enfermas y 500 sanas. Se emplea la función “sample” para organizar la simulación de manera aleatoria.

# Simulación con lote de n=1000 y 500 enfermas
lote <- c(rep("enferma",500),rep("sana",500))
lote <-  sample(lote)

b) Después, se codifica la función estimador que tome las muestras aleatorias y calcule el estimador, considerando los parámetors n -tamaño de la muestra-, m -número de muestras tomadas- y lote -población principal simulada-. Para ello, la función empieza generando un vector vacío donde se almacenará cada resultado de la muestra. Las muestras se generan mediante un ciclo for que toma del lote el número de muestras indicadas en el parámetro m y las almacena, después calcula la proporción de plantas enfermas en la muestra y almacena todo en el vector vacío “resultados”. Con los datos generados, la función dispone todo en un data frame llamado df_resultados que permite visualizar el tamaño de la muestra, el número de repeticiones a la que corresponde cada una y la proporción de plantas enfermas en ella. Finalmente, la función calcula el estimador para el conjunto de muestras generadas, en este caso, el promedio de la proporción de plantas enfermas para el conjunto de muestras.

# Función muestra aleatoria y calculo del estimador de la proporción muestral 
estimador <- function(n, m,lote) {
  
  set.seed(111)
  
    # Crear un vector vacío para almacenar los resultados
  resultados <- vector("list", m)
  
  # Repetir la evaluación m veces
  for (i in 1:m) {
    # Generar la muestra
    muestra <- sample(lote, n)
    # Calcular la proporción de enfermas
    p_enfermas <- sum(muestra == "enferma") / n
    # Almacenar los resultados de esta repetición en la lista
    resultados[[i]] <- c(n = n, Repeticion = i, p_enfermas = p_enfermas)
  }
  
  # Convertir la lista de resultados en un data frame
  df_resultados <- as.data.frame(do.call(rbind, resultados))
  
  #Calculo del estimador P
  calc_estimador <- round(mean(df_resultados$p_enfermas),3)
  
  return(df_resultados)
}

c) La siguiente cuadricula muestra un primer escenario con un tamaño de muestra 50 y 500 repeticiones de la muestra tomada. Al lado izquierdo, un histograma nos permite ver el estimador calculado, para el caso, la proporción de plantas enfermas y cómo se distribuye su cálculo en las 500 muestras tomadas. El gráfico presenta una distribución aproximadamente simétrica, la mayor proporción de plantas enfermas se acumula entre el 45% y 60%. En el lado derecho se presenta un gráfico de normalidad que permite visualizar el comportamiento real del estimador frente a una distribución normal. Para el gráfico se representa una distribución no normal porque los puntos no están alineados con la recta diagonal que representa la coincidencia entre percentiles teóricos normales con los percentiles muestrales.

# Escenario con n = 50 y m = 500

df_resultados <- estimador(50,500,lote)
par(mfrow = c(1, 2), pty = "s", cex=1, cex.axis=0.5, cex.lab=0.5, cex.main=1)
par(mar = c(2, 4, 1, 1))
hist(df_resultados$p_enfermas, freq = FALSE, main = "n = 50", xlab = "p_enfermas")
qqnorm(df_resultados$p_enfermas, main= "n = 50") ; qqline(df_resultados$p_enfermas, col="red")

d) A continuación, se realizan simulaciones con muestras de n = 5, 10, 15, 20, 30, 50, 60, 100, 200, 500, para visualizar su distribución se ejecuta un histograma para cada caso. Según el teorema del límite central, entre mayor el tamaño de la muestra la distribución del estimador se aproximaría cada vez más a una configuración normal. En los 10 escenarios analizados podemos apreciar que esta tendencia se cumple, particularmente desde la muestra n = 60 siendo más evidente en la muestra de n = 500. Lo anterior se constata gracias al p-value arrojado por el Shapiro test que para n = 500 es de 0.4974, valor que supera el nivel de significancia de 0.05 por lo que puede clasificarse como distribución normal.

# Escenario con n = c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500) y m = 500
e_5 <- estimador(5,500,lote)
e_10 <- estimador(10,500,lote)
e_15 <- estimador(15,500,lote)
e_20 <- estimador(20,500,lote)
e_30 <- estimador(30,500,lote)
e_50 <- estimador(50,500,lote)
e_60 <- estimador(60,500,lote)
e_100 <- estimador(100,500,lote)
e_200 <- estimador(200,500,lote)
e_500 <- estimador(500,500,lote)

par(mfrow = c(2, 5), mar = c(2, 4, 3, 1))
hist(e_5$p_enfermas, freq = FALSE, main = paste("n = 5 |", "p-value = ",format(shapiro.test(e_5$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")
hist(e_10$p_enfermas, freq = FALSE, main = paste("n = 10 |", "p-value = ",format(shapiro.test(e_10$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")
hist(e_15$p_enfermas, freq = FALSE, main = paste("n = 15 |", "p-value = ",format(shapiro.test(e_15$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")
hist(e_20$p_enfermas, freq = FALSE, main = paste("n = 20 |", "p-value = ",format(shapiro.test(e_20$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")
hist(e_30$p_enfermas, freq = FALSE, main = paste("n = 30 |", "p-value = ",format(shapiro.test(e_30$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")
hist(e_50$p_enfermas, freq = FALSE, main = paste("n = 50 |", "p-value = ",format(shapiro.test(e_50$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")
hist(e_60$p_enfermas, freq = FALSE, main = paste("n = 60 |", "p-value = ",format(shapiro.test(e_60$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")
hist(e_100$p_enfermas, freq = FALSE, main = paste("n = 100 |", "p-value = ",format(shapiro.test(e_100$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")
hist(e_200$p_enfermas, freq = FALSE, main = paste("n = 200 |", "p-value = ",format(shapiro.test(e_200$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")
hist(e_500$p_enfermas, freq = FALSE, main = paste("n = 500 |", "p-value = ",format(shapiro.test(e_500$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")

Los resultados muestran desde un gráfico de valores no normales hasta un gráfico de puntos completamente alineados con la recta diagonal que representa la coincidencia entre los percentiles teóricos normales con los percentiles muestrales, indicando así una convergencia a la distribución normal.

par(mfrow = c(2, 5), mar = c(5, 4, 2, 1))
qqnorm(e_5$p_enfermas, main= "n = 5") ; qqline(e_5$p_enfermas, col="red")
qqnorm(e_10$p_enfermas, main= "n = 10") ; qqline(e_10$p_enfermas, col="red")
qqnorm(e_15$p_enfermas, main= "n = 15") ; qqline(e_15$p_enfermas, col="red")
qqnorm(e_20$p_enfermas, main= "n = 20") ; qqline(e_20$p_enfermas, col="red")
qqnorm(e_30$p_enfermas, main= "n = 30") ; qqline(e_30$p_enfermas, col="red")
qqnorm(e_50$p_enfermas, main= "n = 50") ; qqline(e_50$p_enfermas, col="red")
qqnorm(e_60$p_enfermas, main= "n = 60") ; qqline(e_60$p_enfermas, col="red")
qqnorm(e_100$p_enfermas, main= "n = 100") ; qqline(e_100$p_enfermas, col="red")
qqnorm(e_200$p_enfermas, main= "n = 200") ; qqline(e_200$p_enfermas, col="red")
qqnorm(e_500$p_enfermas, main= "n = 500") ; qqline(e_500$p_enfermas, col="red")

Finalmente, se presenta una visualización unificada del comportamiento del estimador utilizando gráficos de cajas y bigotes, donde es posible validar que a mayor tamaño de muestra este se aproxima cada vez más al parámetro poblacional siendo su dispersión considerablemente menor, derivado de la reducción del rango intercuartílico.

resultados <- rbind(e_5,e_10,e_15,e_20,e_30,e_50,e_60,e_100,e_200,e_500)

ggplot(resultados, aes(x = as.factor(n), y = p_enfermas)) +
  geom_boxplot(fill = "#FFEE99") +
  geom_hline(yintercept = 0.5, color = "red", lwd = 1) +
  labs(x = "Tamano de Muestra", y = "Proporcion de enfermas",
       title = "Grafica de proporcion de plantas enfermas por tamano de Muestra") +
  theme(axis.text = element_text(color = "black")) +
  theme_minimal()   # Estilo minimalista

e) A partir de esta sección se realiza una simulación de una población diferente, con una proporción de 10% de plantas enfermas, es decir, 900 plantas sanas y 100 enfermas, aplicando el mismo procedimiento descrito desde los numerales a hasta el d.

# Escenario con proporcion de enfermas de 10%
lote <- c(rep("enferma",100),rep("sana",900))
lote <-  sample(lote)

# Escenario con n = 50 y m = 500
df_resultados <- estimador(50,500,lote)
par(mfrow = c(1, 2), pty = "s", cex=1, cex.axis=0.5, cex.lab=0.5, cex.main=1)
par(mar = c(2, 4, 1, 1))
hist(df_resultados$p_enfermas, freq = FALSE, main = "n = 50", xlab = "p_enfermas")
qqnorm(df_resultados$p_enfermas, main= "n = 50") ; qqline(df_resultados$p_enfermas, col="red")

# Escenario con n = c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500) y m = 500
e_5 <- estimador(5,500,lote)
e_10 <- estimador(10,500,lote)
e_15 <- estimador(15,500,lote)
e_20 <- estimador(20,500,lote)
e_30 <- estimador(30,500,lote)
e_50 <- estimador(50,500,lote)
e_60 <- estimador(60,500,lote)
e_100 <- estimador(100,500,lote)
e_200 <- estimador(200,500,lote)
e_500 <- estimador(500,500,lote)

par(mfrow = c(2, 5), mar = c(2, 4, 3, 1))
hist(e_5$p_enfermas, freq = FALSE, main = paste("n = 5 |", "p-value = ",format(shapiro.test(e_5$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")
hist(e_10$p_enfermas, freq = FALSE, main = paste("n = 10 |", "p-value = ",format(shapiro.test(e_10$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")
hist(e_15$p_enfermas, freq = FALSE, main = paste("n = 15 |", "p-value = ",format(shapiro.test(e_15$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")
hist(e_20$p_enfermas, freq = FALSE, main = paste("n = 20 |", "p-value = ",format(shapiro.test(e_20$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")
hist(e_30$p_enfermas, freq = FALSE, main = paste("n = 30 |", "p-value = ",format(shapiro.test(e_30$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")
hist(e_50$p_enfermas, freq = FALSE, main = paste("n = 50 |", "p-value = ",format(shapiro.test(e_50$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")
hist(e_60$p_enfermas, freq = FALSE, main = paste("n = 60 |", "p-value = ",format(shapiro.test(e_60$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")
hist(e_100$p_enfermas, freq = FALSE, main = paste("n = 100 |", "p-value = ",format(shapiro.test(e_100$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")
hist(e_200$p_enfermas, freq = FALSE, main = paste("n = 200 |", "p-value = ",format(shapiro.test(e_200$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")
hist(e_500$p_enfermas, freq = FALSE, main = paste("n = 500 |", "p-value = ",format(shapiro.test(e_500$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")

par(mfrow = c(2, 5), mar = c(5, 4, 2, 1))
qqnorm(e_5$p_enfermas, main= "n = 5") ; qqline(e_5$p_enfermas, col="red")
qqnorm(e_10$p_enfermas, main= "n = 10") ; qqline(e_10$p_enfermas, col="red")
qqnorm(e_15$p_enfermas, main= "n = 15") ; qqline(e_15$p_enfermas, col="red")
qqnorm(e_20$p_enfermas, main= "n = 20") ; qqline(e_20$p_enfermas, col="red")
qqnorm(e_30$p_enfermas, main= "n = 30") ; qqline(e_30$p_enfermas, col="red")
qqnorm(e_50$p_enfermas, main= "n = 50") ; qqline(e_50$p_enfermas, col="red")
qqnorm(e_60$p_enfermas, main= "n = 60") ; qqline(e_60$p_enfermas, col="red")
qqnorm(e_100$p_enfermas, main= "n = 100") ; qqline(e_100$p_enfermas, col="red")
qqnorm(e_200$p_enfermas, main= "n = 200") ; qqline(e_200$p_enfermas, col="red")
qqnorm(e_500$p_enfermas, main= "n = 500") ; qqline(e_500$p_enfermas, col="red")

resultados <- rbind(e_5,e_10,e_15,e_20,e_30,e_50,e_60,e_100,e_200,e_500)

ggplot(resultados, aes(x = as.factor(n), y = p_enfermas)) +
  geom_boxplot(fill = "#FFEE99") +
  geom_hline(yintercept = 0.1, color = "red", lwd = 1) +
  labs(x = "Tamano de Muestra", y = "Proporcion de enfermas",
       title = "Grafica de proporcion de plantas enfermas por tamano de Muestra") +
  theme(axis.text = element_text(color = "black")) +
  theme_minimal()   # Estilo minimalista

A partir de esta sección se realiza una simulación de una población diferente, con una proporción de 90% de plantas enfermas, es decir, 100 plantas sanas y 900 enfermas, aplicando el mismo procedimiento descrito desde los numerales a hasta el d.

# Escenario con proporcion de enfermas de 90%
lote <- c(rep("enferma",900),rep("sana",100))
lote <-  sample(lote)

# Escenario con n = 100 y m = 500
df_resultados <- estimador(50,500,lote)
par(mfrow = c(1, 2), pty = "s", cex=1, cex.axis=0.5, cex.lab=0.5, cex.main=1)
par(mar = c(2, 4, 1, 1))
hist(df_resultados$p_enfermas, freq = FALSE, main = "n = 50", xlab = "p_enfermas")
qqnorm(df_resultados$p_enfermas, main= "n = 50") ; qqline(df_resultados$p_enfermas, col="red")

# Escenario con n = c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500) y m = 500
e_5 <- estimador(5,500,lote)
e_10 <- estimador(10,500,lote)
e_15 <- estimador(15,500,lote)
e_20 <- estimador(20,500,lote)
e_30 <- estimador(30,500,lote)
e_50 <- estimador(50,500,lote)
e_60 <- estimador(60,500,lote)
e_100 <- estimador(100,500,lote)
e_200 <- estimador(200,500,lote)
e_500 <- estimador(500,500,lote)

par(mfrow = c(2, 5), mar = c(2, 4, 3, 1))
hist(e_5$p_enfermas, freq = FALSE, main = paste("n = 5 |", "p-value = ",format(shapiro.test(e_5$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")
hist(e_10$p_enfermas, freq = FALSE, main = paste("n = 10 |", "p-value = ",format(shapiro.test(e_10$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")
hist(e_15$p_enfermas, freq = FALSE, main = paste("n = 15 |", "p-value = ",format(shapiro.test(e_15$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")
hist(e_20$p_enfermas, freq = FALSE, main = paste("n = 20 |", "p-value = ",format(shapiro.test(e_20$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")
hist(e_30$p_enfermas, freq = FALSE, main = paste("n = 30 |", "p-value = ",format(shapiro.test(e_30$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")
hist(e_50$p_enfermas, freq = FALSE, main = paste("n = 50 |", "p-value = ",format(shapiro.test(e_50$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")
hist(e_60$p_enfermas, freq = FALSE, main = paste("n = 60 |", "p-value = ",format(shapiro.test(e_60$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")
hist(e_100$p_enfermas, freq = FALSE, main = paste("n = 100 |", "p-value = ",format(shapiro.test(e_100$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")
hist(e_200$p_enfermas, freq = FALSE, main = paste("n = 200 |", "p-value = ",format(shapiro.test(e_200$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")
hist(e_500$p_enfermas, freq = FALSE, main = paste("n = 500 |", "p-value = ",format(shapiro.test(e_500$p_enfermas)$p.value,scientific = TRUE,digits = 4)), xlab = "p_enfermas")

par(mfrow = c(2, 5), mar = c(5, 4, 2, 1))
qqnorm(e_5$p_enfermas, main= "n = 5") ; qqline(e_5$p_enfermas, col="red")
qqnorm(e_10$p_enfermas, main= "n = 10") ; qqline(e_10$p_enfermas, col="red")
qqnorm(e_15$p_enfermas, main= "n = 15") ; qqline(e_15$p_enfermas, col="red")
qqnorm(e_20$p_enfermas, main= "n = 20") ; qqline(e_20$p_enfermas, col="red")
qqnorm(e_30$p_enfermas, main= "n = 30") ; qqline(e_30$p_enfermas, col="red")
qqnorm(e_50$p_enfermas, main= "n = 50") ; qqline(e_50$p_enfermas, col="red")
qqnorm(e_60$p_enfermas, main= "n = 60") ; qqline(e_60$p_enfermas, col="red")
qqnorm(e_100$p_enfermas, main= "n = 100") ; qqline(e_100$p_enfermas, col="red")
qqnorm(e_200$p_enfermas, main= "n = 200") ; qqline(e_200$p_enfermas, col="red")
qqnorm(e_500$p_enfermas, main= "n = 500") ; qqline(e_500$p_enfermas, col="red")

resultados <- rbind(e_5,e_10,e_15,e_20,e_30,e_50,e_60,e_100,e_200,e_500)

ggplot(resultados, aes(x = as.factor(n), y = p_enfermas)) +
  geom_boxplot(fill = "#FFEE99") +
  geom_hline(yintercept = 0.9, color = "red", lwd = 1) +
  labs(x = "Tamano de Muestra", y = "Proporcion de enfermas",
       title = "Grafica de proporcion de plantas enfermas por tamano de Muestra") +
  theme(axis.text = element_text(color = "black")) +
  theme_minimal()   # Estilo minimalista

De las simulaciones realizadas en la última sección, con poblaciones distintas a la original se denota que a partir de las muestras de tamaño n = 60 la distribución del estimador tiende a ser simétrica. La normalidad apenas se vislumbra en el valor de n = 500, aunque cabe aclarar que el p-value no es suficiente para aceptar la hipótesis nula de que la distribución es normal. Finalmente, puede notarse mayor acumulación de plantas enfermas por encima del parámetro poblacional definido para los casos del 10% y 90% de plantas enfermas.

4. Estimacción boostrap

4.1. Problema

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 % - Para ilustrar el método suponga que coloca los valores de la muestra en una caja y extrae uno al azar. Este correspondería al primer valor de la muestra bootstrap X1. Después de anotado el valor se regresa X1 a la caja y se extrae el valor X2 , regresandolo nuevamente. Este procedimiento se repite hasta completar una muestra de tamaño n , X1 ,X2 ,X3 ,Xn , conformando la muestra bootstrap.

Es necesario extraer un gran número de muestras (suponga k = 1000). Para cada una de las muestra bootstrap obtenidas se calcula la media Xi , obteniéndose un valor para cada muestra. El intervalo de confianza queda conformado por los percentiles P2.5 y P97.5 . Existen dos métodos para estimarlo:

Método 1: (P2.5;P97.5)
Método 2: (2Xprom−P97.5;2Xprom−P2.5)

Construya el intervalo de confianza por los dos métodos y compare los resultados obtenidos. Comente los resultados. ¿Confiaría en estas estimaciones?

4.2. Resultados y Discusiones

Inicialmente se crea el vector x con el registro de las mediciones de eficiencia de combustible en millas/galón de los siete camiones indicados. Después, se realiza un muestreo aleatorio denominado bootstrap de tamaño n = 5000 con reemplazo. A partir de ello se construye una matriz de 5 columnas con 1000 registros, en donde por cada fila se calcula el promedio de los 5 valores que le corresponden.

x <- c( 7.69, 4.97, 4.56, 6.49, 4.34, 6.24, 4.45) 
bootstrap <- sample(x,5000,replace=TRUE)
b <- matrix(bootstrap,nrow=1000,ncol=5)
mx <- apply(b,1,mean)

Después de crear la matriz, se aplican a ella los métodos sugeridos en la consigna para generar los intervalos de confianza del 95%.

Para el primer método, el límite inferior corresponde al percentil 2.5 y el límite superior al percentil 97.5, como se denota a continuación:

ic1 <- quantile(mx, probs=c(0.025, 0.975)) # se calcula IC método 1
ic1
##  2.5% 97.5% 
## 4.532 6.680

En el segundo método, el límite inferior se calcula como dos veces el promedio del promedio de las filas menos el percentil 97.5, y el límite superior es dos veces el promedio del promedio de las filas menos el percentil 2.5, como se muestra a continuación:

ic2 <- c(2*mean(mx)-ic1[2], 2*mean(mx)-ic1[1]) # se calcula IC método 2
ic2
##    97.5%     2.5% 
## 4.354104 6.502104

Para el primer método se obtiene un intervalo entre 4.532 y 6.680. En contraste, para el segundo método el intervalo de confianza está constituido entre 4.354 y 6.502. Con lo anterior, se evidencia que el segundo método realiza una corrección por empleando la media, haciendo que el intervalo se posicione 0.178 por debajo respecto al primer método. También, se observa que en cualquiera de las dos estimaciones se contiene la media de la muestra inicialmente planteada, y como era de suponerse, los dos métodos excluyen el valor mínimo y máximo de la muestra. Según el análisis, ambos métodos parecen ser confiables, ya que los intervalos de confianza del Método 1 y Método 2 se superponen considerablemente, sugiriendo una consistencia en las estimaciones.

hist(mx, las=1, main=" ", ylab = " ", xlab = " ", col="#034A94")
abline(v=ic1, col="#FFC34C",lwd=2)
abline(v=ic2, col="#4CC3FF",lwd=2)
abline(v=mean(x), col = "red", lwd=2)
legend("topright", legend=c("Método 1", "Método 2", "Media muestral"), col=c("#FFC34C", "#4CC3FF", "red"), lwd=2, lty=1)