Este apunte muestra, por simulación en R, tres casos
donde la media muestral \(\bar X\) es estimador
eficiente (alcanza el Cramer–Rao Lower Bound,
CRLB) y dos familias donde \(\bar
X\) es consistente.
Incluye código ejecutable, gráficos y pequeñas guías para
interpretación.
Recordatorio: Un estimador es eficiente si su varianza alcanza el CRLB. Es consistente si converge en probabilidad al verdadero parámetro al aumentar \(n\).
mu <- 2
sigma <- 3 # conocida
n <- 80
R <- 20000
# Estimador candidato: media muestral
xbar <- replicate(R, mean(rnorm(n, mean = mu, sd = sigma)))
# Teoría: I(µ)=1/σ^2 ⇒ CRLB = σ^2/n
crlb_normal <- sigma^2 / n
var_emp <- var(xbar)
ecm_emp <- ECM(xbar, mu)
cat("CRLB teórico: ", round(crlb_normal, 6), "\n")
## CRLB teórico: 0.1125
## Var empírica(X̄): 0.112284
## ECM empírico(X̄): 0.112279
## Var empírica / CRLB = 0.998 (≈ 1 indica eficiencia)
Visualización
par(mfrow = c(1,2))
hist(xbar, breaks = 40, col = "gray", main = "Distribución de X̄",
xlab = "X̄", border = "white")
abline(v = mu, lwd = 2)
plot(density(xbar), main = "Densidad de X̄ vs µ", xlab = "X̄")
abline(v = mu, lwd = 2)
Comparación con un estimador insesgado pero ineficiente (\(X_1\))
x1 <- replicate(R, rnorm(1, mean = mu, sd = sigma))
cat("Var empírica(X1): ", round(var(x1), 6), " (≈ σ^2)\n")
## Var empírica(X1): 9.072962 (≈ σ^2)
lambda <- 4
n <- 120
R <- 20000
xbar_pois <- replicate(R, mean(rpois(n, lambda)))
# Teoría: I(λ)=1/λ ⇒ CRLB = λ/n
crlb_pois <- lambda / n
var_emp_p <- var(xbar_pois)
ecm_emp_p <- ECM(xbar_pois, lambda)
cat("CRLB teórico: ", round(crlb_pois, 6), "\n")
## CRLB teórico: 0.033333
## Var empírica(X̄): 0.033083
## ECM empírico(X̄): 0.033081
## Var empírica / CRLB = 0.992 (≈ 1 indica eficiencia)
par(mfrow = c(1,2))
hist(xbar_pois, breaks = 40, col = "gray", main = "Poisson: X̄",
xlab = "X̄", border = "white")
abline(v = lambda, lwd = 2)
plot(density(xbar_pois), main = "Densidad de X̄ vs λ", xlab = "X̄")
abline(v = lambda, lwd = 2)
p <- 0.35
n <- 200
R <- 20000
xbar_bern <- replicate(R, mean(rbinom(n, size = 1, prob = p)))
# Teoría: I(p)=1/[p(1-p)] ⇒ CRLB = p(1-p)/n
crlb_bern <- p * (1 - p) / n
var_emp_b <- var(xbar_bern)
ecm_emp_b <- ECM(xbar_bern, p)
cat("CRLB teórico: ", round(crlb_bern, 6), "\n")
## CRLB teórico: 0.001137
## Var empírica(X̄): 0.001139
## ECM empírico(X̄): 0.001139
## Var empírica / CRLB = 1.001 (≈ 1 indica eficiencia)
par(mfrow = c(1,2))
hist(xbar_bern, breaks = 40, col = "gray", main = "Bernoulli: X̄",
xlab = "X̄", border = "white")
abline(v = p, lwd = 2)
plot(density(xbar_bern), main = "Densidad de X̄ vs p", xlab = "X̄")
abline(v = p, lwd = 2)
A continuación verificamos la consistencia por simulación comparando \(\bar X\) con \(X_1\) (el primer dato), que no es consistente.
prob_within_epsilon <- function(estimator_values, theta, eps){
mean(abs(estimator_values - theta) <= eps)
}
ECM <- function(estimator_values, theta){
mean((estimator_values - theta)^2)
}
set.seed(123)
mu <- 10
sigma <- 3
n_grid <- c(5,10,20,50,100,200,500,1000)
R <- 5000
epsilon <- 0.3
var_xbar <- var_x1 <- prob_eps_xbar <- prob_eps_x1 <- ecm_xbar <- ecm_x1 <- numeric(length(n_grid))
for(k in seq_along(n_grid)){
n <- n_grid[k]
xbar <- x1 <- numeric(R)
for(r in 1:R){
x <- rnorm(n, mean=mu, sd=sigma)
xbar[r] <- mean(x)
x1[r] <- x[1]
}
var_xbar[k] <- var(xbar)
var_x1[k] <- var(x1)
prob_eps_xbar[k] <- prob_within_epsilon(xbar, mu, epsilon)
prob_eps_x1[k] <- prob_within_epsilon(x1, mu, epsilon)
ecm_xbar[k] <- ECM(xbar, mu)
ecm_x1[k] <- ECM(x1, mu)
}
par(mfrow = c(1,3))
plot(n_grid, var_xbar, type="b", pch=19, xlab="n", ylab="Varianza",
main="Normal: Varianzas de X̄ vs X1"); lines(n_grid, var_x1, type="b", pch=17)
legend("topright", legend=c("Var(X̄)","Var(X1)"), pch=c(19,17))
plot(n_grid, prob_eps_xbar, type="b", pch=19, ylim=c(0,1), xlab="n",
ylab=paste0("P(|Tn - µ| ≤ ", epsilon, ")"),
main="Normal: Probabilidad de estar cerca de µ")
lines(n_grid, prob_eps_x1, type="b", pch=17); abline(h=1, lty=2, col="gray")
legend("bottomright", legend=c("X̄","X1"), pch=c(19,17))
plot(n_grid, ecm_xbar, type="b", pch=19, xlab="n", ylab="ECM",
main="Normal: ECM de X̄ vs X1"); lines(n_grid, ecm_x1, type="b", pch=17)
legend("topright", legend=c("ECM(X̄)","ECM(X1)"), pch=c(19,17))
set.seed(123)
p <- 0.3
n_grid <- c(5,10,20,50,100,200,500,1000)
R <- 5000
epsilon_b <- 0.05
var_xbar_b <- var_x1_b <- prob_eps_xbar_b <- prob_eps_x1_b <- ecm_xbar_b <- ecm_x1_b <- numeric(length(n_grid))
for(k in seq_along(n_grid)){
n <- n_grid[k]
xbar <- x1 <- numeric(R)
for(r in 1:R){
x <- rbinom(n, size=1, prob=p)
xbar[r] <- mean(x)
x1[r] <- x[1]
}
var_xbar_b[k] <- var(xbar)
var_x1_b[k] <- var(x1)
prob_eps_xbar_b[k] <- prob_within_epsilon(xbar, p, epsilon_b)
prob_eps_x1_b[k] <- prob_within_epsilon(x1, p, epsilon_b)
ecm_xbar_b[k] <- ECM(xbar, p)
ecm_x1_b[k] <- ECM(x1, p)
}
par(mfrow = c(1,3))
plot(n_grid, var_xbar_b, type="b", pch=19, xlab="n", ylab="Varianza",
main="Bernoulli: Varianzas de X̄ vs X1"); lines(n_grid, var_x1_b, type="b", pch=17)
legend("topright", legend=c("Var(X̄)","Var(X1)"), pch=c(19,17))
plot(n_grid, prob_eps_xbar_b, type="b", pch=19, ylim=c(0,1), xlab="n",
ylab=paste0("P(|Tn - p| ≤ ", epsilon_b, ")"),
main="Bernoulli: Probabilidad de estar cerca de p")
lines(n_grid, prob_eps_x1_b, type="b", pch=17); abline(h=1, lty=2, col="gray")
legend("bottomright", legend=c("X̄","X1"), pch=c(19,17))
plot(n_grid, ecm_xbar_b, type="b", pch=19, xlab="n", ylab="ECM",
main="Bernoulli: ECM de X̄ vs X1"); lines(n_grid, ecm_x1_b, type="b", pch=17)
legend("topright", legend=c("ECM(X̄)","ECM(X1)"), pch=c(19,17))
## R version 4.3.2 (2023-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 11 x64 (build 26100)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=Spanish_Argentina.utf8 LC_CTYPE=Spanish_Argentina.utf8
## [3] LC_MONETARY=Spanish_Argentina.utf8 LC_NUMERIC=C
## [5] LC_TIME=Spanish_Argentina.utf8
##
## time zone: America/Buenos_Aires
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## loaded via a namespace (and not attached):
## [1] digest_0.6.33 R6_2.5.1 fastmap_1.1.1 xfun_0.52
## [5] cachem_1.0.8 knitr_1.50 htmltools_0.5.7 rmarkdown_2.29
## [9] lifecycle_1.0.4 cli_3.6.2 sass_0.4.8 jquerylib_0.1.4
## [13] compiler_4.3.2 rstudioapi_0.15.0 tools_4.3.2 evaluate_0.23
## [17] bslib_0.6.1 yaml_2.3.8 rlang_1.1.2 jsonlite_1.8.8