GaussRF {RandomFields}R Documentation

Gaussian Random Fields

Description

These functions simulate stationary spatial and spatio-temporal Gaussian random fields using turning bands/layers, circulant embedding, direct methods, and the random coin method.

Usage

GaussRF(x, y=NULL, z=NULL, T=NULL, grid=!missing(gridtriple), model, param,
       trend=NULL, method=NULL, n=1, register=0, gridtriple,
       paired=FALSE, ...)

InitGaussRF(x, y=NULL, z=NULL, T=NULL, grid=!missing(gridtriple), model,
            param, trend=NULL, method=NULL, register=0, gridtriple)

Arguments

x

matrix of coordinates, or vector of x coordinates

y

vector of y coordinates

z

vector of z coordinates

T

vector of time coordinates, may only be given if the random field is defined as an anisotropic random field, i.e. if model=list(list(model=,var=,k=,aniso=),...). T must always be given in the gridtriple format, independently how the spatial part is defined.

grid

logical; determines whether the vectors x, y, and z should be interpreted as a grid definition, see Details. grid does not apply for T.

model

string or list; covariance or variogram model, see CovarianceFct, or type PrintModelList() to get the list of all implemented models; see Details.

param

vector or matrix of parameters or missing, see Details and CovarianceFct; The simplest form is that param is vector of the form param=c(NA,variance,nugget,scale,...), in this order;
The dots ... stand for additional parameters of the model.

trend

trend surface: number (mean) or a vector of length d+1 (linear trend a_0 +a_1 x_1 + ... + a_d x_d), or function; you have the choice of using either x, y, z or X1, X2, X3, ... as spatial variables, as time variable T should be chosen

method

NULL or string; method used for simulating, see RFMethods, or type PrintMethodList() to get all options. If model is given as list then method may not be set if model[[i]]$method, i=1,3,.. is given, and vice versa. However, a global parameter method and specific methods may be given, e.g. list(list(model=..., method="TBM3"), ..., method="ci"); then the specific ones overwrite the global method.

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

paired

logical. If TRUE then the second half of the simulations is obtained by only changing the signs of all the standard Gaussian random variables, on which the first half of the simulations is based. (“Antithetic pairs”.)

...

RFparameters that are locally used only.

Details

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.) See RFMethods for further information. Note that some of the methods do not work for all covariance or variogram models.

GaussRF calls initially 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.

Comments on specific parameters:

Value

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

The object returned GaussRF and DoSimulateRF depends on the parameters n and grid:

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. 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 space and time.

if n > 1 the repetitions make the last dimension

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@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/~schlather

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

References

See RFMethods for the references.

See Also

Covariance, CovarianceFct, DeleteRegister, DoSimulateRF, GetPracticalRange, EmpiricalVariogram, fitvario, 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
 step <- 1    ## nicer, but also time consuming if step <- 0.1
 x <- seq(0, 20, step) 
 y <- seq(0, 20, step)     
 f <- GaussRF(x=x, y=y, model=model, grid=TRUE,
              param=c(mean, variance, nugget, scale, alpha))
 image(x, y, f)


 #############################################################
 ## ... using gridtriple
 step <- 1    ## nicer, but also time consuming if step <- 0.1
 x <- c(0, 20, step)  ## note: vectors of three values, not a 
 y <- c(0, 20, step)  ##       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, Print=5,
              x=x, y=y, z=z, model=model, 
              param=c(mean, variance, nugget, scale, alpha)))


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


 #############################################################
 ## simulating several random fields at once
 step <- 1    ## nicer, but also time consuming if step <- 0.1
 x <- seq(0, 20, step)  # grid
 y <- seq(0, 20, step)
 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])


        

 #############################################################
 ##                                                         ##
 ##      Examples using the extended definition form        ##
 ##                                                         ##
 ##                                                         ##  
 #############################################################


#% library(RandomFields, lib="~/TMP"); RFparameters(Print=6)
x <- (0:100)/10
m <- matrix(c(1,2,3,4),ncol=2)/5
model <- list("$", aniso=m,
              list("*",              
                   list("power", k=2),
                   list("sph"))
              )
