Estimación del factor de frecuencia K_t para la función de distribución de probabilidad gamma incompleta

Selección de las estaciones meteorológicas

Se seleccionaron del programa ERIC III v.2. publicado por el Instituto Mexicano de Tecnología del Agua (IMTA, 2009)Instituto Mexicano de Tecnología del Agua. (2009). ERIC III: Extractor Rápido de Información Climatológica v.2 [CD-ROM]., los valores de precipitación total mensual de seis estaciones meteorológicas de varias regiones de México, con registros que variaron de 31 a 59 años (Tabla 1). Las estaciones seleccionadas están enclavadas en áreas con agricultura de temporal de ambos litorales y del centro de la República Mexicana; con altitudes desde 7 a 2,216 msnm; latitudes desde 14.7714° a 22.500°; y Longitudes desde -89.284° a -105.367° y con precipitaciones totales anuales entre 563.0 a 1889.4 mm.

La base teórica de la FDP Gamma Incompleta y sus parámetros se reproducen en las siguientes líneas, con el único propósito de establecer el procedimiento utilizado en este estudio. La Función de Densidad de la Gamma Incompleta se muestra en la siguiente ecuación (1) $f\left(x,\alpha ,\beta \right)= \frac{1}{{\beta }^{\alpha }\Gamma \left(\alpha \right)} x^{\alpha -1}exp\left(-\frac{x}{\beta }\right)$ y su representación gráfica se reporta en Wilks, (2006)Wilks, D. S. (2006). Statistical methods in the atmospheric sciences (2.a ed.). Academic Press. https://www.sciencedirect.com/book/9780128158234/statistical-methods-in-the-atmospheric-sciences.

Esta ecuación es válida para variables continuas, con valores iguales o mayores a cero, como son los valores de la precipitación.

