Christian David Vera Mendivelso
Paula
Vidal Godoy
Métodos y Simulación estadística
Maestría en Ciencia de Datos
Pontificia Universidad Javeriana de
Cali
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 estimadores propuestos.
Valor esperado y varianza de la función exponencial.
\({E}[X] = \theta\\\) \(\text{Var}(X) = \theta^2\)
Para comprender los conceptos de insesgadez, eficiencia y
consistencia, se calculará el valor esperado y la varianza de
los cuatro estimadores cuando n = 4, los cuales se usaran para
determinar las propiedades de insesgadez y eficiencia junto con
simulaciones de cada uno.
Valor esperado \(\hat{\theta}_1\)
\(\mathbb{E}[\hat{\theta}_1] = \frac{\mathbb{E}[X_1] + \mathbb{E}[X_2]}{6} + \frac{\mathbb{E}[X_3] + \mathbb{E}[X_4]}{3}\)
\(\mathbb{E}[\hat{\theta}_1] = \frac{\theta + \theta}{6} + \frac{\theta + \theta}{3} = \frac{2\theta}{6} + \frac{2\theta}{3} = \frac{\theta}{3} + \frac{2\theta}{3} = \theta\)
Varianza \(\hat{\theta}_1\)
\(\text{Var}(\hat{\theta}_1) = \text{Var}\left(\frac{X_1 + X_2}{6} + \frac{X_3 + X_4}{3}\right)\)
\(\text{Var}(\hat{\theta}_1) = \text{Var}\left(\frac{X_1 + X_2}{6}\right) + \text{Var}\left(\frac{X_3 + X_4}{3}\right)\)
\(\text{Var}(\hat{\theta}_1) = \frac{1}{36} \text{Var}(X_1 + X_2) + \frac{1}{9} \text{Var}(X_3 + X_4)\)
\(\text{Var}(\hat{\theta}_1) = \frac{1}{36} \cdot 2\theta^2 + \frac{1}{9} \cdot 2\theta^2 = \frac{2\theta^2}{36} + \frac{2\theta^2}{9} = \frac{\theta^2}{18} + \frac{2\theta^2}{9} = \frac{5\theta^2}{18}\)
Valor esperado \(\hat{\theta}_2\)
\(\mathbb{E}[\hat{\theta}_2] = \frac{\mathbb{E}[X_1] + 2\mathbb{E}[X_2] + 3\mathbb{E}[X_3] + 4\mathbb{E}[X_4]}{5}\)
\(\mathbb{E}[\hat{\theta}_2] = \frac{\theta + 2\theta + 3\theta + 4\theta}{5} = \frac{10\theta}{5} = 2\theta\)
Varianza \(\hat{\theta}_2\)
\(\text{Var}(\hat{\theta}_2) = \text{Var}\left(\frac{X_1 + 2X_2 + 3X_3 + 4X_4}{5}\right)\)
\(\text{Var}(\hat{\theta}_2) = \frac{1}{25} \text{Var}(X_1 + 2X_2 + 3X_3 + 4X_4)\)
\(\text{Var}(X_1 + 2X_2 + 3X_3 + 4X_4) = \theta^2 + 4\theta^2 + 9\theta^2 + 16\theta^2 = 30\theta^2\)
\(\text{Var}(\hat{\theta}_2) = \frac{30\theta^2}{25} = \frac{6\theta^2}{5}\)
Valor esperado \(\hat{\theta}_3\)
\(\mathbb{E}[\hat{\theta}_3] = \frac{\mathbb{E}[X_1] + \mathbb{E}[X_2] + \mathbb{E}[X_3] + \mathbb{E}[X_4]}{4}\)
\(\mathbb{E}[\hat{\theta}_3] = \frac{\theta + \theta + \theta + \theta}{4} = \frac{4\theta}{4} = \theta\)
Varianza \(\hat{\theta}_3\)
\(\text{Var}(\hat{\theta}_3) = \text{Var}\left(\frac{X_1 + X_2 + X_3 + X_4}{4}\right)\)
\(\text{Var}(\hat{\theta}_3) = \frac{1}{16} \text{Var}(X_1 + X_2 + X_3 + X_4)\)
\(\text{Var}(\hat{\theta}_3) = \frac{4\theta^2}{16} = \frac{\theta^2}{4}\)
Valor esperado y varianza \(\hat{\theta}_4\)
Debido a la dificultad para calcular estos valores de manera analítica, se usaran las simulaciones para estimar estos valores.
Ahora se determinará las propiedades de insesgadez, eficiencia y consistencia para cada uno de los 4 estimadores fijando como parámetro \({\theta}=5\) y generando la semilla 123 para evitar que los resultados varíen.
#Se calcula la funciones de los estimadores y la función con la cual se hacen las simulaciones
# Parámetro y semilla
theta <- 5
set.seed(123)
# Funciones para los estimadores
estimador1 <- function(X) {
return((X[1] + X[2]) / 6 + (X[3] + X[4]) / 3)}
estimador2 <- function(X) {
return((X[1] + 2*X[2] + 3*X[3] + 4*X[4]) / 5)}
estimador3 <- function(X) {
return((X[1] + X[2] + X[3] + X[4]) / 4)}
estimador4 <- function(X) {
return((min(X) + max(X)) / 2)}
# Función para simular datos y calcular estimadores
#cada muestra tiene un tamaño de 4
simular_estimadores <- function(n, theta) {
muestras <- matrix(rexp(4*n, rate = 1 / theta), ncol = 4)
estimadores <- data.frame(
Estimador1 = apply(muestras, 1, estimador1),
Estimador2 = apply(muestras, 1, estimador2),
Estimador3 = apply(muestras, 1, estimador3),
Estimador4 = apply(muestras, 1, estimador4)
)
return(estimadores)
}
# Simular datos para diferentes tamaños
muestras_20 <- simular_estimadores(20, theta)
muestras_50 <- simular_estimadores(50, theta)
muestras_100 <- simular_estimadores(100, theta)
muestras_1000 <- simular_estimadores(1000, theta)Un estimador se considera insesgado si \(E[\hat{\theta}] = \theta\), la Tabla 1 muestra el promedio para cada estimador con 20, 50, 100 y 1.000 simulaciones de una muestra aleatoria de tamaño 4, \(X_1\), \(X_2\), \(X_3\) y \(X_4\). A continuación se menciona la insesgadez de cada estimador con base en la Tabla 1 y el valor esperado calculado para cada uno.
\(\hat{\theta}_1\) : Dado que su valor esperado es igual al parámetro \(\theta\), \(\hat{\theta}_1\) es considerado insesgado.
\(\hat{\theta}_2\) : Dado que su valor esperado difiere al parámetro \(\theta\), \(\hat{\theta}_2\) es considerado sesgado.
\(\hat{\theta}_3\) : Dado que su valor esperado es igual al parámetro \(\theta\), \(\hat{\theta}_3\) es considerado insesgado.
\(\hat{\theta}_4\) : En la Tabla 1 se evidencia que la media de del estimador no se aproxima al valor de \(\theta\) a medida que se incrementa el numero de simulaciones por lo cual es considerado insesgado.
Tabla 1. Media de cada estimador por cantidad de simulaciones.
library(knitr)
library(kableExtra)
library(reshape2)
# Función para calcular la insesgadez
calcular_insesgadez <- function(estimadores, theta) {
colMeans(estimadores)
}
# Obtener la insesgadez para diferentes tamaños de muestra
insesgadez_20 <- calcular_insesgadez(muestras_20, theta)
insesgadez_50 <- calcular_insesgadez(muestras_50, theta)
insesgadez_100 <- calcular_insesgadez(muestras_100, theta)
insesgadez_1000 <- calcular_insesgadez(muestras_1000, theta)
# Crear un data.frame con insesgadez
insesgadez_df <- data.frame(
Tamano = rep(c(20, 50, 100, 1000), each = 4),
Insesgadez = c(insesgadez_20, insesgadez_50, insesgadez_100, insesgadez_1000),
Estimador = rep(c("Estimador 1", "Estimador 2", "Estimador 3", "Estimador 4"), times = 4)
)
# Convertir el data.frame a formato ancho
insesgadez_crosstab <- dcast(insesgadez_df, Tamano ~ Estimador, value.var = "Insesgadez")
# Mostrar la tabla usando kable y kableExtr.
kable(insesgadez_crosstab, format = "html", digits = 6) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = F, position = "center")| Tamano | Estimador 1 | Estimador 2 | Estimador 3 | Estimador 4 |
|---|---|---|---|---|
| 20 | 5.177122 | 10.21979 | 5.085644 | 6.003939 |
| 50 | 5.156130 | 10.32388 | 5.107465 | 5.685169 |
| 100 | 4.989776 | 10.10940 | 4.982349 | 5.844113 |
| 1000 | 5.019668 | 10.02908 | 5.025265 | 5.874587 |
La tabla 1 reafirma la propiedad de insesgadez de \(\hat{\theta}_1\) y \(\hat{\theta}_3\), pues se aprecia que entre mayor sea el número de simulaciones se aproxima más al valor real.Mientras que para \(\hat{\theta}_2\) y \(\hat{\theta}_4\) no se cumple esta afirmación, pues la cantidad de simulaciones no tienen efecto en el promedio calculado.
La eficiencia de un estimador hace referencia a la varianza de su distribución. Se dice que es mas eficiente si tiene una varianza menor respecto a los demás estimadores insesgados. Teniendo en cuenta que \(\hat{\theta}_1\) y \(\hat{\theta}_3\) son estimadores insesgados, se procederá a hallar la varianza de ambos para determinar cual es mas eficiente.
\(\text{Var}(\hat{\theta}_1) = \frac{5\theta^2}{18}\) y \(\text{Var}(\hat{\theta}_3)= \frac{\theta^2}{4}\)
Esto indica que \((\text{Var}(\hat{\theta}_1)\) > \(\text{Var}(\hat{\theta}_3))\), por lo que \((\hat{\theta}_3)\) tiene una menor varianza y, por lo tanto, es más eficiente que \((\hat{\theta}_1)\).Esta afirmación se ve respaldada por los resultados de la Tabla 2 y la Figura 1. Esto tiene relación con el estimador de máxima verosimilitud de una distribución exponencial que es el promedio muestral, puesto que \(\hat{\theta}_3\) utiliza el promedio.
Tabla 2. Varianza de cada estimador por cantidad de simulaciones.
library(knitr)
library(kableExtra)
library(reshape2)
# Función para varianza de los estimadores
calcular_varianza <- function(estimadores) {
apply(estimadores, 2, var)
}
# Se calculan las distintas varianzas
varianza_20 <- calcular_varianza(muestras_20)
varianza_50 <- calcular_varianza(muestras_50)
varianza_100 <- calcular_varianza(muestras_100)
varianza_1000 <- calcular_varianza(muestras_1000)
# Crear un data frame con las varianzas
varianza_df <- data.frame(
Tamano = rep(c(20, 50, 100, 1000), each = 4),
Varianza = c(varianza_20, varianza_50, varianza_100, varianza_1000),
Estimador = rep(c("Estimador 1", "Estimador 2", "Estimador 3", "Estimador 4"), times = 4))
# Convertir el data frame a formato ancho usando reshape2
varianza_crosstab <- dcast(varianza_df, Tamano ~ Estimador, value.var = "Varianza")
# Mostrar la tabla usando kable y kableExtra
kable(varianza_crosstab, format = "html", escape = FALSE, digits = 6) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),full_width = FALSE, position = "center")| Tamano | Estimador 1 | Estimador 2 | Estimador 3 | Estimador 4 |
|---|---|---|---|---|
| 20 | 6.803333 | 21.92193 | 4.964954 | 13.018168 |
| 50 | 5.079508 | 21.18270 | 4.394524 | 6.360973 |
| 100 | 6.791158 | 27.55937 | 6.004551 | 8.874438 |
| 1000 | 7.352848 | 31.80251 | 6.877106 | 10.572378 |
#Se realiza un grafico interactivo usando ggploty de las varianzas por cantidad de simulaciones para ver la eficiencia de manera sencilla.
library(ggplot2)
library(plotly)
# Crear el gráfico con ggplot2
p <- ggplot(varianza_df, aes(x = Tamano, y = Varianza, color = Estimador, shape = Estimador)) +
geom_line() +
geom_point(size = 3) +
scale_x_log10() +
scale_y_log10() +
labs(
x = "Tamaño de Muestra",
y = "Varianza"
) +
theme_minimal() +
theme(
legend.position = "bottom",
axis.title = element_text(size = 9),
axis.text = element_text(size = 9),
plot.title = element_blank()
)
ggplotly(p)Fig 1. Varianza de los estimadores por cantidad de simulaciones.
La propiedad de consistencia de un estimador hace referencia a la propiedad que describe la capacidad de un estimador de aproximarse al valor verdadero del parámetro a medida que el tamaño de la muestra crece. Se dice que un estimador es consistente si a medida que aumenta el tamaño de la muestra tiende al valor del parámetro y su varianza tiende a cero.
A continuación se realizará la generalización de los estimadores iniciales para cualquier tamaño de muestra n, esto con el objetivo de calcular su consistencia.
\(\hat{\theta}_1 = \frac{1}{2(n-1)} \sum_{i=1}^{\frac{n}{2}} X_i + \frac{1}{n-1} \sum_{i=\frac{n}{2} + 1}^{n} X_i\)
\(\hat{\theta}_2 = \frac{\sum_{i=1}^{n} C_i \cdot X_i}{n+1}\)
\(\hat{\theta}_3 = \frac{1}{n} \sum_{i=1}^{n} X_i\)
\(\hat{\theta}_4 = \frac{\text{min}(X_i) + \text{max}(X_i)}{2}\)
#Se realizan las funciones para generalizar los estimadores y ver su consistencia.
# Estimador 1 generalizado
estimador1a <- function(X) {
n <- length(X)
n_mitad <- n / 2
# Calcular las sumas de los bloques
bloque1 <- sum(X[1:n_mitad]) / (2 * (n - 1))
bloque2 <- sum(X[(n_mitad + 1):n]) / (n - 1)
# Sumar ambos bloques para obtener el estimador final
return(bloque1 + bloque2) }
#
estimador2a <- function(X) {
n <- length(X)
coeficientes <- 1:n
return(sum(coeficientes * X) / (1+n))}
# Estimador 3 generalizado (media muestral)
estimador3a <- function(X) {
return(mean(X)) }
# Estimador 4 generalizado (promedio del mínimo y máximo)
estimador4a <- function(X) {
return((min(X) + max(X)) / 2) }
# Función para simular datos y calcular estimadores
# `n` es el número de muestras, `size` es el tamaño de cada muestra
simular_estimadores <- function(n, size, theta) {
muestras <- matrix(rexp(size * n, rate = 1 / theta), ncol = size)
pesos1 <- rep(1, size)
pesos2 <- 1:size
estimadores <- data.frame(
Estimador1 = apply(muestras, 2, estimador1a),
Estimador2 = apply(muestras, 2, estimador2a),
Estimador3 = apply(muestras,2,estimador3a),
Estimador4 = apply(muestras, 1, estimador4a)
)
return(estimadores)
}
# Simular
muestras_20a <- simular_estimadores(20, 10000, theta)
muestras_50a <- simular_estimadores(50, 10000, theta)
muestras_100a <- simular_estimadores(100, 10000, theta)
muestras_1000a <- simular_estimadores(1000, 10000, theta)Tabla 3. Varianza para distintos tamaños de muestra y 10.000 simulaciones.
#Se calculan las varianzas para la generalización de los estimadores para evidenciar la propiedad de consistencia.
calcular_varianza <- function(estimadores) {
apply(estimadores, 2, var)
}
# Obtener varianzas
varianza_20a <- calcular_varianza(muestras_20a)
varianza_50a <- calcular_varianza(muestras_50a)
varianza_100a <- calcular_varianza(muestras_100a)
varianza_1000a <- calcular_varianza(muestras_1000a)
# Crear un data.frame con varianza
varianza_df <- data.frame(
Tamano = rep(c(20, 50, 100, 1000), each = 4),
Varianza = c(varianza_20a, varianza_50a, varianza_100a, varianza_1000a),
Estimador = rep(c("Estimador1", "Estimador2", "Estimador3", "Estimador4"), times = 4)
)
# Convertir el data.frame a formato ancho
varianza_crosstab <- dcast(varianza_df, Tamano ~ Estimador, value.var = "Varianza")
# Mostrar la tabla cruzada en formato kable
kable(varianza_crosstab, format = "html", table.attr = "style='width:50%;'") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
font_size = 12,
position = "center")| Tamano | Estimador1 | Estimador2 | Estimador3 | Estimador4 |
|---|---|---|---|---|
| 20 | 0.8571036 | 162.8809 | 1.2184126 | 12.382892 |
| 50 | 0.3204755 | 406.5021 | 0.4926747 | 12.161861 |
| 100 | 0.1600932 | 819.6182 | 0.2493378 | 14.476345 |
| 1000 | 0.0155735 | 8291.8879 | 0.0247583 | 9.934573 |
Como se evidencia en la Tabla 3, los estimadores consistentes son \(\hat{\theta}_1\) y \(\hat{\theta}_3\) pues su varianza a medida que aumenta el numero de muestras es cercana a cero.
Finalmente la figura 2 muestra como a medida que se aumenta la cantidad de simulaciones, la mediana de los estimadores se sitúa mas cerca del valor real al tiempo que su varianza se ve reducida. Los estimadores que disponen de un mejor sesgo son \(\hat{\theta}_1\) y \(\hat{\theta}_3\), al igual que el que presenta mayor variación
library(ggplot2)
library(plotly)
library(reshape2)
# Crear un data frame que contenga todas las muestras
muestras_totales <- rbind(
data.frame(Tamano = "n = 20", muestras_20),
data.frame(Tamano = "n = 50", muestras_50),
data.frame(Tamano = "n = 100", muestras_100),
data.frame(Tamano = "n = 1000", muestras_1000)
)
# Convertir el data frame a formato largo para ggplot
muestras_largas <- melt(muestras_totales, id.vars = "Tamano", variable.name = "Estimador")
# Cambiar los nombres de las etiquetas de los estimadores
muestras_largas$Estimador <- factor(muestras_largas$Estimador,
levels = c("Estimador1", "Estimador2", "Estimador3", "Estimador4"),
labels = c(expression("Est.1"),
expression("Est.2"),
expression("Est.3"),
expression("Est.4")))
# Crear el boxplot usando ggplot2
p <- ggplot(muestras_largas, aes(x = Estimador, y = value, color = Estimador)) +
geom_boxplot() +
facet_wrap(~ Tamano, scales = "free") +
labs(
x = "Estimadores",
y = "Valor"
) +
theme_minimal() +
theme(
legend.position = "none",
axis.title = element_text(size = 9),
axis.text = element_text(size = 9, angle = 0, hjust = 1),
strip.text = element_text(size = 10),
panel.spacing = unit(2, "lines") # Espacio entre paneles
) +
geom_hline(yintercept = theta, color = "red", linetype = "dashed")
# Convertir a plotly para hacerlo interactivo.
ggplotly(p)Fig.2 Varianzas por estimador por número de simulaciones.