#install.packages("kableExtra")
#install.packages("moments")
#install.packages("gganimate")
#install.packages("gifski")
#install.packages("ggplot2")
#install.packages("magick")
#install.packages("reshape2")
#install.packages("plotly")

Problema 3: Teorema del Límite Central

Se comprobará el teorema del Limite Central a través de una simulación. Primero, se genera una población de \(n = 1000\) individuos (plantas), usando una distribución binomial, donde se medirán como casos de éxito que dichas plantas se encuentren enfermas, cuyo porcentaje será del 50%.

Se escogerá una muestra aleatoria de 50 individuos y luego se proceden a realizar 500 simulaciones. Dicha información se recopila en un histograma para evaluar su simetría y variabilidad.

set.seed(666)
library(moments)

# Generar una población de 1000 individuos, donde 0 = sano y 1 = enfermo usando rbinom()
n_poblacion <- 1000
poblacion <- rbinom(n_poblacion, size = 1, prob = 0.5)

calcular_estimador <- function(poblacion, n_muestra) {
  muestra <- sample(poblacion, size = n_muestra, replace = FALSE)
  p_hat <- mean(muestra)
  return(p_hat)
}

n_simulaciones <- 500
n_muestra <- 50

# Usar sapply para repetir la simulación
estimadores <- sapply(1:n_simulaciones, function(x) calcular_estimador(poblacion, n_muestra))

# Convertir resultados en data.frame para mejor manejo
resultados_df <- data.frame(estimadores)
#print(resultados_df)

# Analizando los resultados
#summary(resultados_df)
kurt_value <- kurtosis(resultados_df$estimadores)
hist(resultados_df$estimadores, main = paste("Distribución del estimador pˆ para n =", n_muestra),
     xlab = "pˆ", breaks = 20, col = "lightblue")
Fig 1. Distribución del estimador \(\hat{p}=0.5\) para \(n=50\) 
# Media y sesgo
media_estimador <- mean(estimadores)
sesgo <- media_estimador - 0.5  # Porque sabemos que la verdadera proporción es 0.5

# Desviación estándar (variabilidad)
desviacion <- sd(estimadores)

cat("<h4><strong>Simetría y variabilidad:</strong></h4><ul>",
    "<li>Media del estimador: ", media_estimador, "</li>",
    "<li>Sesgo: ", sesgo, "</li>",
    "<li>Desviacion estandar: ", desviacion, "</li>",
    "</ul>")

Simetría y variabilidad:


Según se evidencia, la gráfica tiene una figura símetrica que tiende a la de una distribución normal (con forma de campana). Y los valores calculados evidencian que el estimador es bastante cercano al valor real, con un sesgo pequeño y una desviación estandar relativamente pequeña.

Para evaluar la normalidad de los datos, se plantea la \(H_0\): los datos provienen de una población normalmente distribuida. Por tanto, se realizarán simulaciones con distinos tamaños de muestra \((5, 10, 15, 20, 30, 50, 60, 100, 200, 500)\) y se analizarán los \(p\)-\(values\) y los estadísticos de la prueba de normalidad \(Shapiro\)-\(Wilk\), teniendo en cuenta que si el \(p\)-\(value\)\(<0.05\), se rechaza la hipótesis nula.

# Cargar las librerías necesarias
library(moments)
library(magick)

# Generar una población de 1000 individuos, donde 0 = sano y 1 = enfermo usando rbinom()
set.seed(666)
n_poblacion <- 1000
poblacion1 <- rbinom(n_poblacion, size = 1, prob = 0.5)

# Función para calcular el estimador
calcular_estimador <- function(poblacion1, n_muestra) {
  muestra <- sample(poblacion1, size = n_muestra, replace = FALSE)
  p_hat <- mean(muestra)
  return(p_hat)
}

# Parámetros
n_simulaciones <- 500
n_muestras <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)
imagenes_histo05 <- list()  # Para almacenar las imágenes

# Bucle para generar las simulaciones y los gráficos
for (n_muestra in n_muestras) {
  
  # Usar sapply para repetir la simulación con distintos tamaños de muestra
  estimadores <- sapply(1:n_simulaciones, function(x) calcular_estimador(poblacion1, n_muestra))
  
  # Convertir resultados en data.frame para mejor manejo
  resultados_df <- data.frame(estimadores)
  
  # Crear histograma y guardar como archivo PNG
  hist_filename <- paste0("histograma_n_", n_muestra, ".png")
  png(hist_filename)
  hist(resultados_df$estimadores, 
       main = paste("Distribución del estimador pˆ para n =", n_muestra),
       xlab = "pˆ", breaks = 20, col = "lightblue", xlim=range(0,1))
  dev.off()
  
  # Cargar la imagen en el objeto imágenes
  imagenes_histo05[[length(imagenes_histo05) + 1]] <- image_read(hist_filename)
}
# Crear un GIF con las imágenes
gif <- image_animate(image_join(imagenes_histo05), fps = 1)  # 1 cuadro por segundo
image_write(gif, "histogramas05_simulaciones.gif")
Fig. 2 Distribución del estimador $\hat{p}=0.5$ para distintos tamaños de muestra.
Fig. 2 Distribución del estimador \(\hat{p}=0.5\) para distintos tamaños de muestra. 
set.seed(666)

