Use of parameter transformations in nonlinear, discrete flood event models

Journal of Hydrology, 117 (1990) 55-79 Elsevier Science Publishers B.V., Amsterdam - - Printed in The Netherlands

55

I4] U S E OF P A R A M E T E R T R A N S F O R M A T I O N S IN N O N L I N E A R , D I S C R E T E FLOOD E V E N T M O D E L S

BRYSON C. BATES Division of Water Resources, CSIRO, Wembley, W.A. 6014 (Australia) (Received September 5, 1989; accepted for publication September 19, 1989)

ABSTRACT Bates. B.C., 1990. Use of parameter transformations in nonlinear, discrete flood event models. J. Hydrol., 117: 55-79. The use of parameter transformations in calibration and post-calibration studies of conceptual rainfall-runoff models is often justified on the basis of their ability to speed the convergence of automatic parameter estimation algorithms, improve the sampling properties of the parameter estimates, and improve statistical inferences that are based on conventional linear model theory. An examination of the abilities of parameter transformations was undertaken using a nonlinear flood event model, rainfall runoff data from six Australian catchments, a Gauss-Newton optimization algorithm, three measures of statistical nonlinearity, and three types of joint confidence regions (exact, likelihood and linear theory) for the model parameters. Results indicate that the 'expected-value' transformation of Ross and the power transformation of Tsai can lead to improved sampling properties and statistical inferences on model parameters. However, it is clear that the effect of model reparameterization on the convergence of optimization algorithms may only be slight, and that care must be taken when applying parameter transformations to nonlinear flood event models.

INTRODUCTION

The application of statistical parameter estimation techniques to hydrological models has gained considerable popularity in recent years. The literature contains many applications of techniques and descriptions of difficulties encountered during the estimation process (e.g. Ibbitt and O'Donnell, 1974; J o h n s t o n a n d P i l g r i m , 1976; S o r o o s h i a n a n d D r a c u p , 1980; R e s t r e p o - P o s a d a a n d B r a s , 1982; K u c z e r a , 1983; W i l l i a m s a n d Yeh, 1983; B a t e s a n d T o w n l e y , 1988). R e c e n t w o r k b y K u c z e r a (1983), G u p t a a n d S o r o o s h i a n (1983), S o r o o s h i a n a n d G u p t a (1985), C a r r e r a a n d N e u m a n (1986), a n d B a t e s (1988) h a s i n d i c a t e d t h e d e s i r a b i l i t y of u s i n g m o d e l p a r a m e t e r i z a t i o n s w h i c h f a c i l i t a t e t h e c o n v e r g e n c e of p a r a m e t e r e s t i m a t i o n a l g o r i t h m s a n d l e a d t o i m p r o v e m e n t s i n p a r a m e t e r i n f e r e n c e . W h i l e t h e r e s u l t s of t h i s w o r k h a v e b e e n e n c o u r a g i n g , m o s t of t h e s e s t u d i e s h a v e i n v o l v e d a v e r y s m a l l n u m b e r o f s y n t h e t i c o r o b s e r v e d d a t a sets, c o n s i d e r a t i o n of o n l y l o g a r i t h m i c t r a n s f o r m a t i o n s of m o d e l

0022-1694/90/$03.50

© 1990 Elsevier Science Publishers B.V.

56

B C. BATES

NOTATION

f(x,. 4 F@, n g IN

k m n P

PE s U u, l-i ”

V. V.. XI a

2

0

P,

a)

model response function upper a probability point of F distribution with p and n - p degrees of freedom standardized asymptotic third moment of Hougaard (1985) maximum intrinsic curvature parameters to be estimated in case study model number of storm events number of model parameters maximum parameter-effects curvature ‘expected-value’ parameters vector @ x 1) parameter vector @ x 1) starting vector for Gauss-Newton algorithm @ x 1) parameter estimates vector (p x 1) transformed parameters vector Cp x l), eqns. (4H41); power transformed parameters vector after eqn. (41) matrix of first partial derivatives of the model response function with respect to the model parameters (n x p) matrix of second partial derivatives of the model response function with respect to the model parameters (n x p x p) known input vector level of significance random error variance zero vector

parameters, or of transformations that are specific to the given model equations. In many practical applications, there is little available guidance as to which parameter transformations improve the speed of convergence of optimization algorithms or lead to improved parameter inferences. Frequently, the analyst must experiment with many transformations without any guarantee of success. It is the purpose of this paper to build on previous work by examining the application of two general classes of parameter transformations to nonlinear, discrete flood event models and observed rainfall-runoff data. Of particular interest is whether a priori and a posteriori transformations can lead to faster convergence and improved parameter inferences, respectively. The paper begins by setting the parameter estimation and inference problems within a nonlinear regression framework. This is followed by reviews of: the theory on the application of parameter transformations to nonlinear regression models; the measures of statistical nonlinearity used to investigate the sampling properties of the original and transformed parameter estimates; the GaussNewton algorithm used to assess the effect of parameter transformations on the estimation process; and the joint confidence regions for model parameters that are used to assess how well the linear theory confidence regions for the transformed parameters approximate exact regions for the original parameters. A case study involving a nonlinear flood event model and rainfall--runoff data from six Australian catchments is then used to investigate the utility of

57

NONLINEAR, DISCRETE FLOOD EVENT MODELS

p a r a m e t e r t r a n s f o r m a t i o n s . T h r o u g h o u t the paper, p a r t i c u l a r emphasis is placed on the p e r f o r m a n c e of one of the m e a s u r e s of statistical n o n l i n e a r i t y devised by H o u g a a r d (1985), the 'expected-value' p a r a m e t e r t r a n s f o r m a t i o n of Ross (1970), and the p a r a m e t e r power t r a n s f o r m a t i o n of Tsai (1988). BACKGROUND Parameter estimation and inference Suppose t h a t a p p a r a m e t e r , n o n l i n e a r flood e v e n t model is to be c a l i b r a t e d by the fitting of c o m p u t e d and observed h y d r o g r a p h peaks for n discrete storm events. The model m a y be r e p r e s e n t e d by the n o n l i n e a r r e g r e s s i o n equations: qt = f(xt, u) + ~t,

t =

1. . . . .

n

(1)

w h e r e xt is the k n o w n i n p u t v e c t o r (whose c o m p o n e n t s are the o r d i n a t e s of the rainfall or rainfall excess h y e t o g r a p h ) associated with the tth observed peak d i s c h a r g e qt, u is the (p × 1) v e c t o r of u n k n o w n parameters, f is the model response f u n c t i o n w h i c h is assumed to be twice c o n t i n u o u s l y differentiable in u, and the et are r a n d o m errors. If the model response f u n c t i o n is l i n e a r in the p a r a m e t e r s then: P

f(xt, u)

= x~u =

~ x,juj j 1

(2)

where the prime denotes the t r a n s p o s e of a v e c t o r or m a t r i x and xtj denotes the j t h c o m p o n e n t of xt. The most p o p u l a r statistical a p p r o a c h to p a r a m e t e r e s t i m a t i o n in the h y d r o l o g i c l i t e r a t u r e is the m e t h o d of least squares. The least squares estimates fi minimize the sum of squares function: S(u)

=

[q - f(u)]' [q - f(u)]

=

~ [q, - fix,, u)] 2

(3)

t=l

