Taller Simulación en R

Taller Simulación en R

Punto 1: Teorema del Límite Central

Literal a

Simulación de poblacion N=1000 con enfermo 50% y sano 50%

f_lote=function(en,sa){
  lote=c(rep("enfermo",en),rep("sano",sa))
  lote<<-sample(lote)
  return(head(lote))
}

f_lote(500,500)

## [1] "enfermo" "enfermo" "enfermo" "enfermo" "enfermo" "enfermo"

Literal b

Función para muestra aleatoria de la población (para este caso n =7) y cálculo del estimador de la proporción muestral para un tamaño de muestra (para este caso n =7).

estimador = function(n){
  muestra=sample(lote,size=n)
  gorro=sum(muestra == "enfermo")/n
  return(gorro)
}

estimador(7)

## [1] 0.4285714

Literal c

Escenario (b) 500 veces y analisis los resultados en cuanto al comportamiento de los 500 estimadores.

escen_b = function(n,repetic){
  escenarios<<-sapply(rep(n,repetic), estimador)
  return(hist(escenarios))
}

escen_b(7,500)

¿Qué tan simétricos son los datos?

R// Los datos no son simétricos ni demuestran normalidad. Existen probabilidades que no tienen frecuencia genrando vacios y no se presente la campana de Gauss.

¿Son sesgados y qué pasa en cuanto a variabilidad?

R// Los datos no son simetricos, lejos del 50% por lo cual entre cada calculo de una muestra pequeña no se puede determinar estos datos.

Literal d

Se realizan los cálculos para siguientes tamaños de muestra comparando los resultados de los estimadores en cuanto a la normalidad

**_n=5, 10, 15, 20, 30, 50, 60, 100, 200, 500.**

Se puede apreciar que a medida que la muestra de los estimadores aumenta mas ajustado se encuetra la información a la normalidad.

  par(mfrow=c(2,2))
  
  escen_b(5,500)
  title(sub="Muestra_5")
  b1= data.frame(escenarios)
  escen_b(10,500)
  title(sub="Muestra_10")
  b2=data.frame(escenarios)
  escen_b(15,500)
  title(sub="Muestra_15")
  b3=data.frame(escenarios)
  escen_b(20,500)
  title(sub="Muestra_20")

  b4=data.frame(escenarios)
  escen_b(30,500)
  title(sub="Muestra_30")
  b5=data.frame(escenarios)
  escen_b(50,500)
  title(sub="Muestra_50")
  b6=data.frame(escenarios)
  escen_b(60,500)
  title(sub="Muestra_60")
  b7=data.frame(escenarios)
  escen_b(100,500)
  title(sub="Muestra_100")

  b8=data.frame(escenarios)
  escen_b(200,500)
  title(sub="Muestra_200")
  b9=data.frame(escenarios)
  escen_b(500,500)
  title(sub="Muestra_500")
  b10=data.frame(escenarios)

Se corrobora la normalidad a métodos gráficos qq de normalidad

par(mfrow=c(2,2))

qqnorm(b1$escenarios)
qqline(b1$escenarios)
qqnorm(b2$escenarios)
qqline(b2$escenarios) 
qqnorm(b3$escenarios)
qqline(b3$escenarios)
qqnorm(b4$escenarios)
qqline(b4$escenarios)

qqnorm(b5$escenarios)
qqline(b5$escenarios)
qqnorm(b6$escenarios)
qqline(b6$escenarios)
qqnorm(b7$escenarios)
qqline(b7$escenarios)
qqnorm(b8$escenarios)
qqline(b8$escenarios)

qqnorm(b9$escenarios)
qqline(b9$escenarios)
qqnorm(b10$escenarios)
qqline(b10$escenarios)

Para comparar estas gráficas con el resultado se realiza la verificación de las proporciones 50% enfermo:

Como el p-valor es muy alto (1), no se rechaza la hipótesis nula de que p=0.5.

prop.test(x=500,n= 1000, conf.level = 0.9)

## 
##  1-sample proportions test without continuity correction
## 
## data:  500 out of 1000, null probability 0.5
## X-squared = 0, df = 1, p-value = 1
## alternative hypothesis: true p is not equal to 0.5
## 90 percent confidence interval:
##  0.4740277 0.5259723
## sample estimates:
##   p 
## 0.5

Literal e

Se repite toda la simulación (puntos a – d) pero con lote 10% de plantas enfermas.

## Lote de 10% enfermas
f_lote(100,900)

