Parsimonia_AIC

Prompt 1: El problema y la Verosimilitud

“Actúa como mi tutor de estadística para Ingeniería Agrícola. Explícame en máximo 3 viñetas (bullet points) qué es la Verosimilitud (Likelihood) y por qué es necesaria evaluarla cuando varias distribuciones (como Gamma y Lognormal) parecen ajustarse bien a mis datos históricos de caudales máximos anuales. Usa un lenguaje técnico pero claro, ideal para tomar notas a mano.”

El desempate entre distribuciones

Se resuelve con el Principio de Parsimonia (Navaja de Occam), que dicta preferir el modelo más simple cuando dos explican los datos igual de bien, porque un modelo complejo puede sobreajustar ruido en lugar de captar la tendencia real. El Criterio de Información de Akaike (AIC) formaliza esta idea combinando bondad de ajuste y complejidad mediante la fórmula AIC = –2·log(verosimilitud máxima) + 2·k, donde k es el número de parámetros; cuanto menor sea el AIC, mejor es el modelo. En el campo, es como diseñar un canal de riego: si un canal rectilíneo con pendiente uniforme (menos parámetros) conduce el mismo caudal que uno con muchas curvas y ajustes, el rectilíneo gana por simplicidad y facilidad de mantenimiento, y el AIC sería como el presupuesto total que penaliza cada pieza extra. La regla de oro matemática para saber qué distribución gana es calcular ΔAIC = AIC(modelo más alto) – AIC(modelo más bajo); si ΔAIC > 2, el modelo con menor AIC es claramente superior, mientras que si ΔAIC < 2, no hay diferencia significativa y ambos son igualmente plausibles.

# 1. Definición de la función para calcular la verosimilitud de caudales máximos anuales
evaluar_verosimilitud_caudales <- function(caudales) {
  # Filtro para asegurar valores positivos (requerido por Gamma y Lognormal)
  caudales <- caudales[caudales > 0 & !is.na(caudales)]
  n <- length(caudales)
  
  # --- Ajuste Distribución Gamma ---
  # Estimadores de momentos como valores iniciales para optimización numérica
  mean_x <- mean(caudales)
  var_x <- var(caudales)
  shape_init <- mean_x^2 / var_x
  rate_init <- mean_x / var_x
  
  fit_gamma <- optim(
    par = c(shape = shape_init, rate = rate_init),
    fn = function(p) {
      if (p[1] <= 0 || p[2] <= 0) return(1e10)
      -sum(dgamma(caudales, shape = p[1], rate = p[2], log = TRUE))
    }
  )
  
  logLik_gamma <- -fit_gamma$value
  k_gamma <- 2
  aic_gamma <- 2 * k_gamma - 2 * logLik_gamma
  bic_gamma <- k_gamma * log(n) - 2 * logLik_gamma
  
  # --- Ajuste Distribución Lognormal ---
  # Estimadores MLE analíticos exactos
  meanlog_mle <- mean(log(caudales))
  sdlog_mle <- sqrt(mean((log(caudales) - meanlog_mle)^2))
  
  logLik_lognorm <- sum(dlnorm(caudales, meanlog = meanlog_mle, sdlog = sdlog_mle, log = TRUE))
  k_lognorm <- 2
  aic_lognorm <- 2 * k_lognorm - 2 * logLik_lognorm
  bic_lognorm <- k_lognorm * log(n) - 2 * logLik_lognorm
  
  # --- Construcción de la Tabla Comparativa ---
  tabla_resultados <- data.frame(
    Modelo = c("Gamma", "Lognormal"),
    Parametro_1 = c(paste0("shape = ", round(fit_gamma$par[1], 4)), 
                    paste0("meanlog = ", round(meanlog_mle, 4))),
    Parametro_2 = c(paste0("rate = ", round(fit_gamma$par[2], 4)), 
                    paste0("sdlog = ", round(sdlog_mle, 4))),
    Log_Likelihood = round(c(logLik_gamma, logLik_lognorm), 2),
    AIC = round(c(aic_gamma, aic_lognorm), 2),
    BIC = round(c(bic_gamma, bic_lognorm), 2),
    stringsAsFactors = FALSE
  )
  
  return(tabla_resultados)
}

# 2. Ejemplo de uso con caudales máximos anuales ficticios (m³/s)
set.seed(123)
caudales_observados <- rgamma(n = 40, shape = 3.5, rate = 0.01)

# 3. Generación de la tabla
tabla_comparativa <- evaluar_verosimilitud_caudales(caudales_observados)
print(tabla_comparativa)