w h e r e q and f(u) d e n o t e c o l u m n vectors whose tth c o m p o n e n t s are qt and f(xt, u), respectively. F o r models t h a t are l i n e a r in the parameters, the c o n t o u r s of the response surface defined by eqn. (3) form hyperellipsoids and the surface is said to be quadratic. It is possible to obtain unbiased estimates of the p a r a m e t e r s and t h e i r variances, and to m a k e inferences based on these estimates, provided t h r e e assumptions are satisfied: (1) the model response function is c o r r e c t l y specified; (2) the c o m p o n e n t s of the input v e c t o r are n o n s t o c h a s t i c and are m e a s u r e d w i t h o u t error; (3) the errors ~t are i n d e p e n d e n t identically distributed n o r m a l r a n d o m variables with m e a n zero and u n k n o w n c o n s t a n t v a r i a n c e a2 [NID (0, a2)]. If these assumptions hold, the least squares e s t i m a t o r s of u are unbiased, n o r m a l l y distributed, and minimum v a r i a n c e estimators. F o r models t h a t are n o n l i n e a r in the p a r a m e t e r s , the above properties hold only ' a s y m p t o t i c a l l y ' (i.e. as n -~ ~ ) . F o r smaller samples, the response surface

58

B.C. BATES

is n o n q u a d r a t i c and the least squares e s t i m a t o r s m a y be badly biased, nonn o r m a l l y distributed, and h a v e v a r i a n c e s exceeding the m i n i m u m possible variance. M o r t o n (1987) notes t h a t in p r a c t i c e the most i m p o r t a n t f e a t u r e of the d i s t r i b u t i o n of a least squares e s t i m a t o r for a n o n l i n e a r model is its a s y m m e t r y a b o u t its mean. The e x t e n t of this b e h a v i o u r varies widely b e t w e e n different m o d e l / d a t a set combinations, and t h e r e are no guidelines as to how large a sample must be for l i n e a r i t y to be closely a p p r o x i m a t e d (Ratkowsky, 1983, sections 1 2).

Parameter transformations The form of the response surface and the sampling properties of the least squares e s t i m a t o r s are not i n v a r i a n t u n d e r p a r a m e t e r t r a n s f o r m a t i o n s . It may be possible to find a p a r a m e t e r i z a t i o n t h a t gives rise to response surface c o n t o u r s t h a t are more n e a r l y q u a d r a t i c and to estimators t h a t a p p r o x i m a t e the a s y m p t o t i c conditions. Such p a r a m e t e r i z a t i o n s simplify the problem of o b t a i n i n g j o i n t confidence regions for the model parameters. An a l t e r n a t i v e p a r a m e t e r i z a t i o n may be r e p r e s e n t e d by: v

=

T(u)

(4)

implying t h a t the new p a r a m e t e r s v are functions of the original p a r a m e t e r s only and not the input v e c t o r xt. Thus the r e p a r a m e t e r i z e d model gives the same predicted v a l u e for the observed response as the original model. Two general classes of p a r a m e t e r t r a n s f o r m a t i o n s are considered in this paper. The first is the stable p a r a m e t e r or 'expected-value' t r a n s f o r m a t i o n proposed by Ross (1970). H e r e the a n a l y s t chooses a set of i n p u t vectors x* . . . . . x* and uses the expected values vl . . . . . vp where: v~ = f ( x * , u ) ,

i =

1. . . . . p

(5)

as p a r a m e t e r s a f t e r eliminating the original p a r a m e t e r s u from t h e p equations. This type of t r a n s f o r m a t i o n is illustrated in the case study. The second class of t r a n s f o r m a t i o n s is the power family suggested by Tsai (1988): vi = =

uli

ifTi #

0,

i =

1. . . . , p

ln(ui)

ifTi = 0,

i =

1. . . . . p

(6)

The a d v a n t a g e of these t r a n s f o r m a t i o n s is t h a t values for the c o n s t a n t s x* and 7i (i = 1 . . . . . p) can often be found which will improve the shape of the c o n t o u r s of the response surface and the sampling properties of least squares e s t i m a t o r s for a m o d e l / d a t a set combination. One d i s a d v a n t a g e of ~expectedvalue' t r a n s f o r m a t i o n s is the algebraic complexity of t h e i r expressions c o m p a r e d with t h a t of the original model p a r a m e t e r i z a t i o n . The t r a n s f o r m e d p a r a m e t e r s h a v e the t e n d e n c y to a p p e a r m o r e t h a n once in the r e p a r a m e t e r i z e d model response function, and the complexity of these f u n c t i o n s tends to i n c r e a s e as the n u m b e r of model p a r a m e t e r s is increased. A n o t h e r limitation is

59

NONLINEAR. DISCRETE FLOOD EVENT MODELS

that it may not always be possible to solve the equations defined by eqn. (5) to obtain expressions for the original parameters in terms of the transformed parameters alone. Nonlinearity

measures

Several measures of nonlinearity have been proposed for assessing the sampling properties of least squares estimators. These include the bias calculation of Box (1971), the asymptotic bias and the asymptotic third moment of Hougaard (1985), and the measures of intrinsic c u rvat ure and parameter-effects c u r v a tu r e of Bates and Watts (1980). The bias measures of Box and Hougaard are identical for the model specified by eqn. (1), (Cook et al., 1986), and relate to the bias in each individual parameter estimator. Let V. be the (n × p) first derivatives matrix with elements: (V.)ti = ?f(xt, f l ) / ~ u i ,

(t

=

1.....

n; i

=

1.....

p)

(7)

and V.. the (n × p x p) second derivatives matrix with elements: (V..)~ij = ~2f(xt, fi)/<~ui<)uj,

(j =

1. . . . . p)

(8)

Define J = V.'V. and W = V.'V.. where the ( i , j , k)th element of W is given by: W~jk =

~ (V.)ti(V..)t2~,

(k

=

1.....

(9)

p)

t

Hougaard's bias measure is defined as: E(l~li)-

1

Ui -

2

62

~_, WjkrJkr,

(r

=

1,...,p)

(10)

jkrJ~

where E(d~) denotes the expected value of ~ , a term like j i j is the (i, j)th element of the matrix J 1 and the estimate of a 2 is given by: ~2

=

S(fi)/(n

(11)

- p)

An indicator of nonlinear behaviour is the percentage bias, defined here as: %b~ =

1 0 0 ( E ( ~ ) - ui)/$t~,

(i

=

1 .....

p)

(12)

Ratkowsky (1983, sections 2-5) suggests t ha t absolute values of % b i in excess of 1% signal which parameter or parameters may be largely responsible for the overall nonlinearity. However, Ratkowsky (1986) also notes that %bi is location dependent since a high percentage bias will be obtained if u~ ~ 0. Another measure of the nonlinear behaviour of each individual parameter estimator is the asymptotic third moment of Hougaard (1985): E[/t i -

E(/t,)]

3

=

-

o'4~=

d i q g i s g i r ( W q s r + Wsq r + Wrqs)

(13)

qsr

(q = 1. . . . . p; s = 1. . . . . p). Ratkowsky (1990) proposes a standardization of this third moment to give a direct measure of the skewness of u~:

60

gi

B.C. BATES

=

E[tii-

E(tzi)]3/(a2jii) m ,

i

=

1. . . . .

p

(14)

and has found gi to be a v e r y useful i n d i c a t o r of n o n l i n e a r i t y . Intrinsic c u r v a t u r e refers to the shape of the solution locus described by: f(u)

=

[f(xi, u) . . . f(xn, u)]'

