Next Article in Journal
Are Presence/Absence Microbial Tests Appropriate for Monitoring Large Urban Water Supplies in Sub-Saharan Africa?
Next Article in Special Issue
Statistical Analysis of Extreme Events in Precipitation, Stream Discharge, and Groundwater Head Fluctuation: Distribution, Memory, and Correlation
Previous Article in Journal
Integration of Isotopic (2H and 18O) and Geophysical Applications to Define a Groundwater Conceptual Model in Semiarid Regions
Previous Article in Special Issue
Combing Random Forest and Least Square Support Vector Regression for Improving Extreme Rainfall Downscaling
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Third-Order Polynomial Normal Transform Applied to Multivariate Hydrologic Extremes

1
Disaster Prevention and Water Environment Research Center, National Chiao Tung University, Hsinchu 300, Taiwan
2
Department of Civil, Environment, and Architecture, Korea University, Seoul 02841, Korea
*
Author to whom correspondence should be addressed.
Water 2019, 11(3), 490; https://doi.org/10.3390/w11030490
Submission received: 26 December 2018 / Revised: 22 February 2019 / Accepted: 4 March 2019 / Published: 8 March 2019

Abstract

:
Hydro-infrastructural systems (e.g., flood control dams, stormwater detention basins, and seawalls) are designed to protect the public against the adverse impacts of various hydrologic extremes (e.g., floods, droughts, and storm surges). In their design and safety evaluation, the characteristics of concerned hydrologic extremes affecting the hydrosystem performance often are described by several interrelated random variables—not just one—that need to be considered simultaneously. These multiple random variables, in practical problems, have a mixture of non-normal distributions of which the joint distribution function is difficult to establish. To tackle problems involving multivariate non-normal variables, one frequently adopted approach is to transform non-normal variables from their original domain to multivariate normal space under which a large wealth of established theories can be utilized. This study presents a framework for practical normal transform based on the third-order polynomial in the context of a multivariate setting. Especially, the study focuses on multivariate third-order polynomial normal transform (TPNT) with explicit consideration of sampling errors in sample L-moments and correlation coefficients. For illustration, the modeling framework is applied to establish an at-site rainfall intensity–duration-frequency (IDF) relationship. Annual maximum rainfall data analyzed contain seven durations (1–72 h) with 27 years of useable records. Numerical application shows that the proposed modeling framework can produce reasonable rainfall IDF relationships by simultaneously treating several correlated rainfall data series and is a viable tool in dealing with multivariate data with a mixture of non-normal distributions.

1. Introduction

In hydrosystem design, performance evaluation, and simulation, the problems often involve multiple random variables that are correlated with a mixture of non-normal marginal distributions. Under this condition, it is generally difficult, if not impossible, to establish an analytical joint probability distribution for these variables. In comparison with univariate distributions, there are relatively few analytical multivariate distribution functions under special combinations of parametric marginal distributions, and most of them are of the same type, which can be found in [1,2]. Examples of using analytical multivariate distributions in hydrology are bivariate Gamma distribution [3] and bivariate generalized extreme distribution [4]. Their use is somewhat limited to many practical problems because of different marginal distributions.
Due to the difficulty in establishing a truly multivariate joint distribution model for problems involving mixtures of several correlated, non-normal variables, approximated approaches, such as copula or normal transform, are often used by preserving marginal distributions or moments, including the correlation features among the variables. However, one should realize that, unlike using a true multivariate joint distribution function, preservation of the marginal distributions and dependence structure represents the retention of partial information of the concerned multivariate random variables in the analysis [5].
The concept of copula is one type of approximated multivariate approaches that has recently received tremendous attention and applications by researchers in various disciplines, including in hydrology [6]. Some examples of applying copula in multivariate hydrologic modeling can be found in analyzing floods [7,8], droughts [9,10,11,12], dam safety [13], and extreme rainfalls [14,15]. Most of the copula-based applications deal with bivariate problems and some trivariate problems under some restrictive conditions on correlation structures [8]. Applications of copula to higher dimension multivariate problems are rare primarily because there are only a few copula families that are rather restrictive in describing the dependence structure. Recently, the introduction of vine copulas has shown the advantage of overcoming the limitation of currently used copulas in multivariate analysis [16,17,18,19]. A copula-based approach is parametric by nature in that analytical marginal distribution models for the involved variables are specified.
Alternatively, another viable scheme in treating multivariate problems involving correlated non-normal random variables is to apply a NORTA (normal-to-anything) algorithm [20]. By a NORTA algorithm, normal transformation of an individual non-normal variable is made by preserving its marginal probability content in the normal variable domain as Φ ( z ) = F x ( x ) with Φ ( · )   and   F x ( · ) , respectively, being the cumulative distribution functions (CDFs) of the standard normal variable Z and the original variable X. In addition, a relationship must be established to allow the determination of an equivalent correlation coefficient, ρ z j , z k , of a pair of normal transformed variables, Z j   and   Z k , from the correlation coefficient, ρ x j , x k , of the corresponding random variables, X j   and   X k , in the original space. Once the correlation matrix of standard normal variables Z’s, { ρ z j , z k } , is obtained from that of the non-normal variables X’s, { ρ x j , x k } , appropriate orthogonal transformation can be implemented to transform the original correlated variables into uncorrelated standard normal space for analysis.
The determination of ρ z j ,   z k from ρ x j ,   x k is made through the Nataf transform [21], which requires solving an implicit non-linear equation in the form of a double integration involving marginal distributions of a pair of random variables, X j   and   X k , under consideration:
ρ x j , x k = ( x j μ j σ j )   ( x k μ k σ k )   ϕ j k ( z j , z k | ρ z j , z k )   d z j   d z k
where x j = F j 1 [ Φ ( z j ) ] , and ϕ j k ( · ) = bivariate standard normal joint probability density function (PDF). Lebrun and Dutfoy [22] provide an insightful analysis of Nataf transform and uncover that it is a special modeling of dependence structure using Gaussian copula. To facilitate practical engineering applications, a set of empirical equations for 10 commonly used distribution functions has been established to relate ρ z j , z k to ρ x j , x k and their distribution properties [23]. Such empirical relations were applied to reliability analysis of engineering systems [5,24]. Later, computationally more efficient methods based on root finding and linear search [25], the false position method [26], and the artificial neural network method [27] were proposed to solve Equation (1) for ρ z j , z k from the known ρ x j , x k and marginal PDFs of X j   and   X k .
The above mentioned schemes (i.e., copula, NORTA, and Nataf transform) all require the stipulation of marginal PDFs. The stipulation of a distribution function implies knowing the complete statistical information of the random variable, including its moments of all orders. This ideal situation is attainable only when one has a large amount of data, which generally is not the case in practice. Therefore, to relax the information requirement without having to specify the distribution functions, third-order polynomial normal transform (TPNT) can be used. By TPNT, each individual non-normal random variable is related to a 3rd-order polynomial function of the corresponding standard normal variable [28]. The polynomial coefficients are determined by matching the statistical moments or quantiles of the individual random variables. The multivariate version of TPNT was first proposed by Vale and Maurelli [29] to simultaneously consider statistical moments and correlation coefficients. A multivariate TPNT procedure has been applied to different fields including, but not limited to, Monte Carlo simulation for generating multivariate random variates [24,30,31,32,33], wind power modeling [34], load computation in power network planning [35], and reliability analysis [36,37].
It should be noted that the great majority of multivariate TPNT applications are done under the assumption of known marginal statistical moments (i.e., product-moments and L-moments) and correlation coefficients. However, in real-life hydrologic applications, the amount of available data generally is not sufficiently large to reliably ascertain the true marginal probability distribution functions, statistical moments, and correlation coefficients. Therefore, the sample statistical moments and correlation coefficients used could be subject to sampling errors. In this study, a procedure is proposed to (1) optimally estimate multivariate TPNT coefficients by explicitly incorporating sampling errors associated with the sample moments and correlation coefficients, and (2) comply with a one-to-one monotonicity increasing relation between quantiles of the original and normal transformed variables. The procedure is illustrated by analyzing annual maximum rainfall data series involving seven different durations to establish at-site rainfall intensity–duration–frequency (IDF) and depth–duration–frequency (DDF) relationships.

