Introducción

El proceso de simulación es una herramienta poderosa en el campo de la estadística y la cienciad de datos. Se puede emplear para entender relaciones complejas y estimar valores que son difíciles de calcular directamente. La simulación permite a los investigadores crear modelos que representen situaciones del mundo real y experimentar con ellos para obtener información valiosa. Esto puede ser especialmente útil en situaciones donde los datos son escasos o difíciles de obtener, o cuando los modelos matemáticos son demasiado complejos para resolver analíticamente.

Para entender cómo funciona la simulación y cómo se puede aplicar en la práctica, se plantean los siguientes problemas. A través de estos, se ilustrará cómo la simulación puede ser utilizada para obtener estimaciones precisas y confiables, incluso en situaciones complejas.

Librerías usadas:

library(tidyverse)
library(skimr)
library(stringi)
library(highcharter)
library(readxl)
library(heatmaply)
library(plotly)

Problema 1 - Estimación

La siguiente figura sugiere como estimar el valor de \(\pi\) con una simulación. En la figura, un circuito con un área igual a \(\pi\)/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 \(\pi\)/4. Por tanto, se puede estimar el valor de \(\pi\)/4 al contar el número de puntos dentro del círculo, para obtener la estimación de \(\pi\)/4. De este último resultado se encontrar una aproximación para el valor de \(\pi\).

Pasos sugeridos:

\(a.\) Genere \(n\) coordenadas \(x\): \(X_1, \ldots, X_n\). 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)\).

\(b.\) Genere 1000 coordenadas \(y\): \(Y_1,\ldots,Y_n\), utilizando nuevamente la distribución uniforme con valor mínimo de 0 y valor máximo de 1.

\(c.\) Cada punto \((X_i,Y_i)\) se encuentra dentro del círculo si su distancia desde el centro \((0.5,0.5)\) es menor a 0.5. Para cada par \((X_i,Y_i)\) determine si la distancia desde el centro es menor a 0.5. Esto último se puede realizar al calcular el valor \((X_i-0.5)^2+(Y_i-0.5)^2\), que es el cuadrado de la distancia, y al determinar si es menor que 0.25.

Preguntas:

La simulación debe resolver las siguientes preguntas: ¿Cuántos de los puntos están dentro del círculo? ¿Cuál es su estimación de \(\pi\)?

Solución:

Para iniciar se generan \(1000\) valores de \(X\) y \(Y\) con el fin de construir las coordenadas para los puntos \((X_i,Y_i)\) de la siguiente manera:

# Número de puntos
n1 <- 1000

# Generar coordenadas x e y
x <- runif(n1, min = 0, max = 1)
y <- runif(n1, min = 0, max = 1)

tabla <- data.frame(x,y)
head(tabla)
##            x          y
## 1 0.81706316 0.18398733
## 2 0.03993966 0.17478236
## 3 0.51279651 0.68760398
## 4 0.32857806 0.43509626
## 5 0.86529309 0.08041608
## 6 0.09726035 0.69366773

Ahora se calculan los puntos que están dentro del circulo y se grafican

# Calcular la distancia al centro del círculo
distancia <- (x - 0.5)^2 + (y - 0.5)^2

# Determinar si los puntos están dentro del círculo
dentro_circulo <- distancia < 0.25

# Crear un conjunto de datos para los puntos
puntos_data <- data.frame(
  x = x,
  y = y,
  color = ifelse(dentro_circulo, "Dentro", "Fuera")
)

# Crear un gráfico vacío
p <- ggplot() + xlim(0, 1) + ylim(0, 1)

# Agregar un cuadrado de lado 2 centrado en el origen
p <- p + annotate("rect", xmin = 0, xmax = 1, ymin = 0, ymax = 1, fill = NA, color = "black")

# Crear un conjunto de datos para el círculo
circle_data <- data.frame(
  x = 0.5 + 0.5 * cos(seq(0, 2*pi, length.out = 100)),
  y = 0.5 + 0.5 * sin(seq(0, 2*pi, length.out = 100))
)

# Agregar el círculo al gráfico
p <- p + geom_path(data = circle_data, aes(x, y))

# Establecer una relación de aspecto fija entre los ejes x e y
p <- p + coord_fixed()

# Agregar los puntos al gráfico
p <- p + geom_point(data = puntos_data, aes(x, y, color = color))

# Aplicar un tema minimalista al gráfico
p <- p + theme_minimal()

