Simulation and Analysis of Random Fields

Simulation and Analysis of Random Fields [R logo]

Help pages for package ‘RandomFields’ version 1.3.41

CheckAndComplete | Validity of a covariance model
CondSimu | Conditional Simulation
CovarianceFct | Covariance And Variogram Models
DeleteAllRegisters | Deleting Intermediate Results
DeleteRegister | Deleting Intermediate Results
Dev | Choosing the device
DoSimulateRF | Simulation of Random Fields
EmpiricalVariogram | Empirical (Semi-)Variogram
eval.parameters | Interactive menu
FileExists | Files
fitvario | LSQ and Maximum Likelihood Estimation of Random Field Parameters
fitvario.default | LSQ and Maximum Likelihood Estimation of Random Field Parameters
fractal | fractal dimension
fractal.dim | fractal dimension
GaussRF | Gaussian Random Fields
getactionlist | Get input behaviour
getactions | Get input behaviour
GetMethodNames | Simulation Techniques
GetModelInfo | Internal information
GetModelList | Covariance And Variogram Models
GetModelNames | Covariance And Variogram Models
GetPracticalRange | Determination of the practical range
GetRegisterInfo | Internal information
hostname | System calls
Hurst | Hurst coefficient
hurst | Hurst coefficient
InitGaussRF | Gaussian Random Fields
InitMaxStableRF | Max-Stable Random Fields
InitSimulateRF | Simulation of Random Fields
Kriging | Kriging methods
Locator | Graphical Input
LockRemove | Files
long memory dependence | Hurst coefficient
MaxStableRF | Max-Stable Random Fields
mleRF | LSQ and Maximum Likelihood Estimation of Random Field Parameters
parameter.range | Range of model specific parameters
parampositions | Position of the parameters
pid | System calls
PracticalRange | Determination of the practical range
PrintMethodList | Simulation Techniques
PrintModelList | Covariance And Variogram Models
putactions | Set input behaviour
RandomFields | Simulation and Analysis of Random Fields
Readline | Read a Line
regression | Regression plot
RFMethods | Simulation Techniques
RFparameters | Control Parameters
ShowModels | Interactive Choice of Models and Parameters
SimulateRF | Simulation of Random Fields
sleep | Sleep
sleep.milli | Sleep
soil | Soil data of North Bavaria, Germany
useraction | Set input behaviour
Variogram | Covariance And Variogram Models


Internal functions – do not use them directly

Description

CheckXT checks whether the coordinates of the data and related parameters are specified correctly and transforms the coordinates into a standard format

PrepareModel checks whether the parameters of the covariance model and related parameters are specified correctly and transforms the parameters into a standard format

convert.to.readable is the inverse function to PrepareModel; see Details

plotWithCircles displays data values of marked point processes by circles

GetDistributionNames returns the names of the currently available marginal distributions of the random fields

paramextract extracts for some models some parameters from an internal parameter list

Usage

CheckXT(x, y, z, T, grid, gridtriple)
PrepareModel(model, param, timespacedim, trend, method=NULL,
             named=FALSE)
convert.to.readable(l, allowed=c("standard", "nested", "list"))
plotWithCircles(data, factor=1.0, xlim=range(data[,1])+c(-maxr,maxr),
                ylim=range(data[,2])+c(-maxr,maxr),col=1, fill=0, ...)
GetDistributionNames()
paramextract(p, model=c("cutoff"))

Arguments

x | x coordinates
y | y coordinates
z | z coordinates
T | time instances
grid | see GaussRF
gridtriple | see GaussRF
model | see GaussRF
param | see GaussRF
timespacedim | dimension of the random field including the time dimension, if there is any
trend | mean or trend of the random field
method | simulation method
named | logical. If TRUE covnr and param are returned with names
l | list as returned by PrepareModel
allowed | allowed output formats, see CovarianceFct
data | matrix of 3 columns; first two columns give the coordinates, the third the data
factor | enlargement factor for data
xlim | see plot
ylim | see plot
col | border colour of circles
fill | filling colour of circles
... | further graphical parameters
p | internal parameter list; e.g. the columns of CheckAndComplete(...)$param.
model | the name of a covariance model.

Details

convert.to.readable is roughly speaking the inverse function to PrepareModel. convert.to.readable also tries to simplify the model definition, but cannot rediscover the given method for the simulation of the nugget effect in all cases. Due to the simplification in convert.to.readable and the special definition of the nugget effect for nested models, convert.to.readable may return a correct model definition in case of incorrect input, namely if scale is set to 0 in a list definition, see Examples.

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute

See Also

CovarianceFct

Examples

x <- function(...) {
  str(PrepareModel(...))
  cat("--------------------------------\n")
  str(convert.to.readable(PrepareModel(...)))
}
model <- list(list(model="whi", kappa=5, var=2, s=4), "+",
    list(model="whi", kappa=1, var=3, s=0)) ## s=0 should not be used only in
##                             a model definition where the parameters are
##                             are given in a matrix, see the result
x(model=model, ti=1, me="ci")

## since convert.to.readable performs a one-step simplification,
## iterative calls may further simplify the model
xx <- convert.to.readable(PrepareModel(model=model, ti=1, me="ci"))
x(model=xx$mo, pa=xx$pa, ti=1, me=xx$me)

## back to the matrix definition of nested models
str(convert.to.readable(PrepareModel(xx, ti=1), allowed="nested"))

## back to the (correct) list definition
str(convert.to.readable(PrepareModel(xx, ti=1), allowed="list"))
[Package RandomFields version 1.3.41 Index]


Validity of a covariance model

Description

CheckAndComplete checks the validity of a covariance model; further a matrix of internal parameters is returned.

Usage

CheckAndComplete(model, param, dim)

Arguments

model | definition of a covariance function, see CovarianceFct for examples
param | see model
dim | dimension for the space

Value

A list of two parameters is returned:
error a value different from zero signifies an error; to get a readable error information put RFparameters()$PrintLevel to a higher level and rerun.
param (currently for internal purposes only). matrix of the internal parameter list; columns correspond to the simple models and can be extracted using paramextract.

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de, http://www.stochastik.math.uni-goettingen.de/institute

See Also

GaussRF, paramextract, RandomFields, RFMethods.

Examples

 ## for an example see RFMethods


Conditional Simulation

Description

the function returns conditional simulations of a Gaussian random field

Usage

CondSimu(krige.method, x, y=NULL, z=NULL, T=NULL, grid,
         gridtriple=FALSE, model, param, method=NULL, given, data,
         trend, n=1, register=0,
         err.model=NULL, err.param=NULL, err.method=NULL,
         err.register=1, tol=1E-5, pch=".", paired=FALSE, na.rm=FALSE)

Arguments

krige.method | Assumptions on the random field which corresponds to the respective kriging method; currently 'S' (simple kriging) and 'O' (ordinary kriging) are implemented.
x | matrix or vector of x coordinates; points to be kriged.
y | vector of y coordinates.
z | vector of z coordinates.
T | vector in grid triple form for the time coordinates.
grid | logical; determines whether the vectors x, y, and z should be interpreted as a grid definition, see Details.
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.
model | string; covariance model of the random field. See CovarianceFct, or type PrintModelList() to get all options for model. See CovarianceFct for model being a list.
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 covariance functions, see CovarianceFct; the value of mean must be finite in the case of simple kriging, and is ignored otherwise. See CovarianceFct for param being NULL or list.
method | NULL or string; method used for simulating, see RFMethods, or type PrintMethodList() to get all options.
given | matrix or vector of locations where data are available; note that it is not possible to give the points in form of a grid definition.
data | the values measured.
trend | Not programmed yet. (used by universal kriging)
n | number of realisations to generate. If paired=TRUE then n must be even.
register | 0:9; place where intermediate calculations are stored; the numbers are aliases for 10 internal registers; see GaussRF for further details.
err.model | covariance function for the error model. String or list. See model for details.
err.param | parameters for the error model. See also param.
err.method | Only relevant if err.model is not NULL. Then it must be given if and only if method is given; see method for details.
err.register | see register for details.
tol | considered only if grid=TRUE; tolerated distances of a given point to the nearest grid point to be regarded as being zero; see Details.
pch | character. The included kriging procedure can be quite time consuming. The character pch is printed after roughly each 80th part of calculation.
paired | logical. logical. If TRUE then every second simulation is obtained by only changing the signs of the standard Gaussian random variables, the simulation is based on (“antithetic pairs”).
na.rm | logical. If TRUE then NAs are removed from the given data.

Details

The same way as GaussRF the function CondSimu allows for simulating on grids or arbitrary locations. However simulation on a grid is sometimes performed as if the points were at arbitrary locations, what may imply a great reduction in speed. This happens when the given locations do not lay on the specified grid, since in an intermediate step simulation has to be performed simultaneously on both the grid defined by x, y, z, and the locations of given.

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

Value

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 realisations.
* grid=TRUE. An array of dimension d+1, where d is the dimension of the random field as given by x, y, and z, is returned. The last dimension contains the realisations.

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute

References

Chiles, J.-P. and Delfiner, P. (1999) Geostatistics. Modeling Spatial Uncertainty. New York: Wiley.

Cressie, N.A.C. (1993) Statistics for Spatial Data. New York: Wiley.

Goovaerts, P. (1997) Geostatistics for Natural Resources Evaluation. New York: Oxford University Press.

Wackernagel, H. (1998) Multivariate Geostatistics. Berlin: Springer, 2nd edition.

See Also

CovarianceFct, GaussRF, Kriging RandomFields,

Examples

## creating random variables first
## here, a grid is chosen, but any arbitrary points for which
## data are given are fine.  Indeed if the data are given on a
## grid, the grid has to be expanded before calling `CondSimu',
## see below.
## However, locations where values are to be simulated,
## should be given in form of a grid definition whenever
## possible
param <- c(0, 1, 0, 1)
model <- "exponential"

RFparameters(PracticalRange=FALSE)
p <- 1:7
data <- GaussRF(x=p, y=p, grid=TRUE, model=model, param=param)
do.call(getOption("device"), list(height=4,width=4))
do.call(getOption("device"), list(height=4,width=4))
do.call(getOption("device"), list(height=4,width=4))

# another grid, where values are to be simulated
step <- 0.25 # or 0.3
x <-  seq(0, 7, step)

# standardisation of the output
lim <- range( c(x, p) )
zlim <- c(-2.6, 2.6)
colour <- rainbow(100)

## visualise generated spatial data
dev.set(2)
image(p, p, data, xlim=lim, ylim=lim, zlim=zlim, col=colour)

#conditional simulation
krige.method <- "O" ## random field assumption corresponding to
                    ## those of ordinary kriging
cz <- CondSimu(krige.method, x, x,  grid=TRUE,
               model=model, param=param,
               given=expand.grid(p,p),# if data are given on a grid
                                      # then expand the grid first
               data=data)
range(cz)
dev.set(3)
image(x, x, cz, col=colour, xlim=lim, ylim=lim, zlim=zlim)

#conditional simulation with error term
cze <- CondSimu(krige.method, x, x,  grid=TRUE,
                model=model, param=c(0, 1/2, 0, 1),
                err.model="gauss", err.param=c(0, 1/2, 0, 1),
                given=expand.grid(p,p),
                data=data)
range(cze)
dev.set(4)
image(x, x, cze, col=colour, xlim=lim, ylim=lim, zlim=zlim)


Covariance And Variogram Models

Description

CovarianceFct returns the values of a covariance function

Variogram returns the values of a variogram model

PrintModelList prints the list of currently implemented models including the corresponding simulation methods

GetModelList returns a matrix of currently implemented models and their simulation methods

GetModelNames returns a list of currently implemented models

Usage

CovarianceFct(x, model, param, dim=ifelse(is.matrix(x),ncol(x),1),
              fctcall="Covariance")

Variogram(x, model, param, dim=ifelse(is.matrix(x),ncol(x),1))

PrintModelList()

GetModelList(abbr=TRUE)

GetModelNames()

Arguments

x | vector or (n x dim)-matrix. In particular, if the model is isotropic or dim=1 then x is a vector.
model | character or list; if character then name of the covariance function or variogram model - see below, or type PrintModelList() for all options; see Details for the definition of the model by a list.
param | vector or matrix of parameters or missing, see Details and Examples; 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.
dim | dimension of the space in which the model is applied
fctcall | this parameter should not be changed by the user
abbr | logical or numerical. If TRUE the names for the methods are abbreviated. If numerical, abbr gives the number of letters.

Details

