Distribución GVE bivariada (GVEb)

Emil Julius Gumbel estableció a inicios de los años sesenta el llamado modelo bivariado logístico, el cual acepta como distribuciones marginales las de valores extremos: Gumbel, General de Valores Extremos (GVE), Gumbel mixta y TCEV o de valores extremos de dos componentes (^{Escalante & Raynal, 1994}; ^{Ramírez & Aldama, 2000}; ^{Escalante & Reyes, 2004}), su ecuación es:

Donde:

F(x) y F(y) = funciones de distribución de probabilidades (FDP) marginales de las variables aleatorias X y Y.

m = parámetro de asociación, el cual depende de la correlación entre las variables.

Si las marginales son distribuciones GVE o General de Valores Extremos (^{Hosking & Wallis, 1997}; ^{Rao & Hamed, 2000}; ^{Stedinger, 2017}), sus expresiones son:

u, α y k son los parámetros de ubicación, escala y forma. Sustituyendo las ecuaciones 2 y 3 en la 1, se obtiene la FDP de la GVEb (^{Escalante & Raynal, 1994}):

con θ=u1,α1,k1,u2,α2,k2,m ; son los siete parámetros de ajuste y se obtienen con el método estadístico de máxima verosimilitud.

La distribución GVE representa tres modelos probabilísticos: si k < 0 la GVE tiene el tipo Fréchet, sin límite superior; caso aquí estudiado; si k = 0 se tiene el modelo Gumbel y cuando k > 0 se tiene la Weibull, con límite superior. La distribución GVE está definida por el conjunto {x_i: [1 - k(x_i - u)/α] > 0}, lo cual se indica con el signo + afuera del paréntesis.

^{Coles (2001)} destaca que cualquier combinación de parámetros de ajuste que viola la condición anterior de positividad, implica que al menos uno de los puntos observados (x_i), está más allá de los puntos finales de la distribución y entonces la función de verosimilitud es cero y su versión logarítmica -∞.

La solución inversa de las ecuaciones 2 y 3, permiten la estimación de predicciones (x_p, y_p) asociadas a una probabilidad de no excedencia p = F(x) ó p = F(y), son las siguientes:

En este estudio se ajustó la distribución GVE con tres métodos: Sextiles, momentos L y máxima verosimilitud, los cuales se pueden consultar respectivamente en: ^{Campos (2006)}, ^{Hosking & Wallis (1997)} y ^{Rao & Hamed (2000)}.

Restricciones de probabilidad

Las probabilidades de no excedencia univariadas y conjunta de la distribución GVEb, deben cumplir con la restricción siguiente (^{Escalante & Raynal, 1994}; ^{Ramírez & Aldama, 2000}; ^{Escalante, 2004}):

Método de máxima verosimilitud

Para una muestra aleatoria (X₁, X₂, . . . , X_n) de observaciones independientes e idénticamente distribuidas (iid) que siguen una FDP denominada F_θ con parámetros de ajuste θ₁, θ₂, . . . ,θ_q. La probabilidad de obtener un valor X_i, será (^{Rao & Hamed, 2000}; ^{Coles, 2001}; ^{Meylan et al., 2012}):

siendo, f_θ(x_i) la función de densidad de probabilidad. Como los datos son iid, la probabilidad de obtener los n valores X_i, será la probabilidad conjunta o función de verosimilitud, designada L del inglés likelihood; cuya ecuación es:

El método de máxima verosimilitud consiste en encontrar un vector θ^ de parámetros que hagan máxima a L(θ) y por lo tanto, a la probabilidad de obtener la muestra (X₁, X₂, . . . ,X_n ). Con frecuencia, resulta más conveniente tomar logaritmos y trabajar con la función logarítmica de verosimilitud (^{Coles, 2001}), es decir:

Lo anterior, es aceptable debido a que la función logarítmica es monotónica y entonces la función l(θ) alcanza su máximo en el mismo punto que la función L(θ).

