# Simulating Nonhomogeneous Non-Gaussian Field by Using Iterative Rank-Dependent Reordering versus Translation Process-Based Procedure

Authors: H. P. Hong; M. Y. Xiao
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

We compare two commonly used procedures, namely, the iterative rank-dependent reordering (IRDR) procedure and the translation process based procedure, for simulating homogeneous/nonhomogeneous non-Gaussian fields. We identify the limitations and the implicit assumptions of the procedures. We provide a new interpretation of the IRDR procedure, point out that there is no guarantee that the algorithm converges, and suggest modifications in terms of the initial samples, iteration involving decomposition, and convergence requirement to the IRDR procedure for it to become more efficient and robust. The numerical results show that, depending on the prescribed marginal probability distribution, the use of the IRDR procedure may not achieve a prescribed correlation function, a feature that is well-known if the translation process (i.e., Nataf translation system) based procedure is employed. It is shown that the performance of the modified IRDR procedure is comparable to that of the translation process based procedures in terms of limitations and matching the prescribed correlation function. The numerical results also show that the suggested modifications to IRDR in the present study make the algorithm more efficient and robust.


---

## Body

## 1. Introduction

Probabilistic assessment of the responses of engineering systems involves modeling and simulating nonhomogeneous non-Gaussian random vector and field [1–3]. The stochastic modeling could involve random variables with different prescribed non-Gaussian marginal probability distributions and correlation structures. One of the approaches [4–8] takes advantage of the Nataf translation system [9] that approximates a system of non-Gaussian random variables or field by a Gaussian system. The adoption of this approximation facilitates the simulation of the nonhomogeneous and non-Gaussian fields. A key aspect of this approach is to find an underlying or equivalent Gaussian Pearson product-moment correlation coefficient (which will be simply referred to as correlation or Pearson correlation) to a prescribed correlation coefficient for the non-Gaussian correlated variables. This evaluation is based on the solution of a two-dimensional integral equation defining the Nataf translation system for each pair of random variables. There are two well-known possible weaknesses of this approach. The first one is that the underlying Gaussian correlation may not exist for the prescribed combination of non-Gaussian random variables and correlation coefficient. The problem could be dealt with by using the closest approximation to the solution to the integral equation. The second problem is that the correlation matrix (or function) obtained by solving the integral equation may not be positive semidefinite or cannot be decomposed by Cholesky decomposition. This drawback could be overcome by using the nearest correlation coefficient matrix [10, 11]. The use of the nearest correlation coefficient matrix was considered by Kim and Shields [12] in assessing the underlying equivalent Gaussian correlation matrix for the non-Gaussian random field, where the integral equation was solved numerically. The availability of the underlying Gaussian correlation matrix or function simplifies the simulation of the non-Gaussian vector based on the Gaussian copula [13] and the simulation of the non-Gaussian random field.Another approach relies on the rank-dependent reordering (RDR) or shuffling of sampled non-Gaussian random variables to match the prescribed correlation structure [14]. This RDR procedure assumes that the Spearman correlation and Pearson correlation are a good approximation of each other. The use of rank-dependent reordering was considered, among others, by Florian [15]; Olsson et al. [16]; and Yang [17]. Florian [15] emphasized the use of reordering to update Latin hypercube sampling. As the Spearman rank correlation differs from the Pearson correlation, a modification was suggested in Yang [17], which applies a particle swarm optimization algorithm to minimize the difference in the Pearson correlation calculated from the samples to its target. More recently, Zheng et al. [18] proposed simulating the nonhomogeneous non-Gaussian field by iteratively reordering the samples which are initially generated using the marginal distributions; the reordering is carried out according to their ranking such that the correlation matrix of the reordered samples matches the prescribed correlation matrix. Although no link between their procedure and that of Iman and Conover [14] was pointed out, this procedure calls for iteratively applying RDR. Moreover, the limitation of using the iterative RDR (IRDR) procedure to simulate the nonhomogeneous and non-Gaussian field is not thoroughly explored. A comparison of the performance of using the iterative rank-dependent reordering procedure versus using the popular translation process (i.e., Nataf translation system) based procedure to simulate a non-Gaussian field is unavailable in the literature.For completeness, we note that other popular methods to simulate the random field include the spectral representation method (SRM) and the method based on the Karhunen-Loeve (KL) expansion. A random field represented by the KL expansion is controlled by the random coefficients [1, 12, 19–21]. Since the direct use of the KL expansion with random coefficients leads to the Gaussian field, Phoon et al. [19, 20] proposed an iterative procedure that modifies and shuffles the random expansion coefficients based on the ranking to simulate the non-Gaussian process. Kim and Shields [12] used the underlying Gaussian correlation for the non-Gaussian field and the KL expansion to simulate the random field. SRM [22–24] is based on the power spectral density (PSD) function of the random field, where the relation between the PSD function and the covariance function of the field is established based on the Wiener–Khintchine theorem. A discussion of using SRM for conditional simulation of nonstationary processes was given by Cui and Hong [25]. SRM was extended to simulate nonstationary/nonhomogeneous and non-Gaussian fields by many researchers, including Yamazaki and Shinozuka [26]; Grigoriu [27]; Deodatis and Micaletti [28]; Ferrante et al. [29]; Shields and Deodatis [30]; and Hong and Cui [31]. These studies use iterative procedures to find the underlying Gaussian PSD function; the samples of the non-Gaussian field are obtained by transforming the samples generated from the underlying Gaussian PSD function. Without the need to find the underlying Gaussian PSD function, iterative algorithms are presented by Schreiber and Schmitz [32] and Masters and Gurley [33] to generate the surrogate or to simulate the homogeneous non-Gaussian field and by Hong et al. [34–36] and Cui and Hong [37] to simulate the nonhomogeneous and non-Gaussian fields or vectors. These algorithms involve reordering samples such that the PSD function of the sampled field matches the prescribed PSD function. It is acknowledged that the above review of techniques to simulate the non-Gaussian field is not exhausting, considering different complex mathematical modeling and techniques employed on this topic in the literature; a comprehensive review is beyond the scope of the present study.The aim of the present study is twofold. First, we identify whether these two procedures share similar limitations in simulating homogeneous/nonhomogeneous non-Gaussian fields and the implicit assumptions of the procedures. Second, we provide a new interpretation of the IRDR procedure, point out that there is no guarantee that the algorithm converges, and suggest modifications in terms of the initial samples, iteration involving decomposition, and convergence requirement to the IRDR procedure for it to become more efficient and robust. In the following, we first summarize the RDR procedure and provide our interpretation of the IRDR. This serves to identify necessary modifications to the algorithm for it to become more efficient and robust. We then succinctly summarize the properties of the underlying Gaussian correlation based on the Nataf translation system [6, 7] to show the limitation of the Nataf translation system based procedure. This is followed by numerical examples comparing the performance of the IRDR procedure and the translation process based procedure and by conclusions.