The implemented models are in standard notation for a covariance function (variance 1, nugget 0, scale 1) and for positive real arguments x (and t):

  • bessel

    C(x)= 2^a Gamma(a+1)x^(-a) J_a(x)

    The parameter a is greater than or equal to (d-2)/2, where d is the dimension of the random field.

  • Brownian motion
    see fractalB
  • cardinal sine
    see wave
  • cauchy

    C(x)=(1+x^2)^(-a)

    The parameter a is positive. The model possesses two generalisations, the gencauchy model and the hyperbolic model.

  • cauchytbm

    C(x)= (1+(1-b/c)x^a)(1+x^a)^(-b/a-1)

    The parameter a is in (0,2] and b is positive. The model is valid for dimensions d<=c; this has been shown for integer c, but the package allows real values of c.
    It allows for simulating random fields where fractal dimension and Hurst coefficient can be chosen independently. It has negative correlations for b>c and large x.

  • circular

    C(x)=1-2/pi*(x sqrt(1-x^2)+asin(x)) if 0<=x<=1, 0 otherwise

    This isotropic covariance function is valid only for dimensions less than or equal to 2.

  • constant
    Identically constant. Any scale parameter is ignored.
  • cone
    This model is used only for methods based on marked point processes (see RFMethods); it is defined only in two dimensions. The corresponding (boolean) function is a truncated cone with socle. The base has radius 1/2. The model has three parameters, a, b, and c:
    a gives the radius of the top circle of the cone, given as part of the socle radius; a in [0,1).
    b gives the height of the socle.
    c gives the height of the truncated cone.
  • cubic

    C(x)= 1- 7 x^2 + 8.75 x^3 - 3.5 x^5 + 0.75 x^7 if 0<=x<=1, 0 otherwise

    This model is valid only for dimensions less than or equal to 3. It is a 2 times differentiable covariance functions with compact support.

  • cutoff (hypermodel, see also below)

    C(x) = phi(xB), 0 <= xA <= b

    C(x) = b_0 (r - xB)^{2 c}, b <= xA <= bc

    C(x) = 0, bc <= xA

    The cutoff model is a functional of the covariance function phi. Since the model itself is indifferent for scale or anisotropy parameters, the latter must be given only for the submodels. See below for general comments on hypermodels.

    The first parameter, a, gives the number of subsequent models that build phi; b>0, c>0. The parameters r and b_0 are chosen internally such that C is a smooth function. The parameters A and B are the inverse scale parameters for the hypermodel and submodel, respectively. Note that cutoff seldemly works, if A and B are not identical.

    The algorithm that checks the given parameters knows only about some few necessary conditions. Hence it is not ensured that the cutoff-model is a valid covariance function for any choice of phi and the parameters.

    For certain models phi, i.e. stable, whittle and gencauchy, some sufficient conditions are known.

  • dampedcosine

    C(x)= exp(-a x) cos(x)

    This model is valid for dimension 1 iff a>=0, for dimension 2 iff a>=1, and for dimension 3 iff a >= sqrt(3).

  • exponential

    C(x)=exp(-x)

    This model is a special case of the whittlematern model (for a=1/2 there) and the stable class (for a=1).

  • fractalB (fractal Brownian motion)

    gamma(x) = x^a

    The parameter a is in (0,2]. (Implemented for up to three dimensions)

  • FD

    C(k) = (-1)^k G(1-a/2)^2/(G(1-a/2+k)G(1-a/2-k)) for integer k

    and linearly interpolated otherwise. Here, G is the Gamma function. The parameter a is in [-1, 1). The model is defined in 1 dimension only.

    Remark: the fractionally differenced process stems from time series modelling where the grid locations are multiples of the scale parameter.

  • fractgauss

    C(x) = 0.5 (|x+1|^a - 2|x|^a + |x-1|^a)

    This model is the covariance function for the fractional Gaussian noise with Hurst parameter H =a/2, a in (0,2]. In particular, the model is valid only in one dimension.

  • gauss

    C(x)=exp(-x^2)

    This model is a special case of the stable class (for a=2 there). Note that the corresponding function for the random coins method (cf. the methods based on marked point processes in RFMethods) is

    exp(-2 x^2).

    See gneiting for an alternative model that does not have the disadvantages of the Gaussian model.

  • gencauchy (generalised cauchy)

    C(x)= (1+x^a)^(-b/a)

    The parameter a is in (0,2], and b is positive.
    This model allows for simulating random fields where fractal dimension and Hurst coefficient can be chosen independently.

  • gengneiting (generalised gneiting)
    if a=1 then

    C(x)=[1 + (b+1) * x] * (1-x)^(b+1) if 0<=x<=1, 0 otherwise

    if a=2 then

    C(x)= [1 + (b+2) * x + ((b+2)^2-1) * x^2 / 3] * (1-x)^(b+2) if 0<=x<=1, 0 otherwise

    if a=3 then

    C(x)=[1 + (b+3) * x + (2 * (b+3)^2 - 3) * x^2 / 5 + ((b+3)^2 - 4) * (b+3) * x^3 / 15] * (1-x)^(b+3) if 0<=x<=1, 0 otherwise

    The parameter a is a positive integer; here only the cases a=1, 2, 3 are implemented. The parameter b is greater than or equal to (d + 2a +1)/2 where d is the dimension of the random field.

  • gneiting

    C(x)= (1 + 8 s x + 25 s^2 x^2 + 32 s^3 x^3)*(1-s x)^8 if 0 <= x <= 1/s, 0 otherwise

    where s=0.301187465825. This isotropic covariance function is valid only for dimensions less than or equal to 3. It is a 6 times differentiable covariance functions with compact support.
    It is an alternative to the gaussian model since its graph is visually hardly distinguishable from the graph of the Gaussian model, but possesses neither the mathematical and nor the numerical disadvantages of the Gaussian model.
    This model is a special case of gengneiting (for a=3 and b=5 there). Note that, in the original work by Gneiting (1999), s = 10 sqrt(2) / 47 ~= .3008965, a numerical value slightly deviating from the optimal one.

  • gneitingdiff is obsolete, see example below

    C(x)=(1 + 8 x/b + 25 (x/b)^2 + 32 (x/b)^3)*(1-x/b)^8 * 2^{1-a} Gamma(a)^{-1} x^a K_a(x) if 0<=x<=b, 0 otherwise

    This isotropic covariance function is valid only for dimensions less than or equal to 3. The parameters a and b are positive.
    This class of models with compact support allows for smooth parametrisation of the differentiability up to order 6.

  • hyperbolic

    C(x)= c^(-b) (K_b(a*c))^(-1) * (c^2 +x^2)^(0.5 b) * K_b(a sqrt(c^2 + x^2))

    The parameters are such that
    c>=0, a>0 and b>0, or
    c>0 , a>0 and b=0, or
    c>0 , a>=0, and b<0.
    Note that this class is over-parametrised; always one of the three parameters a, c, and scale can be eliminated in the formula. Therefore, one of these parameters should be kept fixed in any simulation study.
    The model contains as special cases the whittlematern model and the cauchy model, for c=0 and a=0, respectively.

  • iacocesare

    C(x, t) = (1+||x||^a+|t|^b)^{-c}

    The parameters a and b take values in [1,2]; the parameters c must be greater than or equal to half the space-time dimension.

  • J-Bessel
    see bessel
  • K-Bessel
    see whittlematern
  • linear with sill
    See power (a=1 there).
  • lgd1 (local-global distinguisher)

    C(x)= 1- b(a+b)^{-1}|x|^a for |x|<= 1 and a(a+b)^{-1}|x|^-b for |x|> 1

    Here b>0 and a is in (0, 1.5-d/2] for dimension d=1,2. The random field has fractal dimension d + 1 - a/2 and Hurst coefficient 1 - b/2 for b in (0,1]

  • mastein (hypermodel for non-separabel space time modelling)

    C(x, t)=Gamma(b + g(t)) Gamma(b + c) / [Gamma(b + g(t) + c) Gamma(b)] W_{b + g(t)}(|x - Vt|)

    Gamma is the Gamma function; g is a variogram; W is the Whittle-Matern model. The first parameter, a, gives the number of subsequent models that build g. Here, the names of covariance models can also be used; the algorithm chooses the corresponding variograms then. The parameter b is the smoothness parameter of the Whittle-Matern model (for t=0) and must be positive. Finally, c must be greater than or equal to half the dimension of x. Instead of the velocity parameter V, the anisotropy matrix for the hyper model is chosen appropriately. Note that the anisotropy matrix must be such that (x,t) is transformed into a purely spatial vector, i.e. the entries in last column of the matrix are all naught. On the other hand, all entries of the anisotropy matrices in the submodels that build gamma is naught except the very last, purely temporal one.

    Note, that for numerical reasons, b + g(t) may not exceed the value 80.0. If exceeded NA is returned or the algorithm fails.

  • matern
    See whittlematern.
  • nsst (Non-Separable Space-Time model)

    C(x,t)= psi(t)^{-f} phi(x / psi(t))

    This model is used for space-time modelling where the spatial component is isotropic.
    phi is the stable model if b=1;
    phi is the whittlematern model if b=2;
    phi is the cauchy model if b=3;
    Here, a is the respective parameter for the model; the restrictions on a are described there.

    The function psi satisfies
    psi^2(t) = (t^c+1)^d if e=1
    psi^2(t) = (d^{-1} t^c+1)/(t^c+1) if e=2
    psi^2(t) = -log(t^c+1/d)/log d if e=3
    The parameter f must be greater than or equal to the genuine spatial dimension of the field. Furthermore, c in (0,2] and d in (0,1). The spatial dimension must be >=1.

  • nsst2

    C(x,t)= psi(t)^{-g} phi(x / psi(t))

    This model is used for space-time modelling where the spatial component is isotropic. Here
    phi is the gencauchy model if c=1.
    The parameters a and b are the respective parameters for the model. The function psi satisfies
    psi^2(t) = (t^d+1)^e if f=1
    psi^2(t) = (e^{-1} t^d+1)/(t^d+1) if f=2
    psi^2(t) = -log(t^d+1/e)/log e if f=3
    The parameter g must be greater than or equal to the genuine spatial dimension of the field. Furthermore, d in (0,2] and e in (0,1]. Necessarily, dim>=2. The spatial dimension must be >=1.

  • nugget

    1(x==0)

    If the model is used in param-definition mode, either param[2], the variance, or param[3], the nugget, must be zero. If the model is used in the list-definition mode, the anisotropy matrix must be given in an anisotropic context, but not the scale parameter in an isotropic context.

  • penta

    C(x)= 1 - 22/3 x^2 +33 x^4 - 77/2 x^5 + 33/2 x^7 - 11/2 x^9 + 5/6 x^11 if 0<=x<=1, 0 otherwise

    valid only for dimensions less than or equal to 3. This is a 4 times differentiable covariance functions with compact support.

  • power

    C(x)= (1-x)^a if 0<=x<=1, 0 otherwise

    This covariance function is valid for dimension d if a >= (d+1)/2. For a=1 we get the well-known triangle (or tent) model, which is valid on the real line, only.

  • powered exponential
    See stable.
  • qexponential

    C(x) = (2 exp(-x)-a exp(-2x))/(2-a)

    The parameter a takes values in [0,1].

  • spherical

    C(x)= 1 - 1.5 x + 0.5 x^3 if 0<=x<=1, 0 otherwise

    This isotropic covariance function is valid only for dimensions less than or equal to 3.

  • stable

    C(x)=exp(-x^a)

    The parameter a is in (0,2]. See exponential and gaussian for special cases.

  • Stein (hypermodel, see also below)

    C(x) = a_0 + a_2 (xB)^2 + phi(xB), 0 <= xA <= b

    C(x) = b_0 (c - xB)^3/(xB), b <= xA <= bc

    C(x) = 0, bc <= xA

    The Stein model is a functional of the covariance function phi. Since the model itself is indifferent for scale or anisotropy parameters, the latter must be given only for the submodels. See below for general comments on hypermodels.

    The first parameter, a, gives the number of subsequent models that build phi; b>0, c>=1. The parameters a_0, a_2 and b_0 are chosen internally such that C becomes a smooth function. The parameters A and B are the inverse scale parameters for the hypermodel and submodel, respectively.

    Note that if A and B are not identical, Stein seldemly works; it may also happen that unsound results are returned without any message of failure.

    The algorithm that checks the given parameters knows only about some few necessary conditions. Hence it is not ensured that the Stein-model is a valid covariance function for any choice of phi and the parameters.

    For certain models phi, i.e. stable, whittle, gencauchy, and the variogram model fractalB some sufficient conditions are known.

  • steinst1 (non-separabel space time model)

    C(x, t) = W_a(y) - <x, z> t W_{a-1}(y) / [(a - 1)(2a + b)]

    Here, W_a is the Whittle-Matern model with smoothness parameter a; b is greater than or equal to the space-time dimension ‘dim’; y = ||(x,t)||. The components of z are given by c, d, ...; the norm of z must less than or equal to 1.

  • symmetric stable
    See stable.
  • tent model
    See power.
  • triangle
    See power.
  • wave

    C(x)=sin(x)/x if x>0 and C(0)=1

    This isotropic covariance function is valid only for dimensions less than or equal to 3. It is a special case of the bessel model (for a=0.5).

  • whittlematern

    C(x)=W_a(x) =2^{1-a} Gamma(a)^{-1} x^a K_a(x),

    The parameter a is positive.
    This is the model of choice if the smoothness of a random field is to be parametrised: if a>=(2m+1)/2 then the graph is m times differentiable.

    The model is a special case of the hyperbolic model (for c=0 there).

Let cov be a model given in standard notation. Then the covariance model applied with arbitrary variance and scale equals

variance * cov( (.)/scale).

For a given covariance function cov the variogram gamma equals

gamma(x) = cov(0) - cov(x).

Note that the value of the covariance function or variogram depends also on RFparameters()$PracticalRange. If the latter is TRUE and the covariance model is isotropic then the covariance function is internally rescaled such that cov(1)~=0.05 for standard parameters (scale=1).

The model and the parameters can be specified by three different forms; the first ‘standard’ form allows for the specification of the covariance model as given above for an isotropic random field. The second form defines isotropic nested models using matrices. The third form allows for defining anisotropic and/or space-time models using lists; here any basic models can arbitrarily be combined by multiplication and summation.

  • model is a string; param is a vector of the form param=c(mean,variance,nugget,scale,...). (These components might be given separately or bound to a simple list passed to model.)

    The first component of param is reserved for the mean of a random field and thus ignored in the evaluation of the covariance function or variogram. The parameters mean, variance, nugget, and scale must be given in this order; additional parameters have to be supplied in case of a parametrised class of models (e.g. hyperbolic, see below), in the order a, b, c.

    Let cov be a model given in standard notation. Then the covariance model applied with arbitrary variance, nugget, and scale equals

    nugget + variance * cov( (.)/scale).

    Some models allow certain parameter combinations only for certain dimensions. As any model valid in d dimensions is also valid in 1 dimension, the default in CovarianceFct and Variogram is dim=1.

  • model is a string; param is a matrix with columns of the form c(variance, scale, ...).

    Except that the entries for the mean and the nugget are missing all explanations given above also apply here. Each column defines a summand of the nested model. A nugget effect is indicated by scale=0; possibly additional parameters are ignored.

  • model is a list as specified below; param is missing.

    model = list(l.1, OP.1, l.2, OP.2, ..., l.n) where n is at most 10 (except cutoff embedding is used, see RFMethods). The lists l.i are all either of the form l.i = list(model=,var=,kappas=,scale=,method=) or of the form l.i = list(model=,var=,kappas=,aniso=,method=). model is a string; var gives the variance; scale is a scalar whereas aniso is a d x d matrix, which is multiplied from the right to the n x d matrix of points; at the transformed points the values of the (isotropic) random field (with scale 1) are calculated. The dimension d of matrix must match the number of columns of x. The models given by l.i can be combined by OP.i="+" or OP.i="*". method is ignored here; it can be set in GaussRF.

  • Hypermodels
    hypermodels are functions or functionals of covariance functions or variograms. The first parameter is always the number of the following covariance models included.

    Important! Hyper models are in an experimental stage:
    (i) the (current) algorithm does not allow for a complete check whether the parameters for a hypermodel are well chosen. So, only use parameter combinations for which you are sure they lead to a positive definite function.
    (ii) behaviour and parameters may change in future version!

Value

CovarianceFct returns a vector of values of the covariance function.
Variogram returns a vector of values of the variogram model.
PrintModelList prints a table of the currently implemented covariance functions and the matching methods. PrintModelList returns NULL.
GetModelNames returns a list of implemented models

Author(s)

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

References

Overviews:

  • Chiles, J.-P. and Delfiner, P. (1999) Geostatistics. Modeling Spatial Uncertainty. New York: Wiley.
  • Gneiting, T. and Schlather, M. (2004) 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.
  • Schlather, M. (2002) Models for stationary max-stable random fields. Extremes 5, 33-44.
  • Yaglom, A.M. (1987) Correlation Theory of Stationary and Related Random Functions I, Basic Results. New York: Springer.
  • Wackernagel, H. (1998) Multivariate Geostatistics. Berlin: Springer, 2nd edition.

Cauchy models, generalisations and extensions

  • Gneiting, T. and Schlather, M. (2004) Stochastic models which separate fractal dimension and Hurst effect. SIAM review 46, 269-282.

Gneiting's models

  • Gneiting, T. (1999) Correlation functions for atmospheric data analysis. Q. J. Roy. Meteor. Soc., Part A 125, 2449-2464.

Holeeffect model

  • Zastavnyi, V.P. (1993) Positive definite functions depending on a norm. Russian Acad. Sci. Dokl. Math. 46, 112-114.

Hyperbolic model

  • Shkarofsky, I.P. (1968) Generalized turbulence space-correlation and wave-number spectrum-function pairs. Can. J. Phys. 46, 2133-2153.

Iaco-Cesare model

  • de Cesare, L., Myers, D.E., and Posa, D. (2002) FORTRAN programs for space-time modeling. Computers & Geosciences 28, 205-212.
  • de Iaco, S.. Myers, D.E., and Posa, D. (2002) Nonseparable space-time covariance models: some parameteric families. Math. Geol. 34, 23-42.

lgd

  • Gneiting, T. and Schlather, M. (2004) Stochastic models which separate fractal dimension and Hurst effect. SIAM review

Ma-Stein model

  • Ma, C. (2003) Spatio-temporal covariance functions generated by mixtures. Math. Geol., 34, 965-975.
  • Stein, M.L. (2005) Space-time covariance functions. JASA, 100, 310-321.

nsst

  • Gneiting, T. (2001) Nonseparable, stationary covariance functions for space-time data, JASA 97, 590-600.
  • Gneiting, T. and Schlather, M. (2001) Space-time covariance models. In El-Shaarawi, A.H. and Piegorsch, W.W.: The Encyclopedia of Environmetrics. Chichester: Wiley.

Power model

  • Golubov, B.I. (1981) On Abel-Poisson type and Riesz means, Analysis Mathematica 7, 161-184.
  • Zastavnyi, V.P. (2000) On positive definiteness of some functions, J. Multiv. Analys. 73, 55-81.

fractalB

  • Stein, M.L. (2002) Fast and exact simulation of fractional Brownian surfaces. J. Comput. Graph. Statist. 11, 587-599.

See Also

EmpiricalVariogram, GetPracticalRange, parameter.range, RandomFields, RFparameters, ShowModels.

Examples

 PrintModelList()
 x <- 0:100

 # the following five model definitions are the same!
 ## (1) very traditional form
 (cv <- CovarianceFct(x, model="bessel", c(NA,2,1,5,0.5)))

 ## (2) traditional form in list notation
 model <- list(model="bessel", param=c(NA,2,1,5,0.5))
 cv - CovarianceFct(x, model=model)

 ## (3) nested model definition
 cv - CovarianceFct(x, model="bessel",
                    param=cbind(c(2, 5, 0.5), c(1, 0, 0)))

 #### most general notation in form of lists
 ## (4) isotropic notation
 model <- list(list(model="bessel", var=2, kappa=0.5, scale=5),
               "+",
               list(model="nugget", var=1))
 cv - CovarianceFct(x, model=model)

 ## (5) anisotropic notation
 model <- list(list(model="bessel", var=2, kappa=0.5, aniso=0.2),
               "+",
               list(model="nugget", var=1, aniso=1))
 cv - CovarianceFct(as.matrix(x), model=model)

 # The model gneitingdiff was defined in RandomFields v1.0.
 # This isotropic covariance function is valid for dimensions less
 # than or equal to 3 and has two positive parameters.
 # It is a class of models with compact support that allows for
 # smooth parametrisation of the differentiability up to order 6.
 # The former model `gneitingdiff' must now be coded as
 gneitingdiff <-  function(p){
   list(list(m="gneiting", v=p[2], s=p[6]*p[4]),
       "*",
       list(m="whittle", k=p[5], v=1.0, s=p[4]),
       "+",
       list(m="nugget", v=p[3]))
 }
 # and then
 param <- c(NA, runif(5,max=10))
 CovarianceFct(0:100,gneitingdiff(param))
 ## instead of formerly CovarianceFct(x,"gneitingdiff",param)

 # definition of a hypermodel is more complex
 model <- list(list(model="mastein", var=1,
                    aniso=c(1, -0.5, 0, 0), kappa=c(1, 0.5, 1.5)),
              "(",
              list(model="exp", var=1, aniso=c(0, 0, 0, 1)))
 CovarianceFct(cbind(0:10, 0:10), model=model)


Choosing the device

Description

Dev chooses between the graphical devices screen, postscript and pdf

Usage

Dev(on, dev, ps=NULL, cur.cex=TRUE, paper="special", width=5, height=5,
    quiet=FALSE, innerwidth, innerheight, mai, horizontal = FALSE, ...)

Arguments

on | logical. Indicates whether dev should be switched on or off
dev | see Details
ps | name of the pdf or postscript file
cur.cex | logical. If TRUE the par parameters of the current device are used, not the standard parameters
paper | kind of paper. Postscript parameter
width | width of figure. Postscript and pdf parameter
height | height of figure. Postscript and pdf parameter
quiet | logical. If FALSE additional information is given.
innerwidth | height of graphic without the margins; overwrites height if given – experimental state
innerheight | width of graphic without the margins; overwrites width if given – experimental state
mai | parameters of function par – experimental state
horizontal | logical. See postscript
... | further parameters for pdf or postscript

Details

The parameter dev might be

logical
If suffix of ps is either "eps", "ps" or "pdf" the respective file is created. Otherwise, a postscript file ending with suffix "eps" is created if dev=TRUE and a pdf file with suffix "pdf" if dev=FALSE
character
A function with name dev is called, and the suffix dev is added to ps
numeric
Dev switches to the device with number dev;
if such a device does not exist, a new X11 device is created.
If dev<2 the device opened by Dev(on=TRUE,...) is
closed by Dev(on=FALSE); otherwise only
par(new=F)
is applied to the device.
The parameters dev, ps, cur.cex, paper, width, height and ... are ignored if on=FALSE.

Value

NULL. Side effect is that the global variable .dev.orig is created when Dev(on=TRUE,...) is called.

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute

Examples

  ## first an eps-file test.eps is created, then a jpeg-file,
  ## finally the figure is plotted on the screen
  dir(pattern="test*")
  dev.list <- list(TRUE, 1)
  if (interactive()) dev.list <- c(dev.list, "jpeg")
  for (dev in dev.list) {
    print(dev)
    size <- if (dev=="jpeg") 450 else 5
    err <- class(try(Dev(TRUE, dev, ps="test", height=size, width=size)))
    if (err!="try-error") {
      plot(0, 0, main=paste("dev=", dev))
      # readline("press return")
      Dev(FALSE)
    }
  }
  dir(pattern="test*")


