## Warning: package 'knitr' was built under R version 4.3.3
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 se encontrar una aproximación para el valor de π.
¿Cuántos de los puntos están dentro del círculo? ¿Cuál es su estimación de π?
En los resultados evidenciados en la tabla se puede observar una relación directamente proporcional de que a medida que el valor de n crece la cantidad de puntos en el circulo es mayor y por ende el cálculo de π también es mas cercano al valor real.
# Establece los números de puntos a generar
n_valores <- c(1000, 10000, 100000)
# Crear una lista para almacenar los resultados
tabla_resultados <- list()
# Función para estimar π y almacenar los resultados en la tabla
estimar_pi <- function(n) {
# Genera coordenadas x y y aleatorias entre 0 y 1
x <- runif(n, 0, 1)
y <- runif(n, 0, 1)
# Define una función para determinar si un punto está dentro del círculo
dentro_del_circulo <- function(x, y) {
distancia_al_centro_cuadrado <- (x - 0.5)^2 + (y - 0.5)^2
return(distancia_al_centro_cuadrado < 0.25)
}
# Cuenta cuántos puntos están dentro del círculo
puntos_dentro <- sum(dentro_del_circulo(x, y))
# Calcula la estimación de π
estimacion_pi <- 4 * (puntos_dentro / n)
# Almacena los resultados en la tabla
tabla_resultados[[as.character(n)]] <<- c(puntos_dentro = puntos_dentro, estimacion_pi = estimacion_pi)
}
# Realiza la estimación de π para cada valor de n y almacena los resultados en la tabla
for (n in n_valores) {
estimar_pi(n)
}
# Crea una tabla a partir de los resultados
tabla <- t(sapply(tabla_resultados, unlist))
# Agrega los nombres de las filas
rownames(tabla) <- n_valores
# Agrega los nombres de las columnas
colnames(tabla) <- c("Puntos dentro del círculo", "Estimación de π")
# Imprime la tabla con kable
knitr::kable(tabla, caption = "Resultados de la estimación de π para diferentes valores de n")| Puntos dentro del círculo | Estimación de π | |
|---|---|---|
| 1000 | 781 | 3.12400 |
| 10000 | 7928 | 3.17120 |
| 1e+05 | 78333 | 3.13332 |
La simulación ayuda a entender y validad las propiedades de los estimadores estadísticos como son. insesgadez, eficiencia y la consistencia principalmente. El siguiente problema permite evidenciar las principales características de un grupo de estimadores propuestos para la estimación de un parámetro asociado a un modelo de probabilidad.
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:
• θ1ˆ= (X1+X2 / 6) + (X3+X4 / 3)
• θ2ˆ= (X1+2X2+3X3+4X4) / 5
• θ3ˆ= ((X1+X2+X3+X4) / 4)
• θ4ˆ= ((min{X1,X2,X3,X4} ) + (max{X1,X2,X3,X4})) / 2
###### Problema N.2 ######
#crear los estimadores
theta_uno = function(x) {((x[1]+x[2])/6) +((x[3]+x[4])/3)}
theta_dos = function(x) {(x[1]+(2*x[2])+(3*x[3])+(4*x[4]))/5}
theta_tres = function(x) {(x[1]+x[2]+x[3]+x[4])/4}
theta_cuatro = function(x) {(min(x[1]+x[2]+x[3]+x[4])+max(x[1]+x[2]+x[3]+x[4]))/2}
#crear muestra para cada estimador
set.seed(123) # Fijar la semilla para reproducibilidad
n = 4 # tamaño de la muestra
theta = 2 # valor asumido de tetha
#Funcion de generacion de estimador, donde se introduce la funcion estimadores
generar_estimadores = function(funcion1,funcion2, funcion3, funcion4){
x = rexp(4,rate = 2)
estimador_uno = funcion1(x)
estimador_dos = funcion2(x)
estimador_tres = funcion3(x)
estimador_cuatro = funcion4(x)
estimadores = c(estimador_uno, estimador_dos, estimador_tres, estimador_cuatro)
return(estimadores)
}
# generacion de muestras n 20, 50, 100, 1000
#codigo para repetir vente veces los estimadores
estimadores_n20 = list()
for (i in 1:20) {
estimadores_n20[[i]]= generar_estimadores(theta_uno,theta_dos,theta_tres,theta_cuatro)
}
# Convertir la lista de resultados en un dataframe
resultados_df20 <- as.data.frame(do.call(rbind, estimadores_n20))
#codigo para repetir cincuenta veces los estimadores
estimadores_n50 = list()
for (i in 1:50) {
estimadores_n50[[i]]= generar_estimadores(theta_uno,theta_dos,theta_tres,theta_cuatro)
}
# Convertir la lista de resultados en un dataframe
resultados_df50 <- as.data.frame(do.call(rbind, estimadores_n50))
#codigo para repetir cien veces los estimadores
estimadores_n100 = list()
for (i in 1:100) {
estimadores_n100[[i]]= generar_estimadores(theta_uno,theta_dos,theta_tres,theta_cuatro)
}
# Convertir la lista de resultados en un dataframe
resultados_df100 <- as.data.frame(do.call(rbind, estimadores_n100))
#codigo para repetir mil veces los estimadores
estimadores_n1000 = list()
for (i in 1:1000) {
estimadores_n1000[[i]]= generar_estimadores(theta_uno,theta_dos,theta_tres,theta_cuatro)
}
# Convertir la lista de resultados en un dataframe
resultados_df1000 <- as.data.frame(do.call(rbind, estimadores_n1000))
# Graficar las muestras obtenidas por cada estimador
par(mfrow=c(2, 2)) # Organizar las gráficas en una cuadrícula de 2x2
# Graficar los estimadores usando un diagrama de caja n=20
boxplot(resultados_df20,
main = "Estimadores para Muestras de tamaño 20",
xlab = "Estimador",
ylab = "Valor Estimado",
col = c("skyblue", "orange", "green", "purple"),
names = c("Uno", "Dos", "Tres", "Cuatro")) # Colores personalizados para cada estimador
# Agregar línea horizontal
abline(h = 2, col = "red")
# Graficar los estimadores usando un diagrama de caja n=50
boxplot(resultados_df50,
main = "Estimadores para Muestras de tamaño 50",
xlab = "Estimador",
ylab = "Valor Estimado",
col = c("skyblue", "orange", "green", "purple"),
names = c("Uno", "Dos", "Tres", "Cuatro")) # Colores personalizados para cada estimador
# Agregar línea horizontal
abline(h = 2, col = "red")
# Graficar los estimadores usando un diagrama de caja n=100
boxplot(resultados_df100,
main = "Estimadores para Muestras de tamaño 100",
xlab = "Estimador",
ylab = "Valor Estimado",
col = c("skyblue", "orange", "green", "purple"),
names = c("Uno", "Dos", "Tres", "Cuatro")) # Colores personalizados para cada estimador
# Agregar línea horizontal
abline(h = 2, col = "red")
# Graficar los estimadores usando un diagrama de caja n=1000
boxplot(resultados_df1000,
main = "Estimadores para Muestras de tamaño 1000",
xlab = "Estimador",
ylab = "Valor Estimado",
col = c("skyblue", "orange", "green", "purple"),
names = c("Uno", "Dos", "Tres", "Cuatro")) # Colores personalizados para cada estimador
# Agregar línea horizontal
abline(h = 2, col = "red")Al comparar los cuatro estimadores propuesto para un parámetro theta desconocido, con diferente tamaño de muestreo, como se visualiza en la grafica anterior, donde la línea roja nos indica el valor asumido de theta, se puede analizar que el estimador número cuatro es el mas cercano para estimar este parámetro, al tener una media cercana a 2, esto se observa en todos los diferentes tamaños de la muestra, pero también se puede ver que tiene una variabilidad alta en todas las muestras comparado con los otros estimadores que presentan desviaciones estándar muy por debajo de 1. Tan solo en la muestra de tamaño 1000 la variabilidad del estimador es baja indicando que a medida que la muestra crece el estimador se va haciendo consistente.
#propiedades
# verificar la propiedad de insesgacion en los estimadores con n=20
medias_20 = colMeans(resultados_df20)
#medias_20
kable(medias_20, caption = "Propiedades del estimador Uno para diferentes tamaños de muestra")| x | |
|---|---|
| V1 | 0.4917774 |
| V2 | 0.9888430 |
| V3 | 0.5085644 |
| V4 | 2.0342578 |
# verificar la propiedad de insesgacion en los estimadores con n=50
medias_50 = colMeans(resultados_df50)
#medias_50
# Tabla para el estimador dos
kable(medias_50, caption = "Propiedades del estimador Dos para diferentes tamaños de muestra")| x | |
|---|---|
| V1 | 0.5106887 |
| V2 | 1.0371388 |
| V3 | 0.5107465 |
| V4 | 2.0429858 |
# verificar la propiedad de insesgacion en los estimadores con n=100
medias_100 = colMeans(resultados_df100)
#medias_100
# Tabla para el estimador tres
kable(medias_100, caption = "Propiedades del estimador Tres para diferentes tamaños de muestra")| x | |
|---|---|
| V1 | 0.4900755 |
| V2 | 0.9808880 |
| V3 | 0.4982349 |
| V4 | 1.9929397 |
# verificar la propiedad de insesgacion en los estimadores con n=1000
medias_1000 = colMeans(resultados_df1000)
#medias_1000
# Tabla para el estimador cuatro
kable(medias_1000, caption = "Propiedades del estimador Cuatro para diferentes tamaños de muestra")| x | |
|---|---|
| V1 | 0.5041606 |
| V2 | 1.0115020 |
| V3 | 0.5025265 |
| V4 | 2.0101062 |
par(mfrow=c(2, 2)) # Organizar las gráficas en una cuadrícula de 2x2
#Graficar media para observar insesgacion
barplot(medias_20,
main = "Media de los estimadores, n=20",
ylab = "Media",
xlab = "Estimadores",
col = "blue",
ylim = range(c(0,3)),
names = c("Uno", "Dos", "Tres", "Cuatro"))
barplot(medias_50,
main = "Media de los estimadores, n=50",
ylab = "Media",
xlab = "Estimadores",
col = "blue",
ylim = range(c(0,3)),
names = c("Uno", "Dos", "Tres", "Cuatro"))
barplot(medias_100,
main = "Media de los estimadores, n=100",
ylab = "Media",
xlab = "Estimadores",
col = "blue",
ylim = range(c(0,3)),
names = c("Uno", "Dos", "Tres", "Cuatro"))
barplot(medias_1000,
main = "Media de los estimadores, n=1000",
ylab = "Media",
xlab = "Estimadores",
col = "blue",
ylim = range(c(0,3)),
names = c("Uno", "Dos", "Tres", "Cuatro"))Según la gráfica de comparación de medias en diferentes tamaños de muestra para todos los estimadores, el estimador con mayor insesgadez es cuatro, donde su media esta muy cerca del parámetro a estimar, que en este caso es 2.
#verificar la propiedad eficiente en los estimadores con la desviacion estandar
# calculo de eficiencia para n=20
sd_20 = apply(resultados_df20,2,sd)
# calculo de eficiencia para n=50
sd_50 = apply(resultados_df50,2,sd)
# calculo de eficiencia para n=100
sd_100 = apply(resultados_df100,2,sd)
# calculo de eficiencia para n=1000
sd_1000 = apply(resultados_df1000,2,sd)
# Tabla para el estimador uno
kable(sd_20, caption = "Propiedades del estimador Uno para diferentes tamaños de muestra")| x | |
|---|---|
| V1 | 0.2254758 |
| V2 | 0.5026344 |
| V3 | 0.2474003 |
| V4 | 0.9896012 |
# Tabla para el estimador dos
kable(sd_50, caption = "Propiedades del estimador Dos para diferentes tamaños de muestra")| x | |
|---|---|
| V1 | 0.2162134 |
| V2 | 0.4488974 |
| V3 | 0.1961480 |
| V4 | 0.7845919 |
# Tabla para el estimador tres
kable(sd_100, caption = "Propiedades del estimador Tres para diferentes tamaños de muestra")| x | |
|---|---|
| V1 | 0.2331738 |
| V2 | 0.4966547 |
| V3 | 0.2281857 |
| V4 | 0.9127430 |
# Tabla para el estimador cuatro
kable(sd_1000, caption = "Propiedades del estimador Cuatro para diferentes tamaños de muestra")| x | |
|---|---|
| V1 | 0.2734736 |
| V2 | 0.5704089 |
| V3 | 0.2502165 |
| V4 | 1.0008661 |
par(mfrow=c(2, 2)) # Organizar las gráficas en una cuadrícula de 2x2
#grafico de la desviacion estandar para determinar la propiedad de eficiencia
barplot(sd_20,
main = "Desviación Estándar para n = 20",
xlab = "Estimador",
ylab = "Desviación Estándar",
names = c("Uno", "Dos", "Tres", "Cuatro"))
barplot(sd_50,
main = "Desviación Estándar para n = 50",
xlab = "Estimador",
ylab = "Desviación Estándar",
names = c("Uno", "Dos", "Tres", "Cuatro"))
barplot(sd_100,
main = "Desviación Estándar para n = 100",
xlab = "Estimador",
ylab = "Desviación Estándar",
names = c("Uno", "Dos", "Tres", "Cuatro"))
barplot(sd_1000,
main = "Desviación Estándar para n = 1000",
xlab = "Estimador",
ylab = "Desviación Estándar",
names = c("Uno", "Dos", "Tres", "Cuatro"))Para la propiedad de eficiencia, se utilizó la desviación estándar para saber su variabilidad, como se observa en la gráfica, se muestra que el estimador cuatro quien presento mayor insesgadez posee una alta variabilidad, deduciendo que tiene una menor eficiencia comparado a otros estimadores. En conclusión, pasando los estimadores por todas las propiedades se obtiene que el estimador cuatro es un estimador INSESGADO CON BAJA EFICIENCIA, PERO CONSISTENTE, a medida que el tamaño de las muestras crece.
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.
A continuación se describen los siguientes pasos para su verificació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%.
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.
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.
Repita los puntos b y c 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.
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.
##### PROBLEMA 3 ######
#simulacion de la poblacion
N = 1000 # tamaño de la población
enfermas = 500 # plantas enfermas
#vector con estado de cada planta (1 para enfermo, 0 para sano)
estado_plantas = c(rep(1, enfermas), rep(0, N - enfermas))
poblacion = sample(estado_plantas) # Barajar la muestra aleatoria de la población
poblacion # poblacion completa## [1] 0 0 1 0 1 1 0 1 1 1 0 0 1 1 0 1 0 1 0 1 0 0 0 1 1 0 1 0 1 1 0 0 1 1 0 0 0
## [38] 0 1 0 1 1 0 1 0 1 1 0 1 1 1 0 1 0 1 0 0 1 0 0 1 0 1 0 1 0 1 1 1 0 1 1 1 1
## [75] 0 1 1 0 1 1 1 1 0 1 0 0 0 1 1 1 0 0 0 1 0 0 1 0 0 0 1 0 0 1 0 0 1 1 1 0 0
## [112] 1 0 1 1 0 0 0 0 1 1 1 1 1 0 0 0 0 1 0 0 0 0 1 1 1 0 1 0 1 0 0 0 1 0 0 1 1
## [149] 1 0 1 1 1 0 0 0 0 1 1 0 0 0 1 0 1 0 0 0 0 0 0 0 1 1 0 1 1 1 0 0 1 0 1 1 0
## [186] 1 0 0 1 0 1 1 1 0 0 0 0 1 0 0 1 1 1 0 0 0 0 0 1 0 0 0 0 1 0 0 1 0 0 1 0 1
## [223] 1 1 1 1 0 0 1 1 1 1 1 1 1 1 0 0 1 0 1 0 1 1 0 1 0 0 1 0 0 0 1 1 0 0 1 0 1
## [260] 1 0 0 0 0 0 1 1 0 1 0 1 1 1 0 1 0 1 1 0 1 1 1 1 1 1 0 0 0 0 0 1 0 1 0 0 0
## [297] 1 1 1 1 1 0 0 0 1 0 0 0 0 1 0 0 0 1 1 1 1 0 1 1 0 0 1 0 1 0 1 0 0 0 1 1 0
## [334] 0 0 0 1 0 0 0 1 0 1 0 0 0 1 0 1 1 0 0 1 1 0 0 1 1 1 1 0 0 1 1 1 0 1 1 1 1
## [371] 1 1 0 1 1 1 0 0 1 0 0 0 0 1 1 1 1 1 0 0 1 1 1 1 0 1 1 0 0 1 1 0 0 1 0 1 0
## [408] 1 0 0 0 1 0 1 0 1 0 1 1 1 1 0 1 1 0 1 0 1 1 0 1 1 0 1 1 1 1 1 0 0 1 0 0 0
## [445] 1 0 0 0 0 0 0 0 1 0 1 1 0 0 1 0 1 1 1 1 0 1 0 1 1 1 1 1 0 0 1 0 0 1 1 1 1
## [482] 1 1 0 1 1 0 1 0 0 0 0 0 1 0 1 1 0 0 1 0 1 1 0 1 1 0 1 1 0 1 1 0 0 1 0 0 1
## [519] 0 0 0 1 1 1 1 0 0 0 1 0 1 1 0 1 1 0 1 1 1 1 1 1 0 1 0 0 1 1 1 0 0 1 1 1 0
## [556] 0 1 0 1 1 0 0 1 0 1 1 1 0 1 0 1 1 0 1 1 1 0 0 1 1 0 0 0 1 1 1 0 0 0 1 1 0
## [593] 0 1 1 0 0 0 1 1 1 0 0 1 1 1 1 0 1 1 1 0 1 0 1 0 1 1 0 0 0 0 1 1 0 1 1 1 1
## [630] 0 0 1 0 0 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 1 1 0 0 0 0 1 1 1 0 1 0 1 0 0 0 0
## [667] 1 1 1 0 1 0 0 1 1 1 0 1 0 1 0 0 1 0 0 1 1 1 1 0 0 0 0 1 1 1 0 0 0 1 0 0 0
## [704] 1 0 1 1 1 1 0 0 0 0 0 1 1 1 0 1 0 0 0 1 1 1 0 1 0 0 0 0 0 1 0 0 1 1 1 1 1
## [741] 1 0 0 0 0 0 0 1 0 1 1 1 1 0 0 0 1 1 0 0 0 1 0 0 0 0 1 1 1 1 1 0 1 1 1 1 0
## [778] 1 0 0 1 0 1 1 1 1 1 0 1 0 1 0 1 1 1 0 0 1 0 0 0 1 1 0 1 1 0 1 1 1 0 1 1 0
## [815] 0 0 0 1 1 1 1 1 0 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 0 0 0 0 0 1 1 1 1 0 1 1 1
## [852] 0 0 0 1 0 1 0 0 0 0 0 1 1 1 0 0 0 0 0 0 1 1 1 0 1 0 0 0 0 0 0 0 1 1 0 1 1
## [889] 1 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 1 1 0 1 1 1 1 1 0 0 1 0 0 1 1
## [926] 1 1 0 0 0 1 0 0 1 1 0 1 0 1 1 1 1 1 1 1 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0
## [963] 1 0 1 0 0 0 0 1 1 0 0 1 0 1 0 0 1 0 0 1 1 1 1 0 0 0 0 0 0 1 0 0 0 0 0 1 1
## [1000] 0
#Funcion de generacion de proporcion muestral para muestra aleatoria
p_muestral = function(poblacion, tamaño){
muestra_aleatoria = sample(poblacion, tamaño) #obtener muestra
plant_enfer_muestra = sum(muestra_aleatoria)
prop_muestral=plant_enfer_muestra/length(muestra_aleatoria)
return(prop_muestral)
}
tamaño = 25
# Funcion para repetir el codigo para repetir 500 veces proporcion muestral
generacion_repeticiones = function(poblacion, tamaño, repeticiones){
p_muestral_rep = vector("list", length = repeticiones)
for (i in 1:repeticiones) {
p_muestral_rep[[i]]= p_muestral(poblacion, tamaño)
}
resultados_df_rep = as.data.frame(do.call(rbind, p_muestral_rep))
return(resultados_df_rep)
}
repeticiones500_25= generacion_repeticiones(poblacion,tamaño = 25,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_25 = mean(repeticiones500_25$V1)
ds_pm500_25 = sd(repeticiones500_25$V1)
resumen = summarytools::descr(repeticiones500_25$V1)
resumen## Descriptive Statistics
## repeticiones500_25$V1
## N: 500
##
## V1
## ----------------- --------
## Mean 0.51
## Std.Dev 0.10
## Min 0.20
## Q1 0.44
## Median 0.52
## Q3 0.56
## Max 0.80
## MAD 0.06
## IQR 0.12
## CV 0.19
## Skewness -0.14
## SE.Skewness 0.11
## Kurtosis 0.19
## N.Valid 500.00
## Pct.Valid 100.00
# Graficar las muestras obtenidas por cada estimador
par(mfrow=c(1, 1)) # Organizar las gráficas en una cuadrícula de 1x1
# Visualizar los resultados
hist(repeticiones500_25$V1,
main = "Histograma de Estimadores de Proporción",
xlab = "Estimador de Proporción",
ylab = "Frecuencia")boxplot(repeticiones500_25$V1,
main = "Proporcion muestral Plantas Enfermas, n=25, 500 rep",
xlab = "Repeticiones",
ylab = "Porcentaje",
col = "green",
names = "Plantas enfermas")En la grafica anterior se puede observar la proporción muestral para el parámetro de plantas enfermas de una muestra de tamaño 25 con 500 repeticiones, donde se puede analizar que se tiene un estimador insesgado cuya media de 0.51 se acerca bastante al valor del parámetro que en este caso es 0.5, cuenta con una desviación estándar de 0.1 indicando que tiene baja variabilidad
# Repeticion para diferentes muestras
# Tamaño de muestra = 5
repeticiones500_5= generacion_repeticiones(poblacion,tamaño=5,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_5 = mean(repeticiones500_5$V1)
ds_pm500_5 = sd(repeticiones500_5$V1)
# Tamaño de muestra = 10
repeticiones500_10= generacion_repeticiones(poblacion,tamaño=10,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_10 = mean(repeticiones500_10$V1)
ds_pm500_10 = sd(repeticiones500_10$V1)
# Tamaño de muestra = 15
repeticiones500_15= generacion_repeticiones(poblacion,tamaño=15,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_15 = mean(repeticiones500_15$V1)
ds_pm500_15 = sd(repeticiones500_15$V1)
# Tamaño de muestra = 20
repeticiones500_20= generacion_repeticiones(poblacion,tamaño=20,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_20 = mean(repeticiones500_20$V1)
ds_pm500_20 = sd(repeticiones500_20$V1)
# Tamaño de muestra = 30
repeticiones500_30= generacion_repeticiones(poblacion,tamaño=30,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_30 = mean(repeticiones500_30$V1)
ds_pm500_30 = sd(repeticiones500_30$V1)
# Tamaño de muestra = 50
repeticiones500_50= generacion_repeticiones(poblacion,tamaño=50,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_50 = mean(repeticiones500_50$V1)
ds_pm500_50 = sd(repeticiones500_50$V1)
# Tamaño de muestra = 60
repeticiones500_60= generacion_repeticiones(poblacion,tamaño=60,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_60 = mean(repeticiones500_60$V1)
ds_pm500_60 = sd(repeticiones500_60$V1)
# Tamaño de muestra = 100
repeticiones500_100= generacion_repeticiones(poblacion,tamaño=100,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_100 = mean(repeticiones500_100$V1)
ds_pm500_100 = sd(repeticiones500_100$V1)
# Tamaño de muestra = 200
repeticiones500_200= generacion_repeticiones(poblacion,tamaño=200,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_200 = mean(repeticiones500_200$V1)
ds_pm500_200 = sd(repeticiones500_200$V1)
# Tamaño de muestra = 500
repeticiones500_500= generacion_repeticiones(poblacion,tamaño=500,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_500 = mean(repeticiones500_500$V1)
ds_pm500_500 = sd(repeticiones500_500$V1)
### meter en tabla
# Test de normalidad de Shapiro-Wilk
#shapiro.test(repeticiones500_25$V1)
#shapiro.test(repeticiones500_5$V1)
#shapiro.test(repeticiones500_10$V1)
#shapiro.test(repeticiones500_15$V1)
#shapiro.test(repeticiones500_20$V1)
#shapiro.test(repeticiones500_30$V1)
#shapiro.test(repeticiones500_50$V1)
#shapiro.test(repeticiones500_60$V1)
#shapiro.test(repeticiones500_100$V1)
#shapiro.test(repeticiones500_200$V1)
#shapiro.test(repeticiones500_500$V1)
# Realizar las pruebas de Shapiro-Wilk
results <- data.frame(
Sample_Size = c(25, 5, 10, 15, 20, 30, 50, 60, 100, 200, 500),
p_value = c(
shapiro.test(repeticiones500_25$V1)$p.value,
shapiro.test(repeticiones500_5$V1)$p.value,
shapiro.test(repeticiones500_10$V1)$p.value,
shapiro.test(repeticiones500_15$V1)$p.value,
shapiro.test(repeticiones500_20$V1)$p.value,
shapiro.test(repeticiones500_30$V1)$p.value,
shapiro.test(repeticiones500_50$V1)$p.value,
shapiro.test(repeticiones500_60$V1)$p.value,
shapiro.test(repeticiones500_100$V1)$p.value,
shapiro.test(repeticiones500_200$V1)$p.value,
shapiro.test(repeticiones500_500$V1)$p.value
)
)
# Aplicar formato a los valores p en notación científica si son muy pequeños
results$p_value <- format(results$p_value, scientific = TRUE)
# Tabla con los resultados
kable(results, caption = "Resultados de la prueba de Shapiro-Wilk para diferentes tamaños de muestra",digits = 4,col.width = "2cm",aling = "c")| Sample_Size | p_value |
|---|---|
| 25 | 6.300701e-06 |
| 5 | 3.072503e-15 |
| 10 | 1.109604e-09 |
| 15 | 1.797807e-07 |
| 20 | 4.053972e-06 |
| 30 | 1.505398e-04 |
| 50 | 4.616669e-04 |
| 60 | 7.755257e-04 |
| 100 | 1.532311e-02 |
| 200 | 3.833420e-01 |
| 500 | 2.650699e-01 |
par(mfrow=c(2, 4)) # Organizar las gráficas en una cuadrícula de 1x1
# grafico de histograma
#hist(repeticiones500_25$V1, main = "n = 25", freq=FALSE, xlab = "Repeticiones(500) n = 25")
hist(repeticiones500_5$V1, main = "n = 5", freq=FALSE, xlab = "Repeticiones(500) n = 5")
hist(repeticiones500_10$V1, main = "n = 10", freq=FALSE, xlab = "Repeticiones(500) n = 10")
#hist(repeticiones500_15$V1, main = "n = 15", freq=FALSE, xlab = "Repeticiones(500) n = 15")
hist(repeticiones500_20$V1, main = "n = 20", freq=FALSE, xlab = "Repeticiones(500) n = 20")
hist(repeticiones500_30$V1, main = "n = 30", freq=FALSE, xlab = "Repeticiones(500) n = 30")
hist(repeticiones500_50$V1, main = "n = 50", freq=FALSE, xlab = "Repeticiones(500) n = 50")
#hist(repeticiones500_60$V1, main = "n = 60", freq=FALSE, xlab = "Repeticiones(500) n = 60")
hist(repeticiones500_100$V1, main = "n = 100", freq=FALSE, xlab = "Repeticiones(500) n = 100")
hist(repeticiones500_200$V1, main = "n = 200", freq=FALSE, xlab = "Repeticiones(500) n = 200")
hist(repeticiones500_500$V1, main = "n = 500", freq=FALSE, xlab = "Repeticiones(500) n = 500")A medida que se va aumentando el tamaño de muestra, en los gráficos anteriores visualmente se observa que hay una convergencia de la distribución de la media muestral a una distribución normal, pero el test de shapiro wilks indica que solo hasta el tamaño de muestra n=200, bajo un nivel de significancia de 0.05, el valor de p supera este umbral y revela que no hay evidencia suficiente para rechazar la hipótesis nula de que los datos siguen una distribución normal con una muestra mayor a 30.
Para el caso de la estimación de la proporción muestral para plantas enfermas con un valor asumido de 10% en una población de 1000 individuos.
######Simulacion para una poblacion N=1000 con 10% plantas enfermas
N = 1000 # tamaño de la población
enfermas10 = 100 # plantas enfermas
#vector con estado de cada planta (1 para enfermo, 0 para sano)
estado_plantas10 = c(rep(1, enfermas10), rep(0, N - enfermas10))
poblacion10E = sample(estado_plantas10) # Barajar la muestra aleatoria de la población
#poblacion10E # poblacion completa
# Repeticion para diferentes muestras
# Tamaño de muestra = 1
repeticiones500_1_10= generacion_repeticiones(poblacion10E,tamaño=1,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_1_10 = mean(repeticiones500_1_10$V1)
ds_pm500_1_10 = sd(repeticiones500_1_10$V1)
# Tamaño de muestra = 5
repeticiones500_5_10= generacion_repeticiones(poblacion10E,tamaño=5,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_5_10 = mean(repeticiones500_5_10$V1)
ds_pm500_5_10 = sd(repeticiones500_5_10$V1)
# Tamaño de muestra = 10
repeticiones500_10_10= generacion_repeticiones(poblacion10E,tamaño=10,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_10_10 = mean(repeticiones500_10_10$V1)
ds_pm500_10_10 = sd(repeticiones500_10_10$V1)
# Tamaño de muestra = 15
repeticiones500_15_10= generacion_repeticiones(poblacion10E,tamaño=15,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_15_10 = mean(repeticiones500_15_10$V1)
ds_pm500_15_10 = sd(repeticiones500_15_10$V1)
# Tamaño de muestra = 20
repeticiones500_20_10= generacion_repeticiones(poblacion10E,tamaño=20,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_20_10 = mean(repeticiones500_20_10$V1)
ds_pm500_20_10 = sd(repeticiones500_20_10$V1)
# Tamaño de muestra = 30
repeticiones500_30_10= generacion_repeticiones(poblacion10E,tamaño=30,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_30_10 = mean(repeticiones500_30_10$V1)
ds_pm500_30_10 = sd(repeticiones500_30_10$V1)
# Tamaño de muestra = 50
repeticiones500_50_10= generacion_repeticiones(poblacion10E,tamaño=50,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_50_10 = mean(repeticiones500_50_10$V1)
ds_pm500_50_10 = sd(repeticiones500_50_10$V1)
# Tamaño de muestra = 60
repeticiones500_60_10= generacion_repeticiones(poblacion10E,tamaño=60,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_60_10 = mean(repeticiones500_60_10$V1)
ds_pm500_60_10 = sd(repeticiones500_60_10$V1)
# Tamaño de muestra = 100
repeticiones500_100_10= generacion_repeticiones(poblacion10E,tamaño=100,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_100_10 = mean(repeticiones500_100_10$V1)
ds_pm500_100_10 = sd(repeticiones500_100_10$V1)
# Tamaño de muestra = 200
repeticiones500_200_10= generacion_repeticiones(poblacion10E,tamaño=200,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_200_10 = mean(repeticiones500_200_10$V1)
ds_pm500_200_10 = sd(repeticiones500_200_10$V1)
# Tamaño de muestra = 500
repeticiones500_500_10= generacion_repeticiones(poblacion10E,tamaño=500,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_500_10 = mean(repeticiones500_500_10$V1)
ds_pm500_500_10 = sd(repeticiones500_500_10$V1)
# Test de normalidad de Shapiro-Wilk
#shapiro.test(repeticiones500_1_10$V1)
#shapiro.test(repeticiones500_5_10$V1)
#shapiro.test(repeticiones500_10_10$V1)
#shapiro.test(repeticiones500_15_10$V1)
#shapiro.test(repeticiones500_20_10$V1)
#shapiro.test(repeticiones500_30_10$V1)
#shapiro.test(repeticiones500_50_10$V1)
#shapiro.test(repeticiones500_60_10$V1)
#shapiro.test(repeticiones500_100_10$V1)
#shapiro.test(repeticiones500_200_10$V1)
#shapiro.test(repeticiones500_500_10$V1)
# Realizar las pruebas de Shapiro-Wilk para los nuevos datos
results_nuevos <- data.frame(
Sample_Size = c(1, 5, 10, 15, 20, 30, 50, 60, 100, 200, 500),
p_value_nuevos = c(
shapiro.test(repeticiones500_1_10$V1)$p.value,
shapiro.test(repeticiones500_5_10$V1)$p.value,
shapiro.test(repeticiones500_10_10$V1)$p.value,
shapiro.test(repeticiones500_15_10$V1)$p.value,
shapiro.test(repeticiones500_20_10$V1)$p.value,
shapiro.test(repeticiones500_30_10$V1)$p.value,
shapiro.test(repeticiones500_50_10$V1)$p.value,
shapiro.test(repeticiones500_60_10$V1)$p.value,
shapiro.test(repeticiones500_100_10$V1)$p.value,
shapiro.test(repeticiones500_200_10$V1)$p.value,
shapiro.test(repeticiones500_500_10$V1)$p.value
)
)
# Aplicar formato a los valores p en notación científica si son muy pequeños
results_nuevos$p_value_nuevos <- format(results_nuevos$p_value_nuevos, scientific = TRUE)
# Tabla con los resultados de los nuevos datos
kable(results_nuevos, caption = "Resultados de la prueba de Shapiro-Wilk para diferentes tamaños de muestra (nuevos datos)", digits = 4, col.width = "2cm")| Sample_Size | p_value_nuevos |
|---|---|
| 1 | 4.016453e-38 |
| 5 | 2.150474e-28 |
| 10 | 4.276654e-23 |
| 15 | 2.821784e-18 |
| 20 | 1.455660e-14 |
| 30 | 1.407972e-12 |
| 50 | 4.612953e-10 |
| 60 | 7.310005e-08 |
| 100 | 1.814241e-06 |
| 200 | 4.619010e-03 |
| 500 | 1.051969e-01 |
par(mfrow=c(2, 4)) # Organizar las gráficas en una cuadrícula de 1x1
# grafico de histograma
hist(repeticiones500_1_10$V1, main = "n = 1", freq=FALSE, xlab = "Repeticiones(500) n = 1")
hist(repeticiones500_5_10$V1, main = "n = 5", freq=FALSE, xlab = "Repeticiones(500) n = 5")
hist(repeticiones500_10_10$V1, main = "n = 10", freq=FALSE, xlab = "Repeticiones(500) n = 10")
#hist(repeticiones500_15_10$V1, main = "n = 15", freq=FALSE, xlab = "Repeticiones(500) n = 15")
#hist(repeticiones500_20_10$V1, main = "n = 20", freq=FALSE, xlab = "Repeticiones(500) n = 20")
hist(repeticiones500_30_10$V1, main = "n = 30", freq=FALSE, xlab = "Repeticiones(500) n = 30")
hist(repeticiones500_50_10$V1, main = "n = 50", freq=FALSE, xlab = "Repeticiones(500) n = 50")
#hist(repeticiones500_60_10$V1, main = "n = 60", freq=FALSE, xlab = "Repeticiones(500) n = 60")
hist(repeticiones500_100_10$V1, main = "n = 100", freq=FALSE, xlab = "Repeticiones(500) n = 100")
hist(repeticiones500_200_10$V1, main = "n = 200", freq=FALSE, xlab = "Repeticiones(500) n = 200")
hist(repeticiones500_500_10$V1, main = "n = 500", freq=FALSE, xlab = "Repeticiones(500) n = 500")Para este caso se tiene que solo en la muestra de tamaño 500 se puede decir que los datos siguen una distribución normal, y se apoya en una media cercana al parámetro 0.1 y una desviación estándar de 0.009, con un resultado de p value = 0.1052 en el test de shapiro wilks.
Para el caso de la estimación de la proporción muestral para plantas enfermas con un valor asumido de 90% en una población de 1000 individuos
####### Simulacion para una poblacion N=1000 con 90% plantas enfermas
N = 1000 # tamaño de la población
enfermas90 = 900 # plantas enfermas
#vector con estado de cada planta (1 para enfermo, 0 para sano)
estado_plantas90 = c(rep(1, enfermas90), rep(0, N - enfermas90))
poblacion90E = sample(estado_plantas90) # Barajar la muestra aleatoria de la población
#poblacion90E # poblacion completa
# Repeticion para diferentes muestras
# Tamaño de muestra = 1
repeticiones500_1_90= generacion_repeticiones(poblacion90E,tamaño=1,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_1_90 = mean(repeticiones500_1_90$V1)
ds_pm500_1_90 = sd(repeticiones500_1_90$V1)
# Tamaño de muestra = 5
repeticiones500_5_90 = generacion_repeticiones(poblacion90E,tamaño=5,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_5_90 = mean(repeticiones500_5_90$V1)
ds_pm500_5_90 = sd(repeticiones500_5_90$V1)
# Tamaño de muestra = 10
repeticiones500_10_90= generacion_repeticiones(poblacion90E,tamaño=10,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_10_90 = mean(repeticiones500_10_90$V1)
ds_pm500_10_90 = sd(repeticiones500_10_90$V1)
# Tamaño de muestra = 15
repeticiones500_15_90= generacion_repeticiones(poblacion90E,tamaño=15,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_15_90 = mean(repeticiones500_15_90$V1)
ds_pm500_15_90 = sd(repeticiones500_15_90$V1)
# Tamaño de muestra = 20
repeticiones500_20_90= generacion_repeticiones(poblacion90E,tamaño=20,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_20_90 = mean(repeticiones500_20_90$V1)
ds_pm500_20_90 = sd(repeticiones500_20_90$V1)
# Tamaño de muestra = 30
repeticiones500_30_90= generacion_repeticiones(poblacion90E,tamaño=30,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_30_90 = mean(repeticiones500_30_90$V1)
ds_pm500_30_90 = sd(repeticiones500_30_90$V1)
# Tamaño de muestra = 50
repeticiones500_50_90= generacion_repeticiones(poblacion90E,tamaño=50,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_50_90 = mean(repeticiones500_50_90$V1)
ds_pm500_50_90 = sd(repeticiones500_50_90$V1)
# Tamaño de muestra = 60
repeticiones500_60_90= generacion_repeticiones(poblacion90E,tamaño=60,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_60_90 = mean(repeticiones500_60_90$V1)
ds_pm500_60_90 = sd(repeticiones500_60_90$V1)
# Tamaño de muestra = 100
repeticiones500_100_90= generacion_repeticiones(poblacion90E,tamaño=100,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_100_90 = mean(repeticiones500_100_90$V1)
ds_pm500_100_90 = sd(repeticiones500_100_90$V1)
# Tamaño de muestra = 200
repeticiones500_200_90= generacion_repeticiones(poblacion90E,tamaño=200,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_200_90 = mean(repeticiones500_200_90$V1)
ds_pm500_200_90 = sd(repeticiones500_200_90$V1)
# Tamaño de muestra = 500
repeticiones500_500_90= generacion_repeticiones(poblacion90E,tamaño=500,repeticiones = 500)
# observacion de propiedades en la proporcion muestral
media_pm500_500_90 = mean(repeticiones500_500_90$V1)
ds_pm500_500_90 = sd(repeticiones500_500_90$V1)
# Test de normalidad de Shapiro-Wilk
#shapiro.test(repeticiones500_1_90$V1)
#shapiro.test(repeticiones500_5_90$V1)
#shapiro.test(repeticiones500_10_90$V1)
#shapiro.test(repeticiones500_15_90$V1)
#shapiro.test(repeticiones500_20_90$V1)
#shapiro.test(repeticiones500_30_90$V1)
#shapiro.test(repeticiones500_50_90$V1)
#shapiro.test(repeticiones500_60_90$V1)
#shapiro.test(repeticiones500_100_90$V1)
#shapiro.test(repeticiones500_200_90$V1)
#shapiro.test(repeticiones500_500_90$V1)
# Realizar las pruebas de Shapiro-Wilk para los nuevos datos
results_nuevos <- data.frame(
Sample_Size = c(1, 5, 10, 15, 20, 30, 50, 60, 100, 200, 500),
p_value_nuevos = c(
shapiro.test(repeticiones500_1_90$V1)$p.value,
shapiro.test(repeticiones500_5_90$V1)$p.value,
shapiro.test(repeticiones500_10_90$V1)$p.value,
shapiro.test(repeticiones500_15_90$V1)$p.value,
shapiro.test(repeticiones500_20_90$V1)$p.value,
shapiro.test(repeticiones500_30_90$V1)$p.value,
shapiro.test(repeticiones500_50_90$V1)$p.value,
shapiro.test(repeticiones500_60_90$V1)$p.value,
shapiro.test(repeticiones500_100_90$V1)$p.value,
shapiro.test(repeticiones500_200_90$V1)$p.value,
shapiro.test(repeticiones500_500_90$V1)$p.value
)
)
# Aplicar formato a los valores p en notación científica si son muy pequeños
results_nuevos$p_value_nuevos <- format(results_nuevos$p_value_nuevos, scientific = TRUE)
# Tabla con los resultados de los nuevos datos
kable(results_nuevos, caption = "Resultados de la prueba de Shapiro-Wilk para diferentes tamaños de muestra (nuevos datos)", digits = 4, col.width = "2cm")| Sample_Size | p_value_nuevos |
|---|---|
| 1 | 5.457684e-39 |
| 5 | 2.329661e-29 |
| 10 | 1.770358e-22 |
| 15 | 1.142050e-17 |
| 20 | 7.595749e-14 |
| 30 | 1.336734e-11 |
| 50 | 1.422300e-07 |
| 60 | 1.437639e-08 |
| 100 | 1.007547e-04 |
| 200 | 6.868262e-03 |
| 500 | 2.751437e-02 |
par(mfrow=c(2, 4)) # Organizar las gráficas en una cuadrícula de 1x1
# grafico de histograma
hist(repeticiones500_1_90$V1, main = "n = 1", freq=FALSE, xlab = "Repeticiones(500) n = 1")
hist(repeticiones500_5_90$V1, main = "n = 5", freq=FALSE, xlab = "Repeticiones(500) n = 5")
hist(repeticiones500_10_90$V1, main = "n = 10", freq=FALSE, xlab = "Repeticiones(500) n = 10")
#hist(repeticiones500_15_90$V1, main = "n = 15", freq=FALSE, xlab = "Repeticiones(500) n = 15")
#hist(repeticiones500_20_90$V1, main = "n = 20", freq=FALSE, xlab = "Repeticiones(500) n = 20")
hist(repeticiones500_30_90$V1, main = "n = 30", freq=FALSE, xlab = "Repeticiones(500) n = 30")
hist(repeticiones500_50_90$V1, main = "n = 50", freq=FALSE, xlab = "Repeticiones(500) n = 50")
#hist(repeticiones500_60_90$V1, main = "n = 60", freq=FALSE, xlab = "Repeticiones(500) n = 60")
hist(repeticiones500_100_90$V1, main = "n = 100", freq=FALSE, xlab = "Repeticiones(500) n = 100")
hist(repeticiones500_200_90$V1, main = "n = 200", freq=FALSE, xlab = "Repeticiones(500) n = 200")
hist(repeticiones500_500_90$V1, main = "n = 500", freq=FALSE, xlab = "Repeticiones(500) n = 500")
Para este caso se tiene que ninguna muestra iterada hasta el tamaño n =
500 siguen una distribución normal, pero se observa que a medida que la
muestra va creciendo su media se va acercando a la del parámetro
estimado que en este caso es 0.9 y su desviación estándar se va haciendo
más pequeña (0.009).
En conclusión, de todo el proceso realizado y de los resultados obtenidos, se puede decir que a medida que la proporción muestral para estimar un parámetro va aumentando su tamaño, esta se va acercando a el valor del parámetro, a la vez que sus datos van siguiendo una distribución normal. También se puede decir que muestras iguales o mayores a un tamaño de 200, pueden tener una convergencia aproximada de los estimadores hacia una distribución normal. Por ultimo se pudo comprobar lo establecido en el teorema de tendencia central que establece que, dado un tamaño de muestra suficientemente grande, la distribución de la media de una muestra aleatoria tomada de cualquier distribución poblacional se aproximará a una distribución normal, independientemente de la forma de la distribución original.
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. Una presentación básica del método se describe a continuación:
El artículo de In-use Emissions from Heavy Duty Dissel Vehicles (J.Yanowitz, 2001) presenta las mediciones de eficiencia de combustible en millas/galón de una muestra de siete camiones. Los datos obtenidos 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∗, regresandolo nuevamente. Este procedimiento se repite hasta completar una muestra de tamaño n: X∗1,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:
Construcción de intervalos de confianza para los dos métodos
# Datos de la muestra original
muestra <- c(7.69, 4.97, 4.56, 6.49, 4.34, 6.24, 4.45)
# Tamaño de la muestra
n <- length(muestra)
# Número de muestras bootstrap
k <- 1000
# Generar muestras bootstrap
muestras_bootstrap <- matrix(NA, nrow = k, ncol = n)
for (i in 1:k) {
muestras_bootstrap[i, ] <- sample(muestra, replace = TRUE)
}
# Calcular las medias para cada muestra bootstrap
medias_bootstrap <- rowMeans(muestras_bootstrap)
# Percentiles P2.5 y P97.5
p2_5 <- quantile(medias_bootstrap, probs = 0.025)
p97_5 <- quantile(medias_bootstrap, probs = 0.975)
# Método 1: Intervalo de confianza (P2.5, P97.5)
intervalo_metodo1 <- c(p2_5, p97_5)
# Método 2: Intervalo de confianza (2X¯ - P97.5, 2X¯ - P2.5)
media_muestra <- mean(muestra)
intervalo_metodo2 <- c(2 * media_muestra - p97_5, 2 * media_muestra - p2_5)
# Resultados
cat("Intervalo de confianza (Método 1): (", round(intervalo_metodo1[1], 2), ", ", round(intervalo_metodo1[2], 2), ")\n")## Intervalo de confianza (Método 1): ( 4.71 , 6.42 )
cat("Intervalo de confianza (Método 2): (", round(intervalo_metodo2[1], 2), ", ", round(intervalo_metodo2[2], 2), ")\n")## Intervalo de confianza (Método 2): ( 4.64 , 6.36 )
Ambos métodos proporcionan intervalos que incluyen la media poblacional. Sin embargo, el método 1 (basado en cuantiles) es más estrecho, mientras que el método 2 (ajustado por la media muestral) es un poco más amplio en el rango inferior y más estrecho en el rango superior.
En general, sí confiaría en estas estimaciones, ya que ambos métodos se acercan al valor real del parámetro (media poblacional) y no difieren significativamente entre sí.
# Crear un nuevo gráfico con mejoras
hist(medias_bootstrap, las = 1, main = "Distribución de Medias Bootstrap",
xlab = "Media de Eficiencia de Combustible", ylab = "Frecuencia",
col = "#8D99AE", border = "white") # Cambio de color de las barras
# Agregar líneas para los intervalos
abline(v = intervalo_metodo1, col = "#C1292E", lwd = 2, lty = 2) # Línea punteada en rojo
abline(v = intervalo_metodo2, col = "#235789", lwd = 2, lty = 2) # Línea punteada en verde
# Leyenda
legend("topright", legend = c("Método 1", "Método 2"),
col = c("#C1292E", "#235789"), lty = 2, title = "Intervalos")
# Etiquetas para los intervalos
text(x = intervalo_metodo1[1], y = 10, labels = paste("(", round(intervalo_metodo1[1], 2), ", ", round(intervalo_metodo1[2], 2), ")", sep = ""),
pos = 1, col = "#C1292E")
text(x = intervalo_metodo2[2], y = 10, labels = paste("(", round(intervalo_metodo2[1], 2), ", ", round(intervalo_metodo2[2], 2), ")", sep = ""),
pos = 3, col = "#235789")