Stochastic Environmental Research and Risk Assessment 15 (2001) 284±309 Ó Springer-Verlag 2001

Estimation of confidence intervals of quantiles for the Weibull distribution J.-H. Heo, J. D. Salas, K.-D. Kim

284 Abstract. Estimation of con®dence limits and intervals for the two- and threeparameter Weibull distributions are presented based on the methods of moment (MOM), probability weighted moments (PWM), and maximum likelihood (ML). The asymptotic variances of the MOM, PWM, and ML quantile estimators are derived as a function of the sample size, return period, and parameters. Such variances can be used for estimating the con®dence limits and con®dence intervals of the population quantiles. Except for the two-parameter Weibull model, the formulas obtained do not have simple forms but can be evaluated numerically. Simulation experiments were performed to verify the applicability of the derived con®dence intervals of quantiles. The results show that overall, the ML method for estimating the con®dence limits performs better than the other two methods in terms of bias and mean square error. This is specially so for c 0:5 even for small sample sizes (e.g. N 10). However, the drawback of the ML method for determining the con®dence limits is that it requires that the shape parameter be bigger than 2. The Weibull model based on the MOM, ML, and PWM estimation methods was applied to ®t the distribution of annual 7-day low ¯ows and 6-h maximum annual rainfall data. The results showed that the differences in the estimated quantiles based on the three methods are not large, generally are less than 10%. However, the differences between the con®dence limits and con®dence intervals obtained by the three estimation methods may be more signi®cant. For instance, for the 7-day low ¯ows the ratio between the estimated con®dence interval to the estimated quantile based on ML is about 17% for T 2 while it is about 30% for estimation based on MOM and PWM methods. J.-H. Heo Department of Civil Engineering, Yonsei University Seoul, Korea J. D. Salas Department of Civil Engineering, Colorado State University Fort Collins, CO 80523, USA K.-D. Kim Korea Infrastructure Safety and Technology Cooperation Korea The research leading to this paper has been sponsored by the US National Science Foundation Grant CMS-9625685 on ``Uncertainty and Risk Analysis under Extreme Hydrologic Events'' and Internal Research Fund of Yonsei University. Acknowledgment is due to two anonymous reviewers who provided important suggestions that improved the paper.

In addition, the analysis of the rainfall data using the three-parameter Weibull showed that while ML parameters can be estimated, the corresponding con®dence limits and intervals could not be found because the shape parameter was smaller than 2.

1 Introduction The Weibull distribution has been quite popular in engineering since it was introduced several decades ago (Weibull, 1939; 1951). In the 1960s and 1970s, many papers were published in the literature concerning parameter estimation of the Weibull distribution based on complete and censored samples. Cohen (1965) described the maximum likelihood estimation method and derived the variance± covariance matrix of the parameters for a two-parameter Weibull distribution based on complete and censored samples. Likewise, Harter and Moore (1965) developed an iterative procedure to ®nd maximum likelihood estimates for a threeparameter Weibull distribution and illustrated numerical examples for the one-, two- and three-parameter Weibull model. Later, Harter and Moore (1967) derived the elements of the information matrix of the maximum likelihood estimates for a three-parameter Weibull distribution. Also, Haan and Beer (1967) developed iterative procedures based on the method of false position and secant method to ®nd the shape parameter for a three-parameter Weibull distribution. In addition, Thoman and Bain (1969) derived the con®dence limits of the parameters for a twoparameter Weibull distribution and Lemon (1975) developed the maximum likelihood estimators for the three-parameter Weibull model based on censored samples and gave the elements of the asymptotic variance±covariance matrix of the parameters based on the method of maximum likelihood. Greenwood et al. (1979) introduced the method of probability weighted moments for the Weibull model. In addition, the Weibull model has been widely used in hydrology for ®tting the frequency distribution of ¯ood and low ¯ow events. For example, Matalas (1963) found the Weibull distribution as an appropriate model for low ¯ow analysis and Kite (1988) showed several estimation techniques for this model. Rao (1981) compared the two- and three-parameter Weibull models for ®tting the frequency distribution of hydrologic data and provided some guidelines for selecting between the two- and the three-parameter Weibull models. Boes et al. (1989) applied a two-parameter Weibull model for regional ¯ood quantile estimation based on the index ¯ood assumption and compared the methods of moments (MOM), probability weighted moments (PWM), and maximum likelihood (ML) by simulation experiments. Also, Heo et al. (1990) derived the asymptotic variances of the quantiles obtained from regional analysis for these three estimation techniques for a two-parameter Weibull model. They also applied the regional Weibull model to ®t the frequency distribution of annual ¯ood data. Likewise, Dodson (1994) gave the con®dence intervals for the shape, scale, location parameters based on the method of maximum likelihood and linear estimation. Other available estimation techniques include those based on the distribution of log X (Johnson and Kotz, 1970), on order statistics (Johnson and Kotz, 1970), and the modi®ed moment estimation (Cohen et al., 1984). The Weibull distribution approximates the normal distribution as the shape parameter is about 3.6 in which case the skewness becomes zero (Johnson and Kotz, 1970). Also, it becomes the exponential and the Rayleigh distributions when the shape parameter is equal to one and two, respectively (Dodson, 1994). Despite its ¯exibility and popularity, not many papers have dealt with estimation of the

285

286

con®dence intervals of quantiles. Most of the results available are limited to the method of moments for the two- and three-parameter Weibull model (see for instance, Kite, 1988). This paper discusses parameter estimation procedures and derives and compares con®dence intervals on population quantiles for the three-parameter Weibull model based on the MOM, PWM, and ML. Simulation experiments are performed to investigate the applicability of derived con®dence intervals for ®nite samples. Examples based on annual 7-day low ¯ows of the Parana River, Argentina and 6-h maximum annual rainfall of Seoul, Korea are included to illustrate and compare the proposed estimation methods.

2 Model description The cumulative distribution function (cdf) and the probability density function (pdf) of the three-parameter Weibull model are de®ned respectively as (Johnson and Kotz, 1970) " # x f b ; xf 1 F x 1 exp a " # b x f b 1 x f b ; exp f x a a a

xf

2

in which a > 0 is the scale parameter, b > 0 is the shape parameter, and f is the location parameter. The three-parameter Weibull distribution is related to the GEV-3 distribution (NERC, 1975). If X is GEV-3 distributed with location parameter f0 , scale parameter a0 , and shape parameter b0 (b0 > 0 for the GEV-3 distribution), then X is Weibull distributed with parameters b 1=b0 ; a a0=b0 and f f0 a0=b0 . Also f 0 for the two-parameter Weibull distribution. The mean and variance of a three-parameter Weibull distribution are given by

l f aC 1 1=b

3

and

r2 a2 C 1 2=b

C2 1 1=b

4

respectively. The skewness coef®cient is given by

c

C 1 3=b

3C 1 2=bC 1 1=b 2C3 1 2=b C 1 2=b

C2 1 1=b3=2

5

where C is the gamma function. The skewness coef®cient c has a lower limit equal to 1:1396 (Gumbel, 1958). Note that the mean, variance, and skewness coef®cient exist for b > 0. For b > 1, the mode can be obtained from Eq. (2). It gives

b mod X f a

1 1=b b

6

287

Fig. 1. Typical pdfs of the Weibull distribution as a function of b for f 0 and a 1

The mode is at zero for 0 < b 1. In addition, the median of the Weibull distribution is given by (Johnson and Kotz, 1970)

med X f alog 21=b

7

Figure 1 shows some examples of the pdf of the Weibull distribution as a function of b for a 1:0 and f 0. Note that the three-parameter Weibull distribution is the exponential distribution if b 1.

3 Estimation of quantiles The quantile estimator X^T of a Weibull distribution can be obtained from Eq. (1) by replacing F x by 1 1=T and solving for x. It gives ^ X^T f^ a^ ln 1=T1=b

8

^ a^ and b^ are the parameter estimators, and T return period which is where f, de®ned here as 1=p (p = exceedance probability). Also, the estimator X^T may be written in terms of the sample mean l^, sample standard deviation r^, and the frequency factor K^T (Chow, 1951) as

X^T l^ K^T r^

9

in which K^T is obtained by combining Eqs. (3), (4), (8) and (9) as ^ ^ ln 1=T1=b C 1 1=b K^T ^ 1=2 ^ C 1 2=b C2 1 1=b

10

Note that the frequency factor K^T of Eq. (10) is a function of the shape parameter b which in turn is a function of the skewness coef®cient. Table 1 gives values of

Table 1. Frequency factor KT for the Weibull distribution Coef®cient Nonexceedance probability (q) of 0.01 0.2 0.5 0.8 skewness c Corresponding return period (T) 1.01 1.25 2 5

288

0.9

0.95

0.98

0.99

0.999

10

20

50

100

1000

1.3310 1.3410 1.3921 1.4492 1.5099 1.5731 1.6372 1.7002 1.7602 1.8152 1.8637 1.9048 1.9379 1.9630 1.9806 1.9913 1.9957 1.9948 1.9893 1.9799 1.9673 1.9521 1.9348 1.9159 1.8956 1.8743 1.8523 1.8298 1.8070 1.7840 1.7609 1.7378 1.7149 1.6921 1.6696 1.6473 1.6253 1.6037 1.5824 1.5615 1.5409 1.5207

1.5488 1.5627 1.6352 1.7181 1.8093 1.9078 2.0124 2.1210 2.2310 2.3399 2.4449 2.5439 2.6354 2.7182 2.7919 2.8564 2.9120 2.9593 2.9988 3.0312 3.0573 3.0777 3.0931 3.1041 3.1111 3.1147 3.1153 3.1133 3.1091 3.1028 3.0949 3.0854 3.0747 3.0629 3.0502 3.0366 3.0223 3.0075 2.9921 2.9763 2.9602 2.9439

1.6825 1.6990 1.7857 1.8864 1.9986 2.1221 2.2559 2.3979 2.5456 2.6959 2.8457 2.9919 3.1321 3.2645 3.3879 3.5015 3.6052 3.6990 3.7834 3.8587 3.9257 3.9850 4.0371 4.0828 4.1225 4.1569 4.1865 4.2117 4.2330 4.2508 4.2653 4.2770 4.2860 4.2927 4.2974 4.3001 4.3011 4.3005 4.2985 4.2953 4.2909 4.2856

2.0164 2.0403 2.1669 2.3178 2.4911 2.6886 2.9110 3.1575 3.4263 3.7142 4.0171 4.3303 4.6494 4.9701 5.2886 5.6019 5.9078 6.2044 6.4906 6.7657 7.0292 7.2811 7.5215 7.7505 7.9686 8.1760 8.3733 8.5608 8.7392 8.9087 9.0700 9.2233 9.3692 9.5081 9.6402 9.7661 9.8860 10.0003 10.1092 10.2132 10.3123 10.4069

)1.14

Lower limit for coef®cient of skewness

)1.04 )1.00 )0.80 )0.60 )0.40 )0.20 0.00 0.20 0.40 0.60 0.80 1.00 1.20 1.40 1.60 1.80 2.00 2.20 2.40 2.60 2.80 3.00 3.20 3.40 3.60 3.80 4.00 4.20 4.40 4.60 4.80 5.00 5.20 5.40 5.60 5.80 6.00 6.20 6.40 6.60 6.80 7.00

)3.0752 )3.0499 )2.9189 )2.7666 )2.5997 )2.4221 )2.2394 )2.0579 )1.8834 )1.7207 )1.5726 )1.4405 )1.3241 )1.2225 )1.1340 )1.0570 )0.9899 )0.9313 )0.8798 )0.8343 )0.7939 )0.7578 )0.7255 )0.6963 )0.6699 )0.6459 )0.6239 )0.6038 )0.5852 )0.5681 )0.5522 )0.5374 )0.5236 )0.5107 )0.4987 )0.4873 )0.4766 )0.4666 )0.4570 )0.4480 )0.4395 )0.4313

)0.7344 )0.7401 )0.7683 )0.7973 )0.8249 )0.8496 )0.8699 )0.8845 )0.8925 )0.8938 )0.8885 )0.8777 )0.8623 )0.8435 )0.8224 )0.7999 )0.7769 )0.7537 )0.7309 )0.7088 )0.6875 )0.6671 )0.6476 )0.6291 )0.6116 )0.5951 )0.5794 )0.5646 )0.5506 )0.5373 )0.5248 )0.5129 )0.5016 )0.4909 )0.4808 )0.4711 )0.4619 )0.4532 )0.4448 )0.4368 )0.4292 )0.4219

0.1550 0.1511 0.1306 0.1056 0.0766 0.0438 0.0076 )0.0308 )0.0705 )0.1098 )0.1477 )0.1831 )0.2151 )0.2436 )0.2682 )0.2892 )0.3069 )0.3214 )0.3332 )0.3426 )0.3500 )0.3557 )0.3598 )0.3627 )0.3646 )0.3656 )0.3659 )0.3656 )0.3648 )0.3635 )0.3619 )0.3601 )0.3580 )0.3557 )0.3533 )0.3508 )0.3482 )0.3455 )0.3428 )0.3400 )0.3373 )0.3345

0.8280 0.8305 0.8422 0.8527 0.8605 0.8646 0.8640 0.8581 0.8463 0.8287 0.8059 0.7785 0.7477 0.7145 0.6799 0.6446 0.6094 0.5749 0.5413 0.5090 0.4782 0.4488 0.4210 0.3947 0.3699 0.3466 0.3246 0.3039 0.2845 0.2662 0.2490 0.2329 0.2177 0.2033 0.1898 0.1771 0.1651 0.1538 0.1430 0.1329 0.1233 0.1142

