https://github.com/cran/fields
Tip revision: 92636eaee56714c89c7b7e7cdc52e0c77aeb79b3 authored by Doug Nychka on 07 February 2005, 00:00:00 UTC
version 2.3
version 2.3
Tip revision: 92636ea
predict.se.Krig.r
"predict.se.Krig" <-
function(object, x, cov.function=NULL, rho, sigma2, weights = NULL, cov
= FALSE, stationary
= TRUE, fixed.mean = TRUE, ...)
{
out<-object # hack for S3
# NOTE: to follow along
# if lambda is an estimate of the random field at a point
# and lambda_hat is an estimate based on linear combination of the data
#
# MSE( lambda - lambda_hat) = VAR( lambda) - 2* COV( lambda, lambda_hat)
# + VAR( lambda_hat)
#
# Perhaps the most confusing aspect of this code is that the MSE
# are found simultanouesly for all the points in x.
# This is to avoid an explicit for loop
#
#
if(missing(x)) x <- out$x
x <- as.matrix(x)
xraw <- x
#### transform x to the scale used to fit model
if(is.null(out$transform)) {
xc <- rep(0, ncol(x))
xs <- rep(1, ncol(x))
}
else {
xc <- out$transform$x.center
xs <- out$transform$x.scale
}
x <- scale(x, xc, xs)
xM <- scale(out$xM, xc, xs)
##
########## if covariance function and parameters are missing
########## extract them from the Krig object
##
## figure out what covariance function to use.
if(!is.null(cov.function)) {
if(is.function(cov.function))
cov.function <- as.character(substitute(cov.function))
}
else{ cov.function<- out$call.name}
if(missing(sigma2)) {
sigma2 <- out$sigma2
}
#
# OK now fix the right value for sigma
#
if(missing(rho)) {
rho <- out$rho
}
#
# set switch if the covariance is stationary
if(!is.null(formals(get(out$call.name))$marginal)) {
stationary <- FALSE
}
#
lambda <- sigma2/rho
if(lambda != 0) {
if((lambda - out$lambda)/lambda > 1e-08) {
warning("lambda value used is different from the one in the Krig object"
)
cat(lambda, out$lambda, fill = TRUE)
}
}
nx <- nrow(xM)
#
# figure out if this is a correlation model and make adustments
#
temp.sd <- 1
var.mean <- 0
if(out$correlation.model) {
# predict standard deviations at xraw values
temp.sd <- c(predict(out$sd.obj, xraw))
if(!fixed.mean) {
# for the very end ... the variance in estimating the mean function.
var.mean <- predict.se(out$mean.obj, xraw)^2
}
}
#
# here is where the real fun begins
#
# wght.vec are the linear combinations of the data ( yM the averages
# if there are replicates) that give the
# correpsonding estimates of the function at the points x
wght.vec <- t(Krig.Amatrix(out, xraw, lambda))
#
# modify weight.vec to give the weighted estimated
#
# Cy is the observed covariance matrix of the data vector
# xM should be in transformed scale and are the unique x values from
# the collapsed data means
#
Cy <- rho * do.call(out$call.name, c(out$args, list(x1 = xM,
x2 = xM))) + sigma2 * diag(1/out$weightsM)
#
# see if the C argument is given, if so utilize it.
C.arg.missing <- is.null(formals(get(out$call.name))$C)
if(cov) {
# this block find the whole covariance matrix
if(C.arg.missing) {
temp <- rho * do.call(out$call.name, c(out$
args, list(x1 = x, x2 = xM))) %*%
wght.vec
}
else {
temp <- rho * do.call(out$call.name, c(out$
args, list(x1 = x, x2 = xM, C =
wght.vec)))
}
# end of if else C.arg.missing stmt
temp <- rho * do.call(out$call.name, c(out$args, list(
x1 = x, x2 = x))) - temp - t(temp)
temp <- temp + t(wght.vec) %*% Cy %*% wght.vec
return(temp.sd * t(t(temp * temp.sd)))
}
#
# this is the second term in the MSE
temp1 <- rho * c(t(wght.vec * do.call(out$call.name, c(out$
args, list(x1 = xM, x2 = x)))) %*% rep(1, nx))
#
# this is the third term in the MSE
temp2 <- c(t(wght.vec * (Cy %*% wght.vec)) %*% rep(1, nx))
#
# VAR of the first term is handled differently depending on
# if this is the MSE at a point or a weighted combination
#
if(stationary) {
x0 <- matrix(0, ncol = ncol(x), nrow = 1)
temp <- rho * do.call(out$call.name, c(out$args, list(
x1 = x0, x2 = x0))) - 2 * temp1 + temp2
return(sqrt(temp * temp.sd^2 + var.mean))
}
#
#
# if covariance is not stationary then loop through each point to get
# the variance of field at that point.
#
#
if(is.null(formals(get(out$call.name))$marginal)) {
temp <- rep(0, nrow(x))
for(k in 1:nrow(x)) {
x0 <- matrix(x[k, ], nrow = 1)
temp[k] <- rho * do.call(out$call.name, c(
out$args, list(x1 = x0, x2 = x0))) -
2 * temp1[k] + temp2[k]
}
}
else {
#
## marginal variances available by single call
#
temp <- rho * do.call(out$call.name, c(out$args, list(
x, marginal = TRUE))) - 2 * temp1 + temp2
}
return(sqrt((temp.sd^2) * temp + var.mean))
}