Populations and sample L-moments

The L-moments are an alternative system to describe the forms of the FDP. Historically they appear as modifications of the weighted probability moments (MPP) developed by ^{Greenwood et al. (1979)}. The MPPs of a random variable x with accumulated FDP F(x) were defined by the following quantities:

Of particular interest is the following special case: β_r = M₁,_r,₀. For a distribution with function of quantiles x(F), equation 1 leads to:

This equation can be contrasted with the ordinary definition of the moments, that is:

The conventional definition of moments involves successive powers of the quantile function x(F), while the MPPs involve successive powers of F and therefore can be considered as integrals of x(F) weighted by polynomials F^r ; hence its name. MPPs substantially improve the properties of sampling, as they are not influenced by scattered values (^{Asquith, 2011}). The L-moments are linear combinations of the MPPs, as follows (^{Hosking and Wallis, 1997}):

In addition, the quotients (τ) of L-moments are defined, starting with L -Cv which is analogous to this coefficient and then those of similarity with the coefficients of asymmetry and kurtosis, these are:

In a sample of size n, with its elements in ascending order (x₁ ≤ x₂ ≤… x_n) the unbiased estimators of β_r are (^{Stedinger et al., 1993}; ^{Hosking and Wallis, 1997}) are:

with the following general expression:

The sampling estimators of λ_r will be l_r being defined by equations 4 a 7 and those of the ratios of L-moments will be t, t₃ and t₄, according to equations 8 to 10.

Fitting with L-moments of the distributions GVE, LOG, and PAG

^{Hosking and Wallis (1997)} highlighted in their “Table 5.1”, that these three FDPs, when their shape parameter is negative (k < 0), have their right tails thicker or denser than all other FDPs commonly used in the AFC; due to this they have gained acceptance in the analysis of frequencies of extreme hydrological data (^{El Adlouni et al., 2008a}). These three PDFs also coincide in having an upper limit when k > 0 and defining functions of two fit parameters known as Gumbel, Logistics and Exponential, when k = 0.

Next, the inverse solution x(F) for the GVE, LOG and PAG distributions is cited, with which the predictions associated with a certain probability of non-exceedance (F) are estimated and the equations that allow us to estimate their three fit parameters (u, a, k) corresponding to the location, scale and form, with the method of L-moments.

Distribution GVE (^{Hosking and Wallis, 1997}): interval of x: u + a/k ≤ x < ∞ si k < 0; - ∞< x < ∞ si k = 0; - ∞ < x ≤ u + a/k si k > 0.

being:

For the evaluation of the Gamma function the Stirling formula (^{Davis, 1972}) was used:

Distribution LOG (^{Hosking and Wallis, 1997}): interval of x, identical to that of the GVE.

For distribution PAG (^{Hosking and Wallis, 1997}) the interval of x is u ≤ x < ∞ if k ≤ 0; u ≤ x ≤ u + a/k if K > 0.

Fitting with L-moments of the distribution GVE₁ and GVE₂

The generalization of the L-moments method proposed by ^{El Adlouni and Ouarda (2008b)} for the fitting of the non-stationary FPD type GVE₁ begins by analyzing the expected value according to the expression:

From this equation it is deduced that an estimator δ^2 of δ₂ can be obtained by a simple linear regression between the variable X and the covariate t. Then a new variable S₁ is defined with the expression:

It is deduced from equation 29, that the new variable is distributed according to a FDP type GVE₀ with parameters δ₁, a and k that are estimated with equations 16 to 19 of the L-moments method, highlighting that u (equation 19) is equal to δ₁, with which the four fitting parameters of the GVE₁ model are estimated. Equation 30 corresponds to one of the first simple approaches suggested to process hydrological records with tendency, consisting in removing first such a deterministic component (^{McCuen and Thomas, 1990}, ^{Campos-Aranda, 2012}). As ^{Mudersbach and Jensen (2010)} have indicated, such an approach is practical but the results of its AFC are valid only in the present and hydraulic works must be safe at the end of their useful life, requiring that the flood of design be estimated in a predetermined future date.

The same approach is used to introduce a quadratic dependence on the location parameter u, therefore, the estimates δ^2 and δ^3 of δ₂ and δ₃ are obtained by means of a multiple linear regression of the variable X against t and t², or, a second-degree polynomial regression of X against t. The new variable S₂ will be:

With the equations 16 to 19 applied to the corrected data sample S₂, the remaining parameters k, a and δ₁ of the GVE₂ model are defined. ^{El Adlouni and Ouarda (2008b)} compared by numerical simulation three procedures for obtaining the four and five fitting parameters of the GVE₁ and GVE₂ models: the maximum likelihood method (^{Coles, 2001}, ^{Nadarajah, 2005}, ^{Katz, 2013}), generalized maximum likelihood (^{Martins and Stedinger, 2000}; ^{El Adlouni et al., 2007}) and their generalization of the L-moments; they conclude that the latter is better than the first, because it has lower bias and mean square error, but it does not exceed the second, especially in records with important asymmetry.

