Gaussian random field simulation by circulant embedding

The covariance \(C(h) = sill\,\rho(|Ah|/range)\) (plus a nugget) on a regular grid is embedded in a torus of g times the grid size whose base block holds C at the wrapped lags; its eigenvalues are the 2-D FFT of that block (Wood and Chan 1994; Dietrich and Newsam 1997). The real part of \(FFT(\sqrt{S/N}(e_1 + i e_2))\) then has covariance exactly C on the grid. g doubles (2, 4, 8) until no eigenvalue is below -1e-8 of the largest; otherwise negatives are clipped and embedding_exact is FALSE. Models: exponential, gaussian, spherical and matern (geoR / fields parameterisation); geometric anisotropy by angle (radians, major axis) and ratio (minor/major). Normals are Philox streams 2k and 2k + 1 of seed (row-major over the torus; nugget stream 2 n_sims + k), identical to the Python arm.

Usage

SpectralGRF(
  coords,
  cov_model = "exponential",
  cov_params = list(),
  n_sims = 1L,
  seed = 42L
)

Arguments

coords

Two-column matrix of regular grid points.

cov_model

"exponential", "gaussian", "matern" or "spherical".

cov_params

List with sill, range, nugget, nu, angle, ratio.

n_sims

Number of realisations.

seed

Philox seed.

Value

List with simulations (n_sims x n), eigenvalues, torus, embedding_exact, min_eigenvalue_ratio.

References

Wood, A. T. A. and Chan, G. (1994). Simulation of stationary Gaussian processes in the unit cube. Journal of Computational and Graphical Statistics 3, 409-432.

Dietrich, C. R. and Newsam, G. N. (1997). Fast and exact simulation of stationary Gaussian processes through circulant embedding of the covariance matrix. SIAM Journal on Scientific Computing 18, 1088-1107.

Examples

g <- as.matrix(expand.grid(y = 0:2 * 0.5, x = 0:3 * 0.5)[, 2:1])
SpectralGRF(g, "matern", list(range = 1, nu = 1.5), seed = 1)$torus
#> [1] 32 24