1.1172 1.1237 1.1566 1.1916 1.2270 1.2614 1.2934 1.3214 1.3442 1.3606 1.3701 1.3727 1.3688 1.3589 1.3440 1.3249 1.3026 1.2778 1.2512 1.2234 1.1948 1.1660 1.1371 1.1085 1.0802 1.0525 1.0253 0.9989 0.9732 0.9483 0.9241 0.9007 0.8781 0.8563 0.8351 0.8148 0.7951 0.7760 0.7577 0.7400 0.7228 0.7063

the frequency factor for given values of the skewness coef®cient and return period. Such a frequency factor can also be used to estimate quantiles. Once the parameter estimates are obtained the mean, standard deviation, and skewness coef®cient of the Weibull model are obtained from Eqs. (3)±(5). Then the quantile for a given return period can be computed by using Eq. (9) and the frequency factor obtained from Table 1. For example, for a given sample data, the mean,

standard deviation, and skewness coef®cient based on the method of probability weighted moments are given by 9133.47, 2698.71 and 0.522, respectively. Then, the 100-year quantile is x^100 9133:47 2698:71 2:637 16250, in which K^100 2:637 is taken from Table 1 for given T 100 and c 0:522. The estimation of quantiles requires the estimation of parameters. Three methods of parameter estimation for the Weibull distribution are considered here. They are: the MOM, method of ML, and method of PWM. While these estimation methods have been widely described in literature for many models, several of the equations presented here are new.

3.1 Method of moments (MOM) ^ a^ and b^ can be obtained by substituting l, r and c in The moment estimators f, Eqs. (3)±(5) by their corresponding sample estimates l^, r^, and c^. Since the skewness coef®cient in Eq. (5) is only a function of the shape parameter b, one can estimate b by solving Eq. (5). A moment estimator of the shape parameter, b^ can be obtained from the approximate equation b^

0:729268

0:338679^ c 4:96077 ^ c 1:14

1:0422

11

0:683609ln ^ c 1:142

which is valid for 1:08 c^ 6:0 0:52 b^ 100 with a determination coef®cient R2 0:9999999 and standard error = 0.008. A better solution of b^ can be obtained by solving Eq. (5) numerically. For instance, to apply the Newton±Raphson method, Eq. (5) is rewritten as

^ G b

^ C 1 3=b

^ ^ 2C3 1 3=b ^ 3C 1 2=bC 1 1=b h i3=2 ^ ^ C 1 2=b C2 1 1=b

c^

12

The ®rst derivative of Eq. (12) with respect to b^ is

^ G0 b

1 h i5=2 ^2 C 1 2=b ^ ^ b C2 1 1=b h ^ 6C0 1 2=bC 1 ^ ^ 3C0 1 3=b 1=b i ^ ^ ^ 2 1 1=b ^ 2=b 6C0 1 1=bC 3C0 1 1=bC 1 h i ^ ^ C 1 2=b C2 1 1=b h i ^ ^ 3C 1 2=bC 1 ^ ^ 2C3 1 1=b C 1 3=b 1=b h i 0 0 ^ ^ ^ 13 3C 1 2=b 3C 1 1=bC 1 1=b

where C0 is the ®rst derivative of the gamma function. Therefore, starting with an initial estimate obtained from Eq. (11) the recursive equation to estimate b^ in the iteration i 1 is

289

b^i1 b^i

G b^i =G0 b^i

The calculations are repeated until satisfying the error criterion

b^ i1 b^i 0

17b

where xj is the order statistic such that x1 x2 xN . From Eq. (16) the ®rst three population PWMs are

A0 f aC 1 1=b h i A1 f a2 1=b C 1 1=b 2 A2 f a3

1=b

C 1 1=b=3

18 19 20

By substituting these three population PWMs by the corresponding sample PWMs, A^0 , A^1 , and A^2 , the PWM estimator of the shape parameter b is a solution of

1 1

3

1=b^

2

1=b^

3A^2 2A^1

A^0 A^0

21

^ By using the Newton±Raphson Equation (21) can be solved numerically for b. method the following equations are de®ned

^ G b

1 1

3

1=b^

2

1=b^

3A^2 2A^1

and

^ G0 b

^2

b 1

1 2

A^0 A^0

h 2 3

22

1=b^

1=b^

1

2

1=b^

ln 3 2

1=b^

1

3

1=b^

291

i ln 2 23

Then, b^ can be determined following a similar procedure as that described in Sect. 3.1 for MOM estimation. Also, in this case, the value of b^ from Eq. (11) can ^ Then, the PWM estimators of the be used as an initial value to solve Eq. (21) for b. parameters a and f may be obtained from Eqs. (18) and (19) as

a^ A^0

2A^1

.h 1

i ^ C 1 1=b

1=b^

2

24

and

f^ A^0

^ a^C 1 1=b

25

In addition, note that using the relationship between the Weibull and the GEV models as indicated before in Sect. 2, i.e. b 1=b0 , the shape parameter can be obtained approximately from the formula suggested by Hosking et al. (1985). Furthermore, the PWM estimators for a two-parameter Weibull distribution (f 0) may be obtained from Eqs. (18) and (19) as

b^ ln 2 ln A^0 =A^1

ln 2

26

and

^ a^ A^0 C 1 1=b

27

3.3 Method of maximum likelihood (ML) The log-likelihood function of a three-parameter Weibull distribution is given by LL x; f; a; b N ln b

Nb ln a b

1

N X i1

ln xi

f

N X xi i1

f

b

a 28

The partial derivatives of LL with respect to f; a; and b are equated to zero yielding

292

N X

N bX xi f b 1 0 a i1 a

oLL of

b

oLL oa

N Nb b X xi f b 0 a a i1 a

oLL N ob b

1

xi

f

1

i1

N ln a

N X

ln xi

29a

29b

f

i1

N X xi i1

a

f b xi ln

f a

0

29c

respectively. These three equations must be solved simultaneously to ®nd the estimators of the parameters f, a, and b. Based on the Newton±Raphson method, one can write

2

3 2 Df o2 LL=of2 4 Da 5 4 o2 LL=oaof Db o2 LL=obof

3 3 12 oLL=of o2 LL=ofob o2 LL=oaob 5 4 oLL=oa 5 oLL=ob o2 LL=ob2

o2 LL=ofoa o2 LL=oa2 o2 LL=oboa

30

where 1 represents the inverse of the matrix and the second partial derivatives of the log-likelihood function of the Weibull distribution are given in Appendix A. One can evaluate the matrices of the right hand side of Eq. (30) for given values of the sample x1 ; . . . ; xN and parameters fi ; ai ; and bi and then obtain the increments Dfi ; Dai ; and Dbi in which the subscript i de®nes the iteration. Thus, the new parameter estimates for iteration i 1 are then determined by

ki1 ki Dki until satisfying the error criterion jDki =ki1 j < e in which k represents any of the parameters f; a; b, and e is a speci®ed relative error. The expected values of the second derivatives of the inverse matrix in Eq. (30) are given in Appendix B. These expected values are also used to obtain the asymptotic variance of the quantile estimator as shown in the next section. The inverse of the square matrix II 1 in Eq. (30) is called the inverse of the information matrix, which may be expressed as

2

II

1

6 6 1 6 6 6 ND 6 6 4

a2 b b 12 a2 h b b 1 abf b 1

a2 h b b 1 a2 a b2 ag

3 abf b 1 7 7 7 7 ag 7 7 7 5 2 bc

31

where

a C 1

2=b1 C00 2

b 1 C00 2

C2 1 2 C0 2 p2 6

c C 1

2=b

f C 1

1=b1

C2 1

1=b1 w 1

32a 32b

1=b

C0 2 w 1

1=b2

32c 1=b

32d

g C 1

2=bC0 2

C2 1

h C 1

1=bf1 C00 2

1=b1 w 1 C0 21 w 1

D bc f 2

1=b

32e

1=bg

32f 32g

00

in which C 2 is the second derivative of the gamma function with argument 2. Note that the coef®cients a, c, and g are de®ned only if b > 2. For a two-parameter Weibull distribution, f 0 in Eqs. (29b) and (29c). Combining these equations gives

N b

N X

Nb ln xi

i1

N P

xbi ln xi

i1 N P

i1

xbi

293

0

33

which is only a function of b. Solving Eq. (33) for b gives the maximum likelihood ^ and then a^ can be obtained from Eq. (29b). estimator b,

4 Confidence intervals on quantiles The 1 d con®dence interval X1 imated by (Kite, 1988) X^1

d

X^T u1

^

d=2 ST

d

on the population quantiles may be approx-

34

where u1 d=2 is the 1 d=2 quantile of the standard normal distribution, X^T is the quantile estimator corresponding to return period T, and S^T is the standard deviation or standard error of X^T . The quantile estimator can be written as

35 X^T g h^1 ; h^2 ; h^3 where h^i denotes in general estimators of either moments or parameters. The asymptotic variance of X^T may be expressed as oXT 2 oXT 2 oXT 2 2 ^ ^ ^ ST Var XT Var h1 Var h2 Var h^3 oh1 oh2 oh3 oXT oXT oXT oXT ^ ^ 2 Cov h1 ; h2 2 Cov h^2 ; h^3 oh1 oh2 oh2 oh3 oXT oXT Cov h^1 ; h^3 2 36 oh1 oh3 The terms in the right hand side of Eq. (36) may be determined depending on the estimation method. In this paper, such standard error can be estimated by the methods of moments, probability weighted moments, and maximum likelihood as described below.

4.1 Standard error by moments By using the ®rst three sample moments for ^ hi in Eq. (36), the variance of X^T can be written as (Kite, 1988)

S2T Var X^T l2 =N 1 KT c1 KT2 c2 oKT =oc1 2c2 2

oKT =oc1 c4

3c21

1=4

6 KT c3

3c1 c3

6c2

6c1 c2 =4 10c1 =4 2 9c1 c2 =4 35c21 =4 9 37

where the cj 's are the so called cumulants which are de®ned as 3=2

294

c1 l3 =l2

38a

c2 l4 =l22

38b

5=2

c3 l5 =l2

38c

c4 l6 =l32

38d

and c1 c is the skewness coef®cient. In turn, the rth central moments, lr are

l2 a2 D2

D21

39a

l3 a3 D3

3D2 D1 2D31

39b

l4 a4 D4

4D3 D1 6D2 D21

l5 a5 D5

5D4 D1 10D3 D21

10D2 D21 4D51

l6 a6 D6

6D5 D1 15D4 D21

20D3 D31 15D2 D41

3D41

39c 39d 5D61

39e

in which Dr C 1 r=b. In addition, the derivative of KT with respect to c1 can be written as

oKT =oc1 oKT =ob ob=oc1

40

The two partial derivatives in the right hand side of Eq. (40) can be determined as follows. From Eq. (10) the derivative of KT with respect to b is

D1 w1 oKT ob

ln BB1=b D2

D21

b2 D2

B1=b

D1

3=2

D2 w2 D21 w1

D21

where wr w 1 r=b C0 1 r=b=C 1 r=b and B from Eq. (5) the derivative of c1 with respect to b becomes

oc1 3 D2 ob D3

D21

2D31 w1 .h 2 D2 w2 D21 w1 b D2

D3 w3 2D2 D1 w2 D2 D1 w1

3D2 D1 2D31

41

ln 1=T. Likewise,

D21 5=2

i

42

4.2 Standard error by probability weighted moments The variance of the PWM quantile estimator X^T for the three-parameter Weibull distribution is obtained by replacing the h^i 's in Eq. (36) by the parameter estimators. It gives

S2T

oXT 2 oXT 2 oXT 2 ^ ^ Var f Var ^ a Var b of oa ob oXT oXT oXT oXT ^ ^ Cov f; a^ 2 Cov ^ a; b 2 of oa oa ob oXT oXT ^ b ^ 2 Cov f; of ob

43

From Eq. (8) the derivatives of XT with respect to the parameters f; a; and b are respectively

oXT 1 of

44a

oXT ln 1=T1=b oa

44b

oXT ob

44c

a ln ln 1=T ln 1=T1=b b2

^ Var(^ ^ Cov(f, ^ a^), Cov(^ ^ and Cov(f, ^ b) ^ in Eq. (43) The terms Var(f), a), Var(b), a, b), are given by (see Appendix C for details)

^ W 2 A00 2W0 W1 A01 W 2 A11 2W0 Wb C1 H Var f 0 1 2W1 Wb C2 H Wb2 CH 2

45a

Var ^ a T02 A00 2T0 T1 A01 T12 A11 2T0 Tb C1 H 2T1 Tb C2 H Tb2 CH 2 45b ^ CH 2 Var b

45c

Cov ^f; ^a W0 T0 A00 W0 T1 A01 W1 T0 A01 W1 T1 A11 W0 Tb C1 H W1 Tb C2 H Wb T0 C1 H Wb T1 C2 H Wb Tb CH 2

45d

^ b ^ W0 C1 H W1 C2 H Wb CH 2 Cov f;

45e

^ T0 C1 H T1 C2 H Tb CH 2 Cov ^ a; b

45f

and the terms in Eqs. (45a) through (45f) are also given in Appendix C. Finally, using B ln 1=T the asymptotic variance of X^T for the three-parameter Weibull distribution becomes

S2T