(15)

This c u r v a t u r e is defined by the m o d e l / d a t a set c o m b i n a t i o n and is u n a l t e r e d by a model r e p a r a m e t e r i z a t i o n . In a l i n e a r model, the solution locus is a h y p e r p l a n e and t h e r e is no intrinsic c u r v a t u r e . T h e g r e a t e r the intrinsic c u r v a t u r e in a n o n l i n e a r model, the g r e a t e r the expected bias in the model predictions and the residual sum of squares. Parameter-effects c u r v a t u r e refers to the b e h a v i o u r of the lines on the solution locus r e p r e s e n t i n g c o n s t a n t values of ui (i = 1 . . . . . p). This c u r v a t u r e is defined by the way in which the ui a p p e a r in the model function, and m a y often be r e d u c e d by a suitable reparameterization. F o r a l i n e a r model, the p a r a m e t e r lines are straight, parallel, and equidist a n t for equal i n c r e m e n t s of u~. If intrinsic and parameter-effects c u r v a t u r e s fall w i t h i n a c c e p t a b l e limits, the least squares estimate fi will be almost unbiased and n o r m a l l y distributed, and the v a r i a n c e of each c o m p o n e n t of will be close to the m i n i m u m possible v a r i a n c e (Ratkowsky, 1983, Sections 1 and 9). For the c u r v a t u r e m e a s u r e s of Bates and Watts (1980), it is assumed t h a t the elements of the V. and V.. m a t r i c e s h a v e been divided by the ' s t a n d a r d radius' 5x/p. C a l c u l a t i o n of the m e a s u r e s involves the (n × p × p) c u r v a t u r e array: A..

-

[Q'][L'V.. L]

(16)

and the Q R decomposition of V. V.

=

QR

(17)

in which Q is an (n × n) o r t h o g o n a l matrix, R is an (n x p) m a t r i x with zeros below the main diagonal, and L = l~-1 w h e r e ~ is the upper (p × p) s u b m a t r i x of R. The s q u a r e b r a c k e t m u l t i p l i c a t i o n in eqn. (16) indicates s u m m a t i o n o v e r the index t, and m u l t i p l i c a t i o n s w i t h i n the b r a c k e t s indicate s u m m a t i o n over index i for p r e m u l t i p l i c a t i o n or index j for postmultiplication. A.. is p a r t i t i o n e d such that: A..

=

APE IA IN

(18)

where A PE is the parameter-effects a c c e l e r a t i o n a r r a y consisting of the first p faces of A.., and A l~ is the intrinsic a c c e l e r a t i o n a r r a y consisting of the last n - p faces. The m a x i m u m intrinsic c u r v a t u r e (IN) and the m a x i m u m parameter-effects c u r v a t u r e (PE) are defined by: IN

= max I[d'A~ Nd [I

(19)

PE

= max I[d'APEd ][

(20)

where d is a u n i t v e c t o r and the double v e r t i c a l bars indicate the norm of a vector.

NONLINEAR, DISCRETE FLOOD EVENT MODELS

61

R a t k o w s k y (1983, p. 18) considers the b e h a v i o u r of a m o d e l / d a t a set combination to be ~close-to-linear' if:

IN

< 1/2~F(p, n - p, a)

(21)

<

(22)

and:

PE

1 / 2 ~ F ( p , n - p, a)

where F(p, n - p, a) denotes the upper a probability point of the F d i s t r i b u t i o n with p and n - p degrees of freedom. In the case study to follow, a v a l u e of a of 0.05 is used to c h e c k the acceptability of I N and PE. I N and P E values t h a t are slightly in excess of the guide values given by the r i g h t - h a n d sides of eqns. (21) and (22) will also be deemed acceptable. A l i m i t a t i o n of P E is t h a t it applies only to the full p a r a m e t e r v e c t o r and does not provide any way of assessing n o n l i n e a r i t y in the individual parameters. It has also been found t h a t P E is, on occasion, an u n r e l i a b l e i n d i c a t o r of n o n l i n e a r i t y (Bates and Watts, 1980; R a t k o w s k y , 1983; Cook and Witmer, 1985). One way of e v a l u a t i n g the n o n l i n e a r i t y m e a s u r e s for a r e p a r a m e t e r i z a t i o n is to express the model f u n c t i o n in terms of the new p a r a m e t e r s and t h e n to r e c a l c u l a t e the first and second derivatives matrices. A more efficient a p p r o a c h is to use the c h a i n rule of differentiation. Let Z. and Z.. be the first and second derivatives m a t r i c e s such that: (Z.)t~ =

~f(xt, v)/0v~,

(t =

1. . . . .

n; i =

1,...,p)

(23)

and:

(Z..)t~j = 0~f(xt, v)/c~vic~vj,

(j =

1 . . . . . p)

(24)

Using the first and second derivatives of the inverse t r a n s f o r m a t i o n u = T l(v), and the original first and second derivatives m a t r i c e s V. and V.., the r e q u i r e d m a t r i c e s are given by: (Z.)ti

(V.)ti(c~ui/~v~)

=

(25)

J

and: (Z..)~i~

(V')tk(~2Uk/~Vi~V~) + Z (V..)tkm(~Um/~Vi)(~Uh/~Vj)

= k

(k

=

1. . . . . p; m =

(26)

krn

1 . . . . . p)

Gauss-Newton optimization algorithm The G a u s s - N e w t o n m e t h o d was used to assess the ability of the p a r a m e t e r t r a n s f o r m a t i o n s r e p r e s e n t e d by eqns. (5) and (6) to speed the c o n v e r g e n c e of a u t o m a t i c p a r a m e t e r e s t i m a t i o n algorithms. This selection was based on the

62

B.C BATES

fact t h a t for a model f u n c t i o n t h a t is l i n e a r in its p a r a m e t e r s , the m e t h o d c o n v e r g e s to fi in a single i t e r a t i o n from a n y s t a r t i n g v e c t o r u0 (Ratkowsky, 1983, p. 41). The m e t h o d also c o n v e r g e s rapidly if the response surface is q u a d r a t i c in the n e i g h b o u r h o o d of the minimum, p a r t i c u l a r l y if this neighbourhood includes u 0 (Ross, 1970). The basic i t e r a t i o n in the G a u s s - N e w t o n m e t h o d is: U (j+l)

=

U (j) ~- J

(27)

1V.'

w h e r e u q) is the v e c t o r of p a r a m e t e r s at the j t h i t e r a t i o n and e = q - f(uO)). T e r m i n a t i o n of the i t e r a t i v e s e a r c h for the m i n i m u m of S(u) was made when the following c r i t e r i a were satisfied simultaneously: ]S(u (j)) - S ( u ~ + l ) ) f / S ( u (j)) < 5s

(28)

](u~ ~ 1 ) - u~))l/]u~)l < 5~,

(29)

(k = 1 , . . . , p )

[((n - p)e'V. J 1V.'e)/pS(u(J))]l/2

< 50

(30)

