Modulo 2 - Actividad 1

Inferencia Estadística y Simulación

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

• Realice una simulación en la cual genere una población de N=1000 (Lote) y además que el porcentaje de individuos (plantas) enfermas sea del 50%.

genPoblación = function(Lote, Enfermas){
  porEnfermas = 0
  while (porEnfermas != Enfermas) {
    Plantas = vector()
    Plantas = sample(x = c(0,1), size = Lote, replace = TRUE)
    porEnfermas = sum(Plantas == 1)/Lote
  }
  
  return(Plantas)
}

Población = genPoblación(Lote = 1000, Enfermas = 0.5)
Población[1:76]

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

• Genere una función que permita obtener una muestra aleatoria de la población y calcule el estimador de la proporción muestral para un tamaño de muestra dado n.

genMuestra = function(limInf, limSup, nRep){
  Población = genPoblación(Lote = 1000, Enfermas = 0.5)
  
  estEnfermas = vector()
  nMuestra = sample(x = c(limInf:limSup), size = 1)
  
  for (i in 1:nRep) {
    Muestra = sample(x = Población, size = nMuestra)
    estEnfermas[i] = sum(Muestra == 1)/nMuestra
  }
  
  my.list = list(estEnfermas, nMuestra)
   
  return(my.list)
}

Estimador = genMuestra(limInf = 30, limSup = 50, nRep = 1)
print(paste("Para un tamaño de muestra de", Estimador[[2]], "se obtuvo un estimador de", Estimador[[1]]))

## [1] "Para un tamaño de muestra de 32 se obtuvo un estimador de 0.46875"

• Repita el escenario anterior (b) 500 veces y analice los resultados en cuanto al comportamiento de los 500 estimadores. ¿Qué tan simétricos son los datos?, ¿Son sesgados y qué pasa en cuanto a variabilidad?

Estimador = genMuestra(limInf = 30, limSup = 50, nRep = 500)
normEstimador = dnorm(x = Estimador[[1]], mean = mean(Estimador[[1]]), sd = sd(Estimador[[1]]))

require(ggplot2)
ggplot(mapping = aes(x = Estimador[[1]])) + labs(x = "Estimador") + ggplot2::geom_histogram(mapping = aes(y = ..density..), bins = 40) + geom_vline(xintercept = mean(Población), col = "red") + geom_line(mapping = aes(y = normEstimador), col = "blue") + geom_text(aes(x = 0.32, label = "Media Parámetro", y = 6), col = "red", angle = 0, text = element_text(size = 9)) + geom_text(aes(x = 0.32, label = "Dist. Normal Estimador", y = 7), col = "blue", angle = 0, text = element_text(size = 9))

Resultados Comparativos Parámetro vs Estimador

Media.Parámetro	SD.Parámetro	Media.Estimador	SD.Estimador	Tamaño.Muestra
0.5	0.5002502	0.4964	0.0920603	30

Como se puede apreciar en la figura anterior, los resultados no son simétricos, además, el estimador no varía continuamente y las desviaciones estándard se encuentran alejadas entre sí. La media del estimador se aproxima en gran proporción a la media de la población con un valor próximo a 0.5.

• Realice los ejercicios completos b y c para tamaños de muestra n=5, 10, 15, 20, 30, 50, 60, 100, 200, 500. Y compare los resultados de los estimadores en cuanto a la normalidad. Investigue y utilice pruebas de bondad y ajuste (shapiro wilks) y métodos gráficos (grafico qq de normalidad).

genMuestraD = function(vMuestra, nRep){
  estEnfermas = data.frame()
  lMuestra = length(vMuestra)
  nMuestra = as.vector(x = vMuestra)
  for (i in 1:lMuestra) {
    for (j in 1:nRep) {
      Muestra = sample(x = Población, size = nMuestra[i])
      estEnfermas[j, i] = sum(Muestra == 1)/nMuestra[i]
    }
  }
  
  return(estEnfermas)
}

resEstimadores = genMuestraD(vMuestra = c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500), nRep = 500)

resEstimadores = data.frame("n 5" = resEstimadores[, 1], "n 10" = resEstimadores[, 2], "n 15" = resEstimadores[, 3], "n 20" = resEstimadores[, 4], "n 30" = resEstimadores[, 5], "n 50" = resEstimadores[, 6], "n 60" = resEstimadores[, 7], "n 100" = resEstimadores[, 8], "n 200" = resEstimadores[, 9], "n 500" = resEstimadores[, 10])