## [1] "enfermo" "enfermo" "enfermo" "enfermo" "enfermo" "enfermo"

## Se inicia un estimador de muestra 7
estimador(7)

## [1] 0

## Se plantea un escenario de muestra 7 unas 500 repeticiones
escen_b(7,500)

## Se realizan los cálculos para siguientes tamaños de muestra n=5, 10, 15, 20, 30, 50, 60, 100, 200 y 500.
par(mfrow=c(2,2))
  
escen_b(5,500)
title(sub="Muestra_5")
b1= data.frame(escenarios)
escen_b(10,500)
title(sub="Muestra_10")
b2=data.frame(escenarios)
escen_b(15,500)
title(sub="Muestra_15")
b3=data.frame(escenarios)
escen_b(20,500)
title(sub="Muestra_20")

b4=data.frame(escenarios)
escen_b(30,500)
title(sub="Muestra_30")
b5=data.frame(escenarios)
escen_b(50,500)
title(sub="Muestra_50")
b6=data.frame(escenarios)
escen_b(60,500)
title(sub="Muestra_60")
b7=data.frame(escenarios)
escen_b(100,500)
title(sub="Muestra_100")

b8=data.frame(escenarios)
escen_b(200,500)
title(sub="Muestra_200")
b9=data.frame(escenarios)
escen_b(500,500)
title(sub="Muestra_500")
b10=data.frame(escenarios)

##Se corrobora la normalidad a métodos gráficos qq de normalidad

par(mfrow=c(2,2))

qqnorm(b1$escenarios)
qqline(b1$escenarios)
qqnorm(b2$escenarios)
qqline(b2$escenarios) 
qqnorm(b3$escenarios)
qqline(b3$escenarios)
qqnorm(b4$escenarios)
qqline(b4$escenarios)

qqnorm(b5$escenarios)
qqline(b5$escenarios)
qqnorm(b6$escenarios)
qqline(b6$escenarios)
qqnorm(b7$escenarios)
qqline(b7$escenarios)
qqnorm(b8$escenarios)
qqline(b8$escenarios)

qqnorm(b9$escenarios)
qqline(b9$escenarios)
qqnorm(b10$escenarios)
qqline(b10$escenarios)

Para comparar estas gráficas con el resultado se realiza la verificación de las proporciones 10% enfermo:

Como el p-valor es muy alto (2.2), no se rechaza la hipótesis nula de que p=0.1.

prop.test(x=100,n= 1000, conf.level = 0.9)

## 
##  1-sample proportions test with continuity correction
## 
## data:  100 out of 1000, null probability 0.5
## X-squared = 638.4, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.5
## 90 percent confidence interval:
##  0.08499444 0.11723306
## sample estimates:
##   p 
## 0.1

Se repite toda la simulación (puntos a – d) pero con lote 90% de plantas enfermas.

## Lote de 90% enfermas
f_lote(900,100)

## [1] "enfermo" "enfermo" "enfermo" "enfermo" "enfermo" "enfermo"

## Se inicia un estimador de muestra 7
estimador(7)

## [1] 0.8571429

## Se plantea un escenario de muestra 7 unas 500 repeticiones
escen_b(7,500)

## Se realizan los cálculos para siguientes tamaños de muestra n=5, 10, 15, 20, 30, 50, 60, 100, 200 y 500.
par(mfrow=c(2,2))
  
escen_b(5,500)
title(sub="Muestra_5")
b1= data.frame(escenarios)
escen_b(10,500)
title(sub="Muestra_10")
b2=data.frame(escenarios)
escen_b(15,500)
title(sub="Muestra_15")
b3=data.frame(escenarios)
escen_b(20,500)
title(sub="Muestra_20")

b4=data.frame(escenarios)
escen_b(30,500)
title(sub="Muestra_30")
b5=data.frame(escenarios)
escen_b(50,500)
title(sub="Muestra_50")
b6=data.frame(escenarios)
escen_b(60,500)
title(sub="Muestra_60")
b7=data.frame(escenarios)
escen_b(100,500)
title(sub="Muestra_100")

b8=data.frame(escenarios)
escen_b(200,500)
title(sub="Muestra_200")
b9=data.frame(escenarios)
escen_b(500,500)
title(sub="Muestra_500")
b10=data.frame(escenarios)

##Se corrobora la normalidad a métodos gráficos qq de normalidad

par(mfrow=c(2,2))