2. Methods

2.1. Third-Order Polynomial Normal Transform (TPNT)

2.1.1. Univariate TPNT

By TPNT, a univariate non-normal random variable, X , is approximated by the standard normal variable, Z , in the form of a 3rd-order polynomial functional relation as [28]
X = TPNT ( Z   |   a 0 ,   a 1 ,   a 2 ,   a 3 ) = a 0 + a 1 Z + a 2 Z 2 + a 3 Z 3
where TPNT ( Z   |   a 0 ,   a 1 ,   a 2 ,   a 3 ) denotes the 3rd-order polynomial transform with a0, a1, a2, and a3 being the transformation coefficients. The TPNT coefficients can be determined by several methods of varying mathematical complexity. By preserving the first four product-moments, the TPNT coefficients are related to the first four product moments of the standardized variable, X = ( X μ x ) / σ x , as [28]
0 = a 0 + a 2
1 = a 1 2 + 6 a 1 a 3 + 2 a 2 2 + 15 a 3 2
γ x = 2 a 2 ( a 1 2 + 24 a 1 a 3 + 105 a 3 2 + 2 )
κ x = 3 + 24 [ a 1 a 3 + a 2 2 ( 1 + a 1 2 + 28 a 1 a 3 ) + a 3 2 ( 12 + 48 a 1 a 3 + 141 a 2 2 + 225 a 3 2 ) ]
in which γ x = skew coefficient; κ x = kurtosis of the original random variable X. Alternatively, the TPNT coefficients in Equation (2) can also be related to the first four L-moments as [38]
a 0 + a 2 = λ 1
0.5642 a 1 + 1.4104 a 3 = λ 2
0.5513 a 2   λ 3
0.0692 a 1 + 0.8078 a 3 = λ 4
in which λ m = the mth order L-moment [39] of the original non-normal random variable, X. Other than the above two moment-matching methods, TPNT coefficients can also be determined by the quantile-based least square method and the Fisher–Cornish asymptotic expansion (FC) method [40]. Chen and Tung [41] investigated the performance of different methods in determining the TPNT coefficients with regard to their accuracy and robustness in capturing the probabilistic features of the random variable X under the condition that the population distribution is known. It was found that, among the various methods for estimating TPNT coefficients, the L-moment based method is computational simplistic and can yield a satisfactory performance under a wide range of distribution conditions. The product-moment method can also yield a satisfactory normal transformation provided that accurate estimations of skew coefficient and kurtosis in Equations (3)–(6) can be made. However, when the statistical moments are to be estimated from finite data, the sample L-moments have been proven to be more stable and robust than those of product-moments [42], especially when the sample size is not large.
By referring to Equations (3)–(6), one also realizes that determining TPNT coefficients based on the product-moments requires solving a system of non-linear equations. It is expected that solving Equations (3)–(6) would be more difficult than solving L-moments based on Equations (7)–(10), which is linear. Sometimes, the solution to the system of non-linear equation may not be attainable. According to Equations (7)–(10), TPNT coefficients can be easily obtained in terms of L-moments as
a 0 = λ 1 1.8138 λ 3
a 1 = 2.2552 λ 2 3.9376 λ 4
a 2 = 1.8138 λ 3
a 3 = 0.1931 λ 2   + 1.5751 λ 4
In the transformation process, it is necessary to preserve probability content in both original space and standard normal space, i.e., F x ( x p ) = Φ ( z p ) = p . This implies that quantiles of the two variables should satisfy the following relationship:
x p = a 0 + a 1 z p + a 2 z p 2 + a 3 z p 3
where x p   and   z p = pth-order quantiles of random variable X and standard normal random variable, Z, respectively, that is, x p = F 1 ( p ) and z p = Φ 1 ( p ) . Furthermore, inherently embedded in Equation (15) is a requirement of one-to-one monotonically increasing relations between x p   and   z p . This, then, requires that TPNT coefficients must comply with the following conditions:
a 3 > 0   and   a 2 2 3 a 1 a 3 < 0 .
It should be noted that the TPNT coefficients obtained from solving Equations (3)–(6), Equations (7)–(10), or other methods mentioned above do not guarantee the compliance of the monotonicity condition stipulated in Equation (16). This is especially a major concern when sample statistics are used in determining TPNT coefficients.

2.1.2. Multivariate TPNT

The TPNT coefficients can be determined by preserving the statistical moments of individual random variables. Specifically, L-moments are used herein to determine the multivariate TPNT coefficients due to simple, linear functional relationships between the TPNT coefficients and the L-moments as shown in Equations (11)–(14). Furthermore, sample L-moments have several desirable sampling properties over the product-moments as proven by Hosking [42]. In the context of fitting the first four L-moments of a total of N correlated variables, Equations (7)–(10) can be re-written as
a 0 j + a 2 j = λ 1 j
0.5642 a 1 j + 1.4104 a 3 j = λ 2 j
0.5513 a 2 j   λ 3 j
0.0692 a 1 j + 0.8078 a 3 j = λ 4 j
in which λ m j = the mth order L-moment of the jth random variable Xj for j = 1, 2, …, N.
In addition to preserving marginal statistical moments of involved variables, multivariate TPNT must also simultaneously preserve the statistical dependence between random variables in the transformation. The correlation coefficient of any two correlated random variables, X j   and   X k , is imbedded in their 2nd-order cross-product moment of which Vale and Maurelli [29] had shown the explicit expressions in terms of TPNT coefficients as
C P j , k ( a j ,   a k ; ρ z j ,   z k ) = E [ X j X k ] = μ j μ k ρ x j ,   x k σ j σ k = ( 6 a 3 j a 3 k )   ρ z j ,   z k 3 + ( 2 a 2 j a 2 k )   ρ z j ,   z k 2 + [ ( a 1 j + 3 a 3 j ) ( a 1 k + 3 a 3 k ) ] ρ z j ,   z k + [ ( a 0 j + a 2 j ) ( a 0 k + a 2 k ) ]
in which ρ x j ,   x k ,   ρ z j ,   z k = correlation coefficient of random variables ( X j , X k ) and its equivalent ( Z j , Z k ) in normal space; μ j   and   σ j = mean and standard deviation of random variable X j , respectively. The correlation coefficient in the original scale, ρ x j ,   x k , is related to its counterpart in the normal space, ρ z j ,   z k , in a 3rd-order polynomial relationship through TPNT coefficients.
Upon the determination of TPNT coefficients for the two concerned random variables, the correlation coefficient in the normal space, ρ z j ,   z k , corresponding to that in the original space, ρ x j ,   x k , can be obtained by finding the real root of Equation (21). The mathematical relations between the two correlation coefficients are [20]
ρ x j ,   x k × ρ z j ,   z k > 0 ;   | ρ x j ,   x k | | ρ z j ,   z k |
Equation (21) is used repeatedly to solve for ρ z j ,   z k for all pairs of correlated random variables to establish the correlation matrix in multivariate normal space.

2.2. Optimization Framework for Determining Multivariate TPNT Coefficients

2.2.1. Objective Function

To determine the multivariate TPNT coefficients that best preserve the known values of L-moments, the least-square criterion is used in the study by which the objective function can be expressed as
M i n i m i z e m = 1 4 j = 1 N δ m j 2
where δ m j = a decision variable defining the deviation between the mth-order TPNT-based L-moments computed by the left-hand side of Equations (17)–(20) and the known values, λ m j , of the jth random variable, Xj. Of course, other forms of objective function, such as minimizing the sum of absolute deviations, can be used.

2.2.2. Constraints

