Sophisticated Models {RandomFields} | R Documentation |
Covariance
returns the values of complex stationary and
nonstationary covariance functions;
see CovarianceFct
for basic isotropic models
Here only the non-isotropic and hyper models are listed;
see CovarianceFct
for basic isotropic models.
The implemented models are in standard notation for a covariance function (variance 1, nugget 0, scale 1) and for positive real arguments h (and t) for the stationary models or parts:
+
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 by
the scale (s=scale
) or
the anisotropy matrix a=anisoT
multiplied from the
right or
the anisotropy matrix A multiplied from the left or
proj
on a lower dimensional space along
the coordinate axis
The parameter scale
is positive, aniso
and A
are matrices, and proj
is a vector indices with between 1 and
the dimension of x.
Note, at most one of the parameters,
anisoT
, A
, proj
may be given at the same time.
The operator $
has 1 submodel.
If the dimension of the field is 1 or aniso
is not given, the
operator allows for derivatives.
ave1
C(h, u) = |E + 2 A h h^t A|^{-1/2} phi(sqrt[|h|^2/2 + (z^t h + u)^2 (1 - 2h^t A (E + 2 A h h^t A)^{-1} A h)])
where E is the identity matrix.
A is a symmetric positive definite (d-1)
x (d-1) and z is a d-1 dimensional vector.
The function φ is normal mixture model, e.g.
whittle
model, see CovarianceFct
and PrintModelList
().
ave2
(nonstationary)
Here C(h) = C_0(h, 0) where C_0 is the
ave1
model.
biWM
(bivariate model)
C_{ij}(h) = c_{ij} W_{ν_{ij}} ( h / s_{ij})
where W_nu is the whittle
model and i,j=1,2.
For (i=j) the constants ν_{ii}, s_{ii}, c_{ii} > 0.
For the offdiagonal elements with have C_{12} = C_{21},
s_{12}=s_{21} > 0,
ν_{12} =ν_{21} = 0.5 (ν_{11} + ν{22}) /
ν_{red}
for some constant ν_{red} \in (0,1].
The scalar c_{12} =c_{21} = ρ_{red} √{f m c_{11} c_{22}}
where
f = Γ(ν_{11} + d/2) * Γ(ν_{22} + d/2) / Γ(ν_{11}) / Γ(ν_{22}) * (Γ(ν_{12}) / Γ(ν_{12}+d/2))^2 * ( s_{12}^{2*ν_{12}} / s_{11}^{ν_{11}} / s_{22}^{ν_{22}} / )^2
and Γ is the Gamma function and d is the dimension of the space. The constant m is the infimum of the function g on [0,∞),
g(t) = (1/s_{12}^2 +t^2)^{2ν_{12} + d} (1/s_{11}^2 + t^2)^{-ν_{11}-d/2} (1/s_{22}^2 + t^2)^{-ν_{22}-d/2}
see the reference below for details on the infimum.
The model now has the parameters
nu
= (nu_{11}, nu_{22})
nured12
=ν_{red}
s
= (s_{11}, s_{22})
s12
= s_{12} = s_{21}\
c
= (c_{11}, c_{22})
rhored
=ρ_{red}
See also parsbiWM
.
constant
This model is designes for the use in fitvario
as a part of a linear model definition.
Its only parameter is a covariance matrix of appropriate size
to match the number of (non-repeated) observations or
the number of columns of parameters X
in model
mixed
, see sophisticated.
coxisham
C(h, u) = |E + u^β D|^{-1/2} φ([ (h - u μ)^t (E + u^β D)^{-1} (h - uμ)]^{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 β is in (0,2] and equals 2 by default.
curlfree
(multivariate)
( - \nabla_x \nabla_x^T ) C_0(x, t)
C_0is 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_0is 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 models divfree
and vector
.
cutoff
C(h)=φ(h), 0≤ h ≤ d
C(h) = b_0 ((dr)^a - h^a)^{2 a}, d ≤ h ≤ dr
C(h) = 0, dr≤ 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
,
whittle
and gencauchy
, some sufficient conditions
are known.
delayeffect
(bivariate)
C_{11}(h) = C_{22}(h) = C_0(h), 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)
( - Δ E + \nabla \nabla^T ) C_0(x, t)
C_0is 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_0is 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 models curlfree
and vector
.
EtAxxA
(auxiliary function)
S(x) = E + R^t A^t x x^t A R, x \in R^3
where E and A are arbitrary 3 x 3 matrices and R is a rotation matrix,
R = matrix(c(cos(alpha x_3), -sin(alpha x_3), 0 sin(alpha x_3), cos(alpha x_3), 0, 0, 0, 1), nc=3, by=TRUE)
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(-γ(h))
where γ is a valid variogram. If a stationary covariance model Cis given in stead of γ, this is automatically turned into a variogram model, i.e. C(h) = \exp(-C(0)+C(h)) .
M
C(h) = M^t φ(h) M
Here phi is a k-variate variogram or covariance, and M is any m x k matrix.
ma1
C(h) = (θ / (1 - (1-θ) * C_0(h)))^α
Here, C_0 is any correlation function, α \in (0,∞) and θ \in (0,1).
ma2
C(h) = (1 - exp(-γ(h)))/γ(h)
Here γ is a variogram model.
mastein
C(h, t)= [ Γ(ν + γ(t))Γ(ν + δ) ] / [ Γ(ν + γ(t) + δ) Γ(ν) ] W_{ν + γ(t)}(\|h - Vt\|)
Γis the Gamma function; γ(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 ν is the smoothness parameter of the Whittle-Matern model (for t=0) and must be positive. Finally, δ 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:
matrix(c(A, -V, 0, 1), nr=2, by=TRUE)
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 g. (Velocity is assumed to be zero in the new coordinates.)
Note, that for numerical reasons, ν+γ+d may not exceed the value 80.0. If exceeded the algorithm fails.
mixed
This model is designed for the use in fitvario
to build up linear regression models with fixed effects, mixed
effects, including geoadditive parts.
The model has two parameters. The first, X
is a matrix of independent variables.
The second, b
, is a vector of regression coefficients.
Furthermore a submodel, covb
, may give the covariance structure
for b
.
Let n
the number of (non-repeated) observations.
The following combinations are allowed:
only X
is given. Then X
is a scalar
or a vector of length n
, and X
defines a known mean.
X
and b
are given. Then X
is a
n x m matrix
where m is the length of the vector b. Then a fixed
effect is defined.
X
and covb
are given.
if covb
is the model constant,
then we have a random model (maybe with preceeding model
$
).
if covb
is any other model then we have
a geoadditive part
The data in the fitvario may contain NA
s, but not X
.
mqam
(multivariate quasi-arithmetic mean)
C_{ij}(h) = ρ_{ij} φ(θ φ^{-1}(C_i(h)) + (1 - θ) φ^{-1}(C_j(h)))
where φ is a completely monotone function and C_i are suitable covariance functions.
The submodel φ is given (by name) as first submodel.
Since φ is completely monotone if and only if
φ( \| \cdot \|^2) is a valid
covariance function for all dimensions, e.g. stable
,
gauss
, exponential
, φ is given
by the name of the corresponding covariance function C,
i.e. φ(\cdot) = C(√(\cdot)).
Warning: RandomFields
cannot check whether the combination
of φ and C_i is valid.
natsc
C(h) = C_0(h / s)
Where C_0 is any stationary and isotropic model. The parameter
s
is chosen by natsc
such that the practical range
(or the mathematical range, if finite) is 1.
nonstWM
C(x, y)=Γ(μ) Γ(ν(x))^{-1/2} Γ(ν(y))^{-1/2} W_{μ} (\|x-y\|)
= 2^{1-μ} Γ(ν(x))^{-1/2} Γ(ν(y))^{-1/2} \|x-y\|^μ K_ν(\|x-y\|)
where μ = [ν(x) + ν(y)]/2 and
ν is a positive function.
If ν is a scalar use the variable nu
.
If ν is a function, use the submodel Nu
.
Note that for Nu
the usual list structure applies
and only the defined covriance models can be used.
nsst
(Non-Separable Space-Time model)
C(h,u)= (ψ(u)+1)^{-δ/2} φ(h /√(ψ(u) +1))
The parameter δ must be greater than or equal to the spatial dimension of the field. φ is normal mixture model and ψ is a variogram.
This model is used for space-time modelling where the spatial
component is isotropic.
nugget
(multivariat model)
diag(1,…,1) 1(h==0)
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_{ν_{ij}} (h / s)
where W_nu is the whittle
model and i,j=1,2.
For (i=j) the constants ν_{ii}, c_{ii} ≥ 0 and s>0.
For the offdiagonal elements with have C_{12} = C_{21}.
Furthermore, ν_{12} =ν_{21} = 0.5 (ν_{11} + ν{22})
and the scalar c_{12} =c_{21} = ρ_{red} √{f m c_{11} c_{22}}
where
f = Γ(ν_{11} + d/2) * Γ(ν_{22} + d/2) / Γ(ν_{11}) / Γ(ν_{22}) * (Γ(ν_{12}) / Γ(ν_{12}+d/2))^2
and Γ is the Gamma function and d is the dimension of the space. The constant m is the infimum of the function g on [0,∞),
g(t) = (1/s_{12}^2 +t^2)^{2ν_{12} + d} (1/s_{11}^2 + t^2)^{-ν_{11}-d/2} (1/s_{22}^2 + t^2)^{-ν_{22}-d/2}
see the reference below for details on the infimum.
The model now has the parameters
nu
= (ν_{11}, ν_{22})
s
= (s_{11}, s_{22})
s12
= s_{12} = s_{21}
c
= (c_{11}, c_{22})
rhored
=ρ_{red}
See also biWM
.
Pow
γ(h) = (γ_0(h))^α
or
C(h) = C_0(0) - [C_0(0) - C_0(h)]^α
where γ_0 is a valid variogram or C_0 is a valid covariance function, and α \in [0, 1].
qam
(Quasi-arithmetic mean)
C(h) = φ(∑_i θ_i φ^{-1}(C_i(h)))
where φ is a completely monotone function and C_i are suitable covariance functions.
The submodel φ is given (by name) as first submodel.
Since φ is completely monotone if and only if
φ( \| \cdot \|^2) is a valid
covariance function for all dimensions, e.g. stable
,
gauss
, exponential
, φ is given
by the name of the corresponding covariance function C,
i.e. φ(\cdot) = C(√(\cdot)).
Warning: RandomFields
cannot check whether the combination
of φ 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 x d matrix and a = (a_0, a_1) is a 2-dimensional vector.
Rotat
(auxiliary function)
S^t(x) = x^t R, x \in R^3
where and R is a rotation matrix,
R = matrix(c(cos(alpha x_3), -sin(alpha x_3), 0 sin(alpha x_3), cos(alpha x_3), 0, 0, 0, 1), nc=3, by=TRUE)
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 + φ(h), 0≤ h ≤ D
C(h)=b_0 (rD - h)^3/(h), r ≤ h ≤ rD
C(h) = 0, rD ≤ 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≥1. The parameters a_0, a_2 and b_0 are chosen internally such that C becomes 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 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(h, t) = W_ν(y) - <h, z> t W_{ν-1}(y) / [ (ν-1) (2ν + d +1) ]
Here, W_{ν} is the Whittle-Matern model with smoothness parameter ν; 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} φ(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 + ξ_2(x) - ξ_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 are
(strictly) positive definite matrices for x \in R^d
an arbitrary d x d matrix
arbitrary
arbitrary d-variate function on R^d
arbitrary univariate function on R^d
a normal mixture model
The model allows for mimicking cyclonic behaviour.
tbm2
C(h) = d/dh int_0^h [ u phi(u) ] / [ sqrt{h^2 - u^2} ] d u
for some stationary and isotropic covariance φ 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, for n=1
reduced to the standard TBM operator
C(h) = d/dh [ h phi(h) ]
for some stationary and isotropic covariance φ
that is valid in at least n+2 dimensions.
n
should be an integer.
This operator is currently only designed for internal use!
vector
(multivariate)
( -0.5 * (a + 1) Δ E + a \nabla \nabla^T ) C_0(x, t)
C_0is 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_0is 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 models divfree
and curlfree
See CovarianceFct
for comments on the use of a
covariance model.
However, for the above sophicated models, the following differences should be considered:
RFparameters
()$PracticalRange
is usually not defined for the above models
only the list notation
can be used, but not the simple model definitions
with model="name"
and
param=c(mean, variance, nugget, scale,...)
.
the use of Covariance
is obligatory if the model is
non-stationary.
the anisotropy matrix belonging to a hypermodel is applied first to the coordinates before any call of the submodels.
To use the above models, a new, very flexible, straight forward
list notation is needed. Background of this notation is that
we have ‘primitives’, i.e. functions that are positive
definite. And we have ‘operators’, i.e. functionals that
make out of given variograms, covariance functions etc. new
models. Examples are "+"
, "*"
, 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
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.
CovarianceFct
and Covariance
return a vector of values of the covariance function.
Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/~schlather
Overviews:
see reference list in CovarianceFct
ave1, ave2
Schlather, M. (2010) On some covariance models based on normal scale mixtures. Bernoulli, 16, 780-797. (Example 13)
biWM, parsbiWM
Gneiting, T., Kleiber, W., Schlather, M. (2011) Matern covariance functions for multivariate random fields JASA
coxisham
Cox, D.R., Isham, V.S. (1988) A simple spatial-temporal model of rainfall. Proc. R. Soc. Lond. A, 415, 317-328.
Schlather, M. (2010) On some covariance models based on normal scale mixtures. Bernoulli, 16, 780-797.
curlfree
see vector
cutoff
Gneiting, T., Sevecikova, H, Percival, D.B., Schlather M., Jiang Y. (2006) Fast and Exact Simulation of Large Gaussian Lattice Systems in $R^2$: Exploring the Limits. J. Comput. Graph. Stat. 15, 483-501.
Stein, M.
delayeffect
Wackernagel, H. (2003) Multivariate Geostatistics. Berlin: Springer, 3nd edition.
divfree
see vector
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.
vector
Fuselier, E.J. (2006) Refined Error Estimates for Matrix-Valued Radial Basis Functions PhD thesis. Texas A&M University
Scheuerer, M. and Schlather, M. (2011) Covariance Models for Random Vector Fields Submitted
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.
ma1/ma2
mixed
Ober, U., Erbe, M., Porcu, E., Schlather, M. and Simianer, H. (2011) Kernel-Based Best Linear Unbiased Prediction with Genomic Data. Submitted.
nonstWM/hyperbolic/cauchy
Stein, M. (2005) Nonstationary Spatial Covariance Functions. Tech. Rep., 2005
nsst
Gneiting, T. (1997) Normal scale mixtures and dual probability densitites, J. Stat. Comput. Simul. 59, 375-384.
Gneiting, T. (2002) 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.
Zastavnyi, V. and Porcu, E. (2011) Caracterization theorems for the Gneiting class space-time covariances. Bernoulli, ??.
Schlather, M. (2010) On some covariance models based on normal scale mixtures. Bernoulli, 16, 780-797.
Quasi-arithmetic means (qam, mqam)
Porcu, E., Mateu, J. & Cchristakos, G. (2007) Quasi-arithmetic means of covariance functions with potential applications to space-time data. Submitted to Journal of Multivariate Analysis.
Paciorek-Stein (steinst1)
Stein, M. (2005) Nonstationary Spatial Covariance Functions. Tech. Rep., 2005
Paciorek, C. (2003) Nonstationary Gaussian Processes for Regression and Spatial Modelling, Carnegie Mellon University, Department of Statistics, PhD thesis.
Stein
Stein, M.
stp
Schlather, M. (2008) On some covariance models based on normal scale mixtures. Submitted
tbm
Gneiting, T. (1999) On the derivatives of radial positive definite function. J. Math. Anal. Appl, 236, 86-99
Matheron, G. (1973). The intrinsic random functions and their applications. Adv . Appl. Probab., 5, 439-468.
CovarianceFct
,
EmpiricalVariogram
,
GetPracticalRange
,
parameter.range
,
RandomFields
,
RFparameters
,
ShowModels
.
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)