Computational and Theoretical Details of Spatial Simulation

The SIM2D Procedure

Computational and Theoretical Details of Spatial Simulation

Introduction

There are a number of approaches to simulating spatial random fields or, more generally, simulating sets of dependent random variables. This includes sequential indicator methods, turning bands, and the Karhunen-Loeve Expansion. Refer to Christakos (1992, Chapter 8) and Duetsch and Journel (1992, Chapter V) for details.

A particularly simple method available for Gaussian spatial random fields is the LU decomposition method. This method is computationally efficient. For a given covariance matrix, the LU = LL^T decomposition is computed once, and the simulation proceeds by repeatedly generating a vector of independent N(0,1) random variables and multiplying by the L matrix.

One problem with this technique is memory requirements; memory is required to hold the full data and grid covariance matrix in core. While this is especially limiting in the three-dimensional case, you can use PROC SIM2D, which handles only two-dimensional data, for moderately sized simulation problems.

Theoretical Development

It is a simple matter to produce an N(0,1) random number, and by stacking k N(0,1) random numbers in a column vector, you can obtain a vector with independent standard normal components $W \sim N_k(0,I)$ . The meaning of the terms independence and randomness in the context of a deterministic algorithm required for the generation of these numbers is a little subtle; refer to Knuth (1981, Vol. 2, Chapter 3) for details.

Rather than $W \sim N_k(0,I)$ , what is required is the generation of a vector ${Z \sim N_k(0, C)}$ , that is,

$Z = [ Z_1\Z_2\\vdots \Z_k]$

with covariance matrix

$C = ( C_{11} & C_{12} & ... & C_{1k}\C_{21} & C_{22} & ... & C_{2k}\& \ddots & & \C_{k1} & C_{k2} & ... & C_{kk})$

If the covariance matrix is symmetric and positive definite, it has a Cholesky root L such that C can be factored as

C = LL^T

where L is lower triangular. Refer to Ralston and Rabinowitz (1978, Chapter 9, Section 3-3) for details. This vector Z can be generated by the transformation Z = LW. Note that this is where the assumption of a Gaussian SRF is crucial. When $W \sim N_k(0,I)$ , then Z = LW is also Gaussian. The mean of Z is

E(Z) = L(E(W)) = 0

and the variance is

Var(Z) = Var(LW) = E(LWW^TL^T) = LE(WW^T)L^T = LL^T = C

Consider now an SRF $Z(s), s \in D \subset \mathcal{R}^2$ , with spatial covariance function C(h). Fix locations s₁,s₂, ... ,s_k, and let Z denote the random vector

$Z = [ Z(s_1)\Z(s_2)\\vdots \Z(s_k)]$

with corresponding covariance matrix

$C_z = ( C(0) & C(s_1-s_2) & ... & C(s_1-s_k)\C(s_2-s_1) & C(0) & ... & C(s_2-s_k)\& \ddots & & \C(s_k-s_1) & C(s_k-s_2) & ... & C(0))$

Since this covariance matrix is symmetric and positive definite, it has a Cholesky root, and the Z(s_i), i = 1, ... ,k can be simulated as described previously. This is how the SIM2D procedure implements unconditional simulation in the zero mean case. More generally,

$Z(s) = \mu(s) + \varepsilon(s)$

with $\mu(s)$ being a quadratic form in the coordinates s = (x,y), and the $\varepsilon(s)$ being an SRF having the same covariance matrix C_z as previously. In this case, the $\mu(s_i), i=1, ... ,k$ is computed once and added to the simulated vector $\varepsilon(s_i), i=1, ... ,k$ for each realization.

For a conditional simulation, this distribution of

$Z = [ Z(s_1)\Z(s_2)\\vdots \Z(s_k)]$

must be conditioned on the observed data. The relevant general result concerning conditional distributions of multivariate normal random variables is the following. Let $X \sim N_m({\mu},{{\Sigma}})$ , where

$X = [ X_1\X_2\]$

${\mu}= [ {\mu}_1\{\mu}_2\]$

and

${{\Sigma}}= ( {{\Sigma}}_{11} & {{\Sigma}}_{12}\{{\Sigma}}_{21} & {{\Sigma}}_{22}\)$

The subvector X₁ is k×1, X₂ is n×1, ${{\Sigma}}_{11}$ is k×k, ${{\Sigma}}_{22}$ is n×n, and ${{\Sigma}}_{12}={{\Sigma}}_{21}^T$ is k×n, with k+n=m. The full vector X is partitioned into two subvectors X₁ and X₂, and ${{\Sigma}}$ is similarly partitioned into covariances and cross covariances.