z <- GaussRF(x=x, y=x, grid=TRUE, model=model, me="TBM3")
Print(c(mean(as.double(z)),var(as.double(z))))
image(z,zlim=c(-3,3))

## to know more what GaussRF does, use Print
## TMB can be very slow. To trace the iteration, use every
##
z <- GaussRF(x=x, y=x, grid=TRUE, model=model, me="TBM3",
             Print=3, every=100)
image(z,zlim=c(-3,3))
## here, GaussRF uses `direct decomp' to simulate on the line
##       and the square root of the covariance matrix is
##       calculated by the Cholesky decomposition



## non-separable space-time model applied for two space dimensions
## note that tbm method works in some special cases.
#%    library(RandomFields, lib="~/TMP")
x <- y <- seq(0, 7, if (interactive()) 0.05 else 0.2) 
T <- c(1,32,1) * 10        ## note necessarily gridtriple definition
model <- list("$", aniso=diag(c(3, 3, 0.02)),
             list("nsst", k1=2,
                  list("gauss"),
                  list("genB", k=c(1, 0.5))
              ))
z <- GaussRF(x=x, y=y, T=T, grid=TRUE, model=model,
             method="ci", CE.strategy=1,
             CE.trials=if (interactive()) 4 else 1)
rl <- function() if (interactive()) readline("Press return")
for (i in 1:dim(z)[3]) { image(z[,,i], zlim=range(z)); rl();}
for (i in 1:dim(z)[2]) { image(z[,i,], zlim=range(z)); rl();}
for (i in 1:dim(z)[1]) { image(z[i,,], zlim=range(z)); rl();}




 #############################################################
 ##                                                         ##
 ##      Example of a 2d random field based on              ##
 ##      covariance functions valid in 1d only              ##
 ##                                                         ##
 #############################################################

x <- seq(0, 10, 1/10)
model <- list("*",
              list("$",  aniso=matrix(nr=2, c(1, 0)),
                   list("fractgauss", k=0.5)),
              list("$",  aniso=matrix(nr=2, c(0, 1)),
                   list("fractgauss", k=0.9))
              )
z <- GaussRF(x, x,  grid=TRUE, gridtriple=FALSE, model=model)
image(x, x, z)




 #############################################################
 ##                                                         ##
 ##                    Brownian motion                      ##
 ##                (using Stein's method)                   ##
 ##                                                         ##  
 #############################################################
# 1d
kappa <- 1   # in [0,2)
z <- GaussRF(x=c(0, 10, 0.001), grid=TRUE, Print=5,
             model=list("fractalB", kappa))
plot(z, type="l")

# 2d
step <- 0.3  ## nicer, but also time consuming if step = 0.1
x <- seq(0, 10, step)
kappa <- 1   # in [0,2)
z <- GaussRF(x=x, y=x, grid=TRUE, model=list("fractalB", kappa))
image(z,zlim=c(-3,3))

# 3d
x <- seq(0, 3, step)
kappa <- 1   # in [0,2)
z <- GaussRF(x=x, y=x, z=x, grid=TRUE,
             model=list("fractalB", kappa)) 
rl <- function() if (interactive()) readline("Press return")
for (i in 1:dim(z)[1]) { image(z[i,,]); rl();}



 #############################################################
 ## This example shows the benefits from stored,            ##
 ## intermediate results: in case of the circulant          ##
 ## embedding method, the speed is doubled in the second    ##
 ## simulation.                                             ##  
 #############################################################

RFparameters(Storing=TRUE)
y <- x <- seq(0, 50, 0.2)
(p <- c(runif(3), runif(1)+1))
ut <- system.time(f <- GaussRF(x=x,y=y,grid=TRUE,model="exponen",
                              method="circ", param=p))
image(x, y, f)


cat("system time (first call)", format(ut,dig=3),"\n")

# second call with the same parameters can be much faster:
ut <- system.time(f <- DoSimulateRF()) 
image(x, y, f)


cat("system time (second call)", format(ut,dig=3),"\n")


 #############################################################
 ##                                                         ##
 ##     Example how the cutoff method can be set            ##
 ##        explicitly using hypermodels                     ##
 ##                                                         ##
 #############################################################

## NOTE: this feature is still in an experimental stage
##       which has not been yet tested intensively;
## further: parameters and algorithms may change in
##       future.

