Actividad 2. Problema 3

PONTIFICIA UNIVERSIDAD JAVERIANA DE CALI

Teorema del Límite Central

A continuación se describen los siguientes pasos para la verificación del teorema:

a. Simulación en una población de n=1000 (Lote), donde el porcentaje de plantas enfermas es de 50%.

set.seed(123)
n=1000
#x=rbinom(n, size=1, prob=0.5)
Lote=c(rep("enferma",500),rep("sana",500))
table(Lote)/n

Lote
enferma    sana 
    0.5     0.5 

b. Generación de una función para: - Obtener una muestra aleatoria de la población y - Calcular el estimador de la proporción muestral pˆpara un tamaño de muestra dado n.

#función
muestra=function(n) {
  m=rbinom(n, size=1, prob=0.5)
  p=sum(m)/n
  return(p)
}

cat("Estimador de la proporción muestral pˆ con n=30")

Estimador de la proporción muestral pˆ con n=30

set.seed(123)
muestra(30)

[1] 0.6333333

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

cat("Estimador de la proporción muestral pˆ con n=500")

Estimador de la proporción muestral pˆ con n=500

set.seed(123)
muestra(500)

[1] 0.47

replica=rep(500,500)
x500=sapply(replica, muestra)
hist(x500, main="Histograma n=500", col="#86AEDF", xlab="Probabilidad de muestras")
line = mean(x500)
abline (v=line, lwd = 4, lty = 2, col="darkblue")

qqnorm(x500)
qqline(x500, col="red")

cat("Prueba de normalidad")

Prueba de normalidad

shapiro.test(x500)

    Shapiro-Wilk normality test

data:  x500
W = 0.99489, p-value = 0.09606

cat("Análisis descriptivo")

Análisis descriptivo

summarytools::descr(x500)

Descriptive Statistics  
x500  
N: 500  

                      x500
----------------- --------
             Mean     0.50
          Std.Dev     0.02
              Min     0.43
               Q1     0.49
           Median     0.50
               Q3     0.51
              Max     0.57
              MAD     0.02
              IQR     0.03
               CV     0.04
         Skewness    -0.02
      SE.Skewness     0.11
         Kurtosis     0.27
          N.Valid   500.00
        Pct.Valid   100.00

Conclusión:

Se puede decir que la distribución es aproximadamente simétrica, teniendo en cuenta el valor del coeficiente de asimetria (-0.02) y que cumple con parámetros de normalidad como se observa en el qqplot y según la prueba de Shaipiro Wilk (p valor= 0.096).

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

set.seed(123)
n=500
m=1000*n
###
X=matrix(rbinom(m, size=1, prob=0.5), ncol=n) 
#generación de muestras
X5=X[,1:5]
X10=X[,1:10]
X15=X[,1:15]
X20=X[,1:20]
X30=X[,1:30]
X50=X[,1:50]
X60=X[,1:60]
X100=X[,1:100]
X200=X[,1:200]
X500=X[,1:500]
#Obención de las medias
Mx5=apply(X5,1,mean)
Mx10=apply(X10,1,mean)
Mx15=apply(X15,1,mean)
Mx20=apply(X20,1,mean)
Mx30=apply(X30,1,mean)
Mx50=apply(X50,1,mean)
Mx60=apply(X60,1,mean)
Mx100=apply(X100,1,mean)
Mx200=apply(X200,1,mean)
Mx500=apply(X500,1,mean)
#generación de densidad empirica
d5=density(Mx5)
d10=density(Mx10)
d15=density(Mx15)
d20=density(Mx20)
d30=density(Mx30)
d50=density(Mx50)
d60=density(Mx60)
d100=density(Mx100)
d200=density(Mx200)
d500=density(Mx500)
#plot
par(cex=0.8, cex.axis=1, cex.lab=.8, cex.main=.8, cex.sub=.8,  mfrow=c(2,2), mai = c(.3, .3, .3, .3))

hist(Mx5, freq = F, xlab="n=5", col = "violet")
hist(Mx10, freq = F, xlab="n=10", col = "violet")
hist(Mx15, freq = F, xlab="n=15", col = "violet")
hist(Mx20, freq = F, xlab="n=20", col = "violet")

hist(Mx30, freq = F, xlab="n=30", col = "violet")
hist(Mx50, freq = F, xlab="n=50", col = "violet")
hist(Mx60, freq = F, xlab="n=60", col = "violet")
hist(Mx100, freq = F, xlab="n=100", col = "violet")