where 5s and (f~ are preselected t o l e r a n c e levels on the r e l a t i v e changes in the sum of squares and the p a r a m e t e r values, respectively, and were set at 10-4. E q u a t i o n (30) is a ' r e l a t i v e offset' c r i t e r i o n proposed by Bates and Watts (1981). It is a m e a s u r e of the o r t h o g o n a l i t y of the residual v e c t o r e to the solution locus and is t h e r e f o r e an absolute m e a s u r e of c o n v e r g e n c e to a local minimum. T h u s 50 is the t o l e r a n c e level used to decide w h e t h e r e is ~sufficiently o r t h o g o n a l ' and was set at 10 -3 as r e c o m m e n d e d by Bates and Watts. Bates and W a t t s (1981, 1988) also suggest t h a t a more n u m e r i c a l l y stable form for the r e l a t i v e offset c r i t e r i o n is:

(HQ;ell/x/p)/(lIQ:~ell/x/(n- p)) < 60

(31)

w h e r e Q1 and Q2 are the first p and the last n - p columns of the Q m a t r i x defined in eqn. (17).

Joint confidence regions for model parameters V a r i o u s m e t h o d s h a v e been proposed for c o m p u t i n g j o i n t confidence regions for the p a r a m e t e r s of n o n l i n e a r models. The usual a p p r o a c h is to assume t h a t the degree of n o n l i n e a r i t y is negligible, and to d e t e r m i n e the 100(1 ~)% l i n e a r t h e o r y confidence region defined by the hyperellipsoid: (u - fi)'P ~(u - ~) ~< pF(p, n - p , ~ )

(32)

w h e r e P = ~2j 1 is the a s y m p t o t i c c o n v a r i a n c e m a t r i x of ft. Such regions are i n e x p e n s i v e to c o n s t r u c t . A l t h o u g h eqn. (32) is exact for linear models, it is not always an a d e q u a t e t e c h n i q u e for a p p r o x i m a t i n g j o i n t confidence regions for the p a r a m e t e r s of a n o n l i n e a r model (Donaldson and Schnabel, 1987). E x a c t confidence regions for u can be o b t a i n e d by the lack-of-fit method:

NONLINEAR, DISCRETE FLOOD EVENT MODELS e'Re/e'(I

-

63

R)e ~< p F ( p , n - p , ~ ) / ( n - p )

(33)

where R = V. j-1V.', V. and e are evaluated at u, and I is the (n × n) identity matrix. Bates and Watts (1988, p. 223) present an alternative form: (llQ;ell2/p)/(llQ~ell2/(n

-

p))

<~ F ( p , n

p , ~)

(34)

The lack-of-fit method is much more computationally expensive than the linear theory method as f(u) and V. must be evaluated for a grid of u values to produce a contour representing the boundary of the region defined by eqn. (33) or eqn. (34). Moreover, the exact regions will contain every local optimum on the surface defined by eqn. (3). An alternative is to construct joint confidence regions based on the likelihood method. These regions are defined by: S(u) ~< S(fi)[1 + p F ( p , n - p , ~ ) / ( n - p)]

(35)

Construction of the boundary of likelihood method confidence regions requires the evaluation off(u) for a grid of u values. Thus they are less computationally expensive than exact regions, particularly in situations where the elements of V. must be evaluated by numerical differentiation. The regions are approximate for nonlinear models in the sense that the confidence coefficient 1 - a is not exact. However, experience indicates that the discrepancy between the actual and nominal confidence levels is usually small (Ratkowsky, 1983, p. 30; Donaldson and Schnabel, 1987). Comparison of linear theory and exact (or likelihood) joint confidence regions facilitates a visual assessment of the utility of parameter transformations. For a reparameterized model/data set combination with low I N and P E , mapping the linear theory region in the transformed parameters into the original parameter space should give a better approximation to the exact (or likelihood) region than the linear approximation based on the original parameters. These improved linear approximation confidence regions are much simpler and much less expensive to produce than likelihood or exact regions. CASE STUDY The runoff routing model (RORB) described by Mein et al. (1974) was chosen for the analysis. Briefly, the model consists of a rainfall excess model and a catchment storage model which routes computed rainfall excess hyetographs to produce a surface runoff hydrograph. Conceptual storage elements are used to represent distributed, temporary storage on the catchment surface. These elements are linked according to the geometry of the channel network. Each element is assumed to have a power function storage~lischarge relation: S

-

kkrQ m

(36)

where S is temporary storage, Q is the discharge from the element, k and m are model parameters applicable to the whole catchment (i.e. u = (k m)'), and kr is

64

B.C.BATES

\ N

ol

20 Indian

Ocean

4o

o3

o5 o6

0 I

50 L

100 km I

J Fig. 1. Location of study catchments in southwestern Australia (for key to numbers see Table 1).

t h e r e l a t i v e d e l a y t i m e a p p l i c a b l e to t h e e l e m e n t . F u r t h e r d e t a i l s o n t h e m o d e l c a n be o b t a i n e d f r o m L a u r e n s o n a n d M e i n (1983). R a i n f a l l - r u n o f f d a t a f r o m six c a t c h m e n t s i n s o u t h w e s t e r n A u s t r a l i a w e r e u s e d for t h e s t u d y as s h o w n i n Fig. 1. A k e y to t h e n u m b e r s s h o w n i n Fig. 1 is g i v e n i n T a b l e 1. D e t a i l s o n c a t c h m e n t c h a r a c t e r i s t i c s a n d r a i n f a l l - r u n o f f d a t a a r e a l s o l i s t e d i n T a b l e 1. R a i n f a l l e x c e s s h y e t o g r a p h s w e r e c a l c u l a t e d b y s u b t r a c t i n g a n i n i t i a l loss f r o m t h e t o t a l r a i n f a l l p r i o r to a p p l y i n g a c o n s t a n t TABLE 1 Characteristics of study catchments and data No.

1 2 3 4 5 6

Catchment

Wungong Bk Dirk Bk at Kentish Farm North Dandalup R. at Scarp Road Little Dandalup R. at Scarp Road Bancell Bk at Wagerup Harvey R. at Dingo Road

Number of storms

Peak surface runoff (m3s 1) Minimum

Maximum

146 35.2

8 13

8.98 0.37

25.5 5.73

153

13

0.34

7,91

39.6

16

0.63

4.69

13.3 148

6 24

1.05 0.07

2.73 18.7

Area (km2)

NONLINEAR. D1SCRETE FLOOD EVENT MODELS

65

proportional loss. During these calculations the volume of rainfall excess was constrained to equal the volume of surface runoff. A more complete description of the study catchments and data is given by Bates and Townley (1988). Although eqn. (36) is imbedded in the calculation of the response variable (peak discharge), it does not necessarily follow that f(xt, u) = kh(xt) m where h(xt) denotes some function of the rainfall input xt. Thus calculation of the V. and V.. matrices involves the use of numerical methods. The finite difference schemes and the step size selection procedure used to compute these matrices are discussed in Bates (1988). Investigations showed that the first and second derivative estimates were stable for step sizes within the interval 0.1A < A < 10A where A denotes the selected step size.

Nonlinearity measures (original parameters) During the course of the investigation, an error was detected in the FORTRAN IV subroutines listed by Ratkowsky (1983, appendix) for the calculation of the Bates and Watts curvature measures. Details of the error are given elsewhere (Bates, 1990). Table 2 reports the revised I N and PE values for the original model parameterization and the six data sets. (The I N values for the Little Dandalup River and Wungong Brook catchments, and the P E value for the Harvey River catchment, are quite different from the values reported earlier by Bates (1988, table 3).) Observe that the intrinsic nonlinearity falls within or is very close to acceptable limits for every model/data set combination. In contrast, the parameter-effects nonlinearity exceeds the guide value 1/2x/F(2 , n - 2, 0.05) in five out of the six cases examined. Three of these P E values also exceed the guide value 1/2x/F(2, n - 2, 0.5). Application of the 1% rule of thumb to the percentage bias values listed in Table 2 indicates that the bias in/~ is relatively small in every case, and that the bias in rh is relatively large in five cases. This suggests that m is largely responsible for the nonlinear behaviour of these model/data set combinations. TABLE 2 Values of nonlinearity Catchment

Bancell Bk Dirk Bk Harvey R Little Dandalup R North Dandalup R Wungong Bk

measures for original model parameterization Curvatures

Guide value a

Percentage

IN

PE

I

II

k

m

gk

g~

0.02 0.26 0.15 0.13 0.19 0.22

0.41 0.40 0.10 0.87 2.0 5.7

0.19 0.25 0.27 0.26 0.25 0.22

0.55 0.58 0.59 0.59 0.58 0.57

- 0.66 - 0.99 0.25 0.32 - 0.64 0.14

1,6 2,4 - 0,32 1,7 4,3 7.2

- 0.45 - 0.35 0.11 - 0.06 - 0.08 0.77

0.80 0.79 - 0.36 0.65 1.2 0.85

a G u i d e v a l u e : I = 1 / 2 x / F ( 2 , n - 2, 0.05); II = 1 / 2 x / F ( 2 , n -

2, 0.5).

bias

g vector

66

B.C.BATES

Consider the values of the c o m p o n e n t s of the g v e c t o r r e p o r t e d in Table 2. E x c e p t for k in the H a r v e y , Little and N o r t h D a n d a l u p cases, k and m show evidence of m a r k e d skewness in t h e i r least squares estimators. It is also a p p a r e n t t h a t the estimates of m are positively skewed in most cases, and t h a t the skew for this p a r a m e t e r is quite s t r o n g for the N o r t h D a n d a l u p c a t c h m e n t . Parameter transformations

Use of the 'expected-value' t r a n s f o r m a t i o n requires the r e p l a c e m e n t of the model p a r a m e t e r s by 'expected-value' parameters. The RORB model does not provide an explicit r e l a t i o n s h i p between rainfall excess and peak discharge. N e v e r t h e l e s s an 'expected-value' t r a n s f o r m a t i o n of the model p a r a m e t e r s k and m in eqn. (36) can be made. Suppose a realistic pair of values Q* and Q* for Q are chosen such t h a t Q* < Q~, and t h a t kr is set equal to its m e d i a n v a l u e for the c a t c h m e n t kr. L e t s = (sl s2)' d e n o t e the 'expected-value' p a r a m e t e r s vector. P a r a m e t e r s s~ and s2 are defined by: sl

=

k £ Q*"

(37)

k k r Q *'~

(38)

and: s2 =

Solving for k and m in terms of sl and s2 yields: k

=

(s]/]~r)(S2/Sl

) - log

