1. Planteamiento

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.

2. Metodología

  1. 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%.

  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 :shapiro.test()) y métodos gráficos (gráfico de normalidad: qqnorm()). Comente en su informe los resultados obtenidos.

  3. 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. Resultados

3.1. Lote con 50% de plantas enfermas.

Se realiza simulación de un lote de 1000 plantas donde el porcentaje de plantas enfermas es del 50% a partir de una generación de números aleatorios de una distrbución binomial con parámetros \(n=1000\) y \(p=0.5\).

# Paso 1

n_lote <- 1000
enfermas <- rbinom(n_lote, size = 1, prob = 0.5)

data_lote <- data.frame(id=1:n_lote, tipo = enfermas) # 1. planta enferma, 0: planta sana

Teniendo ya el lote de plantas generado, se crea una función que permite extraer muestras de algún tamaño deseado donde las salidas de la función es mostrar cuál fue la muestra escogida y cuál fue su respectivo \(\hat{p}\), como se ve a continuación, donde se pone de ejemplo la aplicabilidad de la función extrayendo una muestra de tamaño 20 al lote anteriormente mencionado.

#Paso 2

funcion_muestra <- function(data, n_muestra) {
  muestra <- sample(data_lote$tipo, n_muestra, replace = FALSE)
  p_estimado <- mean(muestra)
  return(list(muestra = muestra, p_estimado = p_estimado))
}

## Aplicamos un ejemplo para ver salida de la función:
ej <- funcion_muestra(data_lote, 20)

cat("La muestra extraída ", ej$muestra, " tiene un p estimado =", ej$p_estimado)
## La muestra extraída  1 1 1 0 1 0 1 1 0 0 1 1 1 1 0 1 1 0 0 0  tiene un p estimado = 0.6

La función anterior se repite 500 veces. Así, ahora se mostrará 500 salidas de muestras extraídas y se calcula a cada una su respectivo \(\hat{p}\). Estas salidas quedarán guardadas en un data frame compuesto por número de la iteración y el \(\hat{p}\).

#Paso 3

n_muestra <- 20
n_repeticiones <- 500

df_estimacion <- data.frame(iteracion = numeric(n_repeticiones), p_estimados = numeric(n_repeticiones))

for (i in 1:n_repeticiones) {
  df_estimacion$iteracion[i] <- i
  muestra_resultado <- funcion_muestra(data_lote, n_muestra)
  df_estimacion$p_estimados[i] <- muestra_resultado$p_estimado
}

Para responder a las preguntas de simetría, variabilidad y sesgo, se graficará un histograma, se calculará la simetría de los datos estimados, el sesgo y el coeficiente de variación.

#Paso 4

# Se realiza gráfico para mirar simetría de la distribución
library(ggplot2)
library(tidyverse)
library(hrbrthemes)

ggplot(df_estimacion, aes(x=p_estimados)) +
  geom_histogram(binwidth=0.05, fill="#d6a692", color="#e9ecef", alpha=0.9) +
  ggtitle(bquote("Distribución de " ~ hat(p) ~ "")) +
  labs(x=(bquote("" ~ hat(p) ~ "")), y="densidad") +
  theme_ipsum() +
  theme(
    plot.title = element_text(size=15)
  )

# Se calcula la simetría de la distribución
library(parameters)
simetria <- skewness(df_estimacion$p_estimados)


# Se calcula el sesgo, asumiendo que p=0.5 que es la probabilidad de que la planta esté enferma.
sesgo <- mean(df_estimacion$p_estimados) - 0.5


# Variabilidad de las estimaciones
cv <- sd(df_estimacion$p_estimados) / mean(df_estimacion$p_estimados)


cat("La simetría es igual a ", simetria$Skewness, "\n")
## La simetría es igual a  -0.1026486
cat("El sesgo es igual a ", sesgo, "\n")
## El sesgo es igual a  0.0223
cat("El coeficiente de variación es igual a ", cv)
## El coeficiente de variación es igual a  0.2107979

Observando los resultados anteriores y teniendo en cuenta que en cada muestra se extrajeron 20 plantas (n_muestra), se infiere que las estimaciones de \(\hat{p}\) presentan un comportamiento que tiende a la simetría, pues visualmente con el histograma los datos se distribuyen de manera Normal y al mirar la simetría calculada, ésta es muy cercana a cero, como también puede observarse con el coeficiente de variación, donde es bajo. Además, el sesgo es mínimo cuando se compara la media de las estimaciones con el \(p\) inicial del ejercicio.