#% library(RandomFields, lib="~/TMP");source("~/R/cran/RandomFields/tests/source.R")
## simuation of the stable model using the cutoff method
#RFparameters(Print=8, Storing=FALSE)
x <- seq(0, 1, 1/24) # nicer pictures with 1/240
scale <- 1.0
model1 <- list("$", scale=scale, list("stable", alpha=1.0))
rs <- get(".Random.seed", envir=.GlobalEnv, inherits = FALSE)
z1 <- GaussRF(x, x,  grid=TRUE, gridtriple=FALSE,
             model=model1, n=1, meth="cutoff", Storing=TRUE)
(size <- GetRegisterInfo(meth=c("cutoff", "circ"))$S$size)
(cut.off.param <-
   GetRegisterInfo(meth=c("cutoff", "circ"), modelname="cutoff")$param)

image(x, x, z1)


## simulation of the same random field using the circulant
## embedding method and defining the hypermodel explicitely
model2 <- list("$", scale = scale,
               list("cutoff", diam=cut.off.param$diam, a=cut.off.param$a,       
                    list("stable", alpha=1.0))
              )
assign(".Random.seed", rs, envir=.GlobalEnv)
z2 <- GaussRF(x, x,  grid=TRUE, gridtriple=FALSE, model=model2,
              meth="circulant", n=1, CE.mmin=size, Storing=TRUE)
image(x, x, z2)
Print(range(z1-z2)) ## essentially no difference between the fields!




#% library(RandomFields)

 #############################################################
 ## The cutoff method simulates on a torus and a (small)    ##
 ## rectangle is taken as the required simulation.          ##
 ##                                                         ##
 ## The following code shows a whole such torus.            ##
 ## The main part of the code sets local.dependent=TRUE and ##
 ## local.mmin to multiples of the basic rectangle lengths  ##
 #############################################################


# definition of the realisation
RFparameters(CE.useprimes=FALSE)
x <- seq(0, 2, len=4) # better 20
y <- seq(0, 1, len=5) # better 40
grid.size <- c(length(x), length(y))
model <- list("$", var=1.1, aniso=matrix(nc=2, c(2, 1, 0.5, 1)),
              list(model="exp"))

# determination of the (minimal) size of the torus
InitGaussRF(x, y, model=model, grid=TRUE, method="cu")
ce.info.size <- GetRegisterInfo(meth=c("cutoff", "circ"))$S$size
blocks <- ceiling(ce.info.size / grid.size / 4) *4
(size <- blocks * grid.size)

# simulation and plot of the torus 
z <- GaussRF(x, y, model=model, grid=TRUE, method="cu", n=prod(blocks) * 2,
             local.dependent=TRUE, local.mmin=size)[,,c(TRUE, FALSE)]
hei <- 8
do.call(getOption("device"),
        list(hei=hei, wid=hei / blocks[2] / diff(range(y)) *
                         blocks[1] * diff(range(x))))

close.screen(close.screen())
sc <- matrix(nc=blocks[1], split.screen(rev(blocks)), byrow=TRUE)
sc <- as.vector(t(sc[nrow(sc):1, ]))

for (i in 1:prod(blocks)) {
  screen(sc[i])
  par(mar=rep(1, 4) * 0.0)
  image(z[,, i], zlim=c(-3, 3), axes=FALSE, col=rainbow(100)) 
}



##############################################################
##	  Simulating with trend (as function)       	    ##
##############################################################

x <- seq(-5,5,0.1)
z <- GaussRF(x=x, y=x, model = "exponential", param=c(1,0,1), grid=TRUE,
             trend=function(x,y) 3*sin(x)*cos(y))
colors=heat.colors(1000)
image(x,x,z,col=colors)




##############################################################
##	  Simulating with linear trend surface    	    ##
##############################################################

x <- seq(-5,5,0.1)
##trend surface: 3x - y
z <- GaussRF(x=x, y=x, model = "cubic", param=c(1,0,2), grid=TRUE,
             trend=c(0,1,-1))
colors=heat.colors(1000)
persp(x,x,z, phi=30, theta=-3)









[Package RandomFields version 2.0.54 Index]