Problema 3
Teorema del Límite Central
Diego Felipe Feria Gómez
Maestría en
Ciencia de Datos
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:
a.
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%.
n=1000
x = rbinom(n,1,0.5)
table(x)/nx
0 1
0.507 0.493 b.
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.
muestra = function(n){
m = rbinom(n,1,0.5)
p = sum(m)/n
return(p)
}
muestra(20)[1] 0.4c.
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.
replica = rep(20,500)
x20=sapply(replica,muestra)
hist(x20)
shapiro.test(x20)
Shapiro-Wilk normality test
data: x20
W = 0.9781, p-value = 7.916e-07para un n de 20, los resultados son simétricos cuando se observan gráficamente , ya que obedecen a una distribución normal. Pero al aplicar el test de normalidad de Shapiro-Wilk, el p-value es muy bajo (<0.05) lo que representa no normalidad, este valor mejoraría si se toma una muestra más grande.
d.
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
n=500
m=1000*n
X= matrix(rbinom(m,1,0.5),ncol=n)
X5=X[ ,1:5] # n=5
X10=X[ ,1:10] # n=10
X15=X[ ,1:15] # n=15
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
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
Mx15=apply(X15,1,mean) # medias de muestras de tamaño n=15
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=500par(cex=1, cex.axis=1, cex.lab=1, cex.main=1, cex.sub=1, mfrow=c(4,3), mai = c(.5, .5, .5, .5))
hist(Mx5, main = "n=5", freq=FALSE, xlab = " ")
hist(Mx10, main ="n=10", freq=FALSE, xlab = " ")
hist(Mx15, main ="n=15", freq=FALSE, xlab = " ")
hist(Mx20, main = "n=20",freq=FALSE, xlab = " ")
hist(Mx30, main = "n=30",freq=FALSE, xlab = " ")
hist(Mx50, main = "n=50",freq=FALSE, xlab = " ")
hist(Mx60, main = "n=60",freq=FALSE, xlab = " ")
hist(Mx100, main = "n=100", freq=FALSE, xlab = " ")
hist(Mx200, main = "n=200", freq=FALSE, xlab = " ")
hist(Mx500, main = "n=500", freq=FALSE, xlab = " ") 
par(cex=1, cex.axis=1, cex.lab=1, cex.main=1, cex.sub=1, mfrow=c(4,3), mai = c(.5, .5, .5, .5))
qqnorm(Mx5,main="n=5", xlab = " ") ; qqline(Mx5, col="blue")
qqnorm(Mx10,main="n=10", xlab = " ") ; qqline(Mx10, col="blue")
qqnorm(Mx15,main="n=15", xlab = " ") ; qqline(Mx15, col="blue")
qqnorm(Mx20,main="n=20", xlab = " ") ; qqline(Mx20, col="blue")
qqnorm(Mx30,main="n=30", xlab = " ") ; qqline(Mx30, col="blue")
qqnorm(Mx50,main="n=50", xlab = " ") ; qqline(Mx50, col="blue")
qqnorm(Mx60,main="n=60", xlab = " ") ; qqline(Mx60, col="blue")
qqnorm(Mx100,main="n=100", xlab = " ") ; qqline(Mx100, col="blue")
qqnorm(Mx200,main="n=200", xlab = " ") ; qqline(Mx200, col="blue")
qqnorm(Mx500,main="n=500", xlab = " ") ; qqline(Mx500, col="blue")
Gráficamente se puede observar una convergencia de la media muestral
a una distribución normal, a medida que el n aumenta, siendo más
ajustada para tamaños mayores a 200.
Lote con 50% de plantas enfermas
P.value5=shapiro.test(Mx5)$p.value
P.value10=shapiro.test(Mx10)$p.value
P.value15=shapiro.test(Mx15)$p.value
P.value20=shapiro.test(Mx20)$p.value
P.value30=shapiro.test(Mx30)$p.value
P.value50=shapiro.test(Mx50)$p.value
P.value60=shapiro.test(Mx60)$p.value
P.value100=shapiro.test(Mx100)$p.value
P.value200=shapiro.test(Mx200)$p.value
P.value500=shapiro.test(Mx500)$p.value
df_p_values <-data.frame(Muestra = c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500),
P_value = c(P.value5, P.value10, P.value15, P.value20, P.value30, P.value50, P.value60, P.value100, P.value200, P.value500))
df_p_values$P_value <- format(df_p_values$P_value, scientific = FALSE)
print(df_p_values) Muestra P_value
1 5 0.000000000000000000005028308
2 10 0.000000000000041153840841872
3 15 0.000000000031609811230227107
4 20 0.000000001991086855555897371
5 30 0.000000180161336932262192137
6 50 0.000049820573857312825220844
7 60 0.000265209829482451720731456
8 100 0.000769491152106604623869823
9 200 0.044365492673859585615048218
10 500 0.075051562685342715353087328Los tamaños de muestra que presentan mayores P valor y por ende
tienden a la normalidad, por contraste de hipótesis, son los de 200
siendo muy cercado y 500.
e.
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.
Lote con 10% de plantas enfermas
n=500
m=1000*n
X= matrix(rbinom(m,1,0.1),ncol=n)
X5=X[ ,1:5] # n=5
X10=X[ ,1:10] # n=10
X15=X[ ,1:15] # n=15
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
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
Mx15=apply(X15,1,mean) # medias de muestras de tamaño n=15
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
P.value5=shapiro.test(Mx5)$p.value
P.value10=shapiro.test(Mx10)$p.value
P.value15=shapiro.test(Mx15)$p.value
P.value20=shapiro.test(Mx20)$p.value
P.value30=shapiro.test(Mx30)$p.value
P.value50=shapiro.test(Mx50)$p.value
P.value60=shapiro.test(Mx60)$p.value
P.value100=shapiro.test(Mx100)$p.value
P.value200=shapiro.test(Mx200)$p.value
P.value500=shapiro.test(Mx500)$p.value
df_p_values <-data.frame(Muestra = c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500),
P_value = c(P.value5, P.value10, P.value15, P.value20, P.value30, P.value50, P.value60, P.value100, P.value200, P.value500))
df_p_values$P_value <- format(df_p_values$P_value, scientific = FALSE)
print(df_p_values) Muestra P_value
1 5 0.000000000000000000000000000000000000003106686
2 10 0.000000000000000000000000000001664407326358767
3 15 0.000000000000000000000000133382751823680738745
4 20 0.000000000000000000000214643624592265301761172
5 30 0.000000000000000000300985791887729155474201503
6 50 0.000000000000690844158941540789674665146158361
7 60 0.000000000004540462388608267315162625443747402
8 100 0.000000002835416876917483193498714655333969858
9 200 0.000024903368841627533912530534365714629529975
10 500 0.010467720689853684395975719212401600088924170Lote con 90% de plantas enfermas
n=500
m=1000*n
X= matrix(rbinom(m,1,0.9),ncol=n)
X5=X[ ,1:5] # n=5
X10=X[ ,1:10] # n=10
X15=X[ ,1:15] # n=15
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
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
Mx15=apply(X15,1,mean) # medias de muestras de tamaño n=15
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
P.value5=shapiro.test(Mx5)$p.value
P.value10=shapiro.test(Mx10)$p.value
P.value15=shapiro.test(Mx15)$p.value
P.value20=shapiro.test(Mx20)$p.value
P.value30=shapiro.test(Mx30)$p.value
P.value50=shapiro.test(Mx50)$p.value
P.value60=shapiro.test(Mx60)$p.value
P.value100=shapiro.test(Mx100)$p.value
P.value200=shapiro.test(Mx200)$p.value
P.value500=shapiro.test(Mx500)$p.value
df_p_values <-data.frame(Muestra = c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500),
P_value = c(P.value5, P.value10, P.value15, P.value20, P.value30, P.value50, P.value60, P.value100, P.value200, P.value500))
df_p_values$P_value <- format(df_p_values$P_value, scientific = FALSE)
print(df_p_values) Muestra P_value
1 5 0.000000000000000000000000000000000000003693462
2 10 0.000000000000000000000000000000721529742010773
3 15 0.000000000000000000000000222352322017886014490
4 20 0.000000000000000000000416533009946646279753751
5 30 0.000000000000000002729489304791072426829498521
6 50 0.000000000000125737586156119481159074147136323
7 60 0.000000000001193638972524367136973655334486466
8 100 0.000000182554042696706438179965470425258899922
9 200 0.000002096048855749475723275407346157805932307
10 500 0.020283231522570110971770418473170138895511627Al trabajar con muestras desbalanceadas (10% y 90%) que presentan
sesgos, es más difícil alcanzar la normalidad, por ello se necesitaría
una muestra aún más grande.