swh:1:snp:33a53053e50f7abe7d281cc0c803be827debf4a3
Tip revision: ae42623708becf61734bc74f16126dad36e203c1 authored by Edzer J. Pebesma on 02 February 2004, 16:39:05 UTC
version 0.9-10
version 0.9-10
Tip revision: ae42623
fit.lmc.R
"fit.lmc" <-
function (v, g, model, fit.ranges = FALSE, fit.lmc = !fit.ranges,
...)
{
posdef = function(X) {
q = eigen(X)
d = q$values
d[d < 0] = 0
q$vectors %*% diag(d) %*% t(q$vectors)
}
if (!inherits(v, "variogram"))
stop("v should be of class variogram")
if (!inherits(g, "gstat"))
stop("g should be of class gstat")
if (!missing(model)) {
if (!inherits(model, "variogram.model"))
stop("model should be of class variogram.model")
}
n = names(g$data)
for (i in 1:length(n)) {
for (j in i:length(n)) {
name = ifelse(i == j, n[i], cross.name(n[i], n[j]))
x = v[v$id == name, ]
if (nrow(x) == 0)
stop(paste("variogram", name, "not present"))
m = g$model[[name]]
if (!missing(model))
m = model
g$model[[name]] = fit.variogram(x, m, fit.ranges = fit.ranges,
...)
}
}
if (fit.lmc) {
m = g$model[[n[1]]]
for (k in 1:nrow(m)) {
psill = matrix(NA, nrow = length(n), ncol = length(n))
for (i in 1:length(n)) {
for (j in i:length(n)) {
name = ifelse(i == j, n[i], cross.name(n[i], n[j]))
psill[i, j] = psill[j, i] = g$model[[name]][k,
"psill"]
}
}
psill = posdef(psill)
for (i in 1:length(n)) {
for (j in i:length(n)) {
name = ifelse(i == j, n[i], cross.name(n[i], n[j]))
g$model[[name]][k, "psill"] = psill[i, j]
}
}
}
}
g
}