PROBLEMA 1: Estimación del valor de \(\pi\)

La siguiente figura sugiere como estimar el valor \(\pi\) de con una simulación. En la figura, un circuito con un área igual a \(\pi\)/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 \(\pi\)/4. Por tanto, se puede estimar el valor de \(\pi\)/4 al contar el número de puntos dentro del círculo, para obtener la estimación de \(\pi\)/4. De este último resultado encontrar una aproximación para el valor de \(\pi\).

Descripción de la imagen

Pasos:

  • Genere n coordenadas x: X1, . . . ,Xn. Utilice la distribución uniforme con valor mínimo de 0 y valor máximo de 1. La distribución uniforme genera variables aleatorias que tienen la misma probabilidad de venir de cualquier parte del intervalo (0,1).

  • Genere 1000 coordenadas y: Y1, . . . ,Yn utilizando nuevamente la distribución uniforme con valor mínimo de 0 y valor máximo de 1.

  • Cada punto (Xi,Yi) se encuentra dentro del círculo si su distancia desde el centro (0.5, 0.5) es menor a 0.5. Para cada par (Xi,Yi) determine si la distancia desde el centro es menor a 0.5. Esto último se puede realizar al calcular el valor \((X_i - 0.5)^2 + (Y_i - 0.5)^2\), que es el cuadrado de la distancia, y al determinar si es menor que 0.25.

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

Solución

El método de Monte Carlo es una técnica estadística que se utiliza para aproximar soluciones numéricas mediante la generación de muestras aleatorias. En este contexto, se realiza una simulación para estimar el valor de π utilizando el método de Monte Carlo.

En la definición de la función que realiza la estimación de π:

  • Generación de coordenadas: En primer lugar, se generan pares de coordenadas aleatorias (x, y) dentro de un cuadrado unitario (de lado 1), utilizando una distribución uniforme (0,1). Estas coordenadas representan puntos distribuidos uniformemente en el cuadrado.

  • Determinación de puntos dentro del círculo: Se calcula si cada par de coordenadas (x, y) generado está dentro de un círculo de radio 0.5 centrado en (0.5, 0.5). Esto se logra calculando la distancia euclidiana desde cada punto al centro del círculo y comprobando si es menor que el radio. Los puntos que cumplen esta condición se consideran que se encuentran “dentro” del círculo.

\((X_i - 0.5)^2 + (Y_i - 0.5)^2\)

  • Estimación de π mediante la proporción de puntos dentro del círculo: Dado que el área del círculo es πr² (donde r es el radio), y el área del cuadrado que contiene al círculo es r², la razón entre el área del círculo y el área del cuadrado es π/4. Por lo tanto, la proporción de puntos generados que caen dentro del círculo se puede utilizar para estimar π multiplicándola por 4.

\[\hat{\pi} = \frac{4 \times \sum_{i=1}^{n} \text{puntos_dentro}}{n}\]

  • Visualización de los resultados: Se visualizan los puntos generados en un gráfico de dispersión utilizando ggplot, donde los puntos dentro del círculo se representan de color “rojo”. Se realiza la simulación para los 3 tamaños de muestra (100, 1.000 y 10.000), mostrando la cantidad de puntos que se encuentran dentro del círculo y el valor estimado de \(\pi\)
set.seed(123)

# Se define la función para estimar pi 

estimacion_pi <- function(n) {
  # Generar coordenadas x e y
  x <- runif(n, 0, 1)
  y <- runif(n, 0, 1)
  
  # Almacenar coordenadas y si están dentro o fuera del círculo
  df <- data.frame(
    x = x,
    y = y,
    puntos_dentro = (x - 0.5)^2 + (y - 0.5)^2 < 0.25 )
  
  # Calcular estimación de pi
  pi_estimado <- 4 * sum(df$puntos_dentro) / n
  
  # Contar la cantidad de puntos dentro del círculo
  conteo_puntos <- sum(df$puntos_dentro)
  
  # Crear el gráfico
  p <- ggplot(df, aes(x, y)) +
    geom_point(aes(color = puntos_dentro)) +
    scale_color_manual(values = c("azure3", "coral2")) +
    labs(x = "X", y = "Y", title = paste("n =", format(n, scientific = FALSE), ", Puntos dentro=", conteo_puntos)) +
    theme_minimal() +
    theme(axis.text = element_text(size = 6), legend.position = "none", aspect.ratio = 1) # Eliminar la leyenda y ajustar el aspecto del gráfico
  
  return(list(p, pi_estimado, conteo_puntos))
}