# Mostrar el gráfico
p

Finalmente se resuelven las preguntas planteadas al inicio:

# Contar cuántos puntos están dentro del círculo
puntos_dentro_circulo <- sum(dentro_circulo)

# Estimar el valor de pi
pi_estimado <- 4 * puntos_dentro_circulo / n1 # Lo mismo que calcular la media mean(dentro_circulo)

pi_symbol <- intToUtf8(0x03C0)

# Mostrar el resultado
cat("Número de puntos dentro del círculo:", puntos_dentro_circulo, "\n")
## Número de puntos dentro del círculo: 786
cat("Estimación de", pi_symbol, "es", pi_estimado)
## Estimación de π es 3.144

Problema 2 - Propiedades de los estimadores

La simulación ayuda a entender y validar las propiedades de los estimadores estadísticos, como son la insesgadez, eficiencia y 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 \(X_1\), \(X_2\), \(X_3\) y \(X_4\), una muestra aleatoria de tamaño \(n=4\) cuya población la conforma una distribución exponencial con parámetro \(\theta\) desconocido. Determine las características de cada uno de los siguientes estimadores propuestos:

Los estimadores propuestos son los siguientes:

Solución:

Se inicia construyendo la función llamada estimadores que permite evaluar los diferentes estimadores planteados para \(4\) valores obtenidos de forma aleatoria a partir de la función de distribución exponencial. Para ello, se hará uso de las funciones rexp(), apply() y function(){}. El código es el siguiente:

### Contrucción de función estimadores

estimadores <- function(repeticiones, theta) {
  X_rep <- lapply(1:repeticiones, function(x) rexp(4, rate=1/theta))
  X_rep <- do.call(rbind, X_rep)
  
  resultados <- t(apply(X_rep, 1, function(x) {
    theta1 <- ((x[1]+x[2])/6) + ((x[3]+x[4])/3)
    theta2 <- (x[1]+(2*x[2])+(3*x[3])+(4*x[4]))/5
    theta3 <- (x[1]+x[2]+x[3]+x[4])/4
    theta4 <- (pmin(x[1],x[2],x[3],x[4])+pmax(x[1],x[2],x[3],x[4]))/2
    c(theta1, theta2, theta3, theta4)
  }))
  
  colnames(resultados) <- c("theta1", "theta2", "theta3", "theta4")
  
  resultados <- as.data.frame(resultados)
  
  return(resultados)
}

Una vez construida la función estimadores, se procede a definir un valor para el parámetro \(\theta\), en este caso será \(1.5\). Luego, se evalúan los estimadores un número diferente de veces, específicamente \(20, 50, 100\) y \(1000\) veces.

### Se define el parámetro θ
theta <- 1.5

### Se evalúa cada uno de los estimadores 20, 50, 100 y 1000 veces 
resultados_20 <- estimadores(20, theta)
resultados_50 <- estimadores(50, theta)
resultados_100 <- estimadores(100, theta)
resultados_1000 <- estimadores(1000, theta)

A continuación, se presentan los resultados obtenidos de la evaluación de los estimadores realizada \(20\) veces. Es importante entender que los resultados de las otras evaluaciones son similares, solo que con un mayor número de observaciones.

resultados_20
##       theta1   theta2    theta3    theta4
## 1  2.6152266 5.185129 2.5844406 2.4920826
## 2  1.7964890 3.839316 1.4440718 1.7074989
## 3  1.6360781 3.151090 1.8766000 1.9387435
## 4  2.5220090 5.039641 2.3443660 2.5329513
## 5  1.0834674 2.055571 1.1936409 1.0254061
## 6  3.1981015 5.897198 2.8231068 2.7804245
## 7  1.7537730 2.861925 1.8805267 2.0545789
## 8  1.3284419 2.475762 1.5234416 1.6993053
## 9  2.1694157 3.516212 2.4567340 2.8798839
## 10 2.2031422 4.698111 2.0640283 2.4813698
## 11 0.5437857 1.189261 0.5591620 0.5130331
## 12 0.7809074 1.595772 0.9997167 1.1562391
## 13 0.6254788 1.123329 0.7306096 0.9445503
## 14 1.9941694 4.317138 2.1601825 2.6582215
## 15 1.0173734 1.830577 1.0768904 1.2554413
## 16 0.5509476 1.254415 0.7122764 1.1962629
## 17 2.8522718 5.183612 3.5622663 4.1546001
## 18 0.9656547 1.797300 1.1213892 1.3006792
## 19 1.9324625 4.163996 1.6286978 2.0726105
## 20 1.3134255 2.653633 1.3146172 1.3110422