V~/|og(Q~/V~)

(39)

and: m

=

(40)

log(s2/s~)/log (Q~/Q*)

which define the inverse t r a n s f o r m a t i o n u = T-1(s). S u b s t i t u t i o n of eqns. (39) and (40) into eqn. (36) gives: S

:

(41)

81 (s2/s 1 )log(Q/Q~)/log(V~/Q~)

F r o m eqn. (41) it follows t h a t S -- sl when Q = Q*, and S = s2 when Q = Q~. Thus eqn. (41) defines the t r a n s f o r m a t i o n s = T(u) and the t r a n s f o r m a t i o n c o n s t a n t s Q* and Q~ can be i n t e r p r e t e d as allowable values of Q. Let v = (vlv2)' d e n o t e the power t r a n s f o r m e d p a r a m e t e r vector. F r o m eqn. (6) it follows t h a t p a r a m e t e r s vl and v2 are defined by: V1

k ~'~

if71 ~

0

In (k)

ifT~ =

0

m :2

if72 ~

0

In (m)

if 72 = 0

(42)

and V2

(43)

It should be n o t e d t h a t it is difficult to provide a physical i n t e r p r e t a t i o n of 71 and 72.

67

NONLINEAR. DISCRETE FLOOD EVENT MODELS

I n i t i a l l y a t t e m p t s w e r e m a d e to find the v a l u e s of t h e t r a n s f o r m a t i o n c o n s t a n t s Q*, Q*, 71, a n d ?2 s u c h t h a t P E was minimized, a n d t h o s e v a l u e s w h i c h g a v e g -- 0 w h e r e 0 is t h e zero vector. No a t t e m p t w a s m a d e to minimize p e r c e n t a g e bias b e c a u s e of its l o c a t i o n d e p e n d e n c e . T h e m i n i m i z a t i o n of P E p r o v e d to be a t i m e c o n s u m i n g process. C o n t o u r s of P E defined by v a r y i n g the Q* (i = 1, 2) for example, w e r e found to be m a r k e d l y n o n - q u a d r a t i c , a n d r e g i o n s w i t h i n the Q* s p a c e were e n c o u n t e r e d w i t h i n w h i c h the i t e r a t i v e p r o c e d u r e devised by B a t e s and W a t t s (1980) for e s t i m a t i n g the m a x i m u m i n t r i n s i c a n d p a r a m e t e r - e f f e c t s c u r v a t u r e s failed to c o n v e r g e (see Fig. 2). T h e s e r e g i o n s w e r e a s s o c i a t e d w i t h r e l a t i v e l y small v a l u e s of P E or Q* v a l u e s w h e r e Q~ ~ Q~. T h u s P E was minimized to the e x t e n t t h a t v a l u e s of t h e t r a n s f o r m a tion c o n s t a n t s w e r e found s u c h t h a t P E was less t h a n or close to the guide v a l u e 1/2x/F(2, n - 2, 0.05) unless c o n v e r g e n c e p r o b l e m s w e r e e n c o u n t e r e d . T h i s m i n i m i z a t i o n was u n d e r t a k e n only w h e n t h e g = 0 c r i t e r i o n failed to p r o d u c e a n a c c e p t a b l e v a l u e of PE. F o r the ' e x p e c t e d - v a l u e ' t r a n s f o r m a t i o n , it was possible to find two points Q* and Q* w h i c h m e t the g = 0 criterion. F i g u r e 3 shows a plot of t h e g vs. Q* c u r v e for the D i r k B r o o k c a t c h m e n t . O b s e r v e t h a t the zeros of g a r e well defined, a n d t h a t t h e s e Q* v a l u e s are r e a l i s t i c in t h e sense t h a t t h e y fall w i t h i n the r a n g e of ' o b s e r v e d ' p e a k s u r f a c e runoff. F o r t h e B a n c e l l a n d W u n g o n g B r o o k c a t c h m e n t s , one of t h e i r zeros of g did not fall w i t h i n t h e i r r a n g e s of ' o b s e r v e d ' p e a k s u r f a c e runoff. H o w e v e r , t h e s e zeros could still be r e g a r d e d as r e a l i s t i c since t h e y w e r e positive a n d were r e a s o n a b l y close to the b o u n d s of lO

~3

E

o o

5

lO

Q,~ ( m 3 s -1 ) Fig. 2. Contours of PE, Dirk Brook catchment. (Shaded areas and the dashed straight line depict the Q~ and Q~ values at which the iterative algorithm of Bates and Watts (1980) failed to converge.) The PE = 0.2 contour is very close to the shaded area, and has been omitted for the sake of clarity.

68

B.c. BATES

Q* (rn3 s-l)

