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.

A continuación se describen los siguientes pasos para su verificación:

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%.
Genere una función que permita:

Obtener una muestra aleatoria de la población y
Calcule el estimador de la proporción muestral pˆ para un tamaño de muestra dado n

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

Solución

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

Para esto se crea una función que permita construir una población n con una probabilidad de éxito p utilizando el método sample

# Establecer la semilla aleatoria para reproducibilidad
set.seed(241)

# Función para generar una población con un tamaño y probabilidad dadas
generar_poblacion <- function(tam, prob) {
  poblacion <- sample(c(rep(1, prob * tam), rep(0, (1 - prob) * tam)))
  return(poblacion)
}

Prueba con 10 elementos y probabilidad del 50%

pob <- generar_poblacion(10, 0.5)
pob

##  [1] 0 0 0 1 0 1 1 1 1 0

table(pob)

## pob
## 0 1 
## 5 5

Prueba con 10 elementos y probabilidad de éxito 10% y de fracaso 90%

pob2 <- generar_poblacion(10, 0.1)
pob2

##  [1] 0 0 0 1 0 0 0 0 0 0

table(pob2)

## pob2
## 0 1 
## 9 1

Genere una función que permita:

Obtener una muestra aleatoria de la población y
Calcule el estimador de la proporción muestral pˆ para un tamaño de muestra dado n

# Función para obtener una muestra aleatoria de la población y calcular el estimador de la proporción muestral
obtener_muestra_y_estimar <- function(poblacion, tam_muestra) {
  muestra <- sample(poblacion, tam_muestra)
  prop_muestral <- mean(muestra)
  return(prop_muestral)
}

#generar población de n = 1000 con el 50% de probabilidad de éxito
poblacion <- generar_poblacion(1000, 0.5)

estimador <- obtener_muestra_y_estimar(poblacion, 20)
print(paste("Estimador:", estimador))

## [1] "Estimador: 0.6"

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.

# Establecer la semilla aleatoria para reproducibilidad
set.seed(241)

# Función para generar una población con un tamaño y probabilidad dadas
generar_poblacion <- function(tam, prob) {
  poblacion <- sample(c(rep(1, prob * tam), rep(0, (1 - prob) * tam)))
  return(poblacion)
}

# Función para obtener una muestra aleatoria de la población y calcular el estimador de la proporción muestral
obtener_muestra_y_estimar <- function(poblacion, tam_muestra) {
  muestra <- sample(poblacion, tam_muestra)
  prop_muestral <- mean(muestra)
  return(prop_muestral)
}

# Simulación para analizar el comportamiento de los estimadores
simulacion <- function(tam_muestra, num_simulaciones, prob) {
  estimaciones <- replicate(num_simulaciones, {
    poblacion <- generar_poblacion(1000, prob)
    estimador <- obtener_muestra_y_estimar(poblacion, tam_muestra)
    return(estimador)
  })
  return(estimaciones)
}

# Ejecutar la simulación para n=20, 500 repeticiones y porcentaje de enfermos del 50%
resultados_20 <- simulacion(20, 500, 0.5)
hist(resultados_20, breaks = 20, main = "Histograma de estimadores (n=20, 50% de plantas enfermas)",
     xlab = "Estimador de proporción muestral", ylab = "Frecuencia")

summarytools::descr(resultados_20)

## Descriptive Statistics  
## resultados_20  
## N: 500  
## 
##                     resultados_20
## ----------------- ---------------
##              Mean            0.50
##           Std.Dev            0.11
##               Min            0.20
##                Q1            0.45
##            Median            0.50
##                Q3            0.60
##               Max            0.85
##               MAD            0.15
##               IQR            0.15
##                CV            0.21
##          Skewness           -0.03
##       SE.Skewness            0.11
##          Kurtosis           -0.18
##           N.Valid          500.00
##         Pct.Valid          100.00

Conclusión y análisis

Esta simulación entrega una estimación con sesgo bajo y distribución simétrica, siendo que solo está considerando una muetra de 20, con el fin de verificar la variabilidad se realizarán cálculos de desviación y se calculará el sesgo vs la media esperada de 0.5.

calcular_sesgo <- function(value, media){
  res <- media - value
  return(res)
}
media_poblacional <-0.5

sesgo <- sapply(resultados_20, calcular_sesgo, media = media_poblacional)
summary(sesgo)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.3500 -0.1000  0.0000 -0.0041  0.0500  0.3000

desviacion <- sd(sesgo)
print(paste("Desviación estándar del estimador:", desviacion))

## [1] "Desviación estándar del estimador: 0.108359002148845"

