R Notebook

Taller de Simulación en R

Autor: Jaime Reinoso

PROBLEMA 1

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.

a. 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%.

El siguente código crea una población con una proporción de plantas enfermas del 50%. Se genera la tabla según tipo de resultado al final.

estadoPlanta = c('Sana','Enferma')
N = 1000
probEnferma = 0.5
poblacionPlantas <- sample(estadoPlanta, 
                           size=N, 
                           replace=TRUE,     # se debe permitir pues si no se acaban los valores a escoger
                           prob=c(1-probEnferma, probEnferma)  # 50% cada tipo, si probEnferma = 0.5
                           )
table(poblacionPlantas)

poblacionPlantas
Enferma    Sana 
    504     496 

Como se observa, la población generada tiene (generalmente) igual cantidad de individuos enfermos como sanos. Por temas de aleatoriedad, normalmente no se obtiene una distribución perfecta 50-50 pero si es bastante cercana siempre.

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

Muestra <- function(n){
              muestra <- sample(poblacionPlantas, size=n, replace=FALSE)  # es un vector de caracteres
           
              #debemos contar cuantos plantas hay Enfermas
              y = length(muestra[muestra== 'Enferma'])
              
              # calculamos la proporción. 
              p = y/n
              return(p)
}

p = Muestra(100)

La función descrita en este punnto toma como parámetro un valor n de la cantidad de elementos que se tomarán de la población para generar la muestra. Se usa el parámetro “replace=FALSE” para evitar que un elemento sea escogido nuevamente, es decir, este solo puede ser escogido una vez.

Como demostración, se ejecuta la función muestra con 100 individuos.

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

# Se crea un vector vacío
t<-c()          

 # Se hace un ciclo 500 veces
for (i in 1:500) {      
  t<-c(t,Muestra(500))   # Se captura el indicador para ese tamaño de muestra y se introduce en el vector c, 500 veces.
}

# Ahora calculemos la prueba de bondad y ajuste (shapiro wilks) de esa serie de indicadores a ver si su distribución se parece a la normal
pruebaST <-shapiro.test(t)

# hay que instalar el paquete nortest
#install.packages("nortest")
library(nortest)


# y la prueba kolmogorov-smirnov
pruebaKS <- lillie.test(t)

stConclusion =""
if (pruebaST$p.value <0.05){
  stConclusion="Sh: Se rechaza Ho"
}else{
  stConclusion="Sh:No se rechaza Ho"
}

ksConclusion =""
if (pruebaKS$p.value <0.05){
  ksConclusion="Kol:Se rechaza Ho"
}else{
  ksConclusion="Kol: No se rechaza Ho"
}

#y hacemos un histograma para ver visualmente la distribución.

hist(t, breaks=100,main=paste("n = 100 Sh=",format(pruebaST$p.value,digits=4),"Kol=",format(pruebaKS$p.value, digits=4)),
     sub=paste(stConclusion," - ",ksConclusion))                 

RESPUESTA: Se observa que la distribución tiende ser simétrica, no sesgada, aunque la prueba Shapiro-Wilks y la de Kolmogorov-Smirnov indican que no sigue una distribución normal. Se entiende que la Shapiro-Wilks debería ser usadas para muestras menores a 50 mientras que la Kolmogorov Smirnov debería ser usada para muestras de tamaño mayor o igual a 50. El ejercicio pide Shapiro-Wilks a pesar que la muestra es pedida de tamaño 500. Por eso se muestran ambas.

d. 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).

vn<-c(5,10,15,20,30,50,60,100,200,500)   # vn es el vector que tiene los diferentes valores de n a probar


#j es el tamaño de muestra que se va a probar

