Definición del problema

La simulación ayuda a entender y validar 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 \(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 \(θ\) desconocido. Determine las características de cada uno de los siguientes estimadores propuestos:

\[θ_1ˆ = (X_1 + X_2)/6 + (X_3 + X_4)/3\] (1) \[θ_2ˆ = (X_1 + 2*X_2 + 3*X_3 + 4*X_4)/5\] (2) \[θ_3ˆ= (X_1 + X_2 + X_3 + X_4)/4\] (3) \[θ_4ˆ = (min(X_1, X_2, X_3, X_4) + max(X_1, X_2 ,X_3 ,X_4))/2\] (4)

Pasos

1. Cargamos la librería necesaria para hacer el ejercicio

library(ggplot2)

2. Vamos a generar las muestras para cada uno de los estimadores con diferentes tamaños de muestra \((n = 20, 50, 100, 1000)\) y suponiendo un valor para el parámetro θ.

Primero, establecemos una semilla para la reproducibilidad, que se refiere a fijar un valor inicial para el generador de números aleatorios. Al establecer una semilla, garantizamos que, si ejecutamos el mismo código nuevamente en el futuro, obtendremos exactamente la misma secuencia de números aleatorios.

set.seed(123)

Suponemos un valor para el parámetro theta \(θ\) y luego generamos muestras para cada estimador y tamaño de muestra.Creamos una lista vacía para almacenar las muestras generadas para cada tamaño \(n\) de muestras \(n = 20, 50, 100, 1000)\)

theta <- 20
n_valores <- c(20,50,100,1000)
result_propiedades <- list()

2. Ahora, vamos a Calcular los estimadores θ₁, θ₂, θ₃ y θ₄ para cada conjunto de muestras y a evaluar las propiedades de insesgadez, eficiencia y consistencia.

Iniciamos un ciclo “for” que recorre los diferentes valores de n que se encuentran en el vector “n_valores” (n_valores = 20,50,100,1000).

Para cada valor de \(n\), se generan muestras aleatorias de una distribución exponencial utilizando la función “rexp()”. El número total de valores generados es \(n * 4\), ya que tenemos cuatro variables aleatorias (X₁, X₂, X₃ y X₄) en cada muestra. El parámetro “rate” se establece en \(1/θ\) para reflejar la distribución exponencial con parámetro \(θ\). Los valores generados se almacenan en una matriz de 4 columnas utilizando la función “matrix()”, donde cada columna representa una de las variables aleatorias.

3. Luego evaluamos las propiedades de los estimadores: insesgadez, eficiencia y consistencia. Para evaluar estas propiedades se deben calcular el sesgo, la varianza, la media, las eficiencias y lasconsistencias de cada estimador y luego visualizar los resultados.

for (n in n_valores) { 
  n_muestras <- n
  muestras <- matrix(rexp(n_muestras * 4, rate = 1/theta), nrow = n_muestras)
  
  estimadores <- list(
    theta1 = function(x) (x[1] + x[2]) / 6 + (x[3] + x[4]) / 3,  #Estimador θ₁
    theta2 = function(x) (x[1] + 2*x[2] + 3*x[3] + 4*x[4]) / 5,  #Estimador θ₂
    theta3 = function(x) mean(x),                                #Estimador θ₃
    theta4 = function(x) (min(x) + max(x))/2                     #Estimador θ₄
  )
  
  resultados <- sapply(estimadores, function(estimador) {
    apply(muestras, 1, estimador)
  })

  #Calcular sesgo y varianza
  
  sesgos <- colMeans(resultados) - theta
  varianzas <- apply(resultados, 2, var)
  medias <- colMeans(resultados)
  eficiencias <- 1 /colMeans((resultados - theta)^2)
  consistencias <- colMeans((resultados - theta)^2)


  result_propiedades[[as.character(n)]]  <- data.frame(Estimador = names(estimadores), n = n, Sesgo = sesgos, Varianza = varianzas, Media = medias, Eficiencia = eficiencias, Consistencia = consistencias)

}

Los resultados de cada propiedad se almacenan en una matriz:

resultados_fin <- do.call(rbind, result_propiedades)
rownames(resultados_fin) <- NULL
print(resultados_fin)
##    Estimador    n       Sesgo  Varianza    Media  Eficiencia Consistencia
## 1     theta1   20  0.70848766 108.85332 20.70849 0.009623471    103.91261
## 2     theta2   20 20.87917566 350.75080 40.87918 0.001300131    769.15323
## 3     theta3   20  0.34257786  79.43927 20.34258 0.013230198     75.58466
## 4     theta4   20  4.01575686 208.29069 24.01576 0.004672843    214.00246
## 5     theta1   50  0.62452044  81.27212 20.62452 0.012494267     80.03671
## 6     theta2   50 21.29552387 338.92327 41.29552 0.001272841    785.64415
## 7     theta3   50  0.42985808  70.31238 20.42986 0.014473683     69.09091
## 8     theta4   50  2.74067440 101.77557 22.74067 0.009323891    107.25136
## 9     theta1  100 -0.04089629 108.65853 19.95910 0.009295960    107.57362
## 10    theta2  100 20.43759565 440.94993 40.43760 0.001170637    854.23575
## 11    theta3  100 -0.07060324  96.07281 19.92940 0.010513360     95.11707
## 12    theta4  100  3.37645298 141.99101 23.37645 0.006580180    151.97154
## 13    theta1 1000  0.07867237 117.64557 20.07867 0.008508168    117.53411
## 14    theta2 1000 20.11631579 508.84010 40.11632 0.001095293    912.99743
## 15    theta3 1000  0.10106150 110.03370 20.10106 0.009096377    109.93388
## 16    theta4 1000  3.49834936 169.15804 23.49835 0.005517932    181.22733

Ahora vamos a ver los gráficos que nos dan una idea mucho más clara de los resultados obtenidos

Primero definimos las variables que incluirán cada gráfica

sesgo_bxp <- ggplot(resultados_fin, aes(Estimador, y = Sesgo)) + geom_boxplot() + theme_minimal()

varianza_bxp <- ggplot(resultados_fin, aes(Estimador, y = Varianza)) + geom_boxplot() + theme_minimal()

eficiencia_bxp <- ggplot(resultados_fin, aes(Estimador, y = Eficiencia)) + geom_boxplot() + theme_minimal()

consistencia_bxp <- ggplot(resultados_fin, aes(Estimador, y = Consistencia)) + geom_boxplot() + theme_minimal()

Y procedemos con los boxplots de las 4 propiedades:

sesgo_bxp

varianza_bxp

eficiencia_bxp

consistencia_bxp