Simulation Techniques

Description

PrintMethodList prints the list of currently implemented methods for simulating random fields

GetMethodNames returns a list of currently implemented methods

Usage

PrintMethodList()

GetMethodNames()

Details

• Boolean functions.
  See marked point processes.
• circulant embedding.
  Introduced by Dietrich & Newsam (1993) and Wood and Chan (1994).

  Circulant embedding is a fast simulation method based on Fourier transformations. It is garantueed to be an exact method for covariance functions with finite support, e.g. the spherical model.

  See also cutoff embedding and intrinsic embedding for variants of the method.

• cutoff embedding.
  Modified circulant embedding method so that exact simulation is garantueed for further covariance models, e.g. the whittle matern model. In fact, the circulant embedding is called with the cutoff hypermodel, see CovarianceFct, and A=B there. cutoff embedding halfens the maximum number of elements models used to define the covariance function of interest (from 10 to 5).

  Here multiplicative models are not allowed (yet).

• direct matrix decomposition.
  This method is based on the well-known method for simulating any multivariate Gaussian distribution, using the square root of the covariance matrix. The method is pretty slow and limited to about 8000 points, i.e. a 20x20x20 grid in three dimensions. This implementation can use the Cholesky decomposition and the singular value decomposition. It allows for arbitrary points and arbitrary grids.
• hyperplane method.
  The method is based on a tessellation of the space by hyperplanes. Each cell takes a spatially constant value of an i.i.d. random variables. The superposition of several such random fields yields approximatively a Gaussian random field.
• intrinsic embedding.
  Modified circulant embedding so that exact simulation is garantueed for further variogram models, e.g. the fractal brownian one. Note that the simulated random field is always non-stationary. In fact, the circulant embedding is called with the Stein hypermodel, see CovarianceFct, and A=B there.

  Here multiplicative models are not allowed (yet).

• Marked point processes.
  Some methods are based on marked point process P = ([x_1,m_1], [x_2,m_2], ...) where the marks m_i are deterministic or i.i.d. random functions on R^d.
  • add.MPP (Random coins).
    Here the functions are elements of the intersection (L1 cap L2) of the Hilbert spaces L1 and L2. A random field Z is obtained by adding the marks:

    Z(.) = sum_i m_i( . - x_i)

    In this package, only stationary Poisson point fields are allowed as underlying unmarked point processes. Thus, if the marks m_i are all indicator functions, we obtain a Poisson random field. If the intensity of the Poisson process is high we obtain an approximate Gaussian random field by the central limit theorem - this is the add.mpp method.

  • max.MPP (Boolean functions).
    If the random functions are multiplied by suitable, independent random values, and then the maximum is taken, a max-stable random field with unit Frechet margins is obtained - this is the max.mpp method.
• nugget.
  The method allows for generating a random field of independent Gaussian random variables. In the isotropic case and if the simple notation of a model (with model and param) is used, this method is called automatically if the nugget effect is positive except the method "circulant embedding" or "direct" have been explicitely.

  The method has been extended to zonal anisotropies, see also parameter nugget.tol in RFparameters.

• Random coins.
  See marked point processes.
• spectral TBM (Spectral turning bands).
  The principle of spectral TBM does not differ from the other turning bands methods. However, line simulations are performed by a spectral technique (Mantoglou and Wilson, 1982); a realisation is given as the cosine with random amplitude and random phase. The implementation allows the simulation of 2-dimensional random fields defined on arbitrary points or arbitrary grids.
• TBM2, TBM3 (Turning bands methods; turning layers).
  It is generally difficult to use the turning bands method (TBM2) directly in the 2-dimensional space. Instead, 2-dimensional random fields are frequently obtained by simulating a 3-dimensional random field (using TBM3) and taking a 2-dimensional cross-section. TBM3 allows for multiplicative models; in case of anisotropy the anisotropy matrices must be multiples of the first matrix or the anisotropy matrix consists of a time component only (i.e. all components are zero except the very last one).
  TBM2 and TBM3 allow for arbitrary points, and arbitrary grids (arbitrary number of points in each direction, arbitrary grid length for each direction).

  Note: Both the precision and the simulation time depend heavily on TBM*.linesimustep and TBM*.linesimufactor that can be set by RFparameters. For covariance models with larger values of the scale parameter, TBM*.linesimufactor=2 is too small.

  The turning layers are used for the simulations with time component. Here, if the model is a multiplicative covariance function then the product may contain matrices with pure time component. All the other matrices must be equal up to a factor and the temporal part of the anisotropy matrix (right column) may contain only zeros, except the very last entry.

Automatic selection algorithm

— details coming soon —

Note