1 a2 ^ ^ B2=b Var ^ a 4 ln B2 B2=b Var b Var f N b h i 2a 1=b 1=b 1=b ^ ^ ^ ^ ln BB Cov f; b B Cov ^ a; b 2B Cov f; a^ b2

46

295

For a two-parameter Weibull distribution, the asymptotic distribution of the sample PWMs becomes bivariate normal by dropping A^2 and the third line and column of the square matrix in Eq. (C3). Following a similar procedure as described for the three-parameter Weibull model, the asymptotic variance has a simpler form as (Heo et al., 1990)

S2T 296

nh i XT2 1 2 21 1=b S 1 21 1=b 4H 1=2 S2 C 1 2=b NC 1 1=b h i o 1 2 211=b S 1 211=b S2 C2 1 1=b 2

47 where S fw 1 1=b

ln ln 1=Tg= ln 0:5.

4.3 Standard error by maximum likelihood The variance of the ML estimator of quantiles for the three-parameter Weibull distribution can be obtained from Eq. (43) where the parameter estimators are now ML estimators. The derivatives of XT with respect to the parameters f; a; and b are the same as Eqs. (44a) through (44c). On the other hand, the variance and covariance terms for ML estimators are the elements of the inverse of information matrix in Eq. (31) (Mood et al., 1974). Hence ^ Var f

a2 b b 12 ND

48a

Var ^ a

a2 a b2 ND

48b

^ Var b

b2 c ND

48c

^ a^ Cov f; ^ Cov ^ a; b

a2 h b b 1ND abf b 1ND

^ b ^ ag Cov f; ND

48d 48e 48f

where the various coef®cients such as a, b, etc. are given by Eq. (32). Thus, substituting Eqs. (44) and (48) into (43) yields

S2T

( a2 b ND b 12

2B1=b B1=b h f ln B 2 a b b 1 b

)

2

2g ln B ln B

49 in which B ln 1=T. Note that for the two-parameter Weibull distribution f 0 and Eq. (43) reduces to

S2T

oXT 2 oXT 2 oXT oXT ^ ^ Cov ^ a; b Var ^ a Var b 2 oa ob oa ob ( ) XT2 C0 2 ln ln 1 q2 1 p2 =6 Nb2

50

5 Simulation experiments Simulation experiments were performed to ®nd out the applicability of the derived con®dence intervals of quantiles based on the methods of moments, probability weighted moments, and maximum likelihood. For this purpose, the location and scale parameters were set to f 0 and a 1, respectively, and the shape parameter b varied as 7.49, 3.62, 2.23, 1.56, and 0.74 which correspond to skewness coef®cients equal to 0:5, 0.0, 0.5, 1.0, and 3.0, respectively. Using a known parameter set (f; a; b), 10,000 sets of data were generated for sample sizes N = 10, 25, 50, and 100. For each generated data set and estimation method, the parameters were estimated and then the con®dence limits X^1 d of Eq. (34) were determined by using the estimates X^T and S^2T obtained from the generated data. These estimates were obtained respectively from Eq. (8) and from either Eqs. (37), (46), or (49) depending on the estimation method. In addition, the theoretical con®dence limits X1 d were obtained also from Eq. (34) but XT and ST were determined based on the assumed (known) parameter set. Then the relative bias RBIASq= BIAS X^1 d =X1 d and the relative root mean square error RRMSE = MSE X^1 d =X1 d were determined. The computations were performed for both the lower and upper con®dence limits. Only the results for the upper con®dence limits for c 0:5, 0.5 and 1.0 are displayed in Figs. 2±5. For N = 10, the RBIAS of the upper con®dence limit generally increases as c increases as shown in Fig. 2. Also it generally increases as T increases except for

Fig. 2. RBIAS vs. return period for the upper con®dence limits for N 10

297

298

Fig. 3. RBIAS vs. sample size for the upper con®dence limits for T 100

Fig. 4. RRMSE vs. return period for the upper con®dence limits for N 10

ML and c 1. MOM gives the smallest RBIAS for skewness coef®cient c 0:5; 0:5 while ML gives the smallest values for c 1. However, RBIAS of PWM rapidly increases as c increases. Furthermore, for T 100 (q 0:99) RBIAS of the upper con®dence limit decreases as the sample size N increases as shown in Fig. 3. RBIAS are about the same for the three estimation methods for c 0:5, MOM gives the smallest RBIAS for c 0:5, while ML gives the smallest value for c 1. In addition, for the PWM method RBIAS becomes consistently

299

Fig. 5. RRMSE vs. sample size for the upper con®dence limits for T 100

larger as c increases. Therefore, the results show that in terms of bias MOM is preferable for negative c while ML is best for positive c. Relative root mean square error (RRMSE) for the upper con®dence limit is shown in Fig. 4 for N 10. RRMSE generally increases as c increases except for the case of MOM for c 0:5. Also RRMSE generally increases as T increases. ML always gives the smallest RRMSE for all ranges of c and PWM gives the largest RRMSE for positive c. For T 100, RRMSE decreases as N increases as shown in Fig. 5. Also RRMSE generally increases as the skewness c increases except for MOM for c 0:5 and N 10. For c 0:5 and N 10, Figs. 2 and 3 show that MOM gives the smallest bias. Yet Figs. 4 and 5 show a large value of RRMSE. This occurred because of the large sample variance of the upper con®dence limit obtained for that particular case. Results for the bias and mean square error for the lower con®dence limits are not shown but can be summarized as follows. For N 10 the bias generally increases as T increases. Also the bias increases as c increases. Generally the biases are small but they can be large for c 1 (e.g. RBIAS is about 25% for the MOM and PWM methods). For T 100 the biases decrease as N increases as expected. The PWM method gives the smallest bias for c 0:5 but as c increases the ML gives the smallest bias. Regarding the mean square error, for N 10 RRMSE generally increases as T increases; this is specially noticeable for c 1. For T 100, RRMSE increases as c increases and it decreases as N increases as expected. In all cases the ML gives the smallest RRMSE. Overall, the ML method for estimating the con®dence limits performs better than the other two methods in terms of bias and mean square error. This is specially so for c 0:5 even for small sample sizes (e.g. N 10).

6 Applications The annual 7-day low ¯ows (cms) of the Parana River at Corrientes, Argentina for the period (1904±1982) and the annual 6-h maximum rainfall data (mm) of Seoul,

Table 2. Parameter estimates of the Weibull distribution for (a) annual 7-day low ¯ows of the Parana River at Corrientes, Argentine, and (b) 6-h maximum annual rainfall data at Seoul, South Korea Methods

Parameters Low ¯ows

300

Rainfall data

f^

a^

b^

f^

MOM

3367.15

6509.15

2.292

33.25

74.70

1.711

ML

3545.67

6313.66

2.238

38.19

68.58 112.58

1.568 2.669

PWM

3562.17

6290.95

2.177

37.90

68.92

1.553

a^

b^

2-Parameter Weibull

South Korea for the period (1925±1987) are used here to illustrate and compare the estimation of the parameters, quantiles, and con®dence intervals assuming the Weibull distribution and the various estimation methods presented herein. The parameters and quantiles for various return periods were estimated by the methods of moments, probability weighted moments, and maximum likelihood as outlined in the previous sections. Table 2 shows the estimated parameters obtained for the referred data. Note that for the rainfall data, the parameters of the 2-parameter Weibull model based on the ML method are also included because in estimating the con®dence intervals using the ML method, the estimated shape parameter should be greater than 2 because the argument of C 1 2=b in Eq. (B1) of Appendix B should be greater than zero. Thus, while the parameter estimates for the 3-parameter Weibull model can be obtained by the ML method, the con®dence intervals cannot be determined because b^ 1:568 < 2. Thus, for the rainfall data, the 2-parameter Weibull model was used for obtaining con®dence intervals of quantiles based on the ML method. The empirical and ®tted cdfs based on the three-parameter Weibull for the two sets of data analyzed are displayed in Figs. 6 and 7. The Weibull cdfs appear to ®t the empirical cdfs quite well. Furthermore, goodness of ®t tests such as the chisquare test and Kolmogorov±Smirnov test did not reject the Weibull model with 5% signi®cance level. The quantiles and 95% con®dence limits, corresponding to return periods T 1:01, 1.05, 1.10, 1.5, 2, 5, 10, 20, 50, 100, and 500 were determined. Table 3 shows the results for the annual 7-day low ¯ows. The estimated quantiles by the three methods are not the same but the differences among them are not too large. Generally they are less than 5% for T 2 (percentage relative to the smaller estimate when comparing between the estimates of any two methods). Somewhat larger differences are observed in the estimated con®dence limits. The maximum differences for both the upper and lower limits are less than 8% for T 500. More notable are the differences in the con®dence intervals as shown for instance in Fig. 8. The ratios of the con®dence interval to the estimated quantile (for a given T) obtained based on the ML method are the narrowest being about 17% for T 2, while the ratios obtained based on the MOM and PWM methods are wider, they are about 30% for T 2. For the rainfall data, the differences of quantiles based on the MOM and PWM methods (for the three-parameter Weibull) are less than 10%. Similar differences are obtained between the estimated con®dence limits by both estimation methods except for the lower con®dence limits for T < 1:11 year. Furthermore, the

301

Fig. 6. Empirical and ®tted cdfs of the 3-parameter Weibull distribution for the annual 7-day low ¯ows of the Parana River at Corrientes, Argentine

Fig. 7. Empirical and ®tted 3-parameter Weibull distribution cdfs for the annual 6-h maximum rainfall data of Seoul, South Korea

con®dence intervals for the three-parameter Weibull model based on the ML method could not be determined because the estimated shape parameter was less than 2. Therefore, for comparison the two-parameter Weibull model was applied based on the ML method (the estimated parameters are given in Table 2). Con®dence intervals of quantiles for the two- and three-parameter Weibull models are shown in Fig. 9. As expected, large differences are observed for the upper con®dence limits obtained for the two- and three-parameter Weibull models.

7 Summary and conclusions Estimation techniques for determining the con®dence intervals for the two- and three-parameter Weibull distributions are presented based on the MOM, PWM, and ML. All three estimation methods require an iterative or numerical solution to estimate the shape parameter. The asymptotic variances of the MOM, PWM, and ML quantile estimators for the Weibull distribution are derived as a function of the sample size, return period, and parameters. Such variances can be used for

Table 3. The estimates of quantiles and 95% con®dence intervals based on the methods of moments, probability weighted moments, and maximum likelihood (Newton±Raphson) for the annual 7-day low ¯ows of the Parana River at Corrientes, Argentine Method

Return period T

Lower con®dence limit

Quantile ^T X

Upper con®dence limit

MOM

1.01 1.05 1.10 1.50 2.00 5.00 10.00 20.00 50.00 100.00 500.00

3005.6 4340.9 5076.1 7113.4 8222.2 10595.7 11780.1 12667.3 13567.6 14117.0 15123.8

4241.6 5109.6 5700.9 7757.0 8914.4 11378.6 12733.8 13873.5 15171.1 16042.0 17813.0

5477.5 5878.2 6325.7 8400.6 9606.4 12161.4 13687.4 15079.8 16774.6 17966.9 20502.2

PWM

1.01 1.05 1.10 1.50 2.00 5.00 10.00 20.00 50.00 100.00 500.00

3358.3 4474.0 5117.2 7084.5 8199.8 10585.9 11799.4 12736.2 13717.4 14332.2 15492.1

4322.4 5133.2 5699.0 7717.7 8878.3 11390.3 12790.2 13975.9 15334.1 16250.1 18123.0

5286.6 5792.3 6280.7 8350.9 9556.8 12194.6 13780.9 15215.6 16950.8 18168.0 20753.9

ML

1.01 1.05 1.10 1.50 2.00 5.00 10.00 20.00 50.00 100.00 500.00

4056.0 4842.4 5382.3 7333.1 8457.8 10805.0 12033.1 13030.4 14132.7 14856.1 16293.8

4353.8 5182.9 5754.1 7763.4 8905.5 11355.5 12711.0 13854.8 15160.5 16038.9 17829.5

4651.6 5523.5 6125.9 8193.8 9353.2 11906.0 13389.0 14679.2 16188.4 17221.6 19365.2

302

estimating the con®dence limits and con®dence intervals of the population quantiles. Except for the two-parameter Weibull model, the formulas obtained do not have simple forms but can be evaluated numerically. Simulation experiments were performed to ®nd out the applicability of the derived con®dence intervals of quantiles based on the three estimation methods. The interest was to see the applicability of the estimations methods for various sample sizes, return periods, and skewness coef®cients. The results show that overall, the ML method for estimating the con®dence limits performs better than the other two methods in terms of bias and mean square error. This is specially so for c 0:5 even for small sample sizes (e.g. N 10). However, the drawback of the ML method for determining the con®dence limits is that it requires that the shape parameter be bigger than 2. Finally, the Weibull model was applied to ®t the distribution of annual 7-day low ¯ows and 6-h maximum annual rainfall data. The purpose was to illustrate

303

Fig. 8. Cumulative distribution of annual 7-day low ¯ows of the Parana River at Corrientes, Argentine and 95% con®dence intervals based on the Weibull-3 model and the MOM, PWM, and ML estimation methods