aquí se puede ver la poca variabilidad del estimador y el sesgo en promedio es del orden de milésimas, incluso con subestimación, esto sugiere que se trata de un buen estimador de la media poblacional, esto será más preciso a medida que se aumente el tamaño de muestra.

qqnorm(resultados_20, main="n=20") ; qqline(resultados_20, col="red")

Conclusión y análisis

A pesar que los puntos forman una escalera, comportamiento esperado ya que se trata de una distribución discreta y con una muestra pequeña, se puede apreciar la relación lineal entre los puntos, se sigue una tendencia ascendente y hay una aparente concentración de los mismos cerca a la diagonal. Esto da una idea de la variabilidad de los datos, el efecto del tamaño de la muestra y sirve como punto de referencia para validar el teorema a medida que se aumente el tamaño de la muestra y se espera que en efecto sigua una distribución normal con mayor fiabildad con cada iteración.

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

# Establecer la semilla aleatoria para reproducibilidad
set.seed(241)

N <- 1000
cant_simulaciones <- 500
prob_exito <- 0.5

generar_poblacion <- function(tam, prob) {
  poblacion <- sample(c(rep(1, prob * tam), rep(0, (1 - prob) * tam)))
  return(poblacion)
}

# Función para obtener una muestra aleatoria de la población y calcular el estimador de la proporción muestral
obtener_muestra_y_estimar <- function(poblacion, tam_muestra) {
  muestra <- sample(poblacion, tam_muestra)
  prop_muestral <- mean(muestra)
  return(prop_muestral)
}

# Simulación para analizar el comportamiento de los estimadores
simulacion <- function(tam_muestra, num_simulaciones, prob) {
  estimaciones <- replicate(num_simulaciones, {
    poblacion <- generar_poblacion(N, prob)
    estimador <- obtener_muestra_y_estimar(poblacion, tam_muestra)
    return(estimador)
  })
  return(estimaciones)
}

# Definir los tamaños de muestra
muestras <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)

# Crear una función para ejecutar la simulación para un tamaño de muestra dado
simulacion_tam <- function(tam_muestra, num_simulaciones, prob) {
  estimaciones <- simulacion(tam_muestra, num_simulaciones, prob)
  hist(estimaciones, main = paste("Histograma (n=", tam_muestra, ", p = 50% )", sep = ""), 
       xlab = "Estimador de proporción muestral", ylab = "Frecuencia")
  
    # Realizar el test de Shapiro-Wilk
  shapiro_test <- shapiro.test(estimaciones)
  cat("\nTest de Shapiro-Wilk para n =", tam_muestra, ":", "\n")
  print(shapiro_test)
  
  # Crear el gráfico de normalidad
  qqnorm(estimaciones, main = paste("Q-Q Plot (n=", tam_muestra, ", p = 50% )", sep = ""))
  qqline(estimaciones)
}

# Configurar el diseño de la cuadrícula de gráficos
par(mfrow = c(2, 2))

# Iterar sobre los diferentes tamaños de muestra y ejecutar la simulación
for (tam in muestras) {
  simulacion_tam(tam, cant_simulaciones, prob_exito)
}

## 
## Test de Shapiro-Wilk para n = 5 : 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.92491, p-value = 4.163e-15

## 
## Test de Shapiro-Wilk para n = 10 : 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.96148, p-value = 3.679e-10

## 
## Test de Shapiro-Wilk para n = 15 : 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.97431, p-value = 1.083e-07

## 
## Test de Shapiro-Wilk para n = 20 : 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.98076, p-value = 3.577e-06

## 
## Test de Shapiro-Wilk para n = 30 : 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.98651, p-value = 0.0001388

## 
## Test de Shapiro-Wilk para n = 50 : 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.98976, p-value = 0.001477

## 
## Test de Shapiro-Wilk para n = 60 : 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.99142, p-value = 0.005431

## 
## Test de Shapiro-Wilk para n = 100 : 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.99355, p-value = 0.03138

## 
## Test de Shapiro-Wilk para n = 200 : 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.98653, p-value = 0.0001402

## 
## Test de Shapiro-Wilk para n = 500 : 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.99662, p-value = 0.3778

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.

# Establecer la semilla aleatoria para reproducibilidad
set.seed(241)

N <- 1000
cant_simulaciones <- 500
prob_exito <- 0.1

generar_poblacion <- function(tam, prob) {
  poblacion <- sample(c(rep(1, prob * tam), rep(0, (1 - prob) * tam)))
  return(poblacion)
}

