PARSIMONIA_AIC

Estadistica

Identificación del estudiante

Estudiante: Luis Perez (1005486307)

URL en RPubs:https://rpubs.com/LuisPerez/parsimonia_aic

Prompt 1: El problema y la Verosimilitud

Prompt 2: El desempate (Parsimonia y AIC)

La Función de Verosimilitud (Likelihood) en la Inferencia Estadística

En la ingeniería de recursos hídricos y el análisis de frecuencias de caudales extremos, la función de verosimilitud (\(\mathcal{L}\)) es el motor matemático central para la estimación de parámetros y la selección de modelos probables.

1. Probabilidad vs Verosimilitud

La diferencia fundamental entre ambos conceptos radica en qué variables son fijas y cuáles son variables:\[\text{Probabilidad: } P(X = x \mid \theta) \quad \text{vs.} \quad \text{Verosimilitud: } \mathcal{L}(\theta \mid X = x)\]

Probabilidad:

Los parámetros del modelo (\(\theta\)) se consideran conocidos y fijos. Se evalúa la probabilidad de observar diferentes conjuntos de datos futuros (\(x\)). La suma/integración sobre el espacio muestral es igual a 1:\[\int P(x \mid \theta) \, dx = 1\]

Verosimilitud:

El conjunto de datos observados (\(X = x\)) es fijo y conocido (ej. tu serie histórica de caudales máximos anuales \(x_1, x_2, \dots, x_n\)). La función evalúa qué tan plausibles son distintos valores de los parámetros (\(\theta\)).

Nota teórica clave:

\(\mathcal{L}(\theta \mid x)\) no es una densidad de probabilidad sobre \(\theta\). Por ende, la superficie bajo la curva de verosimilitud no integra a 1:\[\int \mathcal{L}(\theta \mid x) \, d\theta \neq 1\]

2. Construcción Matemática de la Función

Para una muestra aleatoria independiente e idénticamente distribuida (i.i.d.) de \(n\) mediciones de caudales máximos observados (\(x_1, x_2, \dots, x_n\)) con una función de densidad de probabilidad (PDF) \(f(x \mid \theta)\):\[\mathcal{L}(\theta \mid x_1, x_2, \dots, x_n) = \prod_{i=1}^{n} f(x_i \mid \theta)\]

La Log-Verosimilitud (Log-Likelihood)

Debido a que el producto de probabilidades genera números extremadamente pequeños que causan desbordamiento de flujo bajo (underflow) en cálculos computacionales, se aplica la transformación monotónica del logaritmo natural (\(\ln\)):\[\ell(\theta) = \ln \mathcal{L}(\theta \mid x) = \sum_{i=1}^{n} \ln f(x_i \mid \theta)\]Dado que el logaritmo es una función estrictamente creciente, el valor del parámetro \(\theta\) que maximiza \(\ell(\theta)\) es exactamente el mismo que maximiza \(\mathcal{L}(\theta)\).

3. Estimación de Máxima Verosimilitud (MLE)

El estimador de máxima verosimilitud (\(\hat{\theta}_{\text{MLE}}\)) se define formalmente como:\[\hat{\theta}_{\text{MLE}} = \arg\max_{\theta \in \Theta} \ell(\theta \mid x)\]

Proceso de optimización analítica

Para calcular los parámetros óptimos de una distribución (por ejemplo, los parámetros \(\alpha\) y \(\beta\) de una distribución Gamma, o \(\mu\) y \(\sigma\) de una Lognormal):

1. Derivación de las Ecuaciones de Score:

Se calcula el gradiente del logaritmo de verosimilitud respecto a los parámetros:\[S(\theta) = \nabla_\theta \, \ell(\theta) = 0\]

2. Comprobación de la Matriz Hessiana:

Se evalúa la segunda derivada para garantizar un máximo local (matriz definida negativa):\[\mathbf{H}(\theta) = \nabla_\theta^2 \, \ell(\theta) < 0\]

3. Inversión Numérica::

En distribuciones complejas de hidrología (como Gumbel, Log-Pearson Type III o GEV), el sistema de ecuaciones analítico no tiene solución cerrada y se resuelve numéricamente mediante métodos iterativos como Newton-Raphson o Nelder-Mead.

4. Métricas Hidrológicas Derivadas de la Verosimilitud

# ==============================================================================
# MÉTRICAS HIDROLÓGICAS DERIVADAS DE LA VEROSIMILITUD EN R
# ==============================================================================

# 1. Datos de ejemplo: Caudales máximos anuales (m³/s)
caudales <- c(120, 145, 180, 210, 225, 260, 290, 310, 350, 410, 480, 520)