## 2. Simulation Based on the IRDR Procedure and Suggested Modifications

Consider a one-dimensional random fieldY(x) over the domain of x with prescribed marginal cumulative distribution function, Fyx;xyx;x, and covariance or correlation function. For simplicity and without loss generality, it is considered, in the following, that Y(x) has zero mean and unit variance since the standard deviation and the mean values can be easily considered by scaling and shift.In many practical applications, the data collection and data analysis for a random field are carried out at discrete points with a given sampling interval. Letxp=x1,…xpT represent a vector of x of size p×1 at the discrete points, where the superscript denotes the transpose and x,j represents the j-th element of the vector. It is considered that the correlation function of the field can be adequately represented by the correlation matrix (or covariance matrix) CρY of size p×p.One could sample a matrixU of size p×m with its (i, j)-th element representing the value of the independent standard uniformly distributed random variable. A matrix U of size p×m with the (i, j)-th element zi,j=Φ−1ui,j can be formed. If Y(x) is a Gaussian random field, the simulation of m samples of the field can be obtained using the well-known procedure based on the Cholesky decomposition,(1)Y=L0Z,where each column of Y represents a sample of the random field and L0 denotes the lower triangular matrix obtained by decomposing CρY using the Cholesky decomposition. This decomposition can be carried out only if CρY is a positive semidefinite matrix [38]. The use of equation (1) ensures that the correlation matrix calculated from Y equals the prescribed target CρY.To simulate non-Gaussian correlated random variables with a prescribed correlation matrix and to reduce potential spurious correlation, Iman and Conover [14] proposed a distribution-free RDR procedure in the context of Latin hypercube sampling. The procedure assumes that the Spearman (rank-dependent) correlation could be approximated by the Pearson (product-moment) correlation CρY, which could be decomposed by using the Cholesky decomposition to obtain L0. Multiple sets of the random vector or field Y are sampled according to the marginal probability distribution (Gaussian or non-Gaussian) and Latin hypercube sampling technique. A row-wise matrix R containing the row-wise ranking of Y is formed. The samples of the correlated variables are then formed by reordering Y according to the row-wise ranking of the matrix L0R.The RDRS procedure has been exploited in probabilistic analysis of engineering systems [15–17, 39]. Based on the RDR procedure, Olsson et al. [16] suggested simulating standard multivariate normal variate with a prescribed correlation matrix and then mapping the samples to the non-Gaussian space through the probability transformation. However, Olsson et al. [16] did not elaborate on which correlation matrix should be used. Following Iman and Conover [14], one could assume that the Pearson correlation ρ equals the target Spearman correlation ρs. For bivariate normally distributed random variables, the relation between ρs and ρ is (e.g., [13])(2)ρs=6πarcsinρ2.The absolute difference betweenρs and ρ is less than 0.0181, which occurs for ρ equal to ±0.594. However, since ρs is usually unavailable and equation (2) is not applicable for the non-Gaussian field, it seems that ρ in the non-Gaussian space is employed to carry out the simulation as a practical measure. Based on the above considerations, the RDR procedure followed by Olsson et al. [16] for Latin hypercube sampling may be described as follows:(1)
CalculateL0 from CρY and sample Z based on the standard normal distribution and the Latin hypercube or stratified sampling scheme but without considering correlation.(2)
CalculateAρZ=1/mZZT, and apply the Cholesky decomposition to AρZ to obtain the lower triangular matrix BZ.(3)
Calculate(3)Z∗=L0BZ−1Z.(4)
Determine the matrixR containing the ranks based on Z∗ according to its row-wise ranking and reorder the samples in Z according to R (so a newly ordered Z is obtained).(5)
MapZ to Y by using the probability transformation, where yi,j=Fyxi;xi−1Φzi,j, and yi,j and zi,j denote the (i, j)-th element of Y and Z, respectively.Several observations are worth mentioning. First, we interpretZ and Z∗ in Steps 1 to 3 as projections of another underlying sampled Gaussian field Z0; they are given by Z=BZZ0 and Z∗=L0Z0. From the first projection, we find the underlying sampled Gaussian field Z0, and its use in the second projection leads to equation (3). While the projected Z leads to a covariance matrix that differs from the prescribed covariance matrix, Z∗ leads to the prescribed covariance matrix in the Gaussian space, that is,(4)Z∗Z∗T=L0BZ−1ZL0BZ−1ZT=CρY.Because the row-wise rank ofR, Z∗, and the resulting matrix Z in Step 5 are identical, the Spearman correlation matrix from each of the three matrices is identical by definition. However, this Spearman correlation matrix is likely to be different from the Pearson correlation matrix CρY. In other words, the Pearson correlation from Z is not equal to CρY, although the difference is small if Z is Gaussian. A critical issue with this simulation procedure is that there is no basis to infer if the Pearson correlation matrix obtained from Y equals CρY if Y is strongly non-Gaussian. It can be shown that the steps shown in the above are applicable, whether Z in Step (1) is sampled from independent or correlated normally distributed random variables. In fact, if Z in Step (1) is sampled by considering the correlation matrix CρY, the reordering is unnecessary.Rather than carrying out the sampling in the Gaussian space and mapping it to a non-Gaussian space, Zheng et al. [18] proposed to carry out the initial sampling according to its marginal distribution directly and reordering the samples iteratively. There is no restriction that a stratified simulation is to be considered. More specifically, their procedure consists of generating samples Y based on the marginal probability distribution of Y(x) and iteratively approximating the prescribed correlation matrix by reordering. Although it was not pointed out in Zheng et al. [18], as will be seen, the similarity between their procedure in each iteration and that given by Iman and Conover [14] is striking. The steps of this IRDR procedure are as follows:(1)
CalculateL0 based on the prescribed CρY and sample Y of size p×m according to the prescribed marginal distribution of Y(x).(2)
CalculateAρY=1/mYYT, and apply the Cholesky decomposition to AρY to obtain the lower triangular matrix, BY, where AρY=BYBYT. The calculation of AρY is based on the biased estimate, and the mean and standard deviation of the samples are considered to be equal to zero and unit. This is adequate since m is considered to be very large. Otherwise, the calculation of AρY should be modified accordingly.(3)
CalculateY∗=L0BY−1Y.(4)
DetermineR based on Y∗ according to the row-wise ranking; reorder Y following the row-wise ranking in R (and retain it as new Y). Calculate AρY=1/mYYT.(5)
IfAρY−CρY>εCρY, iterate Steps (3) to (5), where ε is a preselected tolerance value, and the Frobenius norm for the matrix is adopted in the present study.We note from the above that the IRDR procedure, except the iteration and no stratified sampling, is the same as the RDR procedure. Immediately after Step (4) in each iteration, the Spearman correlation matrices fromR, Y∗, and Y are identical by definition. However, there is no guarantee that the algorithm converges for a selected tolerance. We make this assertion based on two factors. First, equation (2) is not applicable for non-normally distributed random variables, and the difference between the Spearman correlation and Pearson correlation, in general, is unknown. Second, the objective of the iterations is to reduce the difference between the Pearson correlation matrix calculated from the samples and the target (i.e., AρY−CρY), but the reordering is based on and aimed at matching Spearman correlation. Similar to equation (4), it is straightforward to show that the Pearson correlation matrix from Y∗ equals the prescribed CρY. However, if Y∗ is used to represent the samples of the field, it may not match the prescribed marginal distribution of the non-Gaussian field; if Y is used to represent the samples of the field, it may not match the prescribed Pearson correlation. One may also use the difference between Y∗ and Y as a convergence measure. For the cases when the adopted convergence criterion in Step (5) is not achieved (examples will be presented shortly), we suggest assigning a maximum iteration number for the procedure and using Y corresponding to the minimum AρY−CρY during the iteration as the samples of the random field.Another problem with the algorithm for practical applications is that the assigned correlation matrix (which could be calculated using the prescribed correlation function)CρY may not be positive semidefinite. This will be elaborated by using numerical examples in the subsequent sections. A matrix must be at least positive semidefinite to be decomposed by the Cholesky decomposition [38]. To ameliorate this problem, we suggest substituting CρY with its corresponding nearest correlation matrix [10, 11]. The same could be done for AρY if it is not positive semidefinite. The use of the nearest correlation matrix makes the algorithm more robust, at least for when CρY is not positive semidefinite; the algorithm given in Qi and Sun [11] is employed for the numerical analysis to be carried out in the following sections.Lastly,Y in Step (1) could be simulated with or without an assigned correlation structure and the prescribed marginal distribution. The use of samples with an assigned correlation matrix could be advantageous in accelerating the convergence. For example, we could sample Y in Step (1) based on Gaussian copula with the prescribed marginal distribution and the assigned correlation matrix equal to CρY. This results in each element of Y in Step (1) being equal to yi,j=Fyxi;xi−1Φzci,j, where zci,j is the (i, j)-th element of Zc and Zc=L0Z.The flowchart describing the modified IRSD procedure based on the above three considerations is shown in Figure1.Figure 1
Flowchart of the modified iterative rank-dependent reordering procedure to simulate the nonhomogeneous and non-Gaussian field.


