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 T1: ",round(var(W$T2), digits = 3)," Eficiente")
}else{paste("Varianza T1: ",round(var(W$T2), digits = 3)," Ineficiente")}
)
print(
if(var(W$T3) <1.3){
paste("Varianza T1: ",round(var(W$T3), digits = 3)," Eficiente")
}else{paste("Varianza T1: ",round(var(W$T3), digits = 3)," Ineficiente")}
)
print(
if(var(W$T4) <1.3){
paste("Varianza T1: ",round(var(W$T4), digits = 3)," Eficiente")
}else{paste("Varianza T1: ",round(var(W$T4), digits = 3)," Ineficiente")}
)
}
funcion(20)
## T1 T2 T3 T4
## Min. :0.5265 Min. :0.9375 Min. :0.638 Min. :0.6898
## 1st Qu.:1.1246 1st Qu.:2.4047 1st Qu.:1.147 1st Qu.:1.2935
## Median :1.7188 Median :3.3108 Median :1.876 Median :2.4182
## Mean :1.9260 Mean :3.7150 Mean :2.041 Mean :2.2740
## 3rd Qu.:2.4387 3rd Qu.:4.4976 3rd Qu.:2.620 3rd Qu.:3.1369
## Max. :4.2273 Max. :8.1654 Max. :4.012 Max. :4.1196
## [1] "-------------------------insesgadez-----------------------------------"
## [1] "Promedio T1: 1.926 Insesgado: cerca a la media"
## [1] "Promedio T2: 3.715 Sesgado: lejos de la media"
## [1] "Promedio T3: 2.041 Insesgado: cerca a la media"
## [1] "Promedio T4: 2.274 Insesgado: cerca a la media"
## [1] "-------------------------Varianzas-----------------------------------"
## [1] "Varianza T1: 1.032 Eficiente"
## [1] "Varianza T1: 3.852 Ineficiente"
## [1] "Varianza T1: 1.143 Eficiente"
## [1] "Varianza T1: 1.252 Eficiente"funcion(50)
## T1 T2 T3 T4
## Min. :0.3126 Min. :0.6135 Min. :0.3192 Min. :0.2994
## 1st Qu.:1.2751 1st Qu.:2.6143 1st Qu.:1.2987 1st Qu.:1.3824
## Median :1.9265 Median :3.8724 Median :1.8882 Median :2.1519
## Mean :2.0467 Mean :4.0993 Mean :1.9861 Mean :2.2743
## 3rd Qu.:2.6344 3rd Qu.:5.3018 3rd Qu.:2.6034 3rd Qu.:2.9724
## Max. :4.5999 Max. :9.4644 Max. :4.8106 Max. :5.3245
## [1] "-------------------------insesgadez-----------------------------------"
## [1] "Promedio T1: 2.047 Insesgado: cerca a la media"
## [1] "Promedio T2: 4.099 Sesgado: lejos de la media"
## [1] "Promedio T3: 1.986 Insesgado: cerca a la media"
## [1] "Promedio T4: 2.274 Insesgado: cerca a la media"
## [1] "-------------------------Varianzas-----------------------------------"
## [1] "Varianza T1: 1.011 Eficiente"
## [1] "Varianza T1: 4.197 Ineficiente"
## [1] "Varianza T1: 0.974 Eficiente"
## [1] "Varianza T1: 1.355 Ineficiente"funcion(100)
## T1 T2 T3 T4
## Min. :0.2978 Min. : 0.5432 Min. :0.3151 Min. :0.367
## 1st Qu.:1.4098 1st Qu.: 2.7743 1st Qu.:1.4012 1st Qu.:1.485
## Median :1.9001 Median : 3.8187 Median :1.8164 Median :2.080
## Mean :2.0380 Mean : 4.0684 Mean :1.9700 Mean :2.356
## 3rd Qu.:2.5370 3rd Qu.: 5.0326 3rd Qu.:2.6094 3rd Qu.:2.986
## Max. :5.3316 Max. :11.7506 Max. :4.5622 Max. :6.729
## [1] "-------------------------insesgadez-----------------------------------"
## [1] "Promedio T1: 2.038 Insesgado: cerca a la media"
## [1] "Promedio T2: 4.068 Sesgado: lejos de la media"
## [1] "Promedio T3: 1.97 Insesgado: cerca a la media"
## [1] "Promedio T4: 2.356 Sesgado: lejos de la media"
## [1] "-------------------------Varianzas-----------------------------------"
## [1] "Varianza T1: 0.906 Eficiente"
## [1] "Varianza T1: 4.034 Ineficiente"
## [1] "Varianza T1: 0.748 Eficiente"
## [1] "Varianza T1: 1.522 Ineficiente"funcion(1000)
## T1 T2 T3 T4
## Min. :0.1968 Min. : 0.384 Min. :0.2631 Min. :0.2689
## 1st Qu.:1.2557 1st Qu.: 2.430 1st Qu.:1.3384 1st Qu.:1.4585
## Median :1.7860 Median : 3.538 Median :1.8675 Median :2.1209
## Mean :1.9753 Mean : 3.946 Mean :1.9920 Mean :2.3401
## 3rd Qu.:2.5112 3rd Qu.: 5.013 3rd Qu.:2.5109 3rd Qu.:2.9633
## Max. :5.9694 Max. :13.717 Max. :5.9838 Max. :7.8954
## [1] "-------------------------insesgadez-----------------------------------"
## [1] "Promedio T1: 1.975 Insesgado: cerca a la media"
## [1] "Promedio T2: 3.946 Sesgado: lejos de la media"
## [1] "Promedio T3: 1.992 Insesgado: cerca a la media"
## [1] "Promedio T4: 2.34 Sesgado: lejos de la media"
## [1] "-------------------------Varianzas-----------------------------------"
## [1] "Varianza T1: 0.973 Eficiente"
## [1] "Varianza T1: 4.387 Ineficiente"
## [1] "Varianza T1: 0.854 Eficiente"
## [1] "Varianza T1: 1.461 Ineficiente"
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