Unconditional and conditional simulation of nonstationary and non-Gaussian vector and field with prescribed marginal and correlation by using iteratively matched correlation

2.1. Estimating the underlying Gaussian correlation and unconditional simulation

Consider a one-dimensional random field $$ U(x) $$ with the marginal cumulative distribution function (CDF) represented by $$ F_{U(x);x}(u(x);x) $$ and the AC function, $$ \rho_Y\left(x_j,x_k\right) $$ , where $$ x_j $$ and $$ x_k $$ represent two points in the field. Its standardized version $$ Y(x) $$ , $$ Y(x)=(U(x)-\mu_{U(x)})/\sigma_{U(x)} $$ , has zero-mean and unit variance, and the marginal CDF, $$ F_{Y(x);x}(y(x);x) $$ equals $$ F_{U(x);x}(u(x);x) $$ . Its AC function is the same as that of $$ U(x), \rho_Y\left(x_j,x_k\right) $$ . Note that the distribution type $$ F_{Y(x);x}(y(x);x) $$ at two different points within the field may differ and $$ \rho_Y\left(x_j,x_k\right) $$ can be used to represent the AC function of a homogeneous/nonhomogeneous field. Once a sample of $$ Y(x), y(x) $$ , is obtained, its corresponding sample $$ u(x) $$ can be calculated using $$ u(x)=\mu_{U(x)}+\sigma_{U(x)}y(x) $$ .

Since, in many practical applications, the data collection and data analysis for a random field is carried out at discrete points with a given sampling interval. Let $$ [x]_p=[x_1,\cdots,x_p]^T $$ represent a vector of $$ x $$ of size $$ p\times1 $$ at the discrete points, where the superscript $$ T $$ denotes the transpose, and $$ x_j $$ represents the $$ j $$ -th element of the vector. It is considered that the AC function of the field at $$ [x]_p $$ can be adequately represented by the correlation coefficient matrix $$ \left[\kern-0.15em\left[ {C_{\rho Y}} \right]\kern-0.15em\right]$$ with the $$ (j,k) $$ -th element equal to $$ \rho_Y\left(x_j,x_k\right) $$ for $$ j, k = 1,\cdots, p $$ . Since $$ F_{Y(x_j);x_j}(y(x_j);x_j) $$ and $$ \left[\kern-0.15em\left[ {C_{\rho Y}} \right]\kern-0.15em\right]$$ can also be used to describe a vector of a non-Gaussian vector with $$ x $$ as a dummy variable (i.e., only the index is required), the proposed simulation procedure for the discretized field in the following is also applicable to a vector of non-Gaussian random variables.

If the marginal CDF is Gaussian and $$ \left[\kern-0.15em\left[ {C_{\rho Y}} \right]\kern-0.15em\right]$$ is positive semi-definite, the simulation can be carried out using the well-known procedure based on the Cholesky decomposition^[37],

where $$ [y]_p=[y(x_1),\cdots,y(x_p)]^T $$ represents a sample of $$ [Y]_p=[Y_1,\cdots,Y_p]^T $$ and $$ Y_j=Y(x_j) $$ , $$ \left[\kern-0.15em\left[ {L} \right]\kern-0.15em\right] $$ denotes the lower triangular matrix obtained by applying Cholesky decomposing to $$ \left[\kern-0.15em\left[ {C_{\rho Y}} \right]\kern-0.15em\right]$$ , $$ [z]_p=[z_1,\cdots,z_p]^T $$ is a sample of $$ [Z]_p $$ of size $$ p\times1 $$ , and $$ Z_j $$ and $$ Z_k $$ are standard independent normal variables for $$ j\neq k $$ .

However, if the prescribed marginal distribution is not Gaussian, $$ [y]_p $$ obtained from Equation (1) is inadequate since the samples from Equation (1) do not follow the prescribed marginal CDF. To overcome this, we note that given the underlying Gaussian correlation, $$ \rho_G $$ , the correlation coefficient $$ Y_j $$ and $$ Y_k $$ , according to the Nataf translation system, $$ c_{jk}(\rho_G) $$ , is given by^{[4, 5, 12, 13]},

where $$ \phi_2(\psi_j,\psi_k;\rho_G) $$ is the joint probability distribution function of two correlated standard normal variables with correlation $$ \rho_G $$ , and $$ \mathrm{\Phi}(\ ) $$ is the standard normal CDF. However, it must be noted that given $$ F_{Y(x);x}(y(x);x_j) $$ and $$ c_{jk}(\rho_G) $$ , the underlying Gaussian correlation $$ \rho_G $$ may not exist^{[4, 13, 16]}. As the mapping shown in Equation (2) is from the Gaussian to non-Gaussian random variables, the properties of $$ c_{jk}(\rho_G) $$ and $$ \rho_G $$ are as follows^{[13, 16, 38, 39]}:

a) $$ \left|c_{jk}(\rho_G)\right|\le\left|\rho_G\right| $$ , which also indicates that $$ c_{jk}(\rho_G)=0 $$ if $$ \rho_G=0 $$ ;

