Covariance returns the values of complex stationary and
nonstationary covariance functions;
see CovarianceFct and Covariance
return a vector of values of the covariance function.+
Operator that adds up at most 10 submodels*
Operator that multiplies at most 10 submodels$
$$C(x, y) = v C(x / s, y / s)$$$$C(x, y) = v C(x a, y a)$$$$C(x, y) = v C(A x, A y)$$$$C(x, y) = v C(p x, p y)$$Operator that modifies the the variance ($v=$var) and the
coordinates or distances byscale) oranisoTmultiplied from the
right orprojon a lower dimensional space along
the coordinate axisscaleis positive,anisoandAare matrices, andprojis a vector indices with between 1 and
the dimension of$x$.
Note, at most one of the parameters,anisoT,A,projmay be given at the same time. The operator$has 1 submodel.
If the dimension of the field is 1 oranisois not given, the
operator allows for derivatives.
ave1% bernoulli 2010 paper example 13$$C(h, u) = |E + 2Ahh^tA|^{-1/2}
\phi(\sqrt(\|h\|^2/ 2 + (z^th + u)^2 (1 - 2h^tA (E+2Ahh^tA)^{-1} Ah)))$$where$E$is the identity matrix.$A$is a symmetric positive definite$(d-1) \times (d-1)$and$z$is a$d-1$dimensional vector.
The function$\phi$is normal mixture model, e.g.whittlemodel, seeave2(nonstationary)
Here$C(h) = C_0(h, 0)$where$C_0$is theave1model.biWM(bivariate model)
$$C_{ij}(h) = c_{ij} W_{\nu_{ij}} ( h / s_{ij})$$where$W_nu$is thewhittlemodel and$i,j=1,2$.
For (i=j) the constants$\nu_{ii}, s_{ii}, c_{ii} > 0$.
For the offdiagonal elements with have$C_{12} = C_{21}$,$s_{12}=s_{21} > 0$,$\nu_{12} =\nu_{21} = 0.5 (\nu_{11} + \nu{22}) /
\nu_{red}$for some constant$\nu_{red} \in (0,1]$.
The scalar$c_{12} =c_{21} = \rho_{red} \sqrt{f m c_{11} c_{22}}$where$$f = \Gamma(\nu_{11} + d/2) * \Gamma(\nu_{22} + d/2) /
\Gamma(\nu_{11}) / \Gamma(\nu_{22}) *
(\Gamma(\nu_{12}) / \Gamma(\nu_{12}+d/2))^2 *
( s_{12}^{2*\nu_{12}} / s_{11}^{\nu_{11}} / s_{22}^{\nu_{22}} / )^2$$and$\Gamma$is the Gamma function and$d$is the dimension
of the space. The constant$m$is
the infimum of the function$g$on$[0,\infty)$,$$g(t) = (1/s_{12}^2 +t^2)^{2\nu_{12} + d}
(1/s_{11}^2 + t^2)^{-\nu_{11}-d/2}
(1/s_{22}^2 + t^2)^{-\nu_{22}-d/2}$$see the reference below for details on the infimum. The model now has the parametersnu$= (nu_{11}, nu_{22})$
nured12$=\nu_{red}$
s$= (s_{11}, s_{22})$
s12$= s_{12} = s_{21}$\c$= (c_{11}, c_{22})$
rhored$=\rho_{red}$See alsoparsbiWM.
constant
This model is designes for the use inXin modelcoxisham$$C(h, u) = |E + u^\beta D|^{-1/2} \phi([
(h - u \mu)^t (E + u^\beta D)^{-1} (h - u\mu)]^{1/2})$$Here $$mu$is vector;$E$is the identity matrix and$D$is a correlation
matrix with$|D| >0$. Currently implementation is done only for$d=2$. The parameter$\beta$is in$(0,2]$and equals 2 by default.curlfree(multivariate)$$( - \nabla_x \nabla_x^T ) C_0(x, t)$$$C_0$is a univariate
covariance model that is motion invariant and at least twice differentiable
in the first component. The operator is applied to the first
component only. The model returns the potential field in the first
component, the corresponding
curlfree field and field of sources and sinks in the last component.
The above formula for the covariance function only gives the part
for the curlfree field. The complete matrix-valued
correlation function, including all
components, is more complicated.$C_0$is either a spatiotemporal model (then$t$is the
time component) or it is an isotropic model. Then, the first$Dspace$coordinates are considered as$x$coordinates and the remaining
ones as$t$coordinates. By default,$Dspace$equals the
dimension of the field (and$t$is identically$0$). See also the modelsdivfreeandvector.
cutoff$$C(h)=\phi(h), 0\le h \le d$$$$C(h) = b_0 ((dr)^a - h^a)^{2 a}, d \le h \le dr$$$$C(h) = 0, dr\le h$$The cutoff model is a functional of the covariance function$\phi$.
Here,$d>0$should be the diameter of the domain on which
simulation is done . The parameter$a>0$has been shown to be
optimal for$a = 1/2$or$a =1$.
The parameters$r$and$b_0$are chosen internally such that$C$is a smooth function. NOTE: 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,whittleandgencauchy, some sufficient conditions
are known.
delayeffect(bivariate)$$C_{11}(h) = C_{22}(h) = C_0(h) \qquad C_{12}(h) = C_0(h +
r), C_{21}(h) = C_0(-h + r)$$Here$r$is a vector of the dimension of the random field,
and$C_0$is a translation invariant, univariate covariance model.divfree(multivariate)$$( - \Delta E + \nabla \nabla^T ) C_0(x, t)$$$C_0$is a univariate
covariance model that is motion invariant and at least twice differentiable
in the first component. The operator is applied to the first
component only. The model returns the potential field in the first
component, the corresponding
divfree field and the field of curl strength in the last component.
The above formula for the covariance function only gives the part
for the divfree field. The complete matrix-valued
correlation function, including all
components, is more complicated.$C_0$is either a spatiotemporal model (then$t$is the
time component) or it is an isotropic model. Then, the first$Dspace$coordinates are considered as$x$coordinates and the remaining
ones as$t$coordinates. By default,$Dspace$equals the
dimension of the field (and$t$is identically$0$). See also the modelscurlfreeandvector.
EtAxxA(auxiliary function)$$S(x) = E + R^t A^t x x^t A R, \qquad x \in R^3$$where$E$and$A$are arbitrary$3 \times 3$matrices
and$R$is a rotation matrix,$$R =\left(
\begin{array}{lll}
\cos (\alpha x_3) & -\sin(\alpha x_3)& 0
\
\sin(\alpha x_3)& \cos (\alpha x_3) & 0
\
0 & 0 & 1
\end{array}
\right)$$This is not a covariance function, but can be used as a submodel for
certain classes of non-stationary covariance functions.Exp$$C(h) = \exp(-\gamma(h))$$where$\gamma$is a valid variogram.
If a stationary covariance model$C$is given in stead of$\gamma$,
this is automatically turned into a variogram model, i.e.$C(h) = \exp(-C(0)+C(h))$.%%%%%%%%%%%%%%%%%%%%%%%%%%%
M$$C(h) = M^t \phi(h) M$$Here$phi$is a$k$-variate variogram or covariance,
and$M$is any$m \times k$matrix.ma1$$C(h) = (\theta / (1 - (1-\theta) * C_0(h)))^\alpha$$Here,$C_0$is any correlation function,$\alpha \in (0,\infty)$and$\theta \in (0,1)$.ma2$$C(h) = (1 - exp(-\gamma(h)))/\gamma(h)$$Here$\gamma$is a variogram model.mastein$$C(h, t)=\frac{\Gamma(\nu + \gamma(t))\Gamma(\nu + \delta)}{
\Gamma(\nu +
\gamma(t) + \delta) \Gamma(\nu)}
W_{\nu + \gamma(t)}(\|h - Vt\|)$$$\Gamma$is the Gamma function;$\gamma(t)$is a variogram on the real axis;$W$is the Whittle-Matern model.
Here, the
names of covariance models can also be used; the algorithm
chooses the corresponding variograms then.
The parameter$\nu$is the smoothness parameter
of the Whittle-Matern model (for$t=0$) and must be positive.
Finally,$\delta$must be greater than or equal to half the
dimension of$h$.
Instead of the velocity parameter$V$in original model
description, a preceeding anisotropy matrix is chosen appropriately:$$\left( \begin{array}{cc} A & -V \ 0 &
1\end{array}\right)$$A is a spatial transformation matrix.
(I.e. (x,t) is multiplied from left on the above matrix and
the first elements of the obtained vector are intepreted as
new spatial components and only these components are used to form
the argument in the Whittle-Matern function.)
The last component in the new coordinates is the time which is
passed to$\gamma$. (Velocity is assumed to be zero in
the new coordinates.)
%Note that
%the anisotropy matrix must be such that \eqn{(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 \eqn{\gamma}{gamma}
%is naught except the very last, purely temporal one.Note, that for numerical reasons,$\nu+\gamma+d$may not exceed the value 80.0. If exceeded the algorithm fails.
mixedThis model is designed for the use inXis a matrix of independent variables.
The second,b, is a vector of regression coefficients.
Furthermore a submodel,covb, may give the covariance structure
forb. Letnthe number of (non-repeated) observations.
The following combinations are allowed:
Xis given. ThenXis a scalar
or a vector of lengthn, andXdefines a known mean.Xandbare given. ThenXis a$(n \times m)$matrix
where$m$is the length of the vector$b$. Then a fixed
effect is defined.Xandcovbare given.covbis the modelconstant,
then we have a random model (maybe with preceeding model$).covbis any other model then we have
a geoadditive partNAs, but notX.mqam(multivariate quasi-arithmetic mean)$$C_{ij}(h) = \rho_{ij} \phi(\theta \phi^{-1}(C_i(h)) + (1 - \theta) \phi^{-1}(C_j(h)))$$where$\phi$is a completely monotone function
and$C_i$are suitable covariance functions.
The submodel$\phi$is given (by name) as first submodel.
Since$\phi$is completely monotone if and only if$\phi(\| . \|^2)$is a valid
covariance function for all dimensions, e.g.stable,gauss,exponential,$\phi$is given
by the name of the corresponding covariance function$C$,
i.e.$phi( . ) = C(sqrt( . ))$. Warning:RandomFieldscannot check whether the combination
of$\phi$and$C_i$is valid.
natsc$$C(h) = C_0(h / s)$$Where$C_0$is any stationary and isotropic model. The parametersis chosen bynatscsuch that the practical range
(or the mathematical range, if finite) is 1.nonstWM$$C(x, y)=\Gamma(\mu) \Gamma(\nu(x))^{-1/2} \Gamma(\nu(y))^{-1/2}
W_{\mu} (\|x-y\|)$$$$= 2^{1-\mu} \Gamma(\nu(x))^{-1/2} \Gamma(\nu(y))^{-1/2}
\|x-y\|^\mu
K_\nu(\|x-y\|)$$where$\mu = [\nu(x) + \nu(y)]/2$and$\nu$is a positive function.
If$\nu$is a scalar use the variablenu.
If$\nu$is a function, use the submodelNu.
Note that forNuthe usual list structure applies
and only the defined covriance models can be used.nsst(Non-Separable Space-Time model)$$C(h,u)= (\psi(u)+1)^{-\delta/2} \phi(h /\sqrt(\psi(u) +1))$$The parameter$\delta$must be greater than or equal to
the spatial dimension of the field.$\phi$is normal mixture model and$\psi$is a variogram.
This model is used for space-time modelling where the spatial
component is isotropic.
nugget(multivariat model)$$C(h)= {\rm diag}(1,\ldots,1) 1_{{0}}(h)$$The components of the multivariate vector are always independent.
The models adapts the multivariate dimension to the calling model.parsbiWM(bivariate model)
$$C_{ij}(h) = c_{ij} W_{\nu_{ij}} (h / s)$$where$W_nu$is thewhittlemodel and$i,j=1,2$.
For (i=j) the constants$\nu_{ii}, c_{ii} \ge 0$and$s>0$.
For the offdiagonal elements with have$C_{12} = C_{21}$.
Furthermore,$\nu_{12} =\nu_{21} = 0.5 (\nu_{11} + \nu{22})$and the scalar$c_{12} =c_{21} = \rho_{red} \sqrt{f m c_{11} c_{22}}$where$$f = \Gamma(\nu_{11} + d/2) * \Gamma(\nu_{22} + d/2) /
\Gamma(\nu_{11}) / \Gamma(\nu_{22}) *
(\Gamma(\nu_{12}) / \Gamma(\nu_{12}+d/2))^2$$and$\Gamma$is the Gamma function and$d$is the dimension
of the space. The constant$m$is
the infimum of the function$g$on$[0,\infty)$,$$g(t) = (1/s_{12}^2 +t^2)^{2\nu_{12} + d}
(1/s_{11}^2 + t^2)^{-\nu_{11}-d/2}
(1/s_{22}^2 + t^2)^{-\nu_{22}-d/2}$$see the reference below for details on the infimum. The model now has the parametersnu$= (\nu_{11}, \nu_{22})$
s$= (s_{11}, s_{22})$
s12$= s_{12} = s_{21}$
c$= (c_{11}, c_{22})$
rhored$=\rho_{red}$See alsobiWM.
Pow$$\gamma(h) = (\gamma_0(h))^\alpha$$or$$C(h) = C_0(0) - [C_0(0) - C_0(h)]^\alpha$$where$\gamma_0$is a valid variogram or$C_0$is a valid covariance function,
and$\alpha \in [0, 1]$.qam(Quasi-arithmetic mean)$$C(h) = \phi(\sum_i \theta_i \phi^{-1}(C_i(h)))$$where$\phi$is a completely monotone function
and$C_i$are suitable covariance functions.
The submodel$\phi$is given (by name) as first submodel.
Since$\phi$is completely monotone if and only if$\phi(\| . \|^2)$is a valid
covariance function for all dimensions, e.g.stable,gauss,exponential,$\phi$is given
by the name of the corresponding covariance function$C$,
i.e.$phi( . ) = C(sqrt( . ))$. Warning:RandomFieldscannot check whether the combination
of$\phi$and$C_i$is valid.
rational(auxiliary)$$S(x) = (a_0 + a_1 * x^t A A^t x)/ (1 + x^t A A^t x)$$where is some$d\times d$matrix
and$a = (a_0, a_1)$is a 2-dimensional vector.Rotat(auxiliary function)$$S^t(x) = x^t R, \qquad x \in R^3$$where
and$R$is a rotation matrix,$$R =\left(
\begin{array}{lll}
\cos (\alpha x_3) & -\sin(\alpha x_3)& 0
\
\sin(\alpha x_3)& \cos \alpha x_3 & 0
\
0 & 0 & 1
\end{array}
\right)$$This is not a covariance function, but can be used a submodel for
certain classes of non-stationary covariance functions.Stein$$C(h)=a_0 + a_2 (h)^2 + \phi(h), 0\le h \le D$$$$C(h)=b_0 (rD - h)^3/(h), r \le h \le rD$$$$C(h) = 0, rD \le h$$The Stein model is a functional of the covariance function$\phi$.
Here,$D>0$should be the diameter of the domain on which
simulation is done,$r\ge1$.
The parameters$a_0$,$a_2$and$b_0$are chosen internally such that$C$becomes a smooth function.
% Note that the scale parameter for the submodel must be \code{double(0)}. NOTE: 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
modelfractalBsome sufficient conditions are known.
steinst1(non-separabel space time model)$$C(h, t) = W_{\nu}(y) -
\frac{\langle h, z \rangle t}{(\nu - 1)(2\nu
+ d)} W_{\nu -1}(y)$$Here,$W_{\nu}$is the Whittle-Matern model with
smoothness parameter$\nu$;$y=\|(h,t)\|$.$z$is a vector whose norm
must less than or equal to 1.stp$$C(x,y) = |S_x|^{1/4} |S_y|^{1/4} |A|^{-1/2}
\phi(Q(x,y)^{1/2})$$where$$Q(x,y) = c^2 - m^2 + h^t (S_x + 2(m + c)M) A^{-1} (A_y + 2
(m-c)M)h,$$$$c = -z^t h + \xi_2(x) - \xi_2(y),$$$$A = S_x + S_y + 4 M h h^t M$$$$m = h^t M h
.$$$$h = H(x) - H(y)$$The parameters aretbm2$$C(h) = \frac{d}{dh}\int_0^h \frac{u\phi(u)}{\sqrt{h^2 - u^2}} d u$$for some stationary and isotropic covariance$\phi$that is valid in at least 2 dimensions.This operator is currently only designed for internal use!
tbm3$$C(h) = \phi(h) + h \phi'(h) / n$$which, forn=1reduced to the standard TBM operator$$C(h) =\frac {d}{d h} h \phi(h)$$for some stationary and isotropic covariance$\phi$that is valid in at least$n+2$dimensions.nshould be an integer.
This operator is currently only designed for internal use!vector(multivariate)$$( -0.5 * (a + 1) \Delta E + a \nabla \nabla^T ) C_0(x, t)$$$C_0$is a univariate
covariance model that is motion invariant and at least twice differentiable
in the first component. The operator is applied to the first
component only. The parameter$a$is in$[-1, 1]$. If$a=-1$then the field is curl free; if$a=1$then the field is
divergence free.$C_0$is either a spatiotemporal model (then$t$is the
time component) or it is an isotropic model. Then, the first$Dspace$coordinates are considered as$x$coordinates and the remaining
ones as$t$coordinates. By default,$Dspace$equals the
dimension of the field (and$t$is identically$0$). See also the modelsdivfreeandcurlfree
()$PracticalRangeis usually not defined for the above modelsmodel="name"andparam=c(mean, variance, nugget, scale,...)."+", "*", or Gneiting's
"nsst". Consequently, we need also an operator, called
"$", that changes the variance and the scale.E.g. a standard exponential model (variance=1, scale=1, nugget=0) is now simply written as $$\hbox{list("exponential")}$$
(And no param must be given!)
Further, a standard exponential model with a nugget effect,
nugget variance 3, is now written as
list("+",
list("exponential"),
list("$", var=3, list("nugget"))
)
Here, only the relevant parameters need to be given; the missing
parameters get standard values whenever standard values exist,
e.g. variance equals 1 if not given.
Further, the parameters can (and must) be called by names, which makes
complex models much more readable.
Submodels, as list("exponential") in the second example above,
can (but need not) be called by name.
CovarianceFctave1, ave2
biWM, parsbiWM
coxisham
curlfree
delayeffect
divfree
vector
Ma-Stein model
mixed
nsst
Quasi-arithmetic means (qam, mqam)
Stein
tbm
RandomFields,
PrintModelList(op=TRUE)
## the subsequent model can be used to model rainfall...
y <- x <- seq(0, 10, len=25) # better 256 -- but will take a while
T <- c(0, 10, 1) # better 0.1
col <- c(topo.colors(300)[1:100], cm.colors(300)[c((1:50) * 2, 101:150)])
model <- list("coxisham", mu=c(1, 1), D=matrix(nr=2, c(1, 0.5, 0.5, 1)),
list("whittle", nu=1)
)
system.time(z <- GaussRF(x, y, T=T, grid =TRUE, spectral.lines=1500,
model = model))
zlim <- range(z)
time <- dim(z)[3]
for (i in 1:time) {
Print(i)
sleep.milli(100)
image(x, y, z[, , i], add=i>1, col=col, zlim=zlim)
}
####################################################
####################################################
# the following five model definitions are the same!
## (1) very traditional form
(cv <- CovarianceFct(x, model="bessel", param=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=rbind(c(2, 5, 0.5), c(1, 0, 0)))
#### most general notation in form of lists
## (4) isotropic notation
model <- list("+",
list("$", var=2, scale=5, list("bessel", 0.5)),
list("nugget"))
cv - CovarianceFct(x, model=model)
## (5) anisotropic notation
model <- list("+",
list("$", var=2, aniso=0.2, list("bessel", 0.5)),
list("nugget"))
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' should now be coded as
gneitingdiff <- function(p){
list("+",
list("$", var=p[3], list("nugget")),
list("$", scale=p[4],
list("*",
list("$", var=p[2], scale=p[6], list("gneiting")),
list("whittle", nu=p[5])
)
)
)
}
# and then
param <- c(NA, runif(5, max=10))
CovarianceFct(0:100, model=gneitingdiff(param))
## instead of formerly CovarianceFct(x,"gneitingdiff",param)Run the code above in your browser using DataCamp Workspace