GaussRF {RandomFields} | R Documentation |
These functions simulate isotropic Gaussian random fields using turning bands, circulant embedding, direct methods, and the random coin method.
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)
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
|
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 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==FALSE)
: 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))
)
model="nugget"
is one possibility to create independent
Gaussian random variables. Without loss of efficiency, any
covariance function with parameter vector
c(mean, 0, nugget, scale, ...)
can also be used. If model="nugget"
is used
the second component of param
must be zero.
param
equals the sill of the variogram.
register
is 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 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 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.
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@uni-bayreuth.de http://www.geo.uni-bayreuth.de/~martin
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.
CovarianceFct
,
DeleteRegister
,
DoSimulateRF
,
GetPracticalRange
,
EmpiricalVariogram
,
mleRF
,
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 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")