c~.

x

E -I Fig. 3. Standardized 'asymptotic' third m o m e n t of H o u g a a r d (1985) as a function of Q*, Dirk Brook catchment.

their respective observed ranges. For the parameter power transformation, values of 71 and 72 which made the components of g zero were easily found by trial and error. Overall, the determination of transformation constants which satisfied the g -- 0 criterion took a small fraction of the time required to find the constant values which minimized PE.

Nonlinearity measures (transformed parameters) Table 3 summarizes the values of the nonlinearity measures for the alternative model reparameterizations and the procedures for optimization of the transformation constants. Percentage biases are not listed since they were below 1% for every reparameterized model/data set combination. Comparison of Tables 2 and 3 reveals t ha t application of both classes of transformations can lead to marked reductions in the PE values obtained for the original parameterization. Seven out of the twelve reparameterized model/data set combinations in Table 3 have acceptable values of PE at ~ = 0.05 when the g = 0 criterion is used. F u r t h e r reductions in PE for the remaining five combinations were usually achieved at the expense of increasing the skewness of the least squares estimators. The only exception is the power transformation/North Dandalup River combination for which a further reduction in PE could not be obtained. However, the attempts for further reductions in PE produced PE values that were acceptable at ~ --- 0.05 in three out of these five cases. A comparison was made of the processor time (CPU) required to calculate the n on lin ear ity measures for the reparameterizations by using the chain rule

69

NONLINEAR, DISCRETE FLOOD EVENT MODELS

i ¢~c5 I

c5c5 I

~k

¢5¢5

c5

I

O

O

¢5 ¢qI

¢4 o¢5

c5 I

c'q II

¢5d

¢~

O

H I

--> 5'q H O

hi

70

B.C. B A T E S

~2 I

I

I

I

o

o%

~g

0

-~ ?~ o~ m, •

,

o

0

,~.N 0

~

0

2

2

o

r..)

~

~

~

~

~

~

~

Z

0

NONLINEAR,

DISCRETE FLOOD EVENT MODELS

Z ¢q 5"q ¢~

(~

I v

qJ

o

~

O

O

.

.

.

.

~9

O

=~o

~-~

o~ oo

m~

0

0 .l~, .~

o~

71

72

B.C.BATES

with t hat required by expressing the model function in terms of the new parameters and calculating the Z. and Z.. matrices. Results showed t hat use of the chain rule gave a one to five hundred-fold reduction in processor time.

Convergence of Gauss-Newton algorithm All of the starting vectors u0 used to investigate the effect of parameter transformations on the rapidity of convergence of the Gauss-Newton algorithm were taken from the boundaries of the linear theory 50, 95 and 99.5% joint confidence regions defined by eqn. (32). Figure 4 illustrates the configuration of the twelve starting vectors when the end-points of the principal axes of the ellipses lie within the feasible region of parameter space. Configurations of this type were used in four out of the six data sets. For the two remaining sets, the position of the starting vectors was compromised according to the parts of the confidence regions having physical significance. Three sets of optimization runs were conducted according to each type of model parameterization. Table 4 lists the original parameter values, the values of the transformation constants which satisfied the g = 0 criterion, the total number of iterations for the four starting vectors on the boundary of each of the three linear theory joint confidence regions, and the processor time required to achieve convergence. Clearly, the use of parameter transformations has led to only a marginal increase in the speed of convergence. In fact the use of an 'expected-value' parameterization has led to a slight reduction in speed in two out of the six cases examined. Also, the values of y, and Y2differ greatly from data set to data set, and Q* and Q~' do not fall at roughly the same locations with respect to the

05

6

q,C -

0.05

I

k

~

*~"

0.50

4

a

major

0

I 1

ax,s

I 2

m

Fig. 4. 100(1 - a)% linear theory joint confidenceregions and configuration of initial parameter estimates for the Gauss-Newton algorithm, Bancell Brook catchment.

73

NONLINEAR, DISCRETE FLOOD EVENT MODELS

b o u n d s of the r a n g e s of ' o b s e r v e d ' p e a k s u r f a c e r u n o f f listed in T a b l e 1. T h e r e f o r e t h e r e is no c o n s i s t e n c y in the o p t i m u m v a l u e s of the t r a n s f o r m a t i o n c o n s t a n t s f r o m one d a t a set to a n o t h e r . I t is s o m e t i m e s s t a t e d t h a t the use of p a r a m e t e r i z a t i o n s w h i c h g r e a t l y r e d u c e the m a g n i t u d e of p a r a m e t e r c o r r e l a t i o n s c a n i m p r o v e t h e r a p i d i t y of convergence. T a b l e 5 lists t h e p a r a m e t e r c o r r e l a t i o n s for the original model para m e t e r i z a t i o n a n d the a l t e r n a t i v e r e p a r a m e t e r i z a t i o n s . T h e c o r r e l a t i o n s for the ' e x p e c t e d - v a l u e ' p a r a m e t e r s are m u c h l o w e r t h a n t h o s e for the o r i g i n a l p a r a m e t e r s in five out of the six cases examined. (The p o w e r t r a n s f o r m a t i o n does not c h a n g e t h e a b s o l u t e v a l u e s of the o r i g i n a l p a r a m e t e r c o r r e l a t i o n s since 7i is i n d e p e n d e n t of ~;~for i ~ j.) Clearly, t h e s e r e d u c t i o n s h a v e n o t led to a m a r k e d i m p r o v e m e n t in the speed of c o n v e r g e n c e of the G a u s s - N e w t o n algorithm. H o w e v e r , T a b l e 4 i n d i c a t e s t h a t the G a u s s - N e w t o n a l g o r i t h m failed to c o n v e r g e for the H a r v e y R i v e r d a t a a n d a n u m b e r of trial s t a r t i n g v e c t o r s w h e n the original model p a r a m e t e r i z a t i o n was used. This d i v e r g e n t b e h a v i o u r was not e x p e r i e n c e d for e i t h e r of the r e p a r a m e t e r i z e d models. Overall, the r e s u l t s a p p e a r to be in a g r e e m e n t w i t h t h o s e of R a t k o w s k y a n d Reedy (1986). T h e i r s i m u l a t i o n e x p e r i m e n t s i n d i c a t e d t h a t while r e p a r a m e t e r i zation h a s little effect on the n u m b e r of i t e r a t i o n s r e q u i r e d to o b t a i n convergence, it m a y d e c r e a s e the p r o b a b i l i t y of divergence.

Joint confidence regions for model parameters B o u n d a r i e s of the 95% l i n e a r theory, likelihood a n d e x a c t j o i n t confidence r e g i o n s for the o r i g i n a l model p a r a m e t e r s were o b t a i n e d for the six d a t a sets. E x p e r i e n c e s h o w e d t h a t a c c u r a t e d e t e r m i n a t i o n of t h e e x a c t confidence regions r e q u i r e d c a r e in the selection of the s t o p p i n g c r i t e r i o n for t h e nonl i n e a r flow r o u t i n g c a l c u l a t i o n s in the R O R B model, p a r t i c u l a r l y for t h a t p a r t of the p a r a m e t e r s p a c e defined by m < 0.5. In e v e r y case t h e use of a c o a r s e c r i t e r i o n affected the s m o o t h n e s s of the model r e s p o n s e f u n c t i o n a n d this led to n u m e r i c a l i n s t a b i l i t i e s in the first d e r i v a t i v e c a l c u l a t i o n s . S u p e r p o s i t i o n of the t h r e e r e g i o n s r e v e a l e d c o n s i d e r a b l e d i s c r e p a n c i e s TABLE 5 Parameter correlations Catchment

k and m

s1 and s2

v1 and v2

Bancell Bk Dirk Bk Harvey R Little Dandalup R North Dandalup R Wungong Bk

-

0.301 0.566 0.822 0.442 0.441 - 0.681

0.876 0.865 0.861 0.952 0.967 0.990

0.876 0.865 0.861 0.952 0.967 0.990

74

B.C. BATES

between the linear approximation and the exact regions for all but the Bancell Brook and Harvey River catchments (e.g. Figs. 5 and 6). Good agreement was found between the likelihood and exact regions in every case. The boundaries of the 95% linear theory joint confidence regions in the 'expected-value' parameters s were mapped through the inverse transformation defined by eqns. (39) and (40) to give approximate confidence regions in u. Superposition of these regions and their corresponding exact regions showed that the 'expected-value' transformation could give confidence regions that were superior to those obtained from the linear approximation in the u parameters (e.g. Fig. 5). However, the back-transformation was not completely successful in several cases. Consider Fig. 6a where the 95% linear theory confidence region in the s parameters crosses the s2 axis. It follows from eqns. (39) and (40) that k --* 0 and m --* ~ as sl ~ 0, and that k and m are undefined for sl < 0. (a)

sl

1

02.5

r

I

3.5

4.5

S2

12

10

k

8

(b) ....-. :

~ Exact .... L i n e a r a D l ~ r o x i m a t i o n

ir~ u

'.

--- Linear approximation

in s

':

6

4

2 0.2

I

0.6

1.0

1.4

1.8

2.2

131

Fig. 5. 95% joint confidence regions for parameter estimates, Dirk Brook catchment: (a) linear approximation region in the s parameters: (b) exact region, and linear approximation regions in the u and back-transformed s parameters.

75

NONLINEAR, DISCRETE FLOOD EVENT MODELS 8

s1

4

O

J 0

20

40

80

s2 40

(b)

~\\ ~\ ~

30

.... Linear aparoxlmatlon - - L i n e a r approximation

in. in s

',. ".. ~ "... "~.

2O

0

-10 0.2

I

I

I

0.6

1.0

1.4

"".)1

