| GaussRF {RandomFields} | R Documentation |
These functions simulate stationary spatial and spatio-temporal Gaussian random fields using turning bands/layers, circulant embedding, direct methods, and the random coin method.
GaussRF(x, y=NULL, z=NULL, T=NULL, grid, model, param, trend,
method=NULL, n=1, register=0, gridtriple=FALSE, paired=FALSE, ...)
InitGaussRF(x, y=NULL, z=NULL, T=NULL, grid, model, param, trend,
method=NULL, register=0, gridtriple=FALSE)
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 |
Not programmed yet. 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(x) |
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. |
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 where
model is the covariance or variogram model and the
parameter is param=c(mean,variance,nugget,scale, ...).
Alternatively the trend can be given (not programmed yet); then
param=c(variance,nugget,scale, ...).
CovarianceFct. If the trend is not given
it is set to 0.
CovarianceFct.
If the trend is not given it is set to 0.
The method may be specified by the global method
or for each model separately, as additional parameter
method for each entry of the list;
note that methods can not be mixed within a multiplicative part.
If model=list(list(model=,var=,k=,aniso=),...) then a time
component might be given. In case of model="nugget",
aniso must still be given as a matrix. Namely if
aniso is a singular matrix then a zonal nugget effect
is obtained.
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:
grid=FALSE : the vectors x, y,
and z are interpreted as vectors of coordinates
(grid=TRUE) && (gridtriple=FALSE) : the vectors
x, y, and z
are increasing sequences with identical lags for each sequence.
A corresponding
grid is created (as given by expand.grid).
(grid=TRUE) && (gridtriple=TRUE) : the vectors
x, y, and z
are triples of the form (start,end,step) defining a grid
(as given by expand.grid(seq(x$start,x$end,x$step),
seq(y$start,y$end,y$step),
seq(z$start,z$end,z$step)))
register is a parameter which may never be
used by most of the users (please let me know if you use it!). In
other words,
the package will work fine if you ignore this parameter.
The parameter register is of interest in the following
situation. Assume you wish to create sequentially
several realisations of two random fields Z1 and
Z2 that have different
specifications of the covariance/variogram models, i.e.
Z1, Z2, Z1, Z2,...
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.
model and method may
be abbreviated as long as the abbreviations match only one
option. See also PrintModelList() and
PrintMethodList()
RFparameters(...).
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:
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 realisations.
* grid=TRUE. An array of dimension
d+1, where d is the dimension of
the random field, is returned. The last
dimension contains the realisations.
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.
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’)
See RFMethods for the references.
CovarianceFct,
DeleteRegister,
DoSimulateRF,
GetPracticalRange,
EmpiricalVariogram,
fitvario,
MaxStableRF,
RFMethods,
RandomFields,
RFparameters,
ShowModels,
#############################################################
## ##
## 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,
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 forms ##
## ##
## ##
#############################################################
## note that the output seems plausible but not checked!!!!
## tbm may also be used for multiplicate models (if they have
## *exactly* the same anisotropy parameters)
x <- (0:100)/10
m <- matrix(c(1,2,3,4),ncol=2)/5
z <- GaussRF(x=x, y=x, grid=TRUE,
model=list(
list(m="power",v=1,k=2,a=m),
"*", list(m="sph", v=1, a=m)
),
me="TBM3", reg=0,n=1)
print(c(mean(as.double(z)),var(as.double(z))))
image(z,zlim=c(-3,3))
## non-separable space-time model applied for two space dimensions
## note that tbm method does not work nicely, but at least
## in some special cases.
x <- y <- (1:32)/2 ## grid definition, but as a sequence
T <- c(1,32,1)*10 ## note necessarily gridtriple definition
aniso <- diag(c(0.5,8,1))
k <- c(1,phi=1,1,0.5,psi=1,dim=2)
model <- list(list(m="nsst", v=1, k=k, a=aniso))
z <- GaussRF(x=x, y=y, T=T, grid=TRUE, model=model)
rl <- function() if (interactive()) readline("Press return")
for (i in 1:dim(z)[3]) { image(z[,,i]); rl();}
for (i in 1:dim(z)[2]) { image(z[,i,]); rl();}
for (i in 1:dim(z)[1]) { image(z[i,,]); rl();}
#############################################################
## ##
## Example of a 2d random field based on ##
## covariance functions valid in 1d only ##
## ##
#############################################################
x <- seq(0, 10, 1/10)
model <- list(list(model="fractgauss", var=1, kappa=0.5,
aniso=c(1, 0, 0, 0)),
"*",
list(model="fractgauss", var=1, kappa=0.5,
aniso=c(0, 0, 0, 1)))
z <- GaussRF(x, x, grid=TRUE, gridtriple=FALSE, model=model)
image(x, x, z)
#############################################################
## ##
## Brownian motion ##
## (using Stein's method) ##
## ##
#############################################################
# 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="fractalB",
param=c(0,1,0,1,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="fractalB",
param=c(0,1,0,1,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. ##
#############################################################
DeleteAllRegisters()
RFparameters(Storing=TRUE, PrintLevel=1)
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)
hist(f)
c( mean(as.vector(f)), var(as.vector(f)) )
cat("system time (first call)", format(ut,dig=3),"\n")
# second call with the *same* parameters is much faster:
ut <- system.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("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.
## simuation of the stable model using the cutoff method
#RFparameters(Print=8, Storing=FALSE)
x <- seq(0, 1, 1/24)
scale <- 1.0
model1 <- list( list(model="stable", var=1, scale=scale, kappa=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()$method[[1]]$mem$new$method[[1]]$mem$size
#str(GetRegisterInfo(), vec=15)
## simulation of the same random field using the circulant
## embedding method and defining the hypermodel explicitely
model2 <- list(list(model="cutoff", var=1, kappa=c(1, sqrt(2), 1),
scale=scale),
"(",
list(model="stable", var=1, scale=scale, kappa=1.0)
)
assign(".Random.seed", rs, envir=.GlobalEnv)
z2 <- GaussRF(x, x, grid=TRUE, gridtriple=FALSE,
model=model2, n=1, meth="circulant", CE.mmin=size)
print(range(z1-z2)) ## essentially no difference between the fields!
#############################################################
## 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
x <- seq(0, 2, len=20)
y <- seq(0, 1, len=40)
grid.size <- c(length(x), length(y))
model <- list(list(model="exp", var=1.1, aniso=c(2,1,0.5,1)))
# determination of the (minimal) size of the torus
DeleteRegister()
InitGaussRF(x, y, model=model, grid=TRUE, method="cu")
ce.info <- GetRegisterInfo()$method[[1]]$mem$new$method[[1]]$mem
blocks <- ceiling(ce.info$size / grid.size)
size <- blocks * grid.size
# simulation and plot of the torus
DeleteRegister()
z <- GaussRF(x, y, model=model, grid=TRUE, method="cu", n=prod(blocks),
local.dependent=TRUE, local.mmin=size)
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())
split.screen(rev(blocks))
k <- 0
for (j in 1:blocks[2]) {
for (i in 1:blocks[1]) {
k <- k + 1
screen(k)
par(mar=rep(1, 4) * 0.02)
image(z[,,(blocks[2]-j) * blocks[1] + i], zlim=c(-3, 3), axes=FALSE)
}
}