# Realizar la simulación para diferentes cantidades de puntos
n_simulacion <- c(1000, 10000, 100000) # Reducir el número de puntos

# Almacenar las estimaciones de pi y la cantidad de puntos dentro del círculo
simulacion_pi_estimado <- numeric(length(n_simulacion))
simulacion_conteo_puntos <- numeric(length(n_simulacion))
plots <- list()

# Llamar a la función para estimar pi y graficar los puntos
for (i in 1:length(n_simulacion)) {
  result <- estimacion_pi(n_simulacion[i])
  plots[[i]] <- result[[1]]
  simulacion_pi_estimado[i] <- result[[2]]
  simulacion_conteo_puntos[i] <- result[[3]]
}

# Mostrar los gráficos en un lienzo múltiple
grid.arrange(grobs = plots, ncol = 3) # No es necesario ajustar el ancho del lienzo

# Crear un marco de datos con los pi estimados y la cantidad de puntos dentro de cada círculo
results_df <- data.frame(
  N_Puntos = format(n_simulacion, scientific = FALSE), 
  Pi_Estimado = simulacion_pi_estimado,
  Puntos_Dentro = simulacion_conteo_puntos
)

names(results_df) <- c("n Puntos", "Pi Estimado", "Puntos Dentro")
#print(results_df)

# Generar una tabla con kable
kable(results_df, format = "markdown")
n Puntos Pi Estimado Puntos Dentro
1000 3.2000 800
10000 3.1528 7882
100000 3.1440 78600
A medida que aumentamos el número de puntos generados (tamaño de muestra), la precisión de la estimación de \(\pi\) tiende a mejorar. Esto se debe a que con más puntos, tenemos una mejor representación de la distribución de puntos dentro del cuadrado unitario, lo que resulta en una estimación más precisa de la proporción de puntos dentro del círculo y, por lo tanto, de \(\pi\).

PROBLEMA 2: Propiedades de los estimadores

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{θ}_1 = \frac{{X_1 + X_2}}{6} + \frac{{X_3 + X_4}}{3}\]

\[\hat{θ}_2 = \frac{{X_1 + 2X_2 + 3X_3 + 4X_4}}{5}\]

\[\hat{θ}_3 = \frac{{X_1 + X_2 + X_3 + X_4}}{4}\]

\[\hat{θ}_4 = \frac{{\min\{X_1, X_2, X_3, X_4\} + \max\{X_1, X_2, X_3, X_4\}}}{2}\] Pasos:

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

  • En cada caso evalue las propiedades de insesgadez, eficiencia y consistencia.

  • Suponga un valor para el parámetro θ.

Solución

Se realiza la simulación de datos muestrales que siguen una distribución exponencial, para calcular y comparar el comportamiento de cuatro estimadores diferentes de un parámetro θ,en función del tamaño de la muestra.

  • Parámetros iniciales: Se definen los tamaños de las muestras, donde se evaluaran los estaimadores para 20, 50, 100 y 1.000 obervaciones, y se establece el parámetro θ en 2.

  • Cálculo y simulación de los estimadores: Definidos los parámetros iniciales, se generan cuatro muestras exponenciales independientes, todas con el mismo parámetro θ, esto es iterado sobre los diferentes tamaños de muestra y son calculador los 4 estimadores para cada uno de ellos.

  • Evaluación de propiedades y estimadores: Para cada estimador se evalua el sesgo, la varianza, la eficiencia y la consistencia.Se devuelve un resumen de las propiedades y los valores estimados, para cada tamaño de muestra y estimador. Finalmente, se crea un boxplot para visualizar y comparar la distribución de los valores estimados de cada estimador en los diferentes tamaños de muestra.

set.seed(123)

# Definir tamaños de las muestras
tamanos_muestra <- c(20, 50, 100, 1000)

# Definir parámetro theta conocido
theta <- 2

