GaussRF {RandomFields}R Documentation

Gaussian Random Fields

Description

These functions simulate isotropic Gaussian random fields using turning bands, circulant embedding, direct methods, and the random coin method.

Usage

GaussRF(x, y=NULL, z=NULL, grid, model, param, method=NULL,
          n=1, register=0, gridtriple=FALSE)

InitGaussRF(x, y=NULL, z=NULL, grid, model, param,
               method=NULL, register=0, gridtriple=FALSE)

Arguments

x matrix of coordinates, or vector of x coordinates
y vector of y coordinates
z vector of z coordinates
grid logical; determines whether the vectors x, y, and z should be interpreted as a grid definition, see Details.
model string; covariance or variogram model, see CovarianceFct, or type PrintModelList() to get all options
param parameter vector: param=c(mean, variance, nugget, scale,...); the parameters must be given in this order; further parameters are to be added in case of a parametrised class of models, see CovarianceFct
method NULL or string; Method used for simulating, see RFMethods, or type PrintMethodList() to get all options
n number of realisations to generate
register 0:9; place where intermediate calculations are stored; the numbers are aliases for 10 internal registers
gridtriple logical. Only relevant if grid==TRUE. If gridtriple==TRUE then x, y, and z are of the form c(start,end,step); if gridtriple==FALSE then x, y, and z must be vectors of ascending values

Details

GaussRF creates an isotropic Gaussian random field with model as covariance function/variogram model and parameters param=c(mean,variance,nugget,scale,...). The sill of the variogram equals variance + nugget.

GaussRF can use different methods for the simulation, i.e., circulant embedding, turning bands, direct methods, and random coin method. If method==NULL then GaussRF searches for a valid method. GaussRF may not find the fastest method neither the most precise one. It just finds any method among the available methods. (However it guesses what is a good choice.) Note that some of the methods do not work for all covariance or variogram models.

GaussRF is split up in an initial InitGaussRF, which does some basic checks on the validity of the parameters. Then, InitGaussRF performs some first calculations, like the first Fourier transform in the circulant embedding method or the matrix decomposition for the direct methods. Random numbers are not involved. GaussRF then calls DoSimulateRF which uses the intermediate results and random numbers to create a simulation.

When InitGaussRF checks the validity of the parameters, it also checks whether the previous simulation has had the same specification of the random field. If so (and if RFparameters()$STORING==TRUE), the stored intermediate results are used instead of being recalculated.

Using InitGaussRF and DoSimulateRF in sequence might be slightly faster than GaussRF (but less convenient).

Comments on specific parameters:

Value

InitGaussRF returns 0 if no error has occured and a positive value if failed.

GaussRF and DoSimulateRF return NULL if an error has occured; otherwise the returned object depends on the parameters n and grid:
n==1:
* grid==FALSE. A vector of simulated values is returned (independent of the dimension of the random field)
* grid==TRUE. An array of the dimension of the random field is returned.

n>1:
* grid==FALSE. A matrix is returned. The columns contain the repetitions.
* grid==TRUE. An array of dimension d+1, where d is the dimension of the random field, is returned. The last dimension contains the repetitions.

Note

The algorithms for all the simulation methods are controlled by additional parameters, see RFparameters(). These parameters have an influence on the speed of the algorithm and the precision of the result. The default parameters are chosen such that the simulations are fine for many models and their parameters. If in doubt modify the example in EmpiricalVariogram() to check the precision.

Author(s)

Martin Schlather, Martin.Schlather@uni-bayreuth.de http://www.geo.uni-bayreuth.de/~martin

References

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

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.

See Also

CovarianceFct, DeleteRegister, DoSimulateRF, GetPracticalRange, EmpiricalVariogram, mleRF, MaxStableRF, RFMethods, RandomFields, RFparameters, ShowModels.

Examples

 #############################################################
 ## Examples using the symmetric stable model, also called  ##
 ## "powered exponential model"                             ## 
 #############################################################
 PrintModelList()    ## the complete list of implemented models
 model <- "stable"   
 mean <- 0
 variance <- 4
 nugget <- 1
 scale <- 10
 alpha <- 1   ## see help("CovarianceFct") for additional
              ## parameters of the covariance functions
 x <- seq(0, 20, 0.1) 
 y <- seq(0, 20, 0.1)     
 f <- GaussRF(x=x, y=y, model=model, grid=TRUE,
              param=c(mean, variance, nugget, scale, alpha))
 image(x, y, f)

 #############################################################
 ## ... using gridtriple
 x <- c(0, 20, 0.1)  ## note: vectors of three values, not a 
 y <- c(0, 20, 0.1)  ##       sequence
 f <- GaussRF(grid=TRUE, gridtriple=TRUE,
               x=x ,y=y, model=model,  
               param=c(mean, variance, nugget, scale, alpha))
 image(seq(x[1],x[2],x[3]), seq(y[1],y[2],y[3]), f)

 #############################################################
 ## arbitrary points
 x <- runif(100, max=20) 
 y <- runif(100, max=20)
 z <- runif(100, max=20) # 100 points in 3 dimensional space
 f <- GaussRF(grid=FALSE,
              x=x, y=y, z=z, model=model, 
              param=c(mean, variance, nugget, scale, alpha))
 f

 #############################################################
 ## usage of a specific method
 ## -- the complete list can be obtained by PrintMethodList()
 x <- runif(100, max=20) # arbitrary points
 y <- runif(100, max=20)
 f <- GaussRF(method="dir",  # direct matrix decomposition
              x=x, y=y, model=model, grid=FALSE, 
              param=c(mean, variance, nugget, scale, alpha))
 f

 #############################################################
 ## simulating several random fields at once
 x <- seq(0, 20, 0.1)  # grid
 y <- seq(0, 20, 0.1)
 f <- GaussRF(n=3,  # three simulations at once
              x=x, y=y, model=model, grid=TRUE,  
              param=c(mean, variance, nugget, scale, alpha))
 image(x, y, f[,,1])
 image(x, y, f[,,2])
 image(x, y, f[,,3])


 #############################################################
 ## This example shows the benefits from stored,            ##
 ## intermediate results: in case of the circulant          ##
 ## embedding method, the speed is doubled in the second    ##
 ## simulation.                                             ##  
 #############################################################
 DeleteAllRegisters()
 RFparameters(Storing=TRUE,PrintLevel=1)
 y <- x <- seq(0, 50, 0.2)
 (p <- c(runif(3), runif(1)+1))
 ut <- unix.time(f <- GaussRF(x=x,y=y,grid=TRUE,model="exponen",
                              method="circ", param=p))
 image(x, y, f)
 hist(f)
 c( mean(as.vector(f)), var(as.vector(f)) )
 cat("unix time (first call)", format(ut,dig=3),"\n")

 # second call with the *same* parameters is much faster:
 ut <- unix.time(f <- GaussRF(x=x,y=y,grid=TRUE,model="exponen",
                              method="circ",param=p)) 
 image(x, y, f)
 hist(f)
 c( mean(as.vector(f)), var(as.vector(f)) )
 cat("unix time (second call)", format(ut,dig=3),"\n")

[Package Contents]