# fields, Tools for spatial data
# Copyright 2004-2011, Institute for Mathematics Applied Geosciences
# University Corporation for Atmospheric Research
# Licensed under the GPL -- www.gpl.org/licenses/gpl.html
"gcv.Krig" <- function(out, lambda.grid = NA, cost = 1,
nstep.cv = 200, rmse = NA, verbose = FALSE, tol = 1e-05,
offset = 0, y = NULL, give.warnings = TRUE, give.warnings.REML = FALSE) {
nt <- out$nt
np <- out$np
N <- out$N
D <- out$matrices$D
# Yet another monster function called by Krig
# but there just many simple steps ...
#
# if a y data vector is not supplied then
# use the one in the Krig object
if (is.null(y)) {
u <- out$matrices$u
shat.pure.error <- out$shat.pure.error
pure.ss <- out$pure.ss
}
else {
#with new data need to update some statistics.
out2 <- Krig.make.u(out, y = y)
u <- out2$u
shat.pure.error <- out2$shat.pure.error
pure.ss <- out2$pure.ss
}
if (verbose) {
cat("u used:", fill = TRUE)
print(u)
}
#
# generate a reasonable grid of lambda based on equally spaced
# effective degrees of freedom
if (is.na(lambda.grid[1])) {
temp.df <- seq(nt, (np - offset) * 0.95, , nstep.cv)
temp.df[1] <- temp.df[1] + 0.001
for (k in 1:nstep.cv) {
lambda.grid[k] <- Krig.df.to.lambda(temp.df[k], D)
}
}
# make sure that the grid is ordered.
lambda.grid <- sort(lambda.grid)
nl <- length(lambda.grid)
nd <- length(D)
V <- V.model <- V.one <- lplike <- trA <- shat <- rep(NA,
nl)
Dl <- rep(NA, nd)
#
# this is small little list used to pass information to the
# objective functions
info <- list(matrices = list(D = D, u = u), N = N, nt = nt,
cost = cost, pure.ss = pure.ss, shat.pure.error = shat.pure.error,
offset = offset)
#
# loop over lambda values for the grid search
for (k in 1:nl) {
#
# all the wonderful things calculated for each lambda
# note the use of the info list.
V[k] <- Krig.fgcv(lambda.grid[k], info)
V.one[k] <- Krig.fgcv.one(lambda.grid[k], info)
V.model[k] <- Krig.fgcv.model(lambda.grid[k], info)
lplike[k] <- Krig.flplike(lambda.grid[k], info)
shat[k] <- sqrt(Krig.fs2hat(lambda.grid[k], info))
trA[k] <- Krig.ftrace(lambda.grid[k], D)
}
#
# clean up as amatrix with all these values.
gcv.grid <- cbind(lambda.grid, trA, V, V.one, V.model, shat,
lplike)
gcv.grid <- as.data.frame(gcv.grid)
names(gcv.grid) <- c("lambda", "trA", "GCV", "GCV.one", "GCV.model",
"shat", "-lnLike Prof")
# at this point the lazy data analyst could stop
# and just pick the minimum over the grid
# but we will now optimze the search producing refined optima
#
if (verbose)
print(info)
if (verbose)
print(gcv.grid)
#
# setup output matrix for refined values
lambda.est <- matrix(ncol = 5, nrow = 6, dimnames = list(c("GCV",
"GCV.model", "GCV.one", "RMSE", "pure error", "REML"),
c("lambda", "trA", "GCV", "shat", "REML")))
#
# now step through the many different ways to find lambda
# This is standard GCV w/o replicates, but with replicates this
# is leave-one-group out.
#
lambda.est[1, 1] <- Krig.find.gcvmin(info, lambda.grid, gcv.grid$GCV,
Krig.fgcv, tol = tol, verbose = verbose, give.warnings = give.warnings)
#
# If replicates then try GCV related to replicate group means
#
if (!is.na(shat.pure.error)) {
temp <- gcv.grid$GCV.model
lambda.est[2, 1] <- Krig.find.gcvmin(info, lambda.grid,
temp, Krig.fgcv.model, tol = tol, verbose = verbose,
give.warnings = give.warnings)
}
#
# Another version of GCV this is really leave an observation out
lambda.est[3, 1] <- Krig.find.gcvmin(info, lambda.grid, gcv.grid$GCV.one,
Krig.fgcv.one, tol = tol, verbose = verbose, give.warnings = give.warnings)
#
# REML
lambda.est[6, 1] <- Krig.find.REML(info, lambda.grid, gcv.grid$"-lnLike Prof",
Krig.flplike, tol = tol, verbose = verbose, give.warnings = give.warnings.REML)
if (verbose) {
cat(" mle estimate", lambda.est[6, 1], fill = TRUE)
}
#
# This is the most obscure way to choose lambda but
# Perry Haaland likes this one ....
# Find the values of lambda so that the
# usual estimate of sigma**2 is equal to either the
# value estimated from replicates and from an external value (rmse) if
# this has been passed.
# in some sense this a moment matching choice for lambda
lambda.rmse <- NA
lambda.pure.error <- NA
# matching external RMSE
if (!is.na(rmse)) {
if (all(gcv.grid$shat < rmse) | all(gcv.grid$shat > rmse)) {
guess <- NA
}
else {
guess <- max(gcv.grid$lambda[gcv.grid$shat < rmse])
}
if (verbose) {
print(rmse)
print(guess)
}
if (!is.na(guess)) {
lambda.rmse <- find.upcross(Krig.fs2hat, info, upcross.level = rmse^2,
guess = guess, tol = tol * rmse^2)
lambda.est[4, 1] <- lambda.rmse
}
else {
if (give.warnings) {
warning("Value of rmse is outside possible range")
}
}
}
#
# matching estimate of sigma from reps.
if (!is.na(shat.pure.error)) {
if (all(gcv.grid$shat < shat.pure.error) | all(gcv.grid$shat >
shat.pure.error)) {
guess <- NA
}
else {
guess <- max(gcv.grid$lambda[gcv.grid$shat < shat.pure.error])
}
if (!is.na(guess) & (guess != -Inf)) {
lambda.pure.error <- find.upcross(Krig.fs2hat, info,
upcross.level = shat.pure.error^2, guess = guess,
tol = tol * shat.pure.error^2)
lambda.est[5, 1] <- lambda.pure.error
}
else {
if (give.warnings) {
warning("Value of pure error estimate is outside possible range")
}
}
}
if (verbose) {
cat(" Done with refining", fill = TRUE)
}
#
# OK done with all six methods.
# fill in return matrix with all the right stuff
# fill in REML results
lam.ml <- lambda.est[6, 1]
lambda.est[6, 2] <- Krig.ftrace(lam.ml, D)
lambda.est[6, 3] <- Krig.fgcv(lam.ml, info)
lambda.est[6, 4] <- sqrt(Krig.fs2hat(lam.ml, info))
lambda.est[6, 5] <- Krig.flplike(lam.ml, info)
# fill in GCV results
for (k in 1:5) {
lam <- lambda.est[k, 1]
if (!is.na(lam)) {
lambda.est[k, 2] <- Krig.ftrace(lam, D)
if (k == 1 | k > 3) {
lambda.est[k, 3] <- Krig.fgcv(lam, info)
lambda.est[k, 5] <- Krig.flplike(lam, info)
}
if (k == 2) {
lambda.est[k, 3] <- Krig.fgcv.model(lam, info)
}
if (k == 3) {
lambda.est[k, 3] <- Krig.fgcv.one(lam, info)
}
lambda.est[k, 4] <- sqrt(Krig.fs2hat(lam, info))
}
if (verbose) {
cat("lam.est ", k, fill = TRUE)
cat(lambda.est[k, ], fill = TRUE)
}
}
# Note that the preferred estimate by default is
# leave-one (obs or group)-out GCV
list(gcv.grid = gcv.grid, lambda.est = lambda.est, lambda.best = lambda.est[1,
1])
}