Empirical (Semi-)Variogram

Description

EmpiricalVariogram calculates the empirical (semi-)variogram of a random field realisation

Usage

EmpiricalVariogram(x, y=NULL, z=NULL, T=NULL, data, grid, bin,
                   gridtriple=FALSE, phi, theta, deltaT)

Arguments

x | vector of x-coordinates, or matrix
y | vector of y-coordinates
z | vector of z-coordinates
T | vector of time components; here T is given in grid format, see GaussRF.
data | vector or matrix of data; if data has a multiple number of components as expected by the definition of the coordinates then it is assumed that the data stem from repeated, independent measurements at the given locations; the empirical variogram is calculated for the repeated data.
grid | logical; if TRUE then x, y, and z define a grid; otherwise x, y, and z are interpreted as points
bin | vector of ascending values giving the bin boundaries
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
phi | vector of two components. First component gives the angle for the first line of midpoints of an angular variogram. The second component gives the number of directions (on the half circle). The spatial dimension must be at least 2.
theta | vector of two components. First component gives the angle for the first line of midpoints of an angular variogram (angle is zero for the xy-plane). The second component gives the number of directions (on the half circle). The spatial dimension must be at least 3.
deltaT | vector of two components. First component gives the largest temporal distance; the second component the grid length, that must be a multiple of T[3].

Details

Comments on specific parameters:

  • data: the number of values must match the number of points (given by x, y, z, grid, and gridtriple). That is, it must equal the number of points or be a multiple of it. In case the number of data equals n times the number of points, the data are interpreted as n independent realisations for the given set of points.
  • (grid=FALSE): the vectors x, y, and z, are interpreted as vectors of coordinates
  • (grid=TRUE) && (gridtriple=FALSE): the vectors code{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)))
  • The bins are left open, right closed intervals, i.e., (b[i],bin[i+1]] for i=1,...,length(bin)-1. Hence, to include zero, bin[1] must be negative.

Value

The function returns a list:

centers | central points of the bins
emp.vario | empirical variogram; vector or matrix or array, depending on the anisotropy definitions. The sequence is distances, phi, theta, Tbins. If phi, theta, or Tbins below are not given, the respective dimensions are missing.
sd | sd of the variogram cloud within each bin
n.bin | number of points within a bin
phi | vector of angles in xy plane
theta | vector of angles in the third dimensions
Tbins | vector of temporal distances

The first four elements are vectors of length (length(bin)-1).

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute

See Also

GaussRF, fitvario, and RandomFields

Examples

  #############################################################
  ## this example checks whether a certain simulation method ##
  ## works well for a specified covariance model and         ##
  ## a configuration of points                               ##
  #############################################################
  x <- seq(0, 10, 0.5)
  y <- seq(0, 10, 0.5)
  gridtriple <- FALSE      ## see help("GaussRF")
  model <- "whittle"       ## whittlematern
  bins <- seq(0, 5, 0.001)
  realisations <- 5 ## by far too small to get reliable results!!
                   ## It should be of order 500,
                   ## but then it will take some time
                   ## to do the simulations
  param <- c(mean=1, variance=10, nugget=5, scale=2, alpha=2)
  f <- GaussRF(x=x, y=y, grid=TRUE, gridtriple=gridtriple,
               model=model, param=param, method="TBM3",
               n=realisations)
  binned <- EmpiricalVariogram(x=x, y=y, data=f, grid=TRUE,
                               gridtriple=gridtriple, bin=bins)
  truevariogram  <- Variogram(binned$c, model, param)
  matplot(binned$c, cbind(truevariogram,binned$e), pch=c("*","e"))
  ##black curve gives the theoretical values


Interactive menu

Description

eval.parameters provides an interactive menu on a X11 graphical device of R

Usage

eval.parameters(variable, entry, update, simulate, dev, create = TRUE,
                col.rect = "red", col.bg = "blue", col.sep = "grey",
                col.left = "red", col.mid = "white", col.right = "white",
                col.line = "red", col.txt = "black",
                cex=0.8, cex.i=cex, sep="-----", ...)

Arguments

variable | string. The name of the variable to be changed. The name consist of $ and [[]] expression pointing to sublists of a list. The complete list must be given by name in ...
entry | A list of lists. See Details.
update | logical. If TRUE then simulate is called after each interactive input.
simulate | function that is called if simulations are to updated. The parameters must equal the variables given by ...; the function must return the complete list indicated by variable.
dev | Before calling eval.parameters() split.screen must have been called. dev gives the screen on which the interactive menu should be plotted.
create | logical. If TRUE missing list elements of variable are created.
col.rect | colour of the button for free input.
col.bg | colour of a interactive bar
col.sep | colour of the separating line
col.left | colour of preceding element
col.mid | colour for the message
col.right | colour of subsequent element
col.line | colour of the marking line in interactive bars of absolute choice.
col.txt | colour of headers
cex | font height of headers
cex.i | font height of text for elements
sep | style of added characters in separating line.
... | The input variables given by name; the names may not start with a dot; There the complete list to which variable refers must be given. Further additional parameters for the function simulate.

Details

eval.parameters shows a menu list on X11. Depending on the mode of the variables the menu bars have a different appearance and behave differently if the user clicks on the bar. Most of the menu bars have a small rectangle on the right hand side. If this rectangle is pressed the input of a variable is expected in the xterm where R is run.

entry is a list of lists. Each list may contain the following elements:

  • name : header for menu button if var is not NULL; otherwise printed instead of a menu button
  • var :
    • NULL : if val=NULL then it is a separating line in colour col.sep; name is surrounded by sep; all other elements of the list are ignored. If val is a string then val is interpreted as a function; a special string is "simulate", which entails the call of the function simulate with the appropriate parameters.
    • string : selected element of the list that is given by variable. A special string "undo" will be installed to undo things.
  • val :
    • function(d, v) gives the update for var. If v is missing, a starting value (for d=1/2) is expected. Otherwise, v is the current value of var and d is the choice of the user, a value in [0,1]
    • TRUE : logical bar.
    • FALSE : logical bar. The value for var is negated before shown. In the menu, the negative value for var is shown.
    • NULL : a string is read from the terminal inot var.
    • character (vector): if var is given then this vector of strings is interpreted as belonging to a categorical variable 1,...,length(val) and var gives the number of the selected elements. If var=NULL then val is interpreted as a function; a special string is "simulate", which entails the call of the function simulate with the appropriate parameters.
  • delta : logical. In the menu bar absolute values are plotted if delta=FALSE and increments otherwise. Only considered if val is a function.
  • cond : The menue points is shown only if the given condition is satisfied. The condition must be expressed by named elements of the list variable, see example.
  • col : colour that overwrites the standard colour for the rectangle or the separating text.
  • update : logical. If not missing, its value overwrites locally the value of the global parameter update.
  • ... : additional parameters for simulate that overwrite the values given in ... in the call of eval.parameters.

Value

The first variable given in “...”, which is a list. To this list the entry .history is added.
If the users enters ‘exit immediately’ at any point, the program stops with an error message.

Note

To the list given by variable the element .history is added. .history is a list that contains the history of the user input. Each element is a list where the first entry is the number of the menu, the second and the third entries are the former and the new value. Exception: for entries with character val, the value of val is returned as second entry. Consequently, the name .history should not be used for other purposes in variable.

Further, any variable name given in ... must start with a letter.

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute

See Also

useraction

Examples

  ## the following lines define a menu that does not make
  ## too much sense, but shows the various kinds of buttons

  quadratic <- function(d, v, a, mini=0, maxi=Inf) {
    d <- pmin(1, pmax(0, d)) - 0.5
    d <- ((d>0) * 2 - 1) * d^2 * a * 4
    if (missing(v)) d else pmax(mini, pmin(maxi, v + d))
  }

  simulate <- function(H, par) { ## not a serious example
    print(c(H$x$var, par, runif(1)))
    return(H)  ## the function must return the first parameter
  }

  entry <- list(
     list(name="Nonsense Menu"),
     list(name="Simulate!", val="simulate", col="blue"),
     list(name="show H", val="str(H)", col="blue"),
     list(name="colx", var="colour",
          val=c("red", "green", "blue", "brown")),
     list(name="open", var="closed", val=FALSE, par=4.5),
     list(name="modifying", var="modify", val=TRUE, par=5),
     list(name="probability", var="probab", delta=FALSE,
          val=function(d, v) pmin(1, pmax(0, d))),
     list(name="variance", var="var", delta=TRUE,
          val=function(d, v) quadratic(d, v, 10)),
     list(name="name", var="name", par=3, cond="modify")
    )

  scr <- split.screen(rbind(c(0, 0.45, 0, 1), c(0.5, 1, 0, 1)))
  ## before proceeding make sure that both the screen and the xterm
  ## are completely visible

  H <- list(modify=5, x=list()) # note that in this example eval.parameters
  ##          will be called by by H$x, hence modify=5 will be left
  ##          unchanged.
  options(locatorBell=FALSE)
  useraction("start.register")  ## registring the user's input
  str(eval.parameters("H$x", entry, simulate, update=TRUE, dev=scr[2],
     H=H, par=17)) # do not forget to call by name
  getactions()

  ## replay the user's input
  useraction("replay")
  str(eval.parameters("H$x", entry, simulate, update=TRUE, dev=scr[2],
     H=H, par=17))


Files

Description

The function FileExists checks whether a file or a lock-file exists

The function LockRemove removes a lock-file

Usage

   FileExists(file, PrintLevel=RFparameters()$Print)

   LockRemove(file)

Arguments

file | name of the data file
PrintLevel | if PrintLevel<=1 no messages are displayed

Details

FileExists checks whether file or file.lock exists. If none of them exists file.lock is created and hostname and PID are written into file.lock. This is useful if several processes use the same directory. Further, it is checked whether another process has tried to create the same file in the same instance. In this case FileExists returns for at least one of the processes that file.lock has already been created.

Value

FileExists returns

1 | if file already exists
2 | if file.lock already exists
3 | if file.lock was tried to be created, but another process inferred and got priority
0 | otherwise, file and file.lock did not exist and file.lock has been created

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute


LSQ and Maximum Likelihood Estimation of Random Field Parameters

Description

The function estimates arbitrary parameters of a random field specification with various methods.

Usage

fitvario(x, y=NULL, z=NULL, T=NULL, data, model, param,
         lower=NULL, upper=NULL, sill=NA, ...)

fitvario.default(x, y=NULL, z=NULL, T=NULL, data, model, param,
         lower=NULL, upper=NULL, sill=NA, trend,
         use.naturalscaling=TRUE, PrintLevel=RFparameters()$Print,
         optim.control=NULL, bins=20, nphi=1, ntheta=1, ntime=20,
         distance.factor=0.5,
         upperbound.scale.factor=10, lowerbound.scale.factor=20,
         lowerbound.scale.LS.factor=5,
         upperbound.var.factor=10, lowerbound.var.factor=100,
         lowerbound.sill=1E-10, scale.max.relative.factor=1000,
         minbounddistance=0.001, minboundreldist=0.02,
         approximate.functioncalls=50, refine.onborder=TRUE,
          pch=RFparameters()$pch,
         var.name="X",
         time.name="T", transform=NULL, standard.style=NULL,
         lsq.methods=c("self", "plain", "sqrt.nr", "sd.inv", "internal"),
         mle.methods=c("ml", "reml"),
         cross.methods=NULL,
         users.guess=NULL, users.beta=NULL, table.format=FALSE)

Arguments