With this notation, the distribution of X₁ conditioned on X₂ = x₂ is $N_k(\tilde{{\mu}},\tilde{{{\Sigma}}})$ ,with

$\tilde{{\mu}}={\mu}_1+{{\Sigma}}_{12}{{\Sigma}}_{22}^{-1}(x_2-{\mu}_2)$

and

$\tilde{{{\Sigma}}} = {{\Sigma}}_{11}-{{\Sigma}}_{12}{{\Sigma}}_{22}^{-1}{{\Sigma}}_{21}$

Refer to Searle (1971, pp. 46 -47) for details. The correspondence with the conditional spatial simulation problem is as follows. Let the coordinates of the observed data points be denoted $\tilde{s}_1,\tilde{s}_2, ... , \tilde{s}_n$ ,with values $\tilde{z}_1,\tilde{z}_2, ... ,\tilde{z}_n$ .Let $\tilde{Z}$ denote the random vector

$\tilde{Z} = [ Z(\tilde{s}_1)\Z(\tilde{s}_2)\\vdots \Z(\tilde{s}_n)]$

The random vector $\tilde{Z}$ corresponds to X₂, while Z corresponds to X₁. Then ${(Z|\tilde{Z}=\tilde{z}) \sim N_k(\tilde{{\mu}},\tilde{C})}$ as in the previous distribution. The matrix

$\tilde{C} = C_{11}-C_{12}C_{22}^{-1}C_{21}$

is again positive definite, so a Cholesky factorization can be performed.

The dimension n for $\tilde{Z}$ is simply the number of nonmissing observations for the VAR= variable; the values $\tilde{z}_1,\tilde{z}_2, ... ,\tilde{z}_n$ are the values of this variable. The coordinates $\tilde{s}_1,\tilde{s}_2, ... , \tilde{s}_n$ are also found in the DATA= data set, with the variables corresponding to the x and y coordinates identified in the COORDINATES statement. Note that all VAR= variables use the same set of conditioning coordinates; this fixes the matrix C₂₂ for all simulations.

The dimension k for Z is the number of grid points specified in the GRID statement. Since there is a single GRID statement, this fixes the matrix C₁₁ for all simulations. Similarly, C₁₂ is fixed.

The Cholesky factorization $\tilde{C}={LL}^T$ is computed once, as is the mean correction

$\tilde{{\mu}}={\mu}_1+{{\Sigma}}_{12}{{\Sigma}}_{22}^{-1}(x_2-{\mu}_2)$

Note that the means ${\mu}_1$ and ${\mu}_2$ are computed using the grid coordinates s₁,s₂, ... ,s_k, the data coordinates $\tilde{s}_1,\tilde{s}_2, ... , \tilde{s}_n$ ,and the quadratic form specification from the MEAN statement. The simulation is now performed exactly as in the unconditional case. A k ×1 vector of independent standard N(0,1) random variables is generated and multiplied by L, and $\tilde{{\mu}}$ is added to the transformed vector. This is repeated N times, where N is the value specified for the NR= option.

Computational Details

In the computation of $\tilde{{\mu}}$ and ${{\Sigma}}$ described in the previous section, the inverse ${{\Sigma}}_{22}^{-1}$ is never actually computed; an equation of the form

${{\Sigma}}_{22} A = B$

is solved for A using a modified Gaussian elimination algorithm that takes advantage of the fact that ${{\Sigma}}_{22}$ is symmetric with constant diagonal C_z(0) that is larger than all off-diagonal elements. The SINGULAR= option pertains to this algorithm. The value specified for the SINGULAR= option is scaled by C_z(0) before comparison with the pivot element.

Memory Usage

For conditional simulations, the largest matrix held in core at any one time depends on the number of grid points and data points. Using the previous notation, the data-data covariance matrix C₂₂ is n ×n, where n is the number of nonmissing observations for the VAR= variable in the DATA= data set. The grid-data cross covariance C₁₂ is n ×k, where k is the number of grid points. The grid-grid covariance C₁₁ is k ×k. The maximum memory required at any one time for storing these matrices is

max(k(k+1),n(n+1)+2(n×k)) × sizeof( double)

There are additional memory requirements that add to the total memory usage, but usually these matrix calculations dominate, especially when the number of grid points is large.

Chapter Contents
Previous
Next
Top