for (n in vn){
  
  t<-c()  # se crea el vector vacío
  
  for (i in 1:500) {      # Se hace un ciclo 500 veces, es decir, se obtienen 500 estadísticos de muestra j
      t<-c(t,Muestra(n))   # Se captura el indicador para ese tamaño de muestra.
  }
  
  # Ahora calculemos la prueba de bondad y ajuste (shapiro wilks) de esa serie de indicadores a ver si su distribución se parece a la normal
  pruebaST <-shapiro.test(t)
  
  stConclusion =""
  if (pruebaST$p.value <0.05){
    stConclusion="Sh: Se rechaza Ho"
  }else{
    stConclusion="Sh:No se rechaza Ho"
  }

  #y hacemos un histograma para ver visualmente la distribución.
  
  hist(t, main=paste("n =",n," Sh=",format(pruebaST$p.value, digits=4)), 
     sub=paste(stConclusion), breaks=40)      
  qqnorm(t, main=paste("n =",n," Sh=",format(pruebaST$p.value, digits=4)))
  
  #hist(t, main = paste("Histograma de n =",j))  # Se publica el resultado.   
}

# nota:  R genera todos los gráficos y estos pueden ser comparados posteriormente.

Como resultado, se observa que a medida que n crece, la distribución del estadístico se acerca a la normal como lo describe el teorema del límite central. Esto también parece confirmarse con la prueba Shapiro-wilk, en donde a medida que se aumenta el tamaño de muestra (n), el p-value va creciendo y por tanto se acerca más y más a no rechazar Ho. Igualmente, al revisasr los gráficos qqnorm, se observa que el gráfica cada vez mas se acerca a la diagonal, indicando que se acerca a una distribución normal.

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

Por facilidad de discusión en este taller, se recrea el código cambiano los valores de n a los pedidos. Para esto, se crea un vector con los valores de n esperados y se incluye un for que revisa cada uno de dichos valores.

estadoPlanta = c('Sana','Enferma')
N = 1000
vProbEnferma = c(0.1, 0.9)   # vamos a hacer el mismo proceso para diferentes valores de probabilidad de plantas enfermas

#Se itera por las diferentes probabilidades de enfermedad
for (probEnferma in vProbEnferma){
  
  #Se crea la población correspondiente
  poblacionPlantas <- sample(estadoPlanta, 
                           size=N, 
                           replace=TRUE,     # se debe permitir pues si no se acaban los valores a escoger
                           prob=c(1-probEnferma, probEnferma)  # 50% cada tipo, si probEnferma = 0.5
                           )
  #Se confirma la distribución obtenida
  table(poblacionPlantas)
  
  #Se define el vector de tamaño de muestra
  vn<-c(5,10,15,20,30,50,60,100,200,500)   # vn es el vector que tiene los diferentes valores de n a probar
  
  #j es el tamaño de muestra a probar
  for (n in vn){
    
    t<-c()  # se crea el vector vacío
    
    for (i in 1:500) {      # Se hace un ciclo 500 veces
      t<-c(t,Muestra(n))   # Se captura el indicador para ese tamaño de muestra.
    }
    
    # Se obtiene la prueba shapiro para t
    pruebaST <-shapiro.test(t)
    
    #Se obtiene el histograma correspondiente.
    hist(t, main = paste(" Porcentaje Enfermas: ", probEnferma, "j =", n, "Shap:",format(pruebaST$p.value, digits=4)))  # Se publica el resultado.  
    qqnorm(t, main = paste(" Porcentaje Enfermas: ", probEnferma, "j =", n, "Shap:",format(pruebaST$p.value, digits=4)))
  }
  
}

NA
NA

En el caso de una distribución del 10% de enfermas, se observa una distribución normal con centro en 10%, a medida que n crece. Esto es consistente con la población usada.

Por otro lado, a medida que n tiende a ser más grande, se confirma que la distribución tiende a normal tanto analizando visualmente el histograma resultante, observando el resultado de la prueba Shapiro aplicada y analizando la grafica qqnorm.

en el caso de la distribución con un 90% de plantas enfermas, ocurre lo mismo, con la diferencia que la muestra tiende a ser normal con una media cercana al 90% como corresponde.

Esto convalida lo expresado por el teorema del límite central.