Ahora, se calculan las estimaciones de los \(4\) estimadores utilizando la media y la varianza para proceder con su respectiva comparación. También, se realiza el cálculo del sesgo.

resultados <- list(resultados_20, resultados_50, resultados_100, resultados_1000)
nombres <- c("20", "50", "100", "1000")

estadisticas <- lapply(resultados, function(x) {
  media <- colMeans(x)
  varianza <- apply(x, 2, var)
  sesgo <- media - theta
  c(media, varianza, sesgo)
})

estadisticas_df <- do.call(rbind, estadisticas)
rownames(estadisticas_df) <- nombres
colnames(estadisticas_df) <- c(paste0("media_theta", 1:4), paste0("varianza_theta", 1:4), paste0("sesgo_theta", 1:4))

Media <- estadisticas_df[,1:4]
Varianza <- estadisticas_df[,5:8]
Sesgo <- estadisticas_df[,9:12]

Con lo anterior, se tienen los resultados para iniciar con la evaluación y comparación de los \(4\) estimadores a partir de las propiedades de Insesgadez, Eficiencia y Consistencia. Se recuerda que el parámetro usado es \(\theta=1.5\)

Insesgadez

Para esta propiedad, se compararán los resultados del sesgo de cada conjunto de datos en términos de su valor absoluto. De esta manera, la comparación se llevará a cabo determinando cuáles son los estimadores que presentan el menor sesgo

# Se muestran los valores obtenidos para cada estimador
sesgos_absolutos <- abs(Sesgo)
sesgos_absolutos
##      sesgo_theta1 sesgo_theta2 sesgo_theta3 sesgo_theta4
## 20     0.14413107     1.691449  0.202838240    0.4077462
## 50     0.04261678     1.542668  0.009350742    0.2937024
## 100    0.06424617     1.624766  0.063561209    0.3025096
## 1000   0.00463893     1.483739  0.009489655    0.2390503
# Para determinar el estimador más insesgado
min_sesgos_pos <- apply(sesgos_absolutos, 1, which.min) # Encuentra la posición del primer sesgo mínimo en cada fila

# Imprime el estimador más insesgado para cada conjunto de resultados
for (i in seq_along(nombres)) {
  cat(paste("Para", nombres[i], "repeticiones, el estimador más insesgado es el estimador", min_sesgos_pos[i], "\n"))
}
## Para 20 repeticiones, el estimador más insesgado es el estimador 1 
## Para 50 repeticiones, el estimador más insesgado es el estimador 3 
## Para 100 repeticiones, el estimador más insesgado es el estimador 3 
## Para 1000 repeticiones, el estimador más insesgado es el estimador 1

Se puede observar que los estimadores \(\hat{\theta}_1\) y \(\hat{\theta}_3\) parecen ser insesgados. Aunque la teoría establece que el estimador \(\hat{\theta}_3\) es, en efecto, un estimador insesgado, resulta interesante que el estimador \(\hat{\theta}_1\) también muestre esta propiedad en ciertos casos. Esto podría ser el resultado de la aleatoriedad.

Eficiencia

Para evaluar esta propiedad, se analizarán los resultados de la varianza para cada conjunto de datos. La comparación se realizará identificando los estimadores que exhiben la varianza más baja. Esto permitirá determinar cuál de los estimadores es el más eficiente.

# Se muestran los valores obtenidos para cada estimador
Varianza
##      varianza_theta1 varianza_theta2 varianza_theta3 varianza_theta4
## 20         0.6306655        2.330016       0.6233868       0.7571243
## 50         0.6554091        2.694804       0.5677870       1.0857602
## 100        0.7365934        3.481193       0.6342152       1.0064951
## 1000       0.6561275        2.821983       0.5890559       0.9726123
# Para determinar el estimador más eficiente (el que tiene la varianza más baja)
min_var_pos <- apply(Varianza, 1, which.min) # Encuentra la posición de la primera varianza mínima en cada fila