n_muestras <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)

resultados <- list()
plots_1 <- list()
promedios_estimaciones <- numeric(length(n_muestras))
imagenes_qq05 <- list()

for (i in seq_along(n_muestras)) {
  n <- n_muestras[i]
  
  estimadores <- sapply(1:n_simulaciones, function(x) calcular_estimador(poblacion1, n))
  
  
  # Prueba de Shapiro-Wilk
  shapiro_result <- shapiro.test(estimadores)
  
  # Determinar si se rechaza o no la hipótesis nula (Ho)
  ho_result <- ifelse(shapiro_result$p.value < 0.05, "Rechazo", "No Rechazo")
  
  # Almacenar resultados en un data.frame
  resultados[[paste("n =", n)]] <- data.frame(
    n = n,
    shapiro_W = shapiro_result$statistic,
    p_value = shapiro_result$p.value,
    Ho = ho_result
  )
  
  # Gráfico QQ plot
  img_filename  <- paste0("qqplot_n_", n, "_simulaciones_", n_simulaciones, ".png")
  png(filename = img_filename)
  qqnorm(estimadores, main = paste("Normal Q-Q Plot para n=", n, "y i =", n_simulaciones), ylim = c(0, 1))
  qqline(estimadores, col="red")
  dev.off()
}

names(promedios_estimaciones) <- paste("n =", n_muestras)

# Crear un data frame con todos los resultados
resultados_df <- do.call(rbind, resultados)

# Mostrar la tabla con kable
kable(resultados_df, format = "html", row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = F) %>%
  add_header_above(c("Resultados para distintos tamaños de muestras (p = 0.5)" = 4)) %>%
  row_spec(0, bold = T) %>%
  column_spec(1:4, border_left = TRUE, border_right = TRUE)
Resultados para distintos tamaños de muestras (p = 0.5)
n shapiro_W p_value Ho
5 0.9296863 0.0000000 Rechazo
10 0.9610802 0.0000000 Rechazo
15 0.9757311 0.0000002 Rechazo
20 0.9806572 0.0000034 Rechazo
30 0.9858755 0.0000898 Rechazo
50 0.9876925 0.0003187 Rechazo
60 0.9918352 0.0075759 Rechazo
100 0.9939519 0.0438787 Rechazo
200 0.9940045 0.0458614 Rechazo
500 0.9962080 0.2792331 No Rechazo 
# Crear un GIF con las imágenes
gif <- image_animate(image_join(imagenes_qq05), fps = 1)  # 1 cuadro por segundo
image_write(gif, "qqplots05_simulaciones.gif")

 

De acuerdo con los resultados de la tabla, se evidencia que el estadístico de la prueba de \(Shapiro\)-\(Wilk\) tiende a 1, lo que indica que al aumentar el tamaño de la muestra, los datos simulados parecen provenir de una distribución normal. Por otro lado, en relación al \(p\)-\(value\), la hipótesis nula se rechaza para todos los tamaños de muestra, exceptuando cuando \(n=500\).

A continuación, se puede observar un gráfico QQ-Plot con las 500 simulaciones para los distintos tamaños de muestras.

Fig. 3 QQ Plots de Estimadores para Diferentes Tamaños de Muestra.
Fig. 3 QQ Plots de Estimadores para Diferentes Tamaños de Muestra.

 

Ahora, se realizará la misma simulación anterior pero asumiendo un \(p=0.9\) y un \(p=0.1\).

P = 0.9

# Generar una población de 1000 individuos, donde 0 = sano y 1 = enfermo usando rbinom()
set.seed(666)
n_poblacion <- 1000
poblacion2 <- rbinom(n_poblacion, size = 1, prob = 0.9)

# Función para calcular el estimador
calcular_estimador <- function(poblacion2, n_muestra) {
  muestra <- sample(poblacion2, size = n_muestra, replace = FALSE)
  p_hat <- mean(muestra)
  return(p_hat)
}

# Parámetros
n_simulaciones <- 500
n_muestras <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)
imagenes_histo09 <- list()  # Para almacenar las imágenes

