Capítulo 3 Ajuste de distribuciones y pruebas de Bondad

3.1 Objetivo de la práctica

Aplicar el análisis de frecuencia hidrológica a una serie de precipitaciones máximas, ajustando los datos a una distribución de probabilidad, para estimar la lluvia esperada en diferentes periodos de retorno y comprobar, con pruebas estadísticas, si el ajuste realizado describe adecuadamente el comportamiento de los datos

3.2 Resultados de la práctica

Aplicar un análisis de frecuencia a datos hidrológicos, realizar el ajuste a distribuciones de probabilidad y evaluar la validez estadística de dicho ajuste

3.3 Instalación de las librerías necesarias

Para este clase necesitaremos 4 librerías que si aún no las has usado es importante que las instales previamente.

librería readxl: Lectura de datos tipo xlsx
librería texmex: Ajuste de distribuciones de probabilidad
librería fitdistrplus: Ajuste de distribuciones, evaluación y diagnóstico de ajustes
librería ggplot2: Producción de gráficos

library("readxl")
library("texmex")
library("fitdistrplus")
library("ggplot2")

3.4 ¿Que datos usaremos en la clase?

Para esta práctica se utilizarán los registros de precipitaciones máximas, con una duración de 1440 minutos, de la estación M0003–Izobamba (Código P16).