b) $$ c_{jk}(\rho_G)\geq0 (\le0) $$ if $$ \rho_G\geq0 (\le0) $$ . $$ c_{jk}(\rho_G) $$ is a non-decreasing function for $$ -1\le\rho_G\le1 $$ and is continuous under mild conditions. $$ c_{jk}(\rho_G) $$ is a strictly increasing function of $$ \rho_G $$ for continuous random variables with finite variance, and

c) The maximum and minimum values of feasible correlation coefficients $$ \hat{\rho}_{N G, j k} $$ and $$ \breve{\rho}_{N G, j k} $$ for $$ Y_j $$ and $$ Y_k $$ are equal to $$ c_{jk}(1) $$ and $$ c_{jk}(-1) $$ , respectively;

$$ c_{jk}(\pm1) $$ (i.e., $$ \hat{\rho}_{N G, j k} $$ and $$ {\breve{\rho}}_{NG,jk}) $$ can be calculated based on $$ m $$ samples of $$ [Z]_p $$ , which is included in a $$ p\times m $$ matrix $$ \left[\kern-0.15em\left[ {z} \right]\kern-0.15em\right] $$ with its $$ j $$ -th column denoted as $$ [z]_{p,j}, j = 1,\cdots, m $$ . Using these samples and Equation (2) results in,

Once $$ \hat{\rho}_{N G, j k} $$ and $$ {\breve{\rho}}_{NG,jk} $$ are calculated, they can be used to modify $$ \left[\kern-0.15em\left[ {C_{\rho Y}} \right]\kern-0.15em\right]$$ to satisfy the maximum and minimum requirements. This results in the modified symmetric matrix $$ \left[\kern-0.15em\left[ {C_{\rho Y}^\ast} \right]\kern-0.15em\right] $$ with its diagonal equal to one and its $$ (j,k) $$ -th element $$ \rho_Y^\ast\left(x_j,x_k\right) $$ , for $$ j\neq k $$ , is given by,

For those elements in $$ \left[\kern-0.15em\left[ {C_{\rho Y}} \right]\kern-0.15em\right]$$ whose elements are greater than $$ \hat{\rho}_{N G, j k} $$ or smaller than $$ {\breve{\rho}}_{NG,jk} $$ , the underlying Gaussian correlation coefficients equal $$ 1 $$ and $$ -1 $$ , respectively. For the elements in $$ \left[\kern-0.15em\left[ {C_{\rho Y} ^\ast} \right]\kern-0.15em\right] $$ that are not equal to $$ \hat{\rho}_{N G, j k} $$ or $$ {\breve{\rho}}_{NG,jk} $$ , the underlying Gaussian correlation coefficient $$ \rho_{G,jk} $$ for the non-Gaussian random variables $$ Y_j $$ and $$ Y_k $$ , for $$ j < k $$ , with prescribed $$ \rho_Y^\ast\left(x_j,x_k\right) $$ can be evaluated by solving $$ c_{jk}(\rho_{G,jk})-\rho_Y^\ast\left(x_j,x_k\right)=0 $$ , where $$ c_{jk}(\rho_{G,jk}) $$ is evaluated by using samples and Equation (2)

We note that $$ \rho_{G,jk} $$ is equal to $$ 1 $$ and $$ -1 $$ for cases where $$ \rho_Y^\ast\left(x_j,x_k\right) $$ equals $$ \hat{\rho}_{N G, j k} $$ and $$ {\breve{\rho}}_{NG,jk} $$ , respectively. For other cases, an iterative numerical root-finding method such as the bisection method and the Newton-Raphson method can be used since $$ \rho_{G,jk} $$ is bounded by $$ \rho_Y^\ast(x_j,x_k)\le\rho_{G,jk}\le1 $$ for $$ \rho_Y^\ast(x_j,x_k)>0 $$ and by $$ -1\le\rho_{G,jk}\le\rho_Y^\ast(x_j,x_k) $$ for $$ \rho_Y^\ast(x_j,x_k)<0 $$ . For the presented study, we implemented the bisection method to find $$ \rho_{G,jk} $$ . We also formulate the problem of finding $$ \rho_{G,jk} $$ as an optimization problem,

A verification that the underlying Gaussian correlation matrix $$ \left[\kern-0.15em\left[ {C_{G}} \right]\kern-0.15em\right] $$ with the off-diagonal element equal to $$ \rho_{G,jk} $$ is a positive semi-definite matrix is to be carried out. If this condition is not satisfied, an algorithm for computing its corresponding nearest correlation matrix $$ \left[\kern-0.15em\left[ {\tilde{C}_{G}} \right]\kern-0.15em\right] $$ ^{[17, 18]} is employed to obtain a proper correlation coefficient matrix. Once the underlying Gaussian correlation matrix $$ \left[\kern-0.15em\left[ {C_{G}} \right]\kern-0.15em\right] $$ (or $$ \left[\kern-0.15em\left[ {\tilde{C}_{G}} \right]\kern-0.15em\right] $$ ) is available, the Gaussian field can be simulated, and the sampled field can be mapped to the non-Gaussian field through the probability distribution transformation. If $$ \left[\kern-0.15em\left[ {\tilde{C}_{G}} \right]\kern-0.15em\right] $$ is used for simulation, the corresponding correlation matrix $$ \left[\kern-0.15em\left[ {\tilde{C}_{\rho Y} ^\ast} \right]\kern-0.15em\right] $$ for the non-Gaussian field can be evaluated using the samples in the non-Gaussian field. Note that $$ \left[\kern-0.15em\left[ {\tilde{C}_{\rho Y} ^\ast} \right]\kern-0.15em\right] $$ may differ from $$ \left[\kern-0.15em\left[ {C_{\rho Y}} \right]\kern-0.15em\right]$$ because of the side effect of using the translation process and the nearest correlation matrix.