## 3. Underlying Gaussian Correlation for Non-Gaussian Field Based on the Nataf Translation System

The Nataf translation system [9] approximates two correlated non-Gaussian variables by bivariate normally distributed random variables, where the underlying Gaussian correlation ρG is determined by solving the following integral equation [4, 6, 40]:(5)cjkρG=∫−∞∞∫−∞∞FYxj;xj−1Φψj;xjFYxk;xk−1Φψk;xkϕ2ψj,ψk;ρGdψjdψk,where cjkρG is the Pearson correlation coefficient between y(xj) and y(xk), cjkρG and ϕ2ψj,ψk;ρG are the joint probability density function of two correlated standard normal variables, and Φ is the standard normal distribution function. Several notable relations between cjkρG and ρG are as follows [6, 7]:(a)
cjkρG≤ρG. cjkρG is a strictly increasing function of ρG for continuous random variables with a finite variance value. cjkρG≥0≤0 if ρG≥0≤0(b)
The maximum and minimum values ofcjkρG, denoted as ρ⌢NG,jk and ρ⌣NG,jk, equal cjk1 and cjk−1, respectivelyAs a consequence of these propertiescjkρG=0 if ρG=0. In order to find the underlying Gaussian symmetric correlation matrix with the diagonal elements equal to one for the non-Gaussian field, we could first determine ρ⌢NG,jk and ρ⌣NG,jk. These values are used to obtain a modified CρY, denoted as CρY∗, with its (j,k)-th element ρY∗xj,xk defined by(6)ρY∗xj,xk=minρYxj,xk,ρ⌢NG,jk,ifρYxj,xk≥0maxρYxj,xk,ρ⌣NG,jk,ifρYxj,xk≤0,where ρYxj,xk denotes the (j,k)-th element of CρY. For cases where ρY∗xj,xk equals ρ⌢NG,jk or ρ⌣NG,jk and j≠k, the underlying Gaussian correlation ρG,jk equals 1 and −1, respectively. For the remaining cases and j≠k, ρG,jk could be obtained by solving equation (5) iteratively with cjkρG,jk=ρY∗xj,xk, noting that the solution to ρG,jk is bracketed by cjkρG,jk≤ρG,jk≤1 if 0<cjkρG, and cjkρG,jk≤ρG,jk<0 if cjkρG,jk<0. If j=k, ρG,jk=1.The underlying Gaussian correlation matrixCG is formed by ρG,jk. If CG is not a positive semidefinite matrix, CG is replaced by its corresponding nearest correlation matrix C˜G [10, 11].The use of the Nataf translation system with the KL expansion to simulate the random field was considered by Kim and Shields [12]. Rather than using the KL expansion, for the field that is represented at p points (i.e., at xp), we could simply sample the Gaussian field based on the Cholesky decomposition of C˜G and then map the sampled Gaussian field to the non-Gaussian field to obtain Y with its (i, j)-th element yi,j=Fyxi;xi−1ΦzGi,j, where zGi,j is the (i, j)-th element of ZG=LGZ, and LG denotes the lower triangular matrix obtained by applying the Cholesky decomposition to CG or C˜G. For simplicity, this procedure is referred to as the Nataf translation system based (NTS) procedure in the following.It must be emphasized that the exposition so far is focused on the one-dimensional random field. Since the simulation of the random field is carried out at discrete points by using the NTS procedure and the IRDR procedure, these procedures can be directly extended to then-dimensional field by viewing and arranging all the discrete points in multidimension as a long vector.


## 4. Comparison of Matched Correlation, Underlying Gaussian Correlation, and Spearman Correlation

### 4.1. Prescribed Conditions for the Random Fields