Several constraints are essential to make sure that multivariate TPNT coefficients obtained are able to preserve the known statistical features and mathematical relationships of concerned random variables.
(a) Preservation of L-moments for the individual variable Xj:
The deviation δ m j in the objective function defining the degree of preserving the known values of the first-four L-moments of individual variable Xj can be written, according to Equations (17)–(20), as
a 0 j + a 2 j + δ 1 j = λ 1 j
0.5642 a 1 j + 1.4104 a 3 j + δ 2 j = λ 2 j
0.5513 a 2 j + δ 3 j = λ 3 j
0.0692 a 1 j + 0.8078 a 3 j + δ 4 j = λ 4 j .
Note that the value of δ m j is unrestricted-in-sign, meaning that its value can be negative, zero, and positive, depending on the relative magnitudes of TPNT-based L-moments and those of the known values.
In reality, statistical properties of a random variable are estimated from a finite number of sample data. Consequently, sample L-moments of random variables, Xj, are subject to uncertainty. In practice, two approaches are used to estimate sample L-moments: plotting position-based estimators and unbiased estimators. This study adopts the latter by which the first four unbiased estimators of L-moments, { λ 1 ,   λ 2 , λ 3 , λ 4 } , can be computed respectively as [42]
1 = ( n 1 ) 1 i = 1 n x i : n
2 = 1 2 ! ( n 2 ) 1 i = 1 n { ( i 1 1 ) ( n i 1 ) } x i : n
3 = 1 3 ! ( n 3 ) 1 i = 1 n { ( i 1 2 ) 2 ( i 1 1 ) ( n i 1 ) + ( n i 2 ) } x i : n
4 = 1 4 ! ( n 4 ) 1   i = 1 n { ( i 1 3 ) 3 ( i 1 2 ) ( n i 1 ) + 3 ( i 1 1 ) ( n i 2 ) ( n i 3 ) } x i : n
in which m = sample estimator of the mth-order L-moment, λ m ; x i : n = the ith ranked sample (ascending order) in a data of size n.
Suppose that the sampling distributions of sample L-moments are derived or approximated. Proper bounds can then be incorporated into Equations (24)–(27) for determining the suitable and probabilistically plausible TPNT coefficients for all N random variables X1, X2, …, XN. Assuming that the lower and upper bounds of the L-moments can be determined from their corresponding sampling distributions, constraint Equations (24)–(27) then can be modified as
1 j ( L ) a 0 j + a 2 j + δ 1 j 1 j ( U )
2 j ( L ) 0.5642 a 1 j + 1.4104 a 3 j + δ 2 j 2 j ( U )
3 j ( L ) 0.5513 a 2 j + δ 3 j 3 j ( U )
4 j ( L ) 0.0692 a 1 j + 0.8078 a 3 j + δ 4 j 4 j ( U )
for j = 1, 2, …, N. In Equations (32)–(35), m j ( U )   and   m j ( L ) are, respectively, the upper and lower bounds containing the unknown population mth-order L-moment of random variable Xj, λ m j . The derivation of bounds for unknown population L-moments is described in Section 2.3.1.
(b) Preservation of the monotonic probability–quantile relationship for individual variable Xj:
a 3 j > 0 ;   a 2 j 2 3 a 1 j a 2 j < 0 ,   for   j = 1 ,   2 ,   ,   N
(c) Preservation of the correlation between all pairs of different variables, Xj and Xk:
Based on Equation (21), any pair of two correlated random variables X j   and   X k must satisfy the following equation.
C P j , k ( a j ,   a k ; ρ z j ,   z k ) ( μ j μ k ρ x j ,   x k σ j σ k ) = 0 ,   for   all   variable   pairs   j k
where C P j , k ( a j ,   a j ; ρ z j ,   z k ) = E ( X j , X k ) = cross-product moment of variables, Xj and Xk, defined in Equation (21), which are functions of the corresponding TPNT coefficients a j = ( a 0 j ,   a 1 j ,   a 2 j ,   a 3 j ) and a k = ( a 0 k ,   a 1 k ,   a 2 k ,   a 3 k ) ; ρ x j ,   x k and ρ z j ,   z k = correlation coefficients of random variables, X j   and   X k , and their normal equivalents, Z j   and   Z k , respectively; μ j , σ j = mean and standard deviation of random variable X j , respectively.
Similarly, constraint Equation (37) on correlation coefficients can be modified as
r x j ,   x k ( L ) C P j , k ( a j ,   a k ; r z j ,   z k ) m j ( a j ) m k ( a k )   s j ( a j ) × s k ( a k ) r x j ,   x k ( U ) ,   for   all   variable   pairs   j   k
in which r x j ,   x k ( L ) ,   r x j ,   x k ( U ) = lower and upper bounds, respectively, of the unknown population correlation coefficient, ρ x j , x k , between the random variables Xj and Xk (see Section 2.3.2); r z j ,   z k = equivalent correlation coefficient of the random variables in the standard normal domain; m j ( a j ) ,   s j ( a j ) = TPNT-based estimation of mean and standard deviation of random variables Xj which can be computed according to Equations (3) and (4) as
m j ( a j ) = a 0 j +   a 2 j
s j 2 ( a j ) = a 1 j 2 + 6 a 1 j a 3 j + 2 a 2 j 2 + 15 a 3 j 2 .
Equation (38) can alternatively be expressed as
C P j , k ( a j ,   a k ; r z j ,   z k ) m j ( a j )   m k ( a k ) r x j ,   x k ( U )   s j ( a j )   s k ( a k ) 0
C P j , k ( a j ,   a k ; r z j ,   z k ) m j ( a j )   m k ( a k ) r x j ,   x k ( L )   s j ( a j )   s k ( a k ) 0
In summary, by considering sampling errors of sample L-moments and correlation coefficients, the optimization model to determine the most plausible TPNT coefficients for establishing multivariate relationships can be summarized as follows:
The objective function is expressed in Equation (23) or its variations, which is subject to the following constraints:
  • Equations (32)–(35) for preserving plausible L-moments (8 × N constraints);
  • Equation (36) for complying with a probability–quantile monotonic relationship (2 × N constraints);
  • Equations (41) and (42) for preserving plausible correlation coefficient (N × (N − 1) constraints); and
  • unrestrictive-in-sign of polynomial coefficients ( a 0 j ,   a 1 j ,   a 2 j ,   a 3 j ) and deviations δ m j .

2.3. Determination of Bounds for L-Moments and Correlation Coefficients

2.3.1. Bounds for L-Moments

To determine the bounds for L-moments, the sampling distributions corresponding to the sample L-moments are needed. For independent random samples of size n from a distribution function F x ( x ) having the mth-order population L-moment λ m , Hosking [39] showed that the statistic n 1 / 2 ( m λ m ) , with m being the sample L-moment of order m, is unbiased, having a sampling distribution asymptotically converge to the normal distribution with the mean zero and variance Λ m m . Therefore, the variance of the mth-order sample L-moment, m , has the variance of σ 2 ( m ) = Λ m m / n . For the first four orders of sample L-moment, the value of Λ m m can be computed by
Λ 11 = x < y ( y x ) 2 d u   d v
Λ 22 = x < y [ ( 3 4 v ) ( 1 4 u ) ] ( y x ) 2 d u   d v
Λ 33 = x < y [ ( 7 24 v + 18 v 2 ) ( 1 12 u + 18 u 2 ) ] ( y x ) 2 d u   d v
Λ 44 = x < y [ ( 13 + 84 v 150 v 2 + 80 v 3 ) ( 1 + 24 u 90 u 2 + 80 u 3 ) ]   ( y x ) 2 d u   d v
in which u = F x ( x )   and   v = F x ( y ) . To estimate the values of Λ m m based on the ranked sample observations, the double integration stated in Equations (43)–(46) can be carried out numerically as
Λ ^ 11 = 1 n ( n 1 ) i = 1 n 1 k = i + 1 n [ x ( k ) x ( i ) ] 2
Λ ^ 22 = 1 n ( n 1 ) i = 1 n 1 k = i + 1 n   { ( 3 4 p k : n ) ( 1 4 p i : n ) [ x ( k ) x ( i ) ] 2 }
Λ ^ 33 = 1 n ( n 1 ) i = 1 n 1 k = i + 1 n   { ( 7 24 p k : n + 18 p k : n 2 ) ( 1 12 p i : n + 18 p i : n 2 ) [ x ( k ) x ( i ) ] 2 }
Λ ^ 44 = 1 n ( n 1 ) i = 1 n 1 k = i + 1 n   { ( 13 + 84 p k : n 150 p k : n 2 + 80 p k : n 3 ) ( 1 + 24 p i : n 90 p i : n 2 +   80 p i : n 3 )   [ x ( k ) x ( i ) ] 2 }
in which x ( i ) = the ith ranked sample in ascending order, i.e., x ( 1 ) < x ( 2 ) < < x ( i ) < < x ( k ) < < x ( n ) ; p i : n = estimated cumulative probability for the ith ranked sample, i.e., Pr [ X x ( i ) ] , by using the well-known Weibull plotting position formula, p i : n = i / ( n + 1 ) . Makkonen [43] has shown that the Weibull plotting position formula [44] provides the best estimate for the underlying non-exceedance probability. The superiority of the Weibull formula gets more pronounced with a decreasing sample size. By adopting the normality distribution assumption, the α-confidence interval for the unknown population λ k can be obtained as
( m ( L )   ,   m ( U )   ) = ( m z α / 2 Λ m m n   ,   m + z α / 2 Λ m m n   )
in which z α / 2 = Φ 1 ( 1 α / 2 ) , a standard normal quantile with an exceedance probability of α / 2 , with Φ 1 ( · ) being the inverse standard normal CDF.

