Taller Simulación Inferencia

Punto 1.

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

## Los 1 = Árboles Enfermos
## Los 0 = Árboles Sanos

lote = c(rep(1, 500), rep(0,500))
lote = sample(lote)

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.

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

Calcula_Estimador(n=200)

## [1] 0.475

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?

En este ejemplo, observamos que cuando se pone a iterar pocas veces, para este ejemplo 50 veces, se tiene poca simetría respecto a un comportamiento normal, pero a medida que se aumentan las iteraciones se va logrando mayor simetría a un comportamiento normal, así como se puede ver en la evolución en 200, 2000 y 10000 iteraciones.

Muestras=sapply(rep(200,50), Calcula_Estimador)
hist(Muestras, prob = TRUE, main = "Histogram with density curve")
line=mean(Muestras)
abline(v=line, col="red", lwd=4)
lines(density(Muestras), col = 2, lwd = 2)

mean(Muestras)-0.5

## [1] -0.0039

sd(Muestras)

## [1] 0.03301654

Muestras=sapply(rep(200,200), Calcula_Estimador)
hist(Muestras, prob = TRUE, main = "Histogram with density curve")
line=mean(Muestras)
abline(v=line, col="red", lwd=4)
lines(density(Muestras), col = 2, lwd = 2)

mean(Muestras)-0.5

## [1] 0.00045

sd(Muestras)

## [1] 0.03150811

Muestras=sapply(rep(200,2000), Calcula_Estimador)
hist(Muestras, prob = TRUE, main = "Histogram with density curve")
line=mean(Muestras)
abline(v=line, col="red", lwd=4)
lines(density(Muestras), col = 2, lwd = 2)

mean(Muestras)-0.5

## [1] 0.001025

sd(Muestras)

## [1] 0.03111926

Muestras=sapply(rep(200,10000), Calcula_Estimador)
hist(Muestras, prob = TRUE, main = "Histogram with density curve")
line=mean(Muestras)
abline(v=line, col="red", lwd=4)
lines(density(Muestras), col = 2, lwd = 2)

mean(Muestras)-0.5

## [1] -0.000197

sd(Muestras)

## [1] 0.03151786

Para el tamaño de muestra n=5 y al hacer 500 iteraciones, nos damos cuenta que no cumple con las pruebas de normalidad, quiere decir que se rechaza la hipótesis nula

n=5

require(MASS)
require(survival)
require(fitdistrplus)

Muestras=sapply(rep(5,500), Calcula_Estimador)
hist(Muestras, prob = TRUE, main = "Histogram with density curve")
line=mean(Muestras)
abline(v=line, col="red", lwd=4)
lines(density(Muestras), col = 2, lwd = 2)

install.packages("MASS")
install.packages("survival")
install.packages("fitdistrplus")


library(MASS)
library(survival)
library(fitdistrplus)

ajuste = fitdist(Muestras,"norm")

ajuste$estimate

##      mean        sd 
## 0.4912000 0.2308302

plot(ajuste)

n=10

require(MASS)
require(survival)
require(fitdistrplus)

Muestras=sapply(rep(10,500), Calcula_Estimador)
hist(Muestras, prob = TRUE, main = "Histogram with density curve")
line=mean(Muestras)
abline(v=line, col="red", lwd=4)
lines(density(Muestras), col = 2, lwd = 2)

install.packages("MASS")
install.packages("survival")
install.packages("fitdistrplus")


library(MASS)
library(survival)
library(fitdistrplus)

ajuste = fitdist(Muestras,"norm")

ajuste$estimate

##      mean        sd 
## 0.5050000 0.1623422

plot(ajuste)

n=15

require(MASS)
require(survival)
require(fitdistrplus)

Muestras=sapply(rep(15,500), Calcula_Estimador)
hist(Muestras, prob = TRUE, main = "Histogram with density curve")
line=mean(Muestras)
abline(v=line, col="red", lwd=4)
lines(density(Muestras), col = 2, lwd = 2)

install.packages("MASS")
install.packages("survival")
install.packages("fitdistrplus")


library(MASS)
library(survival)
library(fitdistrplus)