A set of eight correlation functions were considered by Kim and Shields [12] to test their procedure with the use of the KL expansion to simulate the random field. For each prescribed correlation function, they considered two marginal distributions: one is weakly non-Gaussian, and the other is strongly non-Gaussian. These cases are listed in Table 1. Cases 1 to 4 describe homogeneous non-Gaussian fields, while Cases 5 to 8 describe nonhomogeneous non-Gaussian fields. The consideration of the beta distribution with the parameters given in Table 1 results in weakly non-Gaussian fields. The consideration of the shifted lognormal distribution with the distribution parameters assigned in Table 1 leads to strongly non-Gaussian fields. A simple spectral analysis indicates that Cases 3, 4, 7, and 8 represent narrow-band fields. The use of KL expansion, which was considered in Kim and Shields [12] in connection with the NTS procedure and in Zheng et al. [18] and in Xiao and Hong [41] in connection with the IRDR procedure, is beyond the scope of the present study.Table 1
Definition of the cases (see [12]).
CasePrescribed autocorrelation function for homogeneous cases (1 to 4) or autocovariance function for nonhomogeneous cases (5 to 8)Distribution type and notesBeta distributionLognormal distribution1Cρ1xj,xk=exp−xj−xkFy=ΓpΓq/Γp+q∫0uzp−11−zq−1dzp = 4, q = 2. For Cases 1 to 4, u=y+3.74/1.87+3.74. The field is defined for x within 0 to 2 and discretized using Δx = 0.02.For Cases 5 to 8, u=y−ymin/ymax−ymin, ymin=μbxj−σbxjpp+q+1/q, ymax=μbxj+σbxjqp+q+1/p. The mean μbxj=0 and the standard deviation σbxj is given by the autocovariance function. The field is defined for x within 0 to 1 and discretized using Δx = 0.01Fy=Flny−δ−α/β, β = 1. For Cases 1 to 4, α = −0.7707 and δ = −0.7628. The field is defined for x within 0 to 2 and discretized using Δx = 0.02.For Cases 5 to 8, α(xj) and δ(xj) could be obtained by solving μlxj=δxj+expαxj+β2/2, σl2xj=expβ2−1·exp2αxj+β2. The mean μlxj=0 and the variance σl2xj is given by the autocovariance function. The field is defined for x within 0 to 1 and discretized using Δx = 0.012Cρ2xj,xk=exp−xj−xk23Cρ3xj,xk=exp−xj−xk×cos4πxj−xk4Cρ4xj,xk=exp−xj−xk2×cos4πxj−xk5Cρ5xj,xk=minxj,xk6Cρ6xj,xk=4minxj,xk−xjxk7Cρ7xj,xk=minxj,xk×cos4πxj−xk8Cρ8xj,xk=4minxj,xk−xjxk×cos4πxj−xk


### 4.2. Correlation from Sampled Fields Based on the IRDR Procedure

In this section, the numerical analysis results obtained by iteratively using the IRDR procedure to simulate the non-Gaussian field are presented. For the analysis, unless otherwise indicated, the tolerance valueε is set equal to 10−4, a maximum iteration is set equal to 20, and the field is represented at points with the separation Δx between two adjacent points equal to 0.1 for Cases 1 to 4 and equal to 0.05 for Cases 5 to 8. In all cases, the number of points representing the field equals 101. The value of x for the first point is shifted from 0 to 0.002 for Cases 5 to 8, and the value of x for the last point is shifted from one to 0.998 for Cases 6 and 8 to avoid zero correlation in the diagonal of the correlation matrix. Note that the selection of Δx does not affect the application of the considered two procedures; the selected values of Δx used for the considered examples are sufficiently small to represent the autocorrelation function. The adopted tolerance of ε = 10–4, which is defined as the relative difference in the Frobenius norm for the correlation matrix, for the convergence criterion, is based on the consideration that such a relative difference is deemed sufficiently small. The considered maximum iteration is set equal to 20 since, as will be seen, the lowest relative difference in the Frobenius norm for the correlation matrix during the application of the procedure is achieved in less than 20 iterations.The correlation matrixCρY is calculated by using the correlation function shown in Table 1 for each case. If CρY is not positive semidefinite, it is replaced by its nearest correlation matrix C˜ρY [11]. It was found that CρY is not positive semidefinite only for Cases 2 and 4. The value of the element of C˜ρY−CρY is shown in Figure 2 for Cases 2 and 4. The results presented in the figure indicate that the absolute value of each element in C˜ρY−CρY is less than 10−5.Figure 2
Value of the element ofC˜ρY−CρY at (xj, xk) for Cases 2 and 4.Further analysis results by varyingp from 5 to 100 indicate that the first value of p that the Cholesky decomposition to CρY breaks down equals 17 (i.e., Δx = 0.125) for Case 2 and 25 (i.e., Δx = 0.0833) for Case 4. This shows that the proposed modification to the IRDR procedure by considering the nearest correlation matrix is necessary.By considering a sample size equal to 50,000 (i.e.,m = 50,000), samples of the random field for each combination of cases and the marginal distributions listed in Table 1 are generated. For the sampling, in this section, the initial samples of Y are obtained based on the marginal distribution alone (i.e., Option 1 shown in Figure 1). Typical samples of the random fields by using the IRDR procedure are shown in Figure 3. The Pearson correlation in the non-Gaussian space is calculated by using 50,000 samples for each considered combination. This calculated correlation is denoted as CρYRDS. For easy reference, a list of symbols used to denote the correlation matrices obtained based on different assumptions and approaches is depicted in Table 2. The autocorrelation (AC) function calculated based on CρYRDS is compared with that calculated based on CρY or C˜ρY, as shown in Figure 4 for Cases 1 to 4, in Figure 5 for Cases 5 to 8 by considering the beta distribution, and in Figure 6 for Cases 5 to 8 by considering the lognormal distribution. Also shown in Figures 5 and 6 is the difference CρYRDS−CρY or CρYRDS−C˜ρY.Figure 3
Typical samples generated by using the IRDR procedure.Table 2
List of symbols used for the correlation matrix.
Correlation matrixDefinitionCρYPearson correlation obtained based on the prescribed correlation functionC˜ρYNearest correlation matrix toCρY, if CρY is not positive semidefiniteCρYRDSPearson correlation fromY that are sampled using the IRDR procedureCGRDSPearson correlation in Gaussian space using samples mapped fromY that is simulated using the IRDR procedureCρsYRDSSpearman correlation fromY that is sampled using the IRDR procedureCGUnderlying Gaussian correlation obtained from the prescribed correlation function and equation (5)C˜GNearest correlation matrix toCG, if CG is not positive semidefiniteCρYNTSPearson correlation fromY that are sampled using the NTS procedureCρsYNTSSpearman correlation fromY that is sampled using the NTS procedureFigure 4
Comparison of the autocorrelation function for Cases 1 to 4 (CρY for Cases 2 and 4 is replaced by C˜ρY).Figure 5
Comparison of the correlation matrix for Cases 5 to 8 by considering the beta distribution as the marginal distribution.Figure 6
Comparison of the correlation matrix for Cases 5 to 8 with the lognormal distribution as the marginal distribution function.The results presented in Figures4 and 6 indicate that, for the weakly non-Gaussian cases, CρYRDS is practically equal to its target and that the difference between CρYRDS and its target can be large for strongly non-Gaussian cases. The differences are most noticeable for Cases 3, 4, 7, and 8 if the marginal distribution is the shifted lognormal distribution.To investigate whether the Pearson and Spearman correlations of the sampled fields follow the same trend, we calculated the Spearman correlation from the sampled fieldsCρYRDS. For the cases when the lognormal distribution is used as the marginal distribution, the differences in the autocorrelation function calculated based on CρsYRDS−CρYRDS and CρsYRDS−CρY are presented in Figure 7 for Cases 1 to 4 and in Figure 8 for Cases 5 to 8. The results for the cases when the beta distribution is used as the marginal distribution are not plotted since these differences are small and negligible.Figure 7
Differences in the autocorrelation function calculated based onCρsYRDS−CρYRDS and CρsYRDS−CρY for Cases 1 to 4 (CρY for Cases 2 and 4 is replaced by C˜ρY).Figure 8
Differences in the elements of (a)CρsYRDS−CρYRDS and of (b) CρsYRDS−CρY for Cases 5 to 8.
(a)(b)Results presented in Figures7 and 8 indicate that the Spearman and Pearson correlations from the sampled fields follow the same trend for the considered cases, except for Case 4 with the lognormal distribution as the marginal distribution. This implies that, in general, the accuracy of matching the prescribed Pearson correlation is influenced by the Spearman correlation. It must be emphasized that whether CρYRDS provides an accurate representation of the actual Spearman correlation cannot be inferred since the actual Spearman correlation is not given.An observation on the use of the IRDR procedure for the considered cases deserves comment. It is observed that, for the cases where the convergence is achieved, it is achieved in less than five iterations, as shown in Figure9. The figure shows that the reduction of the error is most significant for the first reordering, implying that the RDR procedure [14] for the Latin hypercube sampling is equally effective for nonstratified sampling. Additional iterations could further reduce the error. However, the figure also shows that convergence is not always achieved, and the algorithm diverges for Cases 2 and 4 where the nearest correlation matrix is used. In such cases, Y corresponding to the minimum ε=AρY−CρY/CρY (CρY is replaced by C˜ρY for Cases 2 and 4) is used as the sampled fields; this adopted criterion is explained in the previous sections. After a few iterations, if ε remains unchanged it implies that there is no additional reordering. In all cases, the procedure is very efficient.Figure 9
Errorε versus iteration by applying the IRDR procedure. Option 1 is used to generate the initial samples.
(a)(b)(c)(d)


