R/predict.R

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

#' @export
predict.CCKriging = function(object, newdata, ...) {
  compute.prediction = function(object, newdata) {
    w = object$x  # mixed inputs used to build the model
    y = object$y
    d = ncol(w)
    n = nrow(w)
    cat.type = object$config$cat.type
    par = object$par  # estimated model parameters
    R = getCorrMatrix(object, ...)
    r0 = getCorrVector(object, x0 = newdata)
    # 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 = as.numeric(mu.hat + t(r0) %*% R.inv %*% (y - ones %*% mu.hat))

    sigma.sq = 1/n * t(y - ones %*% mu.hat) %*% 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
    sd = sqrt(sigma.sq) * sqrt(as.numeric(1 - t(c(1, r0)) %*% solve(tempmat) %*% c(1, r0)))
    ## sd = sigma.sq * (1 - t(r0) %*% R.inv %*% r0 + (1 - t(r0) %*% R.inv %*% ones) %*% solve(t(ones) %*% R.inv %*% ones) %*% t(1 - t(r0) %*% R.inv %*% ones))
    
    return(data.frame(mean = prediction, sd = sd))
  }

  res = sapply(seq_len(nrow(newdata)), function(i) compute.prediction(object, newdata[i, ]))
  res = data.frame(mean = unlist(res[1, ]), sd = unlist(res[2, ]))

  cc.pred = makeCCPred(mean = res$mean, sd = res$sd, x0 = newdata, model = object)
  return(cc.pred)
}