the applicability of the estimation methods to real data and to compare the values of estimated parameters, quantiles, and con®dence intervals based on the MOM, ML, and PWM methods. The results showed that the differences in the estimated quantiles based on the three methods are not large, generally are less than 10%. However, the differences between the con®dence limits and con®dence intervals obtained by the three estimation methods may be more signi®cant. For instance, for the 7-day low ¯ows the ratio between the estimated con®dence interval to the estimated quantile based on ML is about 17% for T 2 while it is about 30% for estimation based on MOM and PWM methods. The analysis of the rainfall data also illustrated the need of having available alternative estimation methods. For example, the analysis of the rainfall data using the three-parameter Weibull showed that while ML parameters can be estimated, the corresponding con®dence limits and intervals could not be found because the shape parameter was smaller than 2.

Appendix A: Second partial derivatives of the log-likelihood function for the Weibull distribution X b 2i hX 2 yi o2 LL=of2 b 1=a2 yi b o2 LL=oaof b=a2 o2 LL=ofob

X xi

X

yib

f

1

A1

1

A2 b=a

X

yib

1

ln yi

1=a

X

yib

1

A3

304

Fig. 9. Cumulative distribution of the annual 6-h maximum rainfall data of Seoul, Korea and 95% con®dence intervals based on the Weibull-2 and -3 models and the MOM, PWM, and ML estimation method

o2 LL=oa2

h b=a2 N

b 1

h o2 LL=oaob 1=a N o2 LL=ob2 N=b2 where yi xi

X

X

f=a and

ybi

X

b

ybi

X

i

i ybi ln yi

ybi ln yi 2

P

A4 A5 A6

represents summation from 1 to N.

Appendix B: Expected values of the second partial derivatives of the log-likelihood function for the Weibull distribution

o2 LL N b 12 E C 1 2=b a2 of2 2 o LL Nb2 E a2 oa2 2 o LL N 2 1 C00 2 E 2 ob b 2 o LL Nb2 C 2 1=b E a2 ofoa

B1 B2 B3 B4

o2 LL E ofob 2 o LL E oaob

N 1 a

1=bC 1

1=b1 w 1

1=b

N 0 C 2 a

B5 B6

where C and w are gamma and digamma functions and C0 2 and C00 2 are the ®rst and second partial derivatives of the gamma function with argument 2, respectively.

Appendix C: Derivation of the variance of the quantile estimator ST2 based on PWM method The quantile estimator of X^T can be written as a function of the ®rst three sample PWMs X^T f A^0 ; A^1 ; A^2

C1

where the sample PWM A^r is de®ned in Eq. (17). The variance S2T can be obtained by

S2T

oXT 2 oXT 2 oXT 2 Var A^0 Var A^1 Var A^2 oA0 oA1 oA2 oXT oXT oXT oXT ^ ^ 2 Cov A0 ; A1 2 Cov A^1 ; A^2 oA0 oA1 oA1 oA2 oXT oXT 2 Cov A^0 ; A^2 oA0 oA2

C2

However, oXT =oA1 in Eq. (C2) is not obtainable for some probability distributions such as the Weibull model, because X^T cannot be written explicitly as a function of the sample PWMs. In other words, the quantile estimator in Eq. (8) cannot be explicitly written as a function of the sample PWMs because the shape parameter estimator is implicitly expressed as a function of the sample PWMs as shown in Eq. (21). Therefore, in this case, the variance±covariance matrix of the sample PWMs should be obtained ®rstly. Then this matrix is transformed to the variance±covariance matrix of the parameter estimators by using a Jacobian transformation. Finally, the variance of the quantile estimator can be obtained by using Eq. (43). Details of the procedure follows. The asymptotic distribution of the sample PWMs can be written as (Chernoff et al., 1967; Hosking et al., 1985)

2

3 0 2 A0 A00 =N A^0 4 A^1 5 [email protected] A1 ; 4 A01 =N A2 A02 =N A^2

A01 =N A11 =N A12 =N

31 A02 =N A12 =N 5A A22 =N

C3

where reads ``is asymptotically distributed as'', TVN is an abbreviation for trivariate normal distribution, and Aij are given by (Heo et al., 1990)

A00 a2 C 1 2=b

C2 1 1=b

C4a

305

h A01 a2 =2 2

2=b

nh A02 a2 =2 3

C 1 2=b 1

2=b

2

2=b

H 1=2

C 1 2=b A11 a2 2

306

21

1=b

1=b

2

i C2 1 1=b

C4b

i

2 3

1=b

o C2 1 1=b

H 1=2C 1 2=b C2 1 1=b h A12 a2 =2 3 2=b H 1=3C 1 2=b 2:6 1=b

C4c

2=b

C4d 2

1=b

C2 1 1=b

i C4e

A22 a2 3

2=b

H 2=3C 1 2=b

C2 1 1=b

C4f

where H z is a hypergeometric function de®ned by

H z F 2=b; 1=b; 1

1=b; z

1 1=b X C n 2=b zn C 2=b n0 n 1=bn!

C5

For the three-parameter Weibull distribution, the asymptotic variance of the PWM estimator of quantile, X^T can be found by using successively the following transformations

A^0 A^1 A^2

A^0 ) A^1 A^2 R

A^0 ) A^1 A^2 b^

f^ ) a^ ) X^T b^

C6

where ) reads ``is transformed as''. As mentioned for Eq. (21), the shape parameter estimator is a function of the sample PWMs, but cannot be written explicitly, thus we need a transformation. In the ®rst transformation, let R 3A^2 A^0 = 2A^1 A^0 which is the right hand side of Eq. (21). Then b^ is ^ ^ given implicitly by 1 3 1=b = 1 2 1b R in the second transformation. Additional details of the transformations follows. (1) 1st Transformation. For the 1st transformation, the Jacobian matrix is given by

2 J1

oA^0 6 oA^^0 6 oA1 6 oA^0 6 oA^ 6 2 4 oA^0 oR oA^0

oA^0 oA^1 oA^1 oA^1 oA^2 oA^1 oR oA^1

oA^0 oA^2 oA^1 oA^2 oA^2 oA^2 oR oA^2

3

2 7 7 6 7 6 74 7 5

1 0 0

3A2 2A1 2A1 A0 2

0 1 0

2 3A2 A0 2A1 A0 2

0 0 1

3 2A1 A0 2A1 A0 2

3 7 7 5

C7

Then the variance±covariance matrix in the 1st transformation, i.e. of the 2nd column terms in (C6) can be obtained by R1 J1 R0 J10 where R0 is the variance± covariance matrix in Eq. (C3). Thus

2

A00 =N 6 A00 =N R1 6 4 A00 =N C1 =N

A01 =N A01 =N A01 =N C2 =N

3 C1 =N C2 =N 7 7 C3 =N 5 C=N

A02 =N A02 =N A02 =N C3 =N

C8

where the elements are de®ned in Eqs. (C4a) through (C4f) and

h C A00 3

1=b

2

1=b

2A01 3

1=b

1 3A02 2

i. 1 M

1=b

C9a h C1 A01 3 h C2 A02 3

1=b

2

1=b

1=b

1=b

2

2A11 3

1 3A12 2

1=b

2A12 3

1=b

1=b

1 3A22 2

i. 1 M

C9b

i. 1 M

1=b

C9c M aC 1 1=b 1

2

1=b

2

C9d

(2) 2nd Transformation. In this case, the Jacobian matrix is

2 oA^

0

J2

^ 6 ooAA^0 6 1 6 oA^0 6 oA^2 6 ^ 4 oA0 ob^ oA^0

oA^0 oA^1 oA^1 oA^1 oA^2 oA^1 ob^ oA^1

where R 1

ob^ oR ln 22

oA^0 oA^2 oA^1 oA^2 oA^2 oA^2 ob^ oA^2

3

1=b 1

1=b

oA^0 oR oA^1 oR oA^2 oR ob^ oR

= 1

3

2 1 7 7 6 7 60 74 0 7 5 0 2

1=b

3 0 0 0 0 7 7 1 0 5 ^ 0 ob=oR

0 1 0 0

C10

and

2 b2 1 2 1=b 3 1=b ln 33

1=b 1

2

C11

1=b

Thus, the variance±covariance matrix in the 2nd transformation is given by

2

A00 =N 6 A00 =N R2 J2 R1 J20 6 4 A00 =N C1 H=N

A01 =N A01 =N A01 =N C2 H=N

A02 =N A02 =N A02 =N C3 H=N

3 C1 H=N C2 H=N 7 7 C3 H=N 5 CH 2 =N

C12

(3) 3rd Transformation. The Jacobian matrix in the 3rd transformation is

2

^ A^0 of=o 4 J3 o^ a=oA^0 0

^ A^1 of=o o^ a=oA^0 0

^ A^2 of=o o^ a=oA^0 0

3 2 ^ b^ W0 of=o o^ a=ob^ 5 4 T0 0 1

W1 T1 0

3 0 Wb 0 Tb 5 0 1 C13

307

where the elements of above matrix are given by

W0 of=oA0

2

1=b

. 1

. W1 of=oA1 2 1 Wb of=ob 308

2

Tb oa=ob h 1 2

1=b

1=b

1=b

.h 2 1

1=b

2

2

C14b 1=b

.h b2 1

2

1=b

i

C14c

i C 1 1=b

1=b

w 1 1=b

C14a

aC 1 1=b ln 22

.h T0 oa=oA0 1 1 T1 oa=oA1

2

C14d

i C 1 1=b

ln 22

1=b

i b2 1

C14e

2

1=b

2

C 1 1=b C14f

Then the variance±covariance matrix can be obtained by

2

Z 1 4 ff 0 Zfa R3 J3 R2 J3 N Z fb

Zfa Zaa Zab

3 Zfb Zab 5 Zbb

C15

where the diagonal and off-diagonal terms are the variances and covariances of estimators, respectively, given by Eqs. (45a) through (45f). (4) 4th Transformation. The quantile estimator of the 3-parameter Weibull distribution for return period T is ^ X^T f^ a^ ln 1=T1=b

C16

and the Jacobian matrix is given by

"

oX^T oX^T oX^T J4 a ob^ of^ o^

#

where each element is given in Eqs. (44a) through (44c). Finally, the S2T of X^T can be obtained by

S2T J4 R3 J40

C17

and the result is given in Eq. (49).

References

Boes DC, Heo JH, Salas JD (1989) Regional ¯ood quantile estimation for a Weibull model. Water Resour. Res., 25: 979±990 Chernoff H, Gastwirth JL, Johns MV (1967) Asymptoti distribution of linear combinations of functions of order statistics with applications to estimation. Annals Math. Statist., 38: 52±72

Chow VT (1951) A general formula for hydrologic frequency analysis. Trans. Geophys. Union, 32: 231±237 Cohen AC (1965) Maximum likelihood estimation in the Weibull distribution based on complete and on censored samples. Techno., 7: 579±588 Cohen AC, Whitten BJ, Ding Y (1984) Modi®ed moment estimation for the three-parameter Weibull distribution. J. Qual. Tech., 16: 159±167 Dodson B (1994) Weibull Analysis. ASQC Press, Wisconsin Greenwood JA, Landwehr JM, Matalas NC, Wallis JR (1979) Probability weighted moments: de®nition and relation to parameters of several distributions expressible in inverse form. Water Resour. Res., 15: 1049±1054 Gumbel EJ (1958) Statistic of Extreme. Columbia University Press, NY Haan CT, Beer CE (1967) Determination of maximum likelihood estimators for the three parameter Weibull distribution. Iowa State J. Science, 42: 37±42. Iowa State University Ames, Iowa Harter HL, Moore AH (1965) Maximum likelihood estimation of the parameters of gamma and Weibull populations from complete and from censored samples. Techno., 7: 639±643 Harter HL, Moore AH (1967) Asymptotic variances and covariances of maximum likelihood estimators, from censored samples, of the parameters of Weibull and gamma populations. Annals Math. Statis., 38: 557±570 Heo JH, Boes DC, Salas JD (1990) Regional ¯ood frequency modeling and estimation. Water Resour., Paper No. 101, Colorado State University, Fort Collins, Colorado Hosking JRM, Wallis JR, Wood EF (1985) Estimation of the generalized extreme-value distribution by the method of probability weighted moments. Techno., 27: 251±261 Jenkinson AF (1969) Statistics of extremes. In: Estimation of Maximum Floods. Technical Report No. 98, WMO Johnson NL, Kotz S (1970) Continuous Univariate Distribution-1. Houghton Mif¯in Company, Boston Kite GW (1988) Frequency and Risk Analysis in Hydrology. Water Resources Publications, Littleton, Colorado Landwehr JM, Matalas NC (1979) Probability weighted moments compared with some traditional techniques in estimating Gumbel parameters and quantiles. Water Resour. Res., 15: 1055±1064 Lemon GH (1975) Maximum likelihood estimation for the three parameter Weibull distribution based on censored samples. Techno., 17: 247±254 Matalas NC (1963) Probability Distribution of Low Flows. USGS Prof. Paper. 434-A Mood AM, Graybill FA, Boes DC (1974) Introduction to the Theory of Statistics, (3rd Edn). McGraw-Hill, NY National Environmental Research Council (1975) Flood Studies Report. vol. 1, London, UK Prescott P, Walden AT (1983) Maximum likelihood estimation of the three-parameter generalized extreme-value distribution from censored samples. J. Statis. Comput. Simul., 16: 241±250 Rao DV (1981) Three-parameter probability distributions. J. Hydraul. Div. Am. Soc. Civ. Eng., 109: 339±358 Thoman DR, Bain LJ (1969) Inferences on the parameters of the Weibull distribution. Techno., 11: 445±460 Weibull W (1939) A statistical theory of the strength of materials. Ing. Vetenskaf Akad. Handl. 151 Weibull W (1951) A statistical distribution function of wide applicability. J. Appl. Mech., Am. Soc. Mech. Eng., 18: 293±297