# Imprime el estimador más eficiente para cada conjunto de resultados
for (i in seq_along(nombres)) {
  cat(paste("Para", nombres[i], "repeticiones, el estimador más eficiente es el estimador", min_var_pos[i], "\n"))
}
## Para 20 repeticiones, el estimador más eficiente es el estimador 3 
## Para 50 repeticiones, el estimador más eficiente es el estimador 3 
## Para 100 repeticiones, el estimador más eficiente es el estimador 3 
## Para 1000 repeticiones, el estimador más eficiente es el estimador 3

Se observa que el estimador \(\hat{\theta}_3\), es en general, el más eficiente.

Consistencia

Para evaluar esta propiedad, es necesario determinar si el incremento en la cantidad de datos contribuye a una estimación más precisa del parámetro. Esto se puede hacer observando el número de repeticiones utilizadas para la estimación de los diferentes estimadores y determinando si un aumento en el tamaño de la muestra mejora la estimación hacia el valor del parámetro, en este caso conocido, \(\theta=1.5\). Además, se puede realizar una comparación entre los estimadores para determinar cuáles cumplen con esta propiedad.

# Se muestran los valores obtenidos para cada estimador
Media
##      media_theta1 media_theta2 media_theta3 media_theta4
## 20       1.644131     3.191449     1.702838     1.907746
## 50       1.542617     3.042668     1.509351     1.793702
## 100      1.564246     3.124766     1.563561     1.802510
## 1000     1.495361     2.983739     1.490510     1.739050

Se observa que el estimador de \(\hat{\theta}_1\) y \(\hat{\theta}_3\) muestran la propiedad de consistencia en general.

Conclusiones

Los estimadores que mejores resultados mostraron en general son:

\(\hat{\theta}_1 = \frac{X_1 + X_2}{6} + \frac{X_3 + X_4}{3}\)

\(\hat{\theta}_3 = \frac{X_1 + X_2 + X_3 + X_4}{4}\)

Sin embargo, el estimador \(\hat{\theta}_3\) es el que, muestra buenos resultados en las tres propiedades en general. Por lo tanto, se concluye que este es el mejor estimador entre los cuatro evaluados.

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 se refiere a la convergencia de los estimadores, como la proporción muestral, a la distribución normal. Algunos autores sostienen que esta aproximación es bastante buena a partir del umbral \(n>30\).

A continuación, se describen los pasos para su verificación:

\(a.\) Realizar una simulación en la que se genere una población de \(n=1000\) que corresponda a una población de plantas sembradas en un cultivo, donde el porcentaje de plantas enfermas sea del \(50\%\).

\(b.\) Generar una función que permita:

\(c.\) Repetir el paso \(b.\) \(500\) veces y analizar los resultados en cuanto al comportamiento de los \(500\) resultados del estimador \(\hat{p}\)

\(d.\) Repetir los puntos \(b.\) y \(c.\) para tamaños de muestra \(n=5, 10, 15, 20, 30, 50, 60, 100, 200, 500\). Comparar los resultados obtenidos para los diferentes tamaños de muestra en cuanto a la normalidad. Utilizar pruebas de bondad y ajuste (Shapiro-Wilks: shapiro.test()) y métodos gráficos (gráfico de normalidad: qqnorm()).

Solución:

Para realizar la simulación se inicia construyendo la población de plantas enfermas:

# Se define el tamaño de la población y la proporción de plantas enfermas
n <- 1000
p <- 0.5

# Se genera la población
enfermas <- rep(c(0, 1), each = n*p) #Esta población son 500 valores de 0 y 500 valores de 1

Ahora se construye una función llamada muestra_estimacion_p que permite obtener una muestra aleatoria de la población anterior y calcular el estimador de la proporción muestral \(\hat{p}\) para un tamaño de muestra dado \(m\).

muestra_estimacion_p <- function(poblacion, m) {
  muestra <- sample(poblacion, m, replace = TRUE)
  p_ <- mean(muestra)
  return(list(muestra_p=muestra, p_estimado = p_))
}

muestra_estimacion_p(enfermas,10) # Se realiza un ejemplo del funcionamiento de la función para 10 valores obtenidos de la población de forma aleatoria y con repetición.
## $muestra_p
##  [1] 0 1 1 1 0 0 1 1 0 1
## 
## $p_estimado
## [1] 0.6

Una vez que se ha verificado el correcto funcionamiento de la función anterior, se procede a repetir su uso 500 veces. Esto se hace con el objetivo de analizar los resultados en términos del comportamiento del estimador \(\hat{p}\).