2.3.2. Bounds for Correlation Coefficients

To quantify the lower and upper bounds of a correlation coefficient, Fisher transform is often used by which the sampling distribution of the inverse hyperbolic tangent function of sample correlation approximately follows a normal distribution as [45,46]
t a n h 1 ( r ) = 1 2   l n ( 1 + r 1 r ) ~ N o r m a l ( 1 2   l n ( 1 + ρ 1 ρ )   ,   1 n 3 )  
in which r , ρ = sample and population correlation coefficients, respectively; n = number of sample pairs. With a specified confidence level α , the corresponding lower and upper bounds for the unknown population coefficient ρ can be obtained as
[ r ( L )   ,   r ( U ) ] = [ t a n h ( t a n h 1 ( r ) z α / 2 n 3 )   , t a n h ( t a n h 1 ( r ) + z α / 2 n 3 )   ]
where the hyperbolic tangent function is defined as t a n h ( θ ) = ( e 2 θ 1 ) / ( e 2 θ + 1 ) .

2.4. Solution Algorithm

A recursive procedure is proposed to solve the above optimization models for determining multivariate TPNT coefficients. The procedure consists of four steps of initialization, optimization, validation, and updating. Solution algorithm for determining multivariate TPNT coefficients considering sampling errors of L-moments and correlation coefficients is detailed below and outlined in Figure 1.
Step (1): Initialization—Since the problem involves a nonlinear optimization model, an initial solution would be needed. One straightforward and sound initial solutions for the TPNT coefficients, a j ( 0 ) , are those provided by Equations (11)–(14) for all variables j = 1, 2, …, N. The initial TPNT coefficients obtained this way will automatically satisfy the constraint Equations (24)–(27). However, they do not necessarily comply with the monotonicity constraint Equation (36). As for the initial normal correlation coefficients, rather than arbitrarily choosing a set of initial correlation coefficients, let { ρ z j ,   z k ( 0 ) } = { r x j ,   x k } in which r x j ,   x k is the sample correlation coefficient between random variables X j   and   X k . Alternatively, obtain a feasible set of initial { ρ z j ,   z k ( 0 ) } by solving the 3rd-order polynomial function of ρ z j ,   z k in Equation (21) according to the initially assumed TPNT coefficients, a j ( 0 ) .
Step (2): Optimization—Based on the initially adopted TPNT coefficients, a j ( 0 ) , and the normal transformed correlation coefficients { ρ z j , z k ( 0 ) } , solve the optimization model with objective function Equation (23) and constraint Equations (32)–(35), (36), and (41)–(42) for the optimal TPNT coefficients a j * = ( a 0 j * ,   a 1 j * ,   a 2 j * ,   a 3 j * ) for random variable X j ,   j = 1 , 2 , , N .
Step (3): Validation—From the optimal feasible TPNT coefficients a j * = ( a 0 j * ,   a 1 j * ,   a 2 j * ,   a 3 j * ) for j = 1, 2, …, N obtained from Step (2), determine the equivalent normal variates corresponding to the sample data by solving the 3rd-order polynomial function:
a 0 j * + a 1 j *   z j , i + a 2 j *   z j , i 2 + a 3 j *   z j , i 3 = x j , i   ,   for   all   j = 1 ~ N ;   i = 1 ~ n
where z j , i = unknown normal variate corresponding to the ith observation of the jth random variable x j , i under the optimal set of TNPT coefficients a j * = ( a 0 j * ,   a 1 j * ,   a 2 j * ,   a 3 j * ) . From the normal-transformed data series of two different variables, z j = ( z j , 1 ,   z j , 2 ,   ,   z j , n ) and z k = ( z k , 1 ,   z k , 2 ,   ,   z k , n ) , the corresponding correlation coefficient, ρ z j ,   z k * , in the normal space is calculated.
Step (4): Updating—Compare the discrepancies between the initialized ρ z j ,   z k ( 0 ) and validated ρ z j ,   z k * for all different pairs of concerned random variables. If the discrepancy in any pair of durations is judged to be significant, update the initial normal correlation as ρ z j ,   z k ( 0 ) = ρ z j ,   z k * and TPNT coefficients a j ( 0 ) = a j * , and the process from Steps (2)–(4) is repeated. Otherwise, the optimal solutions are obtained and the iteration stops.
With regard to the optimization step presented in Step (2), the sequential quadratic programming (SQP) algorithm is implemented [47]. The SQP tackles a nonlinear optimization problem by successively finding the approximated optimum solution to the quadratic programming (QP) representation of the original problem. The approximated solution is improved iteratively by solving the QP problem. Boggs and Tolle [48] elaborated some useful properties of the SQP algorithm. The subroutine “sqp.m” in Matlab is used in this study to solve the optimization model.

3. Numerical Example

In this section, at-site rainfall intensity–duration–frequency (IDF) and depth–duration–frequency (DDF) relations are established to demonstrate the proposed multivariate TPNT method and examine its general performance. Rainfall IDF relations are widely used in the planning, design, and management of hydrosystem infrastructures, such as stormwater sewer systems and detention basins [49,50]. Such relations at a given location involves at-site frequency analysis of annual maximum rainfall intensity (or depth) data of several selected durations. The conventional approach in rainfall frequency analysis chooses a proper parametric probability distribution model to individually fit the observed annual maximum rainfall data of different durations. The choice of a distribution model for the rainfall intensity–frequency relations is largely statistical without much physical justification [51].
By the conventional approach, resulting rainfall intensity–frequency curves of different durations could sometimes intersect within the probability range of practical application. The crossover phenomenon often occurs when data record length is relatively short. According to the physical reality, rainfall intensity–frequency curves of different durations should not crossover or intersect. Porras and Porras [52] attributed the occurrence of crossover of rainfall IDF curves to short record data of questionable representation in which a significant amount of sampling errors existed in the estimated rainfall quantiles by frequency analysis. One other plausible reason for the possible crossover of IDF curves is that frequency analysis of rainfall data is performed separately for each duration without considering the inter-correlations that are intrinsically embedded in rainfall data of different durations. Haktanir [53] earlier pointed out that rainfall frequency analysis of different durations in the process of establishing IDF relationships should not be performed independently of each other, but did not propose a mechanism to handle the correlation directly. Recently, Gräler et al. [54] applied D-vine copula, along with the generalized extreme value distributions, to derive rainfall IDF relationships based on rainfalls of five durations. You and Tung [55], under the TPNT framework, developed a constrained least square model to simultaneously considering rainfall data of seven durations for establishing at-site rainfall IDF relations. However, their model does not explicitly take into account the correlation among rainfall data of different durations.
The multivariate TPNT-based model presented above was applied to establish at-site rainfall IDF relations using annual maximum hourly rainfall data of various durations at a raingauge in Zhongli City of Taoyuan County, Taiwan. Annual maximum rainfall intensity data cover the record period of 1988–2015, but the year 1992 was excluded from the analysis due to long periods of registers with technical issues. Hence, only 27-year data (n = 27) with seven (N = 7) durations (i.e., 1, 2, 6, 12, 24, 48, and 72 h) are used in this illustration (see data in Table 1). The sample values of the mean, standard deviation, and first-four L-moments of rainfall data of different durations are tabulated in Table 2. Furthermore, the standard error values corresponding to the first four sample L-moments, according to Equations (47)–(50), are listed in Table 3. The sample correlation coefficients of all rainfall intensity pairs of different durations in the original and normal-transformed domains are shown in Table 3 and Table 4, respectively. Based on the information given in Table 2 and Table 4, one is able to define the lower and upper bounds for the L-moments and correlation coefficients according to the desired confidence level, α, by Equations (51) and (53), respectively. Table 5 lists the values of correlation coefficients in normal-transformed space, r z j ,   z k , provided by the solution to constraint Equations (41) and (42) in the optimization model.
Under different constraint types and confidence levels for the L-moments and correlation coefficients, the corresponding optimal multivariate TPNT coefficients can vary. With the confidence level of α = 90% for both L-moments and correlation coefficients, Table 6a–d list the optimal TPNT coefficients under four different constraint types, including “LM” for L-moments by Equations (32)–(35), “Mono” for monotonicity by Equation (36), “Corr” for correlation by Equations (41) and (42), and “NC” for no-crossover by Equation (56). Once the optimal TPNT coefficients associated with each rainfall duration are obtained from solving the multivariate TPNT model, the rainfall IDF relations, according to Equation (15), can be established as
i t , T * = a 0 , t * + a 1 , t * z T + a 2 , t * z T 2 + a 3 , t * z T 3
where i t , T * = estimated t-h, T-year rainfall intensity; a 0 , t * , a 1 , t * , a 2 , t * , and a 3 , t * = optimum TPNT coefficients corresponding to rainfall of duration t (h); zT is the standard normal quantile corresponding to return period T-year having an annual exceedance probability of 1 Φ ( z T ) = 1 / T .