# Función para calcular los estimadores
calcular_estimadores <- function(muestra, theta) {
  # Generar las muestras para X1, X2, X3 y X4
  x1 <- rexp(muestra, rate = theta)
  x2 <- rexp(muestra, rate = theta)
  x3 <- rexp(muestra, rate = theta)
  x4 <- rexp(muestra, rate = theta)
  df <- data.frame(x1, x2, x3, x4)
  
  # Calcular estimadores
  theta_1 <- (x1 + x2)/6 + (x3 + x4)/3
  theta_2 <- (x1 + 2*x2 + 3*x3 + 4*x4)/5
  theta_3 <- (x1 + x2 + x3 + x4)/4
  theta_4 <- (apply(df, 1, min) + apply(df, 1, max))/2
  
  # Calcular sesgo, eficiencia y consistencia
  sesgo <- c(mean(theta_1) - theta, mean(theta_2) - theta, mean(theta_3) - theta, mean(theta_4) - theta)
  varianza <- c(var(theta_1), var(theta_2), var(theta_3), var(theta_4))
  eficiencia <- c(theta^2 / varianza[1], theta^2 / varianza[2], theta^2 / varianza[3], theta^2 / varianza[4])
  consistencia <- c(mean((theta_1 - theta)^2), mean((theta_2 - theta)^2), mean((theta_3 - theta)^2), mean((theta_4 - theta)^2))
  
  propiedades <- data.frame(Sesgo = sesgo, Varianza = varianza, Eficiencia = eficiencia, Consistencia = consistencia)
  
  # Crear dataframe con los resultados
  estimadores <- data.frame(Tamaño_Muestra = muestra,
                            Estimador = c(rep("Theta_1", length(theta_1)),
                                          rep("Theta_2", length(theta_2)),
                                          rep("Theta_3", length(theta_3)),
                                          rep("Theta_4", length(theta_4))),
                            Valor = c(theta_1, theta_2, theta_3, theta_4))
  
  # Devolver resultados
  return(list(propiedades = propiedades, estimadores = estimadores))
}

# Realizar simulación para cada tamaño de muestra
resultados_list <- lapply(tamanos_muestra, function(n) calcular_estimadores(n, theta))

# Mostrar los resultados
for (i in seq_along(tamanos_muestra)) {
  cat("Resultados para n =", tamanos_muestra[i], "\n")
  print(resultados_list[[i]]$propiedades)
  cat("\n")
}
## Resultados para n = 20 
##        Sesgo   Varianza Eficiencia Consistencia
## 1 -1.4822878 0.06803333   58.79471     2.261809
## 2 -0.9780206 0.21921925   18.24657     1.164783
## 3 -1.4914356 0.04964954   80.56469     2.271547
## 4 -1.3996061 0.13018168   30.72629     2.082570
## 
## Resultados para n = 50 
##        Sesgo   Varianza Eficiencia Consistencia
## 1 -1.4843870 0.05079508   78.74779     2.253184
## 2 -0.9676119 0.21182705   18.88333     1.143863
## 3 -1.4892535 0.04394524   91.02237     2.260942
## 4 -1.4314831 0.06360973   62.88346     2.111482
## 
## Resultados para n = 100 
##        Sesgo   Varianza Eficiencia Consistencia
## 1 -1.5010224 0.06791158   58.90012     2.320301
## 2 -0.9890601 0.27559371   14.51412     1.251078
## 3 -1.5017651 0.06004551   66.61614     2.314743
## 4 -1.4155887 0.08874438   45.07328     2.091748
## 
## Resultados para n = 1000 
##        Sesgo   Varianza Eficiencia Consistencia
## 1 -1.4980332 0.07352848   54.40069     2.317558
## 2 -0.9970921 0.31802507   12.57762     1.311900
## 3 -1.4974735 0.06877106   58.16400     2.311129
## 4 -1.4125413 0.10572378   37.83444     2.100891
library(knitr)

# Combinar todos los dataframes de propiedades en uno solo
propiedades_combinadas <- do.call(rbind, lapply(seq_along(resultados_list), function(i) {
  prop <- resultados_list[[i]]$propiedades
  prop$Estimador <- rep(c("Theta_1", "Theta_2", "Theta_3", "Theta_4"), each = nrow(prop) / 4)
  prop$Tamaño_Muestra <- rep(tamanos_muestra[i], each = nrow(prop) / length(tamanos_muestra))
  return(prop)
}))

# Reorganizar las columnas
propiedades_combinadas <- propiedades_combinadas[, c("Tamaño_Muestra", "Estimador", names(propiedades_combinadas)[1:(ncol(propiedades_combinadas)-2)])]