La FDP Gamma Incompleta o Gamma de dos parámetros, muestra una gran capacidad para representar varias formas de diferentes distribuciones de probabilidad en función de los valores de sus parámetros de forma, alfa (α) y escala, beta (β). Así, cuando el valor ʼαʼ toma valores iguales o menores a 1.0 se aproxima a la FDP Exponencial. En cambio, cuando el parámetro de forma ʼαʼ se incrementa por encima de 50, la distribución tiende a desplazarse hacia la derecha y se aproxima a la FDP Normal (Arroyo et al., 2014Arroyo, I., Bravo, L. C., Llinás, H. & Muñoz, F. L. (2014). Poisson and gamma distributions: A discrete and continuous relation. Prospectiva, 12(1), 99-107. https://doi.org/10.15665/rp.v12i1.156 y Wilks, 2006Wilks, D. S. (2006). Statistical methods in the atmospheric sciences (2.a ed.). Academic Press. https://www.sciencedirect.com/book/9780128158234/statistical-methods-in-the-atmospheric-sciences).

En el presente estudio, se utilizará la FDP Gamma Incompleta para calcular las probabilidades de excedencia de las precipitaciones mensuales, por las capacidades de representación de dicha FDP antes mencionadas.

El valor de cualquier precipitación en particular (X_i), se puede representar como la suma del valor de la media (μ), más la desviación respecto a la media (∆_XT):

Dónde:

Dónde, ‘σ’ es la desviación estándar (mm) y K_T es un factor de frecuencia correspondiente a un periodo de retorno y/o probabilidad de excedencia que tiene esa precipitación. Resultando:

Para muestras se puede aproximar a:

Esta es la ecuación de la variable reducida que se utiliza para calcular la precipitación a un nivel de probabilidad de excedencia deseado; dónde, X_T es la precipitación total mensual a un cierto periodo de retorno y/o probabilidad de excedencia (mm); P_N es el valor promedio de la precipitación de un conjunto de datos; K_T es el factor de frecuencia y ‘s’ es su correspondiente desviación estándar.

En esta ecuación (5) $X_T=P_N+\left(K_T*s\right)$ , la única incógnita es el valor de K_T , cuyo valor depende del periodo de retorno y/o probabilidad de excedencia del valor de ‘X_T ’ (Chow, 1951Chow, V. T. (1951). A general formula for hydrologic frequency analysis. Eos, Transactions American Geophysical Union, 32(2), 231-237. https://doi.org/10.1029/TR032i002p00231).

Es importante señalar que, los valores de K_T que se estimen, son aplicables para las FDP’s Gamma Incompleta y Log-Gamma de 2 parámetros.

Para encontrar el valor correspondiente del factor de frecuencia K_T para cada valor de la probabilidad selecta de la FDP Gamma Incompleta, en el presente estudio se desarrolló un procedimiento a partir de los valores del parámetro de forma ʼαʻ y del parámetro de escala ʻβʼ de la FDP Gamma Incompleta.

Para otras FDP’s, asimétricas como la Gumbel o la Log-Pearson III (Log-Gamma de tres parámetros), se aplica el mismo modelo de la ecuación (5) $X_T=P_N+\left(K_T*s\right)$ . Solo cambia la manera de obtener el valor del factor de frecuencia K_T , para lo cual existe reportado en la literatura científica modelos que permiten estimar su cálculo, a partir de los valores de ‘Z’ y del coeficiente de asimetría (γ) (Kite, 1977Kite, G. W. (1977). Frequency and risk analyses in hydrology. Water Resources Publications. https://es.scribd.com/document/487410158/Frequency-and-Risk-Analyses-in-Hidrology-G-W-Kite).

Probabilidades de excedencia

Para el cálculo de las probabilidades de excedencia con la FDP Gamma Incompleta, se utilizaron registros mensuales completos de al menos 30 años. Al conjunto de valores de precipitación de un mes en particular, se le determinó las probabilidades de excedencia del 10, 20, 30, 40, 50, 60, 70, 80 y 90%. Se seleccionó la estación meteorológica 27039 de Samaria, Cunduacán, Tabasco, para realizar la calibración y ajuste de los valores de K_T de la ecuación (5) $X_T=P_N+\left(K_T*s\right)$ , debido a que cuenta con un historial de medición de precipitación ininterrumpido y porque es la más representativas del estado. El resto de las cinco estaciones, se utilizaron para la validar los valores de K_T al utilizar dicho modelo en localidades con agricultura de temporal con diferentes climas, para constatar que el procedimiento propuesto puede ser utilizado en otras zonas de México.

Procedimiento para obtener el factor de frecuencia K_T

La FDP Gamma Incompleta o Gamma de dos parámetros (‘α’ y ʻβʼ), se utilizó porque la precipitación total mensual fundamentalmente tiene una distribución asimétrica con sesgo positivo y la FDP Gamma Incompleta es perfecta para representar esos datos (Wilks, 2006Wilks, D. S. (2006). Statistical methods in the atmospheric sciences (2.a ed.). Academic Press. https://www.sciencedirect.com/book/9780128158234/statistical-methods-in-the-atmospheric-sciences). Para la determinación de los parámetros ‘α’ y ʻβʼ, se recomienda el método de máxima verosimilitud, mismo que se muestra en las ecuaciones (6 $D=\ln \left(X_m\right)- \frac{1}{N} {\displaystyle \sum }_{i=1}^{N}ln\left(x_i\right)$ , 7 $\alpha = \frac{1+ \sqrt{1+ \frac{4D}{3}}}{4D}$ y 8 $\beta = \frac{X_m}{\alpha }$ ) reportadas Wilks (2006)Wilks, D. S. (2006). Statistical methods in the atmospheric sciences (2.a ed.). Academic Press. https://www.sciencedirect.com/book/9780128158234/statistical-methods-in-the-atmospheric-sciences y por Arroyo et al. (2014)Arroyo, I., Bravo, L. C., Llinás, H. & Muñoz, F. L. (2014). Poisson and gamma distributions: A discrete and continuous relation. Prospectiva, 12(1), 99-107. https://doi.org/10.15665/rp.v12i1.156.

Dónde; D es la diferencia entre el logaritmo natural de la media aritmética (X_m) y la media de los logaritmos de las observaciones de la precipitación; (X_i). El valor de ‘α’ se calcula con la ecuación (7) $\alpha = \frac{1+ \sqrt{1+ \frac{4D}{3}}}{4D}$ .

Dónde ‘α’ es el parámetro de forma y ‘β’ es el parámetro de escala, que se calcula con la ecuación (8) $\beta = \frac{X_m}{\alpha }$ .

A partir de la ecuación (5) $X_T=P_N+\left(K_T*s\right)$ , se despeja el factor K_T y se obtiene la ecuación (9) $K_T= \frac{\left(X_{T }- P_N\right)}{s}$ siguiente:

Estimación del factor de frecuencia K_T para la FDP Gamma Incompleta

Para obtener el valor del factor K_T en función de los parámetros de la FDP Gamma Incompleta para cualquier probabilidad de excedencia, primero se calculan los parámetros de forma ʻαʼ y de escala ʻβʼ, con las ecuaciones (10) $E\left(X\right)=P_N=\alpha *\beta$ , (11) $P_N=\alpha *\beta$ y (12) $s^2= \alpha *{\beta }^2$ . Para ello, es conveniente recordar que la esperanza matemática o media [E(x)] de la FDP Gamma Incompleta de acuerdo con Arroyo et al. (2014)Arroyo, I., Bravo, L. C., Llinás, H. & Muñoz, F. L. (2014). Poisson and gamma distributions: A discrete and continuous relation. Prospectiva, 12(1), 99-107. https://doi.org/10.15665/rp.v12i1.156 y Wilks (2006)Wilks, D. S. (2006). Statistical methods in the atmospheric sciences (2.a ed.). Academic Press. https://www.sciencedirect.com/book/9780128158234/statistical-methods-in-the-atmospheric-sciences se representa por la siguiente ecuación (10) $E\left(X\right)=P_N=\alpha *\beta$ :

Dónde, P_N es la precipitación promedio o normal climatológica. P_N se puede despejar de la ecuación (14) $K_T= \frac{\left[X_T - \left(\alpha *\beta \right)\right]}{\sqrt{\alpha *{\beta }^2}}$ y se obtiene la ecuación (11) $P_N=\alpha *\beta$ siguiente:

La varianza de la FDP Gamma Incompleta [Var(x)], la reportan Arroyo et al. (2014)Arroyo, I., Bravo, L. C., Llinás, H. & Muñoz, F. L. (2014). Poisson and gamma distributions: A discrete and continuous relation. Prospectiva, 12(1), 99-107. https://doi.org/10.15665/rp.v12i1.156; y Johnson et al. (2018)Johnson, R. A., Miller, I. & Freund, J. E. (2018). Probability and statistics for engineers (9.a ed.). Pearson Education Ltd., como se muestra en la ecuación (12) $s^2= \alpha *{\beta }^2$ :

Así, la desviación estándar ‘s’ será:

Sustituyendo en la ecuación (13) $s= \sqrt{\alpha *{\beta }^2}$ los equivalentes de P _N de la ecuación (11) $P_N=\alpha *\beta$ y ‘s’ de la ecuación (13) $s= \sqrt{\alpha *{\beta }^2}$ , ésta se transforma en la ecuación (14) $K_T= \frac{\left[X_T - \left(\alpha *\beta \right)\right]}{\sqrt{\alpha *{\beta }^2}}$ siguiente:

De esta manera, el factor de frecuencia K_T queda en función de la precipitación con la probabilidad deseada (X_T) y los parámetros ʻαʼ y ʻβʼ de la FDP Gamma Incompleta. Posteriormente, este valor se introduce en la ecuación (5) $X_T=P_N+\left(K_T*s\right)$ para calcular aquella precipitación con cierto nivel de probabilidad de excedencia deseado (X_T).

Cálculo de la precipitación a una determinada probabilidad de excedencia (X_T)

Para estimar el valor de K_T de la ecuación (14) $K_T= \frac{\left[X_T - \left(\alpha *\beta \right)\right]}{\sqrt{\alpha *{\beta }^2}}$ , se calculó el valor de X_T para una determinada probabilidad de excedencia, utilizando el comando de Excel Microsoft ®. DISTR.GAMMA.INV (probabilidad, alfa, beta). Como Excel Microsoft ® calcula automáticamente la probabilidad de No-excedencia, para obtener la probabilidad de excedencia se partió de la siguiente igualdad: Probabilidad de excedencia = 1 - probabilidad de no-excedencia. El valor de X_T resultante fue introducido en la ecuación (14) $K_T= \frac{\left[X_T - \left(\alpha *\beta \right)\right]}{\sqrt{\alpha *{\beta }^2}}$ para obtener el valor de K_T para cada uno de los meses del año. A diferencia de una distribución simétrica, aquí los valores K_T son diferentes para cada uno de los 12 meses y para una probabilidad escogida.

Para obtener un solo valor de K_T para todo el año, se calculó la mediana de esos 12 valores, por ser una mejor representación de la tendencia central, en comparación con la media aritmética, ya que la FDP Gamma Incompleta se distribuye asimétricamente. De esta manera, a cada probabilidad de excedencia se le asoció un único valor de K_T . Así, el valor único del factor K_T resultante para una probabilidad de excedencia escogida, se introdujo en la ecuación (5) $X_T=P_N+\left(K_T*s\right)$ para calcular el valor de la precipitación mensual correspondiente (X_T) a esa probabilidad de excedencia, conocidos su media (P_N) y desviación estándar mensual (s).

Cálculo de la precipitación a una cierta probabilidad de excedencia

En el presente trabajo se calcularon los valores del factor K_T para diferentes probabilidades de excedencia a intervalos de 0.05. Acto seguido, se procedió a realizar una regresión polinómica entre los valores de la probabilidad de excedencia y el valor correspondiente del factor de frecuencia K_T . La ecuación de regresión resultante permite calcular el factor K_T a cualquier valor de probabilidad de excedencia deseado. Además, esto permitió elaborar una tabla donde se asocia el valor de K_T a una determinada probabilidad de excedencia, a intervalos de 0.05 de probabilidad; misma que se presenta en los resultados de este estudio.

Comparación entre valores estimados y observados de la precipitación

Al aplicar el modelo de la ecuación (5) $X_T=P_N+\left(K_T*s\right)$ se obtuvieron los valores correspondientes de precipitación para las diferentes probabilidades de excedencia. Estos valores de precipitación estimados con este primer procedimiento, se compararon mediante un análisis de regresión lineal simple, con los valores obtenidos de la FDP Gamma Incompleta, utilizando el comando de Excel Microsoft ® anterior, para determinar qué tan bueno fue el ajuste al utilizar los valores resultantes de K_T .

Calibración, ajuste y validación de los valores de K_T

Para calibrar los valores de K_T , se utilizaron datos de precipitación total mensual de 52 años de registro, de la estación meteorológica 27039 de Samaria, Cunduacán, Tabasco, México. Para validar los diferentes valores del factor K_T para sus correspondientes probabilidades de excedencia, se utilizaron datos de precipitación mensual de las cinco estaciones meteorológicas restante de la Tabla 1.

Representación de los resultados

A las precipitaciones resultantes para las diferentes probabilidades de excedencia obtenidas con el modelo de la variable reducida y con la FDP Gamma Incompleta, se les realizó un análisis de regresión simple entre ambos datos resultantes. Posteriormente, se elaboraron las gráficas de la distribución temporal de ambas precipitaciones a una determinada probabilidad y las gráficas del análisis de regresión con su correspondiente coeficiente de determinación. También se elaboraron las tablas correspondientes que muestran las ecuaciones de regresión, sus coeficientes de determinación (R²) y la raíz del cuadrado medio del error (RMSE), para cada una de las 6 estaciones selectas.

Índice estadístico para evaluar el desempeño o bondad de ajuste de los modelos

Después del ajuste del modelo se evaluó su desempeño, comparando los valores predichos con el modelo con los datos observados. Para determinar qué tan preciso es el modelo de regresión obtenido para predecir un valor con los modelos de regresión obtenidos, se utilizó el estadístico raíz del error cuadrático medio (RMSE por sus siglas en inglés), que evalúa la relación entre los valores resultantes de la Distribución Gamma Incompleta, con los valores estimados por el modelo propuesto; cuya ecuación la reporta Stöckle et al. (2004)Stöckle, C. O., Kjelgaard, J. & Bellocchi, G. (2004). Evaluation of estimated weather data for calculating Penman-Monteith reference crop evapotranspiration. Irrigation Science, 23(1), 39-46. https://doi.org/10.1007/s00271-004-0091-0.

Dónde, RMSE es la raíz cuadrada del error cuadrático medio (mm); ∑ es el símbolo de sumatoria; N es el número de observaciones; O_i son los valores observados de precipitación (mm); y E_i son los valores estimados de precipitación.

Calibración y ajuste a la FDP Gamma Incompleta

En la Tabla 3, se muestra los valores mensuales de la precipitación calculada con el modelo propuesto y la FDP Gamma Incompleta, para una probabilidad de excedencia de 0.10 (10%), para el área de influencia de la estación meteorológica 27039 de Samaria, Cunduacán, Tabasco, México. Utilizada para realizar la calibración y ajuste del modelo propuesto. Asimismo, en la tabla se muestra el modelo propuesto para calcular la precipitación con una probabilidad de excedencia de 0.10 con un valor del factor K_T de 1.3400.

La Figura 2, muestra de manera gráfica los resultados del comportamiento temporal de los valores de la precipitación obtenidos por ambos procedimientos para una probabilidad de excedencia de 0.10. Puede apreciarse que los resultados con ambas metodologías son muy similares a este nivel de probabilidad y que el modelo propuesto da iguales resultados que la FDP Gamma incompleta.

Los resultados del análisis de regresión entre los valores de precipitación calculados con el modelo propuesto y la función acumulada Gamma Incompleta se muestran en la Figura 3, con un ajuste casi perfecto.

La Tabla 4 muestra los resultados de la precipitación de excedencia con una probabilidad de excedencia de 0.50 (50 %), obtenidos mediante el modelo propuesto y la FDP Gamma Incompleta. Además, muestra el modelo propuesto para calcular la precipitación con una probabilidad de excedencia de 0.50 con un valor del factor K_T de -0.1879.

La Figura 4, se muestra la distribución temporal de los valores de la precipitación con una probabilidad de excedencia de 0.50 obtenidos con ambas metodologías. La Figura 5 muestra que el modelo propuesto da valores muy similares a los obtenidos con la función acumulada Gamma Incompleta.

Los resultados del análisis de regresión entre los valores de precipitación calculados con el modelo propuesto y la FDP Gamma Incompleta se muestran en la Figura 5. En la que se puede apreciar un muy buen ajuste (Coeficiente de Determinación de 0.9974).

Finalmente, y a manera de ejemplo y para la misma estación meteorológica utilizada para calibrar y ajustar el modelo propuesto, en la Tabla 5 se muestran los resultados de la precipitación de excedencia con una probabilidad de 0.80 (80%), obtenidos mediante el modelo propuesto y la función acumulada Gamma Incompleta. En esa misma Tabla 5, se muestra el modelo propuesto para calcular la precipitación con una probabilidad de excedencia de 0.80 con un valor del factor K_T de -0.8457.

En la Figura 6, se muestra la distribución temporal de los valores de la precipitación con una probabilidad de excedencia de 0.80 obtenidos con ambas metodologías, nuevamente se aprecia que el modelo propuesto da valores muy similares de precipitación a los obtenidos con la función acumulada Gamma Incompleta a este nivel de probabilidad de excedencia.

Los resultados del análisis de regresión entre los valores de precipitación calculados con el modelo propuesto y la función acumulada Gamma Incompleta se muestran en la Figura 7. En la que se puede apreciar un muy buen ajuste (Coeficiente de Determinación de 0.9982) a este nivel de probabilidad de excedencia.

Al aplicar los valores del factor K_T de la Tabla 2, a la ecuación (5) $X_T=P_N+\left(K_T*s\right)$ se obtuvieron los valores de precipitación para las probabilidades de excedencia correspondientes a cada 0.10 de intervalo. Estos valores se compararon mediante un análisis de regresión lineal simple con los valores obtenidos de la FDP Gamma Incompleta, para el área de influencia de la estación meteorológica 27039 Samaria, Cunduacán, Tabasco México. Las ecuaciones de regresión y los correspondientes valores de los coeficientes de determinación y el estadístico Raíz del Error Cuadrático Medio (RMSE por sus siglas en inglés), se muestran en el Tabla 6.

Se puede apreciar en la Tabla 6, que el valor de la pendiente es muy cercano al valor de uno. Lo que indica que el valor de precipitación calculado por el modelo, es casi igual al calculado por la función acumulada Gamma Incompleta. Esto se confirma al observar al coeficiente de determinación que siempre fue superior a 0.993, con un error cuadrático medio (RMSE) inferior a los 5.5 mm mensuales, que es insignificante.

Validación del modelo en otras localidades de México

Para validar el modelo general con diferente valor del factor K_T para sus correspondientes probabilidades de excedencia, y constatar la generalización en la aplicación del método propuesto se seleccionaron otras cinco localidades de México, con clima diferente al de la estación meteorológica 27039 Samaria, Cunduacán, Tabasco. Las localidades seleccionadas fueron: Motul de Felipe Carrillo Puerto; Texcoco; Suchiate; Zapopan y Acaponeta, con 45, 31, 35, 48 y 59 años de registro, respectivamente (Tabla 1). Los resultados se muestran en las Tablas 7, 8, 9, 10 y 11.

En la Tabla 7, se puede observar que todos los modelos tuvieron un R² ≥ 0.984, con variación de la RMSE entre 1.1 a 9.1 mm y con pendiente cercana a la unidad (0.959 a 1.00), siendo la menor a la mayor probabilidad. Es decir que la probabilidad de 0.30 presentó el valor más alto de la RMSE (9.1 mm), mientras que la probabilidad de 0.70 presentó el valor más bajo de la RMSE (1.1 mm). En los meses más secos de una localidad disminuye el R² y aumenta la RMSE.

La Tabla 8, muestra que todos los modelos de regresión entre los valores resultantes con el modelo propuesto y la distribución acumulada Gamma incompleta, tuvieron un R² ≥ 0.998; una RMSE que varió entre 0.5 a 2.4 mm y una pendiente cercana a la unidad (0.95 a 1.01) siendo la menor a la mayor probabilidad. La probabilidad de excedencia de 0.30 presentó el valor más alto de la RMSE (2.4 mm), y la probabilidad de excedencia de 0.70 presentó el valor más bajo de la RMSE (0.5 mm).

La Tabla 9, muestra que todos los modelos tuvieron una R² ≥ 0.996; la RMSE varió entre 1.1 a 13.4 mm; una pendiente cercana a la unidad (0.95 a 1.07) siendo la menor a la mayor probabilidad. La probabilidad de 0.40 presentó el valor más alto de la RMSE (13.4mm), mientras que la probabilidad 0.80 presento el valor más bajo de la RMSE (1.1 mm).

La Tabla 10, muestra que todos los modelos tuvieron una R² ≥ 0.996; la RMSE varió entre 0.9 a 7.5 mm; una pendiente cercana a la unidad (0.95 a 1.02) siendo la menor a la mayor probabilidad. La probabilidad de 0.20 presentó el valor más alto de la RMSE (7.5 mm), mientras que la probabilidad 0.80 presentó el valor más bajo de la RMSE (0.9 mm).

La Tabla 11, muestra que todos los modelos tuvieron una R² ≥ 0.995; la RMSE varió entre 1.8 a 12.9 mm; una pendiente cercana a la unidad (0.95 a 1.02) siendo la menor a la mayor probabilidad. Probabilidades iguales o menores a 0.40 presentaron valores de la RMSE ≥10.0 mm, mientras que la probabilidad 0.80 presento el valor más bajo de la RMSE (1.8 mm).