l(θ) de la distribución GVEb

La función logarítmica de verosimilitud de la GVEb para el caso simple de variables aleatorias X, Y con igual amplitud de registro (n), ha sido expuesta por ^{Escalante & Raynal (1994)} y es la siguiente:

Periodos de retorno univariados

La probabilidad de un evento se define como el cociente del número de casos favorables (ncf) entre el número de casos posibles (ncp) a dicho evento y varía de cero a uno. Debido al manejo anual de las variables X, Y, la probabilidad de excedencia F'(x) es igual al recíproco del periodo de retorno (T_X) en años; pues en cada año se tiene, ncf = 1 y ncp = T_X, entonces (^{Yue & Rasmussen, 2002}; ^{Shiau, 2003}):

En las expresiones anteriores, F(x) y F(y) se estiman con las ecuaciones 2 y 3.

Periodos de retorno conjuntos

El primer periodo de retorno conjunto del evento (X, Y) se define bajo la condición OR y la ecuación 14 de la manera siguiente (^{Goel et al., 1998}; ^{Yue, 2000b}; ^{Shiau, 2003}):

F(x,y) es la probabilidad de no excedencia conjunta que se estima con la ecuación 4, previa estimación de sus parámetros de ajuste óptimos u1,α1,k1,u2,α2,k2,m con el método de máxima verosimilitud. Este evento representa el caso en que los límites x ó y, o ambos pueden ser excedidos.

El segundo periodo de retorno conjunto del evento (X, Y) está asociado al caso en que ambos límites son excedidos (X > x, Y > y) o condición AND, su ecuación es (^{Goel et al., 1998}; ^{Aldama, 2000}; ^{Ramírez & Aldama, 2000}; ^{Yue, 2000b} y ^{Shiau, 2003}):

^{Aldama (2000)} obtiene la expresión F’(x,y) de la probabilidad conjunta de excedencia mediante un razonamiento de probabilidades lógico y simple aplicado en el plano cartesiano. En cambio, ^{Yue y Rasmussen (2002)} recurren al plano cartesiano para definir numéricamente un evento bivariado (X, Y), que puede ocurrir en alguno de los cuatro cuadrantes.

Diversos autores (^{Yue, 2000b}; ^{Yue & Rasmussen, 2002}; ^{Shiau, 2003}) han mostrado las gráficas de los dos periodos de retorno conjuntos y han discutido sus diferencias. En la Figura 1 se muestra la gráfica del periodo de retorno conjunto T’(Q,V), construida con los datos de la aplicación numérica que será expuesta.

Figura 1 Gráficas de los cuatro periodos de retorno conjunto T’(Q,V) de diseño, de las crecientes de entrada a la presa Venustiano Carranza (Don Martín), México 

De acuerdo con ^{Yue et al. (1999)}, ^{Yue (2000b)} y ^{Yue & Rasmussen (2002)} existe un tercer tipo de periodo de retorno conjunto, que tiene aplicación en la práctica hidrológica y que se define para un evento X dado que Y ≤ y o para un evento Y dado que X ≤ x y por ello, se designan condicionales. Para tales eventos, sus distribuciones de probabilidad condicional se definen de manera simple, con estas ecuaciones:

Sustituyendo las expresiones anteriores en la ecuación 14, se obtienen las fórmulas del periodo de retorno conjunto condicional:

Eventos críticos del T’(x,y)

^{Volpi y Fiori (2012)} destacan que la gráfica del periodo de retorno conjunto de tipo AND, mostrada como Figura 1, presenta una severa inconsistencia al contener, en un contexto bivariado, umbrales críticos univariados. Debido a lo anterior, tal gráfica se considera integrada por dos porciones, las dos designadas simples (naive part) y la correcta (proper part). Las partes rectas son las colas o rectas asíntotas a la parte curva. La probabilidad de ocurrencia de un evento o pareja de Q y V, es variable en la parte curva y decrece a lo largo de la parte recta, aunque todos los valores definen el mismo periodo de retorno conjunto. En resumen, las parejas de valores de las rectas asíntotas tienen probabilidades de ocurrencia bajas y por ello no deben ser incluidos en los análisis de búsqueda de las crecientes (Q y V) críticas o severas.