qqnorm(b1$escenarios)
qqline(b1$escenarios)
qqnorm(b2$escenarios)
qqline(b2$escenarios) 
qqnorm(b3$escenarios)
qqline(b3$escenarios)
qqnorm(b4$escenarios)
qqline(b4$escenarios)

qqnorm(b5$escenarios)
qqline(b5$escenarios)
qqnorm(b6$escenarios)
qqline(b6$escenarios)
qqnorm(b7$escenarios)
qqline(b7$escenarios)
qqnorm(b8$escenarios)
qqline(b8$escenarios)

qqnorm(b9$escenarios)
qqline(b9$escenarios)
qqnorm(b10$escenarios)
qqline(b10$escenarios)

Para comparar estas gráficas con el resultado se realiza la verificación de las proporciones 90% enfermo:

Como el p-valor es muy alto (2.2), no se rechaza la hipótesis nula de que p=0.9.

prop.test(x=900,n= 1000, conf.level = 0.9)

## 
##  1-sample proportions test with continuity correction
## 
## data:  900 out of 1000, null probability 0.5
## X-squared = 638.4, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.5
## 90 percent confidence interval:
##  0.8827669 0.9150056
## sample estimates:
##   p 
## 0.9

Punto 2: Teorema del Límite Central

Literal a

simulación de dos poblaciones de n1=1000 (Lote1) y n2=1500 (Lote2) con porcentaje de individuos (plantas) enfermas en ambos lotes de 10%.

lote1=c(rep("enfermo",100),rep("sano",900))
n1=sample(lote1)

lote2=c(rep("enfermo",150),rep("sano",1350))
n2=sample(lote2)

Literal b

Función para obtener una muestra aleatoria de los lotes, calcular el estimador de la proporción muestral n1=n2 (para este punto muestra = 500) y sacar la diferencia entre los estimadores p1-p2.

dif_estimador_n = function(n){
  muestra1=sample(lote1,size=n)
  muestra2=sample(lote2,size=n)
  gorro_n1<<-sum(muestra1 == "enfermo")/n
  gorro_n2<<-sum(muestra2 == "enfermo")/n
  gorro_dif = gorro_n1 - gorro_n2
  
  return(gorro_dif)
}

dif_estimador_n(500)

## [1] -0.036

paste("El número anterior es la diferencia entre los estimadores p1 (", gorro_n1,") - (", gorro_n2)

## [1] "El número anterior es la diferencia entre los estimadores p1 ( 0.082 ) - ( 0.118"

Literal c

Escenario 500 veces y analisis de los resultados con 500 estimadores a tráves de un histograma.

escen_b_n = function(n,repetic){
  escenarios_n<<-sapply(rep(n,repetic), dif_estimador_n)
  return(hist(escenarios_n))
}

escen_b_n(500,500)

¿Qué tan simétricos son los datos? Los datos a medida que la muestra es mas grande se establece claramente una simetria en torno al cero, lo cual es ,lo esperado puesto que se tomaron muestras de igual tamaño.

¿Son siempre cero las diferencias?

No siempre, si bien mayoria de los datos estudiados en el histograma anterior tienen a cero hay datos por encima y por debajo de cero en menor proporcion.

Literal d

Presentación de los resultados de los estimadores (p1-p2) para tamaños de muestra n1=n2=5, 10, 15, 20, 30, 50, 60, 100, 200, 500 comparados con la normalidad qqplot

  par(mfrow=c(2,2))
  
  escen_b_n(5,500)
  title(sub="Muestra_5")
  c1= data.frame(escenarios_n)
  escen_b_n(10,500)
  title(sub="Muestra_10")
  c2=data.frame(escenarios_n)
  escen_b_n(15,500)
  title(sub="Muestra_15")
  c3=data.frame(escenarios_n)
  escen_b_n(20,500)
  title(sub="Muestra_20")

  c4=data.frame(escenarios_n)
  escen_b_n(30,500)
  title(sub="Muestra_30")
  c5=data.frame(escenarios_n)
  escen_b_n(50,500)
  title(sub="Muestra_50")
  c6=data.frame(escenarios_n)
  escen_b_n(60,500)
  title(sub="Muestra_60")
  c7=data.frame(escenarios_n)
  escen_b_n(100,500)
  title(sub="Muestra_100")

  c8=data.frame(escenarios_n)
  escen_b_n(200,500)
  title(sub="Muestra_200")
  c9=data.frame(escenarios_n)
  escen_b_n(500,500)
  title(sub="Muestra_500")
  c10=data.frame(escenarios_n)