Análisis.Normalidad = function(nCol, nMuestra){
  
  Gráfica.1 = ggplot(data = resEstimadores, mapping = aes(x = resEstimadores[, nCol])) %>%
    + geom_histogram(mapping = aes(y = ..density..), bins = 40) %>%
    + geom_vline(xintercept = mean(Población), col = "red", size = 1) %>%
    + labs(x = "Estimador", y = "Densidad", title = (paste("Histograma para n =", nMuestra)))
  
  Gráfica.2 = ggplot(data = resEstimadores, mapping = aes(sample = resEstimadores[, nCol])) + stat_qq() + stat_qq_line(col = "red", size = 1) + labs(x = "Theoretical Quantiles", y = "Sample Quantiles", title = (paste("Distribución para n =", nMuestra)))
           
  ShapiroWilks = shapiro.test(resEstimadores[, nCol])
  
  Tabla.1 = data.frame("Tamaño Muestra" = Estimador[[2]], "Media Estimador" = mean(resEstimadores[, nCol]), "P value" = format(ShapiroWilks$p.value, digits = 3))
  
  Gráficas = ggarrange(Gráfica.1, Gráfica.2, ncol = 2, nrow = 1, widths = c(2,3))
  
  my.list = list(Gráficas, Tabla.1)
  
  return(my.list)
  
}

require(ggpubr)

Resultado.Análisis = Análisis.Normalidad(1, 5)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.4908	2.56e-15

Como se observa en la gráfica, los valores del estimador para una muestra n=5 no son continuos, sin embargo la media del estimador se aproxima al valor del parámetro que es igual a 0.5. Por otra parte, el valor de P es muy bajo y apoyándonos en la grafica observamos como los valores se encuentran dispersos alrededor de la línea teórica de la normal.

Resultado.Análisis = Análisis.Normalidad(2, 10)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.4974	2.88e-10

Los resultados obtenidos con una muestra de n = 10 no difieren en mayor medida del análisis realizado para una muestra de n= 5. Aunque el valor de P es mayor que en el punto anterior, aún sigue siendo muy bajo.

Resultado.Análisis = Análisis.Normalidad(3, 15)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.5081333	6.71e-08

Se aprecia como poco a poco los valores se acercan la línea teórica de la normal. Se van generando valores mayores del estimador, aumentando en muy poca medida la continuidad.

Resultado.Análisis = Análisis.Normalidad(4, 20)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.502	3.09e-06

Se observa cierta simetría en relación a una función normal alrededor de la media del parámetro. La proyección es que a medida que la muestra vaya aumentando el valor de P va a seguir acercándose al valor ideal.

Resultado.Análisis = Análisis.Normalidad(5, 30)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.5018	0.000366

Como era de esperarse, el valor de P sigue incrementando poco a poco y el estimador tiene valores más continuos como se observa en el histograma.

Resultado.Análisis = Análisis.Normalidad(6, 50)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.50092	0.000593

El estimador presenta una forma más simétrica, esto debido a los valores están siendo más continuos en el eje x. Ya que la población está restringida al 50% de plantas enfermas, el estimador está constantemente oscilando sobre el 0.5.

Resultado.Análisis = Análisis.Normalidad(7, 60)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.5024	0.0076

Con un número de muestras de 60 sobre la población, se observa una gráfica mucho más continua y uniforme al valor de la media de la población. El valor de P en este punto ha mejorado ampliamente comparándolo con la primera muestra de 5.

Resultado.Análisis = Análisis.Normalidad(8, 100)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.49764	0.0419

Con algunas excepciones, podemos interpretar que la gráfica es simétrica y continua alrededor de la media de la población. Gráficamente se aprecia como los valores se han alineado un poco más con respecto la línea teórica de la normal.

Resultado.Análisis = Análisis.Normalidad(9, 200)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.49709	0.134

Los cambios alrededor del valor de P siguen siendo continuos y los valores del estimador siguen su alineación con la línea teórica de la normal.

Resultado.Análisis = Análisis.Normalidad(10, 500)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.499024	0.283

En este punto vemos una gráfica del estimador con forma simétrica alrededor de la media de población con valores continuos; además, el valor de P a alcanzado un valor relativamente cercano al esperado, con un valor de 0.254.