hist(Mx200, freq = F, xlab="n=200", col = "violet")
hist(Mx500, freq = F, xlab="n=500", col = "violet")

#graficos de densidad
plot(d5, main="n=5", xlab="n=5")
plot(d10, main="n=10", xlab="n=10")

plot(d15, main="n=15", xlab="n=15")
plot(d20, main="n=20", xlab="n=20")
plot(d30, main="n=30", xlab="n=30")
plot(d50, main="n=50", xlab="n=50")

plot(d60, main="n=60", xlab="n=60")
plot(d100, main="n=100", xlab="n=100")
plot(d200, main="n=200", xlab="n=200")
plot(d500, main="n=500", xlab="n=500")

cat("Prueba de normalidad")

Prueba de normalidad

shapiro.test(Mx5)

    Shapiro-Wilk normality test

data:  Mx5
W = 0.92971, p-value < 2.2e-16

shapiro.test(Mx10)

    Shapiro-Wilk normality test

data:  Mx10
W = 0.96536, p-value = 1.106e-14

shapiro.test(Mx20)

    Shapiro-Wilk normality test

data:  Mx20
W = 0.98294, p-value = 2.003e-09

shapiro.test(Mx30)

    Shapiro-Wilk normality test

data:  Mx30
W = 0.98712, p-value = 1.067e-07

shapiro.test(Mx50)

    Shapiro-Wilk normality test

data:  Mx50
W = 0.99202, p-value = 3.128e-05

shapiro.test(Mx60)

    Shapiro-Wilk normality test

data:  Mx60
W = 0.99369, p-value = 0.0003123

shapiro.test(Mx100)

    Shapiro-Wilk normality test

data:  Mx100
W = 0.99477, p-value = 0.001568

shapiro.test(Mx200)

    Shapiro-Wilk normality test

data:  Mx200
W = 0.99723, p-value = 0.08372

shapiro.test(Mx500)

    Shapiro-Wilk normality test

data:  Mx500
W = 0.99785, p-value = 0.2253

# grafico de normalidad
par(cex=0.8, cex.axis=.8, cex.lab=.8, cex.main=.8, cex.sub=.5,  mfrow=c(2,3), mai = c(.3, .3, .3, .3))
qqnorm(Mx5, main="n=5");qqline(Mx5, col="red")
qqnorm(Mx10, main="n=10");qqline(Mx10, col="red")
qqnorm(Mx15, main="n=15");qqline(Mx15, col="red")
qqnorm(Mx20, main="n=20");qqline(Mx20, col="red")
qqnorm(Mx30, main="n=30");qqline(Mx30, col="red")
qqnorm(Mx50, main="n=50");qqline(Mx50, col="red")

qqnorm(Mx60, main="n=60");qqline(Mx60, col="red")
qqnorm(Mx100, main="n=100");qqline(Mx100, col="red")
qqnorm(Mx200, main="n=200");qqline(Mx200, col="red")
qqnorm(Mx500, main="n=500");qqline(Mx500, col="red")

Conclusión: Se observa que al utilizar una distribución binomial y al aumentar el tamaño de la muestra esta se va ir ajustando a una distribución normal, lo cual se puede evidenciar con los graficos de histograma (presentando una distribución simétrica), poligono de frecuencias y el ggplot, con la prueba de Shapiro Wilk, la muestra cumple con el supuesto de normalidad cuando se tiene un n de 200 (p valor=0.08) y un n de 500 (p valor=0.22).

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.

Simulación con 10% de plantas enfermas

set.seed(123)
n=1000
#x=rbinom(n, size=1, prob=0.5)
Lote=c(rep("enferma",100),rep("sana",900))
table(Lote)/n

Lote
enferma    sana 
    0.1     0.9 

#función
muestra=function(n) {
  m=rbinom(n, size=1, prob=0.1)
  p=sum(m)/n
  return(p)
}

cat("**Estimador de la proporción muestral pˆ con n=500**")

**Estimador de la proporción muestral pˆ con n=500**

set.seed(123)
muestra(500)

[1] 0.096

replica=rep(500,500)
x500=sapply(replica, muestra)
hist(x500, main="Histograma n=500", col="#60AEDF", xlab="Probabilidad de muestras")
line = mean(x500)
abline (v=line, lwd = 4, lty = 2, col="darkblue")