Based on the above description, the proposed simulation procedure in the present study for the non-Gaussian vector or field is summarized as follows:

1. Sample the $$ p\times m $$ matrix $$ \left[\kern-0.15em\left[ {Z} \right]\kern-0.15em\right] $$ with its $$ j $$ -th column $$ [z]_{p,j} $$ denoting the $$ j $$ -th sample of $$ [Z]_p $$ ;

2. Evaluate $$ c_{jk}(\pm1) $$ (i.e., $$ \hat{\rho}_{N G, j k} $$ and $$ {\breve{\rho}}_{NG,jk} $$ ) using Equation (3), and calculate $$ \rho_Y^\ast\left(x_j,x_k\right) $$ using Equation (4);

3. Calculate the underlying Gaussian correlation coefficients $$ \rho_{G,jk} $$ by solving $$ c_{jk}(\rho_{G,jk})-\rho_Y^\ast\left(x_j,x_k\right)=0 $$ for $$ j = 1,\cdots,p, $$ and $$ k = j+1,\cdots, p $$ , noting that $$ \rho_{G,jk} $$ equals $$ 1 $$ and $$ -1 $$ for $$ \rho_Y^\ast\left(x_j,x_k\right) $$ equal to $$ \hat{\rho}_{N G, j k} $$ and $$ {\breve{\rho}}_{NG,jk} $$ , respectively;

4. Form the underlying Gaussian correlation matrix $$ \left[\kern-0.15em\left[ {C_{G}} \right]\kern-0.15em\right] $$ and verify it is positive semi-definite. If this condition is not satisfied, compute the nearest correlation matrix $$ \left[\kern-0.15em\left[ {\tilde{C}_{G}} \right]\kern-0.15em\right] $$ ^[18].

5. Apply Cholesky decomposition^[40] to $$ \left[\kern-0.15em\left[ {C_{G}} \right]\kern-0.15em\right] $$ (or $$ \left[\kern-0.15em\left[ {\tilde{C}_{G}} \right]\kern-0.15em\right] $$ if $$ \left[\kern-0.15em\left[ {C_{G}} \right]\kern-0.15em\right] $$ is not positive semi-definite) to obtain the lower diagonal matrix $$ \left[\kern-0.15em\left[ {L} \right]\kern-0.15em\right] $$ , and calculate $$ \left[\kern-0.15em\left[ {\Psi} \right]\kern-0.15em\right]=\left[\kern-0.15em\left[ {L} \right]\kern-0.15em\right] \left[\kern-0.15em\left[ {z} \right]\kern-0.15em\right] $$ . Apply $$ F_{Y(x_j);x_j}^{-1}(\mathrm{\Phi}(\psi_{j,q})) $$ to each element of $$ \left[\kern-0.15em\left[ {\Psi} \right]\kern-0.15em\right] $$ to obtain $$ \left[\kern-0.15em\left[ {y} \right]\kern-0.15em\right] $$ with its $$ j $$ -th column representing the $$ j $$ -th sample of the field. The samples $$ \left[\kern-0.15em\left[ {y} \right]\kern-0.15em\right] $$ can be used to obtain $$ u(x) $$ based on $$ u(x)=\mu_{U(x)}+\sigma_{U(x)}y(x) $$ . If $$ \left[\kern-0.15em\left[ {\tilde{C}_{G}} \right]\kern-0.15em\right] $$ is used for simulation, $$ \left[\kern-0.15em\left[ {y} \right]\kern-0.15em\right] $$ can be used to evaluate $$ \left[\kern-0.15em\left[ {\tilde{C}_{\rho Y} ^\ast} \right]\kern-0.15em\right] $$ for comparison purposes.

A flowchart showing the above steps are presented in Figure 1, where the panel for conditional simulation is to be discussed in the following sections. The procedure can be applied to the homogeneous/nonhomogeneous non-Gaussian process as well as a vector of the non-Gaussian random variables, as mentioned earlier.

Figure 1. Flowchart of the proposed simulation procedure: the left panel represents the unconditional simulation, and the right panel shows the steps for the conditional simulation of the field.

2.2. Sampling of the simulated field by using the KL expansion

In the previous section, it is considered that the marginal CDF and the correlation function are prescribed for the nonhomogeneous non-Gaussian field. Moreover, $$ \left[\kern-0.15em\left[ {C_{\rho Y}} \right]\kern-0.15em\right]$$ (or $$ \left[\kern-0.15em\left[ {\tilde{C}_{\rho Y} ^\ast} \right]\kern-0.15em\right] $$ ) for the field at selected discrete points $$ [x]_p $$ is established, and the underlying Gaussian correlation matrix $$ \left[\kern-0.15em\left[ {C_{G}} \right]\kern-0.15em\right] $$ (or its corresponding nearest correlation matrix $$ \left[\kern-0.15em\left[ {\tilde{C}_{G}} \right]\kern-0.15em\right] $$ ) is obtained at $$ [x]_p $$ . For simplicity and without loss of generality, unless otherwise indicated, the notations $$ \left[\kern-0.15em\left[ {C_{G}} \right]\kern-0.15em\right] $$ and $$ \left[\kern-0.15em\left[ {C_{\rho Y}} \right]\kern-0.15em\right]$$ are replaced by $$ \left[\kern-0.15em\left[ {\tilde{C}_{G}} \right]\kern-0.15em\right] $$ and $$ \left[\kern-0.15em\left[ {\tilde{C}_{\rho Y} ^\ast} \right]\kern-0.15em\right] $$ in the following.

