swh:1:snp:33a53053e50f7abe7d281cc0c803be827debf4a3
Tip revision: ce494d93a27106ca691c205c35c60843717717fa authored by Edzer Pebesma on 14 March 2022, 14:20:05 UTC
version 2.0-9
version 2.0-9
Tip revision: ce494d9
krigeTg.R
# $Id: krigeTg.q,v 1.4 2009-07-07 15:42:39 edzer Exp $
phiInv <- function (x, lambda)
if (lambda==0) log(x) else (x^lambda-1)/lambda
phi <- function(x, lambda)
if (lambda==0) exp(x) else (x*lambda+1)^(1/lambda)
phiPrime <- function (x, lambda)
if (lambda==0) exp(x) else (x*lambda+1)^(1/lambda-1)
phiDouble <- function (x, lambda)
if (lambda==0) exp(x) else
lambda * (1/lambda - 1) * (lambda * x + 1)^(1/lambda-2)
krigeTg <- function(formula, locations, newdata, model = NULL, ...,
nmax = Inf, nmin = 0, maxdist = Inf, block = numeric(0),
nsim = 0, na.action = na.pass, debug.level = 1,
lambda = 1.0)
{
m = model.frame(terms(formula), as.data.frame(locations))
Y = model.extract(m, "response")
if (length(Y) == 0)
stop("no response variable present in formula")
Terms = attr(m, "terms")
X = model.matrix(Terms, m)
has.intercept = attr(Terms, "intercept")
if (ncol(X) > 1)
stop("only formula with intercept allowed, e.g. y ~ 1")
locations$value = phiInv(Y, lambda)
locations$value1 = rep(1, length(locations$value))
OK = krige(value ~ 1, locations, newdata, model,
nmax = nmax, nmin = nmin,
maxdist = maxdist, block = block, nsim = nsim,
na.action = na.action, debug.level = debug.level, ...)
if (nsim > 0) {
OK@data = as.data.frame(phi(OK@data,lambda))
return(OK)
}
# else:
# estimate mu:
g = gstat(formula = value ~ 1, # locations = locations,
data = locations, model = model, nmax = nmax, nmin = nmin,
maxdist = maxdist, ...)
mu = predict(g, newdata = newdata, block = block, nsim = nsim,
na.action = na.action, debug.level = debug.level, BLUE = TRUE)
OK$muhat = mu$var1.pred
SK = krige(value1 ~ 1, locations, newdata,
model = model, beta = 0.0, nmax = nmax, nmin = nmin,
maxdist = maxdist, block = block, nsim = nsim,
na.action = na.action, debug.level = debug.level, ...)
# find m:
OK$m = (OK$var1.var - SK$var1.var)/(1 - SK$var1.pred)
# copy SK output:
OK$var1SK.pred = SK$var1.pred
OK$var1SK.var = SK$var1.var
# compute transGaussian kriging estimate & variance:
OK$var1TG.pred = phi(OK$var1.pred, lambda) +
phiDouble(mu$var1.pred, lambda) * (OK$var1.var/2 - OK$m)
OK$var1TG.var = phiPrime(mu$var1.pred, lambda)^2 * OK$var1.var
OK
}