Data Availability: The data of the annual maximum flow series for selected gauging stations of Polish Rivers have been bought from the Institute of Meteorology and Water Management in Warsaw, Poland, and the authors are obligated not to publish them. Interested researchers can purchase these data by contacting The Expertise Department (historical data) at ekspertyzy@imgw.pl, +48 22 56 94 381 or +48 22 56 94 254. However, the statistical characteristic of these series and other relevant data are within the paper.

Funding: This work was partly financed by the grant of the Polish National Science Centre titled “Modern statistical models for analysis of flood frequency and features of flood waves“, decision nr DEC-2012/05/B/ST10/00482 and by the Polish Ministry of Science and Higher Education under the Grant Iuventus Plus IP 2010 024570 titled “Analysis of the efficiency of estimation methods in flood frequency modelling”. This work was partially supported within statutory activities No 3841/E-41/S/2015 of the Ministry of Science and Higher Education of Poland.

Competing interests: The authors have declared that no competing interests exist.

Introduction

Flood frequency analysis (FFA) provides information about the probable size of flood flows and has been used for the design of civil engineering works over the century. The assessment of the flood (upper) quantiles is required for dimensioning hydraulic structures affected by high waters, such as culverts, dams, bridges, overflow channels, spillways, levees, and others. FFA plays an important role in reducing the flood risk, since the flood quantile estimates are essential in determining the limits of flood zones with varying degree of flood risk as well in estimating the risk of exploitation of floodplains.

The quantile of the order of F ∈ (0,1) is defined as the value xF satisfying the equation: (1) where f is the probability density function (PDF) of the continuous random variable. The flood (upper) quantile means the probable maximum flow of the return period of T years and the relation between the probability of non-exceedance F and return period T has the form: (2) Since the return period T equal to 100, 500, 1000 is usually used, then the probability F is close to 1, i.e. close to its highest value. Equivalently, the probabilities of exceedance p can be applied, where: (3)

The at-site frequency analysis is the most commonly used approach. Then the estimation of flood quantiles refers to the choice of the form of probability density function describing the annual peak flows for the investigated gauging station. The distribution function assumed (also called the model) has a character of statistical hypothesis. Simultaneously, the method of estimation of parameters and, thus, upper quantiles of the assumed distribution is selected. This step is denoted D/E for “distribution and estimation procedure”. To find the best fitting model to the empirical data, the chosen discrimination procedure is applied.

The choice of distribution for fitting the annual maximum flows has attracted considerable interest, e.g. [1–7] and many others. According to the hydrological report of the World Meteorological Organization from 1989 [8], the Gumbel and log-normal distributions were the most commonly used for the description of the peak flow data. In Poland, the Pearson III type distribution has been recommended by Central Office of Water Management for national use [9]. These regulations are still in force, although other models are also applied in practice. Nowadays in many countries around the world, the heavy-tailed distributions are recommend for modelling of extreme flow series, e.g. [10–15]. The heavy-tailed distributions have conventional moments only in a certain range, which decreases with growing moment order. However, the heavy-tailed form of hydrological variables is not sufficiently supported, e.g. [16], [17]. Moreover, the analysis of Polish datasets of annual peak flows in [18] shows that they should be modeled using soft-tailed rather than heavy-tailed distributions.

The characteristics describing properties of the distribution are the summary statistics. Several systems of summary statistics have been developed. Based on different principles they provide, in particular, the measures of location, dispersion, skewness and kurtosis. The summary statistics calculated for a random sample consecutively serve for identifying and fitting PDFs. Among the systems of summary statistics, the most popular are the system of conventional moments (μr) and that of linear moments, called L-moments (λr), presented in Table 1 along with the dimensionless versions of the summary statistic sets in the form of summary statistic ratios (in square brackets). It is convenient to use the dimensionless versions of the summary statistics, since they measure the shape of a distribution independently of its scale of measurement.

As seen from Table 1, the L-moments can be defined by the probability weighted moments of a random variable βr for r = 0,1,2,… [19]. The L-moments create an attractive system because their estimators, in contrast to the classical moments estimators, are not biased and the sampling L-moment ratios have very small biases for moderate and large samples.

For the convenience of the reader, the abbreviations and symbols commonly used in the paper are gathered in Table in S1 Table.