The use of $$ \left[\kern-0.15em\left[ {\tilde{C}_{\rho Y} ^\ast} \right]\kern-0.15em\right] $$ and KL expansion to simulate the non-Gaussian field by using an iterative procedure was presented in several studies^{[32, 33]}. Since $$ \left[\kern-0.15em\left[ {\tilde{C}_{G}} \right]\kern-0.15em\right] $$ is already available, it can be used with the KL expansion to simulate the Gaussian field and map the simulated field to the non-Gaussian field through probability distribution transformation^[35]. The application of the KL expansion has been considered in the literature^[41].

For completeness and in order to provide the background for the conditional simulation based on the KL expansion, the use of the KL expansion is summarized below. In general, the random field (Gaussian or non-Gaussian) $$ Y(x) $$ can be expressed in the KL expansion^{[19, 32]},

where $$ w_q $$ represents independent zero-mean and unit variance random variables. For practical applications, the eigenvalue $$ \lambda_q $$ , and the eigenfunction $$ f_q(x_{,j}) $$ can be determined by solving Fredholm integral equation of the second kind^[42] over the domain of the field $$ \omega $$ ,

By considering that the series shown in Equation (7) is truncated to the first few terms, an approximation to $$ y(x) $$ is obtained.

By taking advantage of the availability of $$ \left[\kern-0.15em\left[ {\tilde{C}_{G}} \right]\kern-0.15em\right] $$ , the practical implementation of simulating the Gaussian implementation by using the KL expansion can be written as^[43],

where $$ [\psi]_p $$ represents a sample of the Gaussian field at $$ [x]_p $$ with unit variance and correlation coefficient $$ \left[\kern-0.15em\left[ {\tilde{C}_{G}} \right]\kern-0.15em\right] $$ , the eigenvalue $$ \lambda_q $$ and eigenvector $$ [\theta]_{p,q} $$ are obtained by using the eigendecomposition or Cholesky decomposition and singular value decomposition to $$ \left[\kern-0.15em\left[ {\tilde{C}_{G}} \right]\kern-0.15em\right] $$ . For efficiency, Equation (9) can be truncated to the first $$ N $$ terms with $$ N \ll p $$ by retaining the first $$ N $$ terms with significant eigenvalues resulting in,

By applying $$ F_{Y(x_j);x_j}^{-1}(\mathrm{\Phi}(\psi_j)) $$ , a sample of the non-Gaussian random field $$ [y]_p $$ is obtained, where $$ \psi_j j = 1,\cdots, p $$ , denotes the $$ j $$ -th element of $$ [\psi]_p $$ . Note that the type of approximation given in Equation (9) forms the basis for sparse sampling or compressive sensing^[44].

2.3. Conditional simulation

Consider that the field to be simulated $$ Y(x) [i.e., U(x)] $$ is constrained by the observed values $$ y(x) $$ at $$ x=x_k $$ where $$ k $$ takes the index in $$ \mathrm{\Omega}_c=[k_1,k_2,\cdots,k_{nc}] $$ . To carry out the conditional simulation of the field $$ y(x) $$ at $$ x_j $$ for $$ j = 1,\cdots,p $$ but $$ j\notin\mathrm{\Omega}_c $$ , we note that the underlying Gaussian correlation matrix $$ \left[\kern-0.15em\left[ {\tilde{C}_{G}} \right]\kern-0.15em\right] $$ the field at $$ [x]_p $$ can be obtained based on the procedure shown in the flowchart depicted in Figure 1. Therefore, we could take advantage of the availability of $$ \left[\kern-0.15em\left[ {\tilde{C}_{G}} \right]\kern-0.15em\right] $$ to carry out conditional simulation for the nonhomogeneous non-Gaussian field. This could be done by first mapping the conditioning $$ y(x_k) $$ , for $$ k\in\mathrm{\Omega}_c $$ , to the Gaussian space using,

Let $$ [\psi_\alpha] $$ denote the vector of size $$ nc\times 1 $$ for all the mapped conditioning in the Gaussian field. Similarly, let $$ [\psi_\beta] $$ denote the vector of size $$ (p-nc)\times 1 $$ containing the values of the (zero-mean and unit variance) Gaussian field for the remaining positions in the field $$ x_j $$ for $$ j $$ within $$ 1 $$ to $$ p $$ but $$ j\notin\mathrm{\Omega}_c $$ . The unknown values of the vector $$ [\psi_\beta] $$ conditioned on $$ [\psi_\alpha] $$ are governed by the following joint Gaussian distribution^[37],

