distMatToPar = function(dist.mat, par) {
par = c(par, dist.mat[lower.tri(dist.mat)])
return(par)
}
## parToQualDistMat: Generate (lower triangular) dist matrix for all combinations
## of the qualitative variables. First d parameters of vector par are ignored
## because they are length scale parameters of the kriging model
parToQualDistMat = function(par, d) {
par = par[-(1:d)]
n.par = length(par)
dim.mat = sqrt(2*n.par + 0.25) + 0.5
if (dim.mat %% 1 != 0) {
stop("Cannot generate quadratic lower triangular matrix from vector par.")
}
distMat = matrix(0, nrow = dim.mat, ncol = dim.mat)
distMat[lower.tri(distMat)] = par
return(distMat)
}
## get the distances of two qualitative vectors
getQualDists = function(x1, x2, dist.mat) {
qual.dists = sapply(seq_along(x1), function(i) {
if (x1[i] == x2[i]) {
return(0)
}
x1.lvl = which(levels(x1) == x1[i])
x2.lvl = which(levels(x1) == x2[i])
if (x1.lvl > x2.lvl) {
return(dist.mat[x1.lvl, x2.lvl])
} else {
return(dist.mat[x2.lvl, x1.lvl])
}
})
return(qual.dists)
}
##### predict
predict.cdKriging = function(object, newdata, ...) {
w = object$x # mixed inputs used to build the model
y = object$y
d = ncol(w)
n = nrow(w)
par = object$par # estimated model parameters
R = getCorrMatrix(x = w, par = par, ...)
qual.dist = parToQualDistMat(par, d = d)
r0 = getCorrVector(x = w, x0 = newdata, par = par, qual.dist = qual.dist, ...)
ones = rep(1, n)
R.inv = solve(R)
mu.hat = solve((t(ones) %*% R.inv %*% ones)) %*% t(ones) %*% R.inv %*% y
prediction = mu.hat + t(r0) %*% R.inv %*% (y - ones %*% mu.hat)
tempmat = matrix(nrow = n + 1, ncol = n + 1)
tempmat[1, ] = c(0, ones)
tempmat[, 1] = c(0, ones)
tempmat[2:(n+1), 2:(n+1)] = R
s.square = 1 - t(c(1, r0)) %*% solve(tempmat) %*% c(1, r0) # FIXME: * sigma^2
return(list(prediction = prediction, s.square = s.square))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.