For two-parameter distributions lower bounded at zero, a basic illustration that provides an intuition to a practitioner to distinguish various distributions is the graph of the relationship between the conventional variation coefficient CV and the conventional skewness coefficient CS or between their linear counterparts, i.e. between the linear variation coefficient LCV and the linear skewness coefficient LCS. These relationships show in what range of CV − CS various distributions can be used, e.g. the log-logistic and log-Gumbel distributions are not proper for modelling the data series of small skewness CS < 1 and average variation CV > 0.5. Both relations, CV − CS and LCV − LCS, are shown in Figs 1 and 2, respectively, for two-parameter distributions commonly used in FFA (lines) plotted with the Polish data of annual peak flows for 38 gauging stations (triangular points). To find the data availability, see the Acknowledgment.

In Figs 1 and 2, if some point lies on the line corresponding to certain distribution or around it, it may indicate that this distribution will be the best fitting to the data series. However, the perfect fit would hold for the asymptotic sample from a given distribution. Due to a limited length of the data series, such graphical analysis is only preliminary and the distribution best fitted to the data is indicated by the discrimination procedures, which will be discussed later in this paper.

The respective measuring sections are listed in Table 2 and illustrated in Fig 3. Most of analyzed data series cover the period 1921–2010.

As seen in Figs 1 and 2, for both conventional and linear moments ratios, there is a range of values taken by the Polish data series and not covered by any distribution. Clearly, there is still room for a new model. The inverse Gaussian (IG) and the generalized exponential (GE) distributions with the scale and shape parameters seem to be a suitable complement (see Figs 4 and 5), i.e. there are many points CV-CS corresponding to the Polish data series, which are on or around the lines of IG and GE distributions.

The GE distribution is used quite effectively to analyze lifetime data in the reliability analysis, being an alternative to the two-parameter gamma, Weibull, Pareto and log-normal distributions [20]. The aim of the study is to assess the usefulness of the generalized exponential distribution in flood frequency analysis for Polish Rivers, as a complementary to the inverse Gaussian distribution, which has proved to be suitable for many Polish data series [21], [22], [18]. In the paper, two-parameter distributions, instead of their three-parameter counterparts, are used for the modelling of relatively large-size samples (i.e. 90 elements), since our studies are intended to be applicable for most of the available observation series, which are much shorter than those investigated here. The short length of the data series hinders the proper selection of the distribution and two-parameter PDFs are usually used for their modelling. To reduce the uncertainty in the estimation of the extreme value distribution quantiles, the multi-model approach proposed by Bogdanowicz [23] is applied.

The paper is organized as follows. After providing some introduction to the topic, the probability distributions analyzed in the paper are shortly presented in second section. Next, the four discrimination procedures used to select the best fitting model are shown. Sequent two sections provide the results of the simulation studies on the probability of correct selection (PCS) among the GE and IG distributions along with the analysis of the asymptotic model error in respect to the upper quantile. In the case study section, fitting the GE and IG distributions to the 90-year series of annual maximum flows is compared for four selected gauging stations of Polish Rivers. Then, the method of aggregated quantiles is proposed for evaluation of upper quantile values. The paper is concluded in the final section.

GE and IG Probability Distributions

The inverse Gaussian (known also under the name of Wald) distribution has several properties analogous to the Gaussian distribution. In fact, the name is misleading, since it is “inverse” only in that, while the Gaussian describes the distribution of distance at a fixed time in Brownian motion, the inverse Gaussian describes PDF of the first passage time for a Brownian motion starting at zero to reach the absorbing barrier at the fixed point [24]. The same function appears in linear flood routing modelling as the impulse response of the semi-infinite channel at a fixed distance for the Froude number equal to zero [25], [26], and the name “linear convective diffusion model” for IG has been used in FFA [27–29]. In the last paper, the similarity between IG and LN distributions was shown by comparison of their first five moment estimates. Moreover, fitting of the two distributions to Polish data was compared there by the likelihood ratio. It indicates the preference of the IG model over the LN model for 27 out of 39 annual peak flow series. The simulation studies on the probability of correct selection among IG and LN have been carried out [21], adopting several discrimination statistics. The discrimination procedures based on the likelihood ratio and the R statistics [22] favor IG over LN, while the discrimination procedure based on the QK statistics [30] favors LN over IG. Investigation of Polish annual maxima datasets by the L-moment ratio diagrams and the test of linearity on log–log plots shows that the inverse Gaussian distribution represents flood frequency characteristics of Polish Rivers quite well, in particular of lowland rivers [18].

The generalized exponential distribution has been developed by Gupta and Kundu [31] and used quite effectively in many situations where a positive skewed distribution is needed. The closeness of GE distribution with gamma, Weibull, and log-normal distributions has been demonstrated [32–35]. The generalized exponential distribution has been applied to analyze lifetime data in the reliability analysis [20]. However, to the best of our knowledge it has not been used in FFA so far but in Poland where the GE model has been introduced for describing random properties of seasonal maximum annual flows [36].