where $$ N(\ ,\ ) $$ denotes the multi-dimensional normal distribution with $$ [\mu_{\beta \mid \alpha}] $$ and $$ \left[\kern-0.15em\left[ {\Sigma_{\beta \mid \alpha}} \right]\kern-0.15em\right] $$ represent the mean vector and covariance matrix, and C with subscript denotes submatrices that form the correlation coefficient matrix of the vector of $$ [[\psi_{\alpha}]^{T},[\psi_{\beta}]^{T}]^{T} $$ ,

This matrix can be constructed by rearranging $$ \left[\kern-0.15em\left[ {\tilde{C}_{G}} \right]\kern-0.15em\right] $$ of the field represented at discretized points of the field at $$ j = 1, \cdots, p. $$

By generating $$ m $$ conditional samples of $$ [\psi_\beta] $$ based on Equation (12) [see also Equation (1)], transforming each sampled value to $$ y(x) $$ by using $$ F_{Y(x);x}^{-1}(\mathrm{\Phi}(\psi)) $$ , and calculating $$ u(x)=\mu_{U(x)}+\sigma_{U(x)}y(x) $$ , samples of the non-Gaussian fields are obtained, which complements the conditioning value to describe the conditional random field. This conditional simulation algorithm for the nonhomogeneous non-Gaussian field is very simple to use [Figure 1] as the conditional Gaussian simulation by using Equation (12) is straightforward.

We also note that a method^[34] was developed to carry out the conditional simulation of the Gaussian field based on the KL expansion. It can be shown that their method for the underlying Gaussian random field at $$ p $$ discretized points $$ [x]_p, [\tilde{\psi}]_p $$ , [i.e., analog to Equation (7) vs. Equation (9)] for the given $$ \left[\kern-0.15em\left[ {\tilde{C}_{G}} \right]\kern-0.15em\right] $$ can be written as,

where the eigenvalue $$ \lambda_q $$ and eigenvector $$ [h]_{p,q} $$ are obtained by applying the eigendecomposition. In Equation (14), $$ {\tilde{\mu}}_q $$ denotes the $$ q $$ -th element of the conditional mean $$ [\tilde{\mu}]q $$ of the size $$ p\times 1 $$ that is given by,

and $$ {\tilde{S}}_{q,j} $$ denotes the $$ (q,j) $$ -the element of the matrix $$ \left[\kern-0.15em\left[ { \tilde{S} } \right]\kern-0.15em\right] $$ of the size $$ p\times p $$ that is given by,

where $$ \left[\kern-0.15em\left[ { I } \right]\kern-0.15em\right] $$ is the $$ p\times p $$ identity matrix, $$ \left[\kern-0.15em\left[ { G } \right]\kern-0.15em\right] $$ is an $$ nc\times p $$ submatrix of $$ \left[\kern-0.15em\left[ { h } \right]\kern-0.15em\right] $$ , in which the $$ j $$ -th column of $$ \left[\kern-0.15em\left[ { G } \right]\kern-0.15em\right] $$ is formed by the $$ k $$ -th elements of the $$ j $$ -th column of $$ \left[\kern-0.15em\left[ { h } \right]\kern-0.15em\right] $$ with the index $$ k\in\mathrm{\Omega}_c $$ .

Once $$ \tilde{\psi}(x) $$ or $$ [\tilde{\psi}]_p $$ is sampled, the conditional field $$ y(x) $$ or $$ [y]_p $$ is obtained by using $$ F_{Y(x_j);x_j}^{-1}(\mathrm{\Phi}({\tilde{\psi}}_j)) $$ . The flowchart for the conditional simulation is included in Figure 1, whether the conditional Gaussian model shown in Equation (12) or by using the KL expansion model shown in Equation (14) is considered. The use of Equation (12) requires the simulation of the correlated Gaussian random variables, which can be done by applying the Cholesky decomposition [Equation (1)]. The use of Equation (14) requires the evaluation of eigenvalues and eigenvectors if they are unavailable. In both cases, the use of probability transformation is required, and the steps are straightforward, as shown in Figure 1.

2.4. Discussion on the extension to the multi-dimensional random field

In this section, we highlight steps to extend the procedure described in the previous section to the $$ n $$ -dimensional case. Consider the case of simulating the $$ n $$ -dimensional zero-mean and unit variance random field that is defined over the domain of x, with the value represented by $$ y $$ (x), where $$ x=[x_1,\cdots,x_n] $$ . The field is defined by the prescribed marginal CDF, $$ F_{y(x);x}(y(x);x) $$ , and the correlation coefficient function. Similar to the one-dimensional field case, it is considered that the field is to be simulated at a set of discrete points that are denoted as $$ x_N=[x_{1,j_1},\cdots,x_{n,j_n}] $$ , where $$ [j]=[j_1,\cdots,j_n] $$ , $$ j_k=1,\cdots,p_k, $$ $$ p_k $$ denotes the total number of points along the $$ x_k $$ -axis, and $$ N $$ represents the total number of points (i.e., combinations of $$ j_k,\; k = 1,\cdots, n $$ , which equals $$ \prod_{k=1}^{n}p_k) $$ within the $$ n $$ -dimensional field to be considered. The correlation coefficient matrix for the random field represented at the discretized points is to be calculated based on the prescribed correlation coefficient function. In other words, we treat a multi-dimensional field by vectorizing it into a long one-dimensional vector, and the procedure described in the previous sections can be directly applied. Perhaps, the manipulation of the $$ n $$ indices in [j] for the points in the $$ n $$ -dimensional space could be tedious. A systematic ordering of the index could facilitate the computer program task.