Test de Wald-Wolfowitz

Esta prueba no paramétrica ha sido utilizada por ^{Bobée & Ashkar (1991)}, ^{Rao & Hamed (2000)} y ^{Meylan et al. (2012)} para probar independencia y estacionariedad en registros de gastos máximos anuales (X_i). Permite verificar su aleatoriedad, cuando su estadístico U no excede de 1.96, en una prueba con nivel de significancia de 5 %.

Estimación de probabilidades empíricas

Las probabilidades de no excedencia empíricas univariadas y bivariadas se estimaron con base en la fórmula de Cunnane, que de acuerdo con ^{Stedinger (2017)} conduce a probabilidades de no excedencia (p) aproximadamente insesgadas, su expresión es:

siendo i el número del dato cuando se ordenan de menor a mayor y n su número total.

Para la estimación de las probabilidades empíricas bivariadas se siguió el mismo principio que aplica para la ecuación 20 (^{Yue et al., 1999}; ^{Yue, 2000b}; ^{Yue & Rasmussen, 2002}), por ello se trabajó en el plano bidimensional, con los datos ordenados en forma progresiva; los gastos pico (Q) en los renglones y los volúmenes (V) en las columnas. El plano formado es un cuadrado de n por n casillas, con n casillas en su diagonal principal, cuando el número de orden del renglón es igual al de la columna. Después cada pareja de datos anuales (Q y V) se localiza en el plano bidimensional y la casilla definida por la intersección del renglón y columna se identifica con el número i que corresponde al año histórico dibujado.

Cuando las n parejas de datos están dibujadas, se busca el año 1 y se define un área rectangular o cuadrada de valores menores de Q y de V, cuyo conteo de casillas numeradas dentro, es NM₁ o combinaciones de Q y V menores. Calculados los n valores de NM_i, se aplica la ecuación 20 para calcular la probabilidad empírica conjunta:

Validación de la distribución GVEb

Esta es la etapa más importante del proceso de ajuste de la GVEb, pues en ella se verifica que tal modelo reproduzca fielmente las probabilidades conjuntas observadas (ecuación 21). ^{Yue (2000a)} indica que la forma más simple de representar las probabilidades conjuntas empíricas y teóricas, consiste en llevar al eje de las abscisas la primera y al eje de las ordenadas la segunda; lógicamente cada pareja de datos define un punto que coincide o se aleja de la recta a 45º.

^{Yue (2000b)} y ^{Yue & Rasmussen (2002)} aplican el test de Kolmogorov-Smirnov con un nivel de significancia (α) de 5 % para aceptar o rechazar la diferencia máxima absoluta (dma) entre las probabilidades conjuntas. Para evaluar la estadística (D_n) del test, se utilizó la expresión que exponen ^{Meylan et al. (2012)}, para α = 5 % es:

n es el número de datos. Si la dma es menor que D_n se acepta la GVEb.

Verificación de la aleatoriedad

A los registros de gasto pico y volumen anual de la Tabla 1 se les aplicó el test de Wald-Wolfowitz, para probar su independencia y estacionariedad. Se encontró que ambas series son aleatorias, con U = 0.726 y U = 1.054.

Ajuste de la distribución GVE

A los registros de la Tabla 1, se les ajustaron distribuciones GVE con los métodos de sextiles, momentos L y máxima verosimilitud. En la Tabla 2 se muestran los estadísticos básicos y los valores de los parámetros de ajuste de cada distribución GVE.

Búsqueda de los parámetros óptimos de la GVE_B

