Actividad 2- Informe 3 - Métodos y simulación Estadística
Cesar Fuentes - Daniel Sanchez
2024-03-18
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, desarrollamos el paso a paso para probar mediante simulación este teorema.
Paso A
Se creó la población de tamaño \(1000\), la cual tiene un porcentaje de plantas enfermas del 50%. Para este caso las plantas enfermas se representan con el número 1.
set.seed(100)
Muestra = list()
p = 0
N= 1000
Plantas_sanas = rep(0, N*0.5)
Plantas_enfermas = rep(1, N*0.5)
cons = c(Plantas_sanas, Plantas_enfermas)
Poblacion = sample(cons)
head(Poblacion,20)## [1] 1 1 0 1 1 1 1 0 1 1 1 1 0 1 0 0 0 1 0 1Paso B
Se desarrolló una función que selecciona muestras aleatorias de la población establecida previamente, según un tamaño especificado. Posteriormente, se calcula el estimador de la proporción muestral \(\hat p\) para el tamaño de muestra especificado, n.
EstimadorP <- function(n){
Muestra = sample(Poblacion, n)
p = sum(Muestra)/n
cat("El estimador de la proporción muestral p̂ es:", p, "para un tamaño de muestra", n )
}Ejemplo:
EstimadorP(45)## El estimador de la proporción muestral p̂ es: 0.4888889 para un tamaño de muestra 45Paso C y D
Se realiza la estimación para 500 repeticiones con los siguientes tamaños de muestra \(n = 5, 10, 15, 20, 30, 50, 60, 100, 200\) y \(500\). Para cada uno de estos tamaños se analizará su simetría, variabilidad y su normalidad.
TLC <- function(n){for (i in 1:500) {
p[i] = sum(sample(Poblacion, n))/n
}
P = data.frame(p)
return(P)
}Simetría y Sesgo
Se muestra a continuación una serie de gráficas de densidad que ilustran la distribución de la muestra a diferentes niveles de n
# Simetría
# Valores de n
valores_n <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)
# Lista para almacenar los resultados de la prueba de Shapiro-Wilk
resultados <- list()
# Crear una ventana gráfica con múltiples paneles
par(mfrow = c(2, 5))
# Iterar sobre los valores de n
for (n in valores_n) {
df <- TLC(n)
# Crear el gráfico de densidad
den<-(density(df$p))
plot(den, frame = FALSE, col = "blue",main = paste("n=",n))
abline(v=0.5)
}
Inicialmente podemos ver en la gráfica anterior que a mayor cantidad de repeticiones la simetría de los datos mejora, viendose un resultado mas notorio a partir de \(n=30\), validandose lo dicho por algunos autores que indican que el umbral es bueno a partir de \(n>30\).
Ahora vemos el comportamiento a diferentes niveles de n en cuanto al sesgo absoluto:
# Sesgo
resultadosses <- list()
for (n in valores_n) {
ses <- abs(mean(TLC(n)$p) - 0.5)
# Almacenar el resultado en la lista
resultadosses[[as.character(n)]] <- c(Sesgo = ses)
}
# Convertir la lista a un dataframe
resultados_df_ses <- do.call(rbind, resultadosses)
resultados_df_ses## Sesgo
## 5 0.0168000000
## 10 0.0012000000
## 15 0.0036000000
## 20 0.0043000000
## 30 0.0007333333
## 50 0.0021600000
## 60 0.0027666667
## 100 0.0005600000
## 200 0.0006700000
## 500 0.0009760000Ahora evaluando el sesgo, vemos que el sesgo absoluto va disminuyendo, acercadose a cero a medida que aumenta el número de repeticiones, presentandose en algunos casos un leve aumento, pero sigue su tendencia al cero.
Por ultimo, se muestran gráficas de boxplot que representan la distribución de los datos a distintos niveles de n. Cada boxplot refleja la simetría de los datos y cualquier indicio de sesgo para el tamaño de muestra correspondiente.
# Crear una ventana gráfica con múltiples paneles
par(mfrow = c(2, 5))
# Iterar sobre los valores de n
for (n in valores_n) {
df <- TLC(n)
# Crear el gráfico de densidad
boxplot(df$p, frame = FALSE,main = paste("n=",n))
abline(h=0.5)
}
Analizando el gráfico de cajas, nos confirma que los datos son simetrico e insesgado a partir de n=30 ya que las cajas se mantienen de un tamaño constante, y su mediana es cercana al valor del parámetro.
Variabilidad
Evaluamos la variabilidad mediante el coeficiente de variación:
# Variabilidad
resultadoscv <- list()
for (n in valores_n) {
cv <- (sd(TLC(n)$p) / mean(TLC(n)$p)) * 100
# Almacenar el resultado en la lista
resultadoscv[[as.character(n)]] <- c(Coeficiente_Variacion = cv)
}
# Convertir la lista a un dataframe
resultados_df_cv <- do.call(rbind, resultadoscv)
resultados_df_cv ## Coeficiente_Variacion
## 5 41.747151
## 10 33.460781
## 15 26.091472
## 20 22.413226
## 30 18.699326
## 50 13.645488
## 60 12.703710
## 100 9.602759
## 200 6.586460
## 500 3.263849Los resultados demuestran que a mayor número de repeticiones la variabilidad disminuye, siguiendo la lógica del teorema del límite central.
Normalidad
Evaluamos la normalidad con la prueba Shapiro-Wilk, dónde la hipótesis nula es que los datos siguen una distribución normal.
# Iterar sobre los valores de n
for (n in valores_n) {
df <- TLC(n)
shapiro_test <- shapiro.test(df$p)
# Almacenar el resultado en la lista
resultados[[as.character(n)]] <- c(p_valor = round(shapiro_test$p.value,5))
}
# Convertir la lista a un dataframe
resultados_df <- do.call(rbind, resultados)
resultados_df <- data.frame(p_valor = resultados_df)
resultados_df$Contraste <- ifelse(resultados_df$p_valor < 0.05,"Rechazo Ho","No rechazo Ho")
resultados_df$n <- as.numeric(rownames(resultados_df))
print(resultados_df)## p_valor Contraste n
## 5 0.00000 Rechazo Ho 5
## 10 0.00000 Rechazo Ho 10
## 15 0.00000 Rechazo Ho 15
## 20 0.00000 Rechazo Ho 20
## 30 0.00013 Rechazo Ho 30
## 50 0.00450 Rechazo Ho 50
## 60 0.00250 Rechazo Ho 60
## 100 0.04410 Rechazo Ho 100
## 200 0.07043 No rechazo Ho 200
## 500 0.10606 No rechazo Ho 500Según los resultados del test, se evidencia que a medida que incrementa el tamaño de muestra los datos tienden a seguir una distribución normal, en este caso a partir de las 200 repeticiones No se rechaza la Hipótesis nula, usando un \(\alpha = 0.05\), lo que sugiere que los datos siguen una distribución normal.
También visualizamos gráficamente mediante qqnorm:
### qqnorm
# Crear una ventana gráfica con múltiples paneles
par(mfrow = c(2, 5))
# Iterar sobre los valores de n
for (n in valores_n) {
df <- TLC(n)
# Crear un gráfico Q-Q normal para la variable p
qqnorm(df$p, main = paste("Q-Q Plot n =", n))
qqline (df$p, col = "RED")
}
Gráficamente, se observa una tendencia de los datos a seguir una distribución normal. Después de 30 repeticiones, los datos están más cercanos a la línea diagonal en el gráfico Q-Q, indicando que los cuantiles de los datos se ajustan bien a una distribución normal. Sin embargo, los resultados del test confirman que es después de 200 repeticiones cuando los datos se distribuyen de manera más cercana a una distribución normal, para esta simulación.
Paso E
Se 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.
Iniciamos cambiando la población para la proporción de enfermos de 10% y 90%:
Vemos un ejemplo de la población con un porcentaje del 10% de enfermos.
set.seed(100)
Muestra = list()
p = 0
N= 1000
Plantas_sanas2 = rep(0, N*0.9)
Plantas_enfermas2 = rep(1, N*0.1)
cons = c(Plantas_sanas2, Plantas_enfermas2)
Poblacion2 = sample(cons)
head(Poblacion2,100)## [1] 0 0 0 0 1 0 1 0 1 0 0 0 0 1 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
## [38] 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
## [75] 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0Vemos un ejemplo de la población con un porcentaje del 90% de enfermos.
Plantas_sanas3 = rep(0, N*0.1)
Plantas_enfermas3 = rep(1, N*0.9)
cons = c(Plantas_sanas3, Plantas_enfermas3)
Poblacion3 = sample(cons)
head(Poblacion3,100)## [1] 0 1 1 1 0 1 1 1 1 1 1 1 1 1 0 1 1 0 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1
## [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [75] 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1Creamos las funciones para cada población que calculan el p 500 veces con una muestra de tamaño n.
TLC10ENF <- function(n){for (i in 1:500) {
p[i] = sum(sample(Poblacion2, n))/n
}
P = data.frame(p)
return(P)
}
TLC90ENF <- function(n){for (i in 1:500) {
p[i] = sum(sample(Poblacion3, n))/n
}
P = data.frame(p)
return(P)
}Simetría y Sesgo
Visualizamos la simetria para la población con 10% de enfermos:
# Simetría
# Valores de n
valores_n <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)
# Lista para almacenar los resultados de la prueba de Shapiro-Wilk
resultados <- list()
# Crear una ventana gráfica con múltiples paneles
par(mfrow = c(2, 5))
# Iterar sobre los valores de n
for (n in valores_n) {
df <- TLC10ENF(n)
# Crear el gráfico de densidad
den<-(density(df$p))
plot(den, frame = FALSE, col = "blue",main = paste("n=",n))
abline(v=0.1)
}
Con la gráfica de densidad podemos ver una simetria a partir de \(n=30\) ya que las colas de los gráficos no se ven mas largas en alguna dirección.
Visualizamos la simetria para la población con 90% de enfermos:
# Crear una ventana gráfica con múltiples paneles
par(mfrow = c(2, 5))
for (n in valores_n) {
df <- TLC90ENF(n)
# Crear el gráfico de densidad
den<-(density(df$p))
plot(den, frame = FALSE, col = "blue",main = paste("n=",n))
abline(v=0.9)
}
Por su parte, para el porcentaje de enfermos igual al 90% vemos un patrón similar al de la gráfica anterior, presentando una simetría a partir de \(n=30\) donde las colas de los gráficos no se ven más largas en alguna dirección.
Ahora, evaluamos el sesgo para 10% enfermos:
# Sesgo
resultadosses10 <- list()
for (n in valores_n) {
ses <- abs(mean(TLC10ENF(n)$p) - 0.1)
# Almacenar el resultado en la lista
resultadosses10[[as.character(n)]] <- c(Sesgo = ses)
}
# Convertir la lista a un dataframe
resultados_df_ses_10_enf <- do.call(rbind, resultadosses10)
resultados_df_ses_10_enf## Sesgo
## 5 0.0056000000
## 10 0.0006000000
## 15 0.0009333333
## 20 0.0013000000
## 30 0.0029333333
## 50 0.0007600000
## 60 0.0009000000
## 100 0.0005200000
## 200 0.0003100000
## 500 0.0000400000# Crear una ventana gráfica con múltiples paneles
par(mfrow = c(2, 5))
# Iterar sobre los valores de n
for (n in valores_n) {
df <- TLC10ENF(n)
# Crear el gráfico de densidad
boxplot(df$p, frame = FALSE,main = paste("n=",n))
abline(h=0.1)
}
Con el resultado del cálculo del sesgo para el porcentaje de plantas enfermas del 10%, vemos que a partir de 30 repeticiones el sesgo empieza a disminuir su valor tendiendo a cero. Con el gráfico de caja podemos validar esta afirmación ya que se observa su simetría en las cajas y que la mediana esta cercana al parámetro.
Ahora, vemos el sesgo con 90% enfermos:
# Sesgo
resultadosses90 <- list()
for (n in valores_n) {
ses <- abs(mean(TLC90ENF(n)$p) - 0.9)
# Almacenar el resultado en la lista
resultadosses90[[as.character(n)]] <- c(Sesgo = ses)
}
# Convertir la lista a un dataframe
resultados_df_ses_90_enf <- do.call(rbind, resultadosses90)
resultados_df_ses_90_enf## Sesgo
## 5 0.000800000
## 10 0.000600000
## 15 0.005466667
## 20 0.000800000
## 30 0.002866667
## 50 0.000560000
## 60 0.002266667
## 100 0.000960000
## 200 0.000820000
## 500 0.000384000# Crear una ventana gráfica con múltiples paneles
par(mfrow = c(2,5 ))
# Iterar sobre los valores de n
for (n in valores_n) {
df <- TLC90ENF(n)
# Crear el gráfico de densidad
boxplot(df$p, frame = FALSE,main = paste("n=",n))
abline(h=0.9)
}
Por su parte, con el resultado del cálculo del sesgo para el porcentaje del plantas enfermas del 90%, vemos un patrón similar al del 10% de plantas enfermas, dónde a partir de 30 repeticiones el sesgo empieza a disminuir su valor tendiendo a cero y en el gráfico de cajas validamos esta afirmación ya que se observa su simetría en las cajas y la mediana cercana al parámetro.
Variabilidad
Evaluamos la variabilidad mediante el coeficiente de variación:
Inicialmente para la población de 10% enfermos:
# Variabilidad
resultadoscv10 <- list()
for (n in valores_n) {
cv <- (sd(TLC10ENF(n)$p) / mean(TLC10ENF(n)$p)) * 100
# Almacenar el resultado en la lista
resultadoscv10[[as.character(n)]] <- c(Coeficiente_Variacion = cv)
}
# Convertir la lista a un dataframe
resultados_df_cv_10enf <- do.call(rbind, resultadoscv10)
resultados_df_cv_10enf## Coeficiente_Variacion
## 5 128.63836
## 10 95.72262
## 15 72.49922
## 20 62.17642
## 30 52.11195
## 50 38.84524
## 60 38.16080
## 100 27.41329
## 200 18.51690
## 500 10.18540Para el 10% de plantas enfermas vemos que tiene una variabilidad alta en las primeras repeticiones, pero va disminuyendo a medida que aumenta el número de repeticiones, llegando con 500 repeticiones hasta el 10% en el coeficiente de variación.
Ahora para la población de 90% enfermos:
# Variabilidad
resultadoscv90 <- list()
for (n in valores_n) {
cv <- (sd(TLC90ENF(n)$p) / mean(TLC90ENF(n)$p)) * 100
# Almacenar el resultado en la lista
resultadoscv90[[as.character(n)]] <- c(Coeficiente_Variacion = cv)
}
# Convertir la lista a un dataframe
resultados_df_cv_90enf <- do.call(rbind, resultadoscv90)
resultados_df_cv_90enf## Coeficiente_Variacion
## 5 14.839166
## 10 11.183830
## 15 8.247376
## 20 7.429268
## 30 5.803508
## 50 4.747461
## 60 4.146956
## 100 3.244822
## 200 2.097762
## 500 1.077853Para el 90% de plantas enfermas vemos que la variabilidad va disminuyendo a medida que aumenta el número de repeticiones, llegando con 500 repeticiones hasta el 1% en el coeficiente de variación.
Si comparamos la variabilidad entre el 10% de plantas enfermas y el 90% vemos que esta segunda siempre tiene un porcentaje de variabilidad menor al primero.
Normalidad
Evaluamos la normalidad con la prueba Shapiro.
Inicialmente para la población de 10% enfermos:
# Iterar sobre los valores de n
resultadosnorm10 <-list()
for (n in valores_n) {
df <- TLC10ENF(n)
shapiro_test <- shapiro.test(df$p)
# Almacenar el resultado en la lista
resultadosnorm10[[as.character(n)]] <- c(p_valor = round(shapiro_test$p.value,5))
}
# Convertir la lista a un dataframe
resultados_df_norm10 <- do.call(rbind, resultadosnorm10)
resultados_df_norm10 <- data.frame(p_valor = resultados_df_norm10)
resultados_df_norm10$contraste <- ifelse(resultados_df_norm10$p_valor < 0.05,"Rechazo Ho","No rechazo Ho")
resultados_df_norm10$n <- as.numeric(rownames(resultados_df_norm10))
print(resultados_df_norm10)## p_valor contraste n
## 5 0.00000 Rechazo Ho 5
## 10 0.00000 Rechazo Ho 10
## 15 0.00000 Rechazo Ho 15
## 20 0.00000 Rechazo Ho 20
## 30 0.00000 Rechazo Ho 30
## 50 0.00000 Rechazo Ho 50
## 60 0.00000 Rechazo Ho 60
## 100 0.00001 Rechazo Ho 100
## 200 0.00045 Rechazo Ho 200
## 500 0.00522 Rechazo Ho 500Según los resultados del test, se evidencia que a medida que incrementa el número de repeticiones el valor p va incrementando, sin embargo, no pasa el umbral de rechazo de la hipótesis nula, por lo tanto, para el número de repeticiones de 500 ó menos, no sigue una distribución normal en esta población.
Ahora para la población de 90% enfermos:
# Iterar sobre los valores de n
resultadosnorm90<-list()
for (n in valores_n) {
df <- TLC90ENF(n)
shapiro_test <- shapiro.test(df$p)
# Almacenar el resultado en la lista
resultadosnorm90[[as.character(n)]] <- c(p_valor = round(shapiro_test$p.value,5))
}
# Convertir la lista a un dataframe
resultados_df_norm90 <- do.call(rbind, resultadosnorm90)
resultados_df_norm90 <- data.frame(p_valor = resultados_df_norm90)
resultados_df_norm90$contraste <- ifelse(resultados_df_norm90$p_valor < 0.05,"Rechazo Ho","No rechazo Ho")
resultados_df_norm90$n <- as.numeric(rownames(resultados_df_norm90))
print(resultados_df_norm90)## p_valor contraste n
## 5 0.00000 Rechazo Ho 5
## 10 0.00000 Rechazo Ho 10
## 15 0.00000 Rechazo Ho 15
## 20 0.00000 Rechazo Ho 20
## 30 0.00000 Rechazo Ho 30
## 50 0.00000 Rechazo Ho 50
## 60 0.00000 Rechazo Ho 60
## 100 0.00000 Rechazo Ho 100
## 200 0.00323 Rechazo Ho 200
## 500 0.05870 No rechazo Ho 500Por su parte, para la población con 90% enfermos, se evidencia que a medida que incrementa el número de repeticiones, el valor p va incrementando también, sin embargo, no pasa el umbral de rechazo de la hipótesis nula hasta las 500 repeticiones, dónde ya sigue una distribución normal.
También visualizamos gráficamente mediante los Q-Q Plot, que con la condición del 10% y 90% solo se acercan a una distribución normal cuando se realizan 500 repeticiones.
Inicialmente para 10% de enfermos:
### qqnorm
# Crear una ventana gráfica con múltiples paneles
par(mfrow = c(2, 5))
# Iterar sobre los valores de n
for (n in valores_n) {
df <- TLC10ENF(n)
# Crear un gráfico Q-Q normal para la variable p
qqnorm(df$p, main = paste("Q-Q Plot n =", n))
qqline (df$p, col = "RED")
}
Posteriormente para 90% de enfermos:
### qqnorm
# Crear una ventana gráfica con múltiples paneles
par(mfrow = c(2, 5))
# Iterar sobre los valores de n
for (n in valores_n) {
df <- TLC90ENF(n)
# Crear un gráfico Q-Q normal para la variable p
qqnorm(df$p, main = paste("Q-Q Plot n =", n))
qqline (df$p, col = "RED")
}
Conclusión final
Evidenciamos que cuando se tiene un parámetro de distribución \(p\) cercano o igual al \(0.5\), los datos tienden a seguir una distribución normal a partir de las 30 repeticiones. Sin embargo, cuando el parámetro de estimación se lleva a los extremos, se requieren más registros para que los datos sigan una distribución normal.
Algo que si se evidencia en este ejemplo es que los datos tienen simetría e insesgadez a medida que se aumenta el número de repeticiones, indiferente al valor del parámetro de distribución \(p\).
En todo el ejercicio evidenciamos que los gráficos son una ayuda visual en la distribución de datos, sin embargo, se deben ejecutar los test para confirmar la distribución de los datos con mayor certeza.