Slope δ₂ according to the Pranab Kumar Sen criterion

To verify numerically the value obtained from δ₂ with equation 34, the ^{Sen criterion (1968)} was used, which utilizes the formula of the trend established in the Kendall test, to estimate its slope, defining it as the median value (MED) of the partial slopes, that is:

X_j and X_i are the data in the j and i times, which represent years: being n the number of data, then we have n(n - 1)/2 partial slopes. δ₂ has units of X/year and its positive sign defines ascending tendencies and negative descending ones. ^{Machiwal and Jha (2012)} indicate that this criterion is resistant or robust to the presence of scattered values (outliers).

Equations of the quadratic regression

The quadratic or parabolic polynomial for the nonlinear trend between the X_i data and the years of the t_i record is (^{Campos-Aranda, 2003}):

whose matrix arrangement of its normal equations is:

with solution:

being, T a square matrix with inverse T^-1, δ a vector column of unknowns and X another vector column of independent terms. The coefficient of determination R² quantifies the degree of correlation of the polynomial fitting, its numerator measures the improvement or reduction of the error due to the regression and its denominator is the dispersion of the dependent variable, that is:

Being:

Since Se² is always smaller than Sx², then R² varies from zero to unity (when Se² = 0) and its square root corresponds to the coefficient of polynomial correlation (^{Campos-Aranda, 2003}).

Fitting with L-moments of the distributions LOG₁, LOG₂, PAG₁, and PAG₂

In ^{Rao and Hamed (2000)} the equations of the expected value of the Generalized Logistics (LOG) and Generalized Pareto (PAG) distributions can be found, similar to the expression 29 of the GVE function. This implies that the generalization of ^{El Adlouni and Ouarda (2008b)} of the L-moments method to fit non-stationary FDP with variable location parameter in time is also applicable in the LOG and PAG models.

Standard error of fit

In the mid-1970s the standard error of fit (EEA) was established as a quantitative statistical indicator, since it evaluates the standard deviation of the differences between the observed values and those estimated with the FDP that is tested; in this study the models: GVE₁, GVE₂, LOG₁, LOG₂, PAG₁ and PAG₂. Their expression is as follows (^{Kite, 1977}, ^{Pandey and Nguyen, 1999}):

in which, n and np are the number of the sample data and of fitting parameters, in this case four and five; X_i are the ordered data from least to greatest and are estimated values with the inverse solution x(F) or quantile function that uses the variable location parameter, for a probability of non-exceedance estimated with the Weibull formula (^{Benson, 1962}):

where, m is the data order number, with 1 for the smallest and n for the largest one.

Records of annual maximum values to be processed

Four records of the specialized literature were selected to illustrate the application of the non-stationary FDPs GVE, LOG and PAG with the generalization method of the L-moments (^{El Adlouni and Ouarda, 2008b}); their values are shown in Table 1. These records have approximate data read in the graphs where they were exposed. The first record was shown by ^{Leclerc and Ouarda (2007)} as an example of a series of flows (m³·s^-1) with a downward trend in a hydrometric station of the Dartmouth River in Canada, with 630 km² of basin area. The second record is the annual daily maximum precipitation (PMD) in millimeters corresponding to the Andong rainfall station in South Korea (^{Park et al., 2011}) and was selected for its strong upward trend.

The third record comes from ^{El Adlouni and Ouarda (2008b)} and they are randomly generated data with a GVE₂ distribution whose values in their fitting parameters were: δ₁ = 10, δ₂ = -0.10, δ₃ = 0.005, a = 1.0 and k = -0.10; corresponds to a PMD series with nonlinear trend. The fourth record also comes from the previous reference and has been used by ^{El Adlouni et al. (2007)}, ^{Aissaoui-Fqayeh et al. (2009}) and ^{Ouarda and El Adlouni (2011)}, corresponds to the PMD recorded in the period from 1952 to 2000 at the Tehachapi rainfall station in Southern California, U.S.A. These values are related to the SOI (Southern Oscillation Index) and, therefore, they will be used as a covariate.

General layout for probabilistic analysis

The graph of data values (X_i ) versus time (t_i ) allows defining if there is a linear or quadratic trend. When the trend is linear, the fitting parameters (δ₂, δ₁, a, k) of the non-stationary FDPs GVE₁, LOG₁ and PAG₁ are estimated, their EEA is quantified with equation 48, to select the FPD that leads to the lowest value of such indicator. When the trend is curved, the non-stationary PDFs the GVE₂, LOG₂ and PAG₂ are tested; for this, their fitting parameters are estimated (δ₂, δ₃, δ₁, a, k), their EEA is quantified, adopting the one that contributed the lowest value. In this process, a non-stationary FDP can be adopted, based on judgments of convenience, for example, the most unfavorable or critical predictions. The above is illustrated in the first numerical application.

