Métodos y Simulación Estadística - Actividad 2
Elaborado por: Carlos Manuel Arcos Ramírez y Harvey Bastidas Caicedo
18 de marzo de 2024
Problema 3: Teorema del Límite Central
El Teorema Central del Límite establece que la distribución de la suma o la media de un gran número de variables aleatorias idénticamente distribuidas tiende a presentar una distribución normal, sin importar la distribución original de las variables. Es uno de los teoremas más importantes en estadística inferencial. Algunos autores afirman que la aproximación es aceptable a partir de un tamaño de muestra \(n>30\).
El siguiente es un procedimiento de simulación que pretende verificar el uso y comportamiento del Teorema del Límite Central.
- Realizaremos una simulación para generar una población de datos de
tamaño \(1000\), en donde se asumirá
que cada dato representa el estado (enferma o no enferma) de cada una de
las plantas dentro de un lote de mil plantas. Para la generación de los
datos utilizaremos la función
rep()que nos permite simular los estados booleanos de individuos enfermos o no enfermos. En esta primera iteración asumimos que el 50% del lote está enfermo.
poblacion_1 = rep(c(0,1),each=500)
kable(table(poblacion_1), "html", escape = FALSE, caption = "Población 1") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed","bordered"), full_width = FALSE)| poblacion_1 | Freq |
|---|---|
| 0 | 500 |
| 1 | 500 |
- Teniendo la población simulada, crearemos una función con la que se obtendrá una muestra aleatoria de la población y se calculará el estimador de la proporción \(\hat p\) de plantas enfermas para un tamaño de muestra \(n\).
muestreo <- function (datos_pob,n_muestra) {
return(sample(datos_pob,n_muestra,replace=TRUE))
}
n_test = 50
muestra = muestreo(poblacion_1,n_test)
kable(table(muestra), "html", escape = FALSE, caption = "Muestra aleatoria n=50") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed","bordered"), full_width = FALSE)| muestra | Freq |
|---|---|
| 0 | 28 |
| 1 | 22 |
phat_test <- function (datos_m) {
p = sum(datos_m)/length(datos_m)
p_table <- data.frame(Proporcion = p)
return(p_table)
}
phat_1 = phat_test(muestra)
kable(phat_1, "html", escape = FALSE, caption = "Proporción plantas enfermas") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed","bordered"), full_width = FALSE)| Proporcion |
|---|
| 0.44 |
- Ahora evaluaremos el escenario del muestreo y del cálculo de la proporción de plantas enfermas en 500 experimientos para \(n=5\) muestras de la población.
phat_grafico = function (pobla, n, m) {
# Contrucción de 500 experimientos a partir de muestras aleatorias
y = matrix(muestreo(pobla,n*m),ncol=n)
# Cálculo de la proporción
phat_funct = function(x){sum(x)/n}
# Evaluación de la proporción en cada experimiento
p_hat = apply(y,1,phat_funct)
# Gráficos
par(mfrow=c(1,3))
hist(p_hat,xlab=NULL,ylab=NULL,main=NULL,axes=F)
par(new=T)
plot(density(p_hat),col=2,lwd=3)
qqnorm(p_hat)
qqline(p_hat, col = "red")
boxplot(p_hat)
}
n = 5
m = 500
phat_grafico(poblacion_1,n,m)
Para \(n=5\) muestras en 500 experimentos los resultados no son muy sesgados y el promedio de la proporción estaría cercano al valor real de 0.5, como se observa en el gráfico Boxplot. También se puede observar que no hay mucha varianza, aunque el histograma no indica completa claridad en la simetría.
- Se realiza el mismo análisis para \(n=10,15,20,30,50,60,100,200,500\).
# Gráfico de histogramas
phat_hist = function (poblacion,n) {
y = matrix(muestreo(poblacion,n*500),ncol=n)
phat_funct = function(x){sum(x)/n}
p_hat = apply(y,1,phat_funct)
hist(p_hat, xlab=NULL,ylab=NULL,main=NULL,axes=F)
par(new=T)
plot(density(p_hat),col=2,lwd=3,main = paste("Histograma para n =",n))
}
# Gráfico de Q-Q Plot
phat_qqnorm = function (poblacion,n) {
y = matrix(muestreo(poblacion,n*500),ncol=n)
phat_funct = function(x){sum(x)/n}
p_hat = apply(y,1,phat_funct)
qqnorm(p_hat, main = paste("Q-Q Plot para n =",n))
qqline(p_hat, col = "red")
}
# Evaluación del test de Shapiro
phat_shaphiro = function(poblacion,n) {
y = matrix(muestreo(poblacion,n*500),ncol=n)
phat_funct = function(x){sum(x)/n}
p_hat = apply(y,1,phat_funct)
return(shapiro.test(p_hat)$p.value)
}
par(mfrow=c(3,3))
phat_hist(poblacion_1,10)
phat_hist(poblacion_1,15)
phat_hist(poblacion_1,20)
phat_hist(poblacion_1,30)
phat_hist(poblacion_1,50)
phat_hist(poblacion_1,60)
phat_hist(poblacion_1,100)
phat_hist(poblacion_1,200)
phat_hist(poblacion_1,500)
A medida que aumenta el número de muestras se percibe mejor la simetría en los histogramas. La media de la proporción de plantas enfermas estaría muy cercana al valor real de 0.5, como se observa en el historama para \(n=5\) muestras.
par(mfrow=c(3,3))
phat_qqnorm(poblacion_1,10)
phat_qqnorm(poblacion_1,15)
phat_qqnorm(poblacion_1,20)
phat_qqnorm(poblacion_1,30)
phat_qqnorm(poblacion_1,50)
phat_qqnorm(poblacion_1,60)
phat_qqnorm(poblacion_1,100)
phat_qqnorm(poblacion_1,200)
phat_qqnorm(poblacion_1,500)
El Q-Q Plot indica cómo los resultados asumen una distribución normal con mayor número de muestras.
p_values = c(phat_shaphiro(poblacion_1,10),
phat_shaphiro(poblacion_1,15),
phat_shaphiro(poblacion_1,20),
phat_shaphiro(poblacion_1,30),
phat_shaphiro(poblacion_1,50),
phat_shaphiro(poblacion_1,60),
phat_shaphiro(poblacion_1,100),
phat_shaphiro(poblacion_1,200),
phat_shaphiro(poblacion_1,500))
shapiro_table <- data.frame(N = c(10,15,20,30,50,60,100,200,500),p_value = p_values)
kable(shapiro_table, "html", escape = FALSE, caption = "Shapiro Test") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed","bordered"), full_width = FALSE)| N | p_value |
|---|---|
| 10 | 0.0000000 |
| 15 | 0.0000001 |
| 20 | 0.0000003 |
| 30 | 0.0000574 |
| 50 | 0.0015803 |
| 60 | 0.0011043 |
| 100 | 0.0040861 |
| 200 | 0.5320705 |
| 500 | 0.2983668 |
Del test de Shapiro-Wilk se concluye que los experimentos con \(n=200\) y \(n=500\) muestras son los que dan evidencia de que existe una distribución normal ya que en esos casos supera el valor de 0.05. Con estos resultados se pudo validar el funcionamiento del Teorema del Límite Central.
- Finalmente repetiremos toda la simulación para los casos en los que se tiene el 10% y el 90% de la población de plantas enfermas.
- Simulación para el caso en que se tiene el 10 % de la población de plantas enfermas.
# Creación de la población con el 10% de plantas enfermas
poblacion_2 = c(rep(1,100),rep(0,900))
kable(table(poblacion_2), "html", escape = FALSE, caption = "Población 2") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed","bordered"), full_width = FALSE)| poblacion_2 | Freq |
|---|---|
| 0 | 900 |
| 1 | 100 |
par(mfrow=c(3,3))
phat_hist(poblacion_2,10)
phat_hist(poblacion_2,15)
phat_hist(poblacion_2,20)
phat_hist(poblacion_2,30)
phat_hist(poblacion_2,50)
phat_hist(poblacion_2,60)
phat_hist(poblacion_2,100)
phat_hist(poblacion_2,200)
phat_hist(poblacion_2,500)
par(mfrow=c(3,3))
phat_qqnorm(poblacion_2,10)
phat_qqnorm(poblacion_2,15)
phat_qqnorm(poblacion_2,20)
phat_qqnorm(poblacion_2,30)
phat_qqnorm(poblacion_2,50)
phat_qqnorm(poblacion_2,60)
phat_qqnorm(poblacion_2,100)
phat_qqnorm(poblacion_2,200)
phat_qqnorm(poblacion_2,500)
p_values_2 = c(phat_shaphiro(poblacion_2,10),
phat_shaphiro(poblacion_2,15),
phat_shaphiro(poblacion_2,20),
phat_shaphiro(poblacion_2,30),
phat_shaphiro(poblacion_2,50),
phat_shaphiro(poblacion_2,60),
phat_shaphiro(poblacion_2,100),
phat_shaphiro(poblacion_2,200),
phat_shaphiro(poblacion_2,500))
shapiro_table_2 <- data.frame(N = c(10,15,20,30,50,60,100,200,500),p_value = p_values_2)
kable(shapiro_table_2, "html", escape = FALSE, caption = "Shapiro Test") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed","bordered"), full_width = FALSE)| N | p_value |
|---|---|
| 10 | 0.0000000 |
| 15 | 0.0000000 |
| 20 | 0.0000000 |
| 30 | 0.0000000 |
| 50 | 0.0000001 |
| 60 | 0.0000000 |
| 100 | 0.0000404 |
| 200 | 0.0000725 |
| 500 | 0.1300931 |
En esta simulación en los histogramas también se observa mejor simetría del estimador de la proporción a medida que aumenta el número de muestras. Para \(n=500\) se tiene simetría alrededor de 0.1, que es el valor de la proporción real. Con el Q-Q Plot podríamos asumir que la distribución es normal para muestras de \(n=200\) y \(n=500\). El test de Shapiro-Wilk indica que únicamente para \(n=500\) se tiene una distribución normal ya que solo en este caso el p-value supera el valor de \(0.05\).
- Simulación para el caso en que se tiene el 90 % de la población de plantas enfermas.
# Creación de la población con el 90% de plantas enfermas
poblacion_3 = c(rep(1,900),rep(0,100))
kable(table(poblacion_3), "html", escape = FALSE, caption = "Población 3") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed","bordered"), full_width = FALSE)| poblacion_3 | Freq |
|---|---|
| 0 | 100 |
| 1 | 900 |
par(mfrow=c(3,3))
phat_hist(poblacion_3,10)
phat_hist(poblacion_3,15)
phat_hist(poblacion_3,20)
phat_hist(poblacion_3,30)
phat_hist(poblacion_3,50)
phat_hist(poblacion_3,60)
phat_hist(poblacion_3,100)
phat_hist(poblacion_3,200)
phat_hist(poblacion_3,500)
par(mfrow=c(3,3))
phat_qqnorm(poblacion_3,10)
phat_qqnorm(poblacion_3,15)
phat_qqnorm(poblacion_3,20)
phat_qqnorm(poblacion_3,30)
phat_qqnorm(poblacion_3,50)
phat_qqnorm(poblacion_3,60)
phat_qqnorm(poblacion_3,100)
phat_qqnorm(poblacion_3,200)
phat_qqnorm(poblacion_3,500)
p_values_3 = c(phat_shaphiro(poblacion_3,10),
phat_shaphiro(poblacion_3,15),
phat_shaphiro(poblacion_3,20),
phat_shaphiro(poblacion_3,30),
phat_shaphiro(poblacion_3,50),
phat_shaphiro(poblacion_3,60),
phat_shaphiro(poblacion_3,100),
phat_shaphiro(poblacion_3,200),
phat_shaphiro(poblacion_3,500))
shapiro_table_3 <- data.frame(N = c(10,15,20,30,50,60,100,200,500),p_value = p_values_3)
kable(shapiro_table_3, "html", escape = FALSE, caption = "Shapiro Test") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed","bordered"), full_width = FALSE)| N | p_value |
|---|---|
| 10 | 0.0000000 |
| 15 | 0.0000000 |
| 20 | 0.0000000 |
| 30 | 0.0000000 |
| 50 | 0.0000000 |
| 60 | 0.0000007 |
| 100 | 0.0000084 |
| 200 | 0.0007041 |
| 500 | 0.3662163 |
Para esta simulación podríamos concluir lo mismo que para la anterior. En el histograma para \(n=500\) se observa que la media es aproximadamente 0.9, lo que demuestra las proporciones de la población estudiada. Al igual que el caso anterior, únicamente para \(n=500\) se tiene una distribución normal, como lo indica la tabla de resultados del test de Shapiro-Wilk.