ajuste = fitdist(Muestras,"norm")

ajuste$estimate

##      mean        sd 
## 0.4977333 0.1292172

plot(ajuste)

n=30

require(MASS)
require(survival)
require(fitdistrplus)

Muestras=sapply(rep(30,500), Calcula_Estimador)
hist(Muestras, prob = TRUE, main = "Histogram with density curve")
line=mean(Muestras)
abline(v=line, col="red", lwd=4)
lines(density(Muestras), col = 2, lwd = 2)

install.packages("MASS")
install.packages("survival")
install.packages("fitdistrplus")


library(MASS)
library(survival)
library(fitdistrplus)

ajuste = fitdist(Muestras,"norm")

ajuste$estimate

##       mean         sd 
## 0.49800000 0.08991601

plot(ajuste)

n=60

require(MASS)
require(survival)
require(fitdistrplus)

Muestras=sapply(rep(60,500), Calcula_Estimador)
hist(Muestras, prob = TRUE, main = "Histogram with density curve")
line=mean(Muestras)
abline(v=line, col="red", lwd=4)
lines(density(Muestras), col = 2, lwd = 2)

install.packages("MASS")
install.packages("survival")
install.packages("fitdistrplus")


library(MASS)
library(survival)
library(fitdistrplus)

ajuste = fitdist(Muestras,"norm")

ajuste$estimate

##       mean         sd 
## 0.49863333 0.06608764

plot(ajuste)

n=200

require(MASS)
require(survival)
require(fitdistrplus)

Muestras=sapply(rep(200,500), Calcula_Estimador)
hist(Muestras, prob = TRUE, main = "Histogram with density curve")
line=mean(Muestras)
abline(v=line, col="red", lwd=4)
lines(density(Muestras), col = 2, lwd = 2)

install.packages("MASS")
install.packages("survival")
install.packages("fitdistrplus")


library(MASS)
library(survival)
library(fitdistrplus)

ajuste = fitdist(Muestras,"norm")

ajuste$estimate

##       mean         sd 
## 0.49872000 0.03304484

plot(ajuste)

n=500

require(MASS)
require(survival)
require(fitdistrplus)

Muestras=sapply(rep(500,500), Calcula_Estimador)
hist(Muestras, prob = TRUE, main = "Histogram with density curve")
line=mean(Muestras)
abline(v=line, col="red", lwd=4)
lines(density(Muestras), col = 2, lwd = 2)

install.packages("MASS")
install.packages("survival")
install.packages("fitdistrplus")


library(MASS)
library(survival)
library(fitdistrplus)

ajuste = fitdist(Muestras,"norm")

ajuste$estimate

##       mean         sd 
## 0.49921200 0.01588745

plot(ajuste)

Conclusión: A mayor iteraciones, mayor es la simetría con la normalidad, menor es el sesgo y muy poca variación.

10% de plantas enfermas:

n=5

## Los 1 = Árboles Enfermos
## Los 0 = Árboles Sanos

lote = c(rep(1, 100), rep(0,900))
lote = sample(lote)

require(MASS)
require(survival)
require(fitdistrplus)

Muestras=sapply(rep(5,500), Calcula_Estimador)
hist(Muestras, prob = TRUE, main = "Histogram with density curve")
line=mean(Muestras)
abline(v=line, col="red", lwd=4)
lines(density(Muestras), col = 2, lwd = 2)

install.packages("MASS")
install.packages("survival")
install.packages("fitdistrplus")


library(MASS)
library(survival)
library(fitdistrplus)

ajuste = fitdist(Muestras,"norm")

ajuste$estimate

##      mean        sd 
## 0.1036000 0.1410923

plot(ajuste)

n=20

require(MASS)
require(survival)
require(fitdistrplus)

Muestras=sapply(rep(20,500), Calcula_Estimador)
hist(Muestras, prob = TRUE, main = "Histogram with density curve")
line=mean(Muestras)
abline(v=line, col="red", lwd=4)
lines(density(Muestras), col = 2, lwd = 2)

install.packages("MASS")
install.packages("survival")
install.packages("fitdistrplus")


library(MASS)
library(survival)
library(fitdistrplus)