x | (n x 2)-matrix of coordinates, or vector of x-coordinates
y | vector of y coordinates
z | vector of z coordinates
T | vector of T coordinates; these coordinates are given in triple notation, see GaussRF
data | vector or matrix of values measured at coord; If also a time component is given, the in the data the indices for the spatial components run the fastest.
model | string or list; covariance model, see CovarianceFct, or type PrintModelList() to get all options. If model is a list, then the parameters with value NA are estimated. Parameters that have value NaN should be explicitely be defined by the function transform. An alternative to define NaN values and the function transform, is to replace the NaN by a real-valued function with solely parameter a list defining a covariance model. In case of the anisotropy matrix, the matrix must be replaced by a list if functions are introduced. Only the list elements variance, scale or anisotropy, and kappas can be used, and not the mean or the trend. Further, the mean or the trend cannot be set by such a function. See also transform below.
param | vector or matrix or NULL. If vector then 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 covariance functions, see CovarianceFct. Any components set to NA are estimated; the others are kept fix. See also model above.
lower | list or vector. Lower bounds for the parameters. If lower and param are vectors and length(lower) < length(param) then lower must match the number of additional parameters a,b,c,.... If param is matrix the length of lower must match the number columns of param or being 2 elements smaller (then lower is filled with NA from the left. The bounds are equally applied to all rows. If lower is a list, then elements that are not given are considered as NA. If lower is not given, or lower contains NA then the missing bounds are generated automatically.
upper | list or vector. Upper bounds for the parameters. See also lower.
sill | If not NA the sill is kept fix. Only used if the standard format for the covariance model is given. See Details.
trend | Not programmed yet. May only be set if missing(param); linear formula : uses X1, X2,... and T as internal parameters for the coordinates; all parameters are estimated matrix : must have the same number of rows as x fixed mean + matrix or linear formula : not possible within this function (just subtract the mean from your data before calling this function)
... | arguments as given in fitvario.default and listed in the following.
use.naturalscaling | logical. Only used if model is given in standard (simple) way. If TRUE then internally, rescaled covariance functions will be used for which cov(1)~=0.05. use.naturalscaling has the advantage that scale and the form parameters of the model get ‘orthogonal’, but use.naturalscaling does not work for all models. See Details.
PrintLevel | level to which messages are shown. See Details.
optim.control | control list for optim, which uses 'L-BFGS-B'. However 'parscale' may not be given.
bins | number of bins of the empirical variogram. See Details.
nphi | scalar or vector of 2 components. If it is a vector then the first component gives the first angle of the xy plane and the second one gives the number of directions on the half circle. If scalar then the first angle is assumed to be zero. Note that a good estimation of the variogramm by LSQ with a anisotropic model a large value for ntheta might be needed (about 20).
ntheta | scalar or vector of 2 components. If it is a vector then the first component gives the first angle in the third direction and the second one gives the number of directions on the half circle. If scalar then the first angle is assumed to be zero. Note that a good estimation of the variogramm by LSQ with a anisotropic model a large value for ntheta might be needed (about 20).
ntime | scalar or vector of 2 components. if ntimes is a vector, then the first component are the maximum time distance (in units of the grid length T[3]) and the second component gives the step size (in units of the grid length T[3]). If scalar then the step size is assumed to 1 (in units of the grid length T[3].
distance.factor | relative right bound for the bins. See Details.
upperbound.scale.factor | relative upper bound for scale in LSQ and MLE. See Details.
lowerbound.scale.factor | relative lower bound for scale in MLE. See Details.
lowerbound.scale.LS.factor | relative lower bound for scale in LSQ. See Details.
upperbound.var.factor | relative upper bound for variance and nugget. See Details.
lowerbound.var.factor | relative lower bound for variance. See Details.
lowerbound.sill | absolute lower bound for variance and nugget. See Details.
scale.max.relative.factor | relative lower bound for scale below which an additional nugget effect is detected. See Details.
minbounddistance | absolute distance to the bounds below which a part of the algorithm is considered as having failed. See Details.
minboundreldist | relative distance to the bounds below which a part of the algorithm is considered as having failed. See Details.
approximate.functioncalls | approximate evaluations of the ML target function on a grid. See Details.
refine.onborder | logical. If refine.onborder=TRUE and if the result of any maximum likelihood method or cross validation method is on a borderline, then the optimisation is redone in a modified way (which takes about double extra time)
pch | character shown before evaluating any method; if pch!="" then one or two additional steps in the MLE methods are marked by “+”. Default: "*".
var.name | basic name for the coordinates in the formula of the trend. Default: ‘X’
time.name | basic name for the time component in the formula of the trend. Default: ‘X’
transform | vector of strings. Essentially, transform allows for the definition of a parameter as a function of other estimated or fixed parameters. All the parameters are supposed to be in a vector called ‘param’ where the positions are given by parampositions. An example of transform is function(param) {param[3] <- 5 - param[1]; param}. Any parameter that is set by transform, should be NaN in the model definition. If it is NA a warning is given. Note that the mean and the trend of the model can be neither set nor used in transform. See also standard.style. Instead of giving transform, in the model definition, all NaN values are replaced by functions whose only parameter is a bare model list, i.e., only the list elements variance, scale or anisotropy, and kappas can be used, and not the mean or the trend. Further, the mean or the trend cannot be set by such a function. Default: NULL
standard.style | logical or NULL. This variable should only be set by the advanced user. If NULL then standard.style will be TRUE if the covariance model allows for a ‘standard’ definition (see convert.to.readable and CovarianceFct) and transform is NULL. If a ‘standard’ definition is given and both the variance and the nugget are either not estimated or do not appear on the right hand side of the transform, then standard.style might be set to TRUE by the user. This accelerates the MLE algorithm. The responsibility is completely left to the user, then.
lsq.methods | variants of the least squares fit of the variogram. See Details.
mle.methods | variants of the maximum likelihood fit of the covariance function. See Details.
cross.methods | Not implemented yet.
users.guess | User's guess of the parameters. All the parameters must be given using the same rules as for either param (except that no NA's should be contained) or model.
users.beta | User's guess of the trend.
table.format | logical. The parameters controls the output. See Value.

Details

The optimisations are performed using optimize if one parameter has to be estimated only and optim, otherwise.

First, by means of various control parameters, see below, the algorithm first tries to estimate the bounds for the parameters to be estimated, if the bounds for the parameters are not given.
Independently whether users.guess is given, the algorithm guesses initial values for the parameters. The automatic guess and the user's guess will be called primitive methods in the following.

Second, the variogram model is fitted by various least squares methods (according to the value of lsq.methods) using the best parameter set among the primitive methods as initial value if the effective number of parameters is greater than 1.

[Remarks: (i) “best” with respect to the target value of the respective lsq method; (ii) the effective number of parameters in the optimisation algorithm can be smaller than the number of estimated parameters, since in some cases, some parameters can be calculated explicitely; relevant for the choice between optimize and optim is the effective number of parameters; (iii) optim needs]

Third, the model is fitted by various maximum likelihood methods (according to the value of mle.methods) using the best parameter set among the primitive methods and the lsq methods as initial value (if the effective number of parameters is greater than 1).

Comments on specific parameters:

  • lower
    The lower bounds are technical bounds that should not restrict really restrict the domaine of the value. However, if these values are too small the optimisation algorithm will frequently run into local minima or get stuck close the border of the parameter domain. It is advised to limit seriously the domain of the additional parameters of the covariance model and/or the total number of parameters to be estimated, if “many” parameters of the covariance model are estimated.

    If the model is given in standard form, the user may supply the lower bounds for the whole parameter vector, or only for the additional form parameters of the model. The lower bound for the mean will be ignored. lower may contain NAs, then these values are generated by the

    If a nested model is given, the bounds may again be supplied for all parameters or only for the additional form parameters of the model. The bounds given apply uniformely to all submodels of the nested model.

    If the model is given in list format, then lower is a list, where components may be missing or NA. These are generated by the algorithm, then.

    If lower is NULL all lower bounds are generated automatically.

  • upper.kappa
    See lower.kappa.
  • sill
    Additionally to estimating nugget and variance separately, they may also be estimated together under the condition that nugget + variance = sill. For the latter a finite value for sill has to be supplied, and nugget and variance are set to NA.

    sill is only used for the standard model.

  • use.naturalscaling
    logical. If TRUE then internally, rescaled covariance functions will be used for which cov(1)~=0.05. However this parameter does not influence the output of fitvario: the parameter vector returned by fitvario refers always to the standard covariance model as given in CovarianceFct. (In contrast to PracticalRange in RFparameters.)
    Advantages if use.naturalscaling=TRUE:
    • scale and the shape parameter of a parameterised covariance model can be estimated better if they are estimated simultaneously.
    • The estimated bounds calculated by means of upperbound.scale.factor and lowerbound.scale.factor, etc. might be more realistic.
    • in case of anisotropic models, the inverse of the elements of the anisotropy matrix should be in the above bounds.

    Disadvantages if use.naturalscaling=TRUE:

    • For some covariance models with additional parameters, the rescaling factor has to be determined numerically. Then, more time is needed to perform fitvario.

    Default: TRUE.

  • PrintLevel
    0 : no message
    1 : error messages
    2 : warnings
    3 : minimum debugging information
    5 : extended debugging information, including graphics
    Default: 0.
  • trace.optim
    see control parameter trace of optim. Default: 0.
  • bins
    vector of explicit boundaries for the bins or the number of bins for the empirical variogram (used in the LSQ target function, which is described at the beginning of the Details). Note that for anisotropic models, the value of bins might be enlarged. Default: 20.
  • distance.factor
    right boundary of the last bin is calculated as distance.factor * (maximum distance between all pairs of points). Only used if bins is a scalar. Default: 0.5.
  • upperbound.scale.factor
    The upper bound for the scale is determined as upperbound.scale.factor * (maximum distance between all pairs of points). Default: 10.
  • lowerbound.scale.factor
    The lower bound for the scale is determined as

    (minimum distance between different pairs of points) / lowerbound.scale.factor.

    Default: 20.

  • lowerbound.scale.LS.factor
    For the LSQ target function a different lower bound for the scale is used. It is determined as

    (minimum distance between different pairs of points) / lowerbound.scale.LS.factor.

    Default: 5.

  • upperbound.var.factor
    The upper bound for the variance and the nugget is determined as

    upperbound.var.factor * var(data).

    Default: 10.

  • lowerbound.var.factor
    The lower bound for the variance and the nugget is determined as

    var(data) / lowerbound.var.factor.

    If a standard model definition is given and either the nugget or the variance is fixed, the parameter to be estimated must also be greater than lowerbound.sill. If a non-standard model definition is given then lowerbound.var.factor is only used for the first model; the other lower bounds for the variance are zero. Default: 100.

  • lowerbound.sill
    See lowerbound.var.factor. Default: 1E-10.
  • scale.max.relative.factor
    If the initial scale value for the ML estimation obtained by the LSQ target function is less than

    (minimum distance between different pairs of points) / scale.max.relative.factor

    a warning is given that probably a nugget effect is present. Note: if scale.max.relative.factor is greater than lowerbound.scale.LS.factor then no warning is given as the scale has the lower bound (minimum distance between different pairs of points) / lowerbound.scale.LS.factor. Default: 1000.

  • minbounddistance
    If any value of the parameter vector returned from the ML estimation is closer than minbounddistance to any of the bounds or if any value has a relative distance smaller than minboundreldist, then it is assumed that the MLE algorithm has dropped into a local minimum, and it will be continued with evaluating the ML target function on a grid, cf. the beginning paragraphs of the Details. Default: 0.001.
  • minboundreldist
    See minbounddistance. Default: 0.02.
  • approximate.functioncalls
    In case the parameter vector is too close to the given bounds, the ML target function is evaluated on a grid to get a new initial value for the ML estimation. The number of points of the grid is approximately approximate.functioncalls. Default: 50.
  • lsq.methods
    Variogram fit by least squares methods; first, a preliminary trend is estimated by a simple regression; second, the variogram is fitted; third, the trend is fitted using the estimated covariance structure.
    • "self" weighted lsq. Weights are the values of the fitted variogram to the power of -2
    • "plain" model fitted by least squares; trends are never taken into account
    • "sqrt.nr" weighted lsq. Weight is the square root of the number of points in the bin
    • "sd.inv" weighted lsq. Weight is the inverse the standard deviation of the variogram cloud within the bin
  • mle.methods
    Model fit by various maximum likelihood methods (according to the value of mle.methods) using the best parameter set among the primitive methods and the lsq methods as initial value (if the effective number of parameters is greater than 1). If the best parameter vector of the MLE found so far is too close to some given bounds, see the specific parameters above, it is assumed that optim ran into a local minimum because of a bad starting value. In this case and if refine.onborder=TRUE the MLE target function is calculated on a grid, the best parameter vector is taken, and the optimisation is restarted with this parameter vector.

    • "ml" maximum likelihood; since ML and REML give the same result if there are not any covariates, ML is performed in that case, independently whether it is given or not.
    • "reml" restricted maximum likelihood

Value

If table.format=TRUE then a matrix is returned where the first rows contain the estimated parameters, followed by the target values of all methods for the given set of parameters; the last rows give the lower and upper bounds used in the estimations. The columns correspond to the various estimation methods for the parameters.
If table.format=FALSE then the function returns a list with the following elements

ev | list returned by EmpiricalVariogram
variogram | list of fitted models; they include the fixed and the estimated parameters; each list element gives the estimated covariance model
values | list of optimal target values

Acknowledgement

Thanks to Paulo Ribeiro for hints and comparing the perliminary versions of fitvario in RandomFields V1.0 to likfit of the package geoR whose homepage is at http://www.est.ufpr.br/geoR/.

Note

This function does not depend on the value of RFparameters()$PracticalRange. The function fitvario always uses the standard specification of the covariance model as given in CovarianceFct.

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute

References

  • Least squares and mle methods

    Cressie, N.A.C. (1993) Statistics for Spatial Data. New York: Wiley.

  • Related software
    Ribeiro, P. and Diggle, P. (2001) Software for geostatistical analysis using R and S-PLUS: geoR and geoS, version 0.6.15. http://www.maths.lancs.ac.uk/~ribeiro/geoR.html.

See Also

CovarianceFct, GetPracticalRange, parampositions RandomFields,

Examples

 model <-"gencauchy"
 param <- c(0, 1, 0, 1, 1, 2)
 estparam <- c(0, NA, 0, NA, NA, 2) ## NA means: "to be estimated"
 ## sequence in `estparam' is
 ## mean, variance, nugget, scale, (+ further model parameters)
 ## So, mean, variance, and scale will be estimated here.
 ## Nugget is fixed and equals zero.
 points <- 100
 x <- runif(points,0,3)
 y <- runif(points,0,3) ## 100 random points in square [0, 3]^2
 d <- GaussRF(x=x, y=y, grid=FALSE, model=model, param=param, n=10)
 str(fitvario(x=cbind(x,y), data=d, model=model, param=estparam,
              lower=c(0.1, 0.1), upper=c(1.9, 5), cross.me=NULL))

#########################################################
## The next two estimations give about the same result.
## For the first the sill is fixed to 1.5. For the second the sill
## is reached if the estimated variance is smaller than 1.5
estparam <-  c(0, NA, NA, NA, NA, NA)
str(fitvario(x=cbind(x,y), data=d, model=model, param=estparam,
    sill=1.5, cross.me=NULL))

estparam <-  c(0, NA, NaN, NA, NA, NA)
parampositions(model=model, param=estparam)
f <- function(param) {
   param[5] <- max(0, 1.5 - param[1])
   return(param)
}
str(fitvario(x=cbind(x,y), data=d, model=model, param=estparam,
    transform=f, cross.me=NULL))

#########################################################
## esimation of coupled parameters (kappa1 = kappa2, here)
estmodel <- list(list(model="gencauchy", var=NA, scale=NA,
                      kappa=list(NA, function(m) m[[1]]$kappa[1])),
                 "+",
                 list(model="nugget", var=NA))
system.time(str(fitvario(x=cbind(x,y), data=d, model=estmodel,
              cross.me=NULL), vec.len=6))

## much faster implementation, but trickier coding
estparam <- c(0, NA, NA, NA, NA, NaN)
parampositions(model=model, param=estparam)
f <- function(param) {param[3] <- param[2]; param}
system.time(str(fitvario(x=cbind(x,y), data=d, model=model,
      param=estparam, transform=f, cross.me=NULL), vec.len=6))

## about as fast as previous coding
## about the same precision as previous coding, or better
## But, a warning is given, since the user may programme
## strange things in this setup, and the programme cannot check it.
system.time(str(fitvario(x=cbind(x,y), data=d, model=model,
      param=estparam, transform=f, cross.me=NULL,
      standard.style=TRUE), vec.len=6))

#########################################################
## estimation in a anisotropic framework
  x <- y <- (1:6)/4
  model <- list(list(m="exp", v=1.5, aniso=c(4,2,-2,1)))
  z <- GaussRF(x=x, y=y, grid=TRUE, model=model, n=10)
  estmodel <- list(list(m="exp", v=NA, aniso=c(NA, NA,-2,1)))
  system.time(str(fitvario(as.matrix(expand.grid(x, y)), data=z,
                           model=estmodel, nphi=20)))


fractal dimension

Description

The function estimates the fractal dimension of a process

Usage

fractal.dim(x, y = NULL, z = NULL, data,
           grid=TRUE, gridtriple = FALSE,
           bin= seq(min(ct$x[3, ]) / 2,
             min((ct$x[2,]-ct$x[1,]) / 4, vario.n * min(ct$x[3,]) + 1),
             min(ct$x[3,])),
           vario.n=5,
           sort=TRUE,
           fft.m = c(65, 86), ## in
           fft.max.length=Inf,
           fft.max.regr=150000,
           fft.shift = 50, # in
           method=c("variogram", "fft"),
           mode=c("plot", "interactive"),
           pch=16, cex=0.2, cex.main=0.85,
           PrintLevel = RFparameters()$Print,
           height=3.5,
           ...)

Arguments

x | matrix of coordinates, or vector of x coordinates; if x is not given a grid with unit grid length is assumed
y | vector of y coordinates
z | vector of z coordinates
data | the values measured.
grid | determines whether the vectors x, y, and z should be interpreted as a grid definition, see Details. grid does not apply for T.
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.
bin | sequence of bin boundaries for the empirical variogram
vario.n | first vario.n value of the empirical variogram are used for the regression fit that are not NA.
sort | If TRUE then the coordinates are permuted such that the largest grid length is in x-direction; this is of interest for algorithms that slice higher dimensional fields into one-dimensional sections.
fft.m | numeric vector of two components; interval of frequencies for which the regression should be calculated; the interval is given in percent of the range of the frequencies in log scale.
fft.max.length | The first dimension of the data is cut into pieces of length fft.max.length. For each piece the FFT is calculated and then the average for all pieces is taken. The pieces may overlap, see the parameter fft.shift.
fft.max.regr | If the fft.m is too large, parts of the regression fit will take a very long time. Therefore, the regression fit is calculated only if the number points given by fft.m is less than fft.max.regr.
fft.shift | This parameter is given in percent [of fft.max.length] and defines the overlap of the pieces defined by fft.max.length. If fft.shift=50 the WOSA estimator is given; if fft.shift=100 no overlap exist.
method | list of implemented methods to calculate the fractal dimension; see Details
mode character. A vector with components 'nographics', 'plot', or 'interactive':
'nographics'
no graphical output
'plot'
the regression line is plotted
'interactive'
the regression domain can be chosen interactively
Usually only one mode is given. Two modes may make sense in the combination c("plot", "interactive") in which case all the results are plotted first, and then the interactive mode is called. In the interactive mode, the regression domain is chosen by two mouse clicks with the left mouse; a right mouse click leaves the plot.
pch | vector or scalar; sign by which data are plotted.
cex | vector or scalar; size of pch.
cex.main | The size of the title in the regression plots.
PrintLevel | integer. If PrintLevel is 0 nothing is printed. If PrintLevel=1 error messages are printed. If PrintLevel=2 warnings and the regression results are given. If PrintLevel>2 tracing information is given.
height | height of the grahics window
... graphical parameters

Details

The function calculates the fractal dimension by various methods:

  • variogram method
  • Fourier transform

Value

The function returns a list with elements vario, fft corresponding to the 2 methods given in the Details.
Each of the elements is itself a list that contains the following elements.

x | the x-coordinates used for the regression fit
y | the y-coordinates used for the regression fit
regr | the return list of the lm
sm | smoothed curve through the (x,y) points
x.u | NULL or the restricted x-coordinates given by the user in the interactive plot
y.u | NULL or y-coordinates according to x.u
regr.u | NULL or the return list of lm for x.u and y.u
D | the fractal dimension
D.u | NULL or the fractal dimension corresponding to the user's regression line

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute

References

variogram method

  • Constantine, A.G. and Hall, P. (1994) Characterizing surface smoothness via estimation of effective fractal dimension. J. R. Statist. Soc. Ser. B 56, 97-113.

fft

  • Chan, Hall and Poskitt (1995)

See Also

CovarianceFct, hurst


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

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 | 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.

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.

  • An isotropic random field is created by 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, ...).
  • Nested models can be defined in the same way as a nested CovarianceFct. If the trend is not given it is set to 0.
  • An anisotropic random field (i.e. zonal anisotropy, geometrical anisotropy, separable models, non-separable space-time models) and a random field based on multiplicative or nested models is defined as in the case of an anisotropic 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.
  • The strings for model and method may be abbreviated as long as the abbreviations match only one option. See also PrintModelList() and PrintMethodList()
  • Further control parameters for the simulation are set by means of RFparameters(...).

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

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/institute

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

References

See RFMethods for the references.

See Also

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,
              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)
  }
}


Get input behaviour

Description

The functions return values stored by useraction, Locator and Readline

Usage

  getactions()
  getactionlist()

Value

getactionlist returns a list of stored values created by Locator and Readline
getactions returns a list of all parameters that influence the behaviour of Locator and Readline.

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute

See Also

Locator, Readline, useraction,

Examples

## see useraction


Determination of the practical range

Description

the function returns the practical range of a covariance model, i.e. the distance for which the model with standard parameters (variance=1, nugget=0, scale=1) values 0.05

Usage

GetPracticalRange(model, kappas=NULL, dim=1)

Arguments

model | string; covariance model, see CovarianceFct, or type PrintModelList() to get all options.
kappas | additional parameters of a parametrised model, see CovarianceFct.
dim | space-time dimension; except for some few anisotropic models, this parameter does not have an effect

Details

In case the practical range has to be determined numerically, only a rough approximation is provided.

The function can only be applied to isotropic models

Value

real number (the practical range)

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute

References

Goovaerts, P. (1997) Geostatistics for Natural Resources Evaluation. New York: Oxford University Press.

See Also

CovarianceFct, RandomFields, RFparameters,

Examples

GetPracticalRange("exponential")
GetPracticalRange("whittlematern",kappa=2)
GetPracticalRange("gengneiting",kappa=c(a=3,b=5))


Internal information

Description

The function returns internal information about the simulation of a random field

Usage

GetRegisterInfo(register=0, ignore.active = FALSE, max.elements=10^6)
GetModelInfo(register=0)

Arguments

register | -2, -1, 0:9; place where intermediate calculations are stored; the numbers 0:9 are aliases for 10 internal registers, see also GaussRF
ignore.active | logical. If FALSE and the register has non-active flag (because of an error or it is deleted) then a list is returned that contains only the element active=FALSE. Otherwise, the full list is always returned.
max.elements | integer; since GetRegisterInfo might be a dump of simulation that needs a huge amount of memory and since all entries are copied, the maximal amount of available memory might be easily exceeded. Therefore, only the size of the critical parts are returned and not the vector or matrix itself, if the number of elements exceeds max.elements.

Details

GetRegisterInfo(register, ignore.active=TRUE) is useful for debugging and specialists' need to control the algorithm, see the examples in RFparameters and GaussRF.

If RFparameters()$Storage=FALSE then values of the internal registers are not kept if GaussRF or DoSimulateRF has been called. Hence GetRegisterInfo cannot provide any information.

GetModelInfo returns the partial information on the internal model structure. It allows for register=-1 returning the model structure for the last use of CovarianceFct or Variogram or similar commands. register=-2 is for internal use only.

Value

List of internal information is returned.

Note

Put Storing=TRUE, see RFparameters if you like to have more internal information in case of an expected failure of an initialisation of a random field simulation.

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute

See Also

GaussRF

Examples

# lets regard at a non-initialised register first:
DeleteAllRegisters()
str(GetRegisterInfo(0))

# now, a register that has successfully initialised
try(GaussRF(1:4, grid=TRUE, model="exp", param=c(1,2,3,4), Storing=TRUE))
str(GetRegisterInfo(0))

