- Email: [email protected]
Using higher probability weighted moments for flood frequency analysis
Journal of Hydrology 194 (1997) 95–106
Using higher probability weighted moments for flood frequency analysis Q.J. Wang Department of Civil and Environmental Engineering, University of Melbourne, Parkville, Vic. 3052, Australia Received 27 February 1996; revised 10 June 1996; accepted 12 June 1996
Abstract When an annual maximum flow series is displayed in a probability plot, the data often exhibit two or more distinct segments. For estimating floods of large return periods, it is probably more appropriate that any interpolation and extrapolation are made by exploring the trend shown by the larger flows in the sample series. Lower bound censored samples may be used to fit distribution functions, but the analytical techniques available for doing so are either inadequate or complicated. In this paper, the use of higher probability weighted moments (PWMs) for distribution fitting is proposed. Higher PWMs give more weight to larger flows. By matching higher PWMs to their sample estimates, the resulting distribution fits more closely to the larger flows in the sample. Examples using actual data and Monte Carlo simulation results show that there is considerable merit in using the method for predicting large return period events. The method involves no more complication than the well-accepted method of using the lowest orders of PWMs.
1. Introduction Estimation of floods of large return periods is often required in engineering structure design. Given an annual maximum flow series, a distribution of assumed function form may be fitted to the data series and then used to calculate flood magnitudes of any return periods. This procedure becomes problematic when an incorrect distribution function form is used. This may be checked, at least within the record length, by displaying both the data and the fitted distribution on a probability plot. Whereas virtually all distribution functions used in practice are smooth in shape, the observed data often exhibit two or more segments which could not all fit into a single smooth curve but could each fit into one smooth curve. Forcing a single smooth function to fit all segments produces a compromised curve, which may lead to serious error in the estimation of floods of large return periods. When the purpose of flood frequency analysis is for estimating floods of large return 0022-1694/97/$17.00 q 1997– Elsevier Science B.V. All rights reserved PII S00 22-1694(96)032 23-4
96
Q.J. Wang/Journal of Hydrology 194 (1997) 95–106
periods, it is probably more appropriate that any interpolation and extrapolation from a record to the required return periods are made by following the trend shown by the larger floods in the record. This can be achieved graphically by drawing a curve following that trend, but different workers who apply the same technique may arrive at different results. To achieve more consistent results, an analytical technique should be employed. This paper examines briefly two existing analytical techniques, namely, the method of maximum likelihood (ML) and the method of partial probability weighted moments (partial PWMs). The method of using higher orders of (ordinary as distinct from partial) PWMs, instead of the lowest orders, to estimate flood distributions is proposed. Formulation of the method for the GEV distribution is given. Examples using actual data and Monte Carlo simulation results are presented to demonstrate the merit of the method.
2. The ML and partial PWM methods The ML method may be employed to fit a distribution function to a lower bound censored sample. With annual maximum flows, the censored sample is formed by discarding those flows which are below a certain threshold. The threshold may be chosen by examining the flow data displayed in a probability plot. The flows above that threshold should at least appear to fit into a single smooth curve. Given a sample of size n, the sample may be ranked in ascending order as x (1) # x (2) # … # x (n). Supposing a total of n 0 flows are below the chosen threshold x 0 a likelihood function may be constructed as L=
x0
−`
f (x)dx
n0
n
∏ f (x(i) ) = [F(x0 )]n0
i = n0 + 1
n
∏ f (x(i) )
i = n0 + 1
(1)
where f(x) is the probability density function and F(x) is the cumulative probability function. The unknown parameters in the distribution function are obtained by maximizing the likelihood function. The maximization may be carried out either by differentiating the likelihood function with respect to each of the unknown parameters in the distribution function and setting the derivatives to zero or by employing a direct numerical search technique. The likelihood function is relatively easy to set up and the maximization can be carried out by using one of the many available computer searching subroutines. However, the method often breaks down before convergence to an optimum is achieved. This is especially so when the distribution is lower or upper bounded. The likelihood function suddenly collapses as the bound places one or more of the sample values outside that bound. For the extreme value distributions, Leese (1973), Prescott and Walden (1983) and Phien and Fang (1989) have presented the details on the use of the method. To overcome the shortcomings of the ML method, Wang (1990a,b, 1996) introduced the method of partial PWMs to fit distribution functions to censored samples. The partial PWMs are defined as (Wang, 1996)
1
b9r
=
x(F)F r dF
F0
1 − F0r + 1
(2)
Q.J. Wang/Journal of Hydrology 194 (1997) 95–106
97
where r = 0,1,2,… denotes the order of the partial PWMs and F 0 is the probability that the threshold x 0 is not exceeded, i.e. F 0 = F(x 0). It should be noted that the above definition differs slightly from that given by Wang (1990a,b). When F 0 = 0, the partial PWMs become the ordinary PWMs. Theoretically speaking, an explicit or implicit analytical expression as a function of the distribution parameters can always be found for the partial PWMs given the distribution function form and F 0 value. Given a sample of length n, the sample may be ranked in ascending order as x (1) # x (2) # … # x (n). The partial PWMs can be estimated by (Wang, 1996) 1 1 n (i − 1)(i − 2)…(i − r) 9 (3) ∑ x(i) b9r = 1 − (n =n)r + 1 n i = 1 (n − 1)(n − 2)…(n − r) 0
where x9(i)
( =
0
if x(i) # x0
x(i)
if x(i) . x0
(4)
and n 0 is the number of events in the sample which do not exceed the threshold x 0. The parameters of the distribution can then be estimated by matching the first few orders of partial PWMs to their sample estimates. Wang (1990a, 1996) presented the formulation for fitting all three types of extreme value distributions. The method of partial PWMs is a useful mathematical tool but its application is likely to be hampered by the complexity involved in the formulation. In the case of the extreme value distributions, it is necessary to calculate values of some special functions (e.g. the Complete and Incomplete Gamma functions, the Exponential Integral function) and to solve a nonlinear equation. (However, there are readily available computer programs for doing so.) The estimation is also sensitive to the value of the threshold nonexceedance probability F 0, which can only be estimated. Moreover, the formulation of the method has so far been derived only for the three types of extreme value distributions.
3. The method of higher PWMS The above two analytical techniques are based on the use of lower bound censored samples to overcome the problem that smaller flows and larger flows in an annual maximum series do not appear to follow one smooth curve. A censored sample is formed by completely discarding those flows which are below a certain threshold. The flows above that threshold should at least appear to fit into a single smooth curve when displayed in a probability plot. By using a lower bound censored sample, floods of large return periods are estimated from the trend shown by the larger flows. In using a lower bound censored sample, the actual values of the smaller flows being censored have no influence on the fitting of the distribution function except that their total number of occurrences is accounted for. However, if the smaller flows are given some (but only very marginal) influence, the final outcome of the fitting will still depend very heavily on the larger flows. The definition of the PWMs makes the method of PWMs the immediate choice for implementing this idea.
98
Q.J. Wang/Journal of Hydrology 194 (1997) 95–106
The (ordinary as distinct from partial) PWMs are defined as (Greenwood et al., 1979) br =
1
x(F)F r dF
(5)
0
where r = 0,1,2,… denotes the order of the PWM. Theoretically speaking, an analytical expression as a function of the distribution parameters can be found for the PWMs for any given distribution function form. Given a ranked sample, x (1) # x (2) # … # x (n), the unbiased estimator of b r is given by (Landwehr et al., 1979a) 1 n (i − 1)(i − 2)…(i − r) (6) x br = ∑ n i = 1 (n − 1)(n − 2)…(n − r) (i) If the distribution function contains p unknown parameters, a total of p equations are needed to solve for them. These p equations may be set up by matching p PWMs of different orders to their sample estimates. All the current practices use the p lowest orders of PWMs, i.e. r = 0,1,…, p − 1. The PWM as defined by Eq. (5) gives more weight to larger x as the order of the PWM becomes higher. The nonexceedance probability F is a monotonic increasing function of x. The probability weight given to the variable x is F r. This power function increases with x at a much higher rate when the power r is higher, thus giving more weight to larger x. In other words, the value of a higher PWM depends more on the larger range of the variable. As discussed above, it is sometimes desirable to give the larger flows in a sample more influence than the smaller flows on the final outcome of the distribution function fitting. When this is so, the method of using higher PWMs should replace the commonly accepted method of using the lowest orders of PWMs. In fact, the orders of the PWMs to be used do not have to be very high because of the nature of the power function F r.
4. Fitting the GEV distribution The extreme value Type I, II and III distributions, also known collectively as the generalized extreme value (GEV) distribution, have the following expression:
8 9 > 1> < = k F(x) = exp − 1 − (x − y) k > > a : ;
1 = exp − exp − (x − y) a
if k Þ 0 (7)
if k = 0
or in inverse form a x(F) = y + [1 − ( − ln F)k ] if k Þ 0 k = y − aln( − ln F) if k = 0
(8)
where y is a location parameter, a is a scale parameter, and k is the important shape parameter. For k = 0, the distribution is unbounded both above and below, and is referred
99
Q.J. Wang/Journal of Hydrology 194 (1997) 95–106 Table 1
Coefficients a 1 and a 2 for eqn (17) and the maximum absolute error ldl in k resulting from using the equation to replace eqn (13) and eqn (16) for −0.5 # k # + 0.5 PWMs used r r r r r
= 0,1,2 = 1,2,3 = 2,3,4 = 3,4,5 = 4,5,6
a1
a2
ldl ,
7.8589 11.9082 15.9316 19.9455 23.9546
2.9534 2.7787 2.7301 2.7072 2.6936
8.8 × 10 −4 3.4 × 10 −4 1.8 × 10 −4 1.1 × 10 −4 7.7 × 10 −5
to as the extreme value Type I or Gumbel distribution. For k , 0, the distribution has a finite lower bound at y + a/k and a thick upper tail, and is referred to as the Type II extreme value distribution. For k . 0, the distribution has a finite upper bound at y + a/k, and is referred to as the Type III extreme value distribution. The PWMs for the GEV distribution have been derived by Hosking et al. (1985) for k Þ 0 as
a G(1 + k) (r + 1)br = y + 1 − k (r + 1)k
(9)
and by Greenwood et al. (1979) for k = 0 as (r + 1)br = y + a[e + ln(r + 1)]
(10)
where G() is the Gamma function, and e = 0.5772156649… is the Euler constant. As there are three unknown parameters, y, a and k, in the GEV distribution, three PWMs are required to set up three equations for solving for the parameters. Let us use the PWMs of orders r = h, h + 1 and h + 2. Writing the expressions for these PWMs using Eq. (9) and rearranging them yields for k Þ 0
(h + 1)bh = y +
a G(1 + k) 1− k (h + 1)k
(h + 2)bh + 1 − (h + 1)bh = a
G(1 + k) [(h + 1) − k − (h + 2) − k ] k
(h + 2)bh + 1 − (h + 1)bh (h + 1) − k − (h + 2) − k = (h + 3)bh + 2 − (h + 1)bh (h + 1) − k − (h + 3) − k
(11)
(12)
(13)
Doing the same for k = 0 using Eq. (10) gives (h + 1)bh = y + a[e + ln(h + 1)]
(14)
(h + 2)bh + 1 − (h + 1)bh = a[ln(h + 2) − ln(h + 1)]
(15)
(h + 2)bh + 1 − (h + 1)bh ln(h + 2) − ln(h + 1) = (h + 3)bh + 2 − (h + 1)bh ln(h + 3) − ln(h + 1)
(16)
100
Q.J. Wang/Journal of Hydrology 194 (1997) 95–106
By replacing bh , bh + 1 and bh + 2 , by their estimates bh , bh + 1 , and b h+2 computed using Eq. (6), the three parameters y, a and k may be found by solving the above equations. Eq. (13) does not give an explicit solution for k and has to be solved numerically using an iterative method. This can be avoided, however, by using the following approximation: k = a 1 c + a2 c 2
(17)
where c=
(h + 2)bh + 1 − (h + 1)bh ln(h + 2) − ln(h + 1) − (h + 3)bh + 2 − (h + 1)bh ln(h + 3) − ln(h + 1)
(18)
The values of coefficients a 1 and a 2 vary with the orders of the PWMs used. Table 1 gives a list of values and the error in k resulting from the use of the approximation within the range −0.5 # k # + 0.5. It should be noted that the variable c has been so defined as in Eq. (18) that Eq. (17) has a zero intercept to satisfy Eq. (16). Once k is obtained, Eq. (12) or Eq. (15) directly gives a solution for a and Eq. (11) or Eq. (14) for y. The Gamma function G() in Eq. (11) and Eq. (12) may be calculated using one of the approximation algorithms. Press et al. (1992) provided a subroutine gammln for calculating the logarithm of the Gamma function. Davis (1965) also listed a number of approximation formulae for calculating the Gamma function. It is clear that the estimation of a probability distribution function by using higher PWMs involves no more complication than using the lowest PWMs. The method is readily applicable to other distributions for which expressions of PWMs have been derived.
5. Examples Three examples are used here to demonstrate the effect of using higher PWMs in estimating the GEV distribution for flood frequency analysis. As a first example, the annual maximum flows at the Bromfleet station of the Albert river in Queensland, Australia, are analysed. The station has a catchment area of 544 km 2. Annual maximum flows for the 50 years from 1918 to 1968 are used, excluding 1926 for which no data are available. The data are plotted as filled circles against the Gumbel reduced variable in the probability plot of Fig. 1. The Gumbel (or extreme value Type I) reduced variable, g, is calculated as g = − ln( − ln F)
(19)
The nonexceedance probability for the ith value in the sample ranked in ascending order is calculated using the Gringorten (1963) formula F(i) =
i − 0:44 n + 0:12
(20)
where n is the sample size and i = 1,2,…,n. Eq. (19) is the Gumbel distribution, in reverse form, of Eq. (8) with k = 0, y = 0, a = 1. When plotted against the Gumbel reduced variable, a sample from the GEV distribution would fall roughly along a straight line if k = 0, a
Q.J. Wang/Journal of Hydrology 194 (1997) 95–106
101
Fig. 1. Fitting of the GEV distribution to annual maximum flows at Bromfleet station of the Albert river, Queensland, Australia.
concave curve if k , 0 and a convex curve if k . 0, corresponding to extreme value distributions Type I, II and III, respectively. The annual maximum flow data shown as filled circles in Fig. 1 appear to form at least two segments, one consisting of the lower flows and having a concave shape, and the other consisting of higher flows and having a convex shape. All the annual maximum flows are used to fit the GEV distribution first using PWMs of r = 0, 1 and 2. The estimated parameter values for y, a and k are 276 m 3 s −1, 353 m 3 s −1 and −0.170, respectively. The fitted distribution is plotted in Fig. 1 as the continuous line. Clearly, the fitting is in serious error, especially for the larger flows. The fitting, however, improves greatly when PWMs of r = 2, 3 and 4 are used. The larger flows are much better fitted by the resulting distribution, shown in Fig. 1 as the dotted line. The estimated parameter values for y, a and k are now 179 m 3 s −1, 691 m 3 s −1 and 0.258, respectively. It is clear that the two lines will give very different estimates of floods of large return periods. The difference will have considerable implication on project design and cost. A second example also demonstrates the advantage of using higher PWMs. The data are from the Wickliffe station, with a catchment area of 1350 km 2, of the Hopkins river, Victoria, Australia. A total of 52 annual maximum flows during the period 1921–1992, excluding the years with missing data, are used. The data are plotted against the Gumbel reduced variable in Fig. 2, showing again that the smaller flows behave differently from the rest. Using PWMs of the lowest three orders r = 0, 1 and 2 to fit the GEV distribution to all the data results in a compromised curve, as shown by the continuous line in Fig. 2. Using PWMs of higher orders r = 2, 3 and 4, on the other hand, produces a good fit to the larger flows, as shown by the dotted line in Fig. 2. The estimated parameter values for y, a and k are respectively 23.2 m 3 s −1, 25.1 m 3 s −1 and 0.061 in the former case, and 21.1 m 3 s −1, 37.2 m 3 s −1 and 0.369 in the latter. A third example is used to show that if the annual maximum flows in a sample can all be reasonably described by a single GEV distribution function, then using PWMs of various
102
Q.J. Wang/Journal of Hydrology 194 (1997) 95–106
Fig. 2. Fitting of the GEV distribution to annual maximum flows at Wickliffe station of the Hopkins river, Victoria, Australia.
orders will produce only minimum change in the resulting distribution curve. The data are from the Savage’s Crossing station, with a catchment area of 10 334 km 2, of the Brisbane river, Queensland, Australia. A total of 60 annual maximum flows during the period 1910– 1969 are available. The data are plotted against the Gumbel reduced variable in Fig. 3. The fitted distributions using PWMs of r = 0, 1 and 2 and PWMs of r = 2, 3 and 4 are also shown in Fig. 3. The difference between the two curves is small, both describing the data reasonably well. The estimated parameter values for y, a and k are respectively 527 m 3 s −1, 579 m 3 s −1 and −0.305 in the former case, and 483 m 3 s −1, 689 m 3 s −1 and −0.224 in the latter.
6. Monte Carlo simulations Monte Carlo simulations have been conducted to study how quantile estimation is affected by using higher PWMs. In practice, the true underlying distribution is never known and often differs from the distribution function actually used. For this reason, the study is carried out by fitting the GEV distribution function to generated Wakeby samples. Houghton (1978) introduced the Wakeby distribution function for flood frequency analysis. It contains five parameters and is thus highly flexible. By varying the parameter values, Landwehr et al. (1979b,c, 1980) constructed six Wakeby distributions, identified as WA-1, WA-2, WA-3, WA-4, WA-5 and WA-6, to represent a wide range of skewness and kurtosis. They then used the six Wakeby distributions as parent distributions in Monte Carlo simulations to assess the performance of various estimation methods and various assumed distribution function forms. The same six Wakeby distributions are used in this study as parent distributions to
Q.J. Wang/Journal of Hydrology 194 (1997) 95–106
103
Fig. 3. Fitting of the GEV distribution to annual maximum flows at Savage’s Crossing station of the Brisbane river, Queensland, Australia.
generate random samples. A total of 100 000 samples of size 30 are used in each case. The samples are fitted to the GEV distribution function using first PWMs of r = 0, 1 and 2 and then PWMs of r = 2, 3 and 4. Quantile estimates are obtained to calculate their bias and mean square error. Fig. 4 and Fig. 5 show respectively bias in x(F = 0.98) and x(F = 0.99) estimates for the first five Wakeby distributions. Using PWMs of r = 2, 3 and 4 almost always leads to much smaller bias than using PWMs of r = 0, 1 and 2 for all five Wakeby distributions. Mean square error has been obtained also for quantiles x(F = 0.98) and x(F0.99). The results for the first five Wakeby distributions are presented in Fig. 6 in terms of relative
Fig. 4. Bias in x(F = 0.98) estimates from fitting the GEV distribution to the Wakeby samples.
104
Q.J. Wang/Journal of Hydrology 194 (1997) 95–106
Fig. 5. Bias in x(F = 0.99) estimates from fitting the GEV distribution to the Wakeby samples.
estimation efficiency, defined as f=
mean square error in estimate using r = 0, 1 and 2 mean square error in estimate using r = 2, 3 and 4
(21)
Fig. 6 shows that the relative estimation efficiency is generally greater than one, implying that using PWMs of r = 2, 3 and 4 almost always leads to much smaller variance than using PWMs of r = 0, 1 and 2 for all five Wakeby distributions. For the sixth Wakeby distribution, WA-6, of Landwehr et al. (1979b,c, 1980), using PWMs of r = 2, 3 and 4 rather than PWMs of r = 0, 1 and 2 results in a significant decrease in efficiency for quantile estimation. It should, however, be pointed out that the WA-6 distribution has a zero product–moment skewness and is not a typical distribution for describing annual maximum streamflows.
Fig. 6. Relative efficiency of x(F = 0.98) and x(F = 0.99) estimates from fitting the GEV distribution to the Wakeby samples.
Q.J. Wang/Journal of Hydrology 194 (1997) 95–106
105
7. Conclusions The use of higher PWMs for estimating flood frequency distributions has been proposed and examined. Higher PWMs give more weight to the larger flows. A distribution estimated by matching higher PWMs to their sample estimates is more closely fitted to the larger flows. Using a distribution so fitted, floods of large return periods are estimated from the trend shown by the larger flows in the recorded sample. This overcomes the problem that the smaller flows may exhibit different trend from the larger flows in a sample, and when a single smooth distribution function, such as the GEV distribution, is forced to fit all the data in such a sample, the fitted curve may have to take up a compromised position and shape, often resulting in poor fit to the larger flows in the sample and serious error in the estimation of floods of large return periods. The examples used in the paper demonstrate the merit of using higher PWMs to fit distributions for estimating floods of large return periods. The orders of the PWMs to be used do not have to be very high. It appears that orders r = 2 and above are sufficient. However, more study is needed to determine the most appropriate orders to be used. As a principle, there is a need to strike a balance between extracting from samples important information on the distribution tail and avoiding the estimation becoming too sensitive to sample outliers. The third example shows that if all the annual maximum flows in a sample can be reasonably described by a single GEV distribution function, then using PWMs of various orders produces only minimum change in the resulting distribution curve. Monte Carlo simulation results from fitting the GEV distribution to samples generated from the first five of the six Wakeby distributions of Landwehr et al. (1979b,c, 1980) show that using PWMs of order r = 2, 3 and 4 instead of PWMs of order r = 0, 1 and 2 almost consistently leads to improved quantile estimation in terms of both bias and mean square error. The method of using higher PWMs involves no more complication than the commonly accepted method of using the lowest PWMs, and the method is readily applicable to all distributions for which expressions of PWMs have been derived. It is suggested that higher PWMs be used as a matter of routine for estimating floods of large return periods.
Acknowledgements The author would like to express his thanks to Professor Tom McMahon and Murray Peel, CRC for Catchment Hydrology, University of Melbourne, for their discussions and comments on the paper.
References Davis, P.J., 1965. Gamma functions and related functions. In: M. Abramowitz and I.A. Stegun (Editors), Handbook of Mathematical Functions. Dover, New York, pp. 253–293. Greenwood, J.A., Landwehr, J.M., Matalas, N.C. and Wallis, J.R., 1979. Probability weighted moments: definition and relation to parameters of distribution expressible in inverse form. Water Resour. Res., 15(5): 1049– 1054.
106
Q.J. Wang/Journal of Hydrology 194 (1997) 95–106
Gringorten, I.I., 1963. A plotting rule for extreme probability paper. J. Geophys. Res., 68(3): 813–814. Hosking, J.R.M., Wallis, J.R. and Wood, E.F., 1985. Estimation of the generalized extreme value distribution by the method of probability weighted moments. Technometrics, 27(3): 251–261. Houghton, J.C., 1978. Birth of a parent: the Wakeby distribution for modeling flood flows. Water Resour. Res., 14(6): 1105–1110. Landwehr, J.M., Matalas, N.C. and Wallis, J.R., 1979a. Probability weighted moments compared with some traditional techniques in estimating Gumbel parameters and quantiles. Water Resour. Res., 15(5): 1055– 1064. Landwehr, J.M., Matalas, N.C. and Wallis, J.R., 1979b. Estimation of parameters and quantiles of Wakeby distributions. 1. Known lower bounds. Water Resour. Res., 15(6): 1361–1372. Landwehr, J.M., Matalas, N.C. and Wallis, J.R., 1979c. Estimation of parameters and quantiles of Wakeby distributions. 2. Unknown lower bounds. Water Resour. Res., 15(6): 1373–1379. Landwehr, J.M., Matalas, N.C. and Wallis, J.R., 1980. Quantile estimation with more or less floodlike distributions. Water Resour. Res., 16(3): 547–555. Leese, M.N., 1973. Use of censored data in the estimation of Gumbel distribution parameters for annual maximum flood series. Water Resour. Res., 9(6): 1534–1542. Phien, H.N. and Fang, T.S.E., 1989. Maximum likelihood estimation of the parameters and quantiles of the generalized extreme-value distribution from censored samples. J. Hydrol., 105: 139–155. Prescott, P. and Walden, A.T., 1983. Maximum likelihood estimation of the parameters of the three-parameter generalized extreme-value distribution from censored samples. J. Stat. Comput. Simul., 16: 241–250. Press, W.H., Teukolsky, S.A., Wetterling, W.T. and Flannery, B.P., 1992. Numerical Recipes in FORTRAN: the Art of Scientific Computing, 2nd edn. Cambridge University Press, New York, 963 pp. Wang, Q.J., 1990a. Estimation of the GEV distribution from censored samples by method of partial probability weighted moments. J. Hydrol., 120: 103–114. Wang, Q.J., 1990b. Unbiased estimation of probability weighted moments and partial probability weighted moments from systematic and historical flood information and their application to estimating the GEV distribution. J. Hydrol., 120: 115–124. Wang, Q.J., 1996. Using partial probability weighted moments to fit the extreme value distributions to censored samples. Water Resour. Res., 32(6): 1767–1771.
Using higher probability weighted moments for flood frequency analysis Q.J. Wang Department of Civil and Environmental Engineering, University of Melbourne, Parkville, Vic. 3052, Australia Received 27 February 1996; revised 10 June 1996; accepted 12 June 1996
Abstract When an annual maximum flow series is displayed in a probability plot, the data often exhibit two or more distinct segments. For estimating floods of large return periods, it is probably more appropriate that any interpolation and extrapolation are made by exploring the trend shown by the larger flows in the sample series. Lower bound censored samples may be used to fit distribution functions, but the analytical techniques available for doing so are either inadequate or complicated. In this paper, the use of higher probability weighted moments (PWMs) for distribution fitting is proposed. Higher PWMs give more weight to larger flows. By matching higher PWMs to their sample estimates, the resulting distribution fits more closely to the larger flows in the sample. Examples using actual data and Monte Carlo simulation results show that there is considerable merit in using the method for predicting large return period events. The method involves no more complication than the well-accepted method of using the lowest orders of PWMs.
1. Introduction Estimation of floods of large return periods is often required in engineering structure design. Given an annual maximum flow series, a distribution of assumed function form may be fitted to the data series and then used to calculate flood magnitudes of any return periods. This procedure becomes problematic when an incorrect distribution function form is used. This may be checked, at least within the record length, by displaying both the data and the fitted distribution on a probability plot. Whereas virtually all distribution functions used in practice are smooth in shape, the observed data often exhibit two or more segments which could not all fit into a single smooth curve but could each fit into one smooth curve. Forcing a single smooth function to fit all segments produces a compromised curve, which may lead to serious error in the estimation of floods of large return periods. When the purpose of flood frequency analysis is for estimating floods of large return 0022-1694/97/$17.00 q 1997– Elsevier Science B.V. All rights reserved PII S00 22-1694(96)032 23-4
96
Q.J. Wang/Journal of Hydrology 194 (1997) 95–106
periods, it is probably more appropriate that any interpolation and extrapolation from a record to the required return periods are made by following the trend shown by the larger floods in the record. This can be achieved graphically by drawing a curve following that trend, but different workers who apply the same technique may arrive at different results. To achieve more consistent results, an analytical technique should be employed. This paper examines briefly two existing analytical techniques, namely, the method of maximum likelihood (ML) and the method of partial probability weighted moments (partial PWMs). The method of using higher orders of (ordinary as distinct from partial) PWMs, instead of the lowest orders, to estimate flood distributions is proposed. Formulation of the method for the GEV distribution is given. Examples using actual data and Monte Carlo simulation results are presented to demonstrate the merit of the method.
2. The ML and partial PWM methods The ML method may be employed to fit a distribution function to a lower bound censored sample. With annual maximum flows, the censored sample is formed by discarding those flows which are below a certain threshold. The threshold may be chosen by examining the flow data displayed in a probability plot. The flows above that threshold should at least appear to fit into a single smooth curve. Given a sample of size n, the sample may be ranked in ascending order as x (1) # x (2) # … # x (n). Supposing a total of n 0 flows are below the chosen threshold x 0 a likelihood function may be constructed as L=
x0
−`
f (x)dx
n0
n
∏ f (x(i) ) = [F(x0 )]n0
i = n0 + 1
n
∏ f (x(i) )
i = n0 + 1
(1)
where f(x) is the probability density function and F(x) is the cumulative probability function. The unknown parameters in the distribution function are obtained by maximizing the likelihood function. The maximization may be carried out either by differentiating the likelihood function with respect to each of the unknown parameters in the distribution function and setting the derivatives to zero or by employing a direct numerical search technique. The likelihood function is relatively easy to set up and the maximization can be carried out by using one of the many available computer searching subroutines. However, the method often breaks down before convergence to an optimum is achieved. This is especially so when the distribution is lower or upper bounded. The likelihood function suddenly collapses as the bound places one or more of the sample values outside that bound. For the extreme value distributions, Leese (1973), Prescott and Walden (1983) and Phien and Fang (1989) have presented the details on the use of the method. To overcome the shortcomings of the ML method, Wang (1990a,b, 1996) introduced the method of partial PWMs to fit distribution functions to censored samples. The partial PWMs are defined as (Wang, 1996)
1
b9r
=
x(F)F r dF
F0
1 − F0r + 1
(2)
Q.J. Wang/Journal of Hydrology 194 (1997) 95–106
97
where r = 0,1,2,… denotes the order of the partial PWMs and F 0 is the probability that the threshold x 0 is not exceeded, i.e. F 0 = F(x 0). It should be noted that the above definition differs slightly from that given by Wang (1990a,b). When F 0 = 0, the partial PWMs become the ordinary PWMs. Theoretically speaking, an explicit or implicit analytical expression as a function of the distribution parameters can always be found for the partial PWMs given the distribution function form and F 0 value. Given a sample of length n, the sample may be ranked in ascending order as x (1) # x (2) # … # x (n). The partial PWMs can be estimated by (Wang, 1996) 1 1 n (i − 1)(i − 2)…(i − r) 9 (3) ∑ x(i) b9r = 1 − (n =n)r + 1 n i = 1 (n − 1)(n − 2)…(n − r) 0
where x9(i)
( =
0
if x(i) # x0
x(i)
if x(i) . x0
(4)
and n 0 is the number of events in the sample which do not exceed the threshold x 0. The parameters of the distribution can then be estimated by matching the first few orders of partial PWMs to their sample estimates. Wang (1990a, 1996) presented the formulation for fitting all three types of extreme value distributions. The method of partial PWMs is a useful mathematical tool but its application is likely to be hampered by the complexity involved in the formulation. In the case of the extreme value distributions, it is necessary to calculate values of some special functions (e.g. the Complete and Incomplete Gamma functions, the Exponential Integral function) and to solve a nonlinear equation. (However, there are readily available computer programs for doing so.) The estimation is also sensitive to the value of the threshold nonexceedance probability F 0, which can only be estimated. Moreover, the formulation of the method has so far been derived only for the three types of extreme value distributions.
3. The method of higher PWMS The above two analytical techniques are based on the use of lower bound censored samples to overcome the problem that smaller flows and larger flows in an annual maximum series do not appear to follow one smooth curve. A censored sample is formed by completely discarding those flows which are below a certain threshold. The flows above that threshold should at least appear to fit into a single smooth curve when displayed in a probability plot. By using a lower bound censored sample, floods of large return periods are estimated from the trend shown by the larger flows. In using a lower bound censored sample, the actual values of the smaller flows being censored have no influence on the fitting of the distribution function except that their total number of occurrences is accounted for. However, if the smaller flows are given some (but only very marginal) influence, the final outcome of the fitting will still depend very heavily on the larger flows. The definition of the PWMs makes the method of PWMs the immediate choice for implementing this idea.
98
Q.J. Wang/Journal of Hydrology 194 (1997) 95–106
The (ordinary as distinct from partial) PWMs are defined as (Greenwood et al., 1979) br =
1
x(F)F r dF
(5)
0
where r = 0,1,2,… denotes the order of the PWM. Theoretically speaking, an analytical expression as a function of the distribution parameters can be found for the PWMs for any given distribution function form. Given a ranked sample, x (1) # x (2) # … # x (n), the unbiased estimator of b r is given by (Landwehr et al., 1979a) 1 n (i − 1)(i − 2)…(i − r) (6) x br = ∑ n i = 1 (n − 1)(n − 2)…(n − r) (i) If the distribution function contains p unknown parameters, a total of p equations are needed to solve for them. These p equations may be set up by matching p PWMs of different orders to their sample estimates. All the current practices use the p lowest orders of PWMs, i.e. r = 0,1,…, p − 1. The PWM as defined by Eq. (5) gives more weight to larger x as the order of the PWM becomes higher. The nonexceedance probability F is a monotonic increasing function of x. The probability weight given to the variable x is F r. This power function increases with x at a much higher rate when the power r is higher, thus giving more weight to larger x. In other words, the value of a higher PWM depends more on the larger range of the variable. As discussed above, it is sometimes desirable to give the larger flows in a sample more influence than the smaller flows on the final outcome of the distribution function fitting. When this is so, the method of using higher PWMs should replace the commonly accepted method of using the lowest orders of PWMs. In fact, the orders of the PWMs to be used do not have to be very high because of the nature of the power function F r.
4. Fitting the GEV distribution The extreme value Type I, II and III distributions, also known collectively as the generalized extreme value (GEV) distribution, have the following expression:
8 9 > 1> < = k F(x) = exp − 1 − (x − y) k > > a : ;
1 = exp − exp − (x − y) a
if k Þ 0 (7)
if k = 0
or in inverse form a x(F) = y + [1 − ( − ln F)k ] if k Þ 0 k = y − aln( − ln F) if k = 0
(8)
where y is a location parameter, a is a scale parameter, and k is the important shape parameter. For k = 0, the distribution is unbounded both above and below, and is referred
99
Q.J. Wang/Journal of Hydrology 194 (1997) 95–106 Table 1
Coefficients a 1 and a 2 for eqn (17) and the maximum absolute error ldl in k resulting from using the equation to replace eqn (13) and eqn (16) for −0.5 # k # + 0.5 PWMs used r r r r r
= 0,1,2 = 1,2,3 = 2,3,4 = 3,4,5 = 4,5,6
a1
a2
ldl ,
7.8589 11.9082 15.9316 19.9455 23.9546
2.9534 2.7787 2.7301 2.7072 2.6936
8.8 × 10 −4 3.4 × 10 −4 1.8 × 10 −4 1.1 × 10 −4 7.7 × 10 −5
to as the extreme value Type I or Gumbel distribution. For k , 0, the distribution has a finite lower bound at y + a/k and a thick upper tail, and is referred to as the Type II extreme value distribution. For k . 0, the distribution has a finite upper bound at y + a/k, and is referred to as the Type III extreme value distribution. The PWMs for the GEV distribution have been derived by Hosking et al. (1985) for k Þ 0 as
a G(1 + k) (r + 1)br = y + 1 − k (r + 1)k
(9)
and by Greenwood et al. (1979) for k = 0 as (r + 1)br = y + a[e + ln(r + 1)]
(10)
where G() is the Gamma function, and e = 0.5772156649… is the Euler constant. As there are three unknown parameters, y, a and k, in the GEV distribution, three PWMs are required to set up three equations for solving for the parameters. Let us use the PWMs of orders r = h, h + 1 and h + 2. Writing the expressions for these PWMs using Eq. (9) and rearranging them yields for k Þ 0
(h + 1)bh = y +
a G(1 + k) 1− k (h + 1)k
(h + 2)bh + 1 − (h + 1)bh = a
G(1 + k) [(h + 1) − k − (h + 2) − k ] k
(h + 2)bh + 1 − (h + 1)bh (h + 1) − k − (h + 2) − k = (h + 3)bh + 2 − (h + 1)bh (h + 1) − k − (h + 3) − k
(11)
(12)
(13)
Doing the same for k = 0 using Eq. (10) gives (h + 1)bh = y + a[e + ln(h + 1)]
(14)
(h + 2)bh + 1 − (h + 1)bh = a[ln(h + 2) − ln(h + 1)]
(15)
(h + 2)bh + 1 − (h + 1)bh ln(h + 2) − ln(h + 1) = (h + 3)bh + 2 − (h + 1)bh ln(h + 3) − ln(h + 1)
(16)
100
Q.J. Wang/Journal of Hydrology 194 (1997) 95–106
By replacing bh , bh + 1 and bh + 2 , by their estimates bh , bh + 1 , and b h+2 computed using Eq. (6), the three parameters y, a and k may be found by solving the above equations. Eq. (13) does not give an explicit solution for k and has to be solved numerically using an iterative method. This can be avoided, however, by using the following approximation: k = a 1 c + a2 c 2
(17)
where c=
(h + 2)bh + 1 − (h + 1)bh ln(h + 2) − ln(h + 1) − (h + 3)bh + 2 − (h + 1)bh ln(h + 3) − ln(h + 1)
(18)
The values of coefficients a 1 and a 2 vary with the orders of the PWMs used. Table 1 gives a list of values and the error in k resulting from the use of the approximation within the range −0.5 # k # + 0.5. It should be noted that the variable c has been so defined as in Eq. (18) that Eq. (17) has a zero intercept to satisfy Eq. (16). Once k is obtained, Eq. (12) or Eq. (15) directly gives a solution for a and Eq. (11) or Eq. (14) for y. The Gamma function G() in Eq. (11) and Eq. (12) may be calculated using one of the approximation algorithms. Press et al. (1992) provided a subroutine gammln for calculating the logarithm of the Gamma function. Davis (1965) also listed a number of approximation formulae for calculating the Gamma function. It is clear that the estimation of a probability distribution function by using higher PWMs involves no more complication than using the lowest PWMs. The method is readily applicable to other distributions for which expressions of PWMs have been derived.
5. Examples Three examples are used here to demonstrate the effect of using higher PWMs in estimating the GEV distribution for flood frequency analysis. As a first example, the annual maximum flows at the Bromfleet station of the Albert river in Queensland, Australia, are analysed. The station has a catchment area of 544 km 2. Annual maximum flows for the 50 years from 1918 to 1968 are used, excluding 1926 for which no data are available. The data are plotted as filled circles against the Gumbel reduced variable in the probability plot of Fig. 1. The Gumbel (or extreme value Type I) reduced variable, g, is calculated as g = − ln( − ln F)
(19)
The nonexceedance probability for the ith value in the sample ranked in ascending order is calculated using the Gringorten (1963) formula F(i) =
i − 0:44 n + 0:12
(20)
where n is the sample size and i = 1,2,…,n. Eq. (19) is the Gumbel distribution, in reverse form, of Eq. (8) with k = 0, y = 0, a = 1. When plotted against the Gumbel reduced variable, a sample from the GEV distribution would fall roughly along a straight line if k = 0, a
Q.J. Wang/Journal of Hydrology 194 (1997) 95–106
101
Fig. 1. Fitting of the GEV distribution to annual maximum flows at Bromfleet station of the Albert river, Queensland, Australia.
concave curve if k , 0 and a convex curve if k . 0, corresponding to extreme value distributions Type I, II and III, respectively. The annual maximum flow data shown as filled circles in Fig. 1 appear to form at least two segments, one consisting of the lower flows and having a concave shape, and the other consisting of higher flows and having a convex shape. All the annual maximum flows are used to fit the GEV distribution first using PWMs of r = 0, 1 and 2. The estimated parameter values for y, a and k are 276 m 3 s −1, 353 m 3 s −1 and −0.170, respectively. The fitted distribution is plotted in Fig. 1 as the continuous line. Clearly, the fitting is in serious error, especially for the larger flows. The fitting, however, improves greatly when PWMs of r = 2, 3 and 4 are used. The larger flows are much better fitted by the resulting distribution, shown in Fig. 1 as the dotted line. The estimated parameter values for y, a and k are now 179 m 3 s −1, 691 m 3 s −1 and 0.258, respectively. It is clear that the two lines will give very different estimates of floods of large return periods. The difference will have considerable implication on project design and cost. A second example also demonstrates the advantage of using higher PWMs. The data are from the Wickliffe station, with a catchment area of 1350 km 2, of the Hopkins river, Victoria, Australia. A total of 52 annual maximum flows during the period 1921–1992, excluding the years with missing data, are used. The data are plotted against the Gumbel reduced variable in Fig. 2, showing again that the smaller flows behave differently from the rest. Using PWMs of the lowest three orders r = 0, 1 and 2 to fit the GEV distribution to all the data results in a compromised curve, as shown by the continuous line in Fig. 2. Using PWMs of higher orders r = 2, 3 and 4, on the other hand, produces a good fit to the larger flows, as shown by the dotted line in Fig. 2. The estimated parameter values for y, a and k are respectively 23.2 m 3 s −1, 25.1 m 3 s −1 and 0.061 in the former case, and 21.1 m 3 s −1, 37.2 m 3 s −1 and 0.369 in the latter. A third example is used to show that if the annual maximum flows in a sample can all be reasonably described by a single GEV distribution function, then using PWMs of various
102
Q.J. Wang/Journal of Hydrology 194 (1997) 95–106
Fig. 2. Fitting of the GEV distribution to annual maximum flows at Wickliffe station of the Hopkins river, Victoria, Australia.
orders will produce only minimum change in the resulting distribution curve. The data are from the Savage’s Crossing station, with a catchment area of 10 334 km 2, of the Brisbane river, Queensland, Australia. A total of 60 annual maximum flows during the period 1910– 1969 are available. The data are plotted against the Gumbel reduced variable in Fig. 3. The fitted distributions using PWMs of r = 0, 1 and 2 and PWMs of r = 2, 3 and 4 are also shown in Fig. 3. The difference between the two curves is small, both describing the data reasonably well. The estimated parameter values for y, a and k are respectively 527 m 3 s −1, 579 m 3 s −1 and −0.305 in the former case, and 483 m 3 s −1, 689 m 3 s −1 and −0.224 in the latter.
6. Monte Carlo simulations Monte Carlo simulations have been conducted to study how quantile estimation is affected by using higher PWMs. In practice, the true underlying distribution is never known and often differs from the distribution function actually used. For this reason, the study is carried out by fitting the GEV distribution function to generated Wakeby samples. Houghton (1978) introduced the Wakeby distribution function for flood frequency analysis. It contains five parameters and is thus highly flexible. By varying the parameter values, Landwehr et al. (1979b,c, 1980) constructed six Wakeby distributions, identified as WA-1, WA-2, WA-3, WA-4, WA-5 and WA-6, to represent a wide range of skewness and kurtosis. They then used the six Wakeby distributions as parent distributions in Monte Carlo simulations to assess the performance of various estimation methods and various assumed distribution function forms. The same six Wakeby distributions are used in this study as parent distributions to
Q.J. Wang/Journal of Hydrology 194 (1997) 95–106
103
Fig. 3. Fitting of the GEV distribution to annual maximum flows at Savage’s Crossing station of the Brisbane river, Queensland, Australia.
generate random samples. A total of 100 000 samples of size 30 are used in each case. The samples are fitted to the GEV distribution function using first PWMs of r = 0, 1 and 2 and then PWMs of r = 2, 3 and 4. Quantile estimates are obtained to calculate their bias and mean square error. Fig. 4 and Fig. 5 show respectively bias in x(F = 0.98) and x(F = 0.99) estimates for the first five Wakeby distributions. Using PWMs of r = 2, 3 and 4 almost always leads to much smaller bias than using PWMs of r = 0, 1 and 2 for all five Wakeby distributions. Mean square error has been obtained also for quantiles x(F = 0.98) and x(F0.99). The results for the first five Wakeby distributions are presented in Fig. 6 in terms of relative
Fig. 4. Bias in x(F = 0.98) estimates from fitting the GEV distribution to the Wakeby samples.
104
Q.J. Wang/Journal of Hydrology 194 (1997) 95–106
Fig. 5. Bias in x(F = 0.99) estimates from fitting the GEV distribution to the Wakeby samples.
estimation efficiency, defined as f=
mean square error in estimate using r = 0, 1 and 2 mean square error in estimate using r = 2, 3 and 4
(21)
Fig. 6 shows that the relative estimation efficiency is generally greater than one, implying that using PWMs of r = 2, 3 and 4 almost always leads to much smaller variance than using PWMs of r = 0, 1 and 2 for all five Wakeby distributions. For the sixth Wakeby distribution, WA-6, of Landwehr et al. (1979b,c, 1980), using PWMs of r = 2, 3 and 4 rather than PWMs of r = 0, 1 and 2 results in a significant decrease in efficiency for quantile estimation. It should, however, be pointed out that the WA-6 distribution has a zero product–moment skewness and is not a typical distribution for describing annual maximum streamflows.
Fig. 6. Relative efficiency of x(F = 0.98) and x(F = 0.99) estimates from fitting the GEV distribution to the Wakeby samples.
Q.J. Wang/Journal of Hydrology 194 (1997) 95–106
105
7. Conclusions The use of higher PWMs for estimating flood frequency distributions has been proposed and examined. Higher PWMs give more weight to the larger flows. A distribution estimated by matching higher PWMs to their sample estimates is more closely fitted to the larger flows. Using a distribution so fitted, floods of large return periods are estimated from the trend shown by the larger flows in the recorded sample. This overcomes the problem that the smaller flows may exhibit different trend from the larger flows in a sample, and when a single smooth distribution function, such as the GEV distribution, is forced to fit all the data in such a sample, the fitted curve may have to take up a compromised position and shape, often resulting in poor fit to the larger flows in the sample and serious error in the estimation of floods of large return periods. The examples used in the paper demonstrate the merit of using higher PWMs to fit distributions for estimating floods of large return periods. The orders of the PWMs to be used do not have to be very high. It appears that orders r = 2 and above are sufficient. However, more study is needed to determine the most appropriate orders to be used. As a principle, there is a need to strike a balance between extracting from samples important information on the distribution tail and avoiding the estimation becoming too sensitive to sample outliers. The third example shows that if all the annual maximum flows in a sample can be reasonably described by a single GEV distribution function, then using PWMs of various orders produces only minimum change in the resulting distribution curve. Monte Carlo simulation results from fitting the GEV distribution to samples generated from the first five of the six Wakeby distributions of Landwehr et al. (1979b,c, 1980) show that using PWMs of order r = 2, 3 and 4 instead of PWMs of order r = 0, 1 and 2 almost consistently leads to improved quantile estimation in terms of both bias and mean square error. The method of using higher PWMs involves no more complication than the commonly accepted method of using the lowest PWMs, and the method is readily applicable to all distributions for which expressions of PWMs have been derived. It is suggested that higher PWMs be used as a matter of routine for estimating floods of large return periods.
Acknowledgements The author would like to express his thanks to Professor Tom McMahon and Murray Peel, CRC for Catchment Hydrology, University of Melbourne, for their discussions and comments on the paper.
References Davis, P.J., 1965. Gamma functions and related functions. In: M. Abramowitz and I.A. Stegun (Editors), Handbook of Mathematical Functions. Dover, New York, pp. 253–293. Greenwood, J.A., Landwehr, J.M., Matalas, N.C. and Wallis, J.R., 1979. Probability weighted moments: definition and relation to parameters of distribution expressible in inverse form. Water Resour. Res., 15(5): 1049– 1054.
106
Q.J. Wang/Journal of Hydrology 194 (1997) 95–106
Gringorten, I.I., 1963. A plotting rule for extreme probability paper. J. Geophys. Res., 68(3): 813–814. Hosking, J.R.M., Wallis, J.R. and Wood, E.F., 1985. Estimation of the generalized extreme value distribution by the method of probability weighted moments. Technometrics, 27(3): 251–261. Houghton, J.C., 1978. Birth of a parent: the Wakeby distribution for modeling flood flows. Water Resour. Res., 14(6): 1105–1110. Landwehr, J.M., Matalas, N.C. and Wallis, J.R., 1979a. Probability weighted moments compared with some traditional techniques in estimating Gumbel parameters and quantiles. Water Resour. Res., 15(5): 1055– 1064. Landwehr, J.M., Matalas, N.C. and Wallis, J.R., 1979b. Estimation of parameters and quantiles of Wakeby distributions. 1. Known lower bounds. Water Resour. Res., 15(6): 1361–1372. Landwehr, J.M., Matalas, N.C. and Wallis, J.R., 1979c. Estimation of parameters and quantiles of Wakeby distributions. 2. Unknown lower bounds. Water Resour. Res., 15(6): 1373–1379. Landwehr, J.M., Matalas, N.C. and Wallis, J.R., 1980. Quantile estimation with more or less floodlike distributions. Water Resour. Res., 16(3): 547–555. Leese, M.N., 1973. Use of censored data in the estimation of Gumbel distribution parameters for annual maximum flood series. Water Resour. Res., 9(6): 1534–1542. Phien, H.N. and Fang, T.S.E., 1989. Maximum likelihood estimation of the parameters and quantiles of the generalized extreme-value distribution from censored samples. J. Hydrol., 105: 139–155. Prescott, P. and Walden, A.T., 1983. Maximum likelihood estimation of the parameters of the three-parameter generalized extreme-value distribution from censored samples. J. Stat. Comput. Simul., 16: 241–250. Press, W.H., Teukolsky, S.A., Wetterling, W.T. and Flannery, B.P., 1992. Numerical Recipes in FORTRAN: the Art of Scientific Computing, 2nd edn. Cambridge University Press, New York, 963 pp. Wang, Q.J., 1990a. Estimation of the GEV distribution from censored samples by method of partial probability weighted moments. J. Hydrol., 120: 103–114. Wang, Q.J., 1990b. Unbiased estimation of probability weighted moments and partial probability weighted moments from systematic and historical flood information and their application to estimating the GEV distribution. J. Hydrol., 120: 115–124. Wang, Q.J., 1996. Using partial probability weighted moments to fit the extreme value distributions to censored samples. Water Resour. Res., 32(6): 1767–1771.