ajuste = fitdist(Muestras,"norm")

ajuste$estimate

##       mean         sd 
## 0.09680000 0.06899101

plot(ajuste)

n=30

require(MASS)
require(survival)
require(fitdistrplus)

Muestras=sapply(rep(30,500), Calcula_Estimador)
hist(Muestras, prob = TRUE, main = "Histogram with density curve")
line=mean(Muestras)
abline(v=line, col="red", lwd=4)
lines(density(Muestras), col = 2, lwd = 2)

install.packages("MASS")
install.packages("survival")
install.packages("fitdistrplus")


library(MASS)
library(survival)
library(fitdistrplus)

ajuste = fitdist(Muestras,"norm")

ajuste$estimate

##       mean         sd 
## 0.09966667 0.05376802

plot(ajuste)

n=50

require(MASS)
require(survival)
require(fitdistrplus)

Muestras=sapply(rep(50,500), Calcula_Estimador)
hist(Muestras, prob = TRUE, main = "Histogram with density curve")
line=mean(Muestras)
abline(v=line, col="red", lwd=4)
lines(density(Muestras), col = 2, lwd = 2)

install.packages("MASS")
install.packages("survival")
install.packages("fitdistrplus")


library(MASS)
library(survival)
library(fitdistrplus)

ajuste = fitdist(Muestras,"norm")

ajuste$estimate

##       mean         sd 
## 0.10348000 0.04137499

plot(ajuste)

n=100

require(MASS)
require(survival)
require(fitdistrplus)

Muestras=sapply(rep(100,500), Calcula_Estimador)
hist(Muestras, prob = TRUE, main = "Histogram with density curve")
line=mean(Muestras)
abline(v=line, col="red", lwd=4)
lines(density(Muestras), col = 2, lwd = 2)

install.packages("MASS")
install.packages("survival")
install.packages("fitdistrplus")


library(MASS)
library(survival)
library(fitdistrplus)

ajuste = fitdist(Muestras,"norm")

ajuste$estimate

##       mean         sd 
## 0.10122000 0.02852914

plot(ajuste)

n=200

require(MASS)
require(survival)
require(fitdistrplus)

Muestras=sapply(rep(200,500), Calcula_Estimador)
hist(Muestras, prob = TRUE, main = "Histogram with density curve")
line=mean(Muestras)
abline(v=line, col="red", lwd=4)
lines(density(Muestras), col = 2, lwd = 2)

install.packages("MASS")
install.packages("survival")
install.packages("fitdistrplus")


library(MASS)
library(survival)
library(fitdistrplus)

ajuste = fitdist(Muestras,"norm")

ajuste$estimate

##      mean        sd 
## 0.1007100 0.0183697

plot(ajuste)

n=500

require(MASS)
require(survival)
require(fitdistrplus)

Muestras=sapply(rep(500,500), Calcula_Estimador)
hist(Muestras, prob = TRUE, main = "Histogram with density curve")
line=mean(Muestras)
abline(v=line, col="red", lwd=4)
lines(density(Muestras), col = 2, lwd = 2)

install.packages("MASS")
install.packages("survival")
install.packages("fitdistrplus")


library(MASS)
library(survival)
library(fitdistrplus)

ajuste = fitdist(Muestras,"norm")

ajuste$estimate

##        mean          sd 
## 0.100432000 0.009530864

plot(ajuste)

90% de plantas enfermas:

n=5

## Los 1 = Árboles Enfermos
## Los 0 = Árboles Sanos

lote = c(rep(1, 900), rep(0,100))
lote = sample(lote)

require(MASS)
require(survival)
require(fitdistrplus)

Muestras=sapply(rep(5,500), Calcula_Estimador)
hist(Muestras, prob = TRUE, main = "Histogram with density curve")
line=mean(Muestras)
abline(v=line, col="red", lwd=4)
lines(density(Muestras), col = 2, lwd = 2)

install.packages("MASS")
install.packages("survival")
install.packages("fitdistrplus")


library(MASS)
library(survival)
library(fitdistrplus)

ajuste = fitdist(Muestras,"norm")

ajuste$estimate