309

Estimation of confidence intervals of quantiles for the Weibull distribution J.-H. Heo, J. D. Salas, K.-D. Kim

284 Abstract. Estimation of con®dence limits and intervals for the two- and threeparameter Weibull distributions are presented based on the methods of moment (MOM), probability weighted moments (PWM), and maximum likelihood (ML). The asymptotic variances of the MOM, PWM, and ML quantile estimators are derived as a function of the sample size, return period, and parameters. Such variances can be used for estimating the con®dence limits and con®dence intervals of the population quantiles. Except for the two-parameter Weibull model, the formulas obtained do not have simple forms but can be evaluated numerically. Simulation experiments were performed to verify the applicability of the derived con®dence intervals of quantiles. The results show that overall, the ML method for estimating the con®dence limits performs better than the other two methods in terms of bias and mean square error. This is specially so for c 0:5 even for small sample sizes (e.g. N 10). However, the drawback of the ML method for determining the con®dence limits is that it requires that the shape parameter be bigger than 2. The Weibull model based on the MOM, ML, and PWM estimation methods was applied to ®t the distribution of annual 7-day low ¯ows and 6-h maximum annual rainfall data. The results showed that the differences in the estimated quantiles based on the three methods are not large, generally are less than 10%. However, the differences between the con®dence limits and con®dence intervals obtained by the three estimation methods may be more signi®cant. For instance, for the 7-day low ¯ows the ratio between the estimated con®dence interval to the estimated quantile based on ML is about 17% for T 2 while it is about 30% for estimation based on MOM and PWM methods. J.-H. Heo Department of Civil Engineering, Yonsei University Seoul, Korea J. D. Salas Department of Civil Engineering, Colorado State University Fort Collins, CO 80523, USA K.-D. Kim Korea Infrastructure Safety and Technology Cooperation Korea The research leading to this paper has been sponsored by the US National Science Foundation Grant CMS-9625685 on ``Uncertainty and Risk Analysis under Extreme Hydrologic Events'' and Internal Research Fund of Yonsei University. Acknowledgment is due to two anonymous reviewers who provided important suggestions that improved the paper.

In addition, the analysis of the rainfall data using the three-parameter Weibull showed that while ML parameters can be estimated, the corresponding con®dence limits and intervals could not be found because the shape parameter was smaller than 2.

1 Introduction The Weibull distribution has been quite popular in engineering since it was introduced several decades ago (Weibull, 1939; 1951). In the 1960s and 1970s, many papers were published in the literature concerning parameter estimation of the Weibull distribution based on complete and censored samples. Cohen (1965) described the maximum likelihood estimation method and derived the variance± covariance matrix of the parameters for a two-parameter Weibull distribution based on complete and censored samples. Likewise, Harter and Moore (1965) developed an iterative procedure to ®nd maximum likelihood estimates for a threeparameter Weibull distribution and illustrated numerical examples for the one-, two- and three-parameter Weibull model. Later, Harter and Moore (1967) derived the elements of the information matrix of the maximum likelihood estimates for a three-parameter Weibull distribution. Also, Haan and Beer (1967) developed iterative procedures based on the method of false position and secant method to ®nd the shape parameter for a three-parameter Weibull distribution. In addition, Thoman and Bain (1969) derived the con®dence limits of the parameters for a twoparameter Weibull distribution and Lemon (1975) developed the maximum likelihood estimators for the three-parameter Weibull model based on censored samples and gave the elements of the asymptotic variance±covariance matrix of the parameters based on the method of maximum likelihood. Greenwood et al. (1979) introduced the method of probability weighted moments for the Weibull model. In addition, the Weibull model has been widely used in hydrology for ®tting the frequency distribution of ¯ood and low ¯ow events. For example, Matalas (1963) found the Weibull distribution as an appropriate model for low ¯ow analysis and Kite (1988) showed several estimation techniques for this model. Rao (1981) compared the two- and three-parameter Weibull models for ®tting the frequency distribution of hydrologic data and provided some guidelines for selecting between the two- and the three-parameter Weibull models. Boes et al. (1989) applied a two-parameter Weibull model for regional ¯ood quantile estimation based on the index ¯ood assumption and compared the methods of moments (MOM), probability weighted moments (PWM), and maximum likelihood (ML) by simulation experiments. Also, Heo et al. (1990) derived the asymptotic variances of the quantiles obtained from regional analysis for these three estimation techniques for a two-parameter Weibull model. They also applied the regional Weibull model to ®t the frequency distribution of annual ¯ood data. Likewise, Dodson (1994) gave the con®dence intervals for the shape, scale, location parameters based on the method of maximum likelihood and linear estimation. Other available estimation techniques include those based on the distribution of log X (Johnson and Kotz, 1970), on order statistics (Johnson and Kotz, 1970), and the modi®ed moment estimation (Cohen et al., 1984). The Weibull distribution approximates the normal distribution as the shape parameter is about 3.6 in which case the skewness becomes zero (Johnson and Kotz, 1970). Also, it becomes the exponential and the Rayleigh distributions when the shape parameter is equal to one and two, respectively (Dodson, 1994). Despite its ¯exibility and popularity, not many papers have dealt with estimation of the

285

286

con®dence intervals of quantiles. Most of the results available are limited to the method of moments for the two- and three-parameter Weibull model (see for instance, Kite, 1988). This paper discusses parameter estimation procedures and derives and compares con®dence intervals on population quantiles for the three-parameter Weibull model based on the MOM, PWM, and ML. Simulation experiments are performed to investigate the applicability of derived con®dence intervals for ®nite samples. Examples based on annual 7-day low ¯ows of the Parana River, Argentina and 6-h maximum annual rainfall of Seoul, Korea are included to illustrate and compare the proposed estimation methods.

2 Model description The cumulative distribution function (cdf) and the probability density function (pdf) of the three-parameter Weibull model are de®ned respectively as (Johnson and Kotz, 1970) " # x f b ; xf 1 F x 1 exp a " # b x f b 1 x f b ; exp f x a a a

xf

2

in which a > 0 is the scale parameter, b > 0 is the shape parameter, and f is the location parameter. The three-parameter Weibull distribution is related to the GEV-3 distribution (NERC, 1975). If X is GEV-3 distributed with location parameter f0 , scale parameter a0 , and shape parameter b0 (b0 > 0 for the GEV-3 distribution), then X is Weibull distributed with parameters b 1=b0 ; a a0=b0 and f f0 a0=b0 . Also f 0 for the two-parameter Weibull distribution. The mean and variance of a three-parameter Weibull distribution are given by

l f aC 1 1=b

3

and

r2 a2 C 1 2=b

C2 1 1=b

4

respectively. The skewness coef®cient is given by

c

C 1 3=b

3C 1 2=bC 1 1=b 2C3 1 2=b C 1 2=b

C2 1 1=b3=2

5

where C is the gamma function. The skewness coef®cient c has a lower limit equal to 1:1396 (Gumbel, 1958). Note that the mean, variance, and skewness coef®cient exist for b > 0. For b > 1, the mode can be obtained from Eq. (2). It gives

b mod X f a

1 1=b b

6

287

Fig. 1. Typical pdfs of the Weibull distribution as a function of b for f 0 and a 1

The mode is at zero for 0 < b 1. In addition, the median of the Weibull distribution is given by (Johnson and Kotz, 1970)

med X f alog 21=b

7

Figure 1 shows some examples of the pdf of the Weibull distribution as a function of b for a 1:0 and f 0. Note that the three-parameter Weibull distribution is the exponential distribution if b 1.

3 Estimation of quantiles The quantile estimator X^T of a Weibull distribution can be obtained from Eq. (1) by replacing F x by 1 1=T and solving for x. It gives ^ X^T f^ a^ ln 1=T1=b

8

^ a^ and b^ are the parameter estimators, and T return period which is where f, de®ned here as 1=p (p = exceedance probability). Also, the estimator X^T may be written in terms of the sample mean l^, sample standard deviation r^, and the frequency factor K^T (Chow, 1951) as

X^T l^ K^T r^

9

in which K^T is obtained by combining Eqs. (3), (4), (8) and (9) as ^ ^ ln 1=T1=b C 1 1=b K^T ^ 1=2 ^ C 1 2=b C2 1 1=b

10

Note that the frequency factor K^T of Eq. (10) is a function of the shape parameter b which in turn is a function of the skewness coef®cient. Table 1 gives values of

Table 1. Frequency factor KT for the Weibull distribution Coef®cient Nonexceedance probability (q) of 0.01 0.2 0.5 0.8 skewness c Corresponding return period (T) 1.01 1.25 2 5

288

0.9

0.95

0.98

0.99

0.999

10

20

50

100

1000

1.3310 1.3410 1.3921 1.4492 1.5099 1.5731 1.6372 1.7002 1.7602 1.8152 1.8637 1.9048 1.9379 1.9630 1.9806 1.9913 1.9957 1.9948 1.9893 1.9799 1.9673 1.9521 1.9348 1.9159 1.8956 1.8743 1.8523 1.8298 1.8070 1.7840 1.7609 1.7378 1.7149 1.6921 1.6696 1.6473 1.6253 1.6037 1.5824 1.5615 1.5409 1.5207

1.5488 1.5627 1.6352 1.7181 1.8093 1.9078 2.0124 2.1210 2.2310 2.3399 2.4449 2.5439 2.6354 2.7182 2.7919 2.8564 2.9120 2.9593 2.9988 3.0312 3.0573 3.0777 3.0931 3.1041 3.1111 3.1147 3.1153 3.1133 3.1091 3.1028 3.0949 3.0854 3.0747 3.0629 3.0502 3.0366 3.0223 3.0075 2.9921 2.9763 2.9602 2.9439

1.6825 1.6990 1.7857 1.8864 1.9986 2.1221 2.2559 2.3979 2.5456 2.6959 2.8457 2.9919 3.1321 3.2645 3.3879 3.5015 3.6052 3.6990 3.7834 3.8587 3.9257 3.9850 4.0371 4.0828 4.1225 4.1569 4.1865 4.2117 4.2330 4.2508 4.2653 4.2770 4.2860 4.2927 4.2974 4.3001 4.3011 4.3005 4.2985 4.2953 4.2909 4.2856

2.0164 2.0403 2.1669 2.3178 2.4911 2.6886 2.9110 3.1575 3.4263 3.7142 4.0171 4.3303 4.6494 4.9701 5.2886 5.6019 5.9078 6.2044 6.4906 6.7657 7.0292 7.2811 7.5215 7.7505 7.9686 8.1760 8.3733 8.5608 8.7392 8.9087 9.0700 9.2233 9.3692 9.5081 9.6402 9.7661 9.8860 10.0003 10.1092 10.2132 10.3123 10.4069

)1.14

Lower limit for coef®cient of skewness

)1.04 )1.00 )0.80 )0.60 )0.40 )0.20 0.00 0.20 0.40 0.60 0.80 1.00 1.20 1.40 1.60 1.80 2.00 2.20 2.40 2.60 2.80 3.00 3.20 3.40 3.60 3.80 4.00 4.20 4.40 4.60 4.80 5.00 5.20 5.40 5.60 5.80 6.00 6.20 6.40 6.60 6.80 7.00

)3.0752 )3.0499 )2.9189 )2.7666 )2.5997 )2.4221 )2.2394 )2.0579 )1.8834 )1.7207 )1.5726 )1.4405 )1.3241 )1.2225 )1.1340 )1.0570 )0.9899 )0.9313 )0.8798 )0.8343 )0.7939 )0.7578 )0.7255 )0.6963 )0.6699 )0.6459 )0.6239 )0.6038 )0.5852 )0.5681 )0.5522 )0.5374 )0.5236 )0.5107 )0.4987 )0.4873 )0.4766 )0.4666 )0.4570 )0.4480 )0.4395 )0.4313

)0.7344 )0.7401 )0.7683 )0.7973 )0.8249 )0.8496 )0.8699 )0.8845 )0.8925 )0.8938 )0.8885 )0.8777 )0.8623 )0.8435 )0.8224 )0.7999 )0.7769 )0.7537 )0.7309 )0.7088 )0.6875 )0.6671 )0.6476 )0.6291 )0.6116 )0.5951 )0.5794 )0.5646 )0.5506 )0.5373 )0.5248 )0.5129 )0.5016 )0.4909 )0.4808 )0.4711 )0.4619 )0.4532 )0.4448 )0.4368 )0.4292 )0.4219

0.1550 0.1511 0.1306 0.1056 0.0766 0.0438 0.0076 )0.0308 )0.0705 )0.1098 )0.1477 )0.1831 )0.2151 )0.2436 )0.2682 )0.2892 )0.3069 )0.3214 )0.3332 )0.3426 )0.3500 )0.3557 )0.3598 )0.3627 )0.3646 )0.3656 )0.3659 )0.3656 )0.3648 )0.3635 )0.3619 )0.3601 )0.3580 )0.3557 )0.3533 )0.3508 )0.3482 )0.3455 )0.3428 )0.3400 )0.3373 )0.3345

0.8280 0.8305 0.8422 0.8527 0.8605 0.8646 0.8640 0.8581 0.8463 0.8287 0.8059 0.7785 0.7477 0.7145 0.6799 0.6446 0.6094 0.5749 0.5413 0.5090 0.4782 0.4488 0.4210 0.3947 0.3699 0.3466 0.3246 0.3039 0.2845 0.2662 0.2490 0.2329 0.2177 0.2033 0.1898 0.1771 0.1651 0.1538 0.1430 0.1329 0.1233 0.1142