Se corrobora la normalidad a métodos gráficos qq de normalidad

par(mfrow=c(2,2))

qqnorm(c1$escenarios_n)
qqline(c1$escenarios_n)
qqnorm(c2$escenarios_n)
qqline(c2$escenarios_n) 
qqnorm(c3$escenarios_n)
qqline(c3$escenarios_n)
qqnorm(c4$escenarios_n)
qqline(c4$escenarios_n)

qqnorm(c5$escenarios_n)
qqline(c5$escenarios_n)
qqnorm(c6$escenarios_n)
qqline(c6$escenarios_n)
qqnorm(c7$escenarios_n)
qqline(c7$escenarios_n)
qqnorm(c8$escenarios_n)
qqline(c8$escenarios_n)

qqnorm(c9$escenarios_n)
qqline(c9$escenarios_n)
qqnorm(c10$escenarios_n)
qqline(c10$escenarios_n)

¿Considera que es más probable concluir que existen diferencias entre los tratamientos con muestras grandes que pequeñas, es decir, cuál considera usted que es el efecto del tamaño de muestra en el caso de la comparación de proporciones?

R// Entre mas pequeña sea una muestra las proporciones no presentan normalidad pero a medida que se amplia la camtidad de la muestra se va perfilando cada vez mas la normalidad de los datos con respecto a la media.

Literal e

Ahora realice nuevamente los puntos a-d bajo un escenario con dos lotes, pero de proporciones de enfermos diferentes (P1=0.1 y P2=0.15), es decir, el tratamiento del lote 1 si presentó un mejor desempeño reduciendo en un 5% el porcentaje de enfermos.

lote3=c(rep("enfermo",100),rep("sano",900))
n3=sample(lote1)

lote4=c(rep("enfermo",225),rep("sano",1275))
n4=sample(lote2)

Bajo este nuevo escenario se compara la distribución de estas diferencias (p1-p2)

dif_estimador_n2 = function(n){
  muestra3=sample(lote3,size=n)
  muestra4=sample(lote4,size=n)
  gorro_n3<<-sum(muestra3 == "enfermo")/n
  gorro_n4<<-sum(muestra4 == "enfermo")/n
  gorro_dif2 = gorro_n3 - gorro_n4
  
  return(gorro_dif2)
}

dif_estimador_n2(500)

## [1] -0.052

paste("El número anterior es la diferencia entre los estimadores p1 (", gorro_n3,") - (", gorro_n4)

## [1] "El número anterior es la diferencia entre los estimadores p1 ( 0.098 ) - ( 0.15"

escen_b_n2 = function(n,repetic){
  escenarios_n2<<-sapply(rep(n,repetic), dif_estimador_n2)
  return(hist(escenarios_n2))
}

escen_b_n2(500,500)

Presentación de los resultados de los estimadores (p1-p2) para tamaños de muestra n1=n2=5, 10, 15, 20, 30, 50, 60, 100, 200, 500

  par(mfrow=c(2,2))
  
  escen_b_n2(5,500)
  title(sub="Muestra_5")
  d1= data.frame(escenarios_n2)
  escen_b_n2(10,500)
  title(sub="Muestra_10")
  d2=data.frame(escenarios_n2)
  escen_b_n2(15,500)
  title(sub="Muestra_15")
  d3=data.frame(escenarios_n2)
  escen_b_n2(20,500)
  title(sub="Muestra_20")

  d4=data.frame(escenarios_n2)
  escen_b_n2(30,500)
  title(sub="Muestra_30")
  d5=data.frame(escenarios_n2)
  escen_b_n2(50,500)
  title(sub="Muestra_50")
  d6=data.frame(escenarios_n2)
  escen_b_n2(60,500)
  title(sub="Muestra_60")
  d7=data.frame(escenarios_n2)
  escen_b_n2(100,500)
  title(sub="Muestra_100")

  d8=data.frame(escenarios_n2)
  escen_b_n2(200,500)
  title(sub="Muestra_200")
  d9=data.frame(escenarios_n2)
  escen_b_n2(500,500)
  title(sub="Muestra_500")
  d10=data.frame(escenarios_n2)

con las observadas bajo igualdad de condiciones en los lotes. ¿Qué puede concluir?

La normalidad presentada mostrando el promedio sobre la diferencia de las poblaciones p1=0,1 y p2=0,15