212.

1.8

m

joint confidence regions for parameter estimates, Wungong Brook catchment: (a) linear approximation region in the s parameters; (b) exact region, and linear approximation regions in the u and back-transformed s parameters. Fig. 6.95%

Thus mapping the linear theory confidence region in the s parameters into the original parameter space gives a confidence region that approaches the m axis asymptotically (see Fig. 6b). A similar result was obtained h)r the North Dandalup River catchment. For the Bancell Brook catchment, the 95% linear theory confidence region in the s parameters almost touches the s2 axis, and hence the back-transformed region has an elongated tail which extends well beyond the boundary of the exact confidence region for the u parameters. Although this behaviour is undesirable when inferences are to be made in

76

B.C. BATES

terms of the u parameters, it does not affect the quality of inferences that are made directly in terms of the s parameters. The boundaries of the 95% linear theory joint confidence regions for the power transformed parameters v were also mapped into the original parameter space. Superposition of the back-transformed regions and their corresponding exact regions showed that the parameter power transformation could give confidence regions that were superior to those obtained from the linear approximation in the u parameters (e.g. Figs. 7 and 8). Although the 95% linear theory regions in the v parameters crossed the s2 axis for the North Dandalup and Wungong Brook cases, the mapping of these regions through the inverse transformation led to truncated banana-shaped regions that were very good approximations to the exact region. The success of these approximate regions is somewhat surprising since the P E values for the North Dandalup River and Wungong Brook data (0.36 and 0.54, respectively) are well in excess of the guide values for a = 0.05 (0.25 and 0.22, respectively). This suggests that the guide value 1/2x//F(2, n - p , 0.05) may be too conservative for these particular cases. Table 4 reveals that the magnitudes of some of the 71 and ~2 values for the Bancell Brook, Harvey River and Little Dandalup River catchments are quite high. A high absolute value of 7i indicates that the skewness of ui is somewhat insensitive to the power transformation. The agreement between the backtransformed confidence regions and their corresponding exact regions may be described as excellent for the Little Dandalup catchment and reasonable for the Bancell Brook and Harvey River catchments. For the Harvey River catchment, however, the linear theory confidence region in the u parameters is in better agreement with the exact region. This result might reflect the marked increase

12 Exact .'"~'~

...... Linear approximation in u nin

10

2 0.2

0.6

1.0

1.4

1.8

v

2.2

m

Fig. 7. 95% j o i n t confidence r e g i o n s for p a r a m e t e r e s t i m a t e s , Dirk B r o o k c a t c h m e n t . L i n e a r a p p r o x i m a t i o n r e g i o n s a r e for t h e u a n d b a c k - t r a n s f o r m e d v p a r a m e t e r s .

N O N L I N E A R , D I S C R E T E FLOOD E V E N T M O D E L S 40

77

Exact ...... Linear approximation

in u

-- Linear approximation

in v

:ii:::iiiiiiiii ;:: 30

20

10 "...~

'>.."2..--_-_-_ •.. '.. -10 0.2

0.6

I.

I

1 0

1.4

"'-:I 1.8

J 2.2

Fig. 8. 95% joint confidence regions for parameter estimates, Wungong Brook catchment. Linear approximation regions are for the u and back-transformed v parameters.

in the maximum parameter-effects curvature incurred by making the power transformation (compare Table 2 and 3). A comparison was also made of the CPU required to compute the exact regions with that required for the back-transformed regions. Results showed that calculation of the exact regions required ten to one hundred times as much processor time. CONCLUSIONS AND RECOMMENDATIONS This paper has investigated the application of the parameter transformations proposed by Ross (1970) and Tsai (1988) to nonlinear, discrete flood event models. The foi|owing conclusions may be presented on the basis of the results reported herein. (1) Marked improvements in the speed of convergence of the Gauss-Newton algorithm cannot always be guaranteed by the use of parameter transformations which make the response surface more nearly quadratic. However, there is some evidence that the use of parameter transformations can reduce the probability of divergence during the estimation process. (2) The standardized asymptotic third moment of Hougaard (1985) is a useful measure of statistical nonlinearity in individual least squares estimators. (3) The 'expected-value' transformation of Ross (1970) and the parameter power transformation of Tsai (1988) can be used advantageously to reduce the asymmetry of least squares estimators to acceptable levels and to improve parameter inferences a posteriori.

78

B.C. BATES