Since the extension to the multi-dimensional field is straightforward, no further consideration for the $$ n $$ -dimensional case is made.

3.1. Examples of using the proposed procedure to evaluate underlying Gaussian correlation and unconditional simulation

Several numerical examples are presented by considering the homogeneous/nonhomogeneous processes with weakly and strongly non-Gaussian processes to illustrate the applicability of the proposed procedures depicted in Figure 1.

The examples in this subsection were considered in the literature^{[32, 35]}. The assumed marginal distributions and the AC function for the homogeneous non-Gaussian cases (i.e., Cases 1 to 8) are listed in Table 1. The use of the beta distribution described in the table leads to a weakly non-Gaussian field, while the use of the lognormal distribution described in the table results in a strongly non-Gaussian field.

First, for each of the considered homogeneous cases, the application of the proposed procedure shown in Figure 1 is employed. For the analysis, samples of the independent standard normally distributed random variables are simulated for a sample size $$ m = 50,000 $$ . This sample size is used throughout the numerical analysis in the present study. The selection of $$ m = 50,000 $$ is based on a sensitivity analysis, indicating that the use of $$ m = 10,000 $$ to $$ 100,000 $$ is adequate for most considered cases. An increased $$ m $$ results in an increased computing time. The selection is also based on balancing the accuracy and computing time. These samples are used to evaluate the bounds $$ \hat{\rho}_{N G, j k} $$ and $$ {\breve{\rho}}_{NG,jk} $$ . Based on the bounds shown in Equations (3) and (4), $$ \left[\kern-0.15em\left[ {C_{\rho Y} ^\ast} \right]\kern-0.15em\right] $$ is obtained, and the corresponding underlying Gaussian matrix $$ \left[\kern-0.15em\left[ {C_{G}} \right]\kern-0.15em\right] $$ is determined by using the bisection method and verified using the optimization approach formulated in Equation (6). In all cases, the obtained $$ \rho_{G,jk} $$ values by both approaches are practically identical. Since the same conclusion is observed for the remaining numerical examples, only the bisection method is used to calculate $$ \rho_{G,jk} $$ for the remaining numerical examples. Note that the use of the bisection method leads to an error (i.e., the last bracketed interval) of less than approximately $$ 1/2^{n_I} $$ , where $$ n_I $$ is the number of iterations. Therefore, an iteration of 7, 10, 13, and 17 leads to the absolute value of the error less than about $$ 1 \times 10^{-2} $$ , $$ 1 \times 10^{-3} $$ , $$ 1 \times 10^{-4} $$ , and $$ 1 \times 10^{-5} $$ . $$ n_I $$ less than 16 is used for all the computations. Moreover, a verification that $$ \left[\kern-0.15em\left[ {C_{G}} \right]\kern-0.15em\right] $$ is positive semi-definite is carried out. $$ \left[\kern-0.15em\left[ {C_{G}} \right]\kern-0.15em\right] $$ is replaced by its corresponding nearest correlation matrix^[18]$$ \left[\kern-0.15em\left[ {\tilde{C}_{G}} \right]\kern-0.15em\right] $$ if such a condition is not satisfied. Finally, the non-Gaussian field is simulated following the steps shown in Figure 1, and the AC function of simulated samples is calculated. The target AC function and the AC functions based on $$ \left[\kern-0.15em\left[ {C_{\rho Y} ^\ast} \right]\kern-0.15em\right] $$ and $$ \left[\kern-0.15em\left[ {C_{G}} \right]\kern-0.15em\right] $$ are presented in Figure 2 and Figure 3. For cases, when $$ \left[\kern-0.15em\left[ {\tilde{C}_{G}} \right]\kern-0.15em\right] $$ and $$ \left[\kern-0.15em\left[ {\tilde{C}_{\rho Y} ^\ast} \right]\kern-0.15em\right] $$ are calculated, they are shown in Figure 2 and Figure 3.

Figure 2. Comparison of the prescribed AC function, estimated AC function based on samples, and underlying Gaussian AC function for Cases 1 to 8 by considering that the marginal distribution is the beta distribution shown in Table 1. The plots for Cases 5 to 8, considered $$ x_k = 0.5 $$ . AC: Autocorrelation.

Figure 3. Comparison of the prescribed AC function, estimated AC function based on samples, and their differences for Cases 5 to 8. The considered marginal distribution is the beta distribution shown in Table 1. AC: Autocorrelation.