Debido a la gran similitud que mostraron los parámetros de ajuste de la distribución GVE en la Tabla 2, con los tres métodos aplicados para ambos registros, se decidió adoptar a los valores inferiores, iniciales y superiores de cada variable de decisión, con base en tales resultados; quedando así: u₁ (160,167,175), α₁ (125,138,148), k₁ (-0.50,-0.52,-0.54), u₂ (38,40,49), α₂ (42,43,47) y k₂ (-0.512,-0.535,-0.545). Además, se definieron cinco intervalos para el parámetro de asociación (m) y se adoptó como inicial el valor medio, tales límites fueron: 1.0 a 1.5, 1.5 a 2.0, 2.0 a 2.5, 2.5 a 3.0 y 3.0 a 3.5. En la Tabla 3 se muestran los resultados principales de las cinco corridas numéricas del algoritmo Complex. Los valores de la última columna se estimaron como se detalla en seguida.

Con base en los parámetros de ajuste óptimos, expuestos en las columnas 4 a 7 de la Tabla 3, se estimaron con la ecuación 4 las probabilidades de no excedencia conjunta teórica F(x,y), utilizando los datos x_i, y_i de la Tabla 1. Con base en la ecuación 21 y su procedimiento gráfico descrito, se calcularon las llamadas probabilidades de no excedencia empíricas bivariadas F_e(x,y), contra las cuales se contrastan los valores de la F(x,y).

La mejor correspondencia o similitud entre ambas probabilidades conjuntas F_e(x,y) y F(x,y), se obtuvo para la quinta aplicación numérica del algoritmo Complex, con un valor del coeficiente de correlación (r_xy) de 0.9959 y diferencias máximas positiva y negativa de 0.0630 y -0.0658, que se indican sombreadas en la Tabla 4.

Ratificación de las marginales

Primeramente, los gastos pico y volúmenes de la Tabla 1 se ordenaron de menor a mayor. Después, se calcularon sus probabilidades de no excedencia teóricas con las ecuaciones 2 y 3 utilizando los parámetros de ajuste óptimos de la Tabla 3 para la quinta corrida del algoritmo Complex. Las probabilidades de no excedencia empíricas de ambas series (Q y V) se estimaron con la ecuación 20.

En las Figuras 3 y 4 se muestra el contraste gráfico de probabilidades para cada serie ordenada. Las diferencias máximas entre probabilidades empíricas y teóricas de los gastos pico y volúmenes fueron 0.0912 y -0.0674; la primera ocurrió en el dato ordenado número 14 y la segunda en el 3, como se muestra en la Tabla 5. Como ambas diferencias absolutas son menores que D_n = 0.1883 (ecuación 22), se acepta que los registros de Q y V de la Tabla 1 tienen las marginales GVE definidas con el algoritmo Complex.

Figura 3 Distribución marginal GVE del gasto pico anual de las crecientes anuales de entrada a la Presa Venustiano Carranza, México 

Figura 4 Distribución marginal GVE del volumen anual de las crecientes anuales de entrada a la Presa Venustiano Carranza, México 

Verificación de las restricciones de probabilidad

Antes de proceder a estimar los periodos de retorno conjuntos de diseño T’(Q,V), es conveniente verificar la ecuación 7, que establece las restricciones de la probabilidad conjunta. Lo anterior se muestra en la Tabla 6 para un número reducido de parejas de datos históricos.

Se observa en la Tabla 6 y en su versión completa, que el valor de la columna 5 es siempre mayor que el de la 4 y menor que el más pequeño de las columnas 6 ó 7. Por lo cual, se cumple con la restricción que impone la ecuación 7.

Predicciones univariadas y periodos de retorno conjuntos de diseño

Aplicando en las ecuaciones 5 y 6 los seis parámetros de ajuste óptimos de la quinta corrida del algoritmo Complex, mostrados en la Tabla 3, se obtuvieron las predicciones mostradas en la Tabla 7.