# Cargar librería para ajuste de distribuciones por MLE
if (!require("fitdistrplus")) install.packages("fitdistrplus")

## Loading required package: fitdistrplus

## Loading required package: MASS

## Loading required package: survival

library(fitdistrplus)

# ------------------------------------------------------------------------------
# MÉTRICA 1 & 2: INFORMACIÓN DE FISHER Y MATRIZ DE COVARIANZA
# ------------------------------------------------------------------------------
# Ajustamos una distribución Gamma por Máxima Verosimilitud
fit_gamma <- fitdist(caudales, "gamma", method = "mle")

# A) Matriz de Covarianza: Var(θ_hat) ≈ I(θ_hat)^(-1)
# En R se extrae directamente con vcov()
matriz_covarianza <- vcov(fit_gamma)

# B) Matriz de Información de Fisher observada: I(θ_hat)
# Se obtiene invirtiendo la matriz de covarianza (Hessiana observada)
informacion_fisher <- solve(matriz_covarianza)

# C) Error Estándar e Intervalos de Confianza al 95% para los parámetros (shape α, rate β)
errores_estandar <- sqrt(diag(matriz_covarianza))
parametros       <- fit_gamma$estimate
ic_inferior      <- parametros - 1.96 * errores_estandar
ic_superior      <- parametros + 1.96 * errores_estandar

tabla_parametros <- data.frame(
  Estimacion = parametros,
  SE         = errores_estandar,
  IC_95_Inf  = ic_inferior,
  IC_95_Sup  = ic_superior
)

cat("\n--- 1. INFORMACIÓN DE FISHER Y MATRIZ DE COVARIANZA ---\n")

## 
## --- 1. INFORMACIÓN DE FISHER Y MATRIZ DE COVARIANZA ---

cat("\nMatriz de Información de Fisher Observada I(θ):\n")

## 
## Matriz de Información de Fisher Observada I(θ):

print(informacion_fisher)

##             shape        rate
## shape    2.378835   -633.7595
## rate  -633.759550 185722.6044

cat("\nMatriz de Covarianza Var(θ_hat):\n")

## 
## Matriz de Covarianza Var(θ_hat):

print(matriz_covarianza)

##            shape         rate
## shape 4.62546398 1.578393e-02
## rate  0.01578393 5.924542e-05

cat("\nParametros Estimados con Intervalos de Confianza (95%):\n")

## 
## Parametros Estimados con Intervalos de Confianza (95%):

print(tabla_parametros)

##       Estimacion          SE   IC_95_Inf  IC_95_Sup
## shape 5.52808252 2.150689188 1.312731709 9.74343332
## rate  0.01895223 0.007697105 0.003865902 0.03403855

# ------------------------------------------------------------------------------
# MÉTRICA 3: PRUEBA DE RAZÓN DE VEROSIMILITUD (LRT)
# ------------------------------------------------------------------------------
# Compara dos modelos anidados: Exponencial (k=1, simplificado) vs Gamma (k=2, completo)
fit_exp <- fitdist(caudales, "exp", method = "mle")

# Extraer Log-Verosimilitudes (l)
loglik_exp   <- fit_exp$loglik     # Modelo nulo (θ_0)
loglik_gamma <- fit_gamma$loglik   # Modelo alternativo (θ)

# Estadística de la Razón de Verosimilitud: Λ = -2 * (l_nulo - l_completo)
lambda_LRT <- -2 * (loglik_exp - loglik_gamma)

# Grados de libertad (diferencia en número de parámetros)
df <- length(fit_gamma$estimate) - length(fit_exp$estimate) # 2 - 1 = 1

# valor p según la distribución Chi-cuadrado (χ²)
p_value <- pchisq(lambda_LRT, df = df, lower.tail = FALSE)

cat("\n--- 2. PRUEBA DE RAZÓN DE VEROSIMILITUD (LRT) ---\n")

## 
## --- 2. PRUEBA DE RAZÓN DE VEROSIMILITUD (LRT) ---

cat(sprintf("Log-Likelihood Exponencial (k=1) : %.4f\n", loglik_exp))

## Log-Likelihood Exponencial (k=1) : -80.1073

cat(sprintf("Log-Likelihood Gamma (k=2)       : %.4f\n", loglik_gamma))

## Log-Likelihood Gamma (k=2)       : -74.1191

cat(sprintf("Estadístico Λ (Lambda)           : %.4f\n", lambda_LRT))

## Estadístico Λ (Lambda)           : 11.9764

cat(sprintf("p-valor (Chi-cuadrado df=%d)      : %.6f\n", df, p_value))

## p-valor (Chi-cuadrado df=1)      : 0.000539