### 4.3. Effect of the Initial Samples on Simulating Random Field by Using the IRDR Procedure

In the previous section, the initial samples of the random fields are based on the marginal distribution alone (i.e., Option 1, see Figure1). The simulation of the random fields is repeated but uses Option 2 to generate initial samples (i.e., using correlated samples defined by the Gaussian copula). To see the impact of Option 2 on the simulated fields, the value of ε obtained in each iteration for all the considered cases is shown in Figure 10. The figure indicates that the error without considering reordering is already small, suggesting that the use of Gaussian copula with the correlation matrix in the Gaussian space equal to that in the non-Gaussian space already provides good approximations for the considered cases. The trends presented in Figure 10 are similar to those shown in Figure 9, except that the iteration number in Figure 10 is downshifted by one as compared to that depicted in Figure 9. The error for some cases decreases but increases for other cases due to reordering. If the error increases with an increased iteration, it indicates that the use of the reordering cannot further improve over the initial samples of the field that are generated based on the Gaussian copula with the correlation matrix equal to CρY (or C˜ρY).Figure 10
Errorε versus iteration by applying the IRDR procedure CρY is replaced by C˜ρY for Cases 2 and 4; Option 2 is used to generate initial samples.
(a)(b)(c)(d)As the obtainedCρYRDS by using Option 2 to generate initial samples for the weakly non-Gaussian cases are practically identical to those shown previously based on Option 1, only the results for the strongly non-Gaussian cases are presented in the following. The AC function based on CρYRDS and its target are presented in Figure 11 for Cases 1 to 4 and in Figure 12 for Cases 5 to 8. A comparison of the results presented in Figures 4 and 11 and in Figures 6 and 12 indicates that CρYRDS is consistent whether Option 1 or Option 2 is used to generate the initial samples. Therefore, it is recommended to use Option 2 to generate initial samples in the IRDR procedure for efficiency.Figure 11
Comparison of autocorrelation function calculated based onCρYRDS and its target (CρY for Cases 2 and 4 is replaced by C˜ρY). Option 2 is used to generate initial samples.Figure 12
Comparison of the elements of the correlation matrix for Cases 5 to 8 with the lognormal distribution as the marginal distribution function. Option 2 is used to generate initial samples.


### 4.4. Comparison of the Results from the IRDR and NTS Procedures

