https://github.com/cran/fields
Tip revision: 6769ffc81115fbf0bf7d9c566cf7ac81be0049dc authored by Doug Nychka on 25 July 2005, 00:00:00 UTC
version 3.04
version 3.04
Tip revision: 6769ffc
predict.se.Krig.r
"predict.se.Krig" <-
function (object, x=NULL, cov = FALSE, verbose=FALSE,...)
{
#
# name of covariance function
call.name<- object$cov.function.name
#
# default is to predict at data x's
if (is.null(x)) {
x <- object$x
}
x <- as.matrix(x)
if (verbose)
print(x)
xraw<- x
# transformations of x values used in Krig
xc <- object$transform$x.center
xs <- object$transform$x.scale
x <- scale(x, xc, xs)
#
# scaled unique observation locations.
xM <- object$xM
# find marginal variance before transforming x.
temp.sd <- 1
if (!is.na( object$sd.obj[1])) {
temp.sd <- c(predict(object$sd.obj, xraw))
}
# NOTE knots are already scaled in Krig object and are used
# Default is to use parameters in best.model
lambda <- object$best.model[1]
rho<- object$best.model[3]
sigma2<- object$best.model[2]
nx <- nrow(xM)
wght.vec <- t(Krig.Amatrix(object, xraw, lambda))
if( verbose) {
cat("wght.vector", fill=TRUE)
print( wght.vec)}
#var( f0 - yhat)= var( f0) - cov( f0,yhat) - cov( yhat, f0) + cov( yhat)
# = temp0 - temp1 - t( temp1) + temp2
if ( !cov){
# find diagonal elements of covariance matrix
# covariance of data
Cov.y <- rho * do.call(
call.name,
c(object$args, list(x1 = xM, x2 = xM))) +
sigma2 * diag(1/object$weightsM)
# now find the three terms.
# note the use of an element by element multiply to only get the
# diagonal elements of the full
# prediction covariance matrix.
#
temp1 <- rho * colSums(
wght.vec* do.call(
call.name, c(object$args,list(x1 = xM, x2 = x) ) )
)
temp2 <- colSums(wght.vec * (Cov.y %*% wght.vec) )
#
# find marginal variances -- trival in the stationary case!
#
temp0 <- rho * do.call(
call.name,
c(object$args, list(x1 = x, marginal=TRUE)) )
#
temp <- temp0 - 2 * temp1 + temp2
#
return(sqrt(temp * temp.sd^2 ))
}
else{
# find full covariance matrix
# covariance of data
Cov.y <- rho *
do.call(
call.name, c(object$args, list(x1 = xM, x2 = xM))) +
sigma2 * diag(1/object$weightsM) # now find the three terms.
#
temp1 <- rho *
t(wght.vec)%*% do.call(
call.name, c(object$args,list(x1 = xM, x2=x)))
#
temp2 <- t(wght.vec) %*% Cov.y %*% wght.vec
#
temp0 <- rho * do.call(
call.name, c(object$args, list(x1 = x, x2 = x)))
#
temp <- temp0 - t(temp1) - temp1 + temp2
temp<- t(t(temp)* temp.sd )* temp.sd
#
return(temp)
}
}