##      Modelo      Parametro_1    Parametro_2 Log_Likelihood    AIC    BIC
## 1     Gamma    shape = 3.753  rate = 0.0126        -254.41 512.81 516.19
## 2 Lognormal meanlog = 5.5582 sdlog = 0.5577        -255.73 515.46 518.83

Prompt 2: El desempate (Parsimonia y AIC)

“Ahora necesito anotar el concepto de desempate. Explícame brevemente el Principio de Parsimonia y el Criterio de Información de Akaike (AIC). Dame una analogía sencilla relacionada con el diseño de una obra en el campo (ej. un canal o un muro) y termina con la ‘regla de oro’ matemática para saber qué distribución gana según el AIC. Estructúralo para copiarlo en mi cuaderno.” **

La verosimilitud (likelihood)

Es una función que mide qué tan plausibles son unos parámetros dados unos datos observados, calculándose como el producto de las densidades de probabilidad (o su logaritmo) para cada valor histórico de caudales máximos anuales; a diferencia de la probabilidad, que predice datos futuros, la verosimilitud evalúa hacia atrás: “dado que vi estos datos, ¿qué parámetros los hacen más probables?”. Es necesaria evaluarla cuando varias distribuciones (como Gamma y Lognormal) parecen ajustarse bien visualmente, porque histogramas o Q-Q plots no cuantifican la incertidumbre de los parámetros; al maximizar la verosimilitud (MLE) se obtienen los parámetros óptimos de cada distribución y, mediante criterios como AIC o BIC que penalizan la complejidad, se decide objetivamente cuál modelo es más verosímil. En la práctica, para caudales extremos, una elección errónea entre estas distribuciones puede subestimar o sobreestimar los caudales de diseño, por lo que la verosimilitud permite detectar, por ejemplo, si la Lognormal asigna baja probabilidad a los valores extremos observados, favoreciendo así a la Gamma u otra distribución con mayor sustento estadístico.

# ============================================================
# DESEMPATE ENTRE DISTRIBUCIONES (Gamma vs Lognormal)
# Principio de Parsimonia y Criterio de Información de Akaike
# Fórmulas incluidas en el código
# ============================================================

# Cargar librería para ajuste por máxima verosimilitud (MLE)
library(MASS)

# Datos de ejemplo: caudales máximos anuales (m³/s)
set.seed(42)  # para reproducibilidad
caudales <- rlnorm(50, meanlog = 3, sdlog = 0.5)  # simulados de Lognormal

# ---- Ajuste de distribuciones por MLE ----

# Ajuste Gamma (parámetros: shape, rate) – k = 2 parámetros
fit_gamma <- fitdistr(caudales, "gamma")
logLik_gamma <- fit_gamma$loglik  # log-verosimilitud máxima ( ln(L) )

# Ajuste Lognormal (parámetros: meanlog, sdlog) – k = 2 parámetros
fit_lnorm <- fitdistr(caudales, "lognormal")
logLik_lnorm <- fit_lnorm$loglik

# ---- Cálculo del AIC con la fórmula ----
# AIC = -2 * log(L) + 2 * k
k <- 2  # ambas distribuciones tienen 2 parámetros

AIC_gamma <- -2 * logLik_gamma + 2 * k
AIC_lnorm <- -2 * logLik_lnorm + 2 * k

cat("AIC (Gamma)    =", AIC_gamma, "\n")

## AIC (Gamma)    = 387.6349

cat("AIC (Lognormal) =", AIC_lnorm, "\n")

## AIC (Lognormal) = 387.8899

# ---- Regla de oro: ΔAIC > 2 ----
# ΔAIC = AIC(modelo más alto) - AIC(modelo más bajo)
delta_AIC <- abs(AIC_gamma - AIC_lnorm)  # valor absoluto
mejor_modelo <- ifelse(AIC_gamma < AIC_lnorm, "Gamma", "Lognormal")

cat("\nΔAIC =", delta_AIC, "\n")

## 
## ΔAIC = 0.2550253

cat("Gana el modelo:", mejor_modelo, "\n")

## Gana el modelo: Gamma

if (delta_AIC > 2) {
  cat("→ Evidencia fuerte: el modelo con menor AIC es claramente superior.\n")
} else {
  cat("→ ΔAIC < 2: no hay diferencia significativa; ambos son igualmente plausibles.\n")
}

## → ΔAIC < 2: no hay diferencia significativa; ambos son igualmente plausibles.

# ---- Nota: fórmula completa de AIC en R ----
# También puedes calcularlo manualmente (sin fitdistr):
# logLik_manual <- sum(dgamma(caudales, shape = fit_gamma$estimate[1], 
#                              rate = fit_gamma$estimate[2], log = TRUE))
# AIC_manual <- -2 * logLik_manual + 2 * k