(#fig:Ubicación_Izobamba)Ubicación de la Estación Izomaba, sector de Cutuglahua, al sur de Quito

(#fig:Información_Izobamba)Información Estación Izobamba

#Carguemos y almacenemos el archivo

max_preci<-read_excel("./Documentos/Precipitacion_maxima_IZOBAMBA.xlsx",
                      skip = 1 )
head(max_preci)
#> # A tibble: 6 × 10
#>    Year    d5   d10   d15   d20   d30   d60  d120  d360
#>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1  1962   9    14.8  12.8  13.5  14.1  16.8  20.4  27  
#> 2  1963   9.2  10.5  12.3  14.3  17.5  23.6  25.6  28.8
#> 3  1964   9.5  12.2  14.7  15.6  20.1  24.3  33.2  37.8
#> 4  1965   8.8  13    16.2  17    20.4  35    39.4  44.4
#> 5  1966   9    11.1  12.9  15    18.5  19.3  22.2  25.8
#> 6  1967   6.5   7.2   8.5  10    14    20.8  24.8  28.2
#> # ℹ 1 more variable: d1440 <dbl>

3.5 Empecemos

El archivo contiene preciptaciones máximas para diferentes duraciones, sin embargo en esta clase evaluaremos únicamente los datos de precipitación máxima con duración de 1440 minutos, por lo que crearemos un vector unicamente con esos datos

preci_dmin <- data.frame(d1440 = max_preci$d1440)

# Creemos un Histograma Sencillo
hist(preci_dmin$d1440,main="", xlab = "Precipitación duración 1440m")

Como puedes observar en el histograma los valores de precipitación se encuentran sesgados, por lo ajustaremos a una función de distribución que permita el análisis de valores extremos. Durante la clase trabajaremos con la distribución de Gumbel (distribución generalizada de valores extremos).

3.6 Ajuste a distribución teórica

Las funciones de distribución de probabilidad se relacionan directamente con la estimación de probabilidades de ocurrencia de eventos hidrológicos. En este contexto, el primer paso consiste en asociar a cada dato observado una probabilidad empírica de excedencia o de no ocurrencia (probabilidad de que un evento no ocurra en un año dado), la cual se vincula con el concepto de periodo de retorno. El periodo de retorno, se define como el intervalo de tiempo promedio entre eventos que igualan o superan una determinada magnitud, y se calcula generalmente como: \(Tr = \\frac{n+1}{m}\) donde n es el número total de datos y m el orden del dato en la serie ordenada

#Ordenemos los datos de mayor a menor
preci_dmin$d1440 <- preci_dmin[order(preci_dmin$d1440, decreasing = TRUE), ]

#Obtengamos el número de datos y una secuencia del 1 a n
#Esto será m el orden de dato
n <- nrow(preci_dmin)
i <- seq_len(n)

#Calculemos el periodo de retorno y la probabilidad
#Observa que la porbabilidad es de no ocurrencia 
preci_dmin$Tr <- round((n + 1) / i,2)
preci_dmin$p <- round(1 - 1 / preci_dmin$Tr,2)

#Analicemos si estamos en lo correcto
head(preci_dmin,10)
#>    d1440    Tr    p
#> 1   82.8 46.00 0.98
#> 2   70.4 23.00 0.96
#> 3   67.2 15.33 0.93
#> 4   57.6 11.50 0.91
#> 5   57.6  9.20 0.89
#> 6   54.7  7.67 0.87
#> 7   50.8  6.57 0.85
#> 8   50.4  5.75 0.83
#> 9   48.3  5.11 0.80
#> 10  48.0  4.60 0.78

Una vez que tenemos que hemos asociado la precipitación a una porbabilidad procederemos a ajustar, la función que ocuparemos en R realiza todo el proceso de ajuste en un paso pero es importante que sepas que hace la función detrás de tus ojos.

Para recordar: Ajustar datos a una función de distribución consiste en encontrar el modelo probabilístico que mejor representa el comportamiento estadístico de un conjunto de observaciones. En otras palabras, buscamos una distribución teórica (como Normal, Log-Normal, Gumbel, etc.) cuyos parámetros permitan describir la variabilidad de los datos observados y realizar inferencias, como la estimación de eventos extremos o el cálculo de probabilidades.

El ajuste permite transformar datos empíricos en un modelo matemático generalizable, que facilita la predicción

(#fig:Diagrama de Flujo)Pasos para ajustar a una distribución Gumble

La función para realizar el ajuste se llama fitdist() y pertenece al paquete fitdistrplus, y trabaja en conjunto con el paquete texmex, para obtener el ajuste son necesarias cuatro variables, la primera los valores a ajustar, en este caso los valores de precipitación máxima con duración de 1440 minutos, el segundo el nombre de la distribución a ajustar, el tercero el media y finalmente la desviación estadar. Es importante mencionar que distribución a distribución puede cambiar un poco la menra de ingresar los datos por lo que les dejo el link del manualde fitdistrplus

G <- fitdist(preci_dmin$d1440, "gumbel", 
                         start=list(mu=mean(preci_dmin$d1440), sigma=sd(preci_dmin$d1440)))
plot(G)

3.7 Pruebas de bondad

Las pruebas de bondad de ajuste son herramientas estadísticas que permiten evaluar si un conjunto de datos observados puede considerarse proveniente de una distribución de probabilidad teórica específica. se plantea una hipótesis nula (Ho), que establece que los datos siguen la distribución propuesta.

El procedimiento consiste en comparar la discrepancia entre la distribución empírica (la que proviene directamente de los datos) y la distribución teórica ajustada. Si esa discrepancia es pequeña y no significativa estadísticamente, se acepta Ho.

Por lo tanto evaluar la validez de los de los modelos de distribución ajustados a los datos es fundamental. En R se emplea con frecuencia el paquete fitdistrplus, que provee herramientas robustas para el ajuste y evaluación de distribuciones de probabilidad. Dentro de este paquete, la función gofstat() permite obtener resultados de las pruebas de bondad de ajuste. Esta función trabaja a partir de un objeto previamente ajustado (por ejemplo, con fitdist, que fue el paso anterior), y calcula de forma automática estadísticos como Kolmogorov–Smirnov, Anderson–Darling, entre otros.

PB_G<- gofstat(G)

#Impresión de resultados de prueba de Kolmogorov–Smirnov
PB_G$ks
#> 1-mle-gumbel 
#>    0.1639417

#Impresion de resumen de resultados completos 
gofstat(G,fitnames = "Gumbel")
#> Goodness-of-fit statistics
#>                                 Gumbel
#> Kolmogorov-Smirnov statistic 0.1639417
#> Cramer-von Mises statistic   0.1485403
#> Anderson-Darling statistic   0.7839431
#> 
#> Goodness-of-fit criteria
#>                                  Gumbel
#> Akaike's Information Criterion 338.5845
#> Bayesian Information Criterion 342.1978

3.8 Cálculo de valores asociados a periodos de retorno

Muy bien, una vez ajustado el modelo y verificada su validez estadística, el siguiente paso consiste en estimar los valores de la variable de interés asociados a distintos periodos de retorno. Para ello, se emplea la función de distribución ajustada, la cual permite calcular cuantiles que representan la magnitud esperada del fenómeno con una determinada probabilidad de excedencia

preci_dmin$pc<- round(qgumbel(preci_dmin$p,mu = mean(preci_dmin$d1440),sigma = sd(preci_dmin$d1440)),2)

head(preci_dmin,10)
#>    d1440    Tr    p    pc
#> 1   82.8 46.00 0.98 87.19
#> 2   70.4 23.00 0.96 79.01
#> 3   67.2 15.33 0.93 72.31
#> 4   57.6 11.50 0.91 69.26
#> 5   57.6  9.20 0.89 66.80
#> 6   54.7  7.67 0.87 64.72
#> 7   50.8  6.57 0.85 62.92
#> 8   50.4  5.75 0.83 61.33
#> 9   48.3  5.11 0.80 59.23
#> 10  48.0  4.60 0.78 57.98

Finalmente, se elaborará un gráfico comparativo entre los valores observados y los estimados en función del periodo de retorno. Para ello se empleará el paquete ggplot2, donde el eje X representará el periodo de retorno y el eje Y la precipitación máxima con una duración de 1440 minutos; los valores observados se representarán mediante puntos, mientras que los valores ajustados se mostrarán mediante una línea.

ggplot(preci_dmin, aes(x = Tr)) + 
  geom_point(aes(y = d1440, colour = "Observados")) +
  geom_line(aes(y = pc, colour = "Calculados")) +
  scale_colour_manual(values = c("Observados" = "#22577a", "Calculados" = "#57cc99")) +
  xlab("Periodo de retorno") + 
  ylab("Precipitación")

3.9 📝 TRABAJO AUTÓNOMO

En esta prueba práctica se determinará la distribución de probabilidad con mejor ajuste para una serie de datos de precipitación máxima a través de pruebas de bondad. El taller se entregará de manera individual y se utilizarán el archivo Datos_Taller_EPMAPS_RUMIPAMBA, que se encuentra en la carpeta Documentos de repositorio de Github, obteniendo los siguientes resultados:

Parte 1: Obtención de Resultados

A partir de las series de datos de precipitación máxima, se deberá calcular para 5, 30 y 1440 min de precipitación lo siguiente: Periodo de retorno, probabilidad de no ocurrencia y precipitación calculada. Para esto se deberá usar las funciones de distribución de Gumbel, Normal y Lognormal (20%)
Gráficos de comparación de densidad (denscomp) y de cuantiles (qqcomp) para las tres distribuciones y cada duración. Escriba un pequeño comentario respecto a lo visualizado en los gráficos (10%)

Parte 2: Análisis

Aplique la prueba de Kolmogorov-Smirnov (K-S) para cada duración y distribución y determine para cada duración cuál de las distribuciones de probabilidad tiene un mejor ajuste. Justifique su respuesta, incluir en la justificación si el análisis visual de los gráficos está acorde a la prueba de bondad aceptada (40%)

Parte 3: Aplicación Simple

¿Cuáles la probabilidad de que en un año cualquiera la precipitación con duración de 5 minutos sea igual o mayor a 8 mm? (De una respuesta por cada distribución de probabilidad) (15%)
Para un periodo de retorno de 50 años ¿Cuál será la precipitación máxima anual de duración 5 min? (15%)

El taller deberá ser entregado el día y hora acordada en formato pdf, a través del aula virtual.

2 Introducción a las distribuciones de probabilidad en R

4 Curvas de Intensidad, Duración y Frecuencia