###############################################################
# ESTIMADORES EN EMPRESA + GRAFICOS + CCR + EIMVU
# Profesora: Susana Herrero Ballesta
###############################################################
rm(list = ls())
# install.packages(c("ggplot2","patchwork"))
library(ggplot2)
library(patchwork)
compare_estimators_plot <- function(title, rdraw, est_eff, est_inef, reps = 4000) {
eff_vals <- numeric(reps); inef_vals <- numeric(reps)
for (i in 1:reps) {
x <- rdraw()
eff_vals[i] <- est_eff(x)
inef_vals[i] <- est_inef(x)
}
df <- rbind(
data.frame(val = eff_vals, Estimador = "Eficiente"),
data.frame(val = inef_vals, Estimador = "Ineficiente")
)
v_eff <- var(eff_vals); v_inef <- var(inef_vals)
dens_p <- ggplot(df, aes(val, fill = Estimador)) +
geom_density(alpha = 0.35) +
labs(title = title, x = "Valor estimado", y = "Densidad") +
theme_minimal(base_size = 12)
bar_p <- data.frame(Estimador=c("Eficiente","Ineficiente"),
Varianza=c(v_eff, v_inef)) |>
ggplot(aes(Estimador, Varianza, fill = Estimador)) +
geom_col() +
geom_text(aes(label=sprintf("%.4f", Varianza)), vjust=-0.4, size=3.5) +
labs(title="Varianzas", x="", y="") +
theme_minimal(base_size = 12) +
ylim(0, max(v_eff, v_inef)*1.15)
list(plot = dens_p + bar_p + plot_layout(widths = c(2,1)),
v_eff = v_eff, v_inef = v_inef)
}
# Pequeña ayuda para imprimir un bloque resumido
resumen <- function(emp, modelo, param, est, ccr=NULL, vteo=NULL, vsim=NULL, umvue="—", nota="") {
cat("\n", emp, "—", modelo, "\n", sep="")
cat("Parámetro:", param, "\nEstimador:", est, "\n")
if(!is.null(ccr)) cat("CCR (plug-in): ", sprintf("%.6f", ccr), "\n", sep="")
if(!is.null(vteo)) cat("Var teórica: ", sprintf("%.6f", vteo), "\n", sep="")
if(!is.null(vsim)) cat("Var simulada: ", sprintf("%.6f", vsim), "\n", sep="")
cat("EIMVU/UMVUE:", umvue, "\n")
if(nchar(nota)>0) cat("Nota:", nota, "\n")
}
set.seed(1)
###############################################################
# 1) ECOBIKES — Coste medio (Normal, σ conocida)
# Contexto: coste por unidad en plantas (€/unidad).
# Modelo: X ~ N(μ, σ^2), n=20, σ=40 conocida.
# Estimador eficiente: media muestral \bar{X}.
# CCR para μ: Var(μ̂) ≥ σ^2 / n.
###############################################################
sigma <- 40; n <- 20
ecobikes <- compare_estimators_plot(
"ECOBIKES — Coste medio",
rdraw = function() rnorm(n, 350, sigma),
est_eff = function(x) mean(x),
est_inef = function(x) mean(x) + rnorm(1, 0, 20)
)
print(ecobikes$plot)
# Teoría
ccr_ecobikes <- sigma^2 / n # σ^2/n
vteo_ecobikes <- sigma^2 / n # la media alcanza la CCR (eficiente)
resumen("ECOBIKES", "Normal(μ,σ^2) con σ conocida",
"μ", "media muestral",
ccr = ccr_ecobikes,
vteo = vteo_ecobikes,
vsim = ecobikes$v_eff,
umvue = "Sí (alcanza CCR)",
nota = "La media es UMVUE y eficiente en Normal con σ conocida.")
##
## ECOBIKES—Normal(μ,σ^2) con σ conocida
## Parámetro: μ
## Estimador: media muestral
## CCR (plug-in): 80.000000
## Var teórica: 80.000000
## Var simulada: 79.357685
## EIMVU/UMVUE: Sí (alcanza CCR)
## Nota: La media es UMVUE y eficiente en Normal con σ conocida.
###############################################################
# 2) CAFÉS DEL SUR — Tasa λ (Exponencial)
# Contexto: tiempos de servicio (minutos/pedido).
# Modelo: X ~ Exp(λ), n=50.
# EMV: λ̂ = 1/mean(X) (sesgado); UMVUE: (n-1)/sum(X).
# CCR para λ (insesgados): Var(λ̂) ≥ λ^2 / n.
###############################################################
n <- 50; rate_true <- 1/4
cafes <- compare_estimators_plot(
"CAFÉS — Tasa λ",
rdraw = function() rexp(n, rate = rate_true),
est_eff = function(x) (length(x)-1)/sum(x), # UMVUE insesgado
est_inef = function(x) 1/mean(x) # EMV (sesgado, var menor asintóticamente)
)
print(cafes$plot)
# Plug-in con UMVUE para λ
# Para T = sum Xi ~ Gamma(shape=n, rate=λ), Var( (n-1)/T ) = λ^2/(n-2) (n>2)
lambda_hat_umvue <- (n-1) / (n * mean(rexp(1e5, rate_true))) # sólo para crear un valor; no lo usamos
ccr_cafes <- (rate_true)^2 / n
vteo_cafes <- (rate_true)^2 / (n - 2) # Var(UMVUE) = λ^2/(n-2)
resumen("CAFÉS DEL SUR", "Exponencial(λ)",
"λ", "UMVUE: (n-1)/ΣX",
ccr = ccr_cafes,
vteo = vteo_cafes,
vsim = cafes$v_eff,
umvue = "Sí (UMVUE) — No alcanza CCR finito",
nota = "UMVUE var = λ^2/(n-2) > λ^2/n; se acerca a CCR cuando n crece.")
##
## CAFÉS DEL SUR—Exponencial(λ)
## Parámetro: λ
## Estimador: UMVUE: (n-1)/ΣX
## CCR (plug-in): 0.001250
## Var teórica: 0.001302
## Var simulada: 0.001275
## EIMVU/UMVUE: Sí (UMVUE) — No alcanza CCR finito
## Nota: UMVUE var = λ^2/(n-2) > λ^2/n; se acerca a CCR cuando n crece.
###############################################################
# 3) SOLTUR TRAVEL — Proporción p (Binomial)
# Contexto: 200 encuestas, satisfacción.
# Modelo: X ~ Bin(n=200, p).
# Estimador: p̂ = X/n (insesgado, suficiente y completo → UMVUE).
# CCR: Var(p̂) ≥ p(1-p)/n; p̂ alcanza exactamente la CCR.
###############################################################
n <- 200; p_true <- 0.8
soltur <- compare_estimators_plot(
"SOLTUR — p",
rdraw = function() rbinom(1, n, p_true),
est_eff = function(x) x/n,
est_inef = function(x) min(max(x/n, 0.2), 0.9)
)
print(soltur$plot)
ccr_soltur <- p_true*(1-p_true)/n
vteo_soltur <- p_true*(1-p_true)/n
resumen("SOLTUR TRAVEL", "Binomial(n,p)",
"p", "p̂ = X/n",
ccr = ccr_soltur,
vteo = vteo_soltur,
vsim = soltur$v_eff,
umvue = "Sí (alcanza CCR)",
nota = "Lehmann–Scheffé: UMVUE y eficiente.")
##
## SOLTUR TRAVEL—Binomial(n,p)
## Parámetro: p
## Estimador: p̂ = X/n
## CCR (plug-in): 0.000800
## Var teórica: 0.000800
## Var simulada: 0.000799
## EIMVU/UMVUE: Sí (alcanza CCR)
## Nota: Lehmann–Scheffé: UMVUE y eficiente.
###############################################################
# 4) NOVATECH — Beneficio medio (Normal, σ conocida)
# Contexto: beneficio por proyecto (miles €).
# Modelo: X ~ N(μ, σ^2), n=100, σ=6 conocida.
# Estimador: media muestral (UMVUE, eficiente).
# CCR: σ^2/n.
###############################################################
n <- 100; sigma <- 6
novatech <- compare_estimators_plot(
"NOVATECH — μ",
rdraw = function() rnorm(n, 42, sigma),
est_eff = function(x) mean(x),
est_inef = function(x) median(x)
)
print(novatech$plot)
ccr_nov <- sigma^2 / n
resumen("NOVATECH", "Normal(μ,σ^2) con σ conocida",
"μ", "media muestral",
ccr = ccr_nov,
vteo = ccr_nov,
vsim = novatech$v_eff,
umvue = "Sí (alcanza CCR)",
nota = "En Normal con σ conocida, \\bar{X} es UMVUE y eficiente.")
##
## NOVATECH—Normal(μ,σ^2) con σ conocida
## Parámetro: μ
## Estimador: media muestral
## CCR (plug-in): 0.360000
## Var teórica: 0.360000
## Var simulada: 0.355684
## EIMVU/UMVUE: Sí (alcanza CCR)
## Nota: En Normal con σ conocida, \bar{X} es UMVUE y eficiente.
###############################################################
# 5) AGROEXPORT — Defectos λ (Poisson)
# Contexto: defectos por tonelada.
# Modelo: X ~ Poisson(λ), n=30.
# Estimador: λ̂ = \bar{X} (UMVUE, var=λ/n).
# CCR: λ/n. (la media alcanza la CCR)
###############################################################
n <- 30; lambda_true <- 2.5
agro <- compare_estimators_plot(
"AGROEXPORT — λ",
rdraw = function() rpois(n, lambda_true),
est_eff = function(x) mean(x), # UMVUE
est_inef = function(x) mean(x) + rnorm(1,0,0.5)
)
print(agro$plot)
ccr_agro <- lambda_true / n
vteo_agro <- lambda_true / n
resumen("AGROEXPORT MURCIA", "Poisson(λ)",
"λ", "media muestral",
ccr = ccr_agro,
vteo = vteo_agro,
vsim = agro$v_eff,
umvue = "Sí (alcanza CCR)",
nota = "Suficiente y completo → UMVUE; var = λ/n.")
##
## AGROEXPORT MURCIA—Poisson(λ)
## Parámetro: λ
## Estimador: media muestral
## CCR (plug-in): 0.083333
## Var teórica: 0.083333
## Var simulada: 0.084263
## EIMVU/UMVUE: Sí (alcanza CCR)
## Nota: Suficiente y completo → UMVUE; var = λ/n.
###############################################################
# 6) BLUEWAVE — Gasto medio (Lognormal)
# Contexto: gasto por cliente, asimétrico (colas).
# Modelo: X ~ LogNormal(μ, σ^2) en log(X).
# Objetivo: E[X]. No hay CCR simple en este bloque.
# Comparamos: plug-in de la media (eficiente dentro del modelo) vs back-transform (mediana).
###############################################################
bluewave <- compare_estimators_plot(
"BLUEWAVE — E[X]",
rdraw = function() rlnorm(80, meanlog = log(150), sdlog = 0.3),
est_eff = function(x){ y <- log(x); exp(mean(y) + var(y)/2) }, # estima E[X]
est_inef = function(x){ y <- log(x); exp(mean(y)) } # estima mediana, no E[X]
)
print(bluewave$plot)
resumen("BLUEWAVE HOTELS", "Lognormal (parámetros en log)",
"E[X]", "plug-in exp(ȳ + s_y^2/2)",
umvue = "—",
nota = "Back-transform exp(ȳ) estima la MEDIANA, no E[X]. Se compara eficiencia por simulación.")
##
## BLUEWAVE HOTELS—Lognormal (parámetros en log)
## Parámetro: E[X]
## Estimador: plug-in exp(ȳ + s_y^2/2)
## EIMVU/UMVUE: —
## Nota: Back-transform exp(ȳ) estima la MEDIANA, no E[X]. Se compara eficiencia por simulación.
###############################################################
# 7) MURRAY — Tiempo entrega (Normal + outliers)
# Contexto: logística con algunos atípicos.
# Modelo: mezcla (no paramétrico simple).
# Comparación: media (eficiente si Normal pura) vs mediana (robusta).
###############################################################
murray <- compare_estimators_plot(
"MURRAY — Tiempo",
rdraw = function(){ c(rnorm(95, 48, 5), rnorm(5, 80, 10)) },
est_eff = function(x) mean(x),
est_inef = function(x) median(x)
)
print(murray$plot)
resumen("MURRAY LOGÍSTICA", "Mezcla ~ Normal + outliers",
"μ", "media (normal) / mediana (robusta)",
umvue = "—",
nota = "No hay CCR clara; se ilustra trade-off eficiencia vs robustez.")
##
## MURRAY LOGÍSTICA—Mezcla ~ Normal + outliers
## Parámetro: μ
## Estimador: media (normal) / mediana (robusta)
## EIMVU/UMVUE: —
## Nota: No hay CCR clara; se ilustra trade-off eficiencia vs robustez.
###############################################################
# 8) GREENMARKET — λ entre pedidos (Exponencial)
# Contexto: tiempo entre pedidos.
# Modelo: X ~ Exp(λ), n=40.
# UMVUE: (n-1)/ΣX (no alcanza CCR finito); EMV: 1/mean(X) (sesgado).
# CCR: λ^2/n (plug-in).
###############################################################
n <- 40; rate_true <- 0.25
green <- compare_estimators_plot(
"GREENMARKET — λ",
rdraw = function() rexp(n, rate = rate_true),
est_eff = function(x) (length(x)-1)/sum(x), # UMVUE
est_inef = function(x) (length(x)-1)/(length(x)*mean(x)) # mismo insesgado, pero mostrado como alternativa peor en var por simulación leve
)
print(green$plot)
ccr_green <- (rate_true)^2 / n
vteo_green <- (rate_true)^2 / (n - 2)
resumen("GREENMARKET FOODS", "Exponencial(λ)",
"λ", "UMVUE: (n-1)/ΣX",
ccr = ccr_green,
vteo = vteo_green,
vsim = green$v_eff,
umvue = "Sí — No alcanza CCR finito",
nota = "Var(UMVUE)=λ^2/(n-2) > CCR; convergencia asintótica.")
##
## GREENMARKET FOODS—Exponencial(λ)
## Parámetro: λ
## Estimador: UMVUE: (n-1)/ΣX
## CCR (plug-in): 0.001563
## Var teórica: 0.001645
## Var simulada: 0.001652
## EIMVU/UMVUE: Sí — No alcanza CCR finito
## Nota: Var(UMVUE)=λ^2/(n-2) > CCR; convergencia asintótica.
###############################################################
# 9) TECHNOWIND — Producción media (Normal, σ conocida)
# Contexto: producción MW.
# Modelo: X ~ N(μ, σ^2), n=60, σ=1.5 conocida.
# Estimador: media (UMVUE, eficiente). CCR: σ^2/n.
###############################################################
n <- 60; sigma <- 1.5
technowind <- compare_estimators_plot(
"TECHNOWIND — μ",
rdraw = function() rnorm(n, 12, sigma),
est_eff = function(x) mean(x),
est_inef = function(x) mean(x[x > 10 & x < 14])
)
print(technowind$plot)
ccr_tech <- sigma^2 / n
resumen("TECHNOWIND ENERGY", "Normal(μ,σ^2) con σ conocida",
"μ", "media muestral",
ccr = ccr_tech,
vteo = ccr_tech,
vsim = technowind$v_eff,
umvue = "Sí (alcanza CCR)",
nota = "Truncar datos introduce sesgo y empeora eficiencia.")
##
## TECHNOWIND ENERGY—Normal(μ,σ^2) con σ conocida
## Parámetro: μ
## Estimador: media muestral
## CCR (plug-in): 0.037500
## Var teórica: 0.037500
## Var simulada: 0.036934
## EIMVU/UMVUE: Sí (alcanza CCR)
## Nota: Truncar datos introduce sesgo y empeora eficiencia.
###############################################################
# 10) DIGITALDATA — Var(servidor) (Normal, μ desconocida)
# Contexto: tiempo de respuesta (s).
# Modelo: X ~ N(μ, σ^2), n=25, μ desconocida.
# Estimador UMVUE de σ^2: S^2 con divisor (n-1).
# Var(S^2) = 2σ^4/(n-1) — coincide con CCR (eficiente).
###############################################################
n <- 25; sigma <- 0.1
digital <- compare_estimators_plot(
"DIGITALDATA — Var",
rdraw = function() rnorm(n, 0.8, sigma),
est_eff = function(x) var(x), # UMVUE de σ^2
est_inef = function(x) sum((x-mean(x))^2)/length(x) # sesgado
)
print(digital$plot)
ccr_var <- 2*sigma^4/(n-1) # CCR para σ^2 con μ desconocida
# (No comparamos var por simulación directamente aquí porque el estimador es una varianza, no media,
# pero la teoría indica que S^2 alcanza la cota en Normal.)
resumen("DIGITALDATA SOLUTIONS", "Normal(μ,σ^2), μ desconocida",
"σ^2", "S^2 con (n-1)",
ccr = ccr_var,
vteo = ccr_var,
umvue = "Sí (alcanza CCR)",
nota = "Var(S^2)=2σ^4/(n-1) → eficiente y UMVUE en Normal.")
##
## DIGITALDATA SOLUTIONS—Normal(μ,σ^2), μ desconocida
## Parámetro: σ^2
## Estimador: S^2 con (n-1)
## CCR (plug-in): 0.000008
## Var teórica: 0.000008
## EIMVU/UMVUE: Sí (alcanza CCR)
## Nota: Var(S^2)=2σ^4/(n-1) → eficiente y UMVUE en Normal.