The results presented in Figure 2 and Figure 3 indicate that the underlying Gaussian AC functions are close to the target AC functions. The closeness is because the use of the beta distribution presented in Table 1 as the marginal distribution leads to only weakly non-Gaussian fields for the considered cases. The weakly non-Gaussian behavior was noted in the literature^[35]. For Cases 2 and 4 (for the considered grid system), it is observed that $$ \left[\kern-0.15em\left[ {C_{G}} \right]\kern-0.15em\right] $$ cannot be decomposed by the Cholesky decomposition. Therefore, its nearest correlation matrix $$ \left[\kern-0.15em\left[ {\tilde{C}_{G}} \right]\kern-0.15em\right] $$ is calculated and used as the basis to simulate the non-Gaussian field. The corresponding $$ \left[\kern-0.15em\left[ {\tilde{C}_{\rho Y} ^\ast} \right]\kern-0.15em\right] $$ is calculated from the simulated non-Gaussian field. The plots shown in Figure 2 indicate that the differences between $$ \left[\kern-0.15em\left[ {C_{G}} \right]\kern-0.15em\right] $$ and $$ \left[\kern-0.15em\left[ {\tilde{C}_{G}} \right]\kern-0.15em\right] $$ and between $$ \left[\kern-0.15em\left[ {C_{\rho Y}} \right]\kern-0.15em\right]$$ and $$ \left[\kern-0.15em\left[ {\tilde{C}_{\rho Y} ^\ast} \right]\kern-0.15em\right] $$ are very small. A further inspection of the numerical results indicates that the absolute value of the difference between the elements in $$ \left[\kern-0.15em\left[ {C_{G}} \right]\kern-0.15em\right] $$ and its corresponding element in $$ \left[\kern-0.15em\left[ {\tilde{C}_{G}} \right]\kern-0.15em\right] $$ is much less than 0.01.

The plots for Cases 5 to 8 shown in Figure 2 are those calculated for $$ x_j = 0.5 $$ and by varying $$ x_k $$ within 0 to 1. Since Cases 5 to 8 are nonhomogeneous, the surfaces for their prescribed AC functions and for $$ \left[\kern-0.15em\left[ {\tilde{C}_{\rho Y} ^\ast} \right]\kern-0.15em\right] $$ are presented in Figure 3, indicating that the latter mimics the former well. In addition, a quantitative comparison of the difference between the prescribed AC function and that based on $$ \left[\kern-0.15em\left[ {\tilde{C}_{\rho Y} ^\ast} \right]\kern-0.15em\right] $$ is presented in Figure 3. The differences are caused by the use of the limited samples (i.e., a sample size of 50000) as well as the use of the nearest correlation matrix. In all cases, the small differences presented in the figure confirm the adequacy of the procedure depicted in Figure 1.

By repeating the analysis that is carried out for the results presented in Figure 2 and Figure 3 but replacing the beta distribution with the lognormal distribution shown in Table 1, the obtained results are shown in Figure 4 and Figure 5. As pointed out in the study^[35], the use of this lognormal distribution results strongly non-Gaussian field. This strongly non-Gaussian behavior is reflected in the calculated underlying Gaussian AC functions.

Figure 4. Comparison of the prescribed AC function, estimated AC function based on samples, and underlying Gaussian AC function for Cases 1 to 8 by considering that the marginal distribution is the lognormal distribution shown in Table 1. The plots for Cases 5 to 8, considered $$ x_k = 0.5 $$ . AC: Autocorrelation.

Figure 5. Comparison of the prescribed AC function, estimated AC function based on samples, and their differences for Cases 5 to 8. The considered marginal distribution is the lognormal distribution shown in Table 1. AC: Autocorrelation.

The plots presented in Figure 4 for Cases 3, 4, 7, and 8 indicate that there are segments where $$ \rho_{G,jk} $$ equals $$ -1 $$ . This is because the calculated $$ \rho_Y(x_j,x_k) $$ is less than $$ {\breve{\rho}}_{NG,jk} $$ [see Equation (4)]. For these cases, $$ \left[\kern-0.15em\left[ {C_{\rho Y} ^\ast} \right]\kern-0.15em\right] $$ and $$ \left[\kern-0.15em\left[ {C_{G}} \right]\kern-0.15em\right] $$ are calculated. Since an assessment indicates that $$ \left[\kern-0.15em\left[ {C_{G}} \right]\kern-0.15em\right] $$ does not form a proper correlation coefficient matrix, its corresponding nearest correlation matrix $$ \left[\kern-0.15em\left[ {\tilde{C}_{G}} \right]\kern-0.15em\right] $$ is then calculated by using the algorithm given in the study^[18]. The non-Gaussian field is then simulated based on $$ \left[\kern-0.15em\left[ {\tilde{C}_{G}} \right]\kern-0.15em\right] $$ using the procedure described in Figure 1 and $$ \left[\kern-0.15em\left[ {\tilde{C}_{\rho Y} ^\ast} \right]\kern-0.15em\right] $$ is evaluated based on the simulated fields.

The plots presented in Figure 4 for Cases 3, 4, 7, and 8 show that $$ \left[\kern-0.15em\left[ {\tilde{C}_{\rho Y} ^\ast} \right]\kern-0.15em\right] $$ can deviate from their target. The differences illustrate the well-known limitation of using the Nataf translation system to represent the non-Gaussian field^{[12-14, 16, 35]} rather than the numerical procedure used in the present study.

Figure 5 shows the surface of the target AC functions and the surface of the underlying Gaussian AC functions. In the same figure, the differences between the surface of the target AC function and the surface of $$ \left[\kern-0.15em\left[ {\tilde{C}_{\rho Y} ^\ast} \right]\kern-0.15em\right] $$ are also depicted. The large differences between the prescribed AC function and that based on $$ \left[\kern-0.15em\left[ {\tilde{C}_{\rho Y} ^\ast} \right]\kern-0.15em\right] $$ for some regions of the field can be explained based on the arguments used to explain the results presented in Figure 4.