1.1172 1.1237 1.1566 1.1916 1.2270 1.2614 1.2934 1.3214 1.3442 1.3606 1.3701 1.3727 1.3688 1.3589 1.3440 1.3249 1.3026 1.2778 1.2512 1.2234 1.1948 1.1660 1.1371 1.1085 1.0802 1.0525 1.0253 0.9989 0.9732 0.9483 0.9241 0.9007 0.8781 0.8563 0.8351 0.8148 0.7951 0.7760 0.7577 0.7400 0.7228 0.7063

the frequency factor for given values of the skewness coef®cient and return period. Such a frequency factor can also be used to estimate quantiles. Once the parameter estimates are obtained the mean, standard deviation, and skewness coef®cient of the Weibull model are obtained from Eqs. (3)±(5). Then the quantile for a given return period can be computed by using Eq. (9) and the frequency factor obtained from Table 1. For example, for a given sample data, the mean,

standard deviation, and skewness coef®cient based on the method of probability weighted moments are given by 9133.47, 2698.71 and 0.522, respectively. Then, the 100-year quantile is x^100 9133:47 2698:71 2:637 16250, in which K^100 2:637 is taken from Table 1 for given T 100 and c 0:522. The estimation of quantiles requires the estimation of parameters. Three methods of parameter estimation for the Weibull distribution are considered here. They are: the MOM, method of ML, and method of PWM. While these estimation methods have been widely described in literature for many models, several of the equations presented here are new.

3.1 Method of moments (MOM) ^ a^ and b^ can be obtained by substituting l, r and c in The moment estimators f, Eqs. (3)±(5) by their corresponding sample estimates l^, r^, and c^. Since the skewness coef®cient in Eq. (5) is only a function of the shape parameter b, one can estimate b by solving Eq. (5). A moment estimator of the shape parameter, b^ can be obtained from the approximate equation b^

0:729268

0:338679^ c 4:96077 ^ c 1:14

1:0422

11

0:683609ln ^ c 1:142

which is valid for 1:08 c^ 6:0 0:52 b^ 100 with a determination coef®cient R2 0:9999999 and standard error = 0.008. A better solution of b^ can be obtained by solving Eq. (5) numerically. For instance, to apply the Newton±Raphson method, Eq. (5) is rewritten as

^ G b

^ C 1 3=b

^ ^ 2C3 1 3=b ^ 3C 1 2=bC 1 1=b h i3=2 ^ ^ C 1 2=b C2 1 1=b

c^

12

The ®rst derivative of Eq. (12) with respect to b^ is

^ G0 b

1 h i5=2 ^2 C 1 2=b ^ ^ b C2 1 1=b h ^ 6C0 1 2=bC 1 ^ ^ 3C0 1 3=b 1=b i ^ ^ ^ 2 1 1=b ^ 2=b 6C0 1 1=bC 3C0 1 1=bC 1 h i ^ ^ C 1 2=b C2 1 1=b h i ^ ^ 3C 1 2=bC 1 ^ ^ 2C3 1 1=b C 1 3=b 1=b h i 0 0 ^ ^ ^ 13 3C 1 2=b 3C 1 1=bC 1 1=b

where C0 is the ®rst derivative of the gamma function. Therefore, starting with an initial estimate obtained from Eq. (11) the recursive equation to estimate b^ in the iteration i 1 is

289

b^i1 b^i

G b^i =G0 b^i

The calculations are repeated until satisfying the error criterion

b^ i1 b^i 0

17b

where xj is the order statistic such that x1 x2 xN . From Eq. (16) the ®rst three population PWMs are

A0 f aC 1 1=b h i A1 f a2 1=b C 1 1=b 2 A2 f a3

1=b

C 1 1=b=3

18 19 20

By substituting these three population PWMs by the corresponding sample PWMs, A^0 , A^1 , and A^2 , the PWM estimator of the shape parameter b is a solution of

1 1

3

1=b^

2

1=b^

3A^2 2A^1

A^0 A^0

21

^ By using the Newton±Raphson Equation (21) can be solved numerically for b. method the following equations are de®ned

^ G b

1 1

3

1=b^

2

1=b^

3A^2 2A^1

and

^ G0 b

^2

b 1

1 2

A^0 A^0

h 2 3

22

1=b^

1=b^

1

2

1=b^

ln 3 2

1=b^

1

3

1=b^

291

i ln 2 23

Then, b^ can be determined following a similar procedure as that described in Sect. 3.1 for MOM estimation. Also, in this case, the value of b^ from Eq. (11) can ^ Then, the PWM estimators of the be used as an initial value to solve Eq. (21) for b. parameters a and f may be obtained from Eqs. (18) and (19) as

a^ A^0

2A^1

.h 1

i ^ C 1 1=b

1=b^

2

24

and

f^ A^0

^ a^C 1 1=b

25

In addition, note that using the relationship between the Weibull and the GEV models as indicated before in Sect. 2, i.e. b 1=b0 , the shape parameter can be obtained approximately from the formula suggested by Hosking et al. (1985). Furthermore, the PWM estimators for a two-parameter Weibull distribution (f 0) may be obtained from Eqs. (18) and (19) as

b^ ln 2 ln A^0 =A^1

ln 2

26

and

^ a^ A^0 C 1 1=b

27

3.3 Method of maximum likelihood (ML) The log-likelihood function of a three-parameter Weibull distribution is given by LL x; f; a; b N ln b

Nb ln a b

1

N X i1

ln xi

f

N X xi i1

f

b

a 28

The partial derivatives of LL with respect to f; a; and b are equated to zero yielding

292

N X

N bX xi f b 1 0 a i1 a

oLL of

b

oLL oa

N Nb b X xi f b 0 a a i1 a

oLL N ob b

1

xi

f

1

i1

N ln a

N X

ln xi

29a

29b

f

i1

N X xi i1

a

f b xi ln

f a

0

29c

respectively. These three equations must be solved simultaneously to ®nd the estimators of the parameters f, a, and b. Based on the Newton±Raphson method, one can write

2

3 2 Df o2 LL=of2 4 Da 5 4 o2 LL=oaof Db o2 LL=obof

3 3 12 oLL=of o2 LL=ofob o2 LL=oaob 5 4 oLL=oa 5 oLL=ob o2 LL=ob2

o2 LL=ofoa o2 LL=oa2 o2 LL=oboa

30

where 1 represents the inverse of the matrix and the second partial derivatives of the log-likelihood function of the Weibull distribution are given in Appendix A. One can evaluate the matrices of the right hand side of Eq. (30) for given values of the sample x1 ; . . . ; xN and parameters fi ; ai ; and bi and then obtain the increments Dfi ; Dai ; and Dbi in which the subscript i de®nes the iteration. Thus, the new parameter estimates for iteration i 1 are then determined by

ki1 ki Dki until satisfying the error criterion jDki =ki1 j < e in which k represents any of the parameters f; a; b, and e is a speci®ed relative error. The expected values of the second derivatives of the inverse matrix in Eq. (30) are given in Appendix B. These expected values are also used to obtain the asymptotic variance of the quantile estimator as shown in the next section. The inverse of the square matrix II 1 in Eq. (30) is called the inverse of the information matrix, which may be expressed as

2

II

1

6 6 1 6 6 6 ND 6 6 4

a2 b b 12 a2 h b b 1 abf b 1

a2 h b b 1 a2 a b2 ag

3 abf b 1 7 7 7 7 ag 7 7 7 5 2 bc

31

where

a C 1

2=b1 C00 2

b 1 C00 2

C2 1 2 C0 2 p2 6

c C 1

2=b

f C 1

1=b1

C2 1

1=b1 w 1

32a 32b

1=b

C0 2 w 1

1=b2

32c 1=b

32d

g C 1

2=bC0 2

C2 1

h C 1

1=bf1 C00 2

1=b1 w 1 C0 21 w 1

D bc f 2

1=b

32e

1=bg

32f 32g

00

in which C 2 is the second derivative of the gamma function with argument 2. Note that the coef®cients a, c, and g are de®ned only if b > 2. For a two-parameter Weibull distribution, f 0 in Eqs. (29b) and (29c). Combining these equations gives

N b

N X

Nb ln xi

i1

N P

xbi ln xi

i1 N P

i1

xbi

293

0

33

which is only a function of b. Solving Eq. (33) for b gives the maximum likelihood ^ and then a^ can be obtained from Eq. (29b). estimator b,

4 Confidence intervals on quantiles The 1 d con®dence interval X1 imated by (Kite, 1988) X^1

d

X^T u1

^

d=2 ST

d

on the population quantiles may be approx-

34

where u1 d=2 is the 1 d=2 quantile of the standard normal distribution, X^T is the quantile estimator corresponding to return period T, and S^T is the standard deviation or standard error of X^T . The quantile estimator can be written as

35 X^T g h^1 ; h^2 ; h^3 where h^i denotes in general estimators of either moments or parameters. The asymptotic variance of X^T may be expressed as oXT 2 oXT 2 oXT 2 2 ^ ^ ^ ST Var XT Var h1 Var h2 Var h^3 oh1 oh2 oh3 oXT oXT oXT oXT ^ ^ 2 Cov h1 ; h2 2 Cov h^2 ; h^3 oh1 oh2 oh2 oh3 oXT oXT Cov h^1 ; h^3 2 36 oh1 oh3 The terms in the right hand side of Eq. (36) may be determined depending on the estimation method. In this paper, such standard error can be estimated by the methods of moments, probability weighted moments, and maximum likelihood as described below.

4.1 Standard error by moments By using the ®rst three sample moments for ^ hi in Eq. (36), the variance of X^T can be written as (Kite, 1988)

S2T Var X^T l2 =N 1 KT c1 KT2 c2 oKT =oc1 2c2 2

oKT =oc1 c4

3c21

1=4

6 KT c3

3c1 c3

6c2

6c1 c2 =4 10c1 =4 2 9c1 c2 =4 35c21 =4 9 37

where the cj 's are the so called cumulants which are de®ned as 3=2

294

c1 l3 =l2

38a

c2 l4 =l22

38b

5=2

c3 l5 =l2

38c

c4 l6 =l32

38d

and c1 c is the skewness coef®cient. In turn, the rth central moments, lr are

l2 a2 D2

D21

39a

l3 a3 D3

3D2 D1 2D31

39b

l4 a4 D4

4D3 D1 6D2 D21

l5 a5 D5

5D4 D1 10D3 D21

10D2 D21 4D51

l6 a6 D6

6D5 D1 15D4 D21

20D3 D31 15D2 D41

3D41

39c 39d 5D61

39e

in which Dr C 1 r=b. In addition, the derivative of KT with respect to c1 can be written as

oKT =oc1 oKT =ob ob=oc1

40

The two partial derivatives in the right hand side of Eq. (40) can be determined as follows. From Eq. (10) the derivative of KT with respect to b is

D1 w1 oKT ob

ln BB1=b D2

D21

b2 D2

B1=b

D1

3=2

D2 w2 D21 w1

D21

where wr w 1 r=b C0 1 r=b=C 1 r=b and B from Eq. (5) the derivative of c1 with respect to b becomes

oc1 3 D2 ob D3

D21

2D31 w1 .h 2 D2 w2 D21 w1 b D2

D3 w3 2D2 D1 w2 D2 D1 w1

3D2 D1 2D31

41

ln 1=T. Likewise,

D21 5=2

i

42

4.2 Standard error by probability weighted moments The variance of the PWM quantile estimator X^T for the three-parameter Weibull distribution is obtained by replacing the h^i 's in Eq. (36) by the parameter estimators. It gives

S2T

oXT 2 oXT 2 oXT 2 ^ ^ Var f Var ^ a Var b of oa ob oXT oXT oXT oXT ^ ^ Cov f; a^ 2 Cov ^ a; b 2 of oa oa ob oXT oXT ^ b ^ 2 Cov f; of ob

43

From Eq. (8) the derivatives of XT with respect to the parameters f; a; and b are respectively

oXT 1 of

44a

oXT ln 1=T1=b oa

44b

oXT ob

44c

a ln ln 1=T ln 1=T1=b b2

^ Var(^ ^ Cov(f, ^ a^), Cov(^ ^ and Cov(f, ^ b) ^ in Eq. (43) The terms Var(f), a), Var(b), a, b), are given by (see Appendix C for details)

^ W 2 A00 2W0 W1 A01 W 2 A11 2W0 Wb C1 H Var f 0 1 2W1 Wb C2 H Wb2 CH 2

45a

Var ^ a T02 A00 2T0 T1 A01 T12 A11 2T0 Tb C1 H 2T1 Tb C2 H Tb2 CH 2 45b ^ CH 2 Var b

45c

Cov ^f; ^a W0 T0 A00 W0 T1 A01 W1 T0 A01 W1 T1 A11 W0 Tb C1 H W1 Tb C2 H Wb T0 C1 H Wb T1 C2 H Wb Tb CH 2

45d

^ b ^ W0 C1 H W1 C2 H Wb CH 2 Cov f;

45e

^ T0 C1 H T1 C2 H Tb CH 2 Cov ^ a; b

45f

and the terms in Eqs. (45a) through (45f) are also given in Appendix C. Finally, using B ln 1=T the asymptotic variance of X^T for the three-parameter Weibull distribution becomes

S2T

1 a2 ^ ^ B2=b Var ^ a 4 ln B2 B2=b Var b Var f N b h i 2a 1=b 1=b 1=b ^ ^ ^ ^ ln BB Cov f; b B Cov ^ a; b 2B Cov f; a^ b2

