attic/CatDist.R

In dominikkirchhoff/CCKriging: Kriging for Mixed Continuous and Categorical Input Variables

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))
}