The functions circulantEmbedding and circulantEmbeddingSetup are more recent fields
functions, more easy to read, and recommended over sim.rf. sim.rf is limited to 2D fields
while circulantEmbedding can handle any number of
dimensions and has some shortcuts to be efficient for the 2D case.
circulantEmbeddingSetup is also the same setup function used
to create covObject for fast multiplication of
a stationary covariance function with values on a regular, 2D grid
(e.g. see stationaryImageCov )
One complication in this code is to take any grid size and enlarge it so that
the grid sizes are powers of 2 and 3 (highly composite). This make the FFT more
efficient. (In fact without any adjustment grid sizes that are prime numbers will
make the FFT quite slow.)
The simulated field has the marginal variance that is determined by
the covariance function for zero distance. Within the fields package the
exponential and Matern functions set this equal to one ( e.g. Matern(0) ==1) so
that one simulates a random field with a marginal variance of one. For
stationary.cov the marginal variance is whatever Covariance(0)
evaluates to and we
recommend that alternative covariance functions also be normalized so
that this is one.
Of course if one requires a Gaussian field with different marginal
variances one can simply scale the result of this function but the process standard deviations.
See the third example below.
Both sim.rf and circulantEmbedding
take an object that includes some preliminary
calculations and so is more efficient for simulating more than one
field from the same covariance.
Adjusting for negative weights
The algorithm using an FFT
known as circulant embedding, may not always work if the correlation
range is large. Specifically the weight function obtained from the FFT of
the covariance field will have some negative values. Currently a lengthy error message
is printed if this is the case.
A simple fix is to increase the size of the domain
so that the correlation scale becomes smaller relative to the extent
of the domain. Increasing the size can be computationally expensive,
however, and so this method has some limitations. But when it works it is
an exact simulation of the random field.
For a stationary model the covariance object ( or list)
for circulantEmbedding should have minimally, the components:
That is
names( obj)
should give
"m", "grid", "M", "wght"
where m is the number of grid points in each dimension, grid is a list
with components giving the grid points in each coordinate.
M is the size of the larger grid that is used for "embedding" and
simulation. Usually M = 2*m and results in an exact
simulation of the stationary Gaussian field. The default if M is not passed
is to find the smallest power of 2 greater than 2*m. wght is an array from
the FFT of the covariance function with dimensions M.
Keep in mind that for the final results only the array
that is within the indices 1: m[i] for each dimension i is retained.
This can give a much larger intermediate array, however, in the computation.
E.g. if m[1] = 100 and m[2]=200 by default then M[1] = 256
and M[2] = 512. A 256 X 512 array
is simluated with to get the 100 by 200 result.
The easiest way to create the
object for simulation is to use circulantEmbeddingSetup.
For the older function sim.rf one uses the image based covariance functions with setup=TRUE to create the list for simulation.
See the example below for this usage.
The classic reference for this algorithm is
Wood, A.T.A. and Chan, G. (1994).
Simulation of Stationary Gaussian Processes in [0,1]^d . Journal of
Computational and Graphical Statistics, 3, 409-432. Micheal Stein and
Tilman Gneiting have also made some additional contributions to the
algortihms and theory.