46

295

For a two-parameter Weibull distribution, the asymptotic distribution of the sample PWMs becomes bivariate normal by dropping A^2 and the third line and column of the square matrix in Eq. (C3). Following a similar procedure as described for the three-parameter Weibull model, the asymptotic variance has a simpler form as (Heo et al., 1990)

S2T 296

nh i XT2 1 2 21 1=b S 1 21 1=b 4H 1=2 S2 C 1 2=b NC 1 1=b h i o 1 2 211=b S 1 211=b S2 C2 1 1=b 2

47 where S fw 1 1=b

ln ln 1=Tg= ln 0:5.

4.3 Standard error by maximum likelihood The variance of the ML estimator of quantiles for the three-parameter Weibull distribution can be obtained from Eq. (43) where the parameter estimators are now ML estimators. The derivatives of XT with respect to the parameters f; a; and b are the same as Eqs. (44a) through (44c). On the other hand, the variance and covariance terms for ML estimators are the elements of the inverse of information matrix in Eq. (31) (Mood et al., 1974). Hence ^ Var f

a2 b b 12 ND

48a

Var ^ a

a2 a b2 ND

48b

^ Var b

b2 c ND

48c

^ a^ Cov f; ^ Cov ^ a; b

a2 h b b 1ND abf b 1ND

^ b ^ ag Cov f; ND

48d 48e 48f

where the various coef®cients such as a, b, etc. are given by Eq. (32). Thus, substituting Eqs. (44) and (48) into (43) yields

S2T

( a2 b ND b 12

2B1=b B1=b h f ln B 2 a b b 1 b

)

2

2g ln B ln B

49 in which B ln 1=T. Note that for the two-parameter Weibull distribution f 0 and Eq. (43) reduces to

S2T

oXT 2 oXT 2 oXT oXT ^ ^ Cov ^ a; b Var ^ a Var b 2 oa ob oa ob ( ) XT2 C0 2 ln ln 1 q2 1 p2 =6 Nb2

50

5 Simulation experiments Simulation experiments were performed to ®nd out the applicability of the derived con®dence intervals of quantiles based on the methods of moments, probability weighted moments, and maximum likelihood. For this purpose, the location and scale parameters were set to f 0 and a 1, respectively, and the shape parameter b varied as 7.49, 3.62, 2.23, 1.56, and 0.74 which correspond to skewness coef®cients equal to 0:5, 0.0, 0.5, 1.0, and 3.0, respectively. Using a known parameter set (f; a; b), 10,000 sets of data were generated for sample sizes N = 10, 25, 50, and 100. For each generated data set and estimation method, the parameters were estimated and then the con®dence limits X^1 d of Eq. (34) were determined by using the estimates X^T and S^2T obtained from the generated data. These estimates were obtained respectively from Eq. (8) and from either Eqs. (37), (46), or (49) depending on the estimation method. In addition, the theoretical con®dence limits X1 d were obtained also from Eq. (34) but XT and ST were determined based on the assumed (known) parameter set. Then the relative bias RBIASq= BIAS X^1 d =X1 d and the relative root mean square error RRMSE = MSE X^1 d =X1 d were determined. The computations were performed for both the lower and upper con®dence limits. Only the results for the upper con®dence limits for c 0:5, 0.5 and 1.0 are displayed in Figs. 2±5. For N = 10, the RBIAS of the upper con®dence limit generally increases as c increases as shown in Fig. 2. Also it generally increases as T increases except for

Fig. 2. RBIAS vs. return period for the upper con®dence limits for N 10

297

298

Fig. 3. RBIAS vs. sample size for the upper con®dence limits for T 100

Fig. 4. RRMSE vs. return period for the upper con®dence limits for N 10

ML and c 1. MOM gives the smallest RBIAS for skewness coef®cient c 0:5; 0:5 while ML gives the smallest values for c 1. However, RBIAS of PWM rapidly increases as c increases. Furthermore, for T 100 (q 0:99) RBIAS of the upper con®dence limit decreases as the sample size N increases as shown in Fig. 3. RBIAS are about the same for the three estimation methods for c 0:5, MOM gives the smallest RBIAS for c 0:5, while ML gives the smallest value for c 1. In addition, for the PWM method RBIAS becomes consistently

299

Fig. 5. RRMSE vs. sample size for the upper con®dence limits for T 100

larger as c increases. Therefore, the results show that in terms of bias MOM is preferable for negative c while ML is best for positive c. Relative root mean square error (RRMSE) for the upper con®dence limit is shown in Fig. 4 for N 10. RRMSE generally increases as c increases except for the case of MOM for c 0:5. Also RRMSE generally increases as T increases. ML always gives the smallest RRMSE for all ranges of c and PWM gives the largest RRMSE for positive c. For T 100, RRMSE decreases as N increases as shown in Fig. 5. Also RRMSE generally increases as the skewness c increases except for MOM for c 0:5 and N 10. For c 0:5 and N 10, Figs. 2 and 3 show that MOM gives the smallest bias. Yet Figs. 4 and 5 show a large value of RRMSE. This occurred because of the large sample variance of the upper con®dence limit obtained for that particular case. Results for the bias and mean square error for the lower con®dence limits are not shown but can be summarized as follows. For N 10 the bias generally increases as T increases. Also the bias increases as c increases. Generally the biases are small but they can be large for c 1 (e.g. RBIAS is about 25% for the MOM and PWM methods). For T 100 the biases decrease as N increases as expected. The PWM method gives the smallest bias for c 0:5 but as c increases the ML gives the smallest bias. Regarding the mean square error, for N 10 RRMSE generally increases as T increases; this is specially noticeable for c 1. For T 100, RRMSE increases as c increases and it decreases as N increases as expected. In all cases the ML gives the smallest RRMSE. Overall, the ML method for estimating the con®dence limits performs better than the other two methods in terms of bias and mean square error. This is specially so for c 0:5 even for small sample sizes (e.g. N 10).

6 Applications The annual 7-day low ¯ows (cms) of the Parana River at Corrientes, Argentina for the period (1904±1982) and the annual 6-h maximum rainfall data (mm) of Seoul,

Table 2. Parameter estimates of the Weibull distribution for (a) annual 7-day low ¯ows of the Parana River at Corrientes, Argentine, and (b) 6-h maximum annual rainfall data at Seoul, South Korea Methods

Parameters Low ¯ows

300

Rainfall data

f^

a^

b^

f^

MOM

3367.15

6509.15

2.292

33.25

74.70

1.711

ML

3545.67

6313.66

2.238

38.19

68.58 112.58

1.568 2.669

PWM

3562.17

6290.95

2.177

37.90

68.92

1.553

a^

b^

2-Parameter Weibull

South Korea for the period (1925±1987) are used here to illustrate and compare the estimation of the parameters, quantiles, and con®dence intervals assuming the Weibull distribution and the various estimation methods presented herein. The parameters and quantiles for various return periods were estimated by the methods of moments, probability weighted moments, and maximum likelihood as outlined in the previous sections. Table 2 shows the estimated parameters obtained for the referred data. Note that for the rainfall data, the parameters of the 2-parameter Weibull model based on the ML method are also included because in estimating the con®dence intervals using the ML method, the estimated shape parameter should be greater than 2 because the argument of C 1 2=b in Eq. (B1) of Appendix B should be greater than zero. Thus, while the parameter estimates for the 3-parameter Weibull model can be obtained by the ML method, the con®dence intervals cannot be determined because b^ 1:568 < 2. Thus, for the rainfall data, the 2-parameter Weibull model was used for obtaining con®dence intervals of quantiles based on the ML method. The empirical and ®tted cdfs based on the three-parameter Weibull for the two sets of data analyzed are displayed in Figs. 6 and 7. The Weibull cdfs appear to ®t the empirical cdfs quite well. Furthermore, goodness of ®t tests such as the chisquare test and Kolmogorov±Smirnov test did not reject the Weibull model with 5% signi®cance level. The quantiles and 95% con®dence limits, corresponding to return periods T 1:01, 1.05, 1.10, 1.5, 2, 5, 10, 20, 50, 100, and 500 were determined. Table 3 shows the results for the annual 7-day low ¯ows. The estimated quantiles by the three methods are not the same but the differences among them are not too large. Generally they are less than 5% for T 2 (percentage relative to the smaller estimate when comparing between the estimates of any two methods). Somewhat larger differences are observed in the estimated con®dence limits. The maximum differences for both the upper and lower limits are less than 8% for T 500. More notable are the differences in the con®dence intervals as shown for instance in Fig. 8. The ratios of the con®dence interval to the estimated quantile (for a given T) obtained based on the ML method are the narrowest being about 17% for T 2, while the ratios obtained based on the MOM and PWM methods are wider, they are about 30% for T 2. For the rainfall data, the differences of quantiles based on the MOM and PWM methods (for the three-parameter Weibull) are less than 10%. Similar differences are obtained between the estimated con®dence limits by both estimation methods except for the lower con®dence limits for T < 1:11 year. Furthermore, the

301

Fig. 6. Empirical and ®tted cdfs of the 3-parameter Weibull distribution for the annual 7-day low ¯ows of the Parana River at Corrientes, Argentine

Fig. 7. Empirical and ®tted 3-parameter Weibull distribution cdfs for the annual 6-h maximum rainfall data of Seoul, South Korea

con®dence intervals for the three-parameter Weibull model based on the ML method could not be determined because the estimated shape parameter was less than 2. Therefore, for comparison the two-parameter Weibull model was applied based on the ML method (the estimated parameters are given in Table 2). Con®dence intervals of quantiles for the two- and three-parameter Weibull models are shown in Fig. 9. As expected, large differences are observed for the upper con®dence limits obtained for the two- and three-parameter Weibull models.

7 Summary and conclusions Estimation techniques for determining the con®dence intervals for the two- and three-parameter Weibull distributions are presented based on the MOM, PWM, and ML. All three estimation methods require an iterative or numerical solution to estimate the shape parameter. The asymptotic variances of the MOM, PWM, and ML quantile estimators for the Weibull distribution are derived as a function of the sample size, return period, and parameters. Such variances can be used for

Table 3. The estimates of quantiles and 95% con®dence intervals based on the methods of moments, probability weighted moments, and maximum likelihood (Newton±Raphson) for the annual 7-day low ¯ows of the Parana River at Corrientes, Argentine Method

Return period T

Lower con®dence limit

Quantile ^T X

Upper con®dence limit

MOM

1.01 1.05 1.10 1.50 2.00 5.00 10.00 20.00 50.00 100.00 500.00

3005.6 4340.9 5076.1 7113.4 8222.2 10595.7 11780.1 12667.3 13567.6 14117.0 15123.8

4241.6 5109.6 5700.9 7757.0 8914.4 11378.6 12733.8 13873.5 15171.1 16042.0 17813.0

5477.5 5878.2 6325.7 8400.6 9606.4 12161.4 13687.4 15079.8 16774.6 17966.9 20502.2

PWM

1.01 1.05 1.10 1.50 2.00 5.00 10.00 20.00 50.00 100.00 500.00

3358.3 4474.0 5117.2 7084.5 8199.8 10585.9 11799.4 12736.2 13717.4 14332.2 15492.1

4322.4 5133.2 5699.0 7717.7 8878.3 11390.3 12790.2 13975.9 15334.1 16250.1 18123.0

5286.6 5792.3 6280.7 8350.9 9556.8 12194.6 13780.9 15215.6 16950.8 18168.0 20753.9

ML

1.01 1.05 1.10 1.50 2.00 5.00 10.00 20.00 50.00 100.00 500.00

4056.0 4842.4 5382.3 7333.1 8457.8 10805.0 12033.1 13030.4 14132.7 14856.1 16293.8

4353.8 5182.9 5754.1 7763.4 8905.5 11355.5 12711.0 13854.8 15160.5 16038.9 17829.5

4651.6 5523.5 6125.9 8193.8 9353.2 11906.0 13389.0 14679.2 16188.4 17221.6 19365.2

302

estimating the con®dence limits and con®dence intervals of the population quantiles. Except for the two-parameter Weibull model, the formulas obtained do not have simple forms but can be evaluated numerically. Simulation experiments were performed to ®nd out the applicability of the derived con®dence intervals of quantiles based on the three estimation methods. The interest was to see the applicability of the estimations methods for various sample sizes, return periods, and skewness coef®cients. The results show that overall, the ML method for estimating the con®dence limits performs better than the other two methods in terms of bias and mean square error. This is specially so for c 0:5 even for small sample sizes (e.g. N 10). However, the drawback of the ML method for determining the con®dence limits is that it requires that the shape parameter be bigger than 2. Finally, the Weibull model was applied to ®t the distribution of annual 7-day low ¯ows and 6-h maximum annual rainfall data. The purpose was to illustrate

303

Fig. 8. Cumulative distribution of annual 7-day low ¯ows of the Parana River at Corrientes, Argentine and 95% con®dence intervals based on the Weibull-3 model and the MOM, PWM, and ML estimation methods

the applicability of the estimation methods to real data and to compare the values of estimated parameters, quantiles, and con®dence intervals based on the MOM, ML, and PWM methods. The results showed that the differences in the estimated quantiles based on the three methods are not large, generally are less than 10%. However, the differences between the con®dence limits and con®dence intervals obtained by the three estimation methods may be more signi®cant. For instance, for the 7-day low ¯ows the ratio between the estimated con®dence interval to the estimated quantile based on ML is about 17% for T 2 while it is about 30% for estimation based on MOM and PWM methods. The analysis of the rainfall data also illustrated the need of having available alternative estimation methods. For example, the analysis of the rainfall data using the three-parameter Weibull showed that while ML parameters can be estimated, the corresponding con®dence limits and intervals could not be found because the shape parameter was smaller than 2.