4. Results and Discussions

By varying the value of zT for different return periods in Equation (55), in conjunction with the optimal TPNT coefficients listed in Table 6a–d, one can establish IDF curves as shown in Figure 2 and Figure 3. Part (a) of Table 6 and Figure 2, Figure 3 and Figure 4 (denoted by “LM”) shows the results from considering only the bounding constraints of L-moments, Equations (32)–(35). In fact, the optimal TPNT coefficients corresponding to each duration can be obtained separately from the exact solutions using sample L-moments in Equations (11)–(14). Note that the TPNT coefficients obtained from each rainfall duration at this stage do not necessarily comply with a one-to-one monotonic increasing relation of rainfall quantile and probability. This can be clearly seen in Table 6a for 1 and 2 h rainfalls for which the two monotonicity constraints are violated (shown by *). Part (b) (denoted by “LM/Mono”) shows the results by considering both L-moment constraints, Equations (32)–(35), and the monotonicity constraints, Equation (36), for each rainfall duration. In this case, both results presented in Parts (a) and (b) in Table 6 and Figure 2, Figure 3 and Figure 4 can be obtained separately by treating rainfall data of different durations without considering their inter-correlations. Results in Part (c), denoted by “LM/Mono/Corr,” were obtained by incorporating correlation constraints of rainfall data with different durations, Equations (41) and (42), in determining the multivariate TPNT coefficients.
To show the degree of goodness-of-fit of normal transformed rainfall data by the proposed multivariate TPNT procedure, a normal probability plot of 24 h rainfall data (after normal transformation) with the fitted line and 95% confidence band are shown in Figure 5 as an example. The goodness-of-fit test shown in Figure 5 was achieved by the Anderson–Darling test [56] by which the test statistic is 0.535 with a p-value of 0.155. Figure 5 represents the worst case among the seven durations considered. The range of p-value varies from 0.155 (for 72 h) to 0.933 (for 2 h), which are higher than the generally adopted significance level of 0.05. This indicates that the normal transform by the proposed multivariate TPNT procedure is quite adequate.
Note that the solution obtained up to this stage does not necessarily comply with the physical reality that rainfall intensity (depth) of a given return period is a decreasing (an increasing) function of duration. In other words, rainfall intensity/depth–frequency curves of different durations should not intersect or crossover each other. However, in the process of establishing rainfall IDF/DDF relationships, one does not know in advance if any two resulting two curves would intersect before the statistical model is developed. Therefore, a special set of intersections avoidance constraints are imposed in establishing the IDF curves:
( a 0 , t j a 0 , t k ) + ( a 1 , t j a 1 , t k ) z T * + ( a 2 , t j a 2 , t k ) z T * 2 + ( a 3 , t j a 3 , t k ) z T * 3 > 0 ,   for durations   t j < t k
where T * = upper limit of selected rainfall return period below which no crossover of IDF curves is permitted to occur; z T * = standard normal quantile obtainable from Φ 1 ( 1 1 T * ) . Hence, additional N − 1 no-crossover (NC) constraints are included in the optimization model to solve for multivariate TPNT coefficients. Part (d) results (denoted by “LM/Mono/Corr/NC”) show the rainfall IDF relations considering the NC constraints.
Figure 2 shows the rainfall intensity–duration curves corresponding to various frequencies. For this particular data set, by only preserving sample L-moments, Figure 2a reveals two unusual features for those curves when return period is high (say, ≥100 years). They are (1) curves that tend to converge together for rainfall duration in the vicinity of 1 h and (2) the relatively pronounced undulation of curves for medium and long duration. These features are indications of possible anomalies that should not appear in a reasonable rainfall IDF relation. The convergence of rainfall intensity–duration curves in Figure 2a, shown in a different form in intensity–frequency relation as Figure 3a, reveal that the 1 h curve (in red) clearly does not satisfy the monotonicity condition according to Equation (36), which requires a rainfall intensity quantile value to increase continuously with a return period (see also Table 4a). In fact, the 2 h intensity–frequency curve (in gold) also mildly violates the monotonicity requirement as the curve starts to bend down for high return periods. The violation of the monotonicity condition can also be observed in the form of the depth–frequency curve for a 1 h duration (see Figure 4a). In this circumstance, the non-monotonicity of the 1 h rainfall intensity–frequency relation produces a crossover with the 2 h curve shown in Figure 3a.
Interestingly, Figure 3a also reveals that 6 and 12 h rainfall intensity–frequency curves have a strong tendency to intersect as rainfall frequency increases. This tendency to intersect could be attributed to a relatively large undulation of intensity–duration curves in the range of 6–12 h when rainfall frequency increases (see Figure 2a). From 6 to 12 h, the gradient of intensity–duration curves flatten out for larger return periods. The empirical results show some evidence of improvement (in terms of a decrease in undulation, for large frequencies) when more constraints are considered. However, the improvement is not significantly enough to remove undulation. In practical engineering applications, the undulation of rainfall IDF curves such as those shown in Figure 2a–d is removed by fitting the estimated rainfall intensity–duration data by an empirical IDF model, such as Sherman’s equation [57].
It is clear that, by considering the monotonicity constraint, Equation (36), the crossover tendency of intensity–duration curves (see Figure 2b) in the vicinity of 1–2 h disappears (see also Table 6b), as does that of the 1 h and 2 h intensity–frequency curves in Figure 3b. Correspondingly, the concave down appearance of the 1 h depth–frequency curve and, to a lesser extent, the 2 h curve is corrected (see Figure 4b).
Notice that joint consideration of complying with L-moments and the monotonicity condition does not truly take into account the inter-correlations of rainfall intensity or depth with different durations. The appearance of undulation in the rainfall intensity–duration curves for medium and long durations (≥6 h), which satisfy the monotonicity condition, is not affected. Hence, the crossover tendency of 6 and 12 h intensity–frequency curves (see Figure 3b) and the actual intersection of 48 and 72 h depth–frequency curves (see Figure 4b) remain unchanged.
With further consideration of inter-correlations of rainfalls of different durations, Equations (41) and (42), the resulting rainfall IDF and DDF curves are shown in Part (c) of Figure 2, Figure 3 and Figure 4. Figure 2c shows that the rainfall intensity–duration curves in the range of short duration for a high return period completely remove the crossover tendency. Both Figure 3c and Figure 4c show that rainfall intensity–frequency and depth–frequency curves for 1 and 2 h are parallel to each other. Still, the 48 and 72 h rainfall depth–frequency curves intersect (see Figure 4c).
For illustration, this application artificially select T * = 5000-year in Equation (56) as the limiting frequency below which rainfall depth–frequency or intensity–frequency curves of any two durations are not allowed to intersect. The obvious results of imposing no-crossover constraint is that the 48 h rainfall depth–frequency curve in Figure 4d would not intersect with the 72 h curve.
As for the effect of confidence level, numerical results indicate that a feasible solution for TPNT coefficients may not exist when the confidence levels for the unknown true L-moments and correlation coefficients are too low. This is expected because the width of confidence interval shrinks toward the sample L-moments and correlation coefficients as the confidence level reduces. At a certain confidence level, the corresponding width of the confidence band might be too restrictive for the optimization model to find feasible TPNT coefficients that simultaneously satisfy the monotonicity constraints. How low the limiting confidence level is depends on the problem. In this numerical example, the limiting confidence level is about 70%, below which no feasible solution can be found for multivariate TPNT coefficients. On the other hand, a reasonable confidence interval allows one to obtain a suitable set of TPNT coefficients to approximate multivariate relations.