qqnorm(x500)
qqline(x500, col="red")

cat("Prueba de normalidad")

Prueba de normalidad

shapiro.test(x500)

    Shapiro-Wilk normality test

data:  x500
W = 0.99315, p-value = 0.02237

cat("Análisis descriptivo")

Análisis descriptivo

summarytools::descr(x500)

Descriptive Statistics  
x500  
N: 500  

                      x500
----------------- --------
             Mean     0.10
          Std.Dev     0.01
              Min     0.06
               Q1     0.09
           Median     0.10
               Q3     0.11
              Max     0.14
              MAD     0.01
              IQR     0.02
               CV     0.14
         Skewness     0.11
      SE.Skewness     0.11
         Kurtosis    -0.37
          N.Valid   500.00
        Pct.Valid   100.00

set.seed(123)
n=500
m=1000*n
###
X=matrix(rbinom(m, size=1, prob=0.1), ncol=n) 
#generación de muestras
X5=X[,1:5]
X10=X[,1:10]
X15=X[,1:15]
X20=X[,1:20]
X30=X[,1:30]
X50=X[,1:50]
X60=X[,1:60]
X100=X[,1:100]
X200=X[,1:200]
X500=X[,1:500]
#Obención de las medias
Mx5=apply(X5,1,mean)
Mx10=apply(X10,1,mean)
Mx15=apply(X15,1,mean)
Mx20=apply(X20,1,mean)
Mx30=apply(X30,1,mean)
Mx50=apply(X50,1,mean)
Mx60=apply(X60,1,mean)
Mx100=apply(X100,1,mean)
Mx200=apply(X200,1,mean)
Mx500=apply(X500,1,mean)
#generación de densidad empirica
d5=density(Mx5)
d10=density(Mx10)
d15=density(Mx15)
d20=density(Mx20)
d30=density(Mx30)
d50=density(Mx50)
d60=density(Mx60)
d100=density(Mx100)
d200=density(Mx200)
d500=density(Mx500)
#plot
par(cex=0.8, cex.axis=1, cex.lab=.8, cex.main=.8, cex.sub=.8,  mfrow=c(2,2), mai = c(.3, .3, .3, .3))

hist(Mx5, freq = F, xlab="n=5", col = "violet")
hist(Mx10, freq = F, xlab="n=10", col = "violet")
hist(Mx15, freq = F, xlab="n=15", col = "violet")
hist(Mx20, freq = F, xlab="n=20", col = "violet")

hist(Mx30, freq = F, xlab="n=30", col = "violet")
hist(Mx50, freq = F, xlab="n=50", col = "violet")
hist(Mx60, freq = F, xlab="n=60", col = "violet")
hist(Mx100, freq = F, xlab="n=100", col = "violet")

hist(Mx200, freq = F, xlab="n=200", col = "violet")
hist(Mx500, freq = F, xlab="n=500", col = "violet")

#graficos de densidad
plot(d5, main="n=5", xlab="n=5")
plot(d10, main="n=10", xlab="n=10")

plot(d15, main="n=15", xlab="n=15")
plot(d20, main="n=20", xlab="n=20")
plot(d30, main="n=30", xlab="n=30")
plot(d50, main="n=50", xlab="n=50")

plot(d60, main="n=60", xlab="n=60")
plot(d100, main="n=100", xlab="n=100")
plot(d200, main="n=200", xlab="n=200")
plot(d500, main="n=500", xlab="n=500")

cat("Prueba de normalidad")

Prueba de normalidad

shapiro.test(Mx5)

    Shapiro-Wilk normality test

data:  Mx5
W = 0.69827, p-value < 2.2e-16

shapiro.test(Mx10)

    Shapiro-Wilk normality test

data:  Mx10
W = 0.83372, p-value < 2.2e-16

shapiro.test(Mx20)

    Shapiro-Wilk normality test

data:  Mx20
W = 0.91931, p-value < 2.2e-16

shapiro.test(Mx30)

    Shapiro-Wilk normality test

data:  Mx30
W = 0.95139, p-value < 2.2e-16

shapiro.test(Mx50)

    Shapiro-Wilk normality test

data:  Mx50
W = 0.97662, p-value = 1.329e-11

shapiro.test(Mx60)

    Shapiro-Wilk normality test

data:  Mx60
W = 0.98039, p-value = 2.355e-10

shapiro.test(Mx100)

    Shapiro-Wilk normality test