# Mostrar los resultados en una tabla con kable
print(kable(propiedades_combinadas))
## 
## 
## | Tamaño_Muestra|Estimador |      Sesgo|  Varianza| Eficiencia| Consistencia|
## |--------------:|:---------|----------:|---------:|----------:|------------:|
## |             20|Theta_1   | -1.4822878| 0.0680333|   58.79471|     2.261809|
## |             20|Theta_2   | -0.9780206| 0.2192192|   18.24657|     1.164783|
## |             20|Theta_3   | -1.4914356| 0.0496495|   80.56469|     2.271547|
## |             20|Theta_4   | -1.3996061| 0.1301817|   30.72629|     2.082570|
## |             50|Theta_1   | -1.4843870| 0.0507951|   78.74779|     2.253184|
## |             50|Theta_2   | -0.9676119| 0.2118270|   18.88333|     1.143863|
## |             50|Theta_3   | -1.4892535| 0.0439452|   91.02237|     2.260943|
## |             50|Theta_4   | -1.4314831| 0.0636097|   62.88346|     2.111482|
## |            100|Theta_1   | -1.5010224| 0.0679116|   58.90012|     2.320301|
## |            100|Theta_2   | -0.9890601| 0.2755937|   14.51412|     1.251078|
## |            100|Theta_3   | -1.5017651| 0.0600455|   66.61614|     2.314743|
## |            100|Theta_4   | -1.4155887| 0.0887444|   45.07328|     2.091748|
## |           1000|Theta_1   | -1.4980332| 0.0735285|   54.40069|     2.317558|
## |           1000|Theta_2   | -0.9970921| 0.3180251|   12.57762|     1.311900|
## |           1000|Theta_3   | -1.4974735| 0.0687711|   58.16400|     2.311129|
## |           1000|Theta_4   | -1.4125413| 0.1057238|   37.83444|     2.100891|
# Combinar resultados en un solo dataframe
resultados_mat <- do.call(rbind, lapply(resultados_list, function(x) x$estimadores))

# Crear boxplot con ggplot
ggplot(resultados_mat, aes(x = factor(Tamaño_Muestra), y = Valor, fill = Estimador)) +
  geom_boxplot() +
  geom_hline(yintercept = 0.5, color= "red")+
  labs(title = "Comparación Estimadores por tamaño de muestra",
       x = "Tamaño de Muestra",
       y = "Valor Estimado",
       fill = "Estimador") +
  theme_minimal()

  1. Sesgo y Varianza:
    • Theta_1, Theta_3 y Theta_4 muestran sesgos cercanos a -1.5 en la mayoría de los casos, mientras que Theta_2 parece estar más cerca de -1. Sin embargo, al considerar que el parámetro conocido es 2, estos sesgos aún indican una subestimación constante.
    • La varianza sigue siendo una medida importante de la dispersión de los estimadores. Theta_2 tiene la varianza más alta en todos los tamaños de muestra, lo que sugiere que es menos consistente en sus estimaciones en comparación con los otros estimadores.
  2. Eficiencia:
    • Los valores de eficiencia más altos indican un mejor rendimiento del estimador. Aunque el sesgo puede estar sesgado hacia abajo, la eficiencia puede variar dependiendo de la varianza. En este caso, Theta_2 muestra la menor eficiencia, lo que indica que aunque tiene un sesgo cercano al verdadero parámetro, su alta varianza reduce su precisión relativa.
  3. Consistencia:
    • La consistencia de Theta_1, Theta_3 y Theta_4 se mantiene, ya que sus sesgos tienden a cero a medida que aumenta el tamaño de la muestra. Sin embargo, para Theta_2, su sesgo parece no converger a cero a medida que aumenta el tamaño de la muestra, lo que podría sugerir problemas de consistencia.
En resumen, aunque los estimadores muestran sesgos y varianzas que pueden no ser óptimos en relación con el parámetro conocido, algunos son más eficientes y consistentes que otros. Theta_2 parece ser el estimador menos deseable en general, mientras que Theta_1, Theta_3 y Theta_4 muestran mejores propiedades en términos de eficiencia y consistencia, aunque pueden tener sesgos hacia abajo consistentes.

PROBLEMA 3: Teorema del Límite Central

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