5. Summary and Conclusions

Statistical modeling and data analysis in hydrosystems engineering often encounter multiple correlated random variables following non-normal distributions. Due to the difficulty in establishing a full joint probability density function for the involved variables, most of the methods tackling multivariate problems preserve the marginal statistical properties (e.g., distributions or moments) of individual variables and their correlation structures. In this study, focus is placed on the third-order polynomial transform (TPNT) procedure, which relies on the preservation of marginal L-moments and correlations among variables. In particular, a general framework is presented to optimally determine multivariate TPNT coefficients incorporating the constraints that (1) preserve the statistical L-moments and correlations with explicit consideration of their associated sampling errors; (2) comply with a one-to-one monotonicity increasing relation between quantiles of the original and normal transformed variables. Other than the above basic constraints required to hold the statistical and mathematical validity of the TPNT method, additional constraints that are relevant to the problem at hand can be incorporated into the modeling framework. In the illustrative example of establishing rainfall intensity–duration–frequency (IDF) relations, the no-intersection constraints for rainfall depth–frequency curves of different durations, Equation (56), are introduced in the model formulation to ensure that the resulting IDF relationships comply with the physical reality. The proposed method not only solves for the suitable multivariate TPNT coefficients that satisfy the monotonicity condition for individual variables, but also produces the correlation coefficients between random variables in the normal space. At this stage, the proposed multivariate TPNT procedure has not gone through a formal mathematical testing for its performance under different scenarios of multivariate distributions, correlation structures, and sample sizes. However, the procedure is based on a good logic with sound statistical and mathematical theory. The results from the empirical application to establish at-site rainfall IDF relationships appear to be quite reasonable.

Author Contributions

Conceptualization: Y.-K.T. and L.W.Y.; formal analysis: Y.-K.T., L.W.Y., and C.S.Y.; investigation: Y.-K.T., L.W.Y., and C.S.Y.; writing-original draft preparation: Y.-K.T. and L.W.Y.; writing-review & editing: Y.-K.T., L.W.Y., and C.S.Y.; project administration: Y.-K.T. and C.S.Y.; funding acquisition: Y.-K.T. and C.S.Y.

Funding

This study was primarily supported by the Joint Cooperative Research Program managed by the National Research Foundation of Korea (2016K2A9A1A06922023) and the Ministry of Science & Technology of Taiwan (MOST 105–2923-E-009-004-MY2). Additional funding was received from the General Research Program of the Ministry of Science & Technology of Taiwan (106-2221-E-009 -067).

Acknowledgments