3.2. Examples using the KL expansion to simulate non-Gaussian field

In this subsection, we consider again all cases shown in Table 1 but use the KL expansion to augment the simulated field. For the analysis, we take advantage of the availability of $$ \left[\kern-0.15em\left[ {C_{G}} \right]\kern-0.15em\right] $$ and $$ \left[\kern-0.15em\left[ {C_{\rho Y}} \right]\kern-0.15em\right]$$ (if they are not proper correlation matrices, we use their corresponding nearest correlation matrices $$ \left[\kern-0.15em\left[ {\tilde{C}_{G}} \right]\kern-0.15em\right] $$ and $$ \left[\kern-0.15em\left[ {\tilde{C}_{\rho Y} ^\ast} \right]\kern-0.15em\right] $$ ), that is discussed in the previous section.

First, the eigendecomposition is applied to $$ \left[\kern-0.15em\left[ {\tilde{C}_{G}} \right]\kern-0.15em\right] $$ (or $$ \left[\kern-0.15em\left[ {C_{G}} \right]\kern-0.15em\right] $$ if it is a proper correlation matrix) for each of Cases 1 to 8. By sorting the eigenvalues in descending order, we sum up the eigenvalues of the first $$ j $$ terms and normalize them with respect to the sum of all the eigenvalues. The normalized value, denoted as $$ R_{j/T} $$ is plotted in Figure 6 for each case, indicating that $$ R_{j/T} $$ is greater than about 0.95. This indicates that the consideration of the first 21 terms of the KL expansion [see Equation (10)] could provide adequate representations of the considered random fields.

Figure 6. Value of $$ R_{j/T} $$ for the considered cases showing the contribution of the eigenvalues from the first $$ j $$ terms.

We use the obtained $$ \left[\kern-0.15em\left[ {\tilde{C}_{G}} \right]\kern-0.15em\right] $$ , considering $$ N = 21 $$ , applying Equation (10) and employing the probability transformation mapping to simulate the non-Gaussian random field at points with the separation between two adjacent points $$ \delta x = 0.1 $$ for Cases 1 to 4, and with $$ \delta x = 0.05 $$ for Cases 5 to 8. An illustration of five samples for each case is depicted in Figure 7. By generating 50, 000 samples for each case, the calculated correlation coefficient is presented in Figure 8 and compared with $$ \left[\kern-0.15em\left[ {\tilde{C}_{\rho Y} ^\ast} \right]\kern-0.15em\right] $$ that is obtained in the previous section (i.e., based on the procedure shown on the left panel in Figure 1). The comparison indicates that the correlation obtained from the samples agrees well with $$ \left[\kern-0.15em\left[ {\tilde{C}_{\rho Y} ^\ast} \right]\kern-0.15em\right] $$ .

Figure 7. Typical samples obtained by applying the KL expansions and probability transformation mapping for each of the cases listed in Table 1. The first two rows are obtained by considering the beta distribution presented in Table 1 and the last two rows are obtained by considering the lognormal distribution presented in Table 1. KL: Karhunen-Loeve.

Figure 8. Comparison of the calculated correlation coefficient from the sampled random field to $$ \left[\kern-0.15em\left[ { \tilde{C}_{\rho Y}^{*} } \right]\kern-0.15em\right] $$ obtained from the sampled random field based on the proposed procedure depicted on the left panel shown in Figure 1. The first two rows are obtained by considering the beta distribution presented in Table 1 and the last two rows are obtained by considering the lognormal distribution presented in Table 1.

3.3. Examples of conditional simulation

For illustrating the conditional simulation, we consider Cases 1 and 5. By assuming that the marginal distribution is the beta or lognormal distribution shown in Table 1, we sample a field at $$ p = 101 $$ points for the considered case. We select 5 points from the sampled field as the conditioning points that are shown in Figure 9. We then apply the procedure described in Figure 1 to sample the random fields conditionally. Five samples for each considered case are presented in Figure 9 by using Equation (12). As shown in the figure, the sampled random fields are constrained by the conditioning points.

Figure 9. Five conditionally sampled fields for Cases 1 and 5. The large dots represent the conditioning points and the lines represent the sampled fields.

Acknowledgments

We gratefully acknowledge the financial support received from the Natural Sciences and Engineering Research Council of Canada (RGPIN-2016-04814, for HPH), the China Scholarship Council (No. 201807980004 for MYX), and the University of Western Ontario.

Authors' contributions

Investigation, Methodology, Formal analysis, Writing - original draft, Writing - review & editing: Xiao MY

Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Writing - original draft, Writing - review & editing: Hong H

Availability of data and materials

Some or all data, models, or code generated or used during the study are available from the author by request.

Financial support and sponsorship

The authors acknowledge the financial support received from the Natural Sciences and Engineering Research Council of Canada (RGPIN-2016-04814, for HPH) and the China Scholarship Council (No.201807980004, for MYX).

Conflicts of interest

All authors declared that there are no conflicts of interest.

Ethical approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Copyright

© The Author(s) 2022.