The numerical analysis by using the NTS procedure to simulate the non-Gaussian field for each combination of cases shown in Table1 is carried out and presented in this section. Again, the field is discretized at 101 points for the simulation. As the results for the weakly non-Gaussian fields are very close to the Gaussian fields and are consistent with those obtained based on the IRDR procedure, only the results for cases that are strongly non-Gaussian are given in the following.First, the underlying Gaussian correlation (or its corresponding nearest correlation) matrixCG is obtained by solving equation (5) for Cases 1 to 8 with the lognormal distribution as the marginal distributions listed in Table 1. Also, since CG for Cases 1 to 8 is not positive semidefinite, its corresponding nearest correlation matrix, C˜G, is calculated and used in the following.Simulation of the non-Gaussian random field, according to the NTS procedure, is carried out for a sample size of 50,000 for each of the considered cases. The Pearson correlationCρYNTS and the Spearman correlation CρsYNTS are calculated by using Y that are sampled based on the NTS procedure. A comparison of CρsYNTS−CρYNTS and CρsYNTS−CρY for Cases 1 to 4 is compared in Figure 13. The difference between CρsYNTS and CρYNTS and between CρsYNTS and CρY for Cases 5 to 8 is shown in Figure 14. A comparison of the results presented in Figures 7 and 13 indicates that while the results presented in Figure 7 for Cases 1 and 2 are mostly less than zero, those shown in Figure 13 are greater than zero. For Cases 3 and 4, although there are differences, the maximum magnitude of the differences is consistent. Also, a comparison of the results presented in Figures 8 and 14 indicates that the trends observed in Figure 14 are consistent with those that can be observed in Figure 8. These observations and the comparison shown in Figures 4 to 6 imply that CρYNTS is consistent with CρYRDR, suggesting that the use of the IRDR procedure and the NTS procedure, at least quantitatively, is equally adequate and also shares similar limitations in matching the prescribed correlation.Figure 13
Differences in the autocorrelation function calculated based onCρsYNTS−CρYNTS and CρsYNTS−CρY for Cases 1 to 4 with the lognormal distribution as the marginal distribution function (CρY for Cases 2 and 4 is replaced by C˜ρY).Figure 14
Differences in the elements of the correlation matrix: (a)CρsYNTS−CρYNTS and (b) CρsYNTS−CρY for Cases 5 to 8.
(a)(b)To quantify the differences in matching the target correlation matrix by using the IRDR procedure and NTS procedure, the smallest and largest values of the elements inCρYRDR−CρY and in CρYNTS−CρY are calculated and shown in Table 3 for all cases by considering the beta distribution or the lognormal distribution as the marginal distribution. The results presented in Table 3 show that the IRDR procedure outperforms the NTS procedure for the weakly non-Gaussian fields; this is reversed for the strongly non-Gaussian fields. Since there is no need to evaluate the underlying Gaussian correlation in the IRDR procedure, this procedure is much more efficient than the NTS procedure.Table 3
The smallest and largest values of the elements inCρYRDR−CρY and in CρYNTS−CρY.
Target correlationCaseBeta distribution as marginalLognormal distribution as marginalCρYRDR−CρYCρYNTS−CρYCρYRDR−CρYCρYNTS−CρYSmallestLargestSmallestLargestSmallestLargestSmallestLargestPositive10.0000.000−0.0450.022−0.0250.014−0.0430.0302−0.0010.001−0.0080.009−0.0630.061−0.0240.0283−0.0010.000−0.0180.019−0.1670.245−0.1320.3064−0.0490.036−0.0150.019−0.5750.547−0.1610.2135−0.0080.001−0.0170.013−0.0360.032−0.0520.02260.0000.000−0.0180.018−0.0310.008−0.0400.0407−0.0010.000−0.0130.017−0.0820.097−0.1250.2958−0.0010.000−0.0140.019−0.0490.086−0.0820.185Negative10.0000.0000.0000.0000.0000.0000.0000.00020.0000.0000.0000.0000.0000.0000.0000.00030.0000.003−0.0190.017−0.0120.474−0.0100.4494−0.0170.040−0.0140.034−0.0200.678−0.0380.58850.0000.0000.0000.0000.0000.0000.0000.00060.0000.0000.0000.0000.0000.0000.0000.00070.0000.006−0.0140.012−0.0210.524−0.0210.50280.0000.002−0.0180.018−0.0130.303−0.0240.284


## 4.1. Prescribed Conditions for the Random Fields

A set of eight correlation functions were considered by Kim and Shields [12] to test their procedure with the use of the KL expansion to simulate the random field. For each prescribed correlation function, they considered two marginal distributions: one is weakly non-Gaussian, and the other is strongly non-Gaussian. These cases are listed in Table 1. Cases 1 to 4 describe homogeneous non-Gaussian fields, while Cases 5 to 8 describe nonhomogeneous non-Gaussian fields. The consideration of the beta distribution with the parameters given in Table 1 results in weakly non-Gaussian fields. The consideration of the shifted lognormal distribution with the distribution parameters assigned in Table 1 leads to strongly non-Gaussian fields. A simple spectral analysis indicates that Cases 3, 4, 7, and 8 represent narrow-band fields. The use of KL expansion, which was considered in Kim and Shields [12] in connection with the NTS procedure and in Zheng et al. [18] and in Xiao and Hong [41] in connection with the IRDR procedure, is beyond the scope of the present study.Table 1
Definition of the cases (see [12]).
CasePrescribed autocorrelation function for homogeneous cases (1 to 4) or autocovariance function for nonhomogeneous cases (5 to 8)Distribution type and notesBeta distributionLognormal distribution1Cρ1xj,xk=exp−xj−xkFy=ΓpΓq/Γp+q∫0uzp−11−zq−1dzp = 4, q = 2. For Cases 1 to 4, u=y+3.74/1.87+3.74. The field is defined for x within 0 to 2 and discretized using Δx = 0.02.For Cases 5 to 8, u=y−ymin/ymax−ymin, ymin=μbxj−σbxjpp+q+1/q, ymax=μbxj+σbxjqp+q+1/p. The mean μbxj=0 and the standard deviation σbxj is given by the autocovariance function. The field is defined for x within 0 to 1 and discretized using Δx = 0.01Fy=Flny−δ−α/β, β = 1. For Cases 1 to 4, α = −0.7707 and δ = −0.7628. The field is defined for x within 0 to 2 and discretized using Δx = 0.02.For Cases 5 to 8, α(xj) and δ(xj) could be obtained by solving μlxj=δxj+expαxj+β2/2, σl2xj=expβ2−1·exp2αxj+β2. The mean μlxj=0 and the variance σl2xj is given by the autocovariance function. The field is defined for x within 0 to 1 and discretized using Δx = 0.012Cρ2xj,xk=exp−xj−xk23Cρ3xj,xk=exp−xj−xk×cos4πxj−xk4Cρ4xj,xk=exp−xj−xk2×cos4πxj−xk5Cρ5xj,xk=minxj,xk6Cρ6xj,xk=4minxj,xk−xjxk7Cρ7xj,xk=minxj,xk×cos4πxj−xk8Cρ8xj,xk=4minxj,xk−xjxk×cos4πxj−xk


## 4.2. Correlation from Sampled Fields Based on the IRDR Procedure