# Bucle para generar las simulaciones y los gráficos
for (n_muestra in n_muestras) {
  
  # Usar sapply para repetir la simulación con distintos tamaños de muestra
  estimadores <- sapply(1:n_simulaciones, function(x) calcular_estimador(poblacion2, n_muestra))
  
  # Convertir resultados en data.frame para mejor manejo
  resultados_df <- data.frame(estimadores)
  
  # Crear histograma y guardar como archivo PNG
  hist_filename <- paste0("histograma_n_", n_muestra, ".png")
  png(hist_filename)
  hist(resultados_df$estimadores, 
       main = paste("Distribución del estimador pˆ para n =", n_muestra),
       xlab = "pˆ", breaks = 20, col = "lightblue", xlim=range(0.5,1))
  dev.off()
  
  # Cargar la imagen en el objeto imágenes
  imagenes_histo09[[length(imagenes_histo09) + 1]] <- image_read(hist_filename)
}
# Crear un GIF con las imágenes
gif <- image_animate(image_join(imagenes_histo09), fps = 1)  # 1 cuadro por segundo
image_write(gif, "histogramas09_simulaciones.gif")
Fig. 4 Distribución del estimador $\hat{p}=0.9$ para distintos tamaños de muestra.
Fig. 4 Distribución del estimador \(\hat{p}=0.9\) para distintos tamaños de muestra. 
set.seed(666)

n_muestras <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)

resultados <- list()
plots_1 <- list()
promedios_estimaciones <- numeric(length(n_muestras))
imagenes_qq09 <- list()

for (i in seq_along(n_muestras)) {
  n <- n_muestras[i]
  
  estimadores <- sapply(1:n_simulaciones, function(x) calcular_estimador(poblacion2, n))
  
  # Prueba de Shapiro-Wilk
  shapiro_result <- shapiro.test(estimadores)
  
  # Determinar si se rechaza o no la hipótesis nula (Ho)
  ho_result <- ifelse(shapiro_result$p.value < 0.05, "Rechazo", "No Rechazo")
  
  # Almacenar resultados en un data.frame
  resultados[[paste("n =", n)]] <- data.frame(
    n = n,
    shapiro_W = shapiro_result$statistic,
    p_value = shapiro_result$p.value,
    Ho = ho_result
  )
  
  # Gráfico QQ plot
  img_filename  <- paste0("qqplot_n_", n, "_simulaciones_", n_simulaciones, ".png")
  png(filename = img_filename)
  qqnorm(estimadores, main = paste("Normal Q-Q Plot para n=", n, "y i =", n_simulaciones), ylim = c(0.5, 1))
  qqline(estimadores, col="red")
  dev.off()
  
  imagenes_qq09[[i]] <- image_read(img_filename)
}

# Crear un data frame con todos los resultados
resultados_df <- do.call(rbind, resultados)

# Mostrar la tabla con kable
kable(resultados_df, format = "html", row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = F) %>%
  add_header_above(c("Resultados para distintos tamaños de muestras (p = 0.9)" = 4)) %>%
  row_spec(0, bold = T) %>%
  column_spec(1:4, border_left = TRUE, border_right = TRUE)
Resultados para distintos tamaños de muestras (p = 0.9)
n shapiro_W p_value Ho
5 0.6932055 0.0000000 Rechazo
10 0.8273326 0.0000000 Rechazo
15 0.8960695 0.0000000 Rechazo
20 0.9243164 0.0000000 Rechazo
30 0.9457173 0.0000000 Rechazo
50 0.9663841 0.0000000 Rechazo
60 0.9707061 0.0000000 Rechazo
100 0.9802342 0.0000026 Rechazo
200 0.9910954 0.0041792 Rechazo
500 0.9914608 0.0055981 Rechazo 
# Crear un GIF con las imágenes
gif <- image_animate(image_join(imagenes_qq09), fps = 1)  # 1 cuadro por segundo
image_write(gif, "qqplots09_simulaciones.gif")
Fig. 5 Normal QQ Plot para Diferentes Tamaños de Muestra.
Fig. 5 Normal QQ Plot para Diferentes Tamaños de Muestra.

 

P = 0.1

# Cargar las librerías necesarias
library(moments)
library(magick)

# Generar una población de 1000 individuos, donde 0 = sano y 1 = enfermo usando rbinom()
set.seed(666)
n_poblacion <- 1000
poblacion3 <- rbinom(n_poblacion, size = 1, prob = 0.1)

# Función para calcular el estimador
calcular_estimador <- function(poblacion3, n_muestra) {
  muestra <- sample(poblacion3, size = n_muestra, replace = FALSE)
  p_hat <- mean(muestra)
  return(p_hat)
}

# Parámetros
n_simulaciones <- 500
n_muestras <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)
imagenes_histo01 <- list()  # Para almacenar las imágenes