data:  Mx100
W = 0.98726, p-value = 1.228e-07

shapiro.test(Mx200)

    Shapiro-Wilk normality test

data:  Mx200
W = 0.99346, p-value = 0.0002248

shapiro.test(Mx500)

    Shapiro-Wilk normality test

data:  Mx500
W = 0.99561, p-value = 0.005827

# grafico de normalidad
par(cex=0.8, cex.axis=.8, cex.lab=.8, cex.main=.8, cex.sub=.5,  mfrow=c(2,3), mai = c(.3, .3, .3, .3))
qqnorm(Mx5, main="n=5");qqline(Mx5, col="red")
qqnorm(Mx10, main="n=10");qqline(Mx10, col="red")
qqnorm(Mx15, main="n=15");qqline(Mx15, col="red")
qqnorm(Mx20, main="n=20");qqline(Mx20, col="red")
qqnorm(Mx30, main="n=30");qqline(Mx30, col="red")
qqnorm(Mx50, main="n=50");qqline(Mx50, col="red")

qqnorm(Mx60, main="n=60");qqline(Mx60, col="red")
qqnorm(Mx100, main="n=100");qqline(Mx100, col="red")
qqnorm(Mx200, main="n=200");qqline(Mx200, col="red")
qqnorm(Mx500, main="n=500");qqline(Mx500, col="red")

Simulación con 90% de plantas enfermas

set.seed(123)
n=1000
#x=rbinom(n, size=1, prob=0.5)
Lote=c(rep("enferma",900),rep("sana",100))
table(Lote)/n

Lote
enferma    sana 
    0.9     0.1 

#función
muestra=function(n) {
  m=rbinom(n, size=1, prob=0.9)
  p=sum(m)/n
  return(p)
}

cat("**Estimador de la proporción muestral pˆ con n=500**")

**Estimador de la proporción muestral pˆ con n=500**

set.seed(123)
muestra(500)

[1] 0.904

replica=rep(500,500)
x500=sapply(replica, muestra)
hist(x500, main="Histograma n=500", col="#10AEDF", xlab="Probabilidad de muestras")
line = mean(x500)
abline (v=line, lwd = 4, lty = 2, col="darkblue")

qqnorm(x500)
qqline(x500, col="red")

cat("Prueba de normalidad")

Prueba de normalidad

shapiro.test(x500)

    Shapiro-Wilk normality test

data:  x500
W = 0.99315, p-value = 0.02237

cat("Análisis descriptivo")

Análisis descriptivo

summarytools::descr(x500)

Descriptive Statistics  
x500  
N: 500  

                      x500
----------------- --------
             Mean     0.90
          Std.Dev     0.01
              Min     0.86
               Q1     0.89
           Median     0.90
               Q3     0.91
              Max     0.94
              MAD     0.01
              IQR     0.02
               CV     0.02
         Skewness    -0.11
      SE.Skewness     0.11
         Kurtosis    -0.37
          N.Valid   500.00
        Pct.Valid   100.00

set.seed(123)
n=500
m=1000*n
###
X=matrix(rbinom(m, size=1, prob=0.9), ncol=n) 
#generación de muestras
X5=X[,1:5]
X10=X[,1:10]
X15=X[,1:15]
X20=X[,1:20]
X30=X[,1:30]
X50=X[,1:50]
X60=X[,1:60]
X100=X[,1:100]
X200=X[,1:200]
X500=X[,1:500]
#Obención de las medias
Mx5=apply(X5,1,mean)
Mx10=apply(X10,1,mean)
Mx15=apply(X15,1,mean)
Mx20=apply(X20,1,mean)
Mx30=apply(X30,1,mean)
Mx50=apply(X50,1,mean)
Mx60=apply(X60,1,mean)
Mx100=apply(X100,1,mean)
Mx200=apply(X200,1,mean)
Mx500=apply(X500,1,mean)
#generación de densidad empirica
d5=density(Mx5)
d10=density(Mx10)
d15=density(Mx15)
d20=density(Mx20)
d30=density(Mx30)
d50=density(Mx50)
d60=density(Mx60)
d100=density(Mx100)
d200=density(Mx200)
d500=density(Mx500)
#plot
par(cex=0.8, cex.axis=1, cex.lab=.8, cex.main=.8, cex.sub=.8,  mfrow=c(2,2), mai = c(.3, .3, .3, .3))