Finalmente, los resultados pueden deberse al tamaño de muestra escogido, puesto que al acercarse a 30, presentan cierta similitud a una distribución Normal.

Ahora, se repetirá el proceso anterior para tamaños de muestra n=5, 10, 15, 20, 30, 50, 60, 100, 200, 500. Para esto, se crea una matriz que guarde 500 muestras generadas por columna y luego calcule el \(\hat{p}\) para cada una de los tamaños de muestra.

#Paso 5

# Generación de matriz
n=500

X <- matrix(unlist(lapply(1:n, function(i) funcion_muestra(enfermas, 100)$muestra)), ncol = n)

#Generación de medias
X5=X[ ,1:5]          # n=5
X10=X[ ,1:10]      # n=10
X20=X[ ,1:20]      # n=20
X30=X[ ,1:30]      # n=30
X50=X[ ,1:50]      # n=50
X60=X[ ,1:60]      # n=60
X100=X[ ,1:100]    # n=100
X200=X[ ,1:200]    # n=200
X500=X[ ,1:500]    # n=500

# Obtención de las medias para diferentes tamaños de muestra.
Mx5=apply(X5,1,mean)       # medias de muestras de tamaño n=5
Mx10=apply(X10,1,mean)      # medias de muestras de tamaño n=10
Mx20=apply(X20,1,mean)      # medias de muestras de tamaño n=20
Mx30=apply(X30,1,mean)      # medias de muestras de tamaño n=30
Mx50=apply(X50,1,mean)      # medias de muestras de tamaño n=50
Mx60=apply(X60,1,mean)      # medias de muestras de tamaño n=60
Mx100=apply(X100,1,mean)    # medias de muestras de tamaño n=100
Mx200=apply(X200,1,mean)    # medias de muestras de tamaño n=200
Mx500=apply(X500,1,mean)    # medias de muestras de tamaño n=500

Para revisar la Normalidad en los datos, se graficarán histogramas y qqnorms de cada uno de los tamaños de muestra, como también se calculará el valor p del test de Shapiro Wilk (SW) para detectar Normalidad en las estimaciones. En esta última, se asume como Hipótesis Nula (\(H_{0}\)) que las estimaciones siguen una distribución Normal, tomando un nivel de significancia del 5%.

# Gráficos de histogramas y qqnorm de cada uno de los tamaños de muestra
par(mfrow=c(3,6), mar = c(4.1, 3.1, 3.1, 1.1)) # divide la gráfica en una matriz 
hist(Mx5,  main = "n=5", freq=FALSE, xlab=expression(hat(p)), ylab="", col="#d6a692")
qqnorm(Mx5, main="n=5", xlab="cuantiles teóricos", ylab="") ; qqline(X5, col="red")
hist(Mx10, main = "n=10",freq=FALSE, xlab=expression(hat(p)), ylab="", col="#d6a692")
qqnorm(Mx10, main ="n=10", xlab="cuantiles teóricos", ylab="") ; qqline(Mx10, col="red")
hist(Mx20, main = "n=20",freq=FALSE, xlab=expression(hat(p)), ylab="", col="#d6a692")
qqnorm(Mx20, main ="n=20", xlab="cuantiles teóricos", ylab="") ; qqline(Mx20, col="red")
hist(Mx30, main = "n=30",freq=FALSE,xlab=expression(hat(p)), ylab="", col="#d6a692")
qqnorm(Mx30, main = "n=30", xlab="cuantiles teóricos", ylab="") ; qqline(Mx30, col="red")
hist(Mx50, main = "n=50",freq=FALSE,xlab=expression(hat(p)), ylab="", col="#d6a692")
qqnorm(Mx50, main = "n=50", xlab="cuantiles teóricos", ylab="") ; qqline(Mx50, col="red")
hist(Mx60, main = "n=60", freq = FALSE,xlab=expression(hat(p)), ylab="", col="#d6a692")
qqnorm(Mx60, main ="n=60", xlab="cuantiles teóricos", ylab="") ; qqline(Mx60, col="red")
hist(Mx100, main = "n=100", freq = FALSE,xlab=expression(hat(p)), ylab="", col="#d6a692")
qqnorm(Mx100, main ="n=100",xlab="cuantiles teóricos" ,ylab="") ; qqline(Mx100, col="red")
hist(Mx200, main = "n=200",freq = FALSE,xlab=expression(hat(p)), ylab="", col="#d6a692")
qqnorm(Mx200, main ="n=200", xlab="cuantiles teóricos", ylab="") ; qqline(Mx200, col="red")
hist(Mx500, main = "n=500", freq = FALSE,xlab=expression(hat(p)), ylab="", col="#d6a692")
qqnorm(Mx500, main="n=500", xlab="cuantiles teóricos", ylab="") ; qqline(Mx500, col="red")

