Simulation of Random Fields – Advanced

Description

This function simulates unconditional random fields:

• univariate and multivariat, spatial and spatio-temporal Gaussian random fields

• stationary Poisson fields

• Chi2 fields

• t fields

• Binary fields

• stationary max-stable random fields.

It also simulates conditional random fields for

• univariate and multivariat, spatial and spatio-temporal Gaussian random fields

For basic simulation of Gaussian random fields, see RFsimulate. See RFsimulate.more.examples and RFsimulate.sophisticated.examples for further examples.

Arguments

Details

RFsimulate simulates different classes of random fields, controlled by the wrapping model.

If the wrapping function of the model argument is a covariance or variogram model (i.e., one of list obtained by RFgetModelNames(type="variogram", group.by="type"), by default, a Gaussian field with the corresponding covariance structure is simulated. By default, the simulation method is chosen automatically through internal algorithms. The simulation method can be set explicitly by enclosing the covariance function with a method specification.

If other than Gaussian fields are to be simulated, the model argument must be enclosed by a function specifying the type of the random field.

There are different possibilities of passing the locations at which the field is to be simulated. If grid=FALSE, all coordinate vectors (except for the time component T) must have the same length and the field is only simulated at the locations given by the rows of x or of cbind(x, y, z). If T is not missing, the field is simulated for all combinations (x[i, ], T[k]) or (x[i], y[i], z[i], T[k]), i=1, ..., nrow(x), k=1, ..., length(T), even if model is not explicitly a space-time model.
If grid=TRUE, the vectors x, y, z and T or the columns of x and T are interpreted as a grid definition, i.e. the field is simulated at all locations (x_i, y_j, z_k, T_l), as given by expand.grid(x, y, z, T). Here, “grid” means “equidistant in each direction”, i.e. all vectors must be equidistant and in ascending order. In case of more than 3 space dimensions, the coordinates must be given in matrix notations. To enable different grid lengths for each direction in combination with the matrix notation, the “gridtriple” notation c(from, stepsize, len) is used: If x, y, z, T or the columns of x are of length 3, they are internally replaced by seq(from=from, to=from+(len-1)*stepsize, by=stepsize) , i.e. the field is simulated at all locations
expand.grid(seq(x$from, length.out=x$len, by=x$stepsize), seq(y$from, length.out=y$len, by=y$stepsize), seq(z$from, length.out=z$len, by=z$stepsize), seq(T$from, length.out=T$len, by=T$stepsize))

If data is passed, conditional simulation is performed.

• if of class RFsp, ncol(data@coords) must equal the dimension of the index space. If data@data contains only a single variable, variable names are optional. If data@data contains more than one variable, variables must be named and model must be given in the tilde notation resp ~ ... (see RFformula) and "resp" must be contained in names(data@data).

• If data is a matrix or a data.frame, either ncol(data) equals (dimension of index space + 1) and the order of the columns is (x, y, z, T, response) or, if data contains more than one response variable (i.e. ncol(data) > (dimension of index space + 1)), colnames(data) must contain colnames(x) or those of "x", "y", "z", "T" that are not missing. The response variable name is matched with model, which must be given in the tilde notation. If "x", "y", "z", "T" are missing and data contains NAs, colnames(data) must contain an element which starts with ‘data’; the corresponding column and those behind it are interpreted as the given data and those before the corresponding column are interpreted as the coordinates.

• if x is missing, RFsimulate searches for NAs in the data and performs a conditional simulation for them.

Specification of err.model: In geostatistics we have two different interpretations of a nugget effect: small scale variability and measurement error. The result of conditional simulation usually does not include the measurement error. Hence the measurement error err.model must be given separately. For sake of generality, any model (and not only the nugget effect) is allowed. Consequently, err.model is ignored when unconditional simulation is performed.

Value

By default, an object of the virtual class RFsp; result is of class RFspatialGridDataFrame if [space-time-dimension > 1] and the coordinates are on a grid, result is of class RFgridDataFrame if [space-time-dimension = 1] and the coordinates are on a grid, result is of class RFspatialPointsDataFrame if [space-time-dimension > 1] and the coordinates are not on a grid, result is of class RFpointsDataFrame if [space-time-dimension = 1] and the coordinates are not on a grid.

The output format can be switched to the "old" array format using RFoptions, either by globally setting RFoptions(spConform=FALSE) or by passing spConform=FALSE in the call of RFsimulate. Then the object returned by RFsimulate depends on the arguments n and grid in the following way:

if vdim > 1 the vdim-variate vector makes the first dimension

if grid=TRUE an array of the dimension of the random field makes the next dimensions. Here, the dimensions are ordered in the sequence x, y, z, T (if given).

Else if no time component is given, then the values are passed as a single vector. Else if the time component is given the next 2 dimensions give the space and the time, respectively.

if n > 1 the repetitions make the last dimension

Note: Conversion between the sp format and the conventional format can be done using the method RFspDataFrame2conventional and the function conventional2RFspDataFrame.

InitRFsimulate returns 0 if no error has occurred and a positive value if failed.

Note

Advanced options are

• spConform (suppressed return of S4 objects)

• practicalrange (forces range of covariances to be one)

• exactness (chooses the simulation method by precision)

• seed (sets .Random.seed locally or globally)

See RFoptions for further options.

Author(s)

Martin Schlather, schlather@math.uni-mannheim.de http://ms.math.uni-mannheim.de/de/publications/software

References

General

• 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.

  Dietrich, C.R. and Newsam, G.N. (1996) A fast and exact method for multidimensional Gaussian stochastic simulations: Extensions to realizations conditioned on direct and indirect measurement Water Resour. Res. 32, 1643-1652.

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

  The code used in RandomFields is based on Dietrich and Newsam (1996).

• 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.

• Markov Gaussian Random Field:

  Rue, H. (2001) Fast sampling of Gaussian Markov random fields. J. R. Statist. Soc., Ser. B, 63 (2), 325-338.

  Rue, H., Held, L. (2005) Gaussian Markov Random Fields: Theory and Applications. Monographs on Statistics and Applied Probability, no 104, Chapman & Hall.

• 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

RFoptions, RMmodel, RFgui, methods for simulating Gaussian random fields, RFfit, RFempiricalvariogram, RFsimulate.more.examples, RFsimulate.sophisticated.examples, RPgauss,

Examples

RFoptions(seed=0) ## *ANY* simulation will have the random seed 0; set
##                   RFoptions(seed=NA) to make them all random again

Results

R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(RandomFields)
Loading required package: sp
Loading required package: RandomFieldsUtils
This is RandomFieldsUtils Version: 0.2.1
This is RandomFields Version: 3.1.16

Attaching package: 'RandomFields'

The following object is masked from 'package:RandomFieldsUtils':

    RFoptions

The following objects are masked from 'package:base':

    abs, acosh, asin, asinh, atan, atan2, atanh, cos, cosh, exp, expm1,
    floor, gamma, lgamma, log, log1p, log2, logb, max, min, round, sin,
    sinh, sqrt, tan, tanh, trunc

> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/RandomFields/RFsimulateAdvanced.Rd_%03d_medium.png", width=480, height=480)
> ### Name: RFsimulateAdvanced
> ### Title: Simulation of Random Fields - Advanced
> ### Aliases: RFsimulateAdvanced
> ### Keywords: spatial
> 
> ### ** Examples
> 
> RFoptions(seed=0) ## *ANY* simulation will have the random seed 0; set
> ##                   RFoptions(seed=NA) to make them all random again
> ## Don't show: 
> StartExample()
> ## End(Don't show)
> ## Don't show: ## Not run: 
> ##D 
> ##D #############################################################
> ##D ##                                                         ##
> ##D ## Example 1: Specification of simulation method           ##
> ##D ##                                                         ##
> ##D #############################################################
> ##D 
> ##D  ## usage of a specific method
> ##D  ## -- the complete list is obtained by RFgetMethodNames()
> ##D  model <- RMstable(alpha=1.5)
> ##D  x <- runif(100, max=20) 
> ##D  y <- runif(100, max=20) # 100 points in 2 dimensional space
> ##D  simulated <- RFsimulate(model = RPdirect(model), x=x, y=y) # cholesky
> ##D  plot(simulated)
> ##D 
> ##D 
> ##D 
> ##D #############################################################
> ##D ##                                                         ##
> ##D ## Example 2: Turnung band with different number of lines  ##
> ##D ##                                                         ##
> ##D #############################################################
> ##D model <- RMstable(alpha=1.5)
> ##D x <- seq(0, 10, 0.01)
> ##D z <- RFsimulate(model = RPtbm(model), x=x, y=x)
> ##D plot(z)
> ##D 
> ##D 
> ##D 
> ##D #############################################################
> ##D ##                                                         ##
> ##D ## Example 3: Shot noise fields (random coins)             ##
> ##D ##                                                         ##
> ##D #############################################################
> ##D 
> ##D  x <- GridTopology(0, .1, 500)
> ##D 
> ##D  z <- RFsimulate(model=RPpoisson(RMgauss()), x=x, mpp.intensity = 100)
> ##D 
> ##D  plot(z)
> ##D  par(mfcol=c(2,1))
> ##D  plot(z@data[,1:min(length(z@data), 1000)], type="l")
> ##D  hist(z@data[,1])
> ##D  
> ##D  
> ##D  z <- RFsimulate(x=x, model=RPpoisson(RMball()), mpp.intensity = 0.1)
> ##D  
> ##D  plot(z)
> ##D  par(mfcol=c(2,1))
> ##D  plot(z@data[,1:min(length(z@data), 1000)], type="l")
> ##D  hist(z@data[,1])
> ##D  
> ##D 
> ##D 
> ##D 
> ##D  #############################################################
> ##D  ##                                                         ##
> ##D  ## Example 4: a 2d random field based on                   ##
> ##D  ## covariance functions valid in 1d only                   ##
> ##D  ##                                                         ##
> ##D  #############################################################
> ##D 
> ##D x <- seq(0, 2, 0.1)
> ##D model <- RMfbm(alpha=0.5, Aniso=matrix(nrow=1, c(1, 0))) + 
> ##D          RMfbm(alpha=0.9, Aniso=matrix(nrow=1, c(0, 1)))
> ##D z <- RFsimulate(x, x, model=model)
> ##D plot(z)
> ##D 
> ##D 
> ##D 
> ##D #############################################################
> ##D ##                                                         ##
> ##D ## Example 5 : Brownian sheet                              ##
> ##D ## (using Stein's method)                                  ##
> ##D ##                                                         ##
> ##D #############################################################
> ##D  
> ##D # 2d
> ##D step <- 0.3 ## nicer, but also time consuming if step = 0.1
> ##D x <- seq(0, 5, step)
> ##D alpha <- 1 # in [0,2)
> ##D z <- RFsimulate(x=x, y=x, model=RMfbm(alpha=alpha))
> ##D plot(z)
> ##D 
> ##D 
> ##D # 3d
> ##D z <- RFsimulate(x=x, y=x, z=x,
> ##D  model=RMfbm(alpha=alpha))
> ##D 
> ##D 
> ##D #############################################################
> ##D ##                                                         ##
> ##D ## Example 5 : Non-Geometric anisotropy                    ##
> ##D ##                                                         ##
> ##D #############################################################
> ##D 
> ##D x <- seq(0.1, 6, 0.12)
> ##D Aniso <- R.c(R.p(1)^2, R.p(2)^1.5)
> ##D z <- RFsimulate(RMexp(Aniso = Aniso) + 10, x, x)
> ##D plot(z) 
> ##D 
> ##D 
> ##D 
> ##D 
> ## End(Not run)## End(Don't show)
> ## Don't show: 
> FinalizeExample()
> ## End(Don't show)
> 
> 
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>