Pasos:

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

  2. Genere una función que permita:

  • Obtener una muestra aleatoria de la población y
  • Calcule el estimador de la proporción muestral \(\hat{p}\) para un tamaño de muestra dado n.

  1. Repita el escenario anterior (b) n = 500 veces y analice los resultados en cuanto al comportamiento de los 500 resultados del estimador \(\hat{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.

  2. 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 denormalidad: qqnorm()). Comente en su informe los resultados obtenidos.

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

Solución

  • Generación de la población inicial: Se crea una población de tamaño 1000 con igual cantidad de plantas sanas y enfermas.

  • Función para obtener una muestra aleatoria y calcular el estimador de proporción muestral: Se define una función llamada muestra_estimador que toma una muestra aleatoria de la población y calcula la proporción de plantas enfermas en esa muestra.

  • Realización del escenario para 500 repeticiones: Se realiza el muestreo aleatorio 500 veces, cada vez con una muestra de tamaño 100, y se almacenan las estimaciones de proporción en un vector llamado resultados.

  • Análisis de los resultados: Se grafican un histograma y un boxplot para visualizar la distribución de las estimaciones de proporción. Se calculan medidas descriptivas como la media, la mediana, la desviación estándar, la asimetría y el coeficiente de variación de las estimaciones.

# Generar población de tamaño N = 1000
set.seed(123)  # Semilla para reproducibilidad
poblacion <- rep(c(0, 1), times = c(500, 500))  # 0 para sanos, 1 para enfermos
tp<-table(poblacion)
names(tp) <- c("Sanos", "Enfermos")
tp
##    Sanos Enfermos 
##      500      500
# Función para obtener muestra aleatoria y calcular estimador de proporción muestral
muestra_estimador <- function(poblacion, n) {
  muestra <- sample(poblacion, size = n, replace = FALSE)
  proporcion_muestral <- mean(muestra)
  return(proporcion_muestral)
}

# Realizar el escenario para 500 repeticiones
n_repeticiones <- 500
resultados <- replicate(n_repeticiones, muestra_estimador(poblacion, n = 100)) # Muestra de tamaño 100

# Análisis de los resultados
# Histograma y boxplot para analizar la simetría y la variabilidad
par(mfrow = c(1, 2))  # Configuramos un lienzo para mostrar múltiples gráficos
hist(resultados, main = "Distribución de las estimaciones p", xlab = "Estimación de proporción", col = "lightblue")
abline(v = mean(resultados), col = "red", lwd = 2) # Línea roja para la media
abline(v = median(resultados), col = "blue", lwd = 2) # Línea azul para la mediana
boxplot(resultados, main = "Distribución de las estimaciones p", ylab = "Estimación de proporción", col = "lightgreen")

# Crear un dataframe con los resultados
tabla_resultados <- data.frame(
  Media = mean(resultados),
  Mediana = median(resultados),
  `Desv Estándar` = sd(resultados),
  Asimetria = skewness(resultados), #Cercana a cero sugiere que es aproximadamente simétrica
  `Coef Variación` = paste(round((sd(resultados) / mean(resultados)) * 100, 6), "%", sep = "")  #Se considera homogéneo cuando es <= 20%
)

# Imprimir la tabla usando kable
kable(tabla_resultados, caption = "Resumen de resultados de los estimadores de proporción")
Resumen de resultados de los estimadores de proporción
Media Mediana Desv.Estándar Asimetria Coef.Variación
0.50002 0.5 0.0475807 0.1010258 9.515769%

Análisis:

  • La media es aproximadamente 0.5, lo que significa que, en promedio las estimaciones de proporción muestral están muy cerca del verdadero valor de proporción en la población, que es 0.5 (mitad de la población enferma y mitad sana).
  • La desviación estándar es aproximadamente 0.0475807, lo que sugiere que las estimaciones tienden a variar alrededor de 0.5, pero no demasiado lejos de este valor.
  • Un valor cercano a cero indica simetría en la distribución. En este caso, la asimetría es 0.1010258, lo que sugiere una ligera asimetría hacia la derecha, pero en general, la distribución es bastante simétrica.
  • El coeficiente de variación es 9.52%, lo que indica que la variabilidad relativa de las estimaciones en relación con su media es relativamente baja, lo que sugiere una cierta homogeneidad en las estimaciones.
  • Simulación de estimaciones para diferentes tamaños de muestra: Se realizan estimaciones de proporción para diferentes tamaños de muestra. Para cada tamaño de muestra, se grafica un QQ plot para evaluar la normalidad de las estimaciones. 
## Simulación de estimaciones para diferentes tamaños de muestra
resultados_lista <- list() # Lista para almacenar los resultados de cada muestra
tamanos_muestra <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500) # Definir tamaños de muestra

par(mfrow = c(2, 5)) # Crear lienzo de 2 filas y 5 columnas para qqnorm

