Actividad 2 - Metodos y simulacion estadistica
Aldair Blanco
2023-09-09
Problema 1
• 1 Estimación del valor de π.
m <- as.integer(c(1000,10000,100000))
for (variable in m) {
n <- variable
x <- runif(n,0,1)
y <- runif(n,0,1)
xy <- data.frame(x,y)
xy <- as.data.frame(xy)
fx = function(x,y){
sqrt(
(x-0.5)^2+(y-0.5)^2
)
}
xy$dist <- mapply(fx,xy$x,xy$y)
puntos_dentro <- xy[xy$dist <= 0.5,]
rad <- 0.5 # Valor del radio
xcenter <- 0.5 # Coordenada en x del centro
ycenter <- 0.5 # Coordenada en y del centro
numero_puntos_dentro <- sum(xy$dist < 0.5)
pii <- numero_puntos_dentro / n
estimar_pi <- 4 * pii
par(mfrow = c(1, 2))
plot(c(0, 1), c(0, 1), type="n", xlab="X", ylab="Y", cex.main=0.7, cex.sub=0.5,
sub=paste0("Para una Muestra: ",n),
main = paste0("Estimación de π: ", round(estimar_pi, digits = 4),
" Puntos dentro del circulo: ",numero_puntos_dentro))
theta <- seq(0, 2 * pi, length = 200)
polygon(x=rad * cos(theta) + xcenter,
y=rad * sin(theta) + ycenter,
lwd=3, col='steelblue1', border='steelblue4')
abline(v = seq(0, 1, 0.1), lty = 2, col = "gray50")
abline(h = seq(0, 1, 0.1), lty = 2, col = "gray50")
points(puntos_dentro$x, puntos_dentro$y, col="red")
puntos_fuera <- xy[xy$dist > 0.5,]
plot(c(0, 1), c(0, 1), type="n", xlab="X", ylab="Y", cex.main=0.7, cex.sub=0.5,
sub=paste0("Para una Muestra: ",n),
main = paste0(" Puntos fuera del circulo: ",
length(puntos_fuera$x)))
theta <- seq(0, 2 * pi, length = 200)
polygon(x=rad * cos(theta) + xcenter,
y=rad * sin(theta) + ycenter,
lwd=3, col='steelblue1', border='steelblue4')
abline(v = seq(0, 1, 0.1), lty = 2, col = "gray50")
abline(h = seq(0, 1, 0.1), lty = 2, col = "gray50")
points(puntos_fuera$x, puntos_fuera$y, col="black")
}
Problema 2
• 2 Propiedades de los estimadores.
m <- c(20,50,100,1000)
#par(mfrow = c(1, 1))
funcion <- function(n) {
n <- n
lambda <- 0.5
x1 <- rexp(n,lambda)
x2 <- rexp(n,lambda)
x3 <- rexp(n,lambda)
x4 <- rexp(n,lambda)
base <- data.frame(x1,x2,x3,x4)
fx1 <- function(x1,x2,x3,x4){
((x1+x2)/6)+((x3+x4)/3)
}
fx2 <- function(x1,x2,x3,x4){
(x1+(2*x2)+(3*x3)+(4*x4))/5
}
fx3 <- function(x1,x2,x3,x4){
(x1+x2+x3+x4)/4
}
fx4 <- function(x1,x2,x3,x4){
(min(x1,x2,x3,x4)+max(x1,x2,x3,x4))/2
}
#fx1(base$x1[1],base$x2[1],base$x3[1],base$x4[1])
T1 <- mapply(fx1,base$x1,base$x2,base$x3,base$x4)
T2 <- mapply(fx2,base$x1,base$x2,base$x3,base$x4)
T3 <- mapply(fx3,base$x1,base$x2,base$x3,base$x4)
T4 <- mapply(fx4,base$x1,base$x2,base$x3,base$x4)
W <- data.frame(T1,T2,T3,T4)
boxplot(W, main=paste0("Numero de muestras n:", n," Lambda:" ,lambda))
abline(h=2, col="red3")
print(
summary(W),
)
print(
"-------------------------insesgadez-----------------------------------"
)
print(
if(mean(W$T1) >1.7 && mean(W$T1) <2.3){
paste("Promedio T1: ",round(mean(W$T1), digits = 3),"Insesgado: cerca a la media")
}else{paste("Promedio T1: ",round(mean(W$T1), digits = 3),"Sesgado: lejos de la media")}
)
print(
if(mean(W$T2) >1.7 && mean(W$T2) <2.3){
paste("Promedio T2: ",round(mean(W$T2), digits = 3),"Insesgado: cerca a la media")
}else{paste("Promedio T2: ",round(mean(W$T2), digits = 3),"Sesgado: lejos de la media")}
)
print(
if(mean(W$T3) >1.7 && mean(W$T3) <2.3){
paste("Promedio T3: ",round(mean(W$T3), digits = 3),"Insesgado: cerca a la media")
}else{paste("Promedio T3: ",round(mean(W$T3), digits = 3),"Sesgado: lejos de la media")}
)
print(
if(mean(W$T4) >1.7 && mean(W$T4) <2.3){
paste("Promedio T4: ",round(mean(W$T4), digits = 3),"Insesgado: cerca a la media")
}else{paste("Promedio T4: ",round(mean(W$T4), digits = 3),"Sesgado: lejos de la media")}
)
print(
"-------------------------Varianzas-----------------------------------"
)
print(
if(var(W$T1) <1.3){
paste("Varianza T1: ",round(var(W$T1), digits = 3)," Eficiente")
}else{paste("Varianza T1: ",round(var(W$T1), digits = 3)," Ineficiente")}
)
print(
if(var(W$T2) <1.3){
paste("Varianza T2: ",round(var(W$T2), digits = 3)," Eficiente")
}else{paste("Varianza T2: ",round(var(W$T2), digits = 3)," Ineficiente")}
)
print(
if(var(W$T3) <1.3){
paste("Varianza T3: ",round(var(W$T3), digits = 3)," Eficiente")
}else{paste("Varianza T3: ",round(var(W$T3), digits = 3)," Ineficiente")}
)
print(
if(var(W$T4) <1.3){
paste("Varianza T4: ",round(var(W$T4), digits = 3)," Eficiente")
}else{paste("Varianza T4: ",round(var(W$T4), digits = 3)," Ineficiente")}
)
}
funcion(20)
## T1 T2 T3 T4
## Min. :0.2065 Min. : 0.441 Min. :0.194 Min. :0.174
## 1st Qu.:1.1174 1st Qu.: 2.351 1st Qu.:1.085 1st Qu.:1.141
## Median :1.6946 Median : 3.047 Median :1.686 Median :1.983
## Mean :2.0717 Mean : 4.089 Mean :2.112 Mean :2.407
## 3rd Qu.:2.4928 3rd Qu.: 5.031 3rd Qu.:2.503 3rd Qu.:3.383
## Max. :7.3305 Max. :14.805 Max. :6.688 Max. :6.970
## [1] "-------------------------insesgadez-----------------------------------"
## [1] "Promedio T1: 2.072 Insesgado: cerca a la media"
## [1] "Promedio T2: 4.089 Sesgado: lejos de la media"
## [1] "Promedio T3: 2.112 Insesgado: cerca a la media"
## [1] "Promedio T4: 2.407 Sesgado: lejos de la media"
## [1] "-------------------------Varianzas-----------------------------------"
## [1] "Varianza T1: 2.521 Ineficiente"
## [1] "Varianza T2: 10.425 Ineficiente"
## [1] "Varianza T3: 2.307 Ineficiente"
## [1] "Varianza T4: 2.654 Ineficiente"funcion(50)
## T1 T2 T3 T4
## Min. :0.515 Min. :0.9771 Min. :0.4368 Min. :0.5472
## 1st Qu.:1.115 1st Qu.:2.1957 1st Qu.:1.1983 1st Qu.:1.3582
## Median :1.551 Median :3.0903 Median :1.6015 Median :2.0228
## Mean :1.725 Mean :3.4121 Mean :1.7939 Mean :2.1609
## 3rd Qu.:2.151 3rd Qu.:4.1707 3rd Qu.:2.3495 3rd Qu.:2.6918
## Max. :3.633 Max. :7.5797 Max. :3.8368 Max. :4.6993
## [1] "-------------------------insesgadez-----------------------------------"
## [1] "Promedio T1: 1.725 Insesgado: cerca a la media"
## [1] "Promedio T2: 3.412 Sesgado: lejos de la media"
## [1] "Promedio T3: 1.794 Insesgado: cerca a la media"
## [1] "Promedio T4: 2.161 Insesgado: cerca a la media"
## [1] "-------------------------Varianzas-----------------------------------"
## [1] "Varianza T1: 0.601 Eficiente"
## [1] "Varianza T2: 2.563 Ineficiente"
## [1] "Varianza T3: 0.69 Eficiente"
## [1] "Varianza T4: 1.06 Eficiente"funcion(100)
## T1 T2 T3 T4
## Min. :0.4473 Min. : 0.8518 Min. :0.650 Min. :0.7189
## 1st Qu.:1.3454 1st Qu.: 2.4906 1st Qu.:1.264 1st Qu.:1.4729
## Median :1.8181 Median : 3.4682 Median :1.913 Median :2.0376
## Mean :2.0582 Mean : 4.1153 Mean :2.048 Mean :2.4010
## 3rd Qu.:2.6402 3rd Qu.: 5.5751 3rd Qu.:2.544 3rd Qu.:3.0650
## Max. :5.1782 Max. :10.6002 Max. :4.398 Max. :6.5441
## [1] "-------------------------insesgadez-----------------------------------"
## [1] "Promedio T1: 2.058 Insesgado: cerca a la media"
## [1] "Promedio T2: 4.115 Sesgado: lejos de la media"
## [1] "Promedio T3: 2.048 Insesgado: cerca a la media"
## [1] "Promedio T4: 2.401 Sesgado: lejos de la media"
## [1] "-------------------------Varianzas-----------------------------------"
## [1] "Varianza T1: 1.025 Eficiente"
## [1] "Varianza T2: 4.496 Ineficiente"
## [1] "Varianza T3: 0.875 Eficiente"
## [1] "Varianza T4: 1.552 Ineficiente"funcion(1000)
## T1 T2 T3 T4
## Min. :0.2061 Min. : 0.4684 Min. :0.2193 Min. :0.2085
## 1st Qu.:1.2244 1st Qu.: 2.3946 1st Qu.:1.2673 1st Qu.:1.4557
## Median :1.8138 Median : 3.6481 Median :1.8586 Median :2.1333
## Mean :2.0265 Mean : 4.0532 Mean :2.0161 Mean :2.3811
## 3rd Qu.:2.6094 3rd Qu.: 5.2548 3rd Qu.:2.6440 3rd Qu.:3.0397
## Max. :7.2557 Max. :15.9063 Max. :6.8337 Max. :8.3720
## [1] "-------------------------insesgadez-----------------------------------"
## [1] "Promedio T1: 2.027 Insesgado: cerca a la media"
## [1] "Promedio T2: 4.053 Sesgado: lejos de la media"
## [1] "Promedio T3: 2.016 Insesgado: cerca a la media"
## [1] "Promedio T4: 2.381 Sesgado: lejos de la media"
## [1] "-------------------------Varianzas-----------------------------------"
## [1] "Varianza T1: 1.128 Eficiente"
## [1] "Varianza T2: 4.889 Ineficiente"
## [1] "Varianza T3: 0.999 Eficiente"
## [1] "Varianza T4: 1.639 Ineficiente"El mejor estimador de θ^ es θ^1 debido a que desde la primera muestra de n=20 su promedio se mantiene cerca de 2 y tiene poca varianza, claro está no dejando al lado el estimador θ^3 que también tuvo un promedio muy cercano a 2 pero cuando la población es n > 40.
Problema 3
• 3 Teorema del Límite Central.
library(moments)
library(nortest)
TeoremaLC <- function(n,porc,m,q){
n<-n
porc <- porc
#a) aun que el punto a me dice que cree una poblacion 1.000, le doy la opcion
#al usuario de crear la poblacion del numero que el quiera para comprobar
#el teorema del limite central, de igaul forma para el porcentaje (porc)
#de plantas enfermas
poblacion <- rbinom(n, 1, porc)
#b) Funcion que calcula el pˆ para n=x
estimador <- function(muestras) {
muestra <- sample(poblacion, muestras, replace = FALSE)
estimar <- mean(muestra)
return(estimar)
}
##Funcion tlc que genera una matrix por el numero de poblacion con
##Respecto al numero de muestras generadas y llama la funcion anterior estimar
##Para estimar la muestra de los resultados, por ultimo guarda los resultados de
##cada muestra y grafica los histograma por cada muestra con su varianza, sesgo, y
##coeficiente de variación de Skewness
tlc <-function(n,m){
resultados <- matrix(nrow = n, ncol = length(m))
for (i in 1:n) {
for (j in 1:length(m)) {
estimadorAux <- estimador(m[j])
resultados[i, j] <- estimadorAux
}
}
resultados <- as.data.frame(resultados)
if (length(m)==1) {
columnas<-1
}else{columnas<-2}
par(cex=0.5, cex.axis=.5,
cex.lab=.5, cex.main=.5,
cex.sub=.5, mfrow=c(as.integer(ceiling(length(m)/2)),columnas),
mai = c(.5, .5, .5, .5))
for (variable in 1:length(m)) {
media <- mean(resultados[,variable])
ds <- sd(resultados[,variable])
varianza <- var(resultados[,variable])
CoeSkewness <- skewness(resultados[,variable])
#CoeFisher <- kurtosis(resultados[,variable])
#if (m[variable]<51) {
prueba <- shapiro.test(resultados[,variable])
nombre <- "Shapiro"
#}else{
#pruebaks <- ks.test(resultados[,variable])
# nombreks <- "Kolmogorov"
#}
prueba <- as.numeric(prueba$p.value)
sesgo <- mean(resultados[,variable]) - porc
varianza <- signif(var(resultados[,variable]), 3)
sesgo <- signif(mean(resultados[,variable]) - porc, 3)
CoeSkewness <- signif(CoeSkewness, 3)
h <-hist(resultados[,variable],
main = paste0("Muestra:",m[variable],"\n",
"Varianza: ", varianza," | ",
"Sesgo: ", sesgo,"\n",
"Asimetria Skewness:", CoeSkewness," | ",
nombre," Test:",signif(prueba, 3)
),
col="darkseagreen2"
)
text(h$mids, h$counts, labels=h$counts, adj=c(0.5,0.25),col="#3b5998")
abline(v = media, col="red", lwd=3, lty=2)
if (q==1) {
qqnorm(resultados[,variable],
main =paste0("Muestra:",m[variable],"\n",
"Varianza: ", varianza," | ",
"Sesgo: ", sesgo,"\n",
"Asimetria Skewness:", CoeSkewness," | ",
nombre," Test:",signif(prueba, 3)
))
qqline(resultados[,variable],col="red")
}
}
}
return(tlc(n,m))
}
#Poblacion: tamaño de poblacion
#%Enferma: Porcentaje de poblacion enferma
#Numero de muestra: Numero de las muestras < poblacion
#q: 1:si quieren el grafico qqnorm 0:No
#TeoremaLC(Poblacion,%Enferma,Numero de muestra,qqnorm)
TeoremaLC(1000,0.5,c(500),0)
Comentaros Para una muestra 500.
Según los resultados obtenidos cuando se ejecuto la función TeoremaLC con los parametros Población 1000, plantas enfermas 0.5(50%) y numero de muestras 500 (TeoremaLC(1000,0.5,c(500))), esto nos genera una variabilidad en los datos casi nula, con un resultado de varianza < 0.005. Genera muy poco sesgo > -0.005 y < 0.005, Cuando le aplicamos el coeficiente de asimetría de Skewness nos da como resultado que la distribución es negativamente sesgada, esto ultimo indica que la cola de la distribución es mas pesada del lado izquierdo que del derecho.
TeoremaLC(1000,0.5,c(5,10,15),1)
TeoremaLC(1000,0.5,c(20,30,50),1)
TeoremaLC(1000,0.5,c(60,100,200),1)
TeoremaLC(1000,0.5,c(500),1)
Comentarios Para una muestra c(5,10,15,20,30,50,60,100,200,500)
Observamos que a medida que va creciendo el numero de muestras no importando que sea aleatorio los datos van tomando una distribución normal esto comprueba el teorema del limite central (TLC) el cual dice que entre mas grande sea la muestra los datos tomaran forma de una distribución normal, vemos de igual forma la prueba de shapiro wilk (shapiro.test) que para las dos muestra que hicimos de 500 no hay evidencia suficiente para rechazar la hipótesis nula (H0) y que los datos siguen una distribución normal.
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.
Para una poblacion con plantas enfermas del 10%
TeoremaLC(1000,0.1,c(5,10,15),1)
TeoremaLC(1000,0.1,c(20,30,50),1)
TeoremaLC(1000,0.1,c(60,100,200),1)
TeoremaLC(1000,0.1,c(500),1)
Para una poblacion con plantas enfermas del 90%
TeoremaLC(1000,0.9,c(5,10,15),1)
TeoremaLC(1000,0.9,c(20,30,50),1)
TeoremaLC(1000,0.9,c(60,100,200),1)
TeoremaLC(1000,0.9,c(500),1)
Para una muestra con plantas del 10% observamos como los datos a medida que va creciendo el numero de la muestra van alcanzando una distribución normal comenzando de izquierda a derecha y para poblacion de plantas enfermas de 90% observamos como sucede lo contrario a medida que crece la muestra los datos tomaran una distribución normal de derecha a izquierda.
Problema 4
• 4 Estimacción boostrap.
library(ggplot2,warn.conflicts = FALSE)## Warning: package 'ggplot2' was built under R version 4.2.3#Semilla para que los datos den igual
set.seed(1)
#Muestra
DatosMP <- c(7.69, 4.97, 4.56, 6.49, 4.34, 6.24, 4.45)
#Muestra con reemplazo de 70000
Muestras <- sample(DatosMP,7000, replace = TRUE)
#DataFrame de datos de 1000*7
DataFrame<-data.frame(matrix(Muestras,nrow=1000,ncol=7))
#Media por cada fila
DataFrame$Media <- apply(DataFrame,1,mean)
#Calculo de metodo 1
IntervaloMetodo1=quantile(DataFrame$Media, probs=c(0.025, 0.975))
#Calculo de metodo 2
IntervaloMetodo2<-c(2*mean(DataFrame$Media)-IntervaloMetodo1[2],
2*mean(DataFrame$Media)-IntervaloMetodo1[1])
IntervaloMetodo1## 2.5% 97.5%
## 4.747000 6.418964IntervaloMetodo2## 97.5% 2.5%
## 4.655576 6.327540##Grafico
ggplot()+
geom_histogram(aes(x=DataFrame$Media),
bins = 12, color="#034A94", fill="#F4CCCC")+
labs(fill="",x="",y="")+
geom_vline(xintercept = IntervaloMetodo1,
linewidth=1, color="#FFFF00")+
geom_vline(xintercept = IntervaloMetodo2,
linewidth=1, color="#6EF68A")
##Medias para estimacion de sesgo
media_original <- mean(DatosMP)
medias_bootstrap <- mean(DataFrame$Media)
sesgo <- mean(medias_bootstrap - media_original)
sesgo <- signif(sesgo,3)
print(paste("Sesgo de la media:", sesgo))## [1] "Sesgo de la media: 0.00298"Resultados:
Debido al ultimo calculo de sesgo 0.00298 el cual es muy bajo y los intervalos de confianzas 2.5% 4.747 97.5% 6.419, podriamos decir que si confiariamos en las estimaciones del metodo de boostrap