# finally, a register that could not be successfully initialised
try(GaussRF(runif(4), grid=FALSE, model="exp", param=c(1,2,3,4),
            me="ci", Storing=TRUE))
str(GetRegisterInfo(0, TRUE))


System calls

Description

The functions hostname and pid return the host name and the PID, respectively.

Usage

hostname()

pid()

Details

If R runs on a unix platform the host name and the PID are returned, otherwise the empty string and naught, respectively.

Value

hostname | returns a string
pid | returns an unsigned integer

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute


Hurst coefficient

Description

The function estimates the Hurst coefficient of a process

Usage

hurst(x, y = NULL, z = NULL, data,
      gridtriple = FALSE, sort = TRUE,
      block.sequ = unique(round(exp(seq(log(3000), log(dim[1]),
                          len=min(100, dim[1]))))),
      fft.m = c(1, min(1000, (fft.len - 1) / 10)),
      fft.max.length = Inf, method=c("dfa", "fft", "var"),
      mode=c("plot", "interactive"), pch=16, cex=0.2, cex.main=0.85,
      PrintLevel=RFparameters()$Print,height=3.5,...)

Arguments

x | matrix of coordinates, or vector of x coordinates
y | vector of y coordinates
z | vector of z coordinates
data | the data
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
sort | logical. If TRUE then the coordinates are permuted such that the largest grid length is in x-direction; this is of interest for algorithms that slice higher dimensional fields into one-dimensional sections.
block.sequ | ascending sequences of block lengths for which the detrended fluctuation analysis and the variance method is performed.
fft.m | vector of 2 integers; lower and upper endpoint of indices for the frequency which are used in the calculation of the regression line for the periodogram near the origin.
fft.max.length | if the number of points in x-direction is larger than fft.max.length then the segments of length fft.max.length are considered, shifted by fft.max.length/2 (WOSA-estimator).
method | list of implemented methods to calculate the Hurst parameter; see Details
mode character. A vector with components 'nographics', 'plot', or 'interactive':
'nographics'
no graphical output
'plot'
the regression line is plotted
'interactive'
the regression domain can be chosen interactively
Usually only one mode is given. Two modes may make sense in the combination c("plot", "interactive") in which case all the results are plotted first, and then the interactive mode is called. In the interactive mode, the regression domain is chosen by two mouse clicks with the left mouse; a right mouse click leaves the plot.
pch | vector or scalar; sign by which data are plotted.
cex | vector or scalar; size of pch.
cex.main | font size for title in regression plot, see regression; only used if mode includes 'plot' or 'interactive'
PrintLevel | integer. If PrintLevel is 0 or 1 nothing is printed. If PrintLevel=2 warnings and the regression results are given. If PrintLevel>2 tracing information is given.
height | height of the graphics window
... graphical parameters

Details

The function is still in development. Several functionalities do not exist - see the code itself for the current stage.

The function calculates the Hurst coefficient by various methods:

  • detrended fluctuation analysis (dfa)
  • aggregated variation (var)
  • periodogram or WOSA estimator (fft)

Value

The function returns a list with elements dfa, varmeth, fft corresponding to the three methods given in the Details.
Each of the elements is itself a list that contains the following elements.

x | the x-coordinates used for the regression fit
y | the y-coordinates used for the regression fit
regr | the coefficients of the lm
sm | smoothed curve through the (x,y) points
x.u | NULL or the restricted x-coordinates given by the user in the interactive plot
y.u | NULL or y-coordinates according to x.u
regr.u | NULL or the coefficients of lm for x.u and y.u
H | the Hurst coefficient
H.u | NULL or the Hurst coefficient corresponding to the user's regression line

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute

References

detrended fluctuation analysis

  • Peng, C.K., Buldyrev, S.V., Havlin, S., Simons, M., Stanley, H.E. and Goldberger, A.L. (1994) Mosaic organization of DNA nucleotides Phys. Rev. E 49, 1685-1689

aggregated variation

  • Taqqu, M.S. and Teverovsky, V. (1998) On estimating the intensity of long range dependence in finite and infinite variance time series. In: Adler, R.J., Feldman, R.E., and Taqqu, M.S. A Practical Guide to Heavy Tails, Statistical Techniques an Applications. Boston: Birkhaeuser
  • Taqqu, M.S. and Teverovsky, V. and Willinger, W. (1995) Estimators for long-range dependence: an empirical study. Fractals 3, 785-798

periodogram

  • Percival, D.B. and Walden, A.T. (1993) Spectral Analysis for Physical Applications: Multitaper and Conventional Univariate Techniques, Cambridge: Cambridge University Press.
  • Welch, P.D. (1967) The use of {F}ast {F}ourier {T}ransform for the estimation of power spectra: a method based on time averaging over short, modified periodograms IEEE Trans. Audio Electroacoustics 15, 70-73.

See Also

CovarianceFct, fractal.dim


Kriging methods

Description

The function allows for different methods of kriging.

Usage

Kriging(krige.method, x, y=NULL, z=NULL, T=NULL, grid,
        gridtriple=FALSE, model, param, given, data, trend, pch=".",
        return.variance=FALSE, internal=FALSE)

Arguments

krige.method | kriging method; currently only 'S' (simple kriging) and 'O' (ordinary kriging) implemented.
x | (n x d) matrix or vector of x coordinates; coordinates of n points to be kriged
y | vector of y coordinates.
z | vector of z coordinates.
T | vector in grid triple form for the time coordinates.
grid | logical; determines whether the vectors x, y, and z should be interpreted as a grid definition, see Details.
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.
model | string; covariance 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 covariance functions, see CovarianceFct. The value of mean must be finite in the case of simple kriging, and is ignored otherwise.
given | matrix or vector of points where data are available.
data | the data values given at given; it might be a vector or a matrix. If a matrix is given multivariate data are assumed which are kriged separately.
trend | not programmed yet (will be used in case of universal kriging)
pch | Kriging procedures are quite time consuming in general. The character pch is printed after roughly each 80th part of calculation.
return.variance | logical. If FALSE the kriged field is returned. If TRUE a list of two elements, estim and var, i.e. the kriged field and the kriging variances, is returned.
internal | FALSE. internal should not be set to TRUE by the user. (In case internal=TRUE, various consistency checks for the input variables are not performed. Further, grid must be FALSE, and model must be given in the output format of PrepareModel.)

Details

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

Value

If variance.return=FALSE Kriging returns a vector or matrix of kriged values corresponding to the specification of x, y, z, and grid, and data.

data: a vector or matrix with one column
* 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 (according to the specification of x, y, and z).

data: a matrix with at least two columns
* grid=FALSE. A matrix with the ncol(data) columns is returned.
* grid=TRUE. An array of dimension d+1, where d is the dimension of the random field, is returned (according to the specification of x, y, and z). The last dimension contains the realisations.
If variance.return=TRUE a list of two elements, estim and var, i.e. the kriged field and the kriging variances, is returned. The format of estim is the same as described above. The format of var is accordingly.

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute

References

Chiles, J.-P. and Delfiner, P. (1999) Geostatistics. Modeling Spatial Uncertainty. New York: Wiley.

Cressie, N.A.C. (1993) Statistics for Spatial Data. New York: Wiley.

Goovaerts, P. (1997) Geostatistics for Natural Resources Evaluation. New York: Oxford University Press.

Wackernagel, H. (1998) Multivariate Geostatistics. Berlin: Springer, 2nd edition.

See Also

CondSimu, CovarianceFct, EmpiricalVariogram, RandomFields,

Examples

## creating random variables first
## here, a grid is chosen, but does not matter
step <- 0.25
x <-  seq(0,7,step)
param <- c(0,1,0,1)
model <- "exponential"
RFparameters(PracticalRange=FALSE)
p <- 1:7
points <- as.matrix(expand.grid(p,p))
data <- GaussRF(points, grid=FALSE, model=model, param=param)

## visualise generated spatial data
zlim <- c(-2.6,2.6)
colour <- rainbow(100)
image(p, p, xlim=range(x), ylim=range(x),
      matrix(data,ncol=length(p)),
      col=colour,zlim=zlim)

## now: kriging
krige.method <- "O" ## ordinary kriging
z <-  Kriging(krige.method=krige.method,
              x=x, y=x, grid=TRUE,
              model=model, param=param,
              given=points, data=data)
image(x,x,z,col=colour,zlim=zlim)


Graphical Input

Description

Reads the position of the graphics cursor when the (first) mouse button is pressed, or replays from storage.

Usage

Locator(n, type="n", info=NULL, ...)

Arguments

n | the maximum number of points to locate; the value of n is ignored if
type | One of '"n"', '"p"', '"l"' or '"o"'. If '"p"' or '"o"' the points are plotted; if '"l"' or '"o"' they are joined by lines.
info | arbitrary object; tracing information that is useful in case the user likes to edit the stored sequence of inputs; the value of info is considered only if useraction has been called by action="start.register" or action="continue.register"
... | additional graphics parameters

Details

The behaviour of Locator depends on the the value of action set by useraction, see there for more information.

See also locator for further information on the parameters.

Value

A list containing 'x' and 'y' components which are the coordinates of the identified points in the user coordinate system, i.e., the one specified by 'par("usr")'.

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute

See Also

getactions, Readline, useraction, userinput

Examples

## see useraction


Max-Stable Random Fields

Description

These functions simulate stationary and isotropic max-stable random fields with unit Frechet margins.

Usage

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

InitMaxStableRF(x, y=NULL, z=NULL, grid, model, param, maxstable,
               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; see CovarianceFct, or type PrintModelList() to get all options; interpretation depends on the value of maxstable, see Details.
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 covariance functions, see CovarianceFct, or be given in one of the extended forms, see Details
maxstable | string. Either 'extremalGauss' or 'BooleanFunction'; see Details.
method | NULL or string; method used for simulating, see RFMethods, or type PrintMethodList() to get all options; interpretation depends on the value of maxstable.
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; if gridtriple=FALSE ascending sequences for the parameters x, y, and z are expected; if gridtriple=TRUE triples of form c(start,end,step) expected; this parameter is used only if grid=TRUE
... | RFparameters that are locally used only.

Details

There are two different kinds of models for max-stable processes implemented:

  • maxstable="extremalGauss"
    Gaussian random fields are multiplied by independent random factors, and the maximum is taken. The random factors are such that the resulting random field has unit Frechet margins; the specification of the random factor is uniquely given by the specification of the random field. The parameter vector param, the model, and the method are interpreted in the same way as for Gaussian random fields, see GaussRF.
  • maxstable="BooleanFunction"
    Deterministic or random, upper semi-continuous L1-functions are randomly centred and multiplied by suitable, independent random factors; the pointwise maximum over all these functions yields a max-stable random field. The simulation technique is related to the random coin method for Gaussian random field simulation, see RFMethods. Hence, only models that are suitable for the random coin method are suitable for this technique, see PrintModelList() for a complete list of suitable covariance models.
    The only value allowed for method is 'max.MPP' (and NULL), see PrintMethodList(). In the parameter list param the first two entries, namely mean and variance, are ignored. If the nugget is positive, for each point an additional independent unit Frechet variable with scale parameter nugget is involved when building the maximum over all functions.

    The model may be defined alternatively in one of the two extended ways as introduced in CovarianceFct and GaussRF. However only a single model may be given! The model may be anisotropic.

Value

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

MaxStableRF and DoSimulateRF return NULL if an error has occurred; otherwise the returned object depends on the parameters:
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.

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute

References

Schlather, M. (2002) Models for stationary max-stable random fields. Extremes 5, 33-44.

See Also

CovarianceFct, GaussRF, RandomFields, RFMethods, RFparameters, DoSimulateRF, .

Examples

n <- 30 ## nicer, but time consuming if n <- 100
x <- y <- 1:n
ms0 <- MaxStableRF(x, y, grid=TRUE, model="exponen",
                 param=c(0,1,0,40), maxstable="extr"
                 ,CE.force = TRUE                 )
image(x,y,ms0)

############################################################
##                                                        ##
##  Plots used in M. Schlather, Extremes 5:1, 33-44, 2002 ##
##                                                        ##
############################################################

pts <- if (interactive()) 512 else 32
x <- (1:pts) / pts * 10
scalegauss <- 1.5
RFparameters(MPP.radius=2*scalegauss, Print=3)
runif(1)
save.seed <- .Random.seed
ms1 <- MaxStableRF(x, x, model="gauss", param=c(0, 1, 0, scalegauss),
                  maxstable="Bool", grid=TRUE)
image(x, x, sqrt(ms1))

scalecone <- 3.7
.Random.seed <- save.seed
ms2 <- MaxStableRF(x, x, model="cone", param=c(0, 1, 0, scalecone,
           0, 0, 1), maxstable="Bool", grid=TRUE)
image(x, x, sqrt(ms2))

.Random.seed <- save.seed
ms3 <- MaxStableRF(x, x, model="exponen", param=c(0, 1, 0, 1),
           maxstable="extr", grid=TRUE, method="ci")
image(x, x, sqrt(ms3))

.Random.seed <- save.seed
ms4 <- MaxStableRF(x, x, model="gauss", param=c(0, 1, 0, 1),
           maxstable="extr", grid=TRUE, method="ci")
image(x, x, sqrt(ms4))


Range of model specific parameters

Description

parameter.range returns the range of the parameters specific to a parametrised variogram model or a covariance function

Usage

parameter.range(model, dim)

Arguments

model | the variogram model or covariance function, see CovarianceFct
dim | the dimension of the random field

Value

If the model is not valid or not specified for the given dimension, NaN is returned.
If the model does not have additional parameters, NULL is returned. Otherwise a list is returned.
If the parameter space is simple, a list of the following elements is returned

theoretical | The theoretical ranges of the parameters. Note that open and closed intervals are not distinguished here.
practical | The ranges of the parameters usually not exceeded in practical applications.

Both elements are matrices with two rows and the number of columns being the number of parameters.

If the parameter space is complex, the parameter space is devided into rectangular areas and theoretical and practical are lists where each element is a matrix as above.

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute

See Also

CovarianceFct

Examples

  parameter.range("exponential", 1) # NULL
  str(parameter.range("power", 1))
  str(parameter.range("power", 2))
  str(parameter.range("gengneiting", 2))
  str(parameter.range("nsst",1))
  str(parameter.range("nsst2",1))
  str(parameter.range("hyper",1))


Position of the parameters

Description

The function returns the internal positions of the model parameters

Usage

  parampositions(model, param, print=TRUE)

Arguments

model | see CovarianceFct
param | see CovarianceFct
print | if FALSE only an invisible list is returned

Value

The model is printed and returned where the values of the parameters are the positions in the internal representation.
An error appears if the model can be substantially simplified.

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute

See Also

CovarianceFct, fitvario

Examples

## compare the output of the following commands
parampositions(model="exp", param=c(0,1,0,NA))
parampositions(model="exp", param=c(0,1,NA,NA))
parampositions(model="exp", param=c(0,0,1,NA))
## that is, the nugget in the standard model is removed if naught!
## the values of the other parameters do not matter. (First value
## of the returned vector refers to the mean position.)

parampositions(model="whi", param=cbind(c(1, 1, 1), c(2, 2, 2)))
parampositions(model="whi", param=cbind(c(1, 1, 1), c(2, 0, 2)))
## second lines, second column defines a nugget effect since scale is 0!

try(parampositions(model="whi", param=cbind(c(1, 1, 1), c(0, 0, 0))))
## leads to an error

try(parampositions(model="whi",
                   param=cbind(c(1, 1, 1), c(2, 0, 0), c(1, 0, 0))))
## leads to an error

try(parampositions(model="whi",
                   param=cbind(c(1, 1, 1), c(NA, 0, 0), c(1, 0, 0))))
## leads to an error

parampositions(model=list(list(model="exp", var=3, scale=6), "+",
                           list(model="whittle", var=2, scale=7,
                                kappa=NA))) ## again the values do not matter


Simulation and Analysis of Random Fields

Description

The package RandomFields allows for simulating various kinds of random fields, including anisotropic processes. Furthermore, algorithms for conditional simulation and simulation of max-stable random fields are provided.

Additionally, the package includes tools for analysing spatial data: Hurst parameter, fractal dimension, empirical variogram, interactive fitting of parameters, LSQ and MLE estimation of parameters. Basic kriging procedures are also provided.

Further extensions of the package are planned.

Details

The following random fields and related functionalities are provided by the package.

  1. stationary and isotropic Gaussian random fields
    • CondSimu : conditional simulation
    • CovarianceFct : covariance functions or variogram models
    • EmpiricalVariogram : empirical variogram
    • GaussRF : simulation of Gaussian random fields; nice examples to get familiar with the simulation features of the package
    • Kriging : simple and ordinary kriging
    • fitvario : variogram/covariance function fit by least squares, maximum likelihood and cross validation techniques
    • PrintMethodList : list of implemented simulation methods
    • ShowModels : interactive, graphical choice of models
    • soil : Soil physical and chemical data; the example gives a simple geostatistical analysis using features of the package
  2. stationary (and isotropic) max-stable random fields
    • CovarianceFct : covariance models for extremal Gaussian random fields
    • MaxStableRF : simulation of max-stable random fields
  3. stationary and isotropic Poisson random fields (not implemented yet)

Functions used in diverse simulation methods:

  • DeleteRegister : deleting internal registers
  • RFparameters : control parameters (advanced settings)

Functions of general purpose:

  • eval.parameters : provides an interactive menu
  • fractal.dim : estimation of the fractal dimension
  • hurst : estimation of the Hurst parameter
  • regression : interactive regression plot
  • Locator, Readline : locator and readline, respectively, with storage and replay functionality
  • sleep.milli : sets the process into sleeping status

Acknowledgement

Many thanks to Martin Maechler, Paulo Ribeiro, and Tilmann Gneiting for proof-reading parts of the code and the help text for the versions 1.0 of this package. Yindeng Jiang jiangyindeng@gmail.com implemented the circulant embedding methods ‘cutoff’ and ‘intrinsic’ in 2004 for the versions 1.2.

The coding of V1.0 has been financially supported by the EU TMR network ERB-FMRX-CT96-0095 on ``Computational and statistical methods for the analysis of spatial data'' in 1999, and by the German Federal Ministry of Research and Technology (BMFT) grant PT BEO 51-0339476C during 2000-03.

Note

From version 1.3 on: Default value of RFparameters()$Storing is FALSE.

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute;

References

This package is announced in and can be cited by:

Schlather, M. (2001) Simulation of stationary and isotropic random fields. R-News 1 (2), 18-20.

Schlather, M. (2004) Simulation of random fields and applications. In preparation.


Read a Line

Description

Readline reads a line from the terminal or from storage

Usage

   Readline(prompt="", info=NULL)

Arguments

prompt | the string printed when prompting the user for input. Should usually end with a space " "
info | arbitrary object; tracing information that is useful in case the user likes to edit the stored sequence of inputs; the value of info is considered only if useraction has been called by action="start.register" or action="continue.register"

Details

The behaviour of Readline depends on the the value of action set by useraction, see there for more information.

The prompt string will be truncated to a maximum allowed length, normally 256 chars (but can be changed in the source code), since the function is based on readline.

Value

A character vector of length one.

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute

See Also

getactions, Locator, readline useraction, userinput,

Examples

## see useraction


Deleting Intermediate Results

Description

DeleteRegister deletes an internal register of the random field simulation

DeleteAllRegisters deletes all internal registers

Usage

DeleteRegister(nr=0)
DeleteAllRegisters()

Arguments

nr number of the register to be deleted; in 0:9

Details

DeleteRegister and DeleteAllRegister should be used if

  1. the simulations are finished and there is not enough memory available to continue the R session, and RFparameters()$Storing has been TRUE.
  2. one of the parameters RFparameters()$Storing or RFparameters()$PracticalRange ((or RFparameters()$TBM.method, once this parameter is available)) has been changed by the user, and if these options should become global.

See GaussRF and RFparameters for further details.

Value

NULL.

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de, http://www.stochastik.math.uni-goettingen.de/institute

See Also

GaussRF, RandomFields and RFparameters.

Examples

 DeleteAllRegisters() ## no effect visible for the user


Regression plot

Description

Interactive regression plot

Usage

regression(x, y, main, scr, value=function(regr) regr$coeff[2],
           variable="var", mode=c("nographics", "plot", "interactive"),
           regr.select=c("coefficients"), averaging=FALSE, cex.main=0.85,
           col.passive = "grey", col.active = "green",
           col.main.passive = "black", col.main.active = "red",
           col.points.chosen = "lightblue", col.abline = "yellow",
           col.abline.chosen = "darkblue", col.smooth = "red", ...)

Arguments

x | x-coordinates of the regression plot
y | y-coordinates
main | main title
scr | split.screen must have been called beforehand; scr gives the screen number
value | function; the characteristic of interest is a function of the regression output.
variable | name of the characteristic of interest, plotted in the graphs
mode A scalar character of one of the following values 'silent', 'plot', or 'interactive':
'silent':
no output except the return of the result
'plot':
the regression line is plotted
'interactive':
the regression domain can be chosen interactively
Any mode includes all the functionalities of the modes that precede in the above enumeration. In the interactive mode, the regression domain is chosen by two mouse clicks with the left mouse; a right mouse click leaves the plot.
regr.select | vector of characters that select elements from the lm return list; these elements are returned by regression.
averaging | logical. If TRUE and y is matrix then the values of y are averaged row-wise beforehand
cex.main | font size of the title
col.passive | colour of the data points if plot is not active
col.active | colour of the data points if plot is active
col.main.passive | colour of the title if plot is not active
col.main.active | colour of the title if plot is active
col.points.chosen | colour of the selected points
col.abline | colour of the regression line for all points
col.abline.chosen | colour of the regression line for the selected points
col.smooth | colour for the smoothed data
... further graphical parameters

Value

The function returns a list with the following components

regr | part of the return list of lm
val | the calculated value of the variable of interest
regr.u | NULL or part of the return list of lm for x.u and y.u
val.u | NULL or the calculated value of the variable of interest for the restricted domain
x.u | NULL or the restricted x-coordinates given by the user in the interactive plot
y.u | NULL or y-coordinates according to x.u
sm | smoothed curve through the (x,y) points

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute


Simulation Techniques

Description

PrintMethodList prints the list of currently implemented methods for simulating random fields

GetMethodNames returns a list of currently implemented methods

Usage

PrintMethodList()

GetMethodNames()

Details

  • Boolean functions.
    See marked point processes.
  • circulant embedding.
    Introduced by Dietrich & Newsam (1993) and Wood and Chan (1994).

    Circulant embedding is a fast simulation method based on Fourier transformations. It is garantueed to be an exact method for covariance functions with finite support, e.g. the spherical model.

    See also cutoff embedding and intrinsic embedding for variants of the method.

  • cutoff embedding.
    Modified circulant embedding method so that exact simulation is garantueed for further covariance models, e.g. the whittle matern model. In fact, the circulant embedding is called with the cutoff hypermodel, see CovarianceFct, and A=B there. cutoff embedding halfens the maximum number of elements models used to define the covariance function of interest (from 10 to 5).

    Here multiplicative models are not allowed (yet).

  • direct matrix decomposition.
    This method is based on the well-known method for simulating any multivariate Gaussian distribution, using the square root of the covariance matrix. The method is pretty slow and limited to about 8000 points, i.e. a 20x20x20 grid in three dimensions. This implementation can use the Cholesky decomposition and the singular value decomposition. It allows for arbitrary points and arbitrary grids.
  • hyperplane method.
    The method is based on a tessellation of the space by hyperplanes. Each cell takes a spatially constant value of an i.i.d. random variables. The superposition of several such random fields yields approximatively a Gaussian random field.
  • intrinsic embedding.
    Modified circulant embedding so that exact simulation is garantueed for further variogram models, e.g. the fractal brownian one. Note that the simulated random field is always non-stationary. In fact, the circulant embedding is called with the Stein hypermodel, see CovarianceFct, and A=B there.

    Here multiplicative models are not allowed (yet).

  • Marked point processes.
    Some methods are based on marked point process P = ([x_1,m_1], [x_2,m_2], ...) where the marks m_i are deterministic or i.i.d. random functions on R^d.
    • add.MPP (Random coins).
      Here the functions are elements of the intersection (L1 cap L2) of the Hilbert spaces L1 and L2. A random field Z is obtained by adding the marks:

      Z(.) = sum_i m_i( . - x_i)

      In this package, only stationary Poisson point fields are allowed as underlying unmarked point processes. Thus, if the marks m_i are all indicator functions, we obtain a Poisson random field. If the intensity of the Poisson process is high we obtain an approximate Gaussian random field by the central limit theorem - this is the add.mpp method.

    • max.MPP (Boolean functions).
      If the random functions are multiplied by suitable, independent random values, and then the maximum is taken, a max-stable random field with unit Frechet margins is obtained - this is the max.mpp method.
  • nugget.
    The method allows for generating a random field of independent Gaussian random variables. In the isotropic case and if the simple notation of a model (with model and param) is used, this method is called automatically if the nugget effect is positive except the method "circulant embedding" or "direct" have been explicitely.

    The method has been extended to zonal anisotropies, see also parameter nugget.tol in RFparameters.

  • Random coins.
    See marked point processes.
  • spectral TBM (Spectral turning bands).
    The principle of spectral TBM does not differ from the other turning bands methods. However, line simulations are performed by a spectral technique (Mantoglou and Wilson, 1982); a realisation is given as the cosine with random amplitude and random phase. The implementation allows the simulation of 2-dimensional random fields defined on arbitrary points or arbitrary grids.
  • TBM2, TBM3 (Turning bands methods; turning layers).
    It is generally difficult to use the turning bands method (TBM2) directly in the 2-dimensional space. Instead, 2-dimensional random fields are frequently obtained by simulating a 3-dimensional random field (using TBM3) and taking a 2-dimensional cross-section. TBM3 allows for multiplicative models; in case of anisotropy the anisotropy matrices must be multiples of the first matrix or the anisotropy matrix consists of a time component only (i.e. all components are zero except the very last one).
    TBM2 and TBM3 allow for arbitrary points, and arbitrary grids (arbitrary number of points in each direction, arbitrary grid length for each direction).

    Note: Both the precision and the simulation time depend heavily on TBM*.linesimustep and TBM*.linesimufactor that can be set by RFparameters. For covariance models with larger values of the scale parameter, TBM*.linesimufactor=2 is too small.

    The turning layers are used for the simulations with time component. Here, if the model is a multiplicative covariance function then the product may contain matrices with pure time component. All the other matrices must be equal up to a factor and the temporal part of the anisotropy matrix (right column) may contain only zeros, except the very last entry.

Automatic selection algorithm

— details coming soon —

Note

Most methods possess additional parameters, see RFparameters() that control the precision of the result. The default parameters are chosen such that the simulations are fine for many models and their parameters. The example in EmpiricalVariogram() shows a way of checking the precision.

Author(s)

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

References

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

Lantuejoul, Ch. (2002) Geostatistical simulation. New York: Springer.

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.

Original work:

  • Circulant embedding:

    Chan, G. and Wood, A.T.A. (1997) An algorithm for simulating stationary Gaussian random fields. J. R. Stat. Soc., Ser. C 46, 171-181.

    Dietrich, C.R. and Newsam, G.N. (1993) A fast and exact method for multidimensional Gaussian stochastic simulations. Water Resour. Res. 29, 2861-2869.

    Wood, A.T.A. and Chan, G. (1994) Simulation of stationary Gaussian processes in [0,1]^d J. Comput. Graph. Stat. 3, 409-432.

  • Intrinsic embedding and Cutoff embedding:

    Stein, M.L. (2002) Fast and exact simulation of fractional Brownian surfaces. J. Comput. Graph. Statist. 11, 587–599.

    Gneiting, T., Sevcikova, H., Percival, D.B., Schlather, M. and Jiang, Y. (2005) Fast and Exact Simulation of Large Gaussian Lattice Systems in R^2: Exploring the Limits J. Comput. Graph. Statist. Submitted.

  • Turning bands method (TBM), turning layers:

    Dietrich, C.R. (1995) A simple and efficient space domain implementation of the turning bands method. Water Resour. Res. 31, 147-156.

    Mantoglou, A. and Wilson, J.L. (1982) The turning bands method for simulation of random fields using line generation by a spectral method. Water. Resour. Res. 18, 1379-1394.

    Matheron, G. (1973) The intrinsic random functions and their applications. Adv. Appl. Probab. 5, 439-468.

    Schlather, M. (2004) Turning layers: A space-time extension of turning bands. Submitted

  • Random coins:

    Matheron, G. (1967) Elements pour une Theorie des Milieux Poreux. Paris: Masson.

See Also

GaussRF, MaxStableRF, PrintModelList, RandomFields.

Examples

 PrintMethodList()

############################################################
##                                                        ##
##                     Figure 1 in                        ##
##           Gneiting, T., Sevcikova, H.,                 ##
##   Percival, D.B., Schlather, M. and Jiang, Y. (2005)   ##
##                                                        ##
############################################################
## the example below shows for which parameter combinations
## of the stable model, the cutoff method works
## note the method 'cutoff' itself has an automatic
## selection algorithm implemented and is called with
## the stable model only.
## Not run:

stabletest <- function(alpha, theta, size=512) {
  RFparameters(CE.trials=1, CE.tolIm = 1e-8, CE.tolRe=0, CE.force = FALSE,
               CE.useprimes=TRUE, CE.strategy=0, Print=2)
  model <- list(list(model="cutoff", var=1, kappa=c(1, theta, 1)),
                "(",
                list(model="stable", var=1, scale=1, kappa=alpha)
                )
  check <- CheckAndComplete(model=model, dim=2)
  if (check$error != 0) {
    cat("error", check$error, "\n")
    return(NA)
  }
  r <- paramextract(check$param[ , 1])$cutoffr
#  x <- c(r/size, r,  r /  size)
  x <- c(0, r, r / (size - 1)) * theta
  RFparameters(CE.mmin=rep(2 * size, 2))
  return(if (InitGaussRF(x, x, grid=TRUE, gridtriple=TRUE,
                         model=model, meth="circulant") == 0)
         r else NA)
}

alphas <- seq(1.52, 2.0, 0.01) # 0.02
thetas <- seq(0.05, 3.5, 0.05) #0.05( 2x)
m <- matrix(NA, nrow=length(thetas), ncol=length(alphas))
for (it in 1:length(thetas)) {
  theta <- thetas[it]
  for (ia in 1:length(alphas)) {
    alpha <- alphas[ia]
    cat("alpha=", alpha, "theta=", theta,"\n")
    print(unix.time(m[it, ia] <- stabletest(alpha=alpha, theta=theta)))
    if (is.na(m[it, ia])) break
  }
  if (any(is.finite(m))) image(thetas, alphas, m, col=rainbow(100))
}

## End(Not run)


Control Parameters

Description

RFparameters sets and returns control parameters for the simulation of random fields

Usage

   RFparameters(..., no.readonly=FALSE)

Arguments

... | arguments in tag = value form, or a list of tagged values.
no.readonly | If RFparameters is called without parameter then all parameters are returned in a list. If no.readonly=TRUE then only rewritable parameters are returned.

Details

The possible parameters are

General options

PracticalRange
logical or 2, 3, 11, 12, 13. If not FALSE the range of the covariance functions is adjusted so that cov(1) is about 0.05 (for scale=1).
  • FALSE : the practical range ajustment is not used.
  • TRUE : PracticalRange is applicable only if the value is known exactly.
  • 2 : PracticalRange is applicable if the value is known pretty well
  • 3 : PracticalRange is applicable if the value is roughly known
  • 11 : if the practical range is not known exactly it is approximated numerically.
  • 12 : if the practical range is not known pretty well it is approximated numerically.
  • 13 : if the practical range is even not known approximately it is approximated numerically.

Note that values beyond FALSE, TRUE, and 11, are only used for specialists' purposes.

Default: FALSE [init].

PrintLevel
If PrintLevel<=0 there is not any output on the screen. The higher the number the more tracing information is given. Default: 1 [init, do].
1 : error messages
2 : messages about partial failures of the algorithm
>2 : additional informations

Note that PrintLevel is also used in other packages as a default, for example in SoPhy (risk.index and create.roots). The changing of PrintLevel here may cause some unexpected effects in these functions. See the documentation there.

General options for simulating
pch
Character or empty string. The character is printed after each performed simulation if more than one simulation is performed at once. If pch='!' then an absolute counter is shown instead of the character. If pch='%' then a counter of percentages is shown instead of the character. Note that also '^H's are printed in the last two cases, which may have undesirable interactions with some few other R functions, e.g. Sweave. Default: '*' [do].
Storing
Logical. If FALSE then the intermediate results are destroyed after the simulation of the random field(s) or if an error had occured.

On the other hand, if Storing=TRUE, then several simulations for the same model parameters are performed faster, see the examples.

Note that in subsequent calls of GaussRF intermediate changes of RFparameters with flag "[init]" do not have any influence on the algorithm, if Storing=TRUE.

See alse CE.several, TBMCE.several and local.several for related parameters.

Default: FALSE [do].

skipchecks
logical. If TRUE, the check whether the given parameter values and the dimension are within the allowed range is skipped. Do not change the value of this variable except you really know what you do.

Default: FALSE [init].

stationary.only
Logical or NA. Used for the automatic choice of methods. See also RFMethods.
  • TRUE: the simulation of non-stationary random fields is refused. In particular, the intrinsic embedding method is excluded and the simulation of Brownian motion is rejected.
  • FALSE: intrinsic embedding is always allowed, actually it's the first one considered in the automatic selection algorithm.
  • NA: the simulation of the Brownian motion allowed, but intrinsic embedding is not used for stationary random fields.

Default: NA [init].

exactness
logical or NA.
  • TRUE: add.MPP, hyperplanes and all turning bands methods are excluded. If the circulant embedding method is considered as badly behaved, then the matrix decomposition methods are preferred.
  • FALSE: if the circulant embedding method is considered as badly behaved or the number of points to be simulated is large, the turning bands methods are rather preferred.
  • NA: approximative non-exact methods are excluded, i.e. TBM2 if the Abel transform of the covariance function cannot be given explicitely.

Default: NA [init].

Options for simulating with the standard circulant embedding method

CE.force
Logical. Circulant embedding does not work if a certain circulant matrix has negative eigenvalues. Sometimes it is convenient to replace all the negative eigenvalues by zero (CE.force=TRUE) after CE.trials number of trials. Default: FALSE [init].
CE.mmin
Scalar or vector, integer if positive. CE.mmin determines the initial size of the circulant matrix. If CE.mmin=0 the minimal starting size is determined automatically according to the dimensions of the grid. If CE.mmin>0 then the absolute starting size is given. If CE.mmin<0 then the automatically determined matrix size is multiplied by |CE.mmin|; here CE.mmin must be smaller than -1; the value -1 takes over the minimal starting size.
Note: in any cases, the initial size might be increased according to CE.useprimes. Default: 0 [init].
CE.strategy
0 : if the circulant matrix has negative eigenvalues then the size in each direction is doubled;
1 : the size is enhanced only in one direction, namely that one where the covariance function has the largest value at the end point of the grid — note that the default value of CE.trials is probably too small in that case.

In some cases CE.strategy=0 works better, in other cases CE.strategy=1. Just try.

Default: 0 [init].

CE.maxmem
maximal total size of the circulant matrix. The total amount of memory needed for the internal calculations is about 16 (=2 * sizeof(double)) times as large as CE.maxmem if RFparameters()$Storing=FALSE and 32 (=4 * sizeof(double)) time as large if Storing=TRUE.

Note that CE.maxmem can be used to control the automatic choice of the simulation algorithm. Namely, in case of huge circulant matrices, other simulation methods (TBM) are faster and might be preferred be the user.

Default: 4096^2 = 16777216 [init].

CE.tolIm
If the modulus of the imaginary part is less than CE.tolIm then the eigenvalue is considered as real. Default: 1E-3 [init].
CE.tolRe
Eigenvalues between CE.tolRe and 0 are considered as 0 and set 0. Default: -1E-7 [init].
CE.trials
A larger circulant matrix is likely to make more eigenvalues non-negative. If at least one of the thresholds CE.tolRe and CE.tolIm are missed then the matrix size is doubled according to CE.strategy, and the matrix is checked again. This procedure is repeated up to CE.trials-1 times. If there are still negative eigenvalues, the simulation method fails if CE.force=FALSE. Default: 3 [init].
CE.several
logical. If FALSE only half the memory is need, but only a single independent realisation can created. Default: TRUE [init].
CE.useprimes
Logical. If FALSE the columns of the circulant matrix have length 2^k for some k. Otherwise the algorithm tries to find a nicely factorizable number close to the size of the given matrix. Default: TRUE [init].
CE.dependent
Logical. If FALSE then independent random fields are created. If TRUE then at least 4 non-overlapping rectangles are taken out of the the expanded grid defined by the circulant matrix. These simulations are dependent. See below for an example. See CE.trials for some more information on the circulant matrix. Default: FALSE [init].

Options for simulating with the local ce methods (cutoff, intrinsic)

local.force
see CE.force above. Default: FALSE [init].
local.mmin
see CE.mmin above.

Difference: if local.mmin=0 the automatic determination of the initial size of the circulant matrix takes into account an expansion factor. This expansion factor is intended to make the circulant matrix positive definite and is either theoretically or numerically known, or guessed.

If the usual strategy of circulant embedding (doubling the grid sizes) should be taken over then local.mmin must be set to -1.

Default: 0 [init].

local.maxmem
see CE.maxmem above. Default: 20000000 [init].
local.tolIm
see CE.tolIm above. Default: 1E-7 [init].
local.tolRe
see CE.tolRe above. Default: -1E-9 [init].
local.several
see CE.several above. Default: 1 [init].
local.useprimes
see CE.useprimes above. Default: TRUE [init].
local.dependent
see CE.dependent above. Default: FALSE [init].

Options for simulating by simple matrix decomposition

direct.bestvariables
integer. When searching for an appropriate simuation method the matrix decomposition method (see ‘direct method’ below) is preferred if the number of variables is less than or equal to direct.bestvariables

Default: 800 [init]

direct.maxvariables
If the number of variables to generate is greater than direct.maxvariables, then any matrix decomposition method is rejected. It is important that this option is set conveniently to avoid great losses of time during the automatic search of a simulation method (method=NULL in GaussRF). Default: 4096 [init]
direct.method
Decomposition of the covariance matrix. If direct.method=1, Cholesky decomposition will not be attempted, but singular value decomposition used instead. Default: 0 [init].
direct.svdtolerance
If SVD decomposition is used for calculating the square root of the covariance matrix then the absolute componentwise difference between the covariance matrix and square of the square root must be less than direct.svdtolerance. No check is performed if direct.svdtolerance is negative. Default: 1e-12 [init].

Options for simulating nugget effects
Simulating a nugget effect seems trivial. It gets complicated and best methods (including direct and circulant embedding!) fail if zonal anisotropies are considered, where sets of points have to be identified that belong to the same subspace of eigenvalue 0 of the anisotropy matrix.

nugget.tol
points at a distance less than or equal to nugget.tol are considered as being identical. This strategy applies to the simulation method and the covariance function itself. Hence, the covariance function is only positive definite if nugget.tol=0.0. However, if the anisotropy matrix does not have full rank and nugget.tol=0.0 then, the simulations are likely to be odd. The value of nugget.tol should be of order 1e-15. Default: 0.0 [init].

Options for simulating with a turning bands method
Currently, there are 3 variants of the turning bands method implemented:

spectral
The spectral turning bands method is implemented for 2 (and 1) dimensions only.
TBM2
It is based on the two dimensional turning bands operator and is applicable for 1 and 2 dimensions. As an additional dimension the time dimension can be added.
TBM3
It is based on the three dimensional turning bands operator and is applicable for 1,2,3 dimensions. As an additional dimension the time dimension can be added.

The following parameters are used.

spectral.grid
Logical. The angle of the lines is random if spectral.grid=FALSE, and k*pi/spectral.lines for k in 1:spectral.lines, otherwise. Default: TRUE [do].
spectral.lines
Spectral turning bands. Number of lines used (in total for all additive components of the covariance function). Default: 500 [do].
TBM.method
character. The preferred method to simulate on the line for TBM2 and TBM3; currently either 'circulant embedding' or 'direct'. If 'direct' then the method is overwritten if the number of points on the grid is larger than direct.maxvariables. If the circulant embedding method is used, then the TBMCE parameters below determine the behaviour of the circulant embedding algorithm.

Default: "circulant embedding" [init].

TBM.center
Scalar or vector. If not NA, the TBM.center is used as the center of the turning bands for TBM2 and TBM3. Otherwise the center is determined automatically such that the line length is minimal. See also TBM.points and the examples below. Default: NA [init].
TBM.points
integer. If greater than 0, TBM.points gives the number of points simulated on the TBM line, hence must be greater than the minimal number of points given by the size of the simulated field and the two paramters TBMx.linesimufactor and TBMx.linesimustep. If TBM.points is not positive the number of points is determined automatically. The use of TBM.center and TBM.points is highlighted in an example below. Default: 0 [init].
TBM2.every
If TBM2.every>0 then every TBM2.everyth iteration is announced. Default: 0 [do].
TBM2.lines
Number of lines used. Default: 60 [do].
TBM2.linesimufactor
TBM2.linesimufactor or TBM2.linesimustep must be non-negative; if TBM2.linesimustep is positive then TBM2.linesimufactor is ignored. If both parameters are naught then TBM.points is used (and must be positive). The grid on the line is TBM2.linesimufactor-times finer than the smallest distance. See also TBM2.linesimustep. Default: 2.0 [init].
TBM2.linesimustep
If TBM2.linesimustep is positive the grid on the line has lag TBM2.linesimustep. See also TBM2.linesimufactor. Default: 0.0 [init].
TBM2.num
Logical. If TRUE then the covariance function on the line is approximated numerically. If FALSE only those models are allowed that have an analytic representation on the line. Default: TRUE [init].
TBM2.layers
Logical. If TRUE then the turning layers are used whenever a time component is given. If FALSE the turning layers are used only when the traditional TBM is not applicable. If negative then turning layers may never be used. If greater than 1 then only turning layers may be used. Default: FALSE [init].
TBM3.every
If TBM3.every>0 then every TBM3.everyth iteration is announced. Default: 0 [do].
TBM3.lines
Number of lines used. Default: 500 [do].
TBM3.linesimufactor
See TBM2.linesimufactor for the meaning. Default: 2.0 [init].
TBM3.layers
See TBM2.layers for the meaning. Default: FALSE [init].
TBM3.linesimus
See TBM2.linesimustep for the meaning. Default: 0.0 [init].
TBMCE.force
see TBM.method and CE.force Default: FALSE [init].
TBMCE.mmin
see TBM.method and CE.mmin. Default: 0 [init].
TBMCE.strategy
see TBM.method and CE.strategy. Default: 0 [init].
TBMCE.maxmem
see TBM.method and CE.maxmem. Default: 10000000 [init].
TBMCE.tolIm
see TBM.method and CE.tolIm. Default: 1E-3 [init].
TBMCE.tolRe
see TBM.method and CE.tolRe. Default: -1E-7 [init].
TBMCE.trials
see TBM.method and CE.trials. Default: 3 [init].
TBMCE.useprimes
see TBM.method and CE.useprimes. Default: TRUE [init].
TBMCE.dependent
see TBM.method and CE.dependent. Default: FALSE [init].

Options for simulating with Poisson point processes

add.MPP.realisations
Number of superposed realisations (to approximate the normal distribution; total number for all (additive) components with same anisotropy); Default: 100 [do].
MPP.approxzero
Functions that do not have compact support are set to zero outside the ball outside for which the function has absolute values less than MPP.approxzero. Default: 0.001 [init].
MPP.radius
In order avoid edge effects, the simulation area is enlarged by a constant r so that all marks have their (supposed) support in the ball with radius r centred at the origin; see also MPP.approxzero. If MPP.radius>0 the true radius r is replaced by MPP.radius. Default: 0.0 [init].

Options for simulating hyperplane tessellations

hyper.superpos
integer. number of superposed hyperplane tessellations. Default: 300 [do].
hyper.maxlines
integer. Maximum number of allowed lines. Default: 1000 [init].
hyper.mar.distr
integer. code for the marginal distribution used in the simulation:
0
uniform distribution
1
Frechet distribution with form parameter hyper.mar.param
2
Bernoulli distribution (Binomial with n=1) with parameter hyper.mar.param

The parameter should not be changed yet. Default: 0 [do].

hyper.mar.param
Parameter used for the marginal distribution. The parameter should not be changed yet. Default: 0 [do].

Options specific to simulating max-stable random fields

maxstable.maxGauss
Max-stable random fields. The simulation of the max-stable process based on random fields uses a stopping rule that necessarily needs a finite upper endpoint of the marginal distribution of the random field. In the case of extremal Gaussian random fields, see MaxStableRF, the upper endpoint is approximated by maxstable.maxGauss. Default: 3.0 [init].

General comments on the options
The following refers to the simulation of Gaussian random fields (InitGaussRF, GaussRF), but most parts also apply for the simulation of max-stable random fields (InitMaxStableRF, MaxStableRF).

Some of the global parameters determine the basic settings of a simulation, e.g. direct.method (which chooses a square root of a positive definite matrix). The values of such parameters are read by InitGaussRF and stored in an internal register. Changing such a parameter between calling InitGaussRF and calling DoSimulateRF or between subsequent calls of GaussRF will not have any effect. These parameters have the flag "[init]".

Parameters like TBM2.lines (which determines the number of i.i.d. processes to be simulated on the line) are only relevant when generating random numbers. These parameters are read by DoSimulateRF (or by the second part of GaussRF), and are marked by "[do]".

Storing has an influence on both, InitGaussRF and DoSimulateRF. InitGaussRF may reserve more memory if Storing=TRUE. DoSimulateRF will free the register if Storing=FALSE, whatever the value of Storing was when InitGaussRF was called.

The distinction between [init] and [do] is also relevant if GaussRF is used and called a second time with the same parameters for the random field and if RFparameters()$Storing=TRUE. Then GaussRF realises that the second call has the same random field parameters, and takes over the stored intermediate results (that have been calculated with the RFparameters() at that time). To prevent the use of stored intermediate results or to take into account intermediate changes of RFparameters set RFparameters(Storing=FALSE) or use DeleteRegister() between calls of GaussRF.

A programme that checks whether the parameters are well adapted to a specific simulation problem is given as an example of EmpiricalVariogram().

For further details on the implemented methods, see RFMethods.

Value

If any parameter has been given RFparameters returns an invisible list of the given parameters in full name.
Otherwise the complete list of parameters is returned. Further the values of the following internal readonly variables are returned:
* covmaxchar: max. name length for variogram/covariance models
* covnr: number of currently implemented variogram/covariance models
* distrmaxchar: max. name length for a distribution
* distrnr: number of currently implemented distributions
* maxdim: maximum number of dimensions for a random field
* maxmodels: maximum number of elementary models in definition of a complex covariance model
* methodmaxchar: max. name length for methods
* methodnr: number of currently implemented simulation methods

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute

References

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

GaussRF, GetPracticalRange, MaxStableRF, RandomFields, and RFMethods.

Examples

 RFparameters(Storing=TRUE)
 str(RFparameters())

############################################################
##                                                        ##
##       use of TBM.points and TBM.center                 ##
##                                                        ##
############################################################
## The following example shows that the same realisation
## can be obtained on different grid geometries (or point
## configurations) using TBM3 (or TBM2)

x1 <- seq(-150,150,1)
y1 <- seq(-15, 15, 1)
x2 <- seq(-50, 50, 1)
model <- "exponential"
param <- c(0, 1, 0, 10)

meth <- "TBM3"
###### simulation of a random field on long thing stripe
runif(1)
rs <- get(".Random.seed", envir=.GlobalEnv, inherits = FALSE)
DeleteAllRegisters()
z1 <- GaussRF(x1, y1, model=model, param=param, grid=TRUE, register=1,
              method=meth, TBM.center=0, Storing=TRUE)
do.call(getOption("device"), list(height=1.55, width=12))
par(mar=c(2.2, 2.2, 0.1, 0.1))
image(x1, y1, z1, col=rainbow(100))
polygon(range(x2)[c(1,2,2,1)], range(y1)[c(1,1,2,2)], border="red", lwd=3)

###### definition of a random field on a square of shorter diagonal
assign(".Random.seed", rs, envir=.GlobalEnv)
z2 <- GaussRF(x2, x2, model=model, param=param, grid=TRUE, register=2,
              method=meth, TBM.center=0,
              TBM.points=length(GetRegisterInfo(1)$method[[1]]$mem$l))
do.call(getOption("device"), list(height=4.3, width=4.3))
par(mar=c(2.2, 2.2, 0.1, 0.1))
image(x2, x2, z2, zlim=range(z1), col=rainbow(100))
polygon(range(x2)[c(1,2,2,1)], range(y1)[c(1,1,2,2)], border="red", lwd=3)

############################################################
##                                                        ##
##                   use of exactness                     ##
##                                                        ##
############################################################
x <- seq(0, 1, 1/30)
model <- list(list(model="stable", var=1, scale=1, kappa=1.0),
              "+",
              list(model="gencauchy", var=1, scale=1, kappa=c(1, 2))
              )

for (exactness in c(NA, FALSE, TRUE)) {
  readline(paste("\n\nexactness: ", exactness, "; press return"))
  DeleteRegister()
  z <- GaussRF(x, x,  grid=TRUE, gridtriple=FALSE,
                model=model, exactness=exactness,
                stationary.only=NA, Print=4, n=1,
                TBM2.linesimustep=1, Storing=TRUE)
   str(lapply(GetRegisterInfo()$method,
              function(x) x[c("name", "covnr")]))
 }

 #############################################################
 ## The following gives a tiny example on the advantage of  ##
 ## local.dependent=TRUE (and CE.dependent=TRUE) if in a    ##
 ## study most of the time is spent with simulating the     ##
 ## Gaussian random fields. Here, the covariance at a pair  ##
 ## of points is estimated.                                 ##
 #############################################################

# In the example below, local.dependent speeds up the simulation
# by about factor 27 at the price of an increased variance of
# factor 1.5

x <- seq(0, 1, len=10)
y <- seq(0, 1, len=10)
grid.size <- c(length(x), length(y))
model <- list(list(model="exp", var=1.1, aniso=c(2,1,0.5,1)))
CovarianceFct(matrix(c(1, -1), ncol=2), model=model) ## true value

RFparameters(Storing=TRUE)
m <- if (interactive()) 1000 else 2

# determine number of non-overlapping realisations on the torus
DeleteRegister()
InitGaussRF(x, y, model=model, grid=TRUE, method="cu")
blocks <- GetRegisterInfo()$method[[1]]$mem$new$method[[1]]$mem$simupos
(n <- prod(blocks) * 1) ## n any multiple of prod(blocks) to avoid
##               dependencies between the m estimated covariance if
##               if local.dep=TRUE; or put RFparameters(Storing=FALSE),
##               but this is slower

# using local.dependent=TRUE...
c1 <- numeric(m)
DeleteRegister()
unix.time(
for (i in 1:m) {
  cat("\n", i)
  z <- GaussRF(x, y, model=model, grid=TRUE, method="cu", n=n,
               local.dependent=TRUE, pch="")
  c1[i] <- cov(z[1,length(y),], z[length(x), 1 , ])
}) # about 35 0.3 35 0  0
var(c1) # about 0.013

# using local.dependent=FALSE...
c2 <- numeric(m)
DeleteRegister()
unix.time(
for (i in 1:m) {
  cat("\n", i)
  z <- GaussRF(x, y, model=model, grid=TRUE, method="cu", n=n,
               local.dependent=FALSE, pch="")
  c2[i] <- cov(z[1,length(y),], z[length(x), 1 , ])
})  # about 950  3 950 0 0
var(c2) # about 0.0087


Interactive Choice of Models and Parameters

Description

ShowModels is an interactive plot for the selection of models and their one- or two-dimensional simulations; it also allows for the fitting of variogram models by eye.

Usage

ShowModels(x, y=NULL,
           covx=ifelse(is.null(empirical), diff(range(x))/5,
                       max(empirical$c)),
           fixed.rs=TRUE, method=NULL, empirical=NULL,
           model=NULL, param=NULL, all.param=NULL,
           legends = TRUE,
           register=0, Mean=NULL, erase=TRUE,
           x.fraction=0.60, cex.names=1, covx.default = 100,
           link.fct=NULL, Zlim=NULL,
           Col.rect="red", Col.bg="blue", Col.sep="grey",
           Col.line="red", Col.txt="black", Col.flash="red",
           Col.vario="blue", Col.main="black",
           Col.model=c("red", "black"),  vario.lty=c(1,2),
           cex.leg = 0.7 * cex.names,  cex.eval = 0.8 * cex.names,
           update=TRUE, screen.new=TRUE, use.outer.RFparameters=FALSE,
           debug=FALSE, ...)

Arguments

x | if NULL simulations are not performed; otherwise it gives the x coordinates of a grid as a sequence of increasing numbers
y | if NULL at most one-dimensional simulations are performed (depending on the value of x); otherwise y gives the y coordinates of a two-dimensional grid (as a sequence of increasing numbers).
covx | if a single value is given, it is the largest distance for which the covariance functions or the variograms are plotted; otherwise the models are plotted for the given values, and the origin.
fixed.rs | if TRUE then the same random seed is used for all simulations until the user clicks on the formula, the title or the subtitles.
method | simulation method, see RFMethods; if NULL then a suitable simulation method is chosen automatically.
empirical | empirical variogram; a list as returned by EmpiricalVariogram. Also empirical variograms with a pair number of anisotropy directions may be passed. Then the first and the middle one are taken.
model covariance model, see CovarianceFct, or type PrintModelList() to get all options. If given, this model is shown at the beginning. Additive or multiplicative models are not allowed.
However, model can also be given by a simple list definition, see CovarianceFct. Then param must not be given. In this case also anisotropic models can be defined.
param | parameter vector: param=c(mean, variance, nugget, scale,...); the parameters must be given in this order; see CovarianceFct for more details. Only considered if model is given. If given, model is initialised by param.
all.param all.param=c(mean, variance, nugget, scale); the parameters must be given in this order; If all.param is given then the parameters of all covariance functions are set to the given values. The values are overwritten for a specific model if model and param are given. Note that it is not possible to set the values of additional (form) parameters of a parametrised class by means of all.param.
In case of an anisotropic model the anisotropy matrix is by default diagonal with both entries equal to 1/all.param[4].
legends | if TRUE then a legend is added to the two-dimensional plot.
register | register where intermediate results of the simulations are stored, see also GaussRF.
Mean | mean of the random field
erase | parameter of split.screen, which is called at the very beginning
x.fraction | the current screen is split into 2 x 2 screens. The parameter x.fraction gives the size of the left screens in the x directions as part of 1. See also the Details.
cex.names | font size for model names
covx.default | if length(cov.x)==1 then [0, cov.x] is covered by covx.default points of equal distance
link.fct NULL or function(values) or "MaxStable". Transformation of the Gaussian random field.

If link.fct="MaxStable" then max-stable random fields are simulated for the given covariance function and the extremal coefficient function is given (up the constant -1) instead of the variogram or the covariance function
Zlim | Vector of two elements or list of two vectors of two elements. Graphical limits for the Gaussian random process (and the transformed field).
Col.rect | colour for interactive plot; see eval.parameters.
Col.bg | colour for interactive plot; see eval.parameters.
Col.sep | colour for interactive plot; see eval.parameters.
Col.line | colour for interactive plot; see eval.parameters.
Col.txt | colour for interactive plot; see eval.parameters.
Col.flash | colour for the previously chosen model
Col.vario | colour for the empirical variogram plot
Col.main | colour for the title of the random field
Col.model | vector of two colours for plotting the variogram of the Gaussian random field and the transformed field
vario.lty | vector of two line types for primary and secondary axis of the variogram
cex.leg | font size used in the legends
cex.eval | font size used in the menue entries
update | logical. If TRUE the plots are updated after each interactive change of the values. Otherwise, the bottom 'simulate' is added in the menu.
screen.new | logical. If FALSE the screen is erased before a simulation and completely rebuild; otherwise the screen is updated. If FALSE flickering appears during the update of the current screen, otherwise it may happen during the reorganisation of any window (and may take quite a lot of time).
use.outer.RFparameters logical. If FALSE the following parameters usually set by RFparameters are internally set
  • PracticalRange=FALSE
  • PrintLevel=1 if debug=FALSE and 5 otherwise.
  • maxstable.maxGauss=2
  • CE.force=TRUE
  • CE.trials=1
  • CE.mmin=-4
  • CE.useprimes=TRUE
debug | logical. If TRUE then internally the RFparameter()$PrintLevel is set to 5.
... additional graphics options for the plot of the one- or two-dimensional simulations, see plot and image.

Details

The interactive plot consists of 3 parts:

  • top left: graph of the covariance function or the variogram. In case empirical is given the empirical variogram is also plotted. If link.fct is given, then also the covariance function or the variogram is plotted. If the correlation model is for a non-stationary random field, the variogram for the transformed random field is not estimated in a primitive way – this is indicated with a star in the legend
  • bottom left: one- or two-dimensional simulation
  • right:

    – list of implemented models; a specific model is chosen by the left mouse button, or:

    – menu for the parameters of the chosen model. The list includes the variance, a nugget effect, the mean and the scale or the anisotropy parameters. Further, some global parameters can be changed. They are the PracticalRange (see RFparameters for details) and the angle of the main variogram direction (or NA, then it follows the angle of the anisotropy). Finally, the user can choose between the plot of the covariance and the corresponding variogram.

The interactive plot is left by clicking any mouse button different from the left when the top right part is active.

Value

list of the last model and its parameters.

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute

See Also

CovarianceFct, eval.parameters, GaussRF, RFMethods, RandomFields.

Examples

# first example: one-dimensional simulations

options(locatorBell=FALSE)
x <-  seq(1,10,0.1);
ShowModels(x=x)

# second example: two-dimensional simulations and
#                 empirical variogram
dx <- runif(300,0,8)
dy <- runif(300,0,8)
dz <- GaussRF(x=dx, y=dy, grid=FALSE, model="gaus",
      param=c(1,2,1,2))
ev <- EmpiricalVariogram(x=dx, y=dy, data=dz, grid=FALSE,
      bin=(-1:20)/4)
x <-  seq(1,5,0.1);
ShowModels(x=x, y=x, empirical=ev)

# third example: two-dimensional anistropic simulations and
#                link function
x <- seq(1,10,0.1)
ShowModels(x=x, y=x, link=function(x) exp(x),
          model=list(list(model="spheric", var=1, aniso=c(1,0,0,5))))

x <-  seq(1,10,0.1)
ShowModels(x=x, link=function(x) exp(x),
           model=list(list(model="spheric",var=1, scale=1)))

x <-  seq(1,10,0.1)
ShowModels(x=x, link="MaxStable", fixed.rs=TRUE,
           model=list(list(model="gauss",var=1, scale=1)), type="l")


Simulation of Random Fields

Description

DoSimulateRF performs an already initialised simulation.

InitSimulateRF internal function; use InitGaussRF and InitMaxStableRF, instead.

Usage

DoSimulateRF(n=1, register=0, paired=FALSE)

InitSimulateRF(x, y=NULL, z=NULL, T=NULL, grid, model, param, trend,
               method=NULL, register=0, gridtriple=FALSE, distribution=NA)

Arguments

x | matrix of coordinates, or vector of x coordinates
y | vector of y coordinates
z | vector of z coordinates
T | time instances
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 | vector or list. param=c(mean, variance, nugget, scale, ...), param=list(c(variance, scale, ...), ..., c(variance,scale,...)), param=matrix(...), or param=list(list(variance, anisotropy, kappa),..., list(variance, anisotropy, kappa)); the parameters must be given in this order; further parameters are to be added in case of a parametrised class of models, see CovarianceFct
trend | Not programmed yet. trend surface: number (mean), vector of length d+1 (linear trend a_0 +a_1 x_1 + ... + a_d x_d), or function
method | NULL or string; Method used for simulating, see RFMethods, or type PrintMethodList() to get all options
register | 0:9; place where intermediate calculations are stored; the numbers are aliases for 10 internal registers
gridtriple | logical; if gridtriple=FALSE ascending sequences for the parameters x, y, and z are expected; if gridtriple=TRUE triples of form c(start,end,step) expected; this parameter is used only if grid=TRUE
distribution | marginal distribution: 'Gauss', 'Poisson', or 'MaxStable'
n | number of realisations to generate; if paired=TRUE then n must be even.
paired | logical. paired may be TRUE only for the simulation of Gaussian random fields. If TRUE then every second simulation is obtained by only changing the signs of the standard Gaussian random variables, the simulation is based on (“antithetic pairs”).

Value

InitSimulateRF returns 0 if no error has occurred during the initialisation process, and a positive value if failed.

DoSimulateRF returns NULL if an error has occurred; 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 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.

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute

See Also

GaussRF, MaxStableRF, RandomFields


Sleep

Description

Process sleeps for a given amount of time

Usage

sleep.milli(milli)

Arguments

milli sleeping time in milliseconds

Value

No value is returned.

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute


Soil data of North Bavaria, Germany

soil(RandomFields) R Documentation

Soil data of North Bavaria Germany

Description

Soil physical and chemical data collected on a field in the Weissenstaedter Becken, Germany

Usage

data(soil)

Format

This data frame contains the following columns:

x.coord
x coordinates given in cm
y.coord
y coordinates given in cm
nr
number of the samples, which were taken in this order
moisture
moisture content [Kg/Kg * 100%]
NO3.N
nitrate nitrogen [mg/Kg]
Total.N
total nitrogen [mg/Kg]
NH4.N
ammonium nitrogen [mg/Kg]
DOC
dissolved organic carbon [mg/Kg]
N20N
nitrous oxide [mg/Kg dried substance]

Details

For technical reasons some of the data were obtained as differences of two measurements (which are not available anymore). Therefore, some of the data have negative values.

Source

The data were collected by Wolfgang Falk, Soil Physics Group, http://www.geo.uni-bayreuth.de/bodenphysik/Welcome.html, University of Bayreuth, Germany.

References

Falk, W. (2000) Kleinskalige raeumliche Variabilitaet von Lachgas und bodenchemischen Parameters [Small Scale Spatial Variability of Nitrous Oxide and Pedo-Chemical Parameters], Master thesis, University of Bayreuth, Germany.

Examples

################################################################
##                                                            ##
##         a geostatistical analysis that demonstrates        ##
##         features of the package `RandomFields'             ##
##                                                            ##
################################################################

data(soil)
names(soil)
pts <- soil[,c(1,2)]
d <- soil$moisture

## define some graphical parameters first
close.screen(close.screen())
maxbin <- max(dist(pts)) / 2
(zlim <- range(d))
cn <- 100
colour <- rainbow(cn)
par(bg="white",cex=1, cex.lab=1.4, cex.axis=1.4, mar=c(4.3,4.3,0.8,0.8))
lu.x <- min(soil$x)
lu.y <- max(soil$y)
y <- x <- seq(min(soil$x), max(soil$x), l=121)

## ... and a certain appearance of the legend
my.legend <- function(lu.x, lu.y, zlim, col, cex=1) {
  ## uses already the legend code of R-1.3.0
  cn <- length(col)
  filler <- vector("character", length=(cn-3)/2)
  legend(lu.x, lu.y, y.i=0.03, x.i=0.1,
         legend=c(format(zlim[2], dig=2), filler,
         format(mean(zlim), dig=2), filler,
         format(zlim[1], dig=2)),
         lty=1, col=rev(col),cex=cex)
}

## plot the data first
plot(pts, col=colour[1+(cn-1)*((d-min(d))/diff(zlim))], pch=16,
     xlab="x [cm]", ylab="y [cm]", cex.axis=1.5, cex.lab=1.5)
my.legend(lu.x, lu.y, zlim=zlim, col=colour, cex=1.5)

## empirical variogram
ev <- EmpiricalVariogram(pts, data=d, grid=FALSE,
                         bin=c(-1,seq(0,maxbin,l=30)))

## show all models,
by.eye <- ShowModels(x=x, y=y, emp=ev, col=colour, Zlim=zlim,
                     Mean=mean(d), me="ci")

## fit parameters of the whittlematern model by MLE
fit <- fitvario(x=pts, data=d, model="whittle", par=rep(NA,5),
                mle.m="ml", cross.m=NULL)$variogram

## plot the fitted model and the empirical variogram
plot(ev$c, ev$emp.var, ylim=c(0,11), ylab="variogram", xlab="lag")
gx <- seq(0.001, max(ev$c), l=100)
if(!is.null(by.eye)) lines(gx, Variogram(gx, model=by.eye))
lines(gx, Variogram(gx, model=fit$ml), col=2)
lines(gx, Variogram(gx, model=fit$plain), col=3)
lines(gx, Variogram(gx, model=fit$sqrt.nr), col=4)
lines(gx, Variogram(gx, model=fit$sd.inv), col=5)
legend(120, 4, c("empirical", "by eye", "ML", "lsq", "sqrt(n) lsq",
               "sd^-1 lsq"),
       lty=c(-1, rep(1, 5)), pch=c(1, rep(-1, 5)),
       col=c(1, 1, 2, 3, 4, 5), cex=1.4)

## map of expected values
k <- Kriging("O", x=x, y=y, grid=TRUE, model=fit$ml, given=pts, data=d)
par(mfrow=c(1,2))
plot(pts, col=colour[1+99*((d-min(d))/diff(zlim))], pch=16, xlab="x [cm]", ylab="y [cm]")
my.legend(lu.x, lu.y, zlim=zlim, col=colour, cex=1)
image(x, y, k, col=colour, zlim=zlim, xlab="x [cm]", ylab="y [cm]")
par(bg="white")

## what is the probability that at no point of the
## grid given by x and y the moisture is greater than 24 percent?
RFparameters(Print=1, CE.force=FALSE, CE.trials=3, CE.useprimes=TRUE)
cs <- CondSimu("O", x=x, y=y, grid=TRUE, model=fit$ml, given=pts,
               data=d, n=10) # better n=100 or n=1000
par(mfrow=c(2,3))
image(x, y, k, col=colour, zlim=zlim, xlab="x [cm]", ylab="y [cm]")
my.legend(lu.x, lu.y, zlim=zlim, col=colour, cex=0.5)
for (i in 1:5)
  image(x, y, cs[, , i], col=colour, zlim=zlim,
        xlab="x [cm]", ylab="y [cm]")
mean(apply(cs<=24, 3, all)) ## about 40 percent ...


Set input behaviour

Description

The functions store values that influence the behaviour of Locator and Readline

Usage

useraction(action=c("none", "start.register", "continue.register",
                    "replay", "endless", "delete"),
                    sleep, wait, PrintLevel, actionlist)
putactions(l)

Arguments

action string of the following values
‘"none"’
the behaviour of Readline and Locator is that of readline and locator, respectively. The values are not stored.
‘"start.register"’
the behaviour of Readline and Locator is that of readline and locator, respectively. The current list of stored values is deleted and new inputs by Readline and Locator are stored.
‘"continue.register"’
the behaviour of Readline and Locator is that of readline and locator, respectively. The current list of stored values is continued with inputs from Readline and Locator.
‘"replay"’
Locator and Readline do not get the data from the mouse or terminal, but read them from the stored list. If the end of storage list is reached, the user is informed and the action is switched to “none”.
‘"endless"’
similar to ‘"replay"’, but storage list is looped.
‘"delete"’
The current list of stored values is deleted and action is set to “none”
actionlist | list of stored values by Readline and Locator
sleep | sleeping time before the blinking starts (replay mode of Locator); value is kept until explicit changing; starting value: 2.5
wait | regulated the blinking speed of the characters (Locator) and of the speed of appearance of characters (Readline); value is kept until explicit changing; starting value: 0.1
PrintLevel | if Printlevel>1 and Readline or Locator is in replay mode then the current pointer position within the storage list and the tracing information (if any) is shown. The value is kept until explicit changing; starting value: 0
l | list of all parameters relevant to the behaviour of Readline and Locator.

Details

The function putactions restores the state taken by getactions.

Concerning the function useraction, the parameter action must be given. If any of the other parameters is not given, its current value is kept. If action="start.register" or action="delete", the value of actionlist is ignored.

If replaying, Readline shows the input values on the screen and Locator shows blinking characters as side effects.

Value

useraction returns NULL.
putactions returns invisibly the state before the new state is set.

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute

See Also

eval.parameters, getactions, Locator, Readline,

Examples

sample.interaction <- function(read.r2=TRUE, type="p", pch=16) {
  r1 <- Readline(prompt="first input: ")
  r2 <- if (read.r2) Readline(prompt="second input: ", info="2nd input")
  else " -- no input  -- "
  cat("Mark some points with the right mouse key, then leave with the left mouse key\n")
  do.call(getOption("device"), list())
  plot(Inf, Inf, xlim=c(-1,1), ylim=c(-1,1), xlab="", ylab="")
  l <- Locator(100, info="locator input", type=type, pch=pch)
  r3 <- Readline(prompt="third input: ", info="last input")
  return(list(r1=r1, r2=r2, l=l, r3=r3))
}

#########################
## user terminal input ##
#########################
{
  useraction("start.register")
  str(sample.interaction())
  str(l <- getactionlist())
}

#################
## just replay ##
#################
useraction("replay", sleep=2, wait=interactive() * 0.2 * 5,
           PrintLevel=2, actionlist=l)
str(sample.interaction(type="l"))

##############################################
## modify first the input list, then replay ##
##############################################
l2 <- l[-2]
l2[[1]] <- list("some other words", info="changed input")
str(l2)
useraction("replay", sleep=1, wait=interactive() * 0.05,
           PrintLevel=0, l2)
str(sample.interaction(read.r2=FALSE, type="o")) # input now from l2
Readline(prompt="new input: ", info="?!") # switch to terminal
##                                                since end of stored list
# str(getactionlist())    # new input has not been not stored ...

####################################################
## use of the two lists, l and l2, in a mixed way ##
####################################################
useraction("replay", sleep=2, wait=interactive() * 0.05,
           PrintLevel=0, actionlist=l)
Readline(prompt="first input of 1: ")
dump <- getactions()
useraction("replay", sleep=0.5, wait=interactive() * 0.05,
           PrintLevel=0, actionlist=l2)
Readline(prompt="first input of 2: ")
dump2 <- putactions(dump)
Readline(prompt="second input of 1: ")
## locator call reading from list 1:
plot(Inf, Inf, xlim=c(-1,1), ylim=c(-1,1), xlab="", ylab="")
Locator(100, info="locator input", type="p")
dump <- putactions(dump2)
##locator call reading from list 2:
Locator(100, info="locator input", type="p", pch=20, col="blue")
Readline(prompt="last input of 2: ")
putactions(dump)
Readline(prompt="last input of 1: ")

#####################################################
## see help("eval.parameters") for another example ##
#####################################################


General input function

Description

The function provides a general functionality to read and reply user input by locator and readline

Usage

userinput(fct, info=NULL, prompt, n,  type="n", pch=par()$pch,
          cex=par()$cex, ...)

Arguments

fct | either ‘locator’ or ‘readline’
info | arbitrary object; tracing information that is useful in case the user likes to edit the stored sequence of inputs
prompt | prompt parameter of readline
n | parameter n of locator
type | parameter type of locator
pch | character type used by locator
cex | character size used by locator
... | additional graphics parameters used by locator

Details

The behaviour of userinput depends on the the value of action set by useraction, see there for more information.

Value

If fct="readline" then a character vector of length one is returned.
If fct="locator" then the function returns a list containing 'x' and 'y' components which are the coordinates of the identified points in the user coordinate system, i.e., the one specified by 'par("usr")'.

Author(s)

Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute

See Also

getactions, Locator, Readline, useraction