GaussRF: 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

matrix of coordinates, or vector of x coordinates

vector of y coordinates

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 Co

method

NULL or string; Method used for simulating, see RFMethods, or type PrintMethodList() to get all options

number of realisations to generate

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

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.

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:

grid==FALSE: the vectorsx,y, andzare interpreted as vectors of coordinates
(grid==TRUE) && (gridtriple==FALSE): the vectorsx,y, andzare increasing sequences with identical lags for each sequence. A corresponding grid is created (as given byexpand.grid).
(grid==TRUE) && (gridtriple==FALSE): the vectorsx,y, andzare triples of the form (start,end,step) defining a grid (as given byexpand.grid(seq(x$start,x$end,x$step), seq(y$start,y$end,y$step), seq(z$start,z$end,z$step)))
model="nugget"is one possibility to create independent Gaussian random variables. Without loss of efficiency, any covariance function with parameter vectorc(mean, 0, nugget, scale, ...)can also be used. Ifmodel="nugget"is used the second component ofparammust be zero. %% this has to be changed in later versions;
The sum of the components variance and nugget in the argumentparamequals the sill of the variogram.
registeris a parameter which may never be used by most users (please let me know if you use it!). In other words, the package will work fine if you ignore this parameter. The parameterregisteris of interest in the following situation. Assume you wish to create sequentially several realisations of two random fields$Z_1$and$Z_2$that have different specifications of the covariance/variogram models, i.e.$Z_1$,$Z_2$,$Z_1$,$Z_2$,... Then, without using different registers, the algorithm will not be able to profit from already calculated intermediate results, as the specifications of the covariance/variogram model change every time. However, using different registers allows for profiting from up to 10 stored intermediate results.
The strings formodelandmethodmay be abbreviated as long as the abbreviations match only one option. See alsoPrintModelList()andPrintMethodList()
Further control parameters for the simulation are set by means ofRFparameters(...).

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.

Examples

Run this code

#############################################################
 ## 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),"")

 # 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),"")

Run the code above in your browser using DataLab