In the analyses that use time (t) as a covariate, based on the inverse solutions (equations 15, 21 and 25) of the GVE, LOG and PAG distributions, predictions with return periods (Tr) of 2, 25, 50 and 100 years were calculated, through the registration period, the location parameter variable u was applied. The first prediction corresponds to the median, since its probability of non-exceedance (F) is 50 % and the following three are calculated for complementary probabilities, to define its superior and inferior value (^{Park et al., 2011}), that is, for the following values: F = 0.96 and F = 0.04 for the Tr of 25 years; F = 0.98 and F = 0.02 for the Tr of 50 years and F = 0.99 and F = 0.01 for the Tr of 100 years. Furthermore, in these analyses, predictions were made for the future, in the years 2020, 2050 and 2100. When the covariate is the SOI, the predictions corresponding to the Tr of 2, 25, 50 and 100 years, fluctuate according to the magnitude of the SOI and then its maximum values are indicated in the extreme magnitude of the SOI. In this case, future predictions are possible using prognostic SOI values.

Predictions at a Dartmouth River gauging station

These data show a linear descent and significant trend, since DS= -2.2931 and DSc = 2.0484 (Figure 1). In this record, the EEAs of the non-stationary FDPs GVE₁, LOG₁ and PAG₁ are 42.0, 43.1 and 41.6 m³·s^-1. Due to the numerical similarity that the EEA have, the results of the function LOG₁ are adopted because they lead to the largest or most critical predictions. Their fitting parameters are δ₂ = -3.125, δ₁ = 227.919, a = 33.318 and k = -0.144; with a linear correlation coefficient of -0.3976 and a slope of -2.8182 as Sen criterion. Table 2 shows some of the predictions within the record and in the future on three pre-fixed dates. This record of 30 data ends in 2003, therefore, the value of time t in 2020 is 47, in 2050 it is 77 and in 2100 it was 127.

Predictions at the Andong pluviometric station

This PMD record has a linear upward and significant trend, because DS = 2.1568 and DSc = 2.0452 (Figure 2). The fitting of the non-stationary PDFs GVE₁, LOG₁ and PAG₁ led to the following EEAs in millimeters: 9.8, 8.8 and 11.9. For this record the LOG₁ function is adopted, for achieving the lowest EEA, its fitting parameters were: δ₂ = 1.1101, δ₁ = 73.805, a = 13.740 and k = -0.113; with a linear correlation coefficient of 0.3718 and a slope of 1.100 according to the Sen criterion. Table 3, lists a part of the predictions within the record and in the future on three pre-fixed dates. As this record covers 31 data and concludes in 2007, the magnitude of time t in 2020 is 44, in 2050 it is 74 and in 2100 it was 124.

Predictions in a pluviometric fictitious station

This PMD record has an upward curve trend (Figure 3). The fitting of the non-stationary PDFs GVE₂, LOG₂ and PAG₂ led these values of the EEA in millimeters: 2.4, 2.3 and 2.5. The LOG₂ function is adopted to achieve the lowest EEA; its fitting parameters were: δ₂ = -0.3367, δ₃ = 0.0109, δ ₁ = 12.828, a = 1.618 and k = -0.264, with a polynomial correlation coefficient of 0.715. Table 4 shows a part of the predictions within the record and in future in three pre-fixed dates. This record of 50 data ends in the year 2010, then the magnitude of time t in 2020 is 60, in 2050 it is 90 and in 2010 it was 140.

Predictions in the pluviometric station Tehachapi

Figure 4 shows the annual PMD record is shown in the ordinates against their corresponding values of the SOI on the abscissa, since it is not possible to select between a linear descending trend or a curve, the two fits will be applied. The non-stationary PDFs GVE₁, LOG₁ and PAG₁ lead to the following values of the EEA in millimeters: 11.6, 11.4 and 12.2. For this linear fitting, the LOG₁ function is adopted for achieving the lowest EEA; its fitting parameters were: δ₂ = -10.2657, δ₁ = 28.376, a = 7.570 and k = -0.104, with a linear correlation coefficient of -0.5680. Table 5 shows a part of the prediction within the record.

The non-stationary FDPs GVE₂, LOG₂ and PAG₂ lead to the following values of the EEA in millimeters: 12.8, 12.7 and 13.3. For this non-linear fitting, the LOG₂ function is adopted, to achieve the lowest EEA value; its fitting parameters were: δ₂ = -8.7385, δ₃ = 2.8432, δ₁ = 25.936, a = 6.942 and k = -0.127, with a polynomial correlation coefficient of 0.6279. Table 6 shows a part of the predictions within the record.

The selection between the two non-stationary FDPs LOG₁ or LOG₂ should not be based on the EEA, but on the best fitting achieved with the X_i and SOI_i data pairs, which was achieved with a location parameter with quadratic variation with respect to the SOI, which defines a 0.6274 as a polynomial correlation coefficient. In addition, such fitting leads to higher predictions in the historical extreme value of the SOI of -3.2, relative to the year of 1986.