hist(Mx5, freq = F, xlab="n=5", col = "violet")
hist(Mx10, freq = F, xlab="n=10", col = "violet")
hist(Mx15, freq = F, xlab="n=15", col = "violet")
hist(Mx20, freq = F, xlab="n=20", col = "violet")

hist(Mx30, freq = F, xlab="n=30", col = "violet")
hist(Mx50, freq = F, xlab="n=50", col = "violet")
hist(Mx60, freq = F, xlab="n=60", col = "violet")
hist(Mx100, freq = F, xlab="n=100", col = "violet")

hist(Mx200, freq = F, xlab="n=200", col = "violet")
hist(Mx500, freq = F, xlab="n=500", col = "violet")

#graficos de densidad
plot(d5, main="n=5", xlab="n=5")
plot(d10, main="n=10", xlab="n=10")

plot(d15, main="n=15", xlab="n=15")
plot(d20, main="n=20", xlab="n=20")
plot(d30, main="n=30", xlab="n=30")
plot(d50, main="n=50", xlab="n=50")

plot(d60, main="n=60", xlab="n=60")
plot(d100, main="n=100", xlab="n=100")
plot(d200, main="n=200", xlab="n=200")
plot(d500, main="n=500", xlab="n=500")

cat("Prueba de normalidad")

Prueba de normalidad

shapiro.test(Mx5)

    Shapiro-Wilk normality test

data:  Mx5
W = 0.69827, p-value < 2.2e-16

shapiro.test(Mx10)

    Shapiro-Wilk normality test

data:  Mx10
W = 0.83372, p-value < 2.2e-16

shapiro.test(Mx20)

    Shapiro-Wilk normality test

data:  Mx20
W = 0.91931, p-value < 2.2e-16

shapiro.test(Mx30)

    Shapiro-Wilk normality test

data:  Mx30
W = 0.95139, p-value < 2.2e-16

shapiro.test(Mx50)

    Shapiro-Wilk normality test

data:  Mx50
W = 0.97662, p-value = 1.329e-11

shapiro.test(Mx60)

    Shapiro-Wilk normality test

data:  Mx60
W = 0.98039, p-value = 2.355e-10

shapiro.test(Mx100)

    Shapiro-Wilk normality test

data:  Mx100
W = 0.98726, p-value = 1.228e-07

shapiro.test(Mx200)

    Shapiro-Wilk normality test

data:  Mx200
W = 0.99346, p-value = 0.0002248

shapiro.test(Mx500)

    Shapiro-Wilk normality test

data:  Mx500
W = 0.99561, p-value = 0.005827

# grafico de normalidad
par(cex=0.8, cex.axis=.8, cex.lab=.8, cex.main=.8, cex.sub=.5,  mfrow=c(2,3), mai = c(.3, .3, .3, .3))
qqnorm(Mx5, main="n=5");qqline(Mx5, col="red")
qqnorm(Mx10, main="n=10");qqline(Mx10, col="red")
qqnorm(Mx15, main="n=15");qqline(Mx15, col="red")
qqnorm(Mx20, main="n=20");qqline(Mx20, col="red")
qqnorm(Mx30, main="n=30");qqline(Mx30, col="red")
qqnorm(Mx50, main="n=50");qqline(Mx50, col="red")

qqnorm(Mx60, main="n=60");qqline(Mx60, col="red")
qqnorm(Mx100, main="n=100");qqline(Mx100, col="red")
qqnorm(Mx200, main="n=200");qqline(Mx200, col="red")
qqnorm(Mx500, main="n=500");qqline(Mx500, col="red")

Conclusión:

Se observa que para lotes con 10% de plantas enfermas, con tamaños de muestras pequeñas estas presentan una distribución asimetrica positiva y a medida que aumenta el tamaño de muestra esta se ajusta a una distribución aproximadamente simetrica lo cual se puede evidenciar con el hitograma y el qqplot, sin embargo con la prueba de Shapiro Wilk las muestras no cumplen con el parámetro de normalidad (p valor<0.05). Para los lotes con 90% de plantas enfermas, con tamaños de muestras pequeñas estas presentan una distribución asimetrica negativa y a medida que aumenta el tamaño de muestra esta se ajusta a una distribución aproximadamente simetrica lo cual se puede evidenciar con el hitograma y el qqplot, sin embargo con la prueba de Shapiro Wilk las muestras no cumplen con el parámetro de normalidad (p valor<0.05).