Theoretical Development :: SAS/STAT(R) 12.1 User's Guide

Theoretical Development

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

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

$\mb {Z} = \left[ \begin{array}{c} Z_1\\ Z_2\\ \vdots \\ Z_ k \end{array} \right]$

with covariance matrix

$\mb {C} = \left( \begin{array}{cccc} C_{11} & C_{12} & \cdots & C_{1k}\\ C_{21} & C_{22} & \cdots & C_{2k}\\ & \ddots & & \\ C_{k1} & C_{k2} & \cdots & C_{kk} \end{array} \right)$

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

$\mb {C}=\mb {LL}’$

where $\mb {L}$ is lower triangular. See Ralston and Rabinowitz (1978, Chapter 9, Section 3-3) for details. This vector $\mb {Z}$ can be generated by the transformation $\mb {Z}=\mb {LW}$ . Here is where the assumption of a Gaussian SRF is crucial. When $\mb {W} \sim N_ k(\Strong{0},\mb {I})$ , then $\mb {Z}=\mb {LW}$ is also Gaussian. The mean of $\mb {Z}$ is

$E(\mb {Z}) = \mb {L}(E(\mb {W})) = \Strong{0}$

and the variance is

$\mbox{Var}(\mb {Z}) = \mbox{Var}(\mb {LW}) = E(\mb {LWW}’\mb {L}’) = \mb {L}E(\mb {WW}’)\mb {L}’ = \mb {LL}’ = \mb {C}$

Consider now an SRF $Z(\bm {s}), \bm {s}\in D \subset \mathcal{R}^2$ , with spatial covariance function $C(\bm {h})$ . Fix locations $\bm {s}_1,\bm {s}_2,\cdots ,\bm {s}_ k$ , and let $\mb {Z}$ denote the random vector

$\mb {Z} = \left[ \begin{array}{c} Z(\bm {s}_1)\\ Z(\bm {s}_2)\\ \vdots \\ Z(\bm {s}_ k) \end{array} \right]$

with corresponding covariance matrix

$\mb {C}_ z = \left( \begin{array}{cccc} C(\Strong{0}) & C(\bm {s}_1-\bm {s}_2) & \cdots & C(\bm {s}_1-\bm {s}_ k)\\ C(\bm {s}_2-\bm {s}_1) & C(\Strong{0}) & \cdots & C(\bm {s}_2-\bm {s}_ k)\\ & \ddots & & \\ C(\bm {s}_ k-\bm {s}_1) & C(\bm {s}_ k-\bm {s}_2) & \cdots & C(\Strong{0}) \end{array} \right)$

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

$Z(\bm {s}) = \mu (\bm {s}) + \varepsilon (\bm {s})$

where $\mu (\bm {s})$ is a quadratic form in the coordinates $\bm {s}=(x,y)$ and the $\varepsilon (\bm {s})$ is an SRF that has the same covariance matrix $\mb {C}_ z$ as previously. In this case, the $\mu (\bm {s}_ i), i=1,\cdots ,k$ , is computed once and added to the simulated vector $\varepsilon (\bm {s}_ i), i=1,\cdots ,k$ , for each realization.

For a conditional simulation, this distribution of

$\mb {Z} = \left[ \begin{array}{c} Z(\bm {s}_1)\\ Z(\bm {s}_2)\\ \vdots \\ Z(\bm {s}_ k) \end{array} \right]$

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

$\mb {X} = \left[ \begin{array}{c} \mb {X}_1\\ \mb {X}_2\\ \end{array} \right]$

$\bmu = \left[ \begin{array}{c} \bmu _1\\ \bmu _2\\ \end{array} \right]$

$\bSigma = \left( \begin{array}{cc} \bSigma _{11} & \bSigma _{12}\\ \bSigma _{21} & \bSigma _{22}\\ \end{array} \right)$

The subvector $\mb {X}_1$ is $k\times 1$ , $\mb {X}_2$ is $n\times 1$ , $\bSigma _{11}$ is $k\times k$ , $\bSigma _{22}$ is $n\times n$ , and $\bSigma _{12}=\bSigma _{21}’$ is $k\times n$ , with . The full vector $\mb {X}$ is partitioned into two subvectors, $\mb {X}_1$ and $\mb {X}_2$ , and $\bSigma$ is similarly partitioned into covariances and cross covariances.

With this notation, the distribution of $\mb {X}_1$ conditioned on $\mb {X}_2=\bm {x}_2$ is $N_ k(\tilde{\bmu },\tilde{\bSigma })$ , with

$\tilde{\bmu }=\bmu _1+\bSigma _{12}\bSigma _{22}^{-1}(\bm {x}_2-\bmu _2)$

and

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

See 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{\bm {s}}_1,\tilde{\bm {s}}_2,\cdots , \tilde{\bm {s}}_ n$ , with values $\tilde{z}_1,\tilde{z}_2,\cdots ,\tilde{z}_ n$ . Let $\tilde{\mb {Z}}$ denote the random vector

$\tilde{\mb {Z}} = \left[ \begin{array}{c} Z(\tilde{\bm {s}}_1)\\ Z(\tilde{\bm {s}}_2)\\ \vdots \\ Z(\tilde{\bm {s}}_ n) \end{array} \right]$

The random vector $\tilde{\mb {Z}}$ corresponds to $\mb {X}_2$ , while $\mb {Z}$ corresponds to $\mb {X}_1$ . Then $\left(\mb {Z}\mid \tilde{\mb {Z}}=\tilde{\mb {z}}\right) \sim N_ k(\tilde{\bmu },\tilde{\mb {C}})$ as in the previous distribution. The matrix

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

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

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

The dimension k for $\mb {Z}$ is the number of grid points specified in the GRID statement. Since there is a single GRID statement, this fixes the matrix $\mb {C}_{11}$ for all simulations. Similarly, $\mb {C}_{12}$ is fixed.

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

$\tilde{\bmu }=\bmu _1+\mb {C}_{12}\mb {C}_{22}^{-1}(\bm {x}_2-\bmu _2)$

The means $\bmu _1$ and $\bmu _2$ are computed using the grid coordinates $\bm {s}_1,\bm {s}_2,\cdots ,\bm {s}_ k$ , the data coordinates $\tilde{\bm {s}}_1,\tilde{\bm {s}}_2,\cdots , \tilde{\bm {s}}_ n$ , and the quadratic form specification from the MEAN statement. The simulation is now performed exactly as in the unconditional case. A $k \times 1$ vector of independent standard random variables is generated and multiplied by $\mb {L}$ , and $\tilde{\bmu }$ is added to the transformed vector. This is repeated N times, where N is the value specified for the NR= option.

The SIM2D Procedure

Theoretical Development