Most methods possess additional parameters, see RFparameters() that control the precision of the result. The default parameters are chosen such that the simulations are fine for many models and their parameters. The example in EmpiricalVariogram() shows a way of checking the precision.

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute

Yindeng Jiang jiangyindeng@gmail.com (circulant embedding methods ‘cutoff’ and ‘intrinsic’)

References

Gneiting, T. and Schlather, M. (2004) Statistical modeling with covariance functions. In preparation.

Lantuejoul, Ch. (2002) Geostatistical simulation. New York: Springer.

Schlather, M. (1999) An introduction to positive definite functions and to unconditional simulation of random fields. Technical report ST 99-10, Dept. of Maths and Statistics, Lancaster University.

Original work:

• Circulant embedding:

  Chan, G. and Wood, A.T.A. (1997) An algorithm for simulating stationary Gaussian random fields. J. R. Stat. Soc., Ser. C 46, 171-181.

  Dietrich, C.R. and Newsam, G.N. (1993) A fast and exact method for multidimensional Gaussian stochastic simulations. Water Resour. Res. 29, 2861-2869.

  Wood, A.T.A. and Chan, G. (1994) Simulation of stationary Gaussian processes in [0,1]^d J. Comput. Graph. Stat. 3, 409-432.

• Intrinsic embedding and Cutoff embedding:

  Stein, M.L. (2002) Fast and exact simulation of fractional Brownian surfaces. J. Comput. Graph. Statist. 11, 587–599.

  Gneiting, T., Sevcikova, H., Percival, D.B., Schlather, M. and Jiang, Y. (2005) Fast and Exact Simulation of Large Gaussian Lattice Systems in R^2: Exploring the Limits J. Comput. Graph. Statist. Submitted.

• Turning bands method (TBM), turning layers:

  Dietrich, C.R. (1995) A simple and efficient space domain implementation of the turning bands method. Water Resour. Res. 31, 147-156.

  Mantoglou, A. and Wilson, J.L. (1982) The turning bands method for simulation of random fields using line generation by a spectral method. Water. Resour. Res. 18, 1379-1394.

  Matheron, G. (1973) The intrinsic random functions and their applications. Adv. Appl. Probab. 5, 439-468.

  Schlather, M. (2004) Turning layers: A space-time extension of turning bands. Submitted

• Random coins:

  Matheron, G. (1967) Elements pour une Theorie des Milieux Poreux. Paris: Masson.

See Also

GaussRF, MaxStableRF, PrintModelList, RandomFields.

Examples

 PrintMethodList()

############################################################
##                                                        ##
##                     Figure 1 in                        ##
##           Gneiting, T., Sevcikova, H.,                 ##
##   Percival, D.B., Schlather, M. and Jiang, Y. (2005)   ##
##                                                        ##
############################################################
## the example below shows for which parameter combinations
## of the stable model, the cutoff method works
## note the method 'cutoff' itself has an automatic
## selection algorithm implemented and is called with
## the stable model only.
## Not run: 

stabletest <- function(alpha, theta, size=512) {
  RFparameters(CE.trials=1, CE.tolIm = 1e-8, CE.tolRe=0, CE.force = FALSE,
               CE.useprimes=TRUE, CE.strategy=0, Print=2)
  model <- list(list(model="cutoff", var=1, kappa=c(1, theta, 1)),
                "(",
                list(model="stable", var=1, scale=1, kappa=alpha)
                )
  check <- CheckAndComplete(model=model, dim=2)
  if (check$error != 0) {
    cat("error", check$error, "\n")
    return(NA)
  }
  r <- paramextract(check$param[ , 1])$cutoffr
#  x <- c(r/size, r,  r /  size)
  x <- c(0, r, r / (size - 1)) * theta
  RFparameters(CE.mmin=rep(2 * size, 2))
  return(if (InitGaussRF(x, x, grid=TRUE, gridtriple=TRUE,
                         model=model, meth="circulant") == 0)
         r else NA)
}

alphas <- seq(1.52, 2.0, 0.01) # 0.02
thetas <- seq(0.05, 3.5, 0.05) #0.05( 2x)
m <- matrix(NA, nrow=length(thetas), ncol=length(alphas))
for (it in 1:length(thetas)) {
  theta <- thetas[it]
  for (ia in 1:length(alphas)) {
    alpha <- alphas[ia]
    cat("alpha=", alpha, "theta=", theta,"\n")
    print(unix.time(m[it, ia] <- stabletest(alpha=alpha, theta=theta)))
    if (is.na(m[it, ia])) break
  }
  if (any(is.finite(m))) image(thetas, alphas, m, col=rainbow(100))
}

## End(Not run)