# Tabla de cálculo de shapiro wilk a cada uno de los tamaños de muestra
sw5 <- shapiro.test(Mx5)
sw10 <- shapiro.test(Mx10)
sw20 <- shapiro.test(Mx20)
sw30 <- shapiro.test(Mx30)
sw50 <- shapiro.test(Mx50)
sw60 <- shapiro.test(Mx60)
sw100 <- shapiro.test(Mx100)
sw200 <- shapiro.test(Mx200)
sw500 <- shapiro.test(Mx500)

muestras <- c("n=5", "n=10", "n=20", "n=30", "n=50", "n=60", "n=100", "n=200", "n=500")
valores_p <- c(sw5$p.value, sw10$p.value, sw20$p.value, sw30$p.value, sw50$p.value, sw60$p.value, sw100$p.value, sw200$p.value, sw500$p.value)
tabla_p <- data.frame(Muestra = muestras, "Valor p" = valores_p)

library(gt)
## Warning: package 'gt' was built under R version 4.3.3
tabla_p %>%
  gt()
Muestra Valor.p
n=5 9.860471e-06
n=10 1.099317e-02
n=20 1.069761e-01
n=30 3.492623e-02
n=50 2.617026e-03
n=60 5.919059e-02
n=100 7.991699e-01
n=200 5.186192e-01
n=500 1.453016e-01

Observando los resultados anteriores, se evidencia que a medida que el tamaño de muestra aumenta, los histogramas tienden a ser más simétricos y que los puntos qqnorm siguen la línea de base. Además el test de SW a partir del tamaño \(n>50\) comienza a no rechazarse \(H_{0}\) al tener un valor p > 5%.

3.2. Lote con 10% de plantas enfermas.

Se realizarán los mismos pasos anteriores pero ahora generando un lote con 10% de plantas enfermas. A continuación se presentarán los histogramas, qqnorm y su respctiva prueba de shapiro wilk.

# Paso 1

n_lote <- 1000
p_10 <- 0.1
enfermas <- rbinom(n_lote, size = 1, prob = p_10)

data_lote <- data.frame(id=1:n_lote, tipo = enfermas) # 1. planta enferma, 0: planta sana

#---------------------------------------------------------------

#Paso 2

funcion_muestra <- function(data, n_muestra) {
  muestra <- sample(data_lote$tipo, n_muestra, replace = FALSE)
  p_estimado <- mean(muestra)
  return(list(muestra = muestra, p_estimado = p_estimado))
}

## Aplicamos un ejemplo para ver salida de la función:
ej <- funcion_muestra(data_lote, 20)

cat("La muestra extraída ", ej$muestra, " tiene un p estimado =", ej$p_estimado)
## La muestra extraída  0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0  tiene un p estimado = 0.1
#---------------------------------------------------------------

# Paso 3
n_muestra <- 20
n_repeticiones <- 500

df_estimacion <- data.frame(iteracion = numeric(n_repeticiones), p_estimados = numeric(n_repeticiones))

for (i in 1:n_repeticiones) {
  df_estimacion$iteracion[i] <- i
  muestra_resultado <- funcion_muestra(data_lote, n_muestra)
  df_estimacion$p_estimados[i] <- muestra_resultado$p_estimado
}

#---------------------------------------------------------------

#Paso 4

# Se realiza gráfico para mirar simetría de la distribución
library(ggplot2)
library(tidyverse)
library(hrbrthemes)