The authors wish to express their gratitude to the two anonymous reviewers for their thorough review and constructive comments.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Kotz, S.; Balakrishnan, N.; Johnson, N.L. Models and Applications. In Continuous Multivariate Distributions, 2nd ed.; Wiley and Sons Inc.: New York, NY, USA, 2005; Volume 1. [Google Scholar]
  2. Hutchinson, T.P.; Lai, C.D. Continuous Bivariate Distributions—Emphasizing Applications; Rumsby Scientific Publishing: Adelaide, South Australia, 1990. [Google Scholar]
  3. Yue, S. A bivariate gamma distribution for use in multivariate flood frequency analysis. Hydrol. Process. 2001, 15, 1033–1045. [Google Scholar] [CrossRef]
  4. Nadarajah, S.; Shiau, J.T. Analysis of extreme flood events for the Pachang River, Taiwan. Water Resour. Manag. 2005, 19, 363–374. [Google Scholar] [CrossRef]
  5. Der Kiureghian, A.; Liu, P.L. Structural reliability under incomplete probability information. J. Eng. Mech. 1986, 112, 85–104. [Google Scholar] [CrossRef]
  6. Genest, C.; Chebana, F. Chapter 30: Coupula Modeling in Hydrologic Frequency Analysis. In Handbook of Applied Hydrology, 2nd ed.; Singh, V.P., Ed.; McGraw-Hill Book Company: New York, NY, USA, 2017. [Google Scholar]
  7. Favre, A.C.; Adlouni, S.E.; Perreault, L.; Thiemonge, N.; Bobee, B. Multivariate hydrological frequency analysis using copulas. Water Resour. Res. 2004, 40, W01101. [Google Scholar] [CrossRef]
  8. Ganguli, P.; Reddy, M.J. Probabilistic assessment of flood risks using trivariate copulas. Theor. Appl. Climatol. 2013, 111, 341–360. [Google Scholar] [CrossRef]
  9. Chen, L.; Singh, V.P.; Guo, S.L.; Mishra, A.K.; Guo, J. Drought analysis using copulas. J. Hydrol. Eng. 2013, 18, 797–808. [Google Scholar] [CrossRef]
  10. Khedun, C.P.; Chowdhary, H.; Mishra, A.K.; Giardino, J.R.; Singh, V.P. Water deficit duration and severity analysis based on runoff derived from Noah land surface model. J. Hydrol. Eng. 2013, 18, 817–833. [Google Scholar] [CrossRef]
  11. Sadri, S.; Burn, D.H. Copula-based polled frequency analysis of droughts in the Canadian Prairies. J. Hydrol. Eng. 2014, 19, 277–289. [Google Scholar] [CrossRef]
  12. Tosunoglu, F.; Kisi, O. Joint modelling of annual maximum drought severity and corresponding duration. J. Hydrol. 2016, 543, 406–422. [Google Scholar] [CrossRef]
  13. Requena, A.I.; Mediero, L.; Garrote, L. A bivariate return period based on copulas for hydrologic dam design: Accounting for reservoir routing in risk estimation. Hydrol. Earth Syst. Sci. 2013, 17, 3023–3038. [Google Scholar] [CrossRef]
  14. Li, C.; Singh, V.P.; Mishra, A.K. A bivariate mixed distribution with a heavy tailed component and its application to single site daily rainfall simulation. Water Resour. Res. 2013, 49, 767–789. [Google Scholar] [CrossRef]
  15. Jun, C.; Qin, X.S.; Gan, T.Y.; Tung, Y.K.; De Michele, C. Bivariate frequency analysis of rainfall intensity and duration for urban stormwater infrastructure design. J. Hydrol. 2017, 553, 374–383. [Google Scholar] [CrossRef]
  16. Aas, K.; Czado, C.; Frigessi, A.; Bakken, H. Pair-copula constructions of multiple dependence. Insur. Math. Econ. 2009, 44, 182–198. [Google Scholar] [CrossRef] [Green Version]
  17. Shafaei, M.; Fard, A.F.; Dinpashoh, Y.; Mirabbasi, R.; Michele, D.C. Modeling flood event characteristics using d-vine structures. Theor. Appl. Climatol. 2017, 130, 713–724. [Google Scholar] [CrossRef]
  18. Vernieuwe, H.; Vandenberghe, S.; De Baets, B.; Verhoest, N.E.C. A continuous rainfall model based on vine copulas. Hydrol. Earth Syst. Sci. 2015, 19, 2685–2699. [Google Scholar] [CrossRef] [Green Version]
  19. Tosunoglu, F.; Singh, V.P. Multivariate Modeling of Annual Instantaneous Maximum Flows Using Copulas. J. Hydrol. Eng. 2018, 23, 04018003. [Google Scholar] [CrossRef]
  20. Qing, X. Generating correlated random vector involving discrete variables. Commun. Stat. Theory Methods 2017, 46, 1594–1605. [Google Scholar]
  21. Nataf, A. Determination des distributions dont les marges sont donnees. Comptes Rendus l’Academie Sciences 1962, 225, 42–43. [Google Scholar]
  22. Lebrun, R.; Dutfoy, A. An innovating analysis of the Nataf transformation from the copula viewpoint. Probab. Eng. Mech. 2009, 24, 312–320. [Google Scholar] [CrossRef]
  23. Liu, P.L.; Der Kiureghian, A. Multivariate distribution models with prescribed marginals and covariances. Probab. Eng. Mech. 1986, 1, 105–112. [Google Scholar] [CrossRef]
  24. Chang, C.H.; Tung, Y.K.; Yang, J.C. Monte Carlo simulation for correlated variables with marginal distributions. J. Hydraul. Eng. 1994, 120, 313–331. [Google Scholar] [CrossRef]
  25. Chen, H.F. Initialization for NORTA: Generation of random vectors with specified marginal and correlations. INFORMS J. Comput. 2001, 13, 312–331. [Google Scholar] [CrossRef]
  26. Li, H.S.; Lu, Z.Z.; Yuan, X.K. Nataf transformation based point estimate method. Chin. Sci. Bull. 2008, 53, 2586–2592. [Google Scholar] [CrossRef]
  27. Niaki, S.T.A.; Abbasi, B. Generating correlation matrices for normal random vectors in NORTA algorithm using artificial neural networks. J. Uncertain Syst. 2008, 2, 192–201. [Google Scholar]
  28. Fleishman, A.L. A method for simulating non-normal distributions. Psychometrika 1978, 43, 521–532. [Google Scholar] [CrossRef]
  29. Vale, C.D.; Maurelli, V.A. Simulating multivariate non-normal distributions. Psychometrika 1983, 48, 465–471. [Google Scholar] [CrossRef]
  30. Headick, T.C.; Sanwilowsky, S.S. Simulating correlated multivariate nonnormal distributions: Extending the Fleishman power method. Psychometrika 1999, 64, 25–35. [Google Scholar] [CrossRef]
  31. Chen, X.Y.; Tung, Y.K. Applications of TPNT in multivariate Monte Carlo simulation. In Water Resources Planning and Management; EWRI/ASCE: Philadelphia, PA, USA, 2003; pp. 23–26. [Google Scholar]
  32. Hodis, F.A. Simulating Univariate and Multivariate Nonnormal Distributions Based on a System of Power Method Distributions. Ph.D. Thesis, Southern Illinois University, Carbondale, IL, USA, 2008. [Google Scholar]
  33. Demirtas, H.; Hedeker, D.; Mermelstein, R.J. Simulation of massive public health data by power polynomials. Stat Med. 2012, 31, 3337–3346. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  34. Yang, H.; Zou, B. The point estimate method using third order polynomial normal transformation technique to solve probabilistic power flow with correlated wind source and load. In Proceedings of the Asia-Pacific Power and Energy Engineering Conference, Shanghai, China, 27–29 March 2012; pp. 1–4. [Google Scholar]
  35. Cai, D.; Shi, D.; Chen, J. Probabilistic load flow computation with polynomial normal transformation and Latin hypercube sampling. IET Gen. Trans. Distrib. 2013, 7, 474–482. [Google Scholar] [CrossRef]
  36. Chen, X.Y. Investigating Third-Order Polynomial Normal Transform and Its Applications to Uncertainty and Reliability Analyses. Master’s Thesis, Hong Kong University of Science and Technology, Hong Kong, China, 2002. [Google Scholar]
  37. Zhao, Y.G.; Lu, Z.H. Fourth moment standardization for structural reliability assessment. J. Struct. Eng. 2007, 133, 916–924. [Google Scholar] [CrossRef]
  38. Tung, Y.K. Polynomial normal transformation in uncertainty analysis. Appl. Probab. Stat. 1999, 1, 167–174. [Google Scholar]
  39. Hosking, J.R.M. The Theory of Probability Weighted Moments; IBM Research Report, RC12210; IBM: Yorktown Heights, NY, USA, 1986. [Google Scholar]
  40. Fisher, R.A.; Cornish, E.A. The percentile points of distribution having known cumulants. Technometrics 1960, 2, 209–225. [Google Scholar] [CrossRef]
  41. Chen, X.Y.; Tung, Y.K. Investigation of polynomial normal transform. Struct. Saf. 2003, 25, 423–445. [Google Scholar] [CrossRef]
  42. Hosking, J.R.M. L-Moments: Analysis and Estimation of Distributions Using Linear Combinations of Order Statistics. J. R. Stat. Soc. Ser. B Method 1990, 52, 105–124. [Google Scholar] [CrossRef]
  43. Makkonen, L.; Pajari, M.; Tikanmaki, M. Discussion on “Plotting positions for fitting distributions and extreme value analysis”. Can. J. Civil Eng. 2013, 40, 130–139. [Google Scholar] [CrossRef]
  44. Weibull, W. A statistical theory of the strength of materials. R. Swed. Inst. Eng. Res. Proc. 1939, 151, 1–45. [Google Scholar]
  45. Fisher, R.A. On the ‘probable error’ of a coefficient of correlation deduced from a small sample. Metron. 1921, 1, 1–32. [Google Scholar]
  46. Kennedy, J.B.; Neville, A.M. Basic Statistical Methods for Engineers and Scientists, 3rd ed.; Happer and Row Publishing: New York, NY, USA, 1986. [Google Scholar]
  47. Wilson, R.B. A simplicial Method for Convex Programming. Ph.D. Thesis, Harvard University, Boston, MA, USA, 1963. [Google Scholar]
  48. Boggs, P.T.; Tolle, J.W. Sequential quadratic programming. Acta Numer. 1995, 4, 1–51. [Google Scholar] [CrossRef]
  49. Akan, A.O.; Houghtalen, R.J. Urban Hydrology, Hydraulics, and Stormwater Quality; Wiley: Hoboken, NJ, USA, 2003. [Google Scholar]
  50. Sun, S.A.; Djordjevic, S.; Khu, S.T. Decision making in flood risk based storm sewer network design. Water Sci Technol. 2010, 64, 247–254. [Google Scholar] [CrossRef]
  51. Singh, V.P.; Strupczewski, W.G. On the status of flood frequency analysis. Hydrol Process. 2002, 16, 3737–3740. [Google Scholar] [CrossRef]
  52. Porras, P.J.S.; Porras, P.J., Jr. New perspective on rainfall frequency curves. J. Hydrol. Eng. 2001, 6, 82–85. [Google Scholar] [CrossRef]
  53. Haktanir, T. Divergence criteria in extreme rainfall series frequency analyses. Hydrol. Sci. J. 2003, 48, 917–937. [Google Scholar] [CrossRef] [Green Version]
  54. Gräler, B.; Fischer, S.; Schumann, A. Joint modeling of annual maximum precipitation across different duration levels. In EGU General Assembly Conference Abstracts; Discussion Paper SFB 823; EGU: Munich, Germany, 2016. [Google Scholar]
  55. You, L.; Tung, Y.K. Derivation of rainfall IDF relations by third-order polynomial normal transform. Stoch. Environ. Res. Risk Assess. 2018, 32, 2309–2324. [Google Scholar] [CrossRef]
  56. D’Agostino, R.B.; Stephens, M.A. Goodness-of-Fit Techniques; Marcel Dekker: New York, NY, USA, 1986. [Google Scholar]
  57. Sherman, C.W. Frequency and intensity of excessive rainfalls at Boston, Massachusetts. Transaction 1931, 95, 951–960. [Google Scholar]