if (p_value < 0.05) {
  cat("Conclusión: Se rechaza el modelo Exponencial. La distribución Gamma justifica la adición del segundo parámetro con significancia estadística.\n")
} else {
  cat("Conclusión: No hay evidencia suficiente para preferir la Gamma sobre la Exponencial.\n")
}

## Conclusión: Se rechaza el modelo Exponencial. La distribución Gamma justifica la adición del segundo parámetro con significancia estadística.

# ------------------------------------------------------------------------------
# MÉTRICA 4: CRITERIO DE INFORMACIÓN DE AKAIKE (AIC)
# ------------------------------------------------------------------------------
# Ajustar modelos alternativos no anidados (Lognormal vs Gamma)
fit_ln <- fitdist(caudales, "lnorm", method = "mle")

# Cálculo manual para verificar la fórmula: AIC = 2k - 2*l(θ)
k_ln     <- length(fit_ln$estimate)
l_ln     <- fit_ln$loglik
aic_calc <- (2 * k_ln) - (2 * l_ln)

# Extracción automática desde R
tabla_aic <- data.frame(
  Modelo     = c("Lognormal", "Gamma", "Exponencial"),
  k          = c(length(fit_ln$estimate), length(fit_gamma$estimate), length(fit_exp$estimate)),
  LogLik     = c(fit_ln$loglik, fit_gamma$loglik, fit_exp$loglik),
  AIC_R      = c(fit_ln$aic, fit_gamma$aic, fit_exp$aic),
  AIC_Manual = c(aic_calc, (2*2 - 2*loglik_gamma), (2*1 - 2*loglik_exp))
)

# Delta AIC respecto al mejor modelo
tabla_aic$Delta_AIC <- tabla_aic$AIC_R - min(tabla_aic$AIC_R)
tabla_aic           <- tabla_aic[order(tabla_aic$AIC_R), ]

cat("\n--- 3. CRITERIO DE INFORMACIÓN DE AKAIKE (AIC) ---\n")

## 
## --- 3. CRITERIO DE INFORMACIÓN DE AKAIKE (AIC) ---

print(tabla_aic)

##        Modelo k    LogLik    AIC_R AIC_Manual Delta_AIC
## 2       Gamma 2 -74.11914 152.2383   152.2383 0.0000000
## 1   Lognormal 2 -74.18202 152.3640   152.3640 0.1257559
## 3 Exponencial 1 -80.10734 162.2147   162.2147 9.9763910

# Dividir el área de dibujo en 2 filas y 2 columnas (4 gráficos)
par(mfrow = c(2, 2))

# Generar los 4 gráficos
plot(1:10, main = "Gráfico 1")
hist(rnorm(100), main = "Gráfico 2", col = "skyblue")
boxplot(rnorm(50), main = "Gráfico 3", col = "lightgreen")
plot(density(rnorm(100)), main = "Gráfico 4", col = "coral")

# Restablecer la vista a 1 solo gráfico
par(mfrow = c(1, 1))

# Instalar y cargar librerías
if (!require("ggplot2")) install.packages("ggplot2")

## Loading required package: ggplot2

if (!require("patchwork")) install.packages("patchwork")

## Loading required package: patchwork

## 
## Attaching package: 'patchwork'

## The following object is masked from 'package:MASS':
## 
##     area

library(ggplot2)
library(patchwork)

# Crear gráficos individuales
g1 <- ggplot(mtcars, aes(x = wt, y = mpg)) + geom_point() + ggtitle("Gráfico A")
g2 <- ggplot(mtcars, aes(x = hp)) + geom_histogram(bins = 10) + ggtitle("Gráfico B")
g3 <- ggplot(mtcars, aes(x = factor(cyl), y = mpg)) + geom_boxplot() + ggtitle("Gráfico C")

# Opción A: Dos al lado y uno abajo
(g1 + g2) / g3

# Opción B: Tres en fila horizontal con título general
g1 + g2 + g3 + plot_annotation(title = "Panel de Evaluación Hidrológica")

# **5. Ejemplo Aplicado: Caudales Máximos Anuales** 
Supón que tienes una muestra de caudales $X = \{120, 180, 210, 350\} \text{ m}^3/\text{s}$.

1. Se asume que los datos siguen una distribución Lognormal con parámetros $\theta_{\text{LN}} = (\mu, \sigma)$ y una distribución Gamma con $\theta_{\text{G}} = (\alpha, \beta)$.

2. Se evalúa la PDF de cada distribución para los $n$ puntos observados.

3. Se calcula el acumulado multiplicativo (o la suma logarítmica) para ambas distribuciones en sus valores óptimos:$$\ell(\hat{\theta}_{\text{LN}}) = -142.5 \quad \text{vs.} \quad \ell(\hat{\theta}_{\text{G}}) = -148.1$$