• Repita toda la simulación (puntos a – d) pero ahora con lotes con 10% y 90% de plantas enfermas. Concluya todo el ejercicio.

10 % de plantas enfermas:

Población = sample(c(rep(1, 100), rep(0, 900)))
resEstimadores = genMuestraD(vMuestra = c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500), nRep = 500)

resEstimadores = data.frame("n 5" = resEstimadores[, 1], "n 10" = resEstimadores[, 2], "n 15" = resEstimadores[, 3], "n 20" = resEstimadores[, 4], "n 30" = resEstimadores[, 5], "n 50" = resEstimadores[, 6], "n 60" = resEstimadores[, 7], "n 100" = resEstimadores[, 8], "n 200" = resEstimadores[, 9], "n 500" = resEstimadores[, 10])

require(ggpubr)

Resultado.Análisis = Análisis.Normalidad(1, 5)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.0984	1.15e-28

Resultado.Análisis = Análisis.Normalidad(2, 10)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.1028	1.21e-21

Resultado.Análisis = Análisis.Normalidad(3, 15)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.1002667	7.29e-18

Resultado.Análisis = Análisis.Normalidad(4, 20)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.0972	3.2e-15

Resultado.Análisis = Análisis.Normalidad(5, 30)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.0972	2.99e-11

Resultado.Análisis = Análisis.Normalidad(6, 50)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.10204	1.01e-08

Resultado.Análisis = Análisis.Normalidad(7, 60)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.0988667	2.55e-07

Resultado.Análisis = Análisis.Normalidad(8, 100)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.10072	4.91e-05

Resultado.Análisis = Análisis.Normalidad(9, 200)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.10004	0.0131

Resultado.Análisis = Análisis.Normalidad(10, 500)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.100388	0.0295

90 % de plantas enfermas:

Población = sample(c(rep(1, 900), rep(0, 100)))
resEstimadores = genMuestraD(vMuestra = c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500), nRep = 500)

resEstimadores = data.frame("n 5" = resEstimadores[, 1], "n 10" = resEstimadores[, 2], "n 15" = resEstimadores[, 3], "n 20" = resEstimadores[, 4], "n 30" = resEstimadores[, 5], "n 50" = resEstimadores[, 6], "n 60" = resEstimadores[, 7], "n 100" = resEstimadores[, 8], "n 200" = resEstimadores[, 9], "n 500" = resEstimadores[, 10])

require(ggpubr)

Resultado.Análisis = Análisis.Normalidad(1, 5)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.8924	5.43e-28

Resultado.Análisis = Análisis.Normalidad(2, 10)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.8984	1e-21

Resultado.Análisis = Análisis.Normalidad(3, 15)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.8989333	1.44e-17

Resultado.Análisis = Análisis.Normalidad(4, 20)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.9037	5.21e-16

Resultado.Análisis = Análisis.Normalidad(5, 30)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.9004667	1.36e-12

Resultado.Análisis = Análisis.Normalidad(6, 50)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.89972	4.23e-08

Resultado.Análisis = Análisis.Normalidad(7, 60)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.9027	2.42e-07

Resultado.Análisis = Análisis.Normalidad(8, 100)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.90078	3.86e-06

Resultado.Análisis = Análisis.Normalidad(9, 200)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.90089	0.00024

Resultado.Análisis = Análisis.Normalidad(10, 500)
Resultado.Análisis[[1]]

Resultado.Análisis[[2]] %>%
  kable(caption = "Comparación de resultados P value", align = 'c') %>%
  kable_paper(bootstrap_options = "striped") 

Comparación de resultados P value

Tamaño.Muestra	Media.Estimador	P.value
30	0.899652	0.0742

Conclusiones Punto 1

• Los valores pequeños de muestras no arrojan simetría alrededor de la media de la población.
• Mientras las muestras sean más pequeñas el estimador tendrá menos valores continuos, lo que se refleja en vanos dentro de una gráfica.
• Lo mismo sucede con el valor de P, en la medida que las muestras sean más pequeñas este valor se aleja mucho más del valor esperado, de las simulaciones encontramos valores de hasta 3.33e-28.
• A medida que aumentan las muestras el valor de P se acerca más a un valor óptimo.
• Gráficamente, los resultados con muestras más pequeñas se dispersan alrededor de la línea teórica de la normal; mientras que cuando las muestras son mayores estos valores se agrupan sobre la línea teórica de la normal.