# Función para obtener una muestra aleatoria de la población y calcular el estimador de la proporción muestral
obtener_muestra_y_estimar <- function(poblacion, tam_muestra) {
  muestra <- sample(poblacion, tam_muestra)
  prop_muestral <- mean(muestra)
  return(prop_muestral)
}

# Simulación para analizar el comportamiento de los estimadores
simulacion <- function(tam_muestra, num_simulaciones, prob) {
  estimaciones <- replicate(num_simulaciones, {
    poblacion <- generar_poblacion(N, prob)
    estimador <- obtener_muestra_y_estimar(poblacion, tam_muestra)
    return(estimador)
  })
  return(estimaciones)
}

# Definir los tamaños de muestra
muestras <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)

# Crear una función para ejecutar la simulación para un tamaño de muestra dado
simulacion_tam <- function(tam_muestra, num_simulaciones, prob) {
  estimaciones <- simulacion(tam_muestra, num_simulaciones, prob)
  hist(estimaciones, main = paste("Histograma (n=", tam_muestra, ", p = 10% )", sep = ""), 
       xlab = "Estimador de proporción muestral", ylab = "Frecuencia")
  
  # Realizar el test de Shapiro-Wilk
  shapiro_test <- shapiro.test(estimaciones)
  cat("\nTest de Shapiro-Wilk para n =", tam_muestra, ":", "\n")
  print(shapiro_test)
  
  # Crear el gráfico de normalidad
  qqnorm(estimaciones, main = paste("Q-Q Plot (n=", tam_muestra, ", p = 10% )", sep = ""))
  qqline(estimaciones)
}

# Configurar el diseño de la cuadrícula de gráficos
par(mfrow = c(2, 2))

# Iterar sobre los diferentes tamaños de muestra y ejecutar la simulación
for (tam in muestras) {
  simulacion_tam(tam, cant_simulaciones, prob_exito)
}

## 
## Test de Shapiro-Wilk para n = 5 : 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.72623, p-value < 2.2e-16

## 
## Test de Shapiro-Wilk para n = 10 : 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.8483, p-value < 2.2e-16

## 
## Test de Shapiro-Wilk para n = 15 : 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.89186, p-value < 2.2e-16

## 
## Test de Shapiro-Wilk para n = 20 : 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.90759, p-value < 2.2e-16

## 
## Test de Shapiro-Wilk para n = 30 : 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.9511, p-value = 8.394e-12

## 
## Test de Shapiro-Wilk para n = 50 : 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.97076, p-value = 1.947e-08

## 
## Test de Shapiro-Wilk para n = 60 : 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.97128, p-value = 2.493e-08

## 
## Test de Shapiro-Wilk para n = 100 : 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.98002, p-value = 2.321e-06

## 
## Test de Shapiro-Wilk para n = 200 : 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.99153, p-value = 0.005902

## 
## Test de Shapiro-Wilk para n = 500 : 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.99275, p-value = 0.01601

Conclusiones y análisis

Con respecto a las simulaciones para una distribución con probabilidad 50%

Se comprueba el teorema del límite central incluso con tamaño de muestra pequeño, desde n=20 se ve una distribución mas o menos simétrica con evidencia de insesgadez y consistencia ya que a medida que se aumentaba el tamaño de la muestra, se puede comprobar mediante el gráfico de normalidad una evidencia que los datos siguen una distribución normal.
Con la prueba de Shapiro-Wilks se ve que la estadistica de prueba cada vez se acerca más a 1 indicando su compatibilidad con una distribución normal, de igual manera ocurre con el p-value que termina estando por encima de 0.05.

Con respecto a las simulaciones para una distribución con probabilidad 10%

Aunque en este caso la prueba de Shapiro-Wilks no es concluyente, se evidencia una tendencia a acercarse a una distribución normal. Por lo tanto también se considera que en este caso se cumple con el teorema de límite central teniendo en cuenta que los histogramas muestran mayor simetría y el gráfico de normalidad tiene buena concentración de valores y linealidad, excepto algunas desviaciones en las colas que puede ser producto de la aleatoriedad de las muestras.
Con respecto al test de Shapiro-Wilks, aunque el W cada vez se acerca más a 1, el p-value nunca estuvo por encima de 0.05, lo que nos permitiría concluir de forma definitiva que se trata de una distribución normal, por lo tanto se atribuye esto a la sensibilidad que tiene el test a desviaciones y tamaños de muestra, y se considera además que un solo test no es concluyente o que incluso podemos tener una muestra robusta que nos permite tener un estimador confiable una vez lo contrastamos con otros indicadores o gráficas, como el histograma y el diagrama de normalidad.

Teorema del límite central

Problema 3 - Metodos y Simulación - PUJC

Julián Alexis Correa

2024-03-16