Puesto que $-142.5 > -148.1$, la distribución Lognormal presenta una mayor verosimilitud ($\mathcal{L}_{\text{LN}} > \mathcal{L}_{\text{G}}$), lo que indica matemáticamente que el proceso lognormal tiene una mayor plausibilidad de haber generado la muestra de caudales máximos analizada.


``` r
# ------------------------------------------------------------------------------
# 1. DATOS HISTÓRICOS DE CAUDALES MÁXIMOS ANUALES (m³/s)
# ------------------------------------------------------------------------------
# Muestra sintética de caudales para el ejemplo
caudales <- c(120, 145, 180, 210, 225, 260, 290, 310, 350, 410, 480, 520)

# ------------------------------------------------------------------------------
# 2. INSTALACIÓN Y CARGA DE LIBRERÍAS
# ------------------------------------------------------------------------------
if (!require("fitdistrplus")) install.packages("fitdistrplus")
library(fitdistrplus)

# ------------------------------------------------------------------------------
# 3. ESTIMACIÓN DE MÁXIMA VEROSIMILITUD (MLE)
# ------------------------------------------------------------------------------
# Ajuste de la distribución Lognormal
fit_ln <- fitdist(caudales, "lnorm", method = "mle")

# Ajuste de la distribución Gamma
fit_gamma <- fitdist(caudales, "gamma", method = "mle")

# Ajuste de la distribución Gumbel (requiere la librería actuar)
if (!require("actuar")) install.packages("actuar")

## Loading required package: actuar

## 
## Attaching package: 'actuar'

## The following objects are masked from 'package:stats':
## 
##     sd, var

## The following object is masked from 'package:grDevices':
## 
##     cm

library(actuar)
fit_gumbel <- fitdist(caudales, "gumbel", method = "mle")

# ------------------------------------------------------------------------------
# 4. EXTRACCIÓN DE LOG-VEROSIMILITUD (logLik) Y AIC
# ------------------------------------------------------------------------------
# Log-Verosimilitud (l)
loglik_ln     <- fit_ln$loglik
loglik_gamma  <- fit_gamma$loglik
loglik_gumbel <- fit_gumbel$loglik

# Criterio de Información de Akaike (AIC)
aic_ln     <- fit_ln$aic
aic_gamma  <- fit_gamma$aic
aic_gumbel <- fit_gumbel$aic

# ------------------------------------------------------------------------------
# 5. TABLA COMPARATIVA Y SELECCIÓN DEL MODELO GANADOR
# ------------------------------------------------------------------------------
resultados <- data.frame(
  Distribucion     = c("Lognormal", "Gamma", "Gumbel"),
  Parametros_k     = c(length(fit_ln$estimate), length(fit_gamma$estimate), length(fit_gumbel$estimate)),
  Log_Likelihood   = c(loglik_ln, loglik_gamma, loglik_gumbel),
  AIC              = c(aic_ln, aic_gamma, aic_gumbel)
)

# Calcular la diferencia de AIC respecto al mínimo (Delta AIC)
resultados$Delta_AIC <- resultados$AIC - min(resultados$AIC)

# Ordenar por el menor AIC (Regla de Oro)
resultados <- resultados[order(resultados$AIC), ]

print("--- TABLA DE SELECCIÓN DE MODELOS ---")

## [1] "--- TABLA DE SELECCIÓN DE MODELOS ---"

print(resultados)

##   Distribucion Parametros_k Log_Likelihood      AIC Delta_AIC
## 2        Gamma            2      -74.11914 152.2383 0.0000000
## 1    Lognormal            2      -74.18202 152.3640 0.1257559
## 3       Gumbel            2      -74.26323 152.5265 0.2881771

# ------------------------------------------------------------------------------
# 6. EVALUACIÓN VISUAL (Gráficos Q-Q y Densidad)
# ------------------------------------------------------------------------------
par(mfrow = c(1, 2)) # Dividir el área de gráfico en 2 columnas

# Gráfico de Densidades Ajustadas
denscomp(list(fit_ln, fit_gamma, fit_gumbel), 
         legendtext = c("Lognormal", "Gamma", "Gumbel"),
         main = "Ajuste de Densidad a Caudales",
         xlab = "Caudal (m³/s)", ylab = "Densidad")

# Gráfico Q-Q (Cuantil-Cuantil)
qqcomp(list(fit_ln, fit_gamma, fit_gumbel), 
       legendtext = c("Lognormal", "Gamma", "Gumbel"),
       main = "Gráfico Q-Q para Desempate")

par(mfrow = c(1, 1)) # Restablecer vista de gráficos