In this section, the numerical analysis results obtained by iteratively using the IRDR procedure to simulate the non-Gaussian field are presented. For the analysis, unless otherwise indicated, the tolerance valueε is set equal to 10−4, a maximum iteration is set equal to 20, and the field is represented at points with the separation Δx between two adjacent points equal to 0.1 for Cases 1 to 4 and equal to 0.05 for Cases 5 to 8. In all cases, the number of points representing the field equals 101. The value of x for the first point is shifted from 0 to 0.002 for Cases 5 to 8, and the value of x for the last point is shifted from one to 0.998 for Cases 6 and 8 to avoid zero correlation in the diagonal of the correlation matrix. Note that the selection of Δx does not affect the application of the considered two procedures; the selected values of Δx used for the considered examples are sufficiently small to represent the autocorrelation function. The adopted tolerance of ε = 10–4, which is defined as the relative difference in the Frobenius norm for the correlation matrix, for the convergence criterion, is based on the consideration that such a relative difference is deemed sufficiently small. The considered maximum iteration is set equal to 20 since, as will be seen, the lowest relative difference in the Frobenius norm for the correlation matrix during the application of the procedure is achieved in less than 20 iterations.The correlation matrixCρY is calculated by using the correlation function shown in Table 1 for each case. If CρY is not positive semidefinite, it is replaced by its nearest correlation matrix C˜ρY [11]. It was found that CρY is not positive semidefinite only for Cases 2 and 4. The value of the element of C˜ρY−CρY is shown in Figure 2 for Cases 2 and 4. The results presented in the figure indicate that the absolute value of each element in C˜ρY−CρY is less than 10−5.Figure 2
Value of the element ofC˜ρY−CρY at (xj, xk) for Cases 2 and 4.Further analysis results by varyingp from 5 to 100 indicate that the first value of p that the Cholesky decomposition to CρY breaks down equals 17 (i.e., Δx = 0.125) for Case 2 and 25 (i.e., Δx = 0.0833) for Case 4. This shows that the proposed modification to the IRDR procedure by considering the nearest correlation matrix is necessary.By considering a sample size equal to 50,000 (i.e.,m = 50,000), samples of the random field for each combination of cases and the marginal distributions listed in Table 1 are generated. For the sampling, in this section, the initial samples of Y are obtained based on the marginal distribution alone (i.e., Option 1 shown in Figure 1). Typical samples of the random fields by using the IRDR procedure are shown in Figure 3. The Pearson correlation in the non-Gaussian space is calculated by using 50,000 samples for each considered combination. This calculated correlation is denoted as CρYRDS. For easy reference, a list of symbols used to denote the correlation matrices obtained based on different assumptions and approaches is depicted in Table 2. The autocorrelation (AC) function calculated based on CρYRDS is compared with that calculated based on CρY or C˜ρY, as shown in Figure 4 for Cases 1 to 4, in Figure 5 for Cases 5 to 8 by considering the beta distribution, and in Figure 6 for Cases 5 to 8 by considering the lognormal distribution. Also shown in Figures 5 and 6 is the difference CρYRDS−CρY or CρYRDS−C˜ρY.Figure 3
Typical samples generated by using the IRDR procedure.Table 2
List of symbols used for the correlation matrix.
Correlation matrixDefinitionCρYPearson correlation obtained based on the prescribed correlation functionC˜ρYNearest correlation matrix toCρY, if CρY is not positive semidefiniteCρYRDSPearson correlation fromY that are sampled using the IRDR procedureCGRDSPearson correlation in Gaussian space using samples mapped fromY that is simulated using the IRDR procedureCρsYRDSSpearman correlation fromY that is sampled using the IRDR procedureCGUnderlying Gaussian correlation obtained from the prescribed correlation function and equation (5)C˜GNearest correlation matrix toCG, if CG is not positive semidefiniteCρYNTSPearson correlation fromY that are sampled using the NTS procedureCρsYNTSSpearman correlation fromY that is sampled using the NTS procedureFigure 4
Comparison of the autocorrelation function for Cases 1 to 4 (CρY for Cases 2 and 4 is replaced by C˜ρY).Figure 5
Comparison of the correlation matrix for Cases 5 to 8 by considering the beta distribution as the marginal distribution.Figure 6
Comparison of the correlation matrix for Cases 5 to 8 with the lognormal distribution as the marginal distribution function.The results presented in Figures4 and 6 indicate that, for the weakly non-Gaussian cases, CρYRDS is practically equal to its target and that the difference between CρYRDS and its target can be large for strongly non-Gaussian cases. The differences are most noticeable for Cases 3, 4, 7, and 8 if the marginal distribution is the shifted lognormal distribution.To investigate whether the Pearson and Spearman correlations of the sampled fields follow the same trend, we calculated the Spearman correlation from the sampled fieldsCρYRDS. For the cases when the lognormal distribution is used as the marginal distribution, the differences in the autocorrelation function calculated based on CρsYRDS−CρYRDS and CρsYRDS−CρY are presented in Figure 7 for Cases 1 to 4 and in Figure 8 for Cases 5 to 8. The results for the cases when the beta distribution is used as the marginal distribution are not plotted since these differences are small and negligible.Figure 7
Differences in the autocorrelation function calculated based onCρsYRDS−CρYRDS and CρsYRDS−CρY for Cases 1 to 4 (CρY for Cases 2 and 4 is replaced by C˜ρY).Figure 8
Differences in the elements of (a)CρsYRDS−CρYRDS and of (b) CρsYRDS−CρY for Cases 5 to 8.
(a)(b)Results presented in Figures7 and 8 indicate that the Spearman and Pearson correlations from the sampled fields follow the same trend for the considered cases, except for Case 4 with the lognormal distribution as the marginal distribution. This implies that, in general, the accuracy of matching the prescribed Pearson correlation is influenced by the Spearman correlation. It must be emphasized that whether CρYRDS provides an accurate representation of the actual Spearman correlation cannot be inferred since the actual Spearman correlation is not given.An observation on the use of the IRDR procedure for the considered cases deserves comment. It is observed that, for the cases where the convergence is achieved, it is achieved in less than five iterations, as shown in Figure9. The figure shows that the reduction of the error is most significant for the first reordering, implying that the RDR procedure [14] for the Latin hypercube sampling is equally effective for nonstratified sampling. Additional iterations could further reduce the error. However, the figure also shows that convergence is not always achieved, and the algorithm diverges for Cases 2 and 4 where the nearest correlation matrix is used. In such cases, Y corresponding to the minimum ε=AρY−CρY/CρY (CρY is replaced by C˜ρY for Cases 2 and 4) is used as the sampled fields; this adopted criterion is explained in the previous sections. After a few iterations, if ε remains unchanged it implies that there is no additional reordering. In all cases, the procedure is very efficient.Figure 9
Errorε versus iteration by applying the IRDR procedure. Option 1 is used to generate the initial samples.
(a)(b)(c)(d)


## 4.3. Effect of the Initial Samples on Simulating Random Field by Using the IRDR Procedure

In the previous section, the initial samples of the random fields are based on the marginal distribution alone (i.e., Option 1, see Figure1). The simulation of the random fields is repeated but uses Option 2 to generate initial samples (i.e., using correlated samples defined by the Gaussian copula). To see the impact of Option 2 on the simulated fields, the value of ε obtained in each iteration for all the considered cases is shown in Figure 10. The figure indicates that the error without considering reordering is already small, suggesting that the use of Gaussian copula with the correlation matrix in the Gaussian space equal to that in the non-Gaussian space already provides good approximations for the considered cases. The trends presented in Figure 10 are similar to those shown in Figure 9, except that the iteration number in Figure 10 is downshifted by one as compared to that depicted in Figure 9. The error for some cases decreases but increases for other cases due to reordering. If the error increases with an increased iteration, it indicates that the use of the reordering cannot further improve over the initial samples of the field that are generated based on the Gaussian copula with the correlation matrix equal to CρY (or C˜ρY).Figure 10
Errorε versus iteration by applying the IRDR procedure CρY is replaced by C˜ρY for Cases 2 and 4; Option 2 is used to generate initial samples.
(a)(b)(c)(d)As the obtainedCρYRDS by using Option 2 to generate initial samples for the weakly non-Gaussian cases are practically identical to those shown previously based on Option 1, only the results for the strongly non-Gaussian cases are presented in the following. The AC function based on CρYRDS and its target are presented in Figure 11 for Cases 1 to 4 and in Figure 12 for Cases 5 to 8. A comparison of the results presented in Figures 4 and 11 and in Figures 6 and 12 indicates that CρYRDS is consistent whether Option 1 or Option 2 is used to generate the initial samples. Therefore, it is recommended to use Option 2 to generate initial samples in the IRDR procedure for efficiency.Figure 11
Comparison of autocorrelation function calculated based onCρYRDS and its target (CρY for Cases 2 and 4 is replaced by C˜ρY). Option 2 is used to generate initial samples.Figure 12
Comparison of the elements of the correlation matrix for Cases 5 to 8 with the lognormal distribution as the marginal distribution function. Option 2 is used to generate initial samples.


