Actividad 2: Probabilidad e Inferencia Estadística
Introducción
El proceso de simulación constituye una herramienta poderosa para la estadística que se puede emplear para entender relaciones complejas y estimar valores difíciles de calcular directamente. Para entenderlo se resolveran 4 problemas.
Cargar Librerías
En primer lugar, se cargan las librerías que utilizaremos durante toda la actividad para solucionar los 4 problemas. Adicionalmente establecemos desde el principio la configuración de las opciones de formato número. En este caso que evite la notación científica y que muestre números en notación decimal estándar.
library(dplyr)
library(tidyverse)
library(summarytools)
library(moments)
library(scales)
options(scipen = 100, digits = 2)Problema 1: Estimación del valor π
La siguiente figura sugiere como estimar el valor de π con una simulación. En la figura, un circuito con un área igual a π/4, está inscrito en un cuadrado cuya área es igual a 1. Se elige de forma aleatoria n puntos dentro del cuadrado. La probabilidad de que un punto esté dentro del círculo es igual a la fracción del área del cuadrado que abarca a este, la cual es π/4. Por tanto, se puede estimar el valor de π/4 al contar el número de puntos dentro del círculo, para obtener la estimación de π/4. De este último resultado encontrar una aproximación para el valor de π.
##
Preparación:
Cada uno de esos puntos se ubica en el cuadrado con coordenadas (x,y). De esos puntos, existirán unos que se encuentren fuera del círculo y para determinar cuáles son, debemos encontrar la distancia que los separa del centro del círculo (0.5,0,5). Por lo tanto, estableceremos el número de puntos a analizar:
nPuntos <- 1000Paso 1: Coordenadas x.
Genere n coordenadas x:X1,.., Xn. Utilice la distribución uniforme con valor mínimo de 0 y valor máximo de 1. La distribución uniforme genera variables aleatorias que tienen la misma probabilidad de venir de cualquier parte del intervalo (0,1).
x <- runif(n= nPuntos, min = 0, max = 1)Paso 2: Coordenadas y.
Genere n coordenadas y: Y1,…,Yn , utilizando nuevamente la distribución uniforme con valor mínimo de 0 y valor máximo de 1.
y <- runif(n= nPuntos, min = 0, max = 1)Paso 3: Distancia.
Cada punto (Xi,Yi) se encuentra dentro del círculo si su distancia desde el centro (0.5,0.5) es menor a 0.5. Para cada par (Xi,Yi). Determine si la distancia desde el centro es menor a 0.5. Esto último se puede realizar al calcular el valor (Xi−0.5)2+(Yi−0.5)2,que es el cuadrado de la distancia, y al determinar si es menor que 0.25.
Se crea el dataframe con el par de coordenadas.
coordenadas <- data.frame(x,y)Se crea una función para calcular la distancia del par de coordenadas al centro del círculo
distancia <- function(xi, yi){
d = sqrt((xi-0.5)^2 + (yi-0.5)^2)
return(d)
}Se aplica la función al dataframe que contiene los puntos con sus respectivas coordenadas
dist_puntos <- apply(coordenadas, 1, function(row){
distancia(row["x"], row["y"])
})Número de puntos dentro del círculo y estimación de π
Se crea un campo numérico que evalúe si las distancias son mayores a 0.5. Si la distancia es mayor a 0.5 evaluar=0 de lo contrario evaluar =1
dist_puntos <-data.frame(dist_puntos) %>%
mutate(evaluar = ifelse(dist_puntos>0.5, 0,1))Paso 4: Respuestas.
¿Cuántos de los puntos están dentro del círculo?
El número de puntos que están dentro del círculo serán los que cumplen la condición de que evaluar=1.
puntosCirculo <- sum(dist_puntos$evaluar)
puntosCirculo## [1] 771¿Cuál es su estimación de π?
La relación entre el área del círculo y el cuadrado es de π/4. El área del cículo vendrían hacer todos los puntos que se encuentran dentro del mismo;mientras que el área del cuadrado serían todos los n puntos.Por lo tanto π es la suma de los puntos que están dentro del círculo por 4 sobre n.
π <- puntosCirculo*4/nPuntos
π## [1] 3.1Aplicar a diferente cantidad de puntos
Teniendo en cuenta todo lo desarrollado anteriormente, se crea una función que contenga los pasos y que nos arroje los dos resultados que necesitamos: El número de puntos que están dentro del cículo y la estimación de π, adicionalmente para hacer una diferenciación agregamos también el número de puntos que estamos analizando.
estimacion_ValorPi <- function(puntos){
x <- runif(n= puntos, min = 0, max = 1)
y <- runif(n= puntos, min = 0, max = 1)
coordenadas <- data.frame(x,y)
dist_puntos <- apply(coordenadas, 1, function(row){
distancia(row["x"], row["y"])
})
dist_puntos <-data.frame(dist_puntos) %>%
mutate(evaluar = ifelse(dist_puntos>0.5, 0,1))
puntosCirculo <- sum(dist_puntos$evaluar)
π <- comma(puntosCirculo*4/puntos, accuracy = 0.001)
resultados <- list(TotalPuntos = puntos, PuntosCirculo = puntosCirculo, Estimacion=π)
return(resultados)
}Los puntos a analizar son 1000, 10000 y 100000. Estos se agregaran en un vector y posteriormente utilizaremos la función sapply para aplicar la función a cada uno de estos puntos de tal forma que nos arroje los resultados.
puntos <- c(1000, 10000, 100000)
resultado <- sapply(puntos, estimacion_ValorPi)
resultado## [,1] [,2] [,3]
## TotalPuntos 1000 10000 100000
## PuntosCirculo 775 7870 78740
## Estimacion "3.100" "3.148" "3.150"Conclusiones.
Como podemos ver, a medida que incrementamos el número de puntos, el valor estimado de π es más cercano a su valor verdadero.
Problema 2: Propiedades de los estimadores
Sean X1, X2, X3 y X4 , una muestra aleatoria de tamaño n=4 cuya población la conforma una distribución exponencial con parámetro θ desconocido. Determine las características de cada uno de los siguientes estimadores propuestos:
Paso 1: Generación de muestras.
Se crea un objeto que tenga el número de valores de la muestra
muestra<- 20Se crean las muestras con distribución exponencial.
lambda <- 2
x1 <- rexp(muestra, lambda)
x2 <- rexp(muestra, lambda)
x3 <- rexp(muestra, lambda)
x4 <- rexp(muestra, lambda)En un solo dataframe dejamos las muestras.
base <- data.frame(x1,x2,x3,x4)Paso 2: Creación de funciones asociadas a cada estimador.
Estimador 1
θ1^=X1+X2/6+X3+X4/3
estimador1 <- function(x1,x2,x3,x4){
est = (x1+x2)/6 + (x3+x4)/3
return(est)
}Estimador 2
θ2^=(X1+2X2+3X3+4X4)/5
estimador2 <- function(x1,x2,x3,x4){
est = (x1+2*x2+3*x3+4*x4)/5
return(est)
}Estimador 3
θ3^=X1+X2+X3+X4/4
estimador3 <- function(x1,x2,x3,x4){
est = (x1+x2+x3+x4)/4
return(est)
}Estimador 4
θ4^=min{X1,X2,X3,X4}+max{X1,X2,X3,X4}/2
estimador4 <- function(x1,x2,x3,x4){
est = (min(x1,x2,x3,x4)+ max(x1,x2,x3,x4))/2
return(est)
}Paso 3: Evaluación cada estimador con la base creada
E1 <- apply(base, 1, function(row){
estimador1(row["x1"], row["x2"], row["x3"], row["x4"])
})
E2 <- apply(base, 1, function(row){
estimador2(row["x1"], row["x2"], row["x3"], row["x4"])
})
E3 <- apply(base, 1, function(row){
estimador3(row["x1"], row["x2"], row["x3"], row["x4"])
})
E4<- apply(base, 1, function(row){
estimador4(row["x1"], row["x2"], row["x3"], row["x4"])
}) Se crea un dataframe que contenga toda la data de los estimadores
estimadores <- data.frame(E1,E2,E3,E4) Paso 4: Crear una función para el problema.
De acuerdo con todo lo desarrollado anteriormente y teniendo en cuenta que vamos analizar diferentes muestras, creamos una sola función que tenga como parámetros la muestra y el valor de lambda.
estimacion_Propiedades <- function(muestras, lambda){
#Creación de Muestras
x1 <- rexp(muestras, lambda)
x2 <- rexp(muestras, lambda)
x3 <- rexp(muestras, lambda)
x4 <- rexp(muestras, lambda)
base <- data.frame(x1,x2,x3,x4)
# Aplicar a la base cada una de las funciones creadas para los estimadores.
E1 <- apply(base, 1, function(row){
estimador1(row["x1"], row["x2"], row["x3"], row["x4"])
})
E2 <- apply(base, 1, function(row){
estimador2(row["x1"], row["x2"], row["x3"], row["x4"])
})
E3 <- apply(base, 1, function(row){
estimador3(row["x1"], row["x2"], row["x3"], row["x4"])
})
E4<- apply(base, 1, function(row){
estimador4(row["x1"], row["x2"], row["x3"], row["x4"])
})
estimadores <- data.frame(E1,E2,E3,E4)
return(estimadores)
}Análisis de resultados.
Se analizan los estimadores a partir de la visualización del gráfico de bigotes. Con la función descr se analiza la media y la varianza. En este caso E[X] = 1/λ Como λ = 2, entonces E[X]= 1/2 , E[X] = 0.5
Muestra: 20
M1 <-estimacion_Propiedades(20,2)
boxplot(M1,col = rainbow(ncol(M1)))
abline(h=0.5, col="red")
descr(M1)## Descriptive Statistics
## M1
## N: 20
##
## E1 E2 E3 E4
## ----------------- -------- -------- -------- --------
## Mean 0.57 1.15 0.56 0.61
## Std.Dev 0.23 0.49 0.20 0.22
## Min 0.18 0.40 0.17 0.20
## Q1 0.42 0.73 0.44 0.50
## Median 0.55 1.15 0.56 0.65
## Q3 0.75 1.56 0.70 0.73
## Max 1.01 1.85 1.01 1.01
## MAD 0.26 0.61 0.18 0.18
## IQR 0.32 0.81 0.25 0.21
## CV 0.40 0.42 0.36 0.36
## Skewness 0.07 -0.03 0.15 -0.07
## SE.Skewness 0.51 0.51 0.51 0.51
## Kurtosis -1.15 -1.55 -0.42 -0.80
## N.Valid 20.00 20.00 20.00 20.00
## Pct.Valid 100.00 100.00 100.00 100.00Muestra: 50
M2 <-estimacion_Propiedades(50,2)
boxplot(M2,col = rainbow(ncol(M2)))
abline(h=0.5, col="red")
descr(M2)## Descriptive Statistics
## M2
## N: 50
##
## E1 E2 E3 E4
## ----------------- -------- -------- -------- --------
## Mean 0.52 1.05 0.52 0.64
## Std.Dev 0.27 0.61 0.27 0.36
## Min 0.09 0.21 0.10 0.13
## Q1 0.27 0.59 0.26 0.32
## Median 0.54 1.04 0.54 0.63
## Q3 0.65 1.30 0.66 0.86
## Max 1.34 3.11 1.46 1.80
## MAD 0.29 0.58 0.25 0.42
## IQR 0.36 0.71 0.39 0.52
## CV 0.52 0.58 0.51 0.56
## Skewness 0.88 1.28 0.82 0.91
## SE.Skewness 0.34 0.34 0.34 0.34
## Kurtosis 0.77 1.96 1.22 0.78
## N.Valid 50.00 50.00 50.00 50.00
## Pct.Valid 100.00 100.00 100.00 100.00Muestra: 100
M3 <-estimacion_Propiedades(100,2)
boxplot(M3,col = rainbow(ncol(M3)))
abline(h=0.5, col="red")
descr(M3)## Descriptive Statistics
## M3
## N: 100
##
## E1 E2 E3 E4
## ----------------- -------- -------- -------- --------
## Mean 0.49 0.97 0.50 0.58
## Std.Dev 0.27 0.53 0.27 0.32
## Min 0.10 0.20 0.09 0.10
## Q1 0.30 0.55 0.31 0.35
## Median 0.43 0.86 0.43 0.49
## Q3 0.62 1.23 0.61 0.76
## Max 1.43 2.73 1.39 1.62
## MAD 0.22 0.49 0.22 0.27
## IQR 0.32 0.67 0.30 0.40
## CV 0.54 0.55 0.54 0.55
## Skewness 1.26 1.09 1.24 1.04
## SE.Skewness 0.24 0.24 0.24 0.24
## Kurtosis 1.65 0.86 1.50 0.53
## N.Valid 100.00 100.00 100.00 100.00
## Pct.Valid 100.00 100.00 100.00 100.00Muestra: 1000
M4 <-estimacion_Propiedades(1000,2)
boxplot(M4,col = rainbow(ncol(M4)))
abline(h=0.5, col="red")
descr(M4)## Descriptive Statistics
## M4
## N: 1000
##
## E1 E2 E3 E4
## ----------------- --------- --------- --------- ---------
## Mean 0.49 0.97 0.49 0.57
## Std.Dev 0.26 0.53 0.24 0.32
## Min 0.02 0.04 0.02 0.02
## Q1 0.30 0.59 0.31 0.34
## Median 0.44 0.88 0.44 0.50
## Q3 0.63 1.25 0.62 0.74
## Max 1.72 3.35 1.51 2.07
## MAD 0.24 0.48 0.23 0.28
## IQR 0.33 0.65 0.31 0.40
## CV 0.53 0.54 0.49 0.55
## Skewness 1.09 1.13 0.87 1.24
## SE.Skewness 0.08 0.08 0.08 0.08
## Kurtosis 1.60 1.62 0.84 2.13
## N.Valid 1000.00 1000.00 1000.00 1000.00
## Pct.Valid 100.00 100.00 100.00 100.00Conclusiones.
Analizando cada uno de los casos encontramos que los mejores estimadores eran el 1 y el 3.
El estimador 1 era el que más se acercaba al valor esperado (Insesgado).
El estimador 3 era el que presentaba menor varianza (Eficiente).
El estimador 4 estaba inicialmente lejos del valor esperado pero a medida que incrementó la muestra, empezó a acercarse al valor esperado (Consistencia).
El estimador 2 es el que presentaba sesgadez y mayor varianza, a medida que incrementó la muestra, empezó acercarse al valor esperado (Consistencia) pero no a la misma velocidad que el estimador 4.
Problema 3: Teorema del límite central
El Teorema del Límite Central es uno de los más importantes en la inferencia estadística y habla sobre la convergencia de los estimadores como la proporción muestral a la distribución normal. Algunos autores afirman que esta aproximación es bastante buena a partir del umbral n>30.
Paso 1: Simulación.
Realice una simulación en la cual genere una población de n=1000(Lote), donde el porcentaje de individuos (supongamos plantas) enfermas sea del 50%.
Se establece N (Población)
N <- 1000Se crea una función que genere población con las características dadas. Los parámetros son: la población que corresponde al número de plantas, la variable 1 y variable 2 que corresponden a la proporción que existe entre plantas sanas y enfermas, en este caso, 50% sanas y el otro 50% enfermas
simulacion <- function(poblacion, var1, var2) {
E = round(poblacion * var1, 0)
S = round(poblacion * var2, 0)
Poblacion = c(rep("Enferma", E), rep("Sana",S))
print(paste("Lote Generado de tamaño: ", poblacion, " Enfermos: ",E,"Sanos: ",S))
return(Poblacion)
}
PoblacionGenerada <- simulacion(N, 0.5,0.5)## [1] "Lote Generado de tamaño: 1000 Enfermos: 500 Sanos: 500"Paso 2: Muestra Aleatoria y Estimador.
Genere una función que permita: Obtener una muestra aleatoria de la población y calcule el estimador de la proporción muestral pˆ para un tamaño de muestra dado n. Se crea una función para obtener una muestra aleatoria de la población, los parámetros son la población que se está analizando y el tamaño que queremos para la muestra.
muestraAleatoria<- function(poblacion, tamaño){
muestra = sample(x = poblacion, size = tamaño)
return(muestra) }Con ayuda de la función obtenemos una muestra aleatoria de la población.
Tamaño de la muestra
n <-100Información de la muestra
data_muestra <- muestraAleatoria(PoblacionGenerada, n)Creamos una función para calcular el estimador de la proporción muestral pˆ. Los parámetros son: la muestra que está asociada a todos los datos de la muestra, el tamaño que corresponde al tamaño de la muestra y la variable a la cual queremos calcular la proporción.
proporcion_muestra <- function(muestra,tamaño,variable){
proporcion = sum(muestra == variable)/tamaño
return(proporcion)
}Como más adelante necesitaremos analizar el estimador, combinaremos las funciones.
CalcularEst <-function(poblacion, tamaño, variable){
muestra = sample(x = poblacion, size = tamaño)
proporcion = sum(muestra == variable)/tamaño
return(proporcion)
}Calculamos el estimador de la proporción.
pestimador <- CalcularEst(PoblacionGenerada,n,"Enferma")Paso 3: Repeticiones.
Repita el escenario anterior (b) n=500 veces y analice los resultados en cuanto al comportamiento de los 500 resultados del estimador pˆ. ¿Qué tan simétricos o sesgados son los resultados obtenidos? y ¿qué se puede observar en cuanto a la variabilidad?.Realice en su informe un comentario sobre los resultados obtenidos.
Se crea una función que permita iteraciones. Los parámetros son la poblaciónq que se está analizando,el tamño de la muestra, la variable que estamos analizando y el número de iteraciones. Creamos un array que tenga el tamaño de las iteraciones y que inicialmente va a tener como valores NA. Cada iteración aplicará la función creada anteriormente que calcula el estimador de la proporción muestral pˆ para un tamaño de muestra dado n.
Calcular_Iter = function(poblacion, tamaño, variable, iteraciones) {
if (tamaño<=length(poblacion)) {
estimadores = array(NA,iteraciones)
for (i in 1:iteraciones) {
estimadores[i] = CalcularEst (poblacion, tamaño, variable)
}
return(estimadores)
}
else {
print("El tamaño de la muestra supera el tamaño de la población.")
}
}Número de iteraciones
iter <- 500
dataEstimadores <-Calcular_Iter(PoblacionGenerada, n, "Enferma", iter)
dataEstimadoresDF <- data.frame(EstimadorP=dataEstimadores)Histograma y análisis de medidas
hist(dataEstimadores, prob = TRUE, main = "Histograma con curva de densidad")
mediaEst <- mean(dataEstimadores)
abline(v=mediaEst, col="#2E86C1", lwd=3)
lines(density(dataEstimadores), col = 2, lwd = 2)
Resumen<-data.frame("Muestra"=0,
"TamañoMuestra"=n,
"Media"=mean(dataEstimadores),
"Mediana"=median(dataEstimadores),
"Desvest"=sd(dataEstimadores),
"Varianza"=var(dataEstimadores),
"Mín."=min(dataEstimadores),
"Máx"=max(dataEstimadores),
"Asimetría"=skewness(dataEstimadores),
"Curtosis"= kurtosis(dataEstimadores))
descr(dataEstimadoresDF)## Descriptive Statistics
## dataEstimadoresDF$EstimadorP
## N: 500
##
## EstimadorP
## ----------------- ------------
## Mean 0.50
## Std.Dev 0.05
## Min 0.35
## Q1 0.47
## Median 0.50
## Q3 0.53
## Max 0.62
## MAD 0.04
## IQR 0.06
## CV 0.10
## Skewness -0.13
## SE.Skewness 0.11
## Kurtosis -0.10
## N.Valid 500.00
## Pct.Valid 100.00¿Qué tan simétricos o sesgados son los resultados obtenidos?
Analizando el histograma y el coeficiente de simetría podemos determinar si los resultados son simétricos o sesgados. En este caso el coeficiente de asimetría es positivo, lo que sugiere que la distribución de los datos es ligeramente sesgada hacia la derecha, es decir, tiene una cola más larga en el lado derecho. Sin embargo, el sesgo es bastante pequeño, por lo que la distribución podría considerarse aproximadamente simétrica. El hecho de que el coeficiente de asimetría sea cercano a cero sugiere que la distribución no está fuertemente sesgada en ninguna dirección y se asemeja a una distribución normal.
¿Qué se puede observar en cuanto a la variabilidad?
En base a las estadísticas descriptivas proporcionadas, los datos tienen una desviación estándar relativamente baja, un rango y un IQR moderados, y un coeficiente de variación del 10%. Esto sugiere que la variabilidad en los datos no es muy alta y que los valores tienden a estar relativamente cerca de la media.
Paso 4: Análsis con diferentes muestras.
Repita los puntos 2 y 3 para tamaños de muestra n=5, 10, 15, 20, 30, 50, 60, 100, 200, 500. Compare los resultados obtenidos para los diferentes tamaños de muestra en cuanto a la normalidad. Utilice pruebas de bondad y ajuste (shapiro wilks :shspiro.test()) y métodos gráficos (gráfico de normalidad: qqnorm()). Comente en su informe los resultados obtenidos.
Establecemos una lista con las muestras indicadas. Utilizamos la función par() para establecer los parámetros en los que presentaremos la información gráficamente. El parámetro será mfrow que controla la disposición de los paneles en una ventana. En este caso, presentaremos en una fila 2 gráficas. La primera gráfica será un histograma y la segunda gráfica será un gráfico Q-Q para evaluar si nuestros datos se ajustan a una distribución normal.
muestras <- c(5,10,15,20,30,50,60,100,200,500)
for(i in 1:length(muestras)){
dataEstimadores <-Calcular_Iter(PoblacionGenerada, muestras[i], "Enferma", iter)
par(mfrow=c(1,2))
hist(dataEstimadores, las=1, ylab = "Frecuencia",
main = paste("Muestra de tamaño:", muestras[i]),col = "#2E4053",prob = TRUE)
abline(v=mean(dataEstimadores), col="#2E86C1", lwd=3)
lines(density(dataEstimadores), col = 2, lwd = 2)
qqnorm(dataEstimadores, xlab="Cuantiles teóricos",
ylab="Cuantiles muestrales",main="qq-normalidad")
qqline(dataEstimadores,col = '#D81B60', lwd = 2, lty = 2)
print(paste('Para una muestra de tamaño: ', muestras[i]))
print(shapiro.test(dataEstimadores))
x=data.frame("Muestra"=i,
"TamañoMuestra"=muestras[i],
"Media"=mean(dataEstimadores),
"Mediana"=median(dataEstimadores),
"Desvest"=sd(dataEstimadores),
"Varianza"=var(dataEstimadores),
"Mín."=min(dataEstimadores),
"Máx"=max(dataEstimadores),
"Asimetría"=skewness(dataEstimadores),
"Curtosis"= kurtosis(dataEstimadores)
)
Resumen=rbind(Resumen,x)
}
## [1] "Para una muestra de tamaño: 5"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 0.9, p-value = 0.000000000000007
## [1] "Para una muestra de tamaño: 10"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 1, p-value = 0.0000000007
## [1] "Para una muestra de tamaño: 15"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 1, p-value = 0.0000001
## [1] "Para una muestra de tamaño: 20"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 1, p-value = 0.000001
## [1] "Para una muestra de tamaño: 30"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 1, p-value = 0.0001
## [1] "Para una muestra de tamaño: 50"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 1, p-value = 0.003
## [1] "Para una muestra de tamaño: 60"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 1, p-value = 0.01
## [1] "Para una muestra de tamaño: 100"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 1, p-value = 0.05
## [1] "Para una muestra de tamaño: 200"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 1, p-value = 0.04
## [1] "Para una muestra de tamaño: 500"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 1, p-value = 0.5Resumen## Muestra TamañoMuestra Media Mediana Desvest Varianza Mín. Máx Asimetría
## 1 0 100 0.50 0.50 0.049 0.00240 0.35 0.62 -0.128
## 2 1 5 0.51 0.40 0.220 0.04847 0.00 1.00 0.096
## 3 2 10 0.51 0.50 0.156 0.02428 0.10 1.00 0.151
## 4 3 15 0.49 0.47 0.127 0.01609 0.13 0.87 0.167
## 5 4 20 0.50 0.50 0.112 0.01260 0.10 0.75 -0.165
## 6 5 30 0.49 0.50 0.093 0.00867 0.27 0.77 0.118
## 7 6 50 0.50 0.50 0.071 0.00500 0.30 0.70 -0.028
## 8 7 60 0.50 0.50 0.063 0.00397 0.30 0.70 -0.045
## 9 8 100 0.50 0.50 0.046 0.00207 0.37 0.62 0.026
## 10 9 200 0.50 0.50 0.032 0.00101 0.41 0.60 0.087
## 11 10 500 0.50 0.50 0.016 0.00027 0.45 0.56 0.031
## Curtosis
## 1 2.9
## 2 2.7
## 3 2.8
## 4 2.9
## 5 2.9
## 6 2.9
## 7 3.1
## 8 2.9
## 9 2.8
## 10 2.8
## 11 3.1Conclusiones
De acuerdo con el anterior ejercicio podemos analizar diferentes puntos:
Convergencia de la media y la mediana: A medida que el tamaño de la muestra aumenta, la media y la mediana de las muestras tienden a aproximarse a 0.50. Esto es coherente con el Teorema del Límite Central, que establece que la media de las muestras converge a la media de la población.
Disminución de la varianza y la desviación estándar: A medida que el tamaño de la muestra aumenta, la varianza y la desviación estándar de las muestras disminuyen. Esto también es coherente con el Teorema del Límite Central, que establece que la varianza de las medias muestrales disminuye a medida que aumenta el tamaño de la muestra.
Asimetría cercana a cero: La asimetría de las muestras tiende a acercarse a cero a medida que el tamaño de la muestra aumenta. Esto sugiere que las muestras tienden a seguir una distribución más simétrica a medida que se vuelven más grandes, lo que es consistente con el Teorema del Límite Central.
Convergencia a la distribución normal: El Teorema del Límite Central implica que, a medida que el tamaño de la muestra aumenta, la distribución de las medias muestrales se aproxima a una distribución normal. Esto se refleja en la disminución de la asimetría y en la estabilización de otras estadísticas a medida que aumenta el tamaño de la muestra.
En resumen, los resultados proporcionados parecen mostrar cómo las estadísticas de las muestras, como la media, la mediana, la varianza y la asimetría, cambian a medida que se aumenta el tamaño de la muestra, y estos cambios son coherentes con las expectativas del Teorema del Límite Central.
- Test de Shapiro: En cuanto al test, podemos ver valores que inician en p = 0.000000000000005 y terminan en p=0.1 a medida que la muestra incrementa. Esto indica que los datos de la muestra inicial de 5 no siguen una distribución normal en absoluto. En este caso, existe evidencia extremadamente fuerte en contra de la hipótesis nula de normalidad.Por otro lado, el valor p de 0.1 en la prueba de Shapiro-Wilk para la muestra de 500 al final, sugiere que los datos de esa muestra tienen menos evidencia en contra de la hipótesis nula de normalidad. Aunque un valor p de 0.1 aún es mayor que el nivel de significancia comúnmente utilizado de 0.05, es menos extremo que el valor muy bajo obtenido en la muestra inicial de 5. Esto podría indicar que, a medida que aumenta el tamaño de la muestra, los datos tienden a seguir una distribución normal con mayor precisión. Es importante recordar que la prueba de normalidad es sensible al tamaño de la muestra, y a medida que aumenta el tamaño de la muestra, es más probable que se obtengan valores p más altos, incluso si los datos tienen algunas desviaciones menores de la normalidad.La decisión de si considerar los datos como aproximadamente normales o no debe basarse en un análisis integral de los resultados de la prueba.
Paso 5: Muestras con diferentes características.
Repita toda la simulación 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.
10% Plantas Enfermas
PoblacionGenerada2 <- simulacion(N, 0.1,0.9)## [1] "Lote Generado de tamaño: 1000 Enfermos: 100 Sanos: 900"for(i in 1:length(muestras)){
dataEstimadores <-Calcular_Iter(PoblacionGenerada2, muestras[i], "Enferma", iter)
par(mfrow=c(1,2))
hist(dataEstimadores, las=1, ylab = "Frecuencia",
main = paste("Muestra de tamaño:", muestras[i]),col = "#2E4053",prob = TRUE)
abline(v=mean(dataEstimadores), col="#2E86C1", lwd=3)
lines(density(dataEstimadores), col = 2, lwd = 2)
qqnorm(dataEstimadores, xlab="Cuantiles teóricos",
ylab="Cuantiles muestrales",main="qq-normalidad")
qqline(dataEstimadores,col = '#D81B60', lwd = 2, lty = 2)
print(paste('Para una muestra de tamaño: ', muestras[i]))
print(shapiro.test(dataEstimadores))
x=data.frame("Muestra"=i,
"TamañoMuestra"=muestras[i],
"Media"=mean(dataEstimadores),
"Mediana"=median(dataEstimadores),
"Desvest"=sd(dataEstimadores),
"Varianza"=var(dataEstimadores),
"Mín."=min(dataEstimadores),
"Máx"=max(dataEstimadores),
"Asimetría"=skewness(dataEstimadores),
"Curtosis"= kurtosis(dataEstimadores)
)
Resumen=rbind(Resumen,x)
}
## [1] "Para una muestra de tamaño: 5"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 0.7, p-value <0.0000000000000002
## [1] "Para una muestra de tamaño: 10"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 0.8, p-value <0.0000000000000002
## [1] "Para una muestra de tamaño: 15"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 0.9, p-value <0.0000000000000002
## [1] "Para una muestra de tamaño: 20"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 0.9, p-value = 0.000000000000002
## [1] "Para una muestra de tamaño: 30"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 1, p-value = 0.000000000008
## [1] "Para una muestra de tamaño: 50"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 1, p-value = 0.000000003
## [1] "Para una muestra de tamaño: 60"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 1, p-value = 0.00000002
## [1] "Para una muestra de tamaño: 100"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 1, p-value = 0.000007
## [1] "Para una muestra de tamaño: 200"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 1, p-value = 0.004
## [1] "Para una muestra de tamaño: 500"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 1, p-value = 0.04Resumen## Muestra TamañoMuestra Media Mediana Desvest Varianza Mín. Máx Asimetría
## 1 0 100 0.498 0.500 0.0490 0.002403 0.350 0.62 -0.128
## 2 1 5 0.505 0.400 0.2202 0.048470 0.000 1.00 0.096
## 3 2 10 0.507 0.500 0.1558 0.024282 0.100 1.00 0.151
## 4 3 15 0.489 0.467 0.1269 0.016093 0.133 0.87 0.167
## 5 4 20 0.499 0.500 0.1123 0.012605 0.100 0.75 -0.165
## 6 5 30 0.494 0.500 0.0931 0.008666 0.267 0.77 0.118
## 7 6 50 0.498 0.500 0.0707 0.004995 0.300 0.70 -0.028
## 8 7 60 0.496 0.500 0.0630 0.003967 0.300 0.70 -0.045
## 9 8 100 0.497 0.500 0.0455 0.002071 0.370 0.62 0.026
## 10 9 200 0.500 0.500 0.0318 0.001011 0.415 0.60 0.087
## 11 10 500 0.499 0.500 0.0164 0.000268 0.454 0.56 0.031
## 12 1 5 0.094 0.000 0.1255 0.015756 0.000 0.60 1.087
## 13 2 10 0.099 0.100 0.0940 0.008837 0.000 0.40 0.759
## 14 3 15 0.101 0.067 0.0770 0.005922 0.000 0.40 0.684
## 15 4 20 0.097 0.100 0.0630 0.003970 0.000 0.30 0.539
## 16 5 30 0.102 0.100 0.0536 0.002869 0.000 0.27 0.460
## 17 6 50 0.103 0.100 0.0434 0.001879 0.000 0.24 0.394
## 18 7 60 0.102 0.100 0.0401 0.001608 0.000 0.28 0.469
## 19 8 100 0.100 0.100 0.0276 0.000761 0.020 0.19 0.257
## 20 9 200 0.098 0.100 0.0184 0.000338 0.045 0.15 0.063
## 21 10 500 0.100 0.100 0.0092 0.000084 0.072 0.13 0.137
## Curtosis
## 1 2.9
## 2 2.7
## 3 2.8
## 4 2.9
## 5 2.9
## 6 2.9
## 7 3.1
## 8 2.9
## 9 2.8
## 10 2.8
## 11 3.1
## 12 3.5
## 13 3.1
## 14 3.3
## 15 3.0
## 16 3.1
## 17 2.8
## 18 3.4
## 19 3.3
## 20 2.7
## 21 3.090% Plantas Enfermas
PoblacionGenerada3 <- simulacion(N, 0.9,0.1)## [1] "Lote Generado de tamaño: 1000 Enfermos: 900 Sanos: 100"for(i in 1:length(muestras)){
dataEstimadores <-Calcular_Iter(PoblacionGenerada3, muestras[i], "Enferma", iter)
par(mfrow=c(1,2))
hist(dataEstimadores, las=1, ylab = "Frecuencia",
main = paste("Muestra de tamaño:", muestras[i]),col = "#2E4053",prob = TRUE)
abline(v=mean(dataEstimadores), col="#2E86C1", lwd=3)
lines(density(dataEstimadores), col = 2, lwd = 2)
qqnorm(dataEstimadores, xlab="Cuantiles teóricos",
ylab="Cuantiles muestrales",main="qq-normalidad")
qqline(dataEstimadores,col = '#D81B60', lwd = 2, lty = 2)
print(paste('Para una muestra de tamaño: ', muestras[i]))
print(shapiro.test(dataEstimadores))
x=data.frame("Muestra"=i,
"TamañoMuestra"=muestras[i],
"Media"=mean(dataEstimadores),
"Mediana"=median(dataEstimadores),
"Desvest"=sd(dataEstimadores),
"Varianza"=var(dataEstimadores),
"Mín."=min(dataEstimadores),
"Máx"=max(dataEstimadores),
"Asimetría"=skewness(dataEstimadores),
"Curtosis"= kurtosis(dataEstimadores)
)
Resumen=rbind(Resumen,x)
}
## [1] "Para una muestra de tamaño: 5"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 0.7, p-value <0.0000000000000002
## [1] "Para una muestra de tamaño: 10"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 0.8, p-value <0.0000000000000002
## [1] "Para una muestra de tamaño: 15"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 0.9, p-value <0.0000000000000002
## [1] "Para una muestra de tamaño: 20"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 0.9, p-value = 0.000000000000003
## [1] "Para una muestra de tamaño: 30"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 0.9, p-value = 0.000000000002
## [1] "Para una muestra de tamaño: 50"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 1, p-value = 0.000000006
## [1] "Para una muestra de tamaño: 60"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 1, p-value = 0.0000002
## [1] "Para una muestra de tamaño: 100"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 1, p-value = 0.00002
## [1] "Para una muestra de tamaño: 200"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 1, p-value = 0.003
## [1] "Para una muestra de tamaño: 500"
##
## Shapiro-Wilk normality test
##
## data: dataEstimadores
## W = 1, p-value = 0.01Resumen## Muestra TamañoMuestra Media Mediana Desvest Varianza Mín. Máx Asimetría
## 1 0 100 0.498 0.500 0.0490 0.002403 0.350 0.62 -0.128
## 2 1 5 0.505 0.400 0.2202 0.048470 0.000 1.00 0.096
## 3 2 10 0.507 0.500 0.1558 0.024282 0.100 1.00 0.151
## 4 3 15 0.489 0.467 0.1269 0.016093 0.133 0.87 0.167
## 5 4 20 0.499 0.500 0.1123 0.012605 0.100 0.75 -0.165
## 6 5 30 0.494 0.500 0.0931 0.008666 0.267 0.77 0.118
## 7 6 50 0.498 0.500 0.0707 0.004995 0.300 0.70 -0.028
## 8 7 60 0.496 0.500 0.0630 0.003967 0.300 0.70 -0.045
## 9 8 100 0.497 0.500 0.0455 0.002071 0.370 0.62 0.026
## 10 9 200 0.500 0.500 0.0318 0.001011 0.415 0.60 0.087
## 11 10 500 0.499 0.500 0.0164 0.000268 0.454 0.56 0.031
## 12 1 5 0.094 0.000 0.1255 0.015756 0.000 0.60 1.087
## 13 2 10 0.099 0.100 0.0940 0.008837 0.000 0.40 0.759
## 14 3 15 0.101 0.067 0.0770 0.005922 0.000 0.40 0.684
## 15 4 20 0.097 0.100 0.0630 0.003970 0.000 0.30 0.539
## 16 5 30 0.102 0.100 0.0536 0.002869 0.000 0.27 0.460
## 17 6 50 0.103 0.100 0.0434 0.001879 0.000 0.24 0.394
## 18 7 60 0.102 0.100 0.0401 0.001608 0.000 0.28 0.469
## 19 8 100 0.100 0.100 0.0276 0.000761 0.020 0.19 0.257
## 20 9 200 0.098 0.100 0.0184 0.000338 0.045 0.15 0.063
## 21 10 500 0.100 0.100 0.0092 0.000084 0.072 0.13 0.137
## 22 1 5 0.906 1.000 0.1311 0.017198 0.400 1.00 -1.199
## 23 2 10 0.897 0.900 0.0944 0.008911 0.500 1.00 -0.807
## 24 3 15 0.898 0.933 0.0767 0.005882 0.667 1.00 -0.481
## 25 4 20 0.899 0.900 0.0674 0.004543 0.600 1.00 -0.604
## 26 5 30 0.897 0.900 0.0515 0.002657 0.700 1.00 -0.511
## 27 6 50 0.900 0.900 0.0410 0.001682 0.740 1.00 -0.447
## 28 7 60 0.899 0.900 0.0384 0.001471 0.767 0.98 -0.295
## 29 8 100 0.899 0.900 0.0294 0.000864 0.790 0.97 -0.308
## 30 9 200 0.901 0.905 0.0186 0.000347 0.825 0.95 -0.109
## 31 10 500 0.900 0.900 0.0089 0.000079 0.872 0.93 -0.188
## Curtosis
## 1 2.9
## 2 2.7
## 3 2.8
## 4 2.9
## 5 2.9
## 6 2.9
## 7 3.1
## 8 2.9
## 9 2.8
## 10 2.8
## 11 3.1
## 12 3.5
## 13 3.1
## 14 3.3
## 15 3.0
## 16 3.1
## 17 2.8
## 18 3.4
## 19 3.3
## 20 2.7
## 21 3.0
## 22 3.8
## 23 3.4
## 24 2.7
## 25 3.4
## 26 3.3
## 27 3.5
## 28 2.8
## 29 3.1
## 30 3.3
## 31 3.1Conclusiones
Realizando el mismo ejercicio pero con características de las muestras diferentes, se ve cómo las estadísticas de las muestras, como la media, la mediana, la varianza y la asimetría, cambian a medida que se aumenta el tamaño de la muestra, y estos cambios son coherentes con las expectativas del Teorema del Límite Central. Pasa igual con el test de Shapiro descrito anteriormente. Lo interesante es que si vemos las gráficas del histograma y el qq-normalidad, la diferencia entre los dos casos es que inician con una distribución sesgada hacia un lado determinado. En el caso de 10% Plantas enfermas, la distribución es sesgada hacia la izquierda, mientras que con el 90% Plantas enfermas, la distribución es sesgada hacia la derecha.
Problema 4: Estimación boostrap
Cuando se extrae una muestra de una población que no es normal y se requiere estimar un intervalo de confianza se pueden utilizar los métodos de estimación bootstrap. Esta metodología supone que se puede reconstruir la población objeto de estudio mediante un muestreo con reemplazo de la muestra que se tiene. Existen varias versiones del método.
Los datos obtenidos para mediciones de eficiencia de combustible son los siguientes: 7.69, 4.97, 4.56, 6.49, 4.34, 6.24 y 4.45. Se supone que es una muestra aleatoria de camiones y que se desea construir un intervalo de confianza del 95 % para la media de la eficiencia de combustible de esta población. No se tiene información de la distribución de los datos. El método bootstrap permite construir intervalos de confianza del 95 %.
Para ilustrar el método suponga que coloca los valores de la muestra en una caja y extrae uno al azar. Este correspondería al primer valor de la muestra bootstrap X∗1. Después de anotado el valor se regresa X∗1 a la caja y se extrae el valor X∗2, regresandolo nuevamente. Este procedimiento se repite hasta completar una muestra de tamaño n, X∗1,X∗2,X∗2,X∗n, conformando la muestra bootstrap.
Es necesario extraer un gran número de muestras (suponga k = 1000). Para cada una de las muestra bootstrap obtenidas se calcula la media X∗i¯, obteniéndose un valor para cada muestra. El intervalo de confianza queda conformado por los percentiles P2.5 y P97.5. Existen dos métodos para estimarlo.
Muestras.
Establecemos un objeto llamado ‘muestras’ que contiene los datos de la eficiencia de combustible de las 7 muestras analizadas.
muestra <- c(7.69, 4.97, 4.56, 6.49, 4.34, 6.24, 4.45)Número de repeticiones bootstrap
k <- 1000Tamaño de la muestra bootstrap
n <- length(muestra)medias_bootstrap <- numeric(k)Método Bootstrap
Teniendo en cuenta las 7 muestras y el número de repeticiones, se hace un bucle que genere muestras y adicionalmente calcule la media de esas muestras.
for (i in 1:k) {
muestra_bootstrap <- sample(muestra, replace = TRUE)
medias_bootstrap[i] <- mean(muestra_bootstrap)
}Se calculan los percentiles inferior y superior
percentil_inferior <- quantile(medias_bootstrap, 0.025)
percentil_superior <- quantile(medias_bootstrap, 0.975)Método 1: (P2.5, P97.5)
Se calcula el intervalo a partir de la definición del método 1.
intervalo_met1 <- c(percentil_inferior, percentil_superior)Método 2: (2X¯ - P97.5, 2X¯ - P2.5)
Se calcula el intervalo a partir de la definición del método 2.
media_muestra <- mean(muestra)
intervalo_met2 <- c(2 * media_muestra - percentil_superior, 2 * media_muestra - percentil_inferior)Resultados:
El Intervalo de confianza para el método 1 es:
cat("Intervalo de confianza (Método 1):", intervalo_met1, "\n")## Intervalo de confianza (Método 1): 4.7 6.4El Intervalo de confianza para el método 2 es:
cat("Intervalo de confianza (Método 2):", intervalo_met2, "\n")## Intervalo de confianza (Método 2): 4.7 6.3