##      mean        sd 
## 0.8908000 0.1368041

plot(ajuste)

n=20

require(MASS)
require(survival)
require(fitdistrplus)

Muestras=sapply(rep(20,500), Calcula_Estimador)
hist(Muestras, prob = TRUE, main = "Histogram with density curve")
line=mean(Muestras)
abline(v=line, col="red", lwd=4)
lines(density(Muestras), col = 2, lwd = 2)

install.packages("MASS")
install.packages("survival")
install.packages("fitdistrplus")


library(MASS)
library(survival)
library(fitdistrplus)

ajuste = fitdist(Muestras,"norm")

ajuste$estimate

##       mean         sd 
## 0.89800000 0.06577233

plot(ajuste)

n=30

require(MASS)
require(survival)
require(fitdistrplus)

Muestras=sapply(rep(30,500), Calcula_Estimador)
hist(Muestras, prob = TRUE, main = "Histogram with density curve")
line=mean(Muestras)
abline(v=line, col="red", lwd=4)
lines(density(Muestras), col = 2, lwd = 2)

install.packages("MASS")
install.packages("survival")
install.packages("fitdistrplus")


library(MASS)
library(survival)
library(fitdistrplus)

ajuste = fitdist(Muestras,"norm")

ajuste$estimate

##       mean         sd 
## 0.89786667 0.05613381

plot(ajuste)

n=50

require(MASS)
require(survival)
require(fitdistrplus)

Muestras=sapply(rep(50,500), Calcula_Estimador)
hist(Muestras, prob = TRUE, main = "Histogram with density curve")
line=mean(Muestras)
abline(v=line, col="red", lwd=4)
lines(density(Muestras), col = 2, lwd = 2)

install.packages("MASS")
install.packages("survival")
install.packages("fitdistrplus")


library(MASS)
library(survival)
library(fitdistrplus)

ajuste = fitdist(Muestras,"norm")

ajuste$estimate

##       mean         sd 
## 0.89776000 0.04234835

plot(ajuste)

n=100

require(MASS)
require(survival)
require(fitdistrplus)

Muestras=sapply(rep(100,500), Calcula_Estimador)
hist(Muestras, prob = TRUE, main = "Histogram with density curve")
line=mean(Muestras)
abline(v=line, col="red", lwd=4)
lines(density(Muestras), col = 2, lwd = 2)

install.packages("MASS")
install.packages("survival")
install.packages("fitdistrplus")


library(MASS)
library(survival)
library(fitdistrplus)

ajuste = fitdist(Muestras,"norm")

ajuste$estimate

##       mean         sd 
## 0.90102000 0.02871863

plot(ajuste)

n=200

require(MASS)
require(survival)
require(fitdistrplus)

Muestras=sapply(rep(200,500), Calcula_Estimador)
hist(Muestras, prob = TRUE, main = "Histogram with density curve")
line=mean(Muestras)
abline(v=line, col="red", lwd=4)
lines(density(Muestras), col = 2, lwd = 2)

install.packages("MASS")
install.packages("survival")
install.packages("fitdistrplus")


library(MASS)
library(survival)
library(fitdistrplus)

ajuste = fitdist(Muestras,"norm")

ajuste$estimate

##       mean         sd 
## 0.89896000 0.01942469

plot(ajuste)

n=500

require(MASS)
require(survival)
require(fitdistrplus)

Muestras=sapply(rep(500,500), Calcula_Estimador)
hist(Muestras, prob = TRUE, main = "Histogram with density curve")
line=mean(Muestras)
abline(v=line, col="red", lwd=4)
lines(density(Muestras), col = 2, lwd = 2)

install.packages("MASS")
install.packages("survival")
install.packages("fitdistrplus")


library(MASS)
library(survival)
library(fitdistrplus)

ajuste = fitdist(Muestras,"norm")

ajuste$estimate

##        mean          sd 
## 0.899668000 0.008947501

plot(ajuste)

Conclusión general:

Siempre el promedio está oscilando entre el valor real del parámetro, para este ejercicio fueron: 50%, 10% y 90%
Entre más iteraciones se realicen, mucho más precisa van a ser los estimadores y los datos se realionan más con un comportamiento normal.
Después de la iteración 30, se observa que la variación permanece constante, quiere decir que después de la iteraciones 30 el modelo se ajusta a la normal.

y b.

## Los 1 = Enfermos
## Los 0 = Sanos

lote_1 = c(rep(1, 100), rep(0,900))
lote_2 = c(rep(1, 150), rep(0,1350))

Calcula_Estimador_2=function(n){
  
  muestra1=sample(lote_1,size=n)
  P1=sum(muestra1)/n
  
  muestra2=sample(lote_2,size=n)
  P2=sum(muestra2)/n
  
  Diferencia=P1-P2
  
  return(Diferencia)
  
  
}

Calcula_Estimador_2(n=200)

## [1] 0.02

Muestras_2=sapply(rep(500,200), Calcula_Estimador_2)

hist(Muestras_2, prob = TRUE, main = "500 Iteraciones")
line=mean(Muestras_2)
abline(v=line, col="red", lwd=4)
lines(density(Muestras_2), col = 2, lwd = 2)

La diferencia de p1 y p2 son simetricos, tiene un comportamiento normal después de 500 iteraciones y su diferencia siempre en cero o muy cercano.

n=5

n=5

Muestras_2=sapply(rep(5,200), Calcula_Estimador_2)

hist(Muestras_2, prob = TRUE, main = "5 Iteraciones")
line=mean(Muestras_2)
abline(v=line, col="red", lwd=4)
lines(density(Muestras_2), col = 2, lwd = 2)

n=10

Muestras_2=sapply(rep(10,200), Calcula_Estimador_2)

hist(Muestras_2, prob = TRUE, main = "10 Iteraciones")
line=mean(Muestras_2)
abline(v=line, col="red", lwd=4)
lines(density(Muestras_2), col = 2, lwd = 2)

n=15

Muestras_2=sapply(rep(15,200), Calcula_Estimador_2)

hist(Muestras_2, prob = TRUE, main = "15 Iteraciones")
line=mean(Muestras_2)
abline(v=line, col="red", lwd=4)
lines(density(Muestras_2), col = 2, lwd = 2)

n=20

Muestras_2=sapply(rep(20,200), Calcula_Estimador_2)

hist(Muestras_2, prob = TRUE, main = "20 Iteraciones")
line=mean(Muestras_2)
abline(v=line, col="red", lwd=4)
lines(density(Muestras_2), col = 2, lwd = 2)

n=30

Muestras_2=sapply(rep(30,200), Calcula_Estimador_2)

hist(Muestras_2, prob = TRUE, main = "30 Iteraciones")
line=mean(Muestras_2)
abline(v=line, col="red", lwd=4)
lines(density(Muestras_2), col = 2, lwd = 2)

n=50

Muestras_2=sapply(rep(50,200), Calcula_Estimador_2)

hist(Muestras_2, prob = TRUE, main = "50 Iteraciones")
line=mean(Muestras_2)
abline(v=line, col="red", lwd=4)
lines(density(Muestras_2), col = 2, lwd = 2)

n=100

Muestras_2=sapply(rep(100,200), Calcula_Estimador_2)

hist(Muestras_2, prob = TRUE, main = "100 Iteraciones")
line=mean(Muestras_2)
abline(v=line, col="red", lwd=4)
lines(density(Muestras_2), col = 2, lwd = 2)

n=200

Muestras_2=sapply(rep(200,200), Calcula_Estimador_2)

hist(Muestras_2, prob = TRUE, main = "200 Iteraciones")
line=mean(Muestras_2)
abline(v=line, col="red", lwd=4)
lines(density(Muestras_2), col = 2, lwd = 2)

n=500

Muestras_2=sapply(rep(500,200), Calcula_Estimador_2)

hist(Muestras_2, prob = TRUE, main = "500 Iteraciones")
line=mean(Muestras_2)
abline(v=line, col="red", lwd=4)
lines(density(Muestras_2), col = 2, lwd = 2)

Conclusión:

A medida que las muestras incrementan de acuerdo a las iteraciones, el grafico de normalidad se va ajustando y después de 100 iteraciones ya se queda simétrico a la distribución normal con su media en 0.