## 4.4. Comparison of the Results from the IRDR and NTS Procedures

The numerical analysis by using the NTS procedure to simulate the non-Gaussian field for each combination of cases shown in Table1 is carried out and presented in this section. Again, the field is discretized at 101 points for the simulation. As the results for the weakly non-Gaussian fields are very close to the Gaussian fields and are consistent with those obtained based on the IRDR procedure, only the results for cases that are strongly non-Gaussian are given in the following.First, the underlying Gaussian correlation (or its corresponding nearest correlation) matrixCG is obtained by solving equation (5) for Cases 1 to 8 with the lognormal distribution as the marginal distributions listed in Table 1. Also, since CG for Cases 1 to 8 is not positive semidefinite, its corresponding nearest correlation matrix, C˜G, is calculated and used in the following.Simulation of the non-Gaussian random field, according to the NTS procedure, is carried out for a sample size of 50,000 for each of the considered cases. The Pearson correlationCρYNTS and the Spearman correlation CρsYNTS are calculated by using Y that are sampled based on the NTS procedure. A comparison of CρsYNTS−CρYNTS and CρsYNTS−CρY for Cases 1 to 4 is compared in Figure 13. The difference between CρsYNTS and CρYNTS and between CρsYNTS and CρY for Cases 5 to 8 is shown in Figure 14. A comparison of the results presented in Figures 7 and 13 indicates that while the results presented in Figure 7 for Cases 1 and 2 are mostly less than zero, those shown in Figure 13 are greater than zero. For Cases 3 and 4, although there are differences, the maximum magnitude of the differences is consistent. Also, a comparison of the results presented in Figures 8 and 14 indicates that the trends observed in Figure 14 are consistent with those that can be observed in Figure 8. These observations and the comparison shown in Figures 4 to 6 imply that CρYNTS is consistent with CρYRDR, suggesting that the use of the IRDR procedure and the NTS procedure, at least quantitatively, is equally adequate and also shares similar limitations in matching the prescribed correlation.Figure 13
Differences in the autocorrelation function calculated based onCρsYNTS−CρYNTS and CρsYNTS−CρY for Cases 1 to 4 with the lognormal distribution as the marginal distribution function (CρY for Cases 2 and 4 is replaced by C˜ρY).Figure 14
Differences in the elements of the correlation matrix: (a)CρsYNTS−CρYNTS and (b) CρsYNTS−CρY for Cases 5 to 8.
(a)(b)To quantify the differences in matching the target correlation matrix by using the IRDR procedure and NTS procedure, the smallest and largest values of the elements inCρYRDR−CρY and in CρYNTS−CρY are calculated and shown in Table 3 for all cases by considering the beta distribution or the lognormal distribution as the marginal distribution. The results presented in Table 3 show that the IRDR procedure outperforms the NTS procedure for the weakly non-Gaussian fields; this is reversed for the strongly non-Gaussian fields. Since there is no need to evaluate the underlying Gaussian correlation in the IRDR procedure, this procedure is much more efficient than the NTS procedure.Table 3
The smallest and largest values of the elements inCρYRDR−CρY and in CρYNTS−CρY.
Target correlationCaseBeta distribution as marginalLognormal distribution as marginalCρYRDR−CρYCρYNTS−CρYCρYRDR−CρYCρYNTS−CρYSmallestLargestSmallestLargestSmallestLargestSmallestLargestPositive10.0000.000−0.0450.022−0.0250.014−0.0430.0302−0.0010.001−0.0080.009−0.0630.061−0.0240.0283−0.0010.000−0.0180.019−0.1670.245−0.1320.3064−0.0490.036−0.0150.019−0.5750.547−0.1610.2135−0.0080.001−0.0170.013−0.0360.032−0.0520.02260.0000.000−0.0180.018−0.0310.008−0.0400.0407−0.0010.000−0.0130.017−0.0820.097−0.1250.2958−0.0010.000−0.0140.019−0.0490.086−0.0820.185Negative10.0000.0000.0000.0000.0000.0000.0000.00020.0000.0000.0000.0000.0000.0000.0000.00030.0000.003−0.0190.017−0.0120.474−0.0100.4494−0.0170.040−0.0140.034−0.0200.678−0.0380.58850.0000.0000.0000.0000.0000.0000.0000.00060.0000.0000.0000.0000.0000.0000.0000.00070.0000.006−0.0140.012−0.0210.524−0.0210.50280.0000.002−0.0180.018−0.0130.303−0.0240.284


## 5. Conclusions

In the present study, we emphasize the underlying assumption for the iterative rank-dependent reordering (IRDR) procedure, indicating that while the objective of simulating the random field is to match the prescribed Pearson correlation, the iterative algorithm reorders the samples based on the Spearman correlation. This mismatching results in that the Pearson correlation of the samples of the random field differs from its target, especially if the random field is highly non-Gaussian. The differences are shown by numerical examples.It is pointed out that the Pearson correlation matrix for a given correlation function may not be positive semidefinite for practical applications. In such cases, we suggest replacing the target Pearson correlation matrix with its nearest correlation matrix. Since the IRDR procedure may not converge for a specified tolerance, we suggest incorporating correlated initial samples and taking a practical approach in assigning the sample fields (see Figure1). In short, we propose modifications to the IRDR procedure in terms of the initial samples, iteration involving decomposition, and convergence requirements to make the IRDR more efficient and robust.Results show that the limitation in matching the prescribed Pearson correlation by using the IRDR procedure to sample the nonhomogeneous and non-Gaussian fields is quantitatively similar to that of using the Nataf translation system. The numerical results for the considered examples also show that the IRDR procedure outperforms the translation process based procedure for the weakly non-Gaussian fields; this trend is reversed for the considered strongly non-Gaussian fields.



---

*Source: 2700498-2022-05-30.xml*