Actividad 2
Guillén León
2024-09-11
Ejercicio 1
Estimacion del valor \(\pi\)
\(\pi\)
n<-1000
x<-runif(n)
y<-runif(n)
d<-(x-0.5)^2+(y-0.5)^2
d2<-as.numeric(d<=0.25)
epi<-(sum(d2)/n)*4 #estimación de pi
error<-pi-epi
p1<-sum((d2)/n)
cat("Proporción de puntos dentro del circulo:",p1,"\n")Proporción de puntos dentro del circulo: 0.803 cat("Estimación de pi:",epi,"\n")Estimación de pi: 3.212 cat("Error de estimación:",error)Error de estimación: -0.07040735n<-10000
x<-runif(n)
y<-runif(n)
d<-(x-0.5)^2+(y-0.5)^2
d2<-as.numeric(d<=0.25)
epi<-(sum(d2)/n)*4 #estimación de pi
error<-pi-epi
p1<-sum((d2)/n)
cat("Proporción de puntos dentro del circulo:",p1,"\n")Proporción de puntos dentro del circulo: 0.7864 cat("Estimación de pi:",epi,"\n")Estimación de pi: 3.1456 cat("Error de estimación:",error)Error de estimación: -0.004007346n<-100000
x<-runif(n)
y<-runif(n)
d<-(x-0.5)^2+(y-0.5)^2
d2<-as.numeric(d<=0.25)
epi<-(sum(d2)/n)*4 #estimación de pi
error<-pi-epi
p1<-sum((d2)/n)
cat("Proporción de puntos dentro del circulo:",p1,"\n")Proporción de puntos dentro del circulo: 0.78669 cat("Estimación de pi:",epi,"\n")Estimación de pi: 3.14676 cat("Error de estimación:",error)Error de estimación: -0.005167346n<-1000000
x<-runif(n)
y<-runif(n)
d<-(x-0.5)^2+(y-0.5)^2
d2<-as.numeric(d<=0.25)
epi<-(sum(d2)/n)*4 #estimación de pi
error<-pi-epi
p1<-sum((d2)/n)
cat("Proporción de puntos dentro del circulo:",p1,"\n")Proporción de puntos dentro del circulo: 0.785247 cat("Estimación de pi:",epi,"\n")Estimación de pi: 3.140988 cat("Error de estimación:",error)Error de estimación: 0.0006046536Conclusión
Como se puede obervar en cuatro casos distintos de n: 1000, 10000, 100000 y 1000000; entre mayor sea la muestra, más se aproxima el valor estimado de π al valor real, cumpliendo con el principio de consistencia de un estimador.
Ejercicio 2
Propiedad de los Estimadores
n=20
Con una muestra de 20 datos, el mejor estimador es e3, dado que se aproxima a 2 y tiene la menor varianza.
set.seed(123)
n=20
l<-2
x1<-rexp(n,1/l)
x2<-rexp(n,1/l)
x3<-rexp(n,1/l)
x4<-rexp(n,1/l)
x1234<-data.frame(x1,x2,x3,x4)
minx<-apply(x1234,1,min)
maxx<-apply(x1234,1,max)
e1234<-data.frame(e1<-(x1+x2)/6+(x3+x4)/3,
e2<-(x1+2*x2+3*x3+4*x4)/5,
e3<-(x1+x2+x3+x4)/4,
e4<-(minx+maxx)/2)
boxplot(e1234, col=c("blue","green","yellow","red"))
abline(h=l,col="pink")
media<-apply(e1234, 2,mean)
varianza<-apply(e1234, 2,var)
sesgo= l-media
rbind(media,varianza) e1.....x1...x2..6....x3...x4..3
media 2.070849
varianza 1.088533
e2.....x1...2...x2...3...x3...4...x4..5 e3.....x1...x2...x3...x4..4
media 4.087918 2.0342578
varianza 3.507508 0.7943927
e4.....minx...maxx..2
media 2.401576
varianza 2.082907print.table(media, varianza, sesgo) e1.....x1...x2..6....x3...x4..3 e2.....x1...2...x2...3...x3...4...x4..5
"2" "4"
e3.....x1...x2...x3...x4..4 e4.....minx...maxx..2
"2" "2" #summary(e1234)n=50
Con una muestra de 50 datos, el mejor estimador es e1, dado que se aproxima más a 2 y tiene la menor varianza, aunque el estimador e3 tambien es un buen estimador pues es menos insesgado y eficiente.
set.seed(123)
n=50
l<-2
x1<-rexp(n,1/l)
x2<-rexp(n,1/l)
x3<-rexp(n,1/l)
x4<-rexp(n,1/l)
x1234<-data.frame(x1,x2,x3,x4)
minx<-apply(x1234,1,min)
maxx<-apply(x1234,1,max)
e1234<-data.frame(e1<-(x1+x2)/6+(x3+x4)/3,
e2<-(x1+2*x2+3*x3+4*x4)/5,
e3<-(x1+x2+x3+x4)/4,
e4<-(minx+maxx)/2)
boxplot(e1234, col=c("blue","green","yellow","red"))
abline(h=l,col="pink")
media<-apply(e1234, 2,mean)
varianza<-apply(e1234, 2,var)
sesgo= l-media
rbind(media,varianza) e1.....x1...x2..6....x3...x4..3
media 1.988809
varianza 1.071557
e2.....x1...2...x2...3...x3...4...x4..5 e3.....x1...x2...x3...x4..4
media 3.958039 2.014466
varianza 4.287298 1.075888
e4.....minx...maxx..2
media 2.315248
varianza 1.972405print.table(media, varianza, sesgo) e1.....x1...x2..6....x3...x4..3 e2.....x1...2...x2...3...x3...4...x4..5
"2" "4"
e3.....x1...x2...x3...x4..4 e4.....minx...maxx..2
"2" "2" n=100
Con una muestra de 100 datos, el mejor estimador es e3, dado que se aproxima más a 2 y tiene la menor varianza, aunque el estimador e1 tambien es un buen estimador pues es menos insesgado y eficiente
set.seed(123)
n=100
l<-2
x1<-rexp(n,1/l)
x2<-rexp(n,1/l)
x3<-rexp(n,1/l)
x4<-rexp(n,1/l)
x1234<-data.frame(x1,x2,x3,x4)
minx<-apply(x1234,1,min)
maxx<-apply(x1234,1,max)
e1234<-data.frame(e1<-(x1+x2)/6+(x3+x4)/3,
e2<-(x1+2*x2+3*x3+4*x4)/5,
e3<-(x1+x2+x3+x4)/4,
e4<-(minx+maxx)/2)
boxplot(e1234, col=c("blue","green","yellow","red"))
abline(h=l,col="pink")
media<-apply(e1234, 2,mean)
varianza<-apply(e1234, 2,var)
sesgo= l-media
rbind(media,varianza) e1.....x1...x2..6....x3...x4..3
media 1.9767021
varianza 0.8288984
e2.....x1...2...x2...3...x3...4...x4..5 e3.....x1...x2...x3...x4..4
media 3.913297 1.9861432
varianza 3.756859 0.7913008
e4.....minx...maxx..2
media 2.293568
varianza 1.316554print.table(media, varianza, sesgo) e1.....x1...x2..6....x3...x4..3 e2.....x1...2...x2...3...x3...4...x4..5
"2" "4"
e3.....x1...x2...x3...x4..4 e4.....minx...maxx..2
"2" "2" n=1000
Con una muestra de 1000 datos, el mejor estimador es e3, dado que se aproxima más a 2 y tiene la menor varianza.
n=1000
n=100
l<-2
x1<-rexp(n,1/l)
x2<-rexp(n,1/l)
x3<-rexp(n,1/l)
x4<-rexp(n,1/l)
x1234<-data.frame(x1,x2,x3,x4)
minx<-apply(x1234,1,min)
maxx<-apply(x1234,1,max)
e1234<-data.frame(e1<-(x1+x2)/6+(x3+x4)/3,
e2<-(x1+2*x2+3*x3+4*x4)/5,
e3<-(x1+x2+x3+x4)/4,
e4<-(minx+maxx)/2)
boxplot(e1234, col=c("blue","green","yellow","red"))
abline(h=l,col="pink")
media<-apply(e1234, 2,mean)
varianza<-apply(e1234, 2,var)
sesgo= l-media
rbind(media,varianza) e1.....x1...x2..6....x3...x4..3
media 2.130337
varianza 1.161209
e2.....x1...2...x2...3...x3...4...x4..5 e3.....x1...x2...x3...x4..4
media 4.279351 2.098959
varianza 5.491821 1.024917
e4.....minx...maxx..2
media 2.467476
varianza 1.874088print.table(media, varianza, sesgo) e1.....x1...x2..6....x3...x4..3 e2.....x1...2...x2...3...x3...4...x4..5
"2" "4"
e3.....x1...x2...x3...x4..4 e4.....minx...maxx..2
"2" "2" Conclusiones
Los estimadores e1, e3 y e4 son insesgados, eficientes y consistentes, ya que su media tiende a aproximarse al valor real del parámetro, su varianza tiende a disminuir y su sesgo tiende a cero a medida que el tamaño de la muestra aumenta. Por otro lado, el estimador e2 no es eficiente, puesto que a pesar de que su varianza disminuya, esta tiende a estabilizarse. Asi mismo, el estimador e2 no es concistente, ya que su media se distancia del valor real del parámetro a medida que el tamaño de la muestra aumenta, y como consecuencia generar un sesgo positivo pues el valor de la media tiende a estar por encima del valor real del parámetro.
La simulación con diferentes tamaños de muestra y el apoyo con graficos permite evidenciar la consistencia y potencia estadistica de las propiedades de los estimadores para tomar decisiones correctas en la selección de los estimadores.
Ejercicio 3
Teorema del Límite Central
A partir de la generación de 1000 (plantas) números aleatorios de una distrbución binomial, y un porcentaje de 50% de plantas enfermas, se realiza una simulación.
planta<- 1000
enfer<- rbinom(planta, size = 1, prob = 0.5)
data_planta <- data.frame(id=1:planta, tipo = enfer) # 1. planta enferma, 0: planta sanaAhora se calcula el estimador de la proporción muestral para un tamaño muestral de 20.
funmuestra <- function(data, n_muestra) {
muestra <- sample(data_planta$tipo, n_muestra, replace = FALSE)
p_estimado <- mean(muestra)
return(list(muestra = muestra, p_estimado = p_estimado))
}
## Se realiza un ejemplo para ver salida de la función:
ej <- funmuestra(data_planta, 20)
cat("La muestra extraída ", ej$muestra, " tiene un p estimado =", ej$p_estimado)La muestra extraída 0 1 1 0 1 0 0 0 0 1 1 1 1 1 1 0 1 1 1 0 tiene un p estimado = 0.6Se repite el ejercicio anterior con n=500 veces y analizan los resultados del estimador pˆ.
n_muestra <- 20
n_repeticiones <- 500
df_estimacion <- data.frame(iteracion = numeric(n_repeticiones), p_estimados = numeric(n_repeticiones))
for (i in 1:n_repeticiones) {
df_estimacion$iteracion[i] <- i
muestra_resultado <- funmuestra(data_planta, n_muestra)
df_estimacion$p_estimados[i] <- muestra_resultado$p_estimado
}Se realiza gráfico para evaluar la simetría de la distribución
library(ggplot2)
library(tidyverse)-- Attaching core tidyverse packages ------------------------ tidyverse 2.0.0 --
v dplyr 1.1.4 v readr 2.1.2
v forcats 1.0.0 v stringr 1.5.1
v lubridate 1.9.3 v tibble 3.2.1
v purrr 1.0.2 v tidyr 1.3.1
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag() masks stats::lag()
i Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(hrbrthemes)
ggplot(df_estimacion, aes(x=p_estimados)) +
geom_histogram(binwidth=0.05, fill="#d6a692", color="#e9ecef", alpha=0.9) +
ggtitle(bquote("Distribución de " ~ hat(p) ~ "")) +
labs(x=(bquote("" ~ hat(p) ~ "")), y="densidad") +
theme_ipsum() +
theme(
plot.title = element_text(size=15)
)Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family
not found in Windows font database
Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family
not found in Windows font database
Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family
not found in Windows font databaseWarning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font database
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font database
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font database
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font database
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font database
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font database
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font database
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font database
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font database
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font database
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font database
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font database
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font database
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font databaseWarning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Se calcula la simetría de la distribución
library(parameters)
simetria <- skewness(df_estimacion$p_estimados)
# Calculo del sesgo, asumiendo p=0.5 probabilidad de que la planta esté enferma.
sesgo <- mean(df_estimacion$p_estimados) - 0.5
# Variabilidad de las estimaciones
cv <- sd(df_estimacion$p_estimados) / mean(df_estimacion$p_estimados)
cat("Simetría es ", simetria$Skewness, "\n")Simetría es 0.1519519 cat("Sesgo es", sesgo, "\n")Sesgo es -0.0044 cat("Coeficiente de variación es", cv)Coeficiente de variación es 0.2131859Segun los resultados de n=20, se infiere que las estimaciones de p^ presentan un comportamiento que tiende a la simetría, y que se distribuyen de manera Normal con una simetría muy cercana a cero y un coeficiente de variación bajo. Además, el sesgo es mínimo cuando se compara la media de las estimaciones con el pp inicial del ejercicio.
Ahora, se repite el proceso anterior con tamaños de muestra n=5, 10, 15, 20, 30, 50, 60, 100, 200, 500, para comparar la normalidad.
n=500
X <- matrix(unlist(lapply(1:n, function(i) funmuestra(enfermas, 100)$muestra)), ncol = n)
#Generación de medias
X5=X[ ,1:5] # n=5
X10=X[ ,1:10] # n=10
X20=X[ ,1:20] # n=20
X30=X[ ,1:30] # n=30
X50=X[ ,1:50] # n=50
X60=X[ ,1:60] # n=60
X100=X[ ,1:100] # n=100
X200=X[ ,1:200] # n=200
X500=X[ ,1:500] # n=500
# Obtención de las medias para diferentes tamaños de muestra.
Mx5=apply(X5,1,mean) # medias de muestras de tamaño n=5
Mx10=apply(X10,1,mean) # medias de muestras de tamaño n=10
Mx20=apply(X20,1,mean) # medias de muestras de tamaño n=20
Mx30=apply(X30,1,mean) # medias de muestras de tamaño n=30
Mx50=apply(X50,1,mean) # medias de muestras de tamaño n=50
Mx60=apply(X60,1,mean) # medias de muestras de tamaño n=60
Mx100=apply(X100,1,mean) # medias de muestras de tamaño n=100
Mx200=apply(X200,1,mean) # medias de muestras de tamaño n=200
Mx500=apply(X500,1,mean) # medias de muestras de tamaño n=500Verificación de normalidad en los datos con histogramas de cada muestra
par(mfrow=c(3,6), mar = c(4.1, 3.1, 3.1, 1.1)) # divide la gráfica en una matriz
hist(Mx5, main = "n=5", freq=FALSE, xlab=expression(hat(p)), ylab="", col="#d6a692")
qqnorm(Mx5, main="n=5", xlab="cuantiles teóricos", ylab="") ; qqline(X5, col="red")
hist(Mx10, main = "n=10",freq=FALSE, xlab=expression(hat(p)), ylab="", col="#d6a692")
qqnorm(Mx10, main ="n=10", xlab="cuantiles teóricos", ylab="") ; qqline(Mx10, col="red")
hist(Mx20, main = "n=20",freq=FALSE, xlab=expression(hat(p)), ylab="", col="#d6a692")
qqnorm(Mx20, main ="n=20", xlab="cuantiles teóricos", ylab="") ; qqline(Mx20, col="red")
hist(Mx30, main = "n=30",freq=FALSE,xlab=expression(hat(p)), ylab="", col="#d6a692")
qqnorm(Mx30, main = "n=30", xlab="cuantiles teóricos", ylab="") ; qqline(Mx30, col="red")
hist(Mx50, main = "n=50",freq=FALSE,xlab=expression(hat(p)), ylab="", col="#d6a692")
qqnorm(Mx50, main = "n=50", xlab="cuantiles teóricos", ylab="") ; qqline(Mx50, col="red")
hist(Mx60, main = "n=60", freq = FALSE,xlab=expression(hat(p)), ylab="", col="#d6a692")
qqnorm(Mx60, main ="n=60", xlab="cuantiles teóricos", ylab="") ; qqline(Mx60, col="red")
hist(Mx100, main = "n=100", freq = FALSE,xlab=expression(hat(p)), ylab="", col="#d6a692")
qqnorm(Mx100, main ="n=100",xlab="cuantiles teóricos" ,ylab="") ; qqline(Mx100, col="red")
hist(Mx200, main = "n=200",freq = FALSE,xlab=expression(hat(p)), ylab="", col="#d6a692")
qqnorm(Mx200, main ="n=200", xlab="cuantiles teóricos", ylab="") ; qqline(Mx200, col="red")
hist(Mx500, main = "n=500", freq = FALSE,xlab=expression(hat(p)), ylab="", col="#d6a692")
qqnorm(Mx500, main="n=500", xlab="cuantiles teóricos", ylab="") ; qqline(Mx500, col="red")
Cálculo de shapiro wilk a cada uno de los tamaños de muestra
sw5 <- shapiro.test(Mx5)
sw10 <- shapiro.test(Mx10)
sw20 <- shapiro.test(Mx20)
sw30 <- shapiro.test(Mx30)
sw50 <- shapiro.test(Mx50)
sw60 <- shapiro.test(Mx60)
sw100 <- shapiro.test(Mx100)
sw200 <- shapiro.test(Mx200)
sw500 <- shapiro.test(Mx500)
muestras <- c("n=5", "n=10", "n=20", "n=30", "n=50", "n=60", "n=100", "n=200", "n=500")
valores_p <- c(sw5$p.value, sw10$p.value, sw20$p.value, sw30$p.value, sw50$p.value, sw60$p.value, sw100$p.value, sw200$p.value, sw500$p.value)
tabla_p <- data.frame(muestra = muestras, "Valor p" = valores_p)
#install.packages("gt")
#library(gt)
#tabla_p %>%
# gt()Observando los resultados anteriores, se evidencia que a medida que el tamaño de muestra aumenta, los histogramas tienden a ser más simétricos y que los puntos qqnorm siguen la línea de base. Además el test de SW a partir del tamaño n>50 comienza a no rechazarse H0 al tener un valor p > 5%.
Ahora se repite toda la simulación con lotes de 10% de plantas enfermas y lotes con un 90% de plantas enfermas.
n_lote <- 1000
p_90 <- 0.9
enfermas <- rbinom(n_lote, size = 1, prob = p_90)
data_lote <- data.frame(id=1:n_lote, tipo = enfermas) # 1. planta enferma, 0: planta sana
#---------------------------------------------------------------
#Paso 2
funcion_muestra <- function(data, n_muestra) {
muestra <- sample(data_lote$tipo, n_muestra, replace = FALSE)
p_estimado <- mean(muestra)
return(list(muestra = muestra, p_estimado = p_estimado))
}
## Aplicamos un ejemplo para ver salida de la función:
ej <- funcion_muestra(data_lote, 20)
cat("La muestra extraída ", ej$muestra, " tiene un p estimado =", ej$p_estimado)La muestra extraída 0 1 1 1 1 0 1 1 1 1 1 1 1 1 1 0 1 1 1 1 tiene un p estimado = 0.85n_muestra <- 20
n_repeticiones <- 500
df_estimacion <- data.frame(iteracion = numeric(n_repeticiones), p_estimados = numeric(n_repeticiones))
for (i in 1:n_repeticiones) {
df_estimacion$iteracion[i] <- i
muestra_resultado <- funcion_muestra(data_lote, n_muestra)
df_estimacion$p_estimados[i] <- muestra_resultado$p_estimado
}
#---------------------------------------------------------------
#Paso 4
# Se realiza gráfico para mirar simetría de la distribución
library(ggplot2)
library(tidyverse)
library(hrbrthemes)
ggplot(df_estimacion, aes(x=p_estimados)) +
geom_histogram(binwidth=0.05, fill="#5191c1", color="#e9ecef", alpha=0.9) +
ggtitle(bquote("Distribución de " ~ hat(p) ~ "")) +
labs(x=(bquote("" ~ hat(p) ~ "")), y="densidad") +
theme_ipsum() +
theme(
plot.title = element_text(size=15)
)Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font database
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font database
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font database
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font database
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font database
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font database
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font database
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font database
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font database
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font database
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font database
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font database
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font database
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font database
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
family not found in Windows font databaseWarning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Ahora se calcula para plantas enfermas 90%
library(parameters)
simetria <- skewness(df_estimacion$p_estimados)
# Se calcula el sesgo, asumiendo que p=0.9 que es la probabilidad de que la planta esté enferma.
sesgo <- mean(df_estimacion$p_estimados) - p_90
# Variabilidad de las estimaciones
cv <- sd(df_estimacion$p_estimados) / mean(df_estimacion$p_estimados)
cat("La simetría es igual a ", simetria$Skewness, "\n")La simetría es igual a -0.6225179 cat("El sesgo es igual a ", sesgo, "\n")El sesgo es igual a -0.0039 cat("El coeficiente de variación es igual a ", cv)El coeficiente de variación es igual a 0.07567887Ahora con n=500
n=500
X <- matrix(unlist(lapply(1:n, function(i) funcion_muestra(enfermas, 100)$muestra)), ncol = n)
#Generación de medias
X5=X[ ,1:5] # n=5
X10=X[ ,1:10] # n=10
X20=X[ ,1:20] # n=20
X30=X[ ,1:30] # n=30
X50=X[ ,1:50] # n=50
X60=X[ ,1:60] # n=60
X100=X[ ,1:100] # n=100
X200=X[ ,1:200] # n=200
X500=X[ ,1:500] # n=500
# Obtención de las medias para diferentes tamaños de muestra.
Mx5=apply(X5,1,mean) # medias de muestras de tamaño n=5
Mx10=apply(X10,1,mean) # medias de muestras de tamaño n=10
Mx20=apply(X20,1,mean) # medias de muestras de tamaño n=20
Mx30=apply(X30,1,mean) # medias de muestras de tamaño n=30
Mx50=apply(X50,1,mean) # medias de muestras de tamaño n=50
Mx60=apply(X60,1,mean) # medias de muestras de tamaño n=60
Mx100=apply(X100,1,mean) # medias de muestras de tamaño n=100
Mx200=apply(X200,1,mean) # medias de muestras de tamaño n=200
Mx500=apply(X500,1,mean) # medias de muestras de tamaño n=500
# Gráficos de histogramas y qqnorm de cada uno de los tamaños de muestra
par(mfrow=c(3,6), mar = c(4.1, 3.1, 3.1, 1.1)) # divide la gráfica en una matriz
hist(Mx5, main = "n=5", freq=FALSE, xlab=expression(hat(p)), ylab="",
col="#5191c1")
qqnorm(Mx5, main="n=5", xlab="cuantiles teóricos", ylab="") ; qqline(X5, col="red")
hist(Mx10, main = "n=10",freq=FALSE, xlab=expression(hat(p)), ylab="",
col="#5191c1")
qqnorm(Mx10, main ="n=10", xlab="cuantiles teóricos", ylab="") ; qqline(Mx10, col="red")
hist(Mx20, main = "n=20",freq=FALSE, xlab=expression(hat(p)), ylab="",
col="#5191c1")
qqnorm(Mx20, main ="n=20", xlab="cuantiles teóricos", ylab="") ; qqline(Mx20, col="red")
hist(Mx30, main = "n=30",freq=FALSE,xlab=expression(hat(p)), ylab="",
col="#5191c1")
qqnorm(Mx30, main = "n=30", xlab="cuantiles teóricos", ylab="") ; qqline(Mx30, col="red")
hist(Mx50, main = "n=50",freq=FALSE,xlab=expression(hat(p)), ylab="",
col="#5191c1")
qqnorm(Mx50, main = "n=50", xlab="cuantiles teóricos", ylab="") ; qqline(Mx50, col="red")
hist(Mx60, main = "n=60", freq = FALSE,xlab=expression(hat(p)), ylab="",
col="#5191c1")
qqnorm(Mx60, main ="n=60", xlab="cuantiles teóricos", ylab="") ; qqline(Mx60, col="red")
hist(Mx100, main = "n=100", freq = FALSE,xlab=expression(hat(p)), ylab="",
col="#5191c1")
qqnorm(Mx100, main ="n=100",xlab="cuantiles teóricos" ,ylab="") ; qqline(Mx100, col="red")
hist(Mx200, main = "n=200",freq = FALSE,xlab=expression(hat(p)), ylab="",
col="#5191c1")
qqnorm(Mx200, main ="n=200", xlab="cuantiles teóricos", ylab="") ; qqline(Mx200, col="red")
hist(Mx500, main = "n=500", freq = FALSE,xlab=expression(hat(p)), ylab="",
col="#5191c1")
qqnorm(Mx500, main="n=500", xlab="cuantiles teóricos", ylab="") ; qqline(Mx500, col="red")
sw5 <- shapiro.test(Mx5)
sw10 <- shapiro.test(Mx10)
sw20 <- shapiro.test(Mx20)
sw30 <- shapiro.test(Mx30)
sw50 <- shapiro.test(Mx50)
sw60 <- shapiro.test(Mx60)
sw100 <- shapiro.test(Mx100)
sw200 <- shapiro.test(Mx200)
sw500 <- shapiro.test(Mx500)
muestras <- c("n=5", "n=10", "n=20", "n=30", "n=50", "n=60", "n=100", "n=200", "n=500")
valores_p <- c(sw5$p.value, sw10$p.value, sw20$p.value, sw30$p.value, sw50$p.value, sw60$p.value, sw100$p.value, sw200$p.value, sw500$p.value)
tabla_p <- data.frame(muestra = muestras, "Valor p" = valores_p)
#library(gt)
#tabla_p %>%
# gt()El primer histograma evidencia que existe una asimetría negativa al estar las estimaciones concentradas a la derecha, la cual puede observarse en el valor de simetría calculada con valor negativo. A pesar de haber asimetría en la distribución, el sesgo expone que el valor de p^ es muy cercano a p con un coeficiente de variación muy cercano a cero. Con los diferentes tamaños de muestra, los histogramas y qqnorms demuestran que se alcanza un comportamiento de normalidad a partir de n>200 y donde el test de SW lo confirma cuando se deja de rechazar H0 a partir de este tamaño de muestra.
Conclusiones
A partir de la información analizada se observa que la tendencia a la normalidad y cumplimiento del teorema del límite central se logra a medida que aumenta el tamaño de la muestra, pues el valor de P estimado, se acerca más al valor de P, disminuye la varianza, y converge hacia una distribución normal. El tener un p=0.5 inicial ayuda a lograr una simetría desde el inicio y con una muestra mayor a 30, p^ se acercaba mucho más al valor real pp. El estimador tiende a acercarse a la proporción real (0.5, 0.10, 0.90) a medida que el tamaño de la muestra aumenta. Respecto al sesgo, en general es bajo para todos los tamaños de muestra, lo que indica que el estimador es no sesgado o de bajo sesgo. Asi mismo, la varianza disminuye a medida que el tamaño de la muestra aumenta. En cuanto al test de Shapiro, para tamaños de muestra pequeños se rechaza la hipótesis nula de normalidad debido a los bajos valores p. Sin embargo, a medida que el tamaño de la muestra crece, se tiende a no rechazar dicha hipótesis, lo que sugiere que las muestras tienden a ser normales como se indico al inicio. Finalmente, para un analisis más certero del comportamiento de la población enferma se requiere de tamaños de muestra más grandes para que la distribución muestral se acerque a la normalidad.
Ejercicio 4
Estimación bootstrap
Con el proposito de saber si la muestra es representativa de la población, se utiliza el método de estimación bootstrap.
x=c(7.69, 4.97, 4.56, 6.49, 4.34, 6.24, 4.45) # datos muestra
boot=sample(x,7000,replace=TRUE) # se extraen n x m muestras
b=matrix(boot,nrow=1000,ncol=7) # se construye matriz de n x m
mx=apply(b,1,mean)
ic1=quantile(mx, probs=c(0.025, 0.975)) # se calcula IC método 1
ic1 2.5% 97.5%
4.709893 6.434357 ic2=c(2*mean(mx)-ic1[2], 2*mean(mx)-ic1[1]) # se calcula IC método 2
ic2 97.5% 2.5%
4.577691 6.302156 hist(mx, las=1, main=" ", ylab = " ", xlab = " ", col="#034A94")
abline(v=ic1, col="#FF7F00",lwd=2)
abline(v=ic2, col="#0EB0C6",lwd=2)
Conclusión
Como se pude observar a partir de una población de 7 muestras, con una replica de 1000 datos para cada muestra, se obtiene unos intervalos de la curva de densidad con valores minimos y maximos para cada estimación empleando los intervalos basados en percentiles para obtener el rango de valores en el que se puede encontrar el parámetro mediante dos metodos. Para un intervalo de confianza del 95% se encuentra que el primer método se basa exclusivamente en los percentiles para construir el intervalo, mientras que el segundo método también incluye la media de todas las observaciones y aporta a la simetría. Si bien, la diferencia entre los intervalos de confianza es muy pequeña en ambos casos, se observa que el metodo 2 se aproxima más al valor 6.5 de eficiencia del combustible. Sin embargo, la muestra es muy poco representativa de la población dado su tamaño, y por tanto las estimaciones efectuadas se deben analizar con cautela, pues el sesgo entre la media de la muestra original y la media de todas las muestras bootstrap es pequeño (0.0052); la desviacion 1.28; varianza 1.65 y coeficientete de variacion 0.23.
Ejercicio 5
Relaciones entre la potencia, el tamaño de los efectos y tamaño de la muestra
Para cada tamaño fijo de los efectos d, se modela la relación entre el tamaño muestral y la potencia (manteniendo constante el nivel de significancia α=0.05). En las siguientes figuras se visualizan los resultados para tamaño de efecto muy pequeño (d=0.1), pequeño (d=0.2), mediano (d=0.5) y grande (d=0.8). Repite el análisis usando 5 valores distintos del nivel de significancia. ¿Cambian los resultados? ¿Qué ocurre cuando el tamaño de muestra de los grupos que se comparan es de 20, 60, 100 y 140? Analiza y compara los resultados.
if(!require(pwr)){install.packages("pwr");library("pwr")}Loading required package: pwr# t-TEST
# Se aplicará power.t.test del paquete stats (ya en R). Calcula la potencia de la prueba t de una o dos muestras, o determina los parámetros para obtener un valor particular de la potencia.
d<-seq(.1,2,by=.1) # 20 tamaños de los efectos
n<-1:150 # Tamaños muestrales
t.test.power.effect <-as.data.frame(do.call("cbind",lapply(1:length(d),function(i)
{
sapply(1:length(n),function(j)
{
power.t.test(n=n[j],d=d[i],sig.level=0.05,power=NULL,type= "two.sample")$power
})
})))
# Si algunas potencias no se pueden calcular, se ajustan a cero:
t.test.power.effect[is.na(t.test.power.effect)] <- 0
colnames(t.test.power.effect)<-paste (d,"effect size")
#Graficando los resultados
prueba <-t.test.power.effect #data frame de 150 X 20 (para graficar)
cuts_num<-c(2,5,8) # cortes
#Cortes basados en: Cohen, J. (1988). Statistical Power Analysis for the Behavioral Sciences (2nd ed.). Hillsdale, NJ: Lawrence Erlbaum Associates, Publishers.
cuts_cat<-c("pequeño","medio","grande")
columnas <- 1:ncol(prueba) #Lista de los valores 1:20
color_linea<-rainbow(length(columnas), alpha=.5) # Lista de 20 colores
grosor_linea=3 # Grosor de la línea
#Para el tipo de línea: (“blank”, “solid”, “dashed”, “dotted”, “dotdash”, “longdash”, “twodash”) ó (0, 1, 2, 3, 4, 5, 6).
#Note que lty = “solid” is idéntica a lty=1.
tipo_linea <- rep(1,length(color_linea)) #Repetir length(color)=20 veces el 1
tipo_linea[cuts_num]<-c(2:(length(cuts_num)+1)) #Asignar 2, 3, 4 en las posiciones 2, 5, 8 de tipo_linea
#Resaltar posiciones importantes
cuts_num<-c(2,5,8) # Cortes
#Cortes basados en: Cohen, J. (1988). Statistical Power Analysis for the Behavioral Sciences (2nd ed.). Hillsdale, NJ: Lawrence Erlbaum Associates, Publishers.
cuts_cat<-c("pequeño","medio","grande")
color_linea[cuts_num]<-c("black")
efecto <- d # Listado de los 20 valores de 20
efecto[cuts_num] <- cuts_cat #Reemplazar en "efecto" las posiciones cuts_num (2, 5, 8) por las categorías de cuts_cat
par(fig=c(0,.8,0,1),new=TRUE)Warning in par(fig = c(0, 0.8, 0, 1), new = TRUE): llamada par(new=TRUE) sin
gráficoplot(1, type="n", #no produce puntos ni líneas
frame.plot=FALSE,
xlab="Tamaño muestral", ylab="Potencia",
xlim=c(1,150), ylim=c(0,1),
main="t-Test", axes = FALSE)
#Editando los ejes, grid, etc.
abline(v=seq(0,150,by=10), col = "lightgray", lty = "dotted") # Grid vertical
abline(h=seq(0,1,by=.05), col = "lightgray", lty = "dotted") # Grid horizontal
axis(1,seq(0,150,by=10)) # Números en eje X
axis(2,seq(0,1,by=.05)) # Números en eje Y
#Plot de las lineas
#columnas <- 1:ncol(prueba) # lista de los valores 1:20
for(i in 1:length(columnas)) #length(columnas)=20
{
lines(1:150,
#prueba (data frame de 150 X 20, para graficar)
#columna <- 1:ncol(prueba) listado de valores 1:20
prueba[,columnas[i]], #filtrar "prueba" para valor de columna
col=color_linea[i], #color_linea[cuts_num]<-c("black")
lwd=grosor_linea, #grosor de cada linea
lty=tipo_linea[i] #tipo_linea[cuts_num]<-c(2:(length(cuts_num)+1))
)
}
#Leyendas
par(fig=c(.65,1,0,1),new=TRUE)
plot.new()
legend("top",legend=efecto, col=color_linea, lwd=3, lty=tipo_linea, title="Tamaño efecto",
bty="n" #Opciones: o (complete box), n (no box), 7, L, C, U
)