| Exponential covariances, radial basis functions and stationary covariances. {fields} | R Documentation |
Given two sets of locations computes the cross covariance matrix for covariances among all pairings.
exp.cov(x1, x2, theta = rep(1, ncol(x1)), p = 1, C = NA, marginal=FALSE)
exp.simple.cov(x1, x2, theta =1, C=NA,marginal=FALSE)
rad.cov(x1, x2, p = 1, with.log = TRUE, with.constant = TRUE,
C=NA,marginal=FALSE)
rad.simple.cov(x1, x2, p=1, with.log = TRUE, with.constant = TRUE,
C = NA, marginal=FALSE)
stationary.cov(x1, x2, Covariance="Exponential", Distance="rdist",
Dist.args=NULL, theta=1.0,C=NA, marginal=FALSE,...)
x1 |
Matrix of first set of locations where each row gives the coordinates of a particular point. |
x2 |
Matrix of second set of locations where each row gives the coordinates of a particular point. If this is missing x1 is used. |
theta |
Range (or scale) parameter. This can be a scalar, vector or matrix.
If a scalar or vector these are expanded to be the diagonal elements of
a linear transformation of the coordinates. In R code the transformation
applied before distances are found is: x1 %*% t(solve(theta)) or
if theta is a scalar: x1/theta.
Default is theta=1. See Details below.
|
C |
A vector with the same length as the number of rows of x2. If specified the covariance matrix will be multiplied by this vector. |
marginal |
If TRUE returns just the diagonal elements of the
covariance matrix using the x1 locations. In this case this is
just 1.0. The marginal argument will trivial for this function is a
required argument and capability for all covariance functions used with
Krig.
|
p |
Exponent in the exponential form. p=1 gives an exponential and p=2 gives a
Gaussian. Default is the exponential form.
For the radial basis function this is the exponent for the distance between locations. |
with.constant |
If TRUE includes complicated constant for radial basis functions.
See the function radbad.constant for more details. The
default is TRUE include the constant. Without the usual constant
the lambda used here will differ by a constant from estimators ( e.g.
cubic smoothing splines) that use the constant. Also a negative value
for the constant may be necessary to make the radial basis positive
definite as opposed to negative definite. |
with.log |
|
Covariance |
|
Distance |
Character string that is the name of the distance
function to use. Choices in fields are rdist, rdist.earth |
Dist.args |
A list of optional arguments to pass to the Distance function. |
... |
Any other arguments that will be passed to the
covariance function. e.g. smoothness for the Matern. |
For purposes of illustration, the function exp.cov.simple is
provided as a simple example and implements the R code discussed below.
It can also serve as a template for creating new covariance functions for the
Krig function. Also see the higher level function
stationary.cov to mix and match different covariance shapes and
distance functions.
Functional Form: If x1 and x2 are matrices where nrow(x1)=m and nrow(x2)=n then this function will return a mXn matrix where the (i,j) element is the covariance between the locations x1[i,] and x2[j,]. The covariance is found as exp( -(D.ij **p)) where D.ij is the Euclidean distance between x1[i,] and x2[j,] but having first been scaled by theta.
Specifically if theta is a matrix to represent a linear
transformation of the coordinates, then let
u= x1%*% t(solve( theta)) and v= x2%*% t(solve(theta)).
Form the mXn distance matrix with elements:
D[i,j] = sqrt( sum( ( u[i,] - v[j,])**2 ) ).
and the cross covariance matrix is found by exp(-D).
Note that if theta is a scalar then this defines an isotropic covariance
function and the functional form is essentially exp(-D/theta).
Implementation:
The function r.dist is a useful FIELDS function that finds
the cross Euclidean distance matrix (D defined above) for two sets of
locations. Thus in compact R code we have
exp(-rdist(u, v)**p)
Note that this function must also support two other kinds of calls:
If marginal is TRUE then just the diagonal elements are returned
(in R code diag( exp(-rdist(u,u)**p) )).
If C is passed then the returned value is
exp(-rdist(u, v)**p) %*% C
Radial basis functions: The functional form is
Constant* rdist(u, v)**p for odd dimensions
and Constant* rdist(u,v)**p * log( rdist(u,v)
For an m th order thin plate spline in d dimensions p= 2*m-d and must
be positive. The constant, depending on m and d, is coded in the fields
function radbas.constant. This form is only a generalized
covariance function – it is only positive definite when restricted to
linear subspace. See rad.simple.cov for a coding of the radial
basis functions in R code.
Stationary covariance: Here the computation is apply the function
Covariance to the distances found by the Distance function.
For example exp.cov(x1,x2, theta=MyTheta) and
stationary.cov( x1,x2, theta=MyTheta, Distance= "rdist",
Covariance="Exponential")
are the same. This also the same as
stationary.cov( x1,x2, theta=MyTheta, Distance= "rdist",
Covariance="Matern",smoothness=.5).
About the FORTRAN: The actual function exp.cov and
rad.cov calls FORTRAN to
make the evaluation more efficient this is especially important when the
C argument is supplied. So unfortunately the actual production code in
exp.cov is not as crisp as the R code sketched above. See
rad.simple.cov for a R coding of the radial basis functions.
If the argument C is NULL the cross covariance matrix. Moreover if x1 is
equal to x2 then this is the covariance matrix for this set of locations.
In general if nrow(x1)=m and nrow(
x2)=n then the returned matrix, Sigma will be mXn.
If C is a vector of length n,
then returned value is the multiplication of the cross covariance matrix
with this vector: Sigma%*%C
Krig, rdist, rdist.earth, gauss.cov, exp.image.cov, Exponential, Matern.
# exponential covariance matrix ( marginal variance =1) for the ozone #locations out<- exp.cov( ozone$x, theta=100) # out is a 20X20 matrix out2<- exp.cov( ozone$x[6:20,],ozone$x[1:2,], theta=100) # out2 is 15X2 matrix # Kriging fit where the nugget variance is found by GCV fit<- Krig( ozone$x, ozone$y, cov.function="exp.cov", theta=100)