(4) O p t i m i z a t i o n of t r a n s f o r m a t i o n c o n s t a n t s by r e d u c i n g the a s y m p t o t i c t h i r d m o m e n t to zero c a n be a n effective m e a n s of r e d u c i n g the a s y m p t o t i c bias of l e a s t s q u a r e s e s t i m a t o r s and m a x i m u m p a r a m e t e r - e f f e c t s c u r v a t u r e . (5) T h e d e t e r m i n a t i o n of e x a c t confidence r e g i o n s for the p a r a m e t e r s of n o n l i n e a r models is quite e x p e n s i v e w h e n first d e r i v a t i v e s of the model r e s p o n s e f u n c t i o n with r e s p e c t to the p a r a m e t e r s h a v e to be o b t a i n e d by n u m e r i c a l methods. In these s i t u a t i o n s , the m a p p i n g of l i n e a r t h e o r y confidence r e g i o n s for the t r a n s f o r m e d p a r a m e t e r s into the o r i g i n a l p a r a m e t e r s p a c e c a n be an efficient m e a n s of o b t a i n i n g reliable j o i n t confidence regions for the model p a r a m e t e r s . (6) T h e ~expected-value' t r a n s f o r m a t i o n p r o p o s e d by Ross (1970) c a n be s u s c e p t i b l e to u n d e s i r a b l e b e h a v i o u r d u r i n g the b a c k - t r a n s f o r m a t i o n process. This is not a c a u s e for c o n c e r n if i n f e r e n c e s a r e m a d e in t e r m s of the ~expectedv a l u e ' r a t h e r t h a n the o r i g i n a l p a r a m e t e r s . T h e findings of this s t u d y are b a s e d on r e s u l t s o b t a i n e d from six c a t c h m e n t s w i t h i n the s a m e g e o g r a p h i c region. F u t u r e r e s e a r c h should c o n c e n t r a t e on the a p p l i c a t i o n of the n o n l i n e a r i t y m e a s u r e s and classes of p a r a m e t e r t r a n s f o r m a tions described in this p a p e r to o t h e r case studies i n v o l v i n g different p a r a m e t e r e s t i m a t i o n a l g o r i t h m s , flood events models and h y d r o l o g i c regimes. ACKNOWLEDGEMENTS I am i n d e b t e d to D a v i d A. R a t k o w s k y , P r o g r a m in Statistics, W a s h i n g t o n S t a t e U n i v e r s i t y , for access to his u n p u b l i s h e d m a n u s c r i p t s on the r e p a r a m e t e r i z a t i o n of n o n l i n e a r r e g r e s s i o n models. T h a n k s a r e also due to G e o r g e K u c z e r a , D e p a r t m e n t of Civil E n g i n e e r i n g and S u r v e y i n g , The U n i v e r s i t y of N e w c a s t l e (Australia), a n d D a v i d R a t k o w s k y for r e v i e w i n g an e a r l i e r v e r s i o n of the paper. The r a i n f a l l - r u n o f f d a t a described in this p a p e r were t a k e n from the a r c h i v e s of the W a t e r A u t h o r i t y of W e s t e r n A u s t r a l i a . Lastly, I am p a r t i c u l a r l y g r a t e f u l to m y c o l l e a g u e Neil R. S u m n e r for his c o m p u t i n g assistance. REFERENCES Bates, B.C., 1988. Nonlinear, discrete flood event models, 2. Assessment of statistical nonlinearity. J. Hydrol., 99: 77-89. Bates, B.C., 1990. Nonlinear, discrete flood event models, 2. Assessment of statistical nonlinearity -- Correction. J. Hydrol. (submitted). Bates, B.C. and Townley, L.R., 1988. Nonlinear, discrete flood event models, 1. Bayesian estimation of parameters. J. Hydrol., 99:61 76. Bates, D.M. and Watts, D.G., 1980. Relative curvature measures of nonlinearity (with Discussion). J.R. Stat. Soc., B, 42(1): 1 25. Bates, D.M. and Watts, D.G., 1981. A relative offset orthogonality convergence criterion for nonlinear least squares. Technometrics, 23(2): 17~183. Bates, D.M. and Watts, D.G., 1988. Nonlinear regression analysis and its applications. Wiley, New York, N.Y., 365pp. Box, M.J., 1971. Bias in nonlinear estimation (with Discussion). J.R. Stat. Soc., B, 33(2): 171 201.

NONLINEAR, DISCRETE FLOODEVENT MODELS

79

Carrera, J. and Neuman, S., 1986. Estimation of aquifer parameters under transient and steady state conditions, Parts 1, 2, and 3. Water Resour. Res., 22(2): 199-242. Cook, R.D. and Witmer, J.A., 1985. A note on parameter-effects curvature. J. Am. Stat. Assoc., 80(392): 872-878. Cook, R.D., Tsai, C.-L. and Wei, B.C., 1986. Bias in nonlinear regression. Biometrika, 73(3): 6154;23. Donaldson, J.R. and Schnabel, R.B., 1987. Computational experience with confidence regions and confidence intervals for nonlinear least squares. Technometrics, 29(1): 67-82 Gupta, V.K. and Sorooshian, S., 1983. Uniqueness and observability of conceptual rainfall runoff model parameters: The percolation process examined. Water Resour. Res., 19(1): 269 276. Hougaard, P., 1985. The appropriateness of the asymptotic distribution in a nonlinear regression model in relation to curvature. J.R. Stat. Soc., B, 47(1): 103-114. [bbitt, R.P. and O'Donnell, T., 1974. Designing conceptual catchment models for automatic fitting methods. In: Mathematical Models in Hydrology Symposium, IAHS-AISH Publ. 2:461 475. Johnston, P.R. and Pilgrim, D.H., 1976. Parameter optimization for watershed models. Water Resour. Res., 12(3): 477~486. Kuczera, G., 1983. Improved parameter inference in catchment models, 1. Evaluating parameter uncertainty. Water Resour. Res., 19(5): 1151-1162. Laurenson, E.M. and Mein, R.G., 1983. RORB - - version 3 runoff routing program user manual. Dep. Civ. Eng., Monash University, Clayton, Vic., 2nd edn. Mein, R.G., Laurenson, E.M. and McMahon, T.A., 1974. Simple nonlinear model for flood estimation. J. Hydraul. Div., Proc. Am. Soc. Civ. Eng., 100(HYll): 1507-1518. Morton, R., 1987. Asymmetry of estimators in nonlinear regression. Biometrika, 74(4): 679~85. Ratkowsky, D.A., 1983. Nonlinear Regression Modeling. Dekker, New York, N.Y., 276 pp. Ratkowsky, D.A., 1986. Statistical properties of alternative parameterizations of the von Bertalanffy growth curve. Can. J. Fish. Aquat. Sci., 43:742 747. Ratkowsky, D.A., 1990. Handbook of Nonlinear Regression Models. Dekker, New York, NY, in press. Ratkowsky, D.A. and Reedy, T.J., 1986. Choosing near-linear parameters in the four-parameter logistic model for radioligand and related assays. Biometrics, 42(3): 575-582. Restrepo-Posada, P.J. and Bras, R.L., 1982. Automatic parameter estimation of a large conceptual rainfall runoff model: A maximum likelihood approach. Ralph M. Parsons Lab., Hydrol. Water Resour. Systems, Dep. Civ. Eng., MIT, Rep. No. 267. Ross, G.J.S., 1970. The efficient use of function minimization in non-linear maximum-likelihood estimation. Appl. Stat. 19(3): 205 221. Sorooshian, S. and Dracup, J.A., 1980. Stochastic parameter estimation procedures for hydrologic rainfall runoff models: Correlated and heteroscedastic error cases. Water Resour Res., 16(2): 430~442. Sorooshian, S. and Gupta, V.K., 1985. The analysis of structural identifiability: Theory and application to conceptual rainfall-runoff models. Water Resour. Res., 21(4): 487495. Tsai, C.-L., 1988. Power transformations and reparameterizations in nonlinear regression models. Technometrics, 30(4): 441448. Williams, B.J. and Yeh, W.W.-G., 1983. Parameter estimation in rainfall runoff models. J. Hydrol., 63:373 393.