for (n in tamanos_muestra) {
  resultados <- replicate(n_repeticiones, muestra_estimador(poblacion, n))
  resultados_lista[[as.character(n)]] <- resultados
  
  shapiro_test <- shapiro.test(resultados)
  
  # QQ Plot
  qqnorm(resultados, main = paste("QQ Plot para n =", n))
  qqline(resultados)
  
  print(paste("Shapiro-Wilk p-value para n =", n, ":", shapiro_test$p.value))
}
## [1] "Shapiro-Wilk p-value para n = 5 : 2.36258457319825e-14"
## [1] "Shapiro-Wilk p-value para n = 10 : 1.14110726364453e-09"
## [1] "Shapiro-Wilk p-value para n = 15 : 8.00760357321995e-08"
## [1] "Shapiro-Wilk p-value para n = 20 : 6.06941777494782e-06"
## [1] "Shapiro-Wilk p-value para n = 30 : 5.64374191954771e-05"
## [1] "Shapiro-Wilk p-value para n = 50 : 0.00114555325539717"
## [1] "Shapiro-Wilk p-value para n = 60 : 0.00675426182657455"
## [1] "Shapiro-Wilk p-value para n = 100 : 0.0190050162511047"
## [1] "Shapiro-Wilk p-value para n = 200 : 0.111690146947713"

## [1] "Shapiro-Wilk p-value para n = 500 : 0.168216156544185"
# Convertir la lista de resultados a un marco de datos
resultados_df <- data.frame(
  Tamaño_Muestra = rep(tamanos_muestra, each = n_repeticiones),
  Estimacion = unlist(resultados_lista)
)

# Crear boxplot comparativo con ggplot
ggplot(resultados_df, aes(x = factor(Tamaño_Muestra), y = Estimacion)) +
  geom_boxplot(fill = "lightblue") +
  labs(title = "Comparación de Estimaciones por Tamaño de Muestra",
       x = "Tamaño de Muestra",
       y = "Estimación de Proporción") +
  theme_minimal()

Análisis:

-Estos resultados indican que las estimaciones de proporción muestral pueden no seguir una distribución normal, especialmente para tamaños de muestra pequeños, pero a medida que aumenta el tamaño de la muestra, los datos tienden a aproximarse más a la normalidad. A medida que aumenta n, los valores p también aumentan gradualmente, pero siguen siendo significativamente pequeños hasta n = 100. Para n = 200 y n = 500, los valores p son mayores que el nivel de significancia de 0.05, lo que sugiere que no hay suficiente evidencia para rechazar la hipótesis nula de normalidad.

  • Simulación de estimaciones para lotes con diferentes proporciones de plantas enfermas: Se realiza una simulación similar a la anterior, pero esta vez con diferentes proporciones de plantas enfermas en la población. 
## Simulación de estimaciones para lotes con 10% y 90% de plantas enfermas
# Definir parámetros
prop_enfermas <- c(0.1, 0.9)
resultados_lista <- list()

par(mfrow = c(2, 5)) # Configurar el lienzo para mostrar múltiples gráficos
contador <- 1 # Contador para controlar la posición del gráfico en la cuadrícula

for (prop in prop_enfermas) {
  set.seed(123)  # Semilla para reproducibilidad
  poblacion <- rep(c(0, 1), times = c(1000 * (1 - prop), 1000 * prop))  # Ajuste de proporción de plantas enfermas
  
  resultados_prop <- list()
  
  for (n in tamanos_muestra) {
    resultados <- replicate(n_repeticiones, muestra_estimador(poblacion, n))
    resultados_prop[[as.character(n)]] <- resultados
    
    shapiro_test <- shapiro.test(resultados)
    
    # QQ Plot
    # Determinar posición en la cuadrícula
    fila <- ceiling(contador / 5)
    columna <- contador %% 5
    if (columna == 0) columna <- 5
    
    # Ajustar marco de gráfico
    par(mar = c(5, 4, 4, 2))
    
    # Generar QQ plot en la posición especificada
    qqnorm(resultados, main = paste("n =", n, "y proporción =", prop),
           xlab = "Cuantiles teóricos", ylab = "Cuantiles observados")
    qqline(resultados)
    
    # Incrementar el contador
    contador <- contador + 1
    
    print(paste("Shapiro-Wilk p-value para n =", n, "y proporción de enfermos =", prop, ":", shapiro_test$p.value))
  }
  
  resultados_lista[[as.character(prop)]] <- resultados_prop
}
## [1] "Shapiro-Wilk p-value para n = 5 y proporción de enfermos = 0.1 : 3.59812956053355e-28"
## [1] "Shapiro-Wilk p-value para n = 10 y proporción de enfermos = 0.1 : 1.61347519086642e-22"
## [1] "Shapiro-Wilk p-value para n = 15 y proporción de enfermos = 0.1 : 1.98228169949802e-17"
## [1] "Shapiro-Wilk p-value para n = 20 y proporción de enfermos = 0.1 : 1.00272386857994e-14"
## [1] "Shapiro-Wilk p-value para n = 30 y proporción de enfermos = 0.1 : 1.81714179168535e-11"
## [1] "Shapiro-Wilk p-value para n = 50 y proporción de enfermos = 0.1 : 4.03922222645246e-08"
## [1] "Shapiro-Wilk p-value para n = 60 y proporción de enfermos = 0.1 : 2.59095689627499e-08"
## [1] "Shapiro-Wilk p-value para n = 100 y proporción de enfermos = 0.1 : 0.00014125987500612"
## [1] "Shapiro-Wilk p-value para n = 200 y proporción de enfermos = 0.1 : 0.00200107337274328"

