Generate Gaussian process with squared exponential covariance function

Question

In a (stationary) Gaussian Process, values which are closeby are more similar than values far away from each other. The correlation function tends to zero as distance increases. Often, one models the decaying correlation functon $C$ as:

$C(x_i, x_j) = \theta \, e^{-||x_i - x_j||^2}$

I believe this model also underpins the Kriging method of interpolation.

However, how does one generate (i.e. simulate on a computer) a random field with such a property? You may, for simplicity, assume it's a one dimensional function $x(t)$ with mean $\mu = 0$ and standard deviation $\sigma = 1$.

Answer 1

A Gaussian process is a probability distribution over functions, parameterized by a mean function $\mu(x)$ and covariance function $C(x, x')$. For any set of points $\{x_1, \dots, x_n\}$, the corresponding function values $y = [f(x_1), \dots, f(x_n)]^T$ have a joint Gaussian distribution with mean $m = [\mu(x_1), \dots, \mu(x_n)]^T$ and covariance matrix $K$, where $K_{ij} = C(x_i, x_j)$:

$$p(y \mid m, K) = \text{det}(2 \pi K)^{-\frac{1}{2}} \exp \left[ -\frac{1}{2} (y-m)^T K^{-1} (y-m) \right]$$

To generate a Gaussian process, you would simply pick the mean and covariance functions. To sample from this Gaussian process, you would first pick the points $\{x_1, \dots, x_n\}$ at which the function is to be evaluated. Compute $m$ and $K$ as above. Then, generate the function values by sampling vector $y$ from a Gaussian distribution with mean $m$ and covariance matrix $K$.

Answer 2

user20160 is right. Just to add a little technical help: you simulate data from a Gaussian distribution with covariance matrix $K$ by calculating the Cholesky-Decomposition $K = L L^\top$ and generating a vector of independent Gaussian RVs $Z$. Then $Y = L Z$ has the wanted distribution.

In R, you can do this (or directly simulate a Gaussian process) with the help of the packages RandomFields and RandomFieldsUtils.