The basic statistical characteristics of both IG and GE distributions are presented in Table 3.

The polygamma functions are defined as the logarithmic derivative of the gamma function [37]: (4) For real positive arguments z, digamma function ψ(z) is a concave increasing function of z which satisfies the following relation [37], [38]: (5) Differentiating Eq 5 appropriate number of times, one gets the evaluations of polygamma functions that can be used for numerical calculations instead of analytical formulas.

Only the first two linear moments of GE distribution have been derived so far [20]. The formula for its third linear moment (λ3) and thus for the linear skewness coefficient (LCS) has been derived by the authors (see Appendix) and presented in Table 3. Since the linear moments of IG distribution have no analytical form, their integral formulas are applied for computational calculations and the trapezoidal rule is used for approximation of the definite integral [39]. The details concerning the derivation of the formula for the quantile corresponding to probabilities of non-exceedance F (xF) for IG distribution are shown in [27].

As shown in Fig 4, for the variation coefficient CV equal to 0.41, the skewness coefficients CS of both GE and IG distributions are the same and amount to 1.23. As the two basic characteristics for the two-parameter distributions are equal, the shapes of distribution density functions are almost identical; see solid lines in Fig 6. However, the PDFs are not identical, since the values of kurtosis CK = μ4 / μ2 vary and equal to 5.67 and 3.68 for GE and IG distributions, respectively. As you move away from CV = 0.41, the differences in the values of CS of both distributions increase (Fig 4); therefore, the shapes of their density functions differ from each other. This is exemplified by CV = 0.8 and corresponds to the values CS = 1.69 and CS = 2.4 for GE and IG distributions, respectively; see solid lines in Fig 6.

Discrimination Procedures

The main disadvantage of using the wrong form of distribution for a flood series is that it over- or under-designs the hydraulic structures. Even if the sample size is not sufficiently large for making a proper choice among alternative distribution functions, a selection method is still needed and moreover all available information should be utilized for it. To find the best fitting model to empirical data from the set of competing models, a discrimination procedure is required. It must define a test statistics as well as a decision rule indicating the action to be taken for the sample under consideration. One can also prioritize all competing models according to the values of the selection criterion. However, the use of a discrimination procedure without the knowledge of its performance for the considered set of PDFs may be “a foolhardy gamble” [40] and may lead to erroneous conclusions [21]. To increase the efficiency of the model selection techniques in FFA, the use of several discrimination procedures along with the knowledge of their efficiency for a particular case is advisable.

K procedure

The K procedure [41], [42] of model selection is based on the likelihood functions for i = 1,…,k and k is the number of considered distributions expressed by their density functions fi. In fact, the K procedure is equivalent to the Akaike information criterion (AIC) [43] for distributions with the same number of parameters. The procedure points out the model with the highest value of the logarithm of the likelihood function as the true or the closest to the true model among all competing models, i.e.: (6) where is a set of distribution parameters evaluated by any estimation method. In this study, three methods of the assessment of parameters and, thus, of flood quantiles are applied, i.e. the method of moments (MOM) (e.g. [44]), the method of linear moments (LMM) (e.g. [19]), and the maximum likelihood method (MLM) (e.g. [45]). These methods were applied for the IG and GE distributions in [27] and [46], respectively. The accuracy of the estimates of large quantiles obtained from these three methods for the two- and three-parameter log-normal and GEV distributions have been analyzed in [47] both in the case of true and false hypothetical models, while the asymptotic bias of a quantile caused by the wrong distributional assumption has been analytically derived for a wide set of two-parameter distributions in [48], [49] and [17].

QK procedure

The QK discrimination procedure bases on the statistics that is invariant under scale transformation of the data [30]: (7) where N is the sample size and fi is the probability density function with scale parameter equal to one for k alternative models, i = 1,…,k. The unknown shape parameter of each of the considered distributions is estimated by the MLM method and substituted into Eq 7. As the selection rule, Quesenberry and Kent [50] proposed to choose the model which corresponds to the highest value of the Si statistics among competing PDFs. They showed that the QK discrimination procedure minimizes the sum of the probabilities of selecting the incorrect families of distribution. In practice the logarithm of the selection statistics Si instead of the statistics itself is usually applied: (8)