## [1] "Shapiro-Wilk p-value para n = 500 y proporción de enfermos = 0.1 : 0.00616971540326558"
## [1] "Shapiro-Wilk p-value para n = 5 y proporción de enfermos = 0.9 : 5.41820993822321e-28"
## [1] "Shapiro-Wilk p-value para n = 10 y proporción de enfermos = 0.9 : 5.34996180890682e-22"
## [1] "Shapiro-Wilk p-value para n = 15 y proporción de enfermos = 0.9 : 1.4660813448635e-17"
## [1] "Shapiro-Wilk p-value para n = 20 y proporción de enfermos = 0.9 : 3.17294926773443e-16"
## [1] "Shapiro-Wilk p-value para n = 30 y proporción de enfermos = 0.9 : 7.77839631959109e-12"
## [1] "Shapiro-Wilk p-value para n = 50 y proporción de enfermos = 0.9 : 3.64015232331105e-10"
## [1] "Shapiro-Wilk p-value para n = 60 y proporción de enfermos = 0.9 : 2.13823235869883e-07"
## [1] "Shapiro-Wilk p-value para n = 100 y proporción de enfermos = 0.9 : 2.87228898635745e-06"
## [1] "Shapiro-Wilk p-value para n = 200 y proporción de enfermos = 0.9 : 0.00868569234736011"

## [1] "Shapiro-Wilk p-value para n = 500 y proporción de enfermos = 0.9 : 0.0348777628589641"
# Reiniciar el parámetro mfrow
par(mfrow = c(1, 1))

# Convertir la lista de resultados a un marco de datos
resultados_df <- data.frame(
  Proporcion_Enfermos = rep(prop_enfermas, each = length(tamanos_muestra) * n_repeticiones),
  Tamaño_Muestra = rep(tamanos_muestra, times = length(prop_enfermas) * n_repeticiones),
  Estimacion = unlist(unlist(resultados_lista))
)

# Crear boxplots comparativos para cada proporción de enfermos
ggplot(resultados_df, aes(x = factor(Tamaño_Muestra), y = Estimacion, fill = factor(Proporcion_Enfermos))) +
  geom_boxplot() +
  labs(title = "Comparación de Estimaciones por Tamaño de Muestra y Proporción de Plantas Enfermas",
       x = "Tamaño de Muestra",
       y = "Estimación de Proporción",
       fill = "Proporción de Enfermos") +
  theme_minimal()

Para proporción de enfermos = 0.1:

  • Los valores p son pequeños en todos los tamaños de muestra, lo que indica evidencia significativa para rechazar la hipótesis nula de normalidad. Por ejemplo, para n = 5, el valor p es aproximadamente 3.598e-28.

  • Estos resultados sugieren que las estimaciones de proporción muestral para una proporción de enfermos del 10% no siguen una distribución normal, independientemente del tamaño de la muestra.

Para proporción de enfermos = 0.9:

  • Los valores p también son muy pequeños en general, pero no tan pequeños como en el caso anterior. A medida que aumenta el tamaño de la muestra, los valores p tienden a aumentar, pero siguen siendo significativamente pequeños.

  • Sin embargo, para algunos tamaños de muestra más grandes (n = 200, n = 500), los valores p son mayores que el nivel de significancia de 0.05, lo que sugiere que no hay suficiente evidencia para rechazar la hipótesis nula de normalidad.

  • Estos resultados sugieren que las estimaciones de proporción muestral para una proporción de enfermos del 90% no siguen una distribución normal para tamaños de muestra más pequeños, pero pueden aproximarse a una distribución normal para tamaños de muestra más grandes.

En resumen, estos resultados indican que la normalidad de las estimaciones de proporción muestral depende tanto del tamaño de la muestra como de la proporción de plantas enfermas en la población, siendo más probable que las estimaciones sean no normales para proporciones de enfermos más extremas y tamaños de muestra más pequeños.