# Bucle para generar las simulaciones y los gráficos
for (n_muestra in n_muestras) {
  
  # Usar sapply para repetir la simulación con distintos tamaños de muestra
  estimadores <- sapply(1:n_simulaciones, function(x) calcular_estimador(poblacion3, n_muestra))
  
  # Convertir resultados en data.frame para mejor manejo
  resultados_df <- data.frame(estimadores)
  
  # Crear histograma y guardar como archivo PNG
  hist_filename <- paste0("histograma_n_", n_muestra, ".png")
  png(hist_filename)
  hist(resultados_df$estimadores, 
       main = paste("Distribución del estimador pˆ para n =", n_muestra),
       xlab = "pˆ", breaks = 20, col = "lightblue", xlim=range(0,0.5))
  dev.off()
  
  # Cargar la imagen en el objeto imágenes
  imagenes_histo01[[length(imagenes_histo01) + 1]] <- image_read(hist_filename)
}
# Crear un GIF con las imágenes
gif <- image_animate(image_join(imagenes_histo01), fps = 1)  # 1 cuadro por segundo
image_write(gif, "histogramas01_simulaciones.gif")
Fig. 6 Distribución del estimador $\hat{p}=0.1$ para distintos tamaños de muestra.
Fig. 6 Distribución del estimador \(\hat{p}=0.1\) para distintos tamaños de muestra. 
set.seed(666)

n_muestras <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)

resultados <- list()
plots_1 <- list()
promedios_estimaciones <- numeric(length(n_muestras))
imagenes_qq01 <- list()

for (i in seq_along(n_muestras)) {
  n <- n_muestras[i]
  
  estimadores <- sapply(1:n_simulaciones, function(x) calcular_estimador(poblacion3, n))
  
  # Prueba de Shapiro-Wilk
  shapiro_result <- shapiro.test(estimadores)
  
  # Determinar si se rechaza o no la hipótesis nula (Ho)
  ho_result <- ifelse(shapiro_result$p.value < 0.05, "Rechazo", "No Rechazo")
  
  # Almacenar resultados en un data.frame
  resultados[[paste("n =", n)]] <- data.frame(
    n = n,
    shapiro_W = shapiro_result$statistic,
    p_value = shapiro_result$p.value,
    Ho = ho_result
  )
  
  # Gráfico QQ plot
  img_filename  <- paste0("qqplot_n_", n, "_simulaciones_", n_simulaciones, ".png")
  png(filename = img_filename)
  qqnorm(estimadores, main = paste("Normal Q-Q Plot para n=", n, "y i =", n_simulaciones), ylim = c(0, 0.5))
  qqline(estimadores, col="red")
  dev.off()
  
  imagenes_qq01[[i]] <- image_read(img_filename)
}

# Crear un data frame con todos los resultados
resultados_df <- do.call(rbind, resultados)

# Mostrar la tabla con kable
kable(resultados_df, format = "html", row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = F) %>%
  add_header_above(c("Resultados para distintos tamaños de muestras (p = 0.1)" = 4)) %>%
  row_spec(0, bold = T) %>%
  column_spec(1:4, border_left = TRUE, border_right = TRUE)
Resultados para distintos tamaños de muestras (p = 0.1)
n shapiro_W p_value Ho
5 0.6932055 0.0000000 Rechazo
10 0.8273326 0.0000000 Rechazo
15 0.8960695 0.0000000 Rechazo
20 0.9243164 0.0000000 Rechazo
30 0.9457173 0.0000000 Rechazo
50 0.9663841 0.0000000 Rechazo
60 0.9707061 0.0000000 Rechazo
100 0.9802342 0.0000026 Rechazo
200 0.9910954 0.0041792 Rechazo
500 0.9914608 0.0055981 Rechazo 
# Crear un GIF con las imágenes
gif <- image_animate(image_join(imagenes_qq01), fps = 1)  # 1 cuadro por segundo
image_write(gif, "qqplots01_simulaciones.gif")
Fig. 7 Normal QQ Plot para Diferentes Tamaños de Muestra.
Fig. 7 Normal QQ Plot para Diferentes Tamaños de Muestra.

 

Cuando las probabilidades son \(p=0.9\) y \(p=0.1\), se evidencia que a pesar de que los estimadores de la prueba de normalidad \(Shapiro\)-\(Wilk\) tienden a \(1\) a medida que el tamaño de la muestra aumenta, para ambas distribuciones, se terminó rechazando la \(H_0\) independiente del tamaño de muestra. Lo que demuestra que a pesar de que los datos parecen ajustarse bien a una distribución normal, a través del \(p\)-\(value\) no se logra afirmar que los datos realmente sigan dicha distribución.