En esta aplicación numérica, las predicciones univariadas y los periodos de retorno conjuntos de diseño, fueron establecidos con apego estricto a sus ocurrencias reales. Por lo anterior, los datos de Q y V de la Tabla 1, se ordenaron de menor a mayor, para detectar sus valores extremos o dispersos (outliers); es decir, los que se apartan de manera notable de la tendencia general. Para el gasto pico se definió uno de 4320.7 m³/s y para el volumen se obtuvieron dos: uno de 529.11 y otro de 983.02 Mm³.

Debido a la amplitud del registro conjunto, de 52 años, las primeras predicciones por revisar de la Tabla 7 son del periodo de retorno de 50 años, cuyos valores de 1671 m³/s y 561 Mm³ se consideran acordes con los datos extremos observados. Para los periodos de retorno mayores, resultan aceptables, como máximo, las de 500 años, cuya predicción de gasto pico de 5499 m³/s puede aceptarse factible de ocurrir; en cambio, la del volumen es casi del doble del máximo observado. Por lo anterior, los periodos de retorno conjuntos de diseño serán los cuatro siguientes: 50, 100, 250 y 500 años.

Eventos de diseño obtenidos con regresión

El diagrama de dispersión de las 52 parejas de datos originales (Tabla 1) mostró una nube de puntos con tendencia lineal con un coeficiente de correlación lineal (r_xy) de 0.9188 La ecuación de regresión lineal que la representa es la siguiente (^{Campos, 2003}):

De la Tabla 7 se obtienen las predicciones siguientes para el gasto pico (Q) y los cuatro periodos de retorno conjuntos de diseño: 1671, 2406, 3861 y 5499 m³/s. Con base en la ecuación 23 se definen los volúmenes anuales (V) siguientes: 433, 617, 980 y 1389 Mm³, para los eventos de diseño buscados. Estas cuatro parejas de valores de Q y V se han dibujado en la Figura 1 y se indican con la letra “r”.

Eventos de diseño condicionales tipo T(Q|V)

Están definidos por las ecuaciones 18 y 19, cuya aplicación emplea las expresiones 2 a 4. Para los cuatro periodos de retorno conjuntos de diseño definidos, se obtienen de la Tabla 7 los siguientes cuatro gastos pico: 1671, 2406, 3861 y 5499 m³/s. Adoptando tales gastos como valores condicionantes (X ≤ x), se procedió por tanteos del volumen (y) a estimar, con la ecuación 19, el periodo de retorno condicional que debe igualar al del gasto pico. Los volúmenes estimados fueron: 386, 571, 940 y 1360 Mm³. Estas cuatro parejas de valores de Q y V se han dibujado en la Figura 1 y se indican con la letra “c”.

Gráficas del periodo de retorno conjunto T’(Q,V)

Los periodos de retorno conjuntos de tipo AND se estiman con base en la ecuación 15. Para cada periodo de retorno conjunto de diseño se seleccionan arbitrariamente volúmenes y gastos pico, para obtener sus probabilidades de no excedencia marginales (ecuaciones 2 y 3) y conjunta (ecuación 4). En la Tabla 8 se muestran algunos resultados para definir las cuatro gráficas de la Figura 1.

En la Figura 1 o en la Tabla 8 se pueden seleccionar infinitas parejas de Q y V, que satisfacen el periodo de retorno conjunto de diseño y que definen como subgrupo de parejas críticas las que están dentro de la porción curva de cada gráfica de T’(Q,V), fuera de las rectas asíntotas (^{Volpi & Fiori, 2012}).

Las combinaciones de gasto pico y volumen que tienen el mismo periodo de retorno conjunto, establecen crecientes o hidrogramas que producirán diferentes efectos en el embalse que se diseña o revisa; adoptando por seguridad, el que genera las condiciones más críticas, severas o desfavorables. Lo anterior, está incorporando en el diseño hidrológico las características físicas del vertedor y almacenamiento o vaso del embalse en proyecto o bajo revisión. Para formar cada hidrograma de diseño, existen métodos teóricos y empíricos (^{Aldama, 2000}; ^{Ramírez & Aldama, 2000}).