Muestras_2=sapply(rep(500,200), Calcula_Estimador_2)

hist(Muestras_2, prob = TRUE, main = "500 Iteraciones")
line=mean(Muestras_2)
abline(v=line, col="red", lwd=4)
lines(density(Muestras_2), col = 2, lwd = 2)

La diferencia de p1 y p2 son simetricos, tiene un comportamiento normal después de 500 iteraciones y su diferencia siempre en cero o muy cercano.

## Los 1 = Enfermos
## Los 0 = Sanos

lote_1 = c(rep(1, 100), rep(0,900))
lote_2 = c(rep(1, 225), rep(0,1275))

Calcula_Estimador_2=function(n){
  
  muestra1=sample(lote_1,size=n)
  P1=sum(muestra1)/n
  
  muestra2=sample(lote_2,size=n)
  P2=sum(muestra2)/n
  
  Diferencia=P1-P2
  
  return(Diferencia)
  
  
}

Calcula_Estimador_2(n=200)

## [1] 0.005

n=5

Muestras_2=sapply(rep(5,200), Calcula_Estimador_2)

hist(Muestras_2, prob = TRUE, main = "5 Iteraciones")
line=mean(Muestras_2)
abline(v=line, col="red", lwd=4)
lines(density(Muestras_2), col = 2, lwd = 2)

n=10

Muestras_2=sapply(rep(10,200), Calcula_Estimador_2)

hist(Muestras_2, prob = TRUE, main = "10 Iteraciones")
line=mean(Muestras_2)
abline(v=line, col="red", lwd=4)
lines(density(Muestras_2), col = 2, lwd = 2)

n=15

Muestras_2=sapply(rep(15,200), Calcula_Estimador_2)

hist(Muestras_2, prob = TRUE, main = "15 Iteraciones")
line=mean(Muestras_2)
abline(v=line, col="red", lwd=4)
lines(density(Muestras_2), col = 2, lwd = 2)

n=20

Muestras_2=sapply(rep(20,200), Calcula_Estimador_2)

hist(Muestras_2, prob = TRUE, main = "20 Iteraciones")
line=mean(Muestras_2)
abline(v=line, col="red", lwd=4)
lines(density(Muestras_2), col = 2, lwd = 2)

n=30

Muestras_2=sapply(rep(30,200), Calcula_Estimador_2)

hist(Muestras_2, prob = TRUE, main = "30 Iteraciones")
line=mean(Muestras_2)
abline(v=line, col="red", lwd=4)
lines(density(Muestras_2), col = 2, lwd = 2)

n=50

Muestras_2=sapply(rep(50,200), Calcula_Estimador_2)

hist(Muestras_2, prob = TRUE, main = "50 Iteraciones")
line=mean(Muestras_2)
abline(v=line, col="red", lwd=4)
lines(density(Muestras_2), col = 2, lwd = 2)

n=100

Muestras_2=sapply(rep(100,200), Calcula_Estimador_2)

hist(Muestras_2, prob = TRUE, main = "100 Iteraciones")
line=mean(Muestras_2)
abline(v=line, col="red", lwd=4)
lines(density(Muestras_2), col = 2, lwd = 2)

n=200

Muestras_2=sapply(rep(200,200), Calcula_Estimador_2)

hist(Muestras_2, prob = TRUE, main = "200 Iteraciones")
line=mean(Muestras_2)
abline(v=line, col="red", lwd=4)
lines(density(Muestras_2), col = 2, lwd = 2)

n=500

Muestras_2=sapply(rep(500,200), Calcula_Estimador_2)

hist(Muestras_2, prob = TRUE, main = "500 Iteraciones")
line=mean(Muestras_2)
abline(v=line, col="red", lwd=4)
lines(density(Muestras_2), col = 2, lwd = 2)

Conclusión:

Al cambiar el % de simetría, la diferencia de p, va oscilando entre -0.1 y 0 hasta 50 iteraciones, de 100 en adelante la diferencia pierde su simetría y arroja diferencias muy variables.

Taller Simulación Inferencia

Santiago Garcia Cifuentes

19 Marzo 2022