https://github.com/cran/fields
Tip revision: 328f151e07864e976578f9134025cbfdb72eb2f3 authored by Doug Nychka on 19 July 2013, 23:45:29 UTC
version 6.8
version 6.8
Tip revision: 328f151
sim.Krig.R
# 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
"sim.Krig" <- function(object, xp, M = 1,
verbose = FALSE) {
sigma2 <- object$best.model[2]
rho <- object$best.model[3]
#
# check for unique rows of xp
if (any(duplicated(xp))) {
stop(" predictions locations should be unique")
}
#
# set up various sizes of arrays
m <- nrow(xp)
n <- nrow(object$xM)
N <- length(object$y)
if (verbose) {
cat(" m,n,N", m, n, N, fill = TRUE)
}
#transform the new points
xc <- object$transform$x.center
xs <- object$transform$x.scale
xpM <- scale(xp, xc, xs)
# complete set of points for prediction.
# check for replicates and adjust
x <- rbind(object$xM, xpM)
if (verbose) {
cat("full x ", fill = TRUE)
print(x)
}
#
# find indices of all rows of xp that correspond to rows of
# xM and then collapse x to unique rows.
rep.x.info <- fields.duplicated.matrix(x)
x <- as.matrix(x[!duplicated(rep.x.info), ])
if (verbose) {
cat("full x without duplicates ", fill = TRUE)
print(x)
}
N.full <- nrow(x)
if (verbose) {
cat("N.full", N.full, fill = TRUE)
}
# these give locations in x matrix to reconstruct xp matrix
xp.ind <- rep.x.info[(1:m) + n]
if (verbose) {
print(N.full)
print(x)
}
if (verbose) {
cat("reconstruction of xp from collapsed locations",
fill = TRUE)
print(x[xp.ind, ])
}
#
# Sigma is full covariance at the data locations and at prediction points.
#
Sigma <- rho * do.call(object$cov.function.name, c(object$args,
list(x1 = x, x2 = x)))
#
# square root of Sigma for simulating field
# Cholesky is fast but not very stable.
#
# the following code line is similar to chol(Sigma)-> Scol
# but adds possible additional arguments controlling the Cholesky
# from the Krig object.
#
Schol <- do.call("chol", c(list(x = Sigma), object$chol.args))
#
# output matrix to hold results
N.full <- nrow(x)
out <- matrix(NA, ncol = m, nrow = M)
#
# find conditional mean field from initial fit
# don't multiply by sd or add mean if this is
# a correlation model fit.
# (these are added at the predict step).
h.hat <- predict(object, xp)
# marginal standard deviation of field.
temp.sd <- 1
#
#
# this is not 1 if Krig object is a corelation model.
if (object$correlation.model) {
if (!is.na(object$sd.obj[1])) {
temp.sd <- predict(object$sd.obj, x)
}
}
#
# Define W2i for simulating errors.
#
W2i <- Krig.make.Wi(object)$W2i
for (k in 1:M) {
# simulate full field
h <- t(Schol) %*% rnorm(N.full)
# value of simulated field at observations
#
# NOTE: fixed part of model (null space) need not be simulated
# because the estimator is unbiased for this part.
# the variability is still captured because the fixed part
# is still estimated as part of the predict step below
h.data <- h[1:n]
# expand the values according to the replicate pattern
h.data <- h.data[object$rep.info]
# create synthetic data
y.synthetic <- h.data + sqrt(sigma2) * W2i %d*% rnorm(N)
# predict at xp using these data
# and subtract from 'true' value
# note that true values of field have to be expanded in the
# case of common locations between xM and xp.
h.true <- (h[xp.ind])
temp.error <- predict(object, xp, y = y.synthetic, eval.correlation.model = FALSE) -
h.true
# add the error to the actual estimate (conditional mean)
# and adjust by marginal standard deviation
out[k, ] <- h.hat + temp.error * temp.sd
}
out
}