# Número de repeticiones
n_repetitions <- 500

# Crear data frame vacío para almacenar resultados
results <- data.frame(repeticion = numeric(n_repetitions), estimacion = numeric(n_repetitions))


# Obtener estimaciones y guardar en data frame
for (i in 1:n_repetitions) {
  results$repeticion[i] <- i
  results$estimacion[i] <- muestra_estimacion_p(enfermas,55)$p_estimado
}

Se responden las preguntas formuladas en el apartado \(c.\)

\(1.\) ¿Qué tan simétricos o sesgados son los resultados obtenidos?

# Calculando el sesgo
Sesgo <- mean(results$estimacion) - p
abs(Sesgo)
## [1] 0

En relación al sesgo, los valores que se obtienen son muy bajos, por lo que se podría asumir que la estimación para \(\hat{p}\) es insesgada.

# Histograma
ggplot(results, aes(x=estimacion)) +
  geom_histogram(aes(y=after_stat(count)), binwidth=0.1, colour="black", fill="white") +
  geom_vline(aes(xintercept=0.5), color="red", linetype="dashed") + # Añade una línea vertical en 0.5
  annotate("text", x = 0.5, y = Inf, label = "Parámetro p", vjust = 2, hjust=-1, color = "red") + # Añade una etiqueta para la línea
  labs(title=bquote("Histograma de los valores de" ~ hat(p)), x=expression(hat(p)), y="Frecuencia") +
  theme(plot.title = element_text(hjust = 0.5)) # Centra el título

En cuanto a la simetría, el histograma muestra que las 500 estimaciones de \(\hat{p}\) se distribuyen de manera simétrica.

\(2.\) ¿Qué se puede observar en cuanto a variabilidad?

# Calcula el coeficiente de variación
cv <- (sd(results$estimacion) / mean(results$estimacion)) * 100
cv
## [1] 13.18968

El resultado del coeficiente de variación es aceptable, pero probablemente se deba a que se están calculando los \(\hat{p}\) con muestras de tamaño 55, teniendo en cuenta que \(n>30\). Por lo tanto, se espera que al aumentar el tamaño de la muestra, los resultados sean cada vez menos dispersos, como menciona la teoría.

Finalmente, teniendo en cuenta los tamaños de muestra \(n=5, 10, 15, 20, 30, 50, 60, 100, 200, 500\), se utiliza la función creada anteriormente, muestra_estimacion_p() que sirve para extraer las muestras de la población y estimar \(\hat{p}\), para crear una matriz que servirá de insumo para almacenar los valores simulados y finalmente comparar los resultados obtenidos respecto a su simetría y normalidad.

El siguiente código sirve de insumo para construir luego las gráficas y comparaciones:

# Genera una matriz donde cada columna es una muestra extraída de la población
m <- 1000 # Número de muestras
X <- matrix(unlist(lapply(1:m, function(i) muestra_estimacion_p(enfermas, 100)$muestra_p)), ncol=m)

# Generación de muestras
X5=X[ ,1:5]          # n=5
X10=X[ ,1:10]        # n=10
X15=X[ ,1:15]        # n=15
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

# Cálculo de medias
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
Mx15=apply(X15,1,mean)      # medias de muestras de tamaño n=15
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

La prueba de Shapiro-Wilk se utiliza para determinar si una muestra de datos proviene de una distribución normal. Las hipótesis para esta prueba son las siguientes:

Hipótesis nula \((H_0)\): Los datos siguen una distribución normal.

Hipótesis alternativa \((H_1)\): Los datos no siguen una distribución normal.

Para verificar si las estimaciones de \(\hat{p}\) siguen una distribución normal, se examinará si los valores p rechazan o no la hipótesis nula \(H_0\).

# Crear una lista con todas las medias
medias_list <- list(Mx5, Mx10, Mx15, Mx20, Mx30, Mx50, Mx60, Mx100, Mx200, Mx500)
names(medias_list) <- c("n=5", "n=10", "n=15", "n=20", "n=30", "n=50", "n=60", "n=100", "n=200", "n=500")

# Aplicar la prueba de Shapiro-Wilk a cada conjunto de medias
resultados_shapiro <- lapply(medias_list, shapiro.test)

# Extraer los valores p de los resultados
valores_p <- sapply(resultados_shapiro, function(resultado) resultado$p.value)

# Desactivar la notación científica
options(scipen=999)