ggplot(df_estimacion, aes(x=p_estimados)) +
  geom_histogram(binwidth=0.05, fill="#3b8952", color="#e9ecef", alpha=0.9) +
  ggtitle(bquote("Distribución de " ~ hat(p) ~ "")) +
  labs(x=(bquote("" ~ hat(p) ~ "")), y="densidad") +
  theme_ipsum() +
  theme(
    plot.title = element_text(size=15)
  )

# Se calcula la simetría de la distribución
library(parameters)
simetria <- skewness(df_estimacion$p_estimados)


# Se calcula el sesgo, asumiendo que p=0.1 que es la probabilidad de que la planta esté enferma.
sesgo <- mean(df_estimacion$p_estimados) - p_10


# Variabilidad de las estimaciones
cv <- sd(df_estimacion$p_estimados) / mean(df_estimacion$p_estimados)


cat("La simetría es igual a ", simetria$Skewness, "\n")
## La simetría es igual a  0.8120182
cat("El sesgo es igual a ", sesgo, "\n")
## El sesgo es igual a  -0.004
cat("El coeficiente de variación es igual a ", cv)
## El coeficiente de variación es igual a  0.6776812
#---------------------------------------------------------------

#Paso 5

# Generación de matriz
n=500

X <- matrix(unlist(lapply(1:n, function(i) funcion_muestra(enfermas, 100)$muestra)), ncol = n)

#Generación de medias
X5=X[ ,1:5]          # n=5
X10=X[ ,1:10]      # n=10
X20=X[ ,1:20]      # n=20
X30=X[ ,1:30]      # n=30
X50=X[ ,1:50]      # n=50
X60=X[ ,1:60]      # n=60
X100=X[ ,1:100]    # n=100
X200=X[ ,1:200]    # n=200
X500=X[ ,1:500]    # n=500

# Obtención de las medias para diferentes tamaños de muestra.
Mx5=apply(X5,1,mean)       # medias de muestras de tamaño n=5
Mx10=apply(X10,1,mean)      # medias de muestras de tamaño n=10
Mx20=apply(X20,1,mean)      # medias de muestras de tamaño n=20
Mx30=apply(X30,1,mean)      # medias de muestras de tamaño n=30
Mx50=apply(X50,1,mean)      # medias de muestras de tamaño n=50
Mx60=apply(X60,1,mean)      # medias de muestras de tamaño n=60
Mx100=apply(X100,1,mean)    # medias de muestras de tamaño n=100
Mx200=apply(X200,1,mean)    # medias de muestras de tamaño n=200
Mx500=apply(X500,1,mean)    # medias de muestras de tamaño n=500

# Gráficos de histogramas y qqnorm de cada uno de los tamaños de muestra
par(mfrow=c(3,6), mar = c(4.1, 3.1, 3.1, 1.1)) # divide la gráfica en una matriz 
hist(Mx5,  main = "n=5", freq=FALSE, xlab=expression(hat(p)), ylab="", 
     col="#3b8952")
qqnorm(Mx5, main="n=5", xlab="cuantiles teóricos", ylab="") ; qqline(X5, col="red")
hist(Mx10, main = "n=10",freq=FALSE, xlab=expression(hat(p)), ylab="", 
     col="#3b8952")
qqnorm(Mx10, main ="n=10", xlab="cuantiles teóricos", ylab="") ; qqline(Mx10, col="red")
hist(Mx20, main = "n=20",freq=FALSE, xlab=expression(hat(p)), ylab="", 
     col="#3b8952")
qqnorm(Mx20, main ="n=20", xlab="cuantiles teóricos", ylab="") ; qqline(Mx20, col="red")
hist(Mx30, main = "n=30",freq=FALSE,xlab=expression(hat(p)), ylab="", 
     col="#3b8952")
qqnorm(Mx30, main = "n=30", xlab="cuantiles teóricos", ylab="") ; qqline(Mx30, col="red")
hist(Mx50, main = "n=50",freq=FALSE,xlab=expression(hat(p)), ylab="", 
     col="#3b8952")
qqnorm(Mx50, main = "n=50", xlab="cuantiles teóricos", ylab="") ; qqline(Mx50, col="red")
hist(Mx60, main = "n=60", freq = FALSE,xlab=expression(hat(p)), ylab="", 
     col="#3b8952")
qqnorm(Mx60, main ="n=60", xlab="cuantiles teóricos", ylab="") ; qqline(Mx60, col="red")
hist(Mx100, main = "n=100", freq = FALSE,xlab=expression(hat(p)), ylab="", 
     col="#3b8952")
qqnorm(Mx100, main ="n=100",xlab="cuantiles teóricos" ,ylab="") ; qqline(Mx100, col="red")
hist(Mx200, main = "n=200",freq = FALSE,xlab=expression(hat(p)), ylab="", 
     col="#3b8952")
qqnorm(Mx200, main ="n=200", xlab="cuantiles teóricos", ylab="") ; qqline(Mx200, col="red")
hist(Mx500, main = "n=500", freq = FALSE,xlab=expression(hat(p)), ylab="", 
     col="#3b8952")
qqnorm(Mx500, main="n=500", xlab="cuantiles teóricos", ylab="") ; qqline(Mx500, col="red")

# Tabla de cálculo de shapiro wilk a cada uno de los tamaños de muestra
sw5 <- shapiro.test(Mx5)
sw10 <- shapiro.test(Mx10)
sw20 <- shapiro.test(Mx20)
sw30 <- shapiro.test(Mx30)
sw50 <- shapiro.test(Mx50)
sw60 <- shapiro.test(Mx60)
sw100 <- shapiro.test(Mx100)
sw200 <- shapiro.test(Mx200)
sw500 <- shapiro.test(Mx500)

muestras <- c("n=5", "n=10", "n=20", "n=30", "n=50", "n=60", "n=100", "n=200", "n=500")
valores_p <- c(sw5$p.value, sw10$p.value, sw20$p.value, sw30$p.value, sw50$p.value, sw60$p.value, sw100$p.value, sw200$p.value, sw500$p.value)
tabla_p <- data.frame(Muestra = muestras, "Valor p" = valores_p)

library(gt)
tabla_p %>%
  gt()
Muestra Valor.p
n=5 3.725016e-13
n=10 6.787448e-10
n=20 5.623429e-07
n=30 1.093030e-03
n=50 7.067330e-03
n=60 5.580949e-04
n=100 1.146693e-02
n=200 1.055419e-01
n=500 2.970885e-02

Analizando el primer histograma de los \(\hat{p}\), no se observa que los datos sean simétricos, pues están concentrados a la izquierda del gráfico (asimetría positiva) y se puede identificar en el cálculo de la simetría, donde su valor es mayor que cero. Acompañado está el sesgo, que es relativamente bajo, significando que está muy cerca al valor real del parámetro \(p\). Sin embargo, la variabilidad de los datos refleja con el coeficiente de variación que ésta es un poco dispersa.

Luego, en cuanto a los análisis de los diferentes tamaños de muestra, la simetría en las distribuciones de los histogramas se alcanzan con un tamaño mayor que 200 (n>200), donde el qqnorm también se ajusta de mejor manera; inferencia que viene justificada junto a la prueba de SW, en donde el cumplimiento de \(H_{0}\) se da a partir del mencionado tamaño de muestra (\(valor p > 0.05\)).

3.3. Lote con 90% de plantas enfermas.

Ahora se hará generando un lote con 90% de plantas enfermas. A continuación se presentarán los histogramas, qqnorm y su respctiva prueba de shapiro wilk.

# Paso 1

n_lote <- 1000
p_90 <- 0.9
enfermas <- rbinom(n_lote, size = 1, prob = p_90)

data_lote <- data.frame(id=1:n_lote, tipo = enfermas) # 1. planta enferma, 0: planta sana

#---------------------------------------------------------------

#Paso 2

funcion_muestra <- function(data, n_muestra) {
  muestra <- sample(data_lote$tipo, n_muestra, replace = FALSE)
  p_estimado <- mean(muestra)
  return(list(muestra = muestra, p_estimado = p_estimado))
}

## Aplicamos un ejemplo para ver salida de la función:
ej <- funcion_muestra(data_lote, 20)

cat("La muestra extraída ", ej$muestra, " tiene un p estimado =", ej$p_estimado)
## La muestra extraída  0 1 1 1 1 1 1 0 1 1 1 1 1 0 1 0 1 1 1 1  tiene un p estimado = 0.8
#---------------------------------------------------------------

# Paso 3
n_muestra <- 20
n_repeticiones <- 500

df_estimacion <- data.frame(iteracion = numeric(n_repeticiones), p_estimados = numeric(n_repeticiones))

for (i in 1:n_repeticiones) {
  df_estimacion$iteracion[i] <- i
  muestra_resultado <- funcion_muestra(data_lote, n_muestra)
  df_estimacion$p_estimados[i] <- muestra_resultado$p_estimado
}

#---------------------------------------------------------------

#Paso 4

# Se realiza gráfico para mirar simetría de la distribución
library(ggplot2)
library(tidyverse)
library(hrbrthemes)

ggplot(df_estimacion, aes(x=p_estimados)) +
  geom_histogram(binwidth=0.05, fill="#5191c1", color="#e9ecef", alpha=0.9) +
  ggtitle(bquote("Distribución de " ~ hat(p) ~ "")) +
  labs(x=(bquote("" ~ hat(p) ~ "")), y="densidad") +
  theme_ipsum() +
  theme(
    plot.title = element_text(size=15)
  )

# Se calcula la simetría de la distribución
library(parameters)
simetria <- skewness(df_estimacion$p_estimados)


# Se calcula el sesgo, asumiendo que p=0.9 que es la probabilidad de que la planta esté enferma.
sesgo <- mean(df_estimacion$p_estimados) - p_90


# Variabilidad de las estimaciones
cv <- sd(df_estimacion$p_estimados) / mean(df_estimacion$p_estimados)


cat("La simetría es igual a ", simetria$Skewness, "\n")
## La simetría es igual a  -0.233768
cat("El sesgo es igual a ", sesgo, "\n")
## El sesgo es igual a  -0.0297
cat("El coeficiente de variación es igual a ", cv)
## El coeficiente de variación es igual a  0.08004219
#---------------------------------------------------------------

#Paso 5

# Generación de matriz
n=500

X <- matrix(unlist(lapply(1:n, function(i) funcion_muestra(enfermas, 100)$muestra)), ncol = n)

#Generación de medias
X5=X[ ,1:5]          # n=5
X10=X[ ,1:10]      # n=10
X20=X[ ,1:20]      # n=20
X30=X[ ,1:30]      # n=30
X50=X[ ,1:50]      # n=50
X60=X[ ,1:60]      # n=60
X100=X[ ,1:100]    # n=100
X200=X[ ,1:200]    # n=200
X500=X[ ,1:500]    # n=500

# Obtención de las medias para diferentes tamaños de muestra.
Mx5=apply(X5,1,mean)       # medias de muestras de tamaño n=5
Mx10=apply(X10,1,mean)      # medias de muestras de tamaño n=10
Mx20=apply(X20,1,mean)      # medias de muestras de tamaño n=20
Mx30=apply(X30,1,mean)      # medias de muestras de tamaño n=30
Mx50=apply(X50,1,mean)      # medias de muestras de tamaño n=50
Mx60=apply(X60,1,mean)      # medias de muestras de tamaño n=60
Mx100=apply(X100,1,mean)    # medias de muestras de tamaño n=100
Mx200=apply(X200,1,mean)    # medias de muestras de tamaño n=200
Mx500=apply(X500,1,mean)    # medias de muestras de tamaño n=500

# Gráficos de histogramas y qqnorm de cada uno de los tamaños de muestra
par(mfrow=c(3,6), mar = c(4.1, 3.1, 3.1, 1.1)) # divide la gráfica en una matriz 
hist(Mx5,  main = "n=5", freq=FALSE, xlab=expression(hat(p)), ylab="", 
     col="#5191c1")
qqnorm(Mx5, main="n=5", xlab="cuantiles teóricos", ylab="") ; qqline(X5, col="red")
hist(Mx10, main = "n=10",freq=FALSE, xlab=expression(hat(p)), ylab="", 
     col="#5191c1")
qqnorm(Mx10, main ="n=10", xlab="cuantiles teóricos", ylab="") ; qqline(Mx10, col="red")
hist(Mx20, main = "n=20",freq=FALSE, xlab=expression(hat(p)), ylab="", 
     col="#5191c1")
qqnorm(Mx20, main ="n=20", xlab="cuantiles teóricos", ylab="") ; qqline(Mx20, col="red")
hist(Mx30, main = "n=30",freq=FALSE,xlab=expression(hat(p)), ylab="", 
     col="#5191c1")
