Derivatives of sample paths:
In d dimensions a Gaussian process with Matern covariance and smoothness
parameter
nu will have nu + d/2 derivatives that exist in a mean square sense.
The marginal variance is given by the parameter scale. When this
parameter is set to one this yields a correlation
function. There several different ways to parameterize the Matern family and the
reader is referred to Stein's book page 49 for discussion. In terms
of the more geostatistical terminology, we note that out theta is also the
"range" and the scale is also the "sill" if there is no nugget variance
included in the covariance. We caution that the range for this function
gives a qualitatively different scaling as one varies the smoothness.
Functional Form: If x1 and x2 are matrices where nrow(x1)=m and nrow(
x2)=n and each row are the coordinates of a location, then this
function should 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 H( D.ij) where D.ij is the Euclidean
distance between x1[i,] and x2[j,] but having first been scaled by theta.
H is proportional to a modified Bessel function
of third kind using denoted by K.nu . In our parameterization we take
smoothness = nu and H is normalized so that H(0)=1. (See the function matern for a succinct definition as R code.)
The reader is referred
to Stein's book, page 31 for more details. Note that we do not use Stein's
normalization, however, as it seems more useful to normalize the
covariance so that .
Definition of the distance matrix:
D.ij = sqrt( sum.k (( x1[i,k] - x2[j,k]) /theta[k])**2 ).
Note that if theta is a scalar then this defines an isotropic covariance
function.
Implementation: The function rdist is a useful FIELDS function that finds
the cross
distance matrix ( D defined above) for two sets of locations. Thus in
compact S code we have
u <- t(solve(theta)v <- t(solve(theta)H(-rdist(u, v))
where solve(theta) is the (matrix) inverse for theta.
A simple modification of this function for the user would be to substitute
rdist.earth for rdist to give a distance metric that makes sense for
a small region in lon/lat coordinates.