1. Introducción

El teorema del límite central es un teorema fundamental de probabilidad y estadística. El teorema describe la distribución de la media de una muestra aleatoria proveniente de una población con varianza finita. Cuando el tamaño de la muestra es lo suficientemente grande, la distribución de las medias sigue aproximadamente una distribución normal. El teorema se aplica independientemente de la forma de la distribución de la población, mediante simulación en R se propone denostrar la validez del teorema.

2. Construcción del modelo

  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.
  • Calcule el estimador de la proporción muestral 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 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. Estimaciones

3.1 Muestra 50% plantas sanas 50% 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

Se generan numeros aleatorios medinate una distribución de probabilidad tipo binomial, se crea una función que permite extraer muestras de tamaño deseado donde el retorno de la función es el valor del estimador p^.

#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 0 1 0 0 0 1 0 0 0 1 1 0 0 1 1 0 1 0 1  tiene un p estimado = 0.45

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 p^ . Estas salidas quedarán guardadas en un data frame compuesto por número de la iteración y el 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
}

Generamos un histograma para observar la distribción, se calcula la simetría de los datos estimados, el sesgo y el coeficiente de variación.

#Paso 4

# Revisamos la distribución
library(ggplot2)
library(tidyverse)
library(hrbrthemes)

ggplot(df_estimacion, aes(x=p_estimados)) +
  geom_histogram(binwidth=0.05, fill="#d6a", 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)

# Calculamos 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.03281369
cat("El sesgo es igual a ", sesgo, "\n")
## El sesgo es igual a  -0.0066
cat("El coeficiente de variación es igual a ", cv)
## El coeficiente de variación es igual a  0.2291397

Observamos los resultados para 20 plantas, analizamos que las estimaciones de p^ presentan un comportamiento que tiende a la simetría, se deduce que 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.

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

Se repite 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 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

Verificamos la Normalidad de los datos, graficamos histogramas y qqnorms para cada uno de los tamaños de muestra, como también se calculará el valor p del test de Shapiro Wilk (SW) para verificar Normalidad en las estimaciones. En esta última, se asume como Hipótesis Nula (H0) 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="blue")
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="blue")
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="blue")
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="blue")
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="blue")
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="blue")
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="blue")
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="blue")
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="blue")
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.622913e-05
n=10 2.322422e-04
n=20 5.767686e-03
n=30 1.326152e-01
n=50 1.152676e-01
n=60 6.644464e-01
n=100 3.238166e-01
n=200 3.203707e-01
n=500 2.810447e-01

Revisando los resultados, 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 Shapiro wilk a partir del tamaño n>50 comienza a no rechazarse H0 al tener un valor p > 5%.

3.2 Muestra 90% plantas sanas 10% plantas enfermas

Se genera ahora un lote con 10% de plantas enfermas. A continuación se presentarán los histogramas, qqnorm y la 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 0 1 0 0 0 0 0 0 0 0 0  tiene un p estimado = 0.05
#---------------------------------------------------------------

# 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="grey", 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.5278716
cat("El sesgo es igual a ", sesgo, "\n")
## El sesgo es igual a  0.0087
cat("El coeficiente de variación es igual a ", cv)
## El coeficiente de variación es igual a  0.6265609
#---------------------------------------------------------------

#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="#3b8")
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="#3b8")
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="#3b8")
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="#3b8")
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="#3b8")
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="#3b8")
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="#3b8")
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="#3b8")
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="#3b8")
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 4.895364e-11
n=10 4.517152e-09
n=20 2.020044e-05
n=30 3.245669e-03
n=50 1.585995e-02
n=60 3.418930e-02
n=100 4.728603e-02
n=200 1.150513e-01
n=500 3.608452e-01

Analizando el primer histograma de los 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 H0 se da a partir del mencionado tamaño de muestra (valorp>0.05).

3.3 Muestra 10% plantas sanas 90% plantas enfermas

# 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  1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1  tiene un p estimado = 0.85
#---------------------------------------------------------------

# 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="grey", 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.4586448
cat("El sesgo es igual a ", sesgo, "\n")
## El sesgo es igual a  -0.0084
cat("El coeficiente de variación es igual a ", cv)
## El coeficiente de variación es igual a  0.07720928
#---------------------------------------------------------------

#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="#519")
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="#519")
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="#519")
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="#519")
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="#519")
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="#519")
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="#519")
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="#519")
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="#519")
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 2.525197e-12
n=10 1.026341e-08
n=20 3.833152e-05
n=30 4.578833e-03
n=50 1.847509e-02
n=60 2.804416e-02
n=100 2.359133e-01
n=200 6.119903e-01
n=500 5.984941e-01

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 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 H0 a partir de este tamaño de muestra.

3. 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 po 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, 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.