Figure 1. Flow diagram showing the procedure of multivariate TPNT modeling considering sampling errors of sample statistics.
Figure 1. Flow diagram showing the procedure of multivariate TPNT modeling considering sampling errors of sample statistics.
Water 11 00490 g001
Figure 2. Multivariate TPNT modeling of rainfall intensity–duration relationships of varying return periods under different constraints which consider: (a) L-moments only; (b) L-moments and monotonicity; (c) L-moments, monotonicity and correlation; (d) L-moments, monotonicity, correlation and no crossover.
Figure 2. Multivariate TPNT modeling of rainfall intensity–duration relationships of varying return periods under different constraints which consider: (a) L-moments only; (b) L-moments and monotonicity; (c) L-moments, monotonicity and correlation; (d) L-moments, monotonicity, correlation and no crossover.
Water 11 00490 g002
Figure 3. Multivariate TPNT modeling of rainfall intensity–frequency relationships of varying durations under different constraints which consider: (a) L-moments only; (b) L-moments and monotonicity; (c) L-moments, monotonicity and correlation; (d) L-moments, monotonicity, correlation and no crossover.
Figure 3. Multivariate TPNT modeling of rainfall intensity–frequency relationships of varying durations under different constraints which consider: (a) L-moments only; (b) L-moments and monotonicity; (c) L-moments, monotonicity and correlation; (d) L-moments, monotonicity, correlation and no crossover.
Water 11 00490 g003
Figure 4. Multivariate TPNT modeling of rainfall depth–frequency relationships of varying durations under different constraints which consider: (a) L-moments only; (b) L-moments and monotonicity; (c) L-moments, monotonicity and correlation; (d) L-moments, monotonicity, correlation and no crossover.
Figure 4. Multivariate TPNT modeling of rainfall depth–frequency relationships of varying durations under different constraints which consider: (a) L-moments only; (b) L-moments and monotonicity; (c) L-moments, monotonicity and correlation; (d) L-moments, monotonicity, correlation and no crossover.
Water 11 00490 g004
Figure 5. Probability plot of normalized 24 h rainfall data.
Figure 5. Probability plot of normalized 24 h rainfall data.
Water 11 00490 g005
Table 1. Annual maximum rainfall intensities (mm/h) at Zhongli Station, Taiwan.
Table 1. Annual maximum rainfall intensities (mm/h) at Zhongli Station, Taiwan.
Year1 h2 h6 h12 h24 h48 h72 h
198848.041.018.19.04.62.92.2
198983.569.538.219.910.06.14.1
199089.051.325.114.08.15.73.8
199152.536.313.16.84.43.92.9
199373.050.525.815.58.14.33.0
199444.531.019.111.78.75.03.3
199587.559.027.113.77.04.73.5
199680.540.314.010.38.04.63.1
199735.522.816.411.06.33.72.7
199870.535.321.211.17.84.53.0
199958.529.312.27.75.52.92.0
200043.525.819.814.011.28.25.8
200149.543.528.021.514.312.48.7
200235.028.312.46.95.02.71.8
200330.520.59.27.14.72.61.7
200451.535.319.413.39.66.14.1
200538.030.314.89.37.64.63.7
200648.532.319.713.87.75.34.6
200765.560.524.912.68.46.55.3
200842.032.313.59.57.15.54.8
200933.029.314.910.76.23.22.6
201047.535.013.69.37.03.62.4
201171.548.821.012.67.13.72.6
201272.047.843.434.317.68.85.9
201352.046.025.814.57.66.74.5
201441.026.315.011.37.94.93.6
201535.525.515.610.35.63.22.1
Table 2. Sample moments (in mm/h) of rainfall intensity data.
Table 2. Sample moments (in mm/h) of rainfall intensity data.
Moments1 h2 h6 h12 h24 h48 h72 h
μ ^ 54.8038.2620.0412.657.895.053.63
σ ^ 17.812.47.95.62.92.21.5
1 54.8038.2620.0412.657.895.053.63
2 10.217.024.252.671.461.130.83
3 1.6731.3931.1100.8790.3890.3170.200
4 0.3220.7150.7060.8630.4600.2460.148
Table 3. Standard errors (in mm/h) of sample L-moments of rainfall intensity data.
Table 3. Standard errors (in mm/h) of sample L-moments of rainfall intensity data.
Std. Error1 h2 h6 h12 h24 h48 h72 h
s e ( 1 ) 3.4232.3861.5191.0710.5540.4150.296
s e ( 2 ) 1.1100.9520.7500.6820.3240.2260.148
s e ( 3 ) 1.0050.6600.4470.3890.1590.1270.089
s e ( 4 ) 0.7130.4040.2090.1820.0880.0660.050
Table 4. Sample correlation coefficients between rainfall intensity of different durations.
Table 4. Sample correlation coefficients between rainfall intensity of different durations.
Duration1 h2 h6 h12 h24 h48 h72 h
1 h1.000
2 h0.8261.000
6 h0.6240.7601.000
12 h0.4210.4960.9101.000
24 h0.3060.3350.7620.9151.000
48 h0.2100.3380.6350.7320.8691.000
72 h0.1610.3240.5760.6710.8120.9761.000
Table 5. Correlation coefficients of normal-transformed rainfall intensity of different durations.
Table 5. Correlation coefficients of normal-transformed rainfall intensity of different durations.
Duration1 h2 h6 h12 h24 h48 h72 h
1 h1.000
2 h0.8491.000
6 h0.6540.7681.000
12 h0.5110.5690.9181.000
24 h0.4300.4400.7360.8701.000
48 h0.3960.4620.7130.7960.9031.000
72 h0.3230.4380.6700.7550.8540.9791.000
Table 6. Multivariate TPNT coefficients obtained under different constraints with α = 0.90.
Table 6. Multivariate TPNT coefficients obtained under different constraints with α = 0.90.
(a) Constraints: L-moments only (LM)
TPNT Coefficients1 h2 h6 h12 h24 h48 h72 h
a 0 51.7635.7318.0311.067.194.473.27
a 1 21.7613.026.802.621.471.581.29
a 2 3.032.532.011.590.700.580.36
a 3 −1.465−0.2300.2920.8430.4440.1700.072
a 3 > 0 −1.46 *−0.23 *0.290.840.440.170.07
a 2 2 3 a 1 a 2 < 0 104.8 *15.4 *−1.90−4.09−1.46−0.47−0.15
(b) Constraints: L-moments and Monotonicity (LM/Mono)
1 h2 h6 h12 h24 h48 h72 h
a 0 51.9935.8118.0311.067.194.473.27
a 1 17.5211.986.802.621.471.581.29
a 2 2.802.452.011.590.700.580.36
a 3 0.1500.1670.2920.8430.4440.1700.072
a 3 > 0 0.150.170.290.840.440.170.07
a 2 2 3 a 1 a 2 < 0 −0.01−0.01−1.9−4.1−1.5−0.5−0.1
(c) Constraints: L-moments, Monotonicity, and Correlation (LM/Mono/Corr)
1 h2 h6 h12 h24 h48 h72 h
a 0 51.7635.8118.0311.067.194.473.27
a 1 17.1411.986.802.621.471.581.29
a 2 3.032.452.011.590.700.580.36
a 3 0.3820.1670.2920.8430.4440.1700.072
a 3 > 0 0.380.170.290.840.440.170.07
a 2 2 3 a 1 a 2 < 0 −10.4−0.01−1.9−4.1−1.5−0.5−0.1
(d) Constraints: L-moments, Monotonicity, Correlation, and No Crossover (LM/Mono/Corr/NC)
1 h2 h6 h12 h24 h48 h72 h
a 0 51.9935.8118.0111.067.194.483.26
a 1 17.5211.986.802.621.471.611.25
a 2 2.802.452.031.590.710.570.37
a 3 0.1500.1670.2920.8430.4440.1590.089
a 3 > 0 0.150.170.290.840.440.160.09
a 2 2 3 a 1 a 2 < 0 −0.01−0.01−1.8−4.1−1.5−0.4−0.2
Note: * indicates a violation of monotonicity condition.

Share and Cite

MDPI and ACS Style

Tung, Y.-K.; You, L.; Yoo, C. Third-Order Polynomial Normal Transform Applied to Multivariate Hydrologic Extremes. Water 2019, 11, 490. https://doi.org/10.3390/w11030490

AMA Style

Tung Y-K, You L, Yoo C. Third-Order Polynomial Normal Transform Applied to Multivariate Hydrologic Extremes. Water. 2019; 11(3):490. https://doi.org/10.3390/w11030490

Chicago/Turabian Style

Tung, Yeou-Koung, Lingwan You, and Chulsang Yoo. 2019. "Third-Order Polynomial Normal Transform Applied to Multivariate Hydrologic Extremes" Water 11, no. 3: 490. https://doi.org/10.3390/w11030490

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop