|
|
3.4.2 Defining New Covariance Functions
Any covariance function can be passed into the Krig function provided
it has the right format. For example, suppose I wish to create a new
covariance function. First, I would define a function named, say, foo.cov.S.
Thus,
foo.cov.S <- function( x1, x2, range) {
exp( -(rdist(x1, x2)/range)**2)
} # end of foo.cov fcn (Note that foo.cov.S is the Gaussian covariance fcn.)
To use this covariance
function on the ozone data (using range = 10), and store it in an object called
foo.fit, simply use the following command.
foo.fit <- Krig( ozone$x, ozone$y, cov.function=foo.cov.S, range=10)
Notice that range is a parameter defined in my foo.cov.S function, and
it is passed into Krig as if it were a Krig parameter. The same kinds of
summary information can be obtained as before by using the usual summary
command.
summary( foo.fit)
Call:
Krig(x = ozone$x, Y = ozone$y, cov.function = foo.cov.S, range = 10)
Number of Observations: 20
Number of unique points: 20
Degree of polynomial null space ( base model): 1
Number of parameters in the null space 3
Effective degrees of freedom: 3.2
Residual degrees of freedom: 16.8
MLE sigma 4.402
GCV est. sigma 4.402
MLE rho 0.2765
Scale used for covariance (rho) 0.2765
Scale used for nugget (sigma^2) 19.38
lambda (sigma2/rho) 70.08
Cost in GCV 1
GCV Minimum 23.06
Residuals:
min 1st Q median 3rd Q max
-7.802 -2.736 -0.3941 2.757 7.472
REMARKS
Covariance function: foo.cov.S
One caveat to creating new covariance functions is that they must not contain an argument called "C" unless the C argument is used in a special.
Utilizing the C argument
The fastidious eye will have noticed that two exponential covariance
functions are included with the Fields software package. Indeed, this is the
case with many of the provided functions. There is no difference in the results
between using exp.cov and
exp.cov.S. However, for large data sets exp.cov
is more efficient than exp.cov.S if the actual
covariance matrix, Sigma, is not needed, but instead simply the result of right
multiplying the Sigma by a vector, say C.
For example, if Sigma = my.cov( x1, x2)
is the covariance matrix, and y is a vector of length =
nrow( x2) one requires the result
Sigma %*% y,
This particular computation can be invoked using the C argument. That is,
my.cov( x1, x2, C=y)
The advantage is that Krig will not create the covariance
matrix, Sigma. Instead, it will perform the calculations of Sigma %*% C using, in most
cases (true for all Fields provided covariance functions), a Fortran subroutine and save
only the result of the multiplication.
By contrast, the Fields function exp.cov.S (which
is the same as exp.cov except that it does not utilize the C
option) will create the matrix Sigma (thereby using up more memory) and then it performs
the calculation Sigma %*% C. Thus, when using large data sets,
exp.cov, is faster and can be more efficient than exp.cov.S.
|