qqnorm(Mx30, main = "n=30", xlab="cuantiles teóricos", ylab="") ; qqline(Mx30, col="red")
hist(Mx50, main = "n=50",freq=FALSE,xlab=expression(hat(p)), ylab="", 
     col="#5191c1")
qqnorm(Mx50, main = "n=50", xlab="cuantiles teóricos", ylab="") ; qqline(Mx50, col="red")
hist(Mx60, main = "n=60", freq = FALSE,xlab=expression(hat(p)), ylab="", 
     col="#5191c1")
qqnorm(Mx60, main ="n=60", xlab="cuantiles teóricos", ylab="") ; qqline(Mx60, col="red")
hist(Mx100, main = "n=100", freq = FALSE,xlab=expression(hat(p)), ylab="", 
     col="#5191c1")
qqnorm(Mx100, main ="n=100",xlab="cuantiles teóricos" ,ylab="") ; qqline(Mx100, col="red")
hist(Mx200, main = "n=200",freq = FALSE,xlab=expression(hat(p)), ylab="", 
     col="#5191c1")
qqnorm(Mx200, main ="n=200", xlab="cuantiles teóricos", ylab="") ; qqline(Mx200, col="red")
hist(Mx500, main = "n=500", freq = FALSE,xlab=expression(hat(p)), ylab="", 
     col="#5191c1")
qqnorm(Mx500, main="n=500", xlab="cuantiles teóricos", ylab="") ; qqline(Mx500, col="red")

# Tabla de cálculo de shapiro wilk a cada uno de los tamaños de muestra
sw5 <- shapiro.test(Mx5)
sw10 <- shapiro.test(Mx10)
sw20 <- shapiro.test(Mx20)
sw30 <- shapiro.test(Mx30)
sw50 <- shapiro.test(Mx50)
sw60 <- shapiro.test(Mx60)
sw100 <- shapiro.test(Mx100)
sw200 <- shapiro.test(Mx200)
sw500 <- shapiro.test(Mx500)

muestras <- c("n=5", "n=10", "n=20", "n=30", "n=50", "n=60", "n=100", "n=200", "n=500")
valores_p <- c(sw5$p.value, sw10$p.value, sw20$p.value, sw30$p.value, sw50$p.value, sw60$p.value, sw100$p.value, sw200$p.value, sw500$p.value)
tabla_p <- data.frame(Muestra = muestras, "Valor p" = valores_p)

library(gt)
tabla_p %>%
  gt()
Muestra Valor.p
n=5 1.839000e-10
n=10 5.010504e-07
n=20 4.772843e-04
n=30 4.622600e-03
n=50 8.829914e-03
n=60 5.587520e-02
n=100 4.836092e-02
n=200 7.928454e-02
n=500 6.538411e-02

El primer histograma evidencia que existe una asimetría negativa al estar las estimaciones concentradas a la derecha, la cual puede observarse en el valor de simetría calculada que da negativo. A pesar de haber asimetría en la distribución, el sesgo expone que el valor de \(\hat{p}\) es muy cercano a \(p\) y donde el coeficiente de variación es muy cercano a cero.

Con los diferentes tamaños de muestra, los histogramas y qqnorms demuestran que se alcanza un comportamiento de Normalidad a partir de \(n>200\) y donde el test de SW lo confirma cuando se deja de rechazar \(H_{0}\) a partir de este tamaño de muestra.

4. Conclusiones

Con los casos expuestos se infiere que la convergencia a la Normalidad y cumplimiento del Teorema del Límite Central viene dado tanto por el \(p\) o parámetro de la población como también de los tamaños de muestra a escoger. Se observó que al tener desde un principio un \(p=0.5\) llegar a la convergencia fue mucho más rápido, pues se lograba simetría desde un principio y por tanto con un tamaño de muestra mayor a 30, \(\hat{p}\) se acercaba mucho más al valor real \(p\). En cambio, esto mismo no sucedió cuando se tenían \(p\) extremos, como fueron los casos de \(p=0.1\) y \(p=0.9\), donde la simetría no se cumplía en un inicio y llegar a una convergencia necesitaría mayor tamaño de muestra.