Appendix A: Second partial derivatives of the log-likelihood function for the Weibull distribution X b 2i hX 2 yi o2 LL=of2 b 1=a2 yi b o2 LL=oaof b=a2 o2 LL=ofob

X xi

X

yib

f

1

A1

1

A2 b=a

X

yib

1

ln yi

1=a

X

yib

1

A3

304

Fig. 9. Cumulative distribution of the annual 6-h maximum rainfall data of Seoul, Korea and 95% con®dence intervals based on the Weibull-2 and -3 models and the MOM, PWM, and ML estimation method

o2 LL=oa2

h b=a2 N

b 1

h o2 LL=oaob 1=a N o2 LL=ob2 N=b2 where yi xi

X

X

f=a and

ybi

X

b

ybi

X

i

i ybi ln yi

ybi ln yi 2

P

A4 A5 A6

represents summation from 1 to N.

Appendix B: Expected values of the second partial derivatives of the log-likelihood function for the Weibull distribution

o2 LL N b 12 E C 1 2=b a2 of2 2 o LL Nb2 E a2 oa2 2 o LL N 2 1 C00 2 E 2 ob b 2 o LL Nb2 C 2 1=b E a2 ofoa

B1 B2 B3 B4

o2 LL E ofob 2 o LL E oaob

N 1 a

1=bC 1

1=b1 w 1

1=b

N 0 C 2 a

B5 B6

where C and w are gamma and digamma functions and C0 2 and C00 2 are the ®rst and second partial derivatives of the gamma function with argument 2, respectively.

Appendix C: Derivation of the variance of the quantile estimator ST2 based on PWM method The quantile estimator of X^T can be written as a function of the ®rst three sample PWMs X^T f A^0 ; A^1 ; A^2

C1

where the sample PWM A^r is de®ned in Eq. (17). The variance S2T can be obtained by

S2T

oXT 2 oXT 2 oXT 2 Var A^0 Var A^1 Var A^2 oA0 oA1 oA2 oXT oXT oXT oXT ^ ^ 2 Cov A0 ; A1 2 Cov A^1 ; A^2 oA0 oA1 oA1 oA2 oXT oXT 2 Cov A^0 ; A^2 oA0 oA2

C2

However, oXT =oA1 in Eq. (C2) is not obtainable for some probability distributions such as the Weibull model, because X^T cannot be written explicitly as a function of the sample PWMs. In other words, the quantile estimator in Eq. (8) cannot be explicitly written as a function of the sample PWMs because the shape parameter estimator is implicitly expressed as a function of the sample PWMs as shown in Eq. (21). Therefore, in this case, the variance±covariance matrix of the sample PWMs should be obtained ®rstly. Then this matrix is transformed to the variance±covariance matrix of the parameter estimators by using a Jacobian transformation. Finally, the variance of the quantile estimator can be obtained by using Eq. (43). Details of the procedure follows. The asymptotic distribution of the sample PWMs can be written as (Chernoff et al., 1967; Hosking et al., 1985)

2

3 0 2 A0 A00 =N A^0 4 A^1 5 [email protected] A1 ; 4 A01 =N A2 A02 =N A^2

A01 =N A11 =N A12 =N

31 A02 =N A12 =N 5A A22 =N

C3

where reads ``is asymptotically distributed as'', TVN is an abbreviation for trivariate normal distribution, and Aij are given by (Heo et al., 1990)

A00 a2 C 1 2=b

C2 1 1=b

C4a

305

h A01 a2 =2 2

2=b

nh A02 a2 =2 3

C 1 2=b 1

2=b

2

2=b

H 1=2

C 1 2=b A11 a2 2

306

21

1=b

1=b

2

i C2 1 1=b

C4b

i

2 3

1=b

o C2 1 1=b

H 1=2C 1 2=b C2 1 1=b h A12 a2 =2 3 2=b H 1=3C 1 2=b 2:6 1=b

C4c

2=b

C4d 2

1=b

C2 1 1=b

i C4e

A22 a2 3

2=b

H 2=3C 1 2=b

C2 1 1=b

C4f

where H z is a hypergeometric function de®ned by

H z F 2=b; 1=b; 1

1=b; z

1 1=b X C n 2=b zn C 2=b n0 n 1=bn!

C5

For the three-parameter Weibull distribution, the asymptotic variance of the PWM estimator of quantile, X^T can be found by using successively the following transformations

A^0 A^1 A^2

A^0 ) A^1 A^2 R

A^0 ) A^1 A^2 b^

f^ ) a^ ) X^T b^

C6

where ) reads ``is transformed as''. As mentioned for Eq. (21), the shape parameter estimator is a function of the sample PWMs, but cannot be written explicitly, thus we need a transformation. In the ®rst transformation, let R 3A^2 A^0 = 2A^1 A^0 which is the right hand side of Eq. (21). Then b^ is ^ ^ given implicitly by 1 3 1=b = 1 2 1b R in the second transformation. Additional details of the transformations follows. (1) 1st Transformation. For the 1st transformation, the Jacobian matrix is given by

2 J1

oA^0 6 oA^^0 6 oA1 6 oA^0 6 oA^ 6 2 4 oA^0 oR oA^0

oA^0 oA^1 oA^1 oA^1 oA^2 oA^1 oR oA^1

oA^0 oA^2 oA^1 oA^2 oA^2 oA^2 oR oA^2

3

2 7 7 6 7 6 74 7 5

1 0 0

3A2 2A1 2A1 A0 2

0 1 0

2 3A2 A0 2A1 A0 2

0 0 1

3 2A1 A0 2A1 A0 2

3 7 7 5

C7

Then the variance±covariance matrix in the 1st transformation, i.e. of the 2nd column terms in (C6) can be obtained by R1 J1 R0 J10 where R0 is the variance± covariance matrix in Eq. (C3). Thus

2

A00 =N 6 A00 =N R1 6 4 A00 =N C1 =N

A01 =N A01 =N A01 =N C2 =N

3 C1 =N C2 =N 7 7 C3 =N 5 C=N

A02 =N A02 =N A02 =N C3 =N

C8

where the elements are de®ned in Eqs. (C4a) through (C4f) and

h C A00 3

1=b

2

1=b

2A01 3

1=b

1 3A02 2

i. 1 M

1=b

C9a h C1 A01 3 h C2 A02 3

1=b

2

1=b

1=b

1=b

2

2A11 3

1 3A12 2

1=b

2A12 3

1=b

1=b

1 3A22 2

i. 1 M

C9b

i. 1 M

1=b

C9c M aC 1 1=b 1

2

1=b

2

C9d

(2) 2nd Transformation. In this case, the Jacobian matrix is

2 oA^

0

J2

^ 6 ooAA^0 6 1 6 oA^0 6 oA^2 6 ^ 4 oA0 ob^ oA^0

oA^0 oA^1 oA^1 oA^1 oA^2 oA^1 ob^ oA^1

where R 1

ob^ oR ln 22

oA^0 oA^2 oA^1 oA^2 oA^2 oA^2 ob^ oA^2

3

1=b 1

1=b

oA^0 oR oA^1 oR oA^2 oR ob^ oR

= 1

3

2 1 7 7 6 7 60 74 0 7 5 0 2

1=b

3 0 0 0 0 7 7 1 0 5 ^ 0 ob=oR

0 1 0 0

C10

and

2 b2 1 2 1=b 3 1=b ln 33

1=b 1

2

C11

1=b

Thus, the variance±covariance matrix in the 2nd transformation is given by

2

A00 =N 6 A00 =N R2 J2 R1 J20 6 4 A00 =N C1 H=N

A01 =N A01 =N A01 =N C2 H=N

A02 =N A02 =N A02 =N C3 H=N

3 C1 H=N C2 H=N 7 7 C3 H=N 5 CH 2 =N

C12

(3) 3rd Transformation. The Jacobian matrix in the 3rd transformation is

2

^ A^0 of=o 4 J3 o^ a=oA^0 0

^ A^1 of=o o^ a=oA^0 0

^ A^2 of=o o^ a=oA^0 0

3 2 ^ b^ W0 of=o o^ a=ob^ 5 4 T0 0 1

W1 T1 0

3 0 Wb 0 Tb 5 0 1 C13

307

where the elements of above matrix are given by

W0 of=oA0

2

1=b

. 1

. W1 of=oA1 2 1 Wb of=ob 308

2

Tb oa=ob h 1 2

1=b

1=b

1=b

.h 2 1

1=b

2

2

C14b 1=b

.h b2 1

2

1=b

i

C14c

i C 1 1=b

1=b

w 1 1=b

C14a

aC 1 1=b ln 22

.h T0 oa=oA0 1 1 T1 oa=oA1

2

C14d

i C 1 1=b

ln 22

1=b

i b2 1

C14e

2

1=b

2

C 1 1=b C14f

Then the variance±covariance matrix can be obtained by

2

Z 1 4 ff 0 Zfa R3 J3 R2 J3 N Z fb

Zfa Zaa Zab

3 Zfb Zab 5 Zbb

C15

where the diagonal and off-diagonal terms are the variances and covariances of estimators, respectively, given by Eqs. (45a) through (45f). (4) 4th Transformation. The quantile estimator of the 3-parameter Weibull distribution for return period T is ^ X^T f^ a^ ln 1=T1=b

C16

and the Jacobian matrix is given by

"

oX^T oX^T oX^T J4 a ob^ of^ o^

#

where each element is given in Eqs. (44a) through (44c). Finally, the S2T of X^T can be obtained by

S2T J4 R3 J40

C17

and the result is given in Eq. (49).

References

Boes DC, Heo JH, Salas JD (1989) Regional ¯ood quantile estimation for a Weibull model. Water Resour. Res., 25: 979±990 Chernoff H, Gastwirth JL, Johns MV (1967) Asymptoti distribution of linear combinations of functions of order statistics with applications to estimation. Annals Math. Statist., 38: 52±72

Chow VT (1951) A general formula for hydrologic frequency analysis. Trans. Geophys. Union, 32: 231±237 Cohen AC (1965) Maximum likelihood estimation in the Weibull distribution based on complete and on censored samples. Techno., 7: 579±588 Cohen AC, Whitten BJ, Ding Y (1984) Modi®ed moment estimation for the three-parameter Weibull distribution. J. Qual. Tech., 16: 159±167 Dodson B (1994) Weibull Analysis. ASQC Press, Wisconsin Greenwood JA, Landwehr JM, Matalas NC, Wallis JR (1979) Probability weighted moments: de®nition and relation to parameters of several distributions expressible in inverse form. Water Resour. Res., 15: 1049±1054 Gumbel EJ (1958) Statistic of Extreme. Columbia University Press, NY Haan CT, Beer CE (1967) Determination of maximum likelihood estimators for the three parameter Weibull distribution. Iowa State J. Science, 42: 37±42. Iowa State University Ames, Iowa Harter HL, Moore AH (1965) Maximum likelihood estimation of the parameters of gamma and Weibull populations from complete and from censored samples. Techno., 7: 639±643 Harter HL, Moore AH (1967) Asymptotic variances and covariances of maximum likelihood estimators, from censored samples, of the parameters of Weibull and gamma populations. Annals Math. Statis., 38: 557±570 Heo JH, Boes DC, Salas JD (1990) Regional ¯ood frequency modeling and estimation. Water Resour., Paper No. 101, Colorado State University, Fort Collins, Colorado Hosking JRM, Wallis JR, Wood EF (1985) Estimation of the generalized extreme-value distribution by the method of probability weighted moments. Techno., 27: 251±261 Jenkinson AF (1969) Statistics of extremes. In: Estimation of Maximum Floods. Technical Report No. 98, WMO Johnson NL, Kotz S (1970) Continuous Univariate Distribution-1. Houghton Mif¯in Company, Boston Kite GW (1988) Frequency and Risk Analysis in Hydrology. Water Resources Publications, Littleton, Colorado Landwehr JM, Matalas NC (1979) Probability weighted moments compared with some traditional techniques in estimating Gumbel parameters and quantiles. Water Resour. Res., 15: 1055±1064 Lemon GH (1975) Maximum likelihood estimation for the three parameter Weibull distribution based on censored samples. Techno., 17: 247±254 Matalas NC (1963) Probability Distribution of Low Flows. USGS Prof. Paper. 434-A Mood AM, Graybill FA, Boes DC (1974) Introduction to the Theory of Statistics, (3rd Edn). McGraw-Hill, NY National Environmental Research Council (1975) Flood Studies Report. vol. 1, London, UK Prescott P, Walden AT (1983) Maximum likelihood estimation of the three-parameter generalized extreme-value distribution from censored samples. J. Statis. Comput. Simul., 16: 241±250 Rao DV (1981) Three-parameter probability distributions. J. Hydraul. Div. Am. Soc. Civ. Eng., 109: 339±358 Thoman DR, Bain LJ (1969) Inferences on the parameters of the Weibull distribution. Techno., 11: 445±460 Weibull W (1939) A statistical theory of the strength of materials. Ing. Vetenskaf Akad. Handl. 151 Weibull W (1951) A statistical distribution function of wide applicability. J. Appl. Mech., Am. Soc. Mech. Eng., 18: 293±297

309