# valores p
valores_p
##           n=5          n=10          n=15          n=20          n=30 
## 0.00006982919 0.00905469756 0.00636954312 0.07564326846 0.01165964001 
##          n=50          n=60         n=100         n=200         n=500 
## 0.05326189483 0.60112750550 0.56465651645 0.22266541979 0.19902039498
# Rechazo si valor-p > 0.05

valores_p < 0.05 # Se rechaza la hipotesis nula de la prueba Shapiro wilks
##   n=5  n=10  n=15  n=20  n=30  n=50  n=60 n=100 n=200 n=500 
##  TRUE  TRUE  TRUE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE

Los resultados muestran que con valores de muestra más grandes se empieza a no rechazar \(H_0\).

Ahora la idea es mostrar mediante un histograma el comportamiento de las estimaciones de \(\hat{p}\) en la medida de que aumentan los tamaños muestrales.

par(mfrow=c(2,5), mar = c(4.1, 3.1, 3.1, 1.1)) # divide la gráfica en una matriz 

hist(X5,  main = "n=5", freq=FALSE, xlab="", ylab="")
hist(Mx10, main = "n=10",freq=FALSE, xlab="", ylab="")
hist(Mx15, main = "n=15",freq=FALSE, xlab="", ylab="")
hist(Mx20, main = "n=20",freq=FALSE, xlab="", ylab="")
hist(Mx30, main = "n=30",freq=FALSE,xlab="", ylab="") 
hist(Mx50, main = "n=50",freq=FALSE,xlab="", ylab="") 
hist(Mx60, main = "n=60", freq = FALSE,xlab="", ylab="")
hist(Mx100, main = "n=100", freq = FALSE,xlab=expression(hat(p)), ylab="") 
hist(Mx200, main = "n=200",freq = FALSE,xlab="", ylab="")
hist(Mx500, main = "n=500", freq = FALSE,xlab="", ylab="")

Se puede observar fácilmente cómo, a medida que aumentan los tamaños de las muestras en las que se calcula \(\hat{p}\), la distribución adquiere cada vez más una forma simétrica, lo cual es característico de los datos que siguen una distribución normal.

Finalmente, se realiza una prueba gráfica conocida como qqplot(). La idea es que los puntos generados por los cuantiles por cada muestra, a medida que aumentan, se ajusten cada vez mejor a la curva qqline. Esto demuestra que existe un ajuste hacia la distribución normal.

par(mfrow=c(2,5), mar = c(4.1, 3.1, 3.1, 1.1))  # divide la gráfica en una matriz

qqnorm(X5, main="n=5", xlab="", ylab="") ; qqline(X5, col="red")
qqnorm(Mx10, main ="n=10", xlab="", ylab="") ; qqline(Mx10, col="red")
qqnorm(Mx15, main ="n=15", xlab="", ylab="") ; qqline(Mx15, col="red")
qqnorm(Mx20, main ="n=20", xlab="", ylab="") ; qqline(Mx20, col="red")
qqnorm(Mx30, main = "n=30", xlab="", ylab="") ; qqline(Mx30, col="red")
qqnorm(Mx50, main = "n=50", xlab="", ylab="") ; qqline(Mx50, col="red")
qqnorm(Mx60, main ="n=60", xlab="", ylab="") ; qqline(Mx60, col="red")
qqnorm(Mx100, main ="n=100",xlab="Cuantiles Teorícos" ,ylab="") ; qqline(Mx100, col="red")
qqnorm(Mx200, main ="n=200", xlab="", ylab="") ; qqline(Mx200, col="red")
qqnorm(Mx500, main="n=500", xlab="", ylab="") ; qqline(Mx500, col="red")

Se puede observar que a medida que aumentan los tamaños de las muestras, el ajuste hacia la qqline mejora. Esto demuestra una convergencia entre los cuantiles teóricos normales y los cuantiles muestrales, lo que sugiere que estas estimaciones tienden a seguir una distribución normal.

Problema 4 - Estimación Boostrap

Cuando se extrae una muestra de una población que no sigue una distribución normal y se necesita estimar un intervalo de confianza, los métodos de estimación bootstrap pueden ser útiles. Esta metodología asume que es posible reconstruir la población objeto de estudio mediante un muestreo con reemplazo de la muestra disponible. Existen varias versiones del método. A continuación, se presenta una descripción básica del método:

El artículo “In-use Emissions from Heavy Duty Diesel Vehicles” (J.Yanowitz, 2001) presenta las mediciones de eficiencia de combustible en millas por 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, 4.45\). Se supone que es una muestra aleatoria de camiones y 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 sobre la distribución de los datos.

El método bootstrap permite construir intervalos de confianza del \(95\%\). Para ilustrar el método, supongamos que colocamos los valores de la muestra en una caja y extraemos uno al azar. Este correspondería al primer valor de la muestra bootstrap \(X^*_1\). Después de anotar el valor, se regresa \(X^*_1\) a la caja y se extrae el valor \(X^*_2\), devolviéndolo 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 (supongamos \(k = 1000\)). Para cada una de las muestras bootstrap obtenidas, se calcula la media \(\bar{X}^*_i\), obteniéndose un valor para cada muestra. El intervalo de confianza queda conformado por los percentiles \(P_{2.5}\) y \(P_{97.5}\).

Los dos métodos para estimar el intervalo de confianza son:

Método 1: \((P_{2.5}; P_{97.5})\)

Método 2: \((2\bar{X} - P_{97.5}; 2\bar{X} - P_{2.5})\)

Solución:

Para la construcción del bootstrap, se extraen \(1000\) muestras del vector llamado datos. A cada muestra extraída se le calcula su media y se guarda en otro vector llamado medias.

# Datos de la muestra
datos <- c(7.69, 4.97, 4.56, 6.49, 4.34, 6.24, 4.45)

# Número de muestras bootstrap
k <- 10000

# Tamaño de la muestra
n <- length(datos)

# Vector para almacenar las medias de las muestras bootstrap
medias <- numeric(k)

# Generar k muestras bootstrap y calcular la media de cada una
for (i in 1:k) {
  muestra_bootstrap <- sample(datos, n, replace = TRUE)
  medias[i] <- mean(muestra_bootstrap)
}

Ahora, se procede a construir los intervalos de confianza al \(95\%\) utilizando ambos métodos. Para el primer método, utilizaremos la función quantile. Para el segundo método, se puede hacer uso del primer método, ya que consiste en usar los valores anteriores y ajustarlos o corregirlos, usando dos veces el valor de la media del vector datos.

# Calcular los percentiles P2.5 y P97.5
percentiles <- quantile(medias, c(0.025, 0.975))

# Intervalo de confianza por el método 1
IC1 <- percentiles

# Intervalo de confianza por el método 2
media_muestra <- mean(datos)
IC2 <- c(2 * media_muestra - percentiles[2], 2 * media_muestra - percentiles[1])

# Mostrar los resultados
cat("Intervalo de confianza por el método 1: (", IC1[1], ";", IC1[2], ")\n")
## Intervalo de confianza por el método 1: ( 4.721429 ; 6.485714 )
cat("Intervalo de confianza por el método 2: (", IC2[1], ";", IC2[2], ")\n")
## Intervalo de confianza por el método 2: ( 4.582857 ; 6.347143 )

Los resultados son muy similares entre sí. Los intervalos con un nivel de confianza del \(95\%\), al parecer poseen la misma diferencia entre sus intervalos superior e inferior.

A continuación, se construye un histograma para visualizar este comportamiento en relación con la media del vector datos.

# Graficar el histograma de las medias bootstrap con ggplot
ggplot(data.frame(medias = medias), aes(x = medias)) +
  geom_histogram(bins = 20, color = "black", fill = "lightgray") +
  labs(title = "Distribución de las medias bootstrap", x = "Media", y = "Frecuencia") +
  geom_vline(xintercept = IC1[1], linetype = 2, color = "red") +
  geom_vline(xintercept = IC1[2], linetype = 2, color = "red") +
  geom_vline(xintercept = IC2[1], linetype = 2, color = "blue") +
  geom_vline(xintercept = IC2[2], linetype = 2, color = "blue") +
  geom_vline(xintercept = mean(datos), linetype = 1, color = "green") + theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

La media de los datos iniciales es de \(5.53\), por lo que está dentro de ambos intervalos de confianza. Esto significa que, con un \(95\%\) de confianza, se puede concluir que la media poblacional es igual o similar a la media de la muestra. Al comparar ambos métodos, se puede entender el segundo método como un ajuste o corrección del primer método. Ambos métodos son confiables para estimar un intervalo de confianza cuando se desconoce la distribución de probabilidad de la que provienen los datos.