PROBLEMA 4: Estimacció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 %. 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 \(X^{*}_{1}\). Después de anotado el valor se regresa \(X^{*}_{1}\) a la caja y se extrae el valor \(X^{*}_{2}\), regresándolo nuevamente. Este procedimiento se repite hasta completar una muestra de tamaño n = 5, \(X^{*}_{1}\), \(X^{*}_{2}\), \(X^{*}_{2}\), \(X^{*}_{n}\), 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 \(\bar{X}^{*}_{i}\), obteniéndose un valor para cada muestra. El intervalo de confianza queda conformado por los percentiles \(P_{2.5}\) y \(P_{97.5}\). Existen dos métodos para estimarlo:

Método 1: (\(P_{2.5}\); \(P_{97.5}\))
Método 2: (2\(\bar{X}\) - \(P_{97.5}\); 2\(\bar{X}\) – \(P_{2.5}\))

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

Solución

Se realiza un análisis de remuestreo para estimar el intervalo de confianza de la media de una muestra de datos de eficiencia de combustible.

  • Definición de datos: Se establece un vector que contiene los valores de eficiencia de combustible. En la variable n se almacena el tamaño de la muestra y se define k como el número de muestras bootstrap que se generarán.

  • Función para generar una muestra bootstrap: La función generar_muestra_bootstrap toma una muestra aleatoria con reemplazo de los datos y devuelve la media de esa muestra.

  • Generación de muestras bootstrap: Se utilizan muestras bootstrap para estimar la distribución de la media muestral. Esto se realiza mediante la función replicate, que genera k muestras utilizando la función generar_muestra_bootstrap.

  • Cálculo del intervalo de confianza utilizando dos métodos:

    • Método 1: Se calculan los cuantiles 2.5% y 97.5% de las muestras bootstrap.
    • Método 2: Se calcula el intervalo de confianza a partir de la media bootstrap y los cuantiles obtenidos en el método 1.
# Datos de eficiencia de combustible
datos <- c(7.69, 4.97, 4.56, 6.49, 4.34, 6.24, 4.45)

# Tamaño de la muestra
n <- length(datos)

# Número de muestras bootstrap
k <- 1000

# Función para generar una muestra bootstrap
generar_muestra_bootstrap <- function(datos) {
  muestra_bootstrap <- sample(datos, replace = TRUE)
  return(mean(muestra_bootstrap))
}

# Generar k muestras bootstrap
muestras_bootstrap <- replicate(k, generar_muestra_bootstrap(datos))

# Calcular intervalo de confianza utilizando el método 1
intervalo_metodo1 <- quantile(muestras_bootstrap, c(0.025, 0.975))

# Calcular intervalo de confianza utilizando el método 2
media_bootstrap <- mean(muestras_bootstrap)
intervalo_metodo2 <- c(2 * media_bootstrap - intervalo_metodo1[2], 
                       2 * media_bootstrap - intervalo_metodo1[1])

# Mostrar resultados
cat("Intervalo de confianza Método 1:", intervalo_metodo1, "\n")
## Intervalo de confianza Método 1: 4.721429 6.440107
cat("Intervalo de confianza Método 2:", intervalo_metodo2, "\n")
## Intervalo de confianza Método 2: 4.599264 6.317943
  • Se imprimen en la consola los intervalos de confianza obtenidos por ambos métodos y se construye un histograma que muestra la distribución de las medias bootstrap y los intervalos de confianza calculados por ambos métodos son superpuestos.
summary(muestras_bootstrap)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   4.434   5.225   5.514   5.520   5.803   6.959
# Crear histograma
hist(muestras_bootstrap, breaks = 30, main = "Distribución de Medias Bootstrap",
     xlab = "Media Bootstrap", ylab = "Frecuencia", col = "lightblue")

# Agregar intervalos de confianza al histograma
abline(v = intervalo_metodo1, col = "red", lty = 2)
abline(v = intervalo_metodo2, col = "blue", lty = 2)

# Agregar leyenda
legend("topright", legend = c("Método 1", "Método 2"), 
       fill = c("red", "blue"), bty = "n")

- Ambos intervalos de confianza son similares, con ligeras diferencias en sus extremos.

  • Ambos intervalos incluyen la media de las medias bootstrap 5.54.Los valores mínimo y máximo de las medias bootstrap también están dentro de los intervalos de confianza calculados.

  • El Método 2 parece proporcionar un intervalo de confianza ligeramente más estrecho en comparación con el Método 1.Esto sugiere ser mas consistente y estable en múltiples muestras bootstrap.

  • Se puede inferir que los intervalos son confiables en el sentido de que proporcionan una estimación precisa de la variabilidad de la media muestral.