The analytical formula for the logarithms of the Si statistics of the inverse Gaussian distribution has been derived and published with small editorial error in [21]. Therefore, its corrected form is given below: (9) where Kν is the modified Bessel function of the second kind (e.g. [37]). Since we failed to get the analytical QK formula for the generalized exponential distribution, the selection statistics ln SGE has been calculated numerically from the definition Eq 7 using the trapezoidal rule for approximation of the definite integral (e.g. [39]).

R procedure

Since no simple statistical model can reproduce the dataset in its entire range of variability, it seems to be a right idea that the shape of the distribution tail should be a leading statistics when choosing a hypothetical distribution. Some guidelines and procedures for selecting the class of distributions that provides the best fit to the sample extremes are presented, for example, in [54–55]. However, there is a problem with a small number of data from the scope of the tail.

The R procedure follows the thought of fitting the model to the data in the range of the distribution tail. The parametric methods of the estimation of a model density function are asymptotically unbiased, which means that the assessment of any model parameter tends to its exact value for the sample withdrawn from the population of known distribution function. Then, in particular, the estimate of any quantile converges to its true value. Basing on this rule, the differences between the estimates obtained from various methods have been used to assess model fitting to the sample. The procedure of model discrimination has been explicitly proposed in [22] and based on the difference between 1% quantile assessment provided by the method of moments and the maximum likelihood method. The 1% quantile assessment is the most commonly used design value and corresponds to a probability of exceedance p = 0.01 (i.e. F = 0.99), expressed as a percentage. According to the relation between the probability and return period (Eqs 2 and 3), the determines the probable maximum flow which appears, on average, once in 100 years.

Here, two discrimination statistics, and , are proposed, for i = 1,…,k and k being the number of competing PDFs: (13)(14) where are the 1% quantiles estimated by the moment method, the linear moments method and the maximum likelihood method, respectively. For both discrimination statistics (Eqs 13 and 14), the best model is the one with their lowest value: (15)(16)

Note that for normal distribution the MOM and MLM methods are equivalent. Therefore, if the normal distribution is among the alternative distributions, it would be chosen by the R procedure. What's more, all three estimation methods also give the same estimate of the mean for gamma and the two-parameter inverse Gaussian distributions. This gives for these distributions the similarity of the estimates of quantiles for these three methods in the range of the main probability mass and, to a certain extent, for higher quantiles as well. Therefore, to use the R procedure properly, the knowledge of its performance for the considered set of PDFs is required.

Evaluation of Efficiency of Discrimination Procedures

Each discrimination procedure is considered to be of universal use, i.e. can be applied for model selection among any set of alternative PDFs, regardless of the sample size. However, the real drawback appears when for a small or medium sample size the discrimination procedures tend to favour some alternative distributions.

Discrimination between the generalized exponential and other two-parameter distributions has been already investigated in respect to the gamma [34], Weilbull [32] and log-normal [35] distributions. The ratio of the maximized likelihood functions has been used there to determine the probability of correct selection. Additionally, the selection among the We, LN and GE distributions has been studied in [56]. Here, the discrimination between the generalized exponential and the inverse Gaussian and vice versa is the subject of investigation. The efficiency of four procedures of discrimination has been evaluated using simulated data with GE as true (T) model and IG as a false (F) one and vice versa. S = 10,000 pseudo-random samples have been generated from the GE and IG distributions, respectively, for sample sizes N = 20, 50 and 100. To unify the distributions with respect to parameters, the original parameters were replaced by the mean μ and the variation coefficient CV. Without a loss of generality, the mean equal to one is assumed. The variation coefficient varies from 0.3 to 1.0; this covers the range of CV values for Polish data (see Fig 1). The results available in the literature [21], [22] indicate significant differences in the values of the PCS obtained by the K and QK discrimination procedures for different pairs of distributions and small sample sizes generated. A similar result was expected for the pair of distributions IG and GE. However, to our surprise, the values of PCS according to the K discrimination procedure (Fig 7) are almost identical to the values from the QK procedure (Fig 8). The differences between these two procedures are shown in Fig 9.

The values of the probability of inconsistent selection (PIS) in Fig 9 mean that for a single sample generated from the assumed PDF, one of the procedures, K or QK, points out the right PDF (correct selection), while, at the same time, the latter procedure points out the wrong PDF (incorrect selection). In other words, the values of PIS are the percentages of inconsistency of the two procedures.

We have detected the identity of the K and QK procedures for the pairs of the inverse Gaussian with log-normal or gamma distributions. This issue will be further investigated. However, it seems that the K and QK procedures of discrimination are equivalent when IG is one of the alternative distributions, which is a unique feature of this distribution.

Finally, the results obtained from the KS and both variants of R discrimination procedures, i.e. R1 and R2, are presented in Figs 10–12, respectively.

PCS for the pair of GE and IG distributions

It is quite clear from the above figures that the PCS increases with increasing sample size, i.e. the probability of correct selection is the smallest for 20-element samples and the highest for 100-element samples. The exception is the R2 procedure applied for the true GE distribution with IG as an alternative within the range of variation coefficient CV from 0.3 to 0.45 (Fig 12). It is also clear that as CV moves away from the value around 0.4, the PCS increases. For all considered discrimination procedures, if the variation coefficient is about 0.4, a sharp decline of the PCS value is visible. Only for the KS procedure the decrease of PCS value is more evident when IG is the true sample distribution (Fig 10), while for the other procedures this effect is stronger when GE is the true sample distribution (Figs 7, 8, 11 and 12). For example, for the K and QK procedures, when the data are drawn from the IG distribution, the PCS is about 63% at the minimum point for all the considered sizes of the sample, while if the data are drawn from the GE distribution, the PCS decreases up to about 34% for the sample size N = 20, i.e. then 66% of samples generated from GE distribution will be wrongly recognized as originated from IG parent distribution. The lowest value of the PCS among the GE and IG models for CV ∼ 0.4 is related with the fact that for the variation coefficient equal to 0.41, the skewness coefficients CS of both distributions are the same and amount to 1.23. Hence, for the range of CV around 0.4, the investigated distributions have a similar shape, as shown in Fig 6.

Similar results as for the procedures K and QK are obtained for the procedure R (Figs 11 and 12). However, note that for the variant R2 (Fig 12), the minimum of the PCS obtained for the case of T = GE and F = IG decreases below 30%. In general, if GE is the true sample distribution and CV is lower than 0.5, it does not make sense to use any variant of R discrimination procedure, since the probability of the correct selection of the distribution is lower than 50%. Then the decision based on “head and tail” rule is more efficient and easier to use. The same applies to the procedures K and QK and CV about 0.4, depending on the sample size N.

For all four discrimination procedures, the generalized exponential model is better recognizable than the inverse Gaussian, for moderate and large CV values, i.e. CV > 0.5, while for small CV values, i.e. CV < 0.5, the GE model is favoured only by KS procedure (Fig 10) and IG model is favoured by K, QK and R procedures (Figs 7, 8, 11 and 12).

For the range of CV from 0.6 to 0.8, which covers most data of Polish Rivers, the probability of correct selection among GE and IG distributions is quite large. Except of some cases of 20-element series drawn both from the IG and GE distributions, the PCS is higher than 70%. However, it should be remembered that the above experiment relates to a special theoretical case when one of the two competing distributions is the true one. If there are more alternative models, the probability of the selection of the true distribution significantly decreases. Similarly, the PCS is smaller if a set of alternative distributions consists of the PDFs which are similar in type to the true distribution. Note, the exponential distribution is a special case of the gamma, Weibull and generalized exponential models, if the shape parameter is equal to one, being a special case of Pareto, if the shape parameter is equal to zero. Then the PCS among any pair of the distributions from the set above would be lower than the PCS among GE and IG distributions.

Asymptotic Model Error in Respect to the Upper Quantile

In flood frequency analysis, the assumed (hypothetical) model is treated as the correct (true) model and any assessment of the accuracy of the estimation of its parameters and quantiles is usually made assuming that the considered random sample is derived from that probability distribution. In this way, the error of the choice of false distribution, i.e. the model error, is omitted, although this error can have a significant impact on the accuracy of the quantiles estimation. For a given estimation method, the total bias of quantile estimate consists of a sampling bias which asymptotically converges to zero and a model bias caused by wrong distributional selection. Those biases can be of opposite signs. The theoretical background for the asymptotic bias caused by false distributional assumption for various estimation methods has been presented in [48] followed by derivations for various pairs of (True, False) distributions in [49], [17].

Here, the set of competing PDFs involves the generalized exponential and inverse Gaussian distributions and the interest is in the derivation of the asymptotic model bias of the 1% quantile estimate; when the GE is the true (T) population model, then the IG is falsely (F) adopted for the hypothetical PDF, and vice versa. Three estimation methods presented in section 3.1 are used as approximation method of T distribution by F distribution. The relative model bias is defined for each approximation method as: (17)

To present a unified treatment for both distributions and three estimation methods, the mean of the T distribution equal to 1.0 and the variation coefficient CV varying from 0.3 to 1.0 are considered in the experiment. The relative asymptotic bias of is determined analytically by the methods of MOM and LMM, while to use MLM, the samples of size N