Nothing
R/buildKrigingForrester.R
In SPOT: Sequential Parameter Optimization Toolbox
Defines functions
print.kriging
predictKrigingReinterpolation
predict.kriging
krigingLikelihood
normalizeMatrix
buildKriging
Documented in
buildKriging
krigingLikelihood
normalizeMatrix
predict.kriging
predictKrigingReinterpolation
print.kriging
#' @title Build Kriging Model
#'
#' @description This function builds a Kriging model based on code by Forrester et al..
#' By default exponents (p) are fixed at a value of two, and a nugget (or regularization constant) is used.
#' To correct the uncertainty estimates in case of nugget, re-interpolation is also by default turned on.
#'
#' @details The model uses a Gaussian kernel: \code{k(x,z)=exp(-sum(theta_i * |x_i-z_i|^p_i))}. By default, \code{p_i = 2}.
#' Note that if dimension \code{x_i} is a factor variable (see parameter \code{types}), Hamming distance will be used
#' instead of \code{|x_i-z_i|}.
#'
#' @param x design matrix (sample locations)
#' @param y vector of observations at \code{x}
#' @param control (list), with the options for the model building procedure. Note: This can also be changed after the model has been built, by manipulating the respective \code{object$target} value.
#' \describe{
#' \item{\code{types}}{a character vector giving the data type of each variable.
#' All but "factor" will be handled as numeric, "factor" (categorical) variables
#' will be subject to the hamming distance.}
#' \item{\code{thetaLower}}{lower boundary for theta, default is \code{1e-4}}
#' \item{\code{thetaUpper}}{upper boundary for theta, default is \code{1e2}}
#' \item{\code{algTheta}}{algorithm used to find theta, default is \code{optimDE}.}
#' \item{\code{budgetAlgTheta}}{budget for the above mentioned algorithm,
#' default is \code{200}. The value will be multiplied with the length of
#' the model parameter vector to be optimized.}
#' \item{\code{optimizeP}}{boolean that specifies whether the exponents (\code{p}) should be optimized. Else they will be set to two. Default is \code{FALSE}.}
#' \item{\code{useLambda}}{whether or not to use the regularization constant lambda
#' (nugget effect). Default is \code{TRUE}.}
#' \item{\code{lambdaLower}}{lower boundary for log10{lambda}, default is \code{-6}}
#' \item{\code{lambdaUpper}}{upper boundary for log10{lambda}, default is \code{0}}
#' \item{\code{startTheta}}{optional start value for theta optimization, default is \code{NULL}}
#' \item{\code{reinterpolate}}{whether (\code{TRUE},default) or not (\code{FALSE}) reinterpolation
#' should be performed.}
#' \item{\code{target}}{values of the prediction, a vector of strings.
#' Each string specifies a value to be predicted, e.g., "y" for mean, "s" for standard deviation, "ei" for expected improvement.
#' See also \code{\link{predict.kriging}}.}
#' }
#'
#' @return an object of class \code{kriging}. Basically a list, with the options
#' and found parameters for the model which has to be passed to the predictor function:\cr
#' \code{x} sample locations (scaled to values between 0 and 1)\cr
#' \code{y} observations at sample locations (see parameters)\cr
#' \code{thetaLower} lower boundary for theta (see parameters)\cr
#' \code{thetaUpper} upper boundary for theta (see parameters)\cr
#' \code{algTheta} algorithm to find theta (see parameters)\cr
#' \code{budgetAlgTheta} budget for the above mentioned algorithm (see parameters)\cr
#' \code{optimizeP} boolean that specifies whether the exponents (\code{p}) were optimized (see parameters)\cr
#' \code{normalizeymin} minimum in normalized space\cr
#' \code{normalizeymax} maximum in normalized space\cr
#' \code{normalizexmin} minimum in input space\cr
#' \code{normalizexmax} maximum in input space\cr
#' \code{dmodeltheta} vector of activity parameters\cr
#' \code{Theta} log_10 vector of activity parameters (i.e. \code{log10(dmodeltheta)})\cr
#' \code{dmodellambda} regularization constant (nugget) \cr
#' \code{Lambda} log_10 of regularization constant (nugget) (i.e. \code{log10(dmodellambda)})\cr
#' \code{yonemu} \code{Ay-ones*mu} \cr
#' \code{ssq} sigma square\cr
#' \code{mu} mean mu\cr
#' \code{Psi} matrix large Psi\cr
#' \code{Psinv} inverse of Psi\cr
#' \code{nevals} number of Likelihood evaluations during MLE
#'
#' @importFrom stats predict
#' @importFrom stats dist
#' @importFrom stats var
#'
#' @export
#' @seealso \code{\link{predict.kriging}}
#' @references Forrester, Alexander I.J.; Sobester, Andras; Keane, Andy J. (2008). Engineering Design via Surrogate Modelling - A Practical Guide. John Wiley & Sons.
#'
#' @examples
#' ## Create design points
#' set.seed(1)
#' x <- cbind(runif(20)*15-5,runif(20)*15)
#' y <- funBranin(x)
#' ## Create model with default settings
#' fit <- buildKriging(x,y,control = list(algTheta=optimLHD))
#' ## Print model parameters
#' print(fit)
#' ## Predict at new location
#' predict(fit,cbind(1,2))
#' ## True value at location
#' funBranin(matrix(c(1,2), 1))
#' ##
#' ## Next Example: Handling factor variables
#' \donttest{
#' ## create a test function:
#' braninFunctionFactor <- function (x) {
#' y <- (x[2] - 5.1/(4 * pi^2) * (x[1] ^2) + 5/pi * x[1] - 6)^2 +
#' 10 * (1 - 1/(8 * pi)) * cos(x[1] ) + 10
#' if(x[3]==1)
#' y <- y +1
#' else if(x[3]==2)
#' y <- y -1
#' y
#' }
#' ## create training data
#' set.seed(1)
#' x <- cbind(runif(50)*15-5,runif(50)*15,sample(1:3,50,replace=TRUE))
#' y <- as.matrix(apply(x,1,braninFunctionFactor))
#' ## fit the model (default: assume all variables are numeric)
#' fitDefault <- buildKriging(x,y,control = list(algTheta=optimDE))
#' ## fit the model (give information about the factor variable)
#' fitFactor <- buildKriging(x,y,control =
#' list(algTheta=optimDE,types=c("numeric","numeric","factor")))
#' ## create test data
#' xtest <- cbind(runif(200)*15-5,runif(200)*15,sample(1:3,200,replace=TRUE))
#' ytest <- as.matrix(apply(xtest,1,braninFunctionFactor))
#' ## Predict test data with both models, and compute error
#' ypredDef <- predict(fitDefault,xtest)$y
#' ypredFact <- predict(fitFactor,xtest)$y
#' mean((ypredDef-ytest)^2)
#' mean((ypredFact-ytest)^2)
#' }
buildKriging <- function(x, y, control = list()) {
## (K-1) Set Parameters
x <-
data.matrix(x)
# data.matrix is better than as.matrix,
# because it always converts to numeric, not character! (in case of mixed data types.)
k <- ncol(x) # dimension
n <- nrow(x) # number of observations
con <- list(
thetaLower = 1e-4,
thetaUpper = 1e2,
types = rep("numeric", k),
algTheta = optimDE,
budgetAlgTheta = 200,
optimizeP = FALSE,
useLambda = TRUE,
lambdaLower = -6,
lambdaUpper = 0,
startTheta = NULL,
reinterpolate = TRUE,
target = "y"
)
con[names(control)] <- control
control <- con
fit <- control
fit$x <- x
fit$y <- y
LowerTheta <- rep(1, k) * log10(fit$thetaLower)
UpperTheta <- rep(1, k) * log10(fit$thetaUpper)
## (K-2) Normalize input data
ymin <- 0
ymax <- 1
fit$normalizeymin <- ymin
fit$normalizeymax <- ymax
# normalization in each column:
res <- normalizeMatrix(fit$x, ymin, ymax)
fit$scaledx <- res$y
fit$normalizexmin <- res$xmin
fit$normalizexmax <- res$xmax
## (K-3) Prepare distance/correlation matrix
A <- matrix(0, k, n * n)
for (i in 1:k) {
if (control$types[i] != "factor") {
A[i, ] <-
as.numeric(as.matrix(dist(fit$scaledx[, i]))) #euclidean distance
} else {
tmp <-
outer(fit$scaledx[, i], fit$scaledx[, i], '!=') #hamming distance
class(tmp) <- "numeric"
A[i, ] <- tmp
}
}
## (K-4) Prepare starting points, search bounds and penalty value for MLE optimization
## 4.1 theta
if (is.null(fit$startTheta))
x1 <- rep(n / (100 * k), k) # start point for theta
else
x1 <- fit$startTheta
## 4.2 p
if (fit$optimizeP) {
LowerTheta <- c(LowerTheta, rep(1, k) * 0.01)
UpperTheta <- c(UpperTheta, rep(1, k) * 2)
x3 <- rep(1, k) * 1.9 #start values for p
x0 <- c(x1, x3)
} else{
# p is fixed to 2 and the array A is completely precalculated
A <- A ^ 2
x0 <- c(x1)
}
## 4.3 lambda
if (fit$useLambda) {
# start value for lambda:
x2 <- (fit$lambdaUpper + fit$lambdaLower) / 2
x0 <- c(x0, x2)
#append regression constant lambda (nugget)
LowerTheta <- c(LowerTheta, fit$lambdaLower)
UpperTheta <- c(UpperTheta, fit$lambdaUpper)
}
#force x0 into bounds
x0 <- pmin(x0, UpperTheta)
x0 <- pmax(x0, LowerTheta)
x0 <- matrix(x0, 1) #matrix with one row
opts <- list(funEvals = fit$budgetAlgTheta * ncol(x0))
## 4.4 penalty
#determine a good penalty value (based on number of samples and variance of y)
penval <- n * log(var(y)) + 1e4
## (K-5) MLE objective function
# Wrapper for optimizing theta based on krigingLikelihood:
fitFun <-
function (x, fX, fy, optimizeP, useLambda, penval) {
#todo vectorize, at least for cma_es with active vectorize?
as.numeric(krigingLikelihood(x, fX, fy, optimizeP, useLambda, penval)$NegLnLike)
}
## (K-7) MLE optimization
res <- fit$algTheta(
x = x0,
fun =
function(x, fX, fy, optimizeP, useLambda, penval) {
if (!is.matrix(x)) {
fitFun(x, fX, fy, optimizeP, useLambda, penval)
}
else{
apply(x, 1, fitFun, fX, fy, optimizeP, useLambda, penval)
}
},
lower = LowerTheta,
upper = UpperTheta,
control = opts,
fX = A,
fy = fit$y,
optimizeP = fit$optimizeP,
useLambda = fit$useLambda,
penval = penval
)
if (is.null(res$xbest))
res$xbest <- x0
## (K-8): compile results from MLE optimization (to fit object)
Params <- res$xbest
nevals <- as.numeric(res$count[[1]])
fit$Theta <- Params[1:k]
fit$dmodeltheta <- 10 ^ Params[1:k]
if (fit$optimizeP) {
fit$P <- Params[(k + 1):(2 * k)]
}
if (fit$useLambda) {
fit$Lambda <- Params[length(Params)]
fit$dmodellambda <- 10 ^ Params[length(Params)]
} else{
fit$Lambda <- -Inf
fit$dmodellambda <- 0
}
## (K-9) Evaluate with optimized parameters
res <-
krigingLikelihood(c(fit$Theta, fit$P, fit$Lambda),
A,
fit$y,
fit$optimizeP,
fit$useLambda)
if (is.na(res$Psinv[1]))
stop(
"buildKriging failed to produce a valid correlation matrix. Consider activating the nugget/regularization constant via useLambda=TRUE in the control list."
)
## (K-10) Add results from MLE evaluation to fit object
fit$yonemu <- res$yonemu
fit$ssq <- as.numeric(res$ssq)
fit$mu <- res$mu
fit$Psi <- res$Psi
fit$Psinv <- res$Psinv
fit$nevals <- nevals
fit$like <- res$NegLnLike
fit$returnCrossCor <- FALSE
## (K-11) Calculate observed minimum
xlist <- split(x, 1:nrow(x))
uniquex <- unique(xlist)
ymean <- NULL
for (xi in uniquex) {
ind <- xlist %in% list(xi)
ymean <- c(ymean, mean(y[ind]))
}
fit$min <- min(ymean)
class(fit) <- "kriging"
return(fit)
}
#' Normalize design 2
#'
#' Normalize design with given maximum and minimum in input space.
#' Supportive function for Kriging model, not to be used directly.
#'
#' @param x design matrix in input space (n rows for each point, k columns for each parameter)
#' @param ymin minimum vector of normalized space
#' @param ymax maximum vector of normalized space
#' @param xmin minimum vector of input space
#' @param xmax maximum vector of input space
#'
#' @return normalized design matrix
#' @seealso \code{\link{buildKriging}}
#' @export
normalizeMatrix2 <- function (x, ymin, ymax, xmin, xmax) {
rangex <- xmax - xmin
rangey <- ymax - ymin
s <- dim(x)[1]
rangey * (x - matrix(rep(xmin, s), nrow = s, byrow = TRUE)) / matrix(rep(rangex, s), nrow =
s, byrow = TRUE) + ymin
}
#' @title Normalize design matrix
#'
#' @description Normalize design by using minimum and maximum of the design values
#' for input space. Each column has entries in the range from \code{ymin} to \code{ymax}.
#'
#' @param x design matrix in input space
#' @param ymin minimum vector of normalized space
#' @param ymax maximum vector of normalized space
#' @param MARGIN a vector giving the subscripts which the function will be applied over.
#' 1 indicates rows, 2 indicates columns. Default: \code{2}.
#'
#' @return \code{list} with the following entries:
#' \describe{
#' \item{\code{y}}{normalized design matrix in the range [ymin, ymax]}
#' \item{\code{xmin}}{min in each column}
#' \item{\code{xmax}}{max in each column}
#' }
#'
#' @seealso \code{\link{buildKriging}}
#' @examples
#' set.seed(1)
#' x <- matrix(c(rep(1,3), rep(2,3),rep(3,3), rep(4,3)),3,4)
#' ## columnwise:
#' normalizeMatrix(x, ymin=0, ymax=1)
#' ## rowwise
#' normalizeMatrix(x, ymin=0, ymax=1, MARGIN=1)
#' # rows with identical values are mapped to the mean:
#' x <- matrix(rep(0,4),2,2)
#' normalizeMatrix(x, ymin=0, ymax=1)
#'
#' @export
normalizeMatrix <- function(x, ymin, ymax, MARGIN = 2) {
if (MARGIN == 1)
x <- t(x)
# Return the maximum from each row:
xmax <- apply(x, 2, max)
# Return the minimum from each row:
xmin <- apply(x, 2, min)
s <- dim(x)[1]
rangex <- xmax - xmin
rangey <- ymax - ymin
xmin[rangex == 0] <- xmin[rangex == 0] - 0.5
xmax[rangex == 0] <- xmax[rangex == 0] + 0.5
rangex[rangex == 0] <- 1
y <-
rangey * (x - matrix(rep(xmin, s), nrow = s, byrow = TRUE)) / matrix(rep(rangex, s), nrow =
s, byrow = TRUE) + ymin
if (MARGIN == 1)
y <- t(y)
list(y = y, xmin = xmin, xmax = xmax)
}
#' Calculate negative log-likelihood
#'
#' Used to determine theta/lambda values for the Kriging model in \code{\link{buildKriging}}.
#' Supportive function for Kriging model, not to be used directly.
#'
#' @param x vector, containing parameters log10(theta), log10(lambda) and p.
#' @param AX 3 dimensional array, constructed by buildKriging from the sample locations
#' @param Ay vector of observations at sample locations
#' @param optimizeP boolean, whether to optimize parameter p (exponents) or fix at two.
#' @param useLambda boolean, whether to use nugget
#' @param penval a penalty value which affects the value returned for invalid correlation matrices / configurations
#'
#' @return list with elements\cr
#' \code{NegLnLike} concentrated log-likelihood *-1 for minimising \cr
#' \code{Psi} correlation matrix\cr
#' \code{Psinv} inverse of correlation matrix (to save computation time during prediction)\cr
#' \code{mu} \cr
#' \code{ssq}
#' @seealso \code{\link{buildKriging}}
#' @export
#' @keywords internal
krigingLikelihood <-
function(x,
AX,
Ay,
optimizeP = FALSE,
useLambda = TRUE,
penval = 1e8) {
## note: this penalty value should not be a hard constant.
## the scale of the likelihood (n*log(SigmaSqr) + LnDetPsi)
## at least depends on log(var) and the number of samples
## hence, large number of samples may lead to cases where the
## penality is lower than the likelihood of most valid parameterizations
## suggested penalty:
#penval <- n*log(var(y)) + 1e4
## currently, this better penalty is set in the buildKriging function,
## when calling krigingLikelihood
## (L-1) Starting Points
nx <- nrow(AX)
theta <- 10 ^ x[1:nx]
if (optimizeP) {
AX <- abs(AX) ^ (x[(nx + 1):(2 * nx)])
}
lambda <- 0
if (useLambda) {
lambda <- 10 ^ x[length(x)]
}
if (any(theta == 0) ||
any(is.infinite(c(theta, lambda)))) {
#for instance caused by bound violation
return(list(
NegLnLike = 1e4,
Psi = NA,
Psinv = NA,
mu = NA,
ssq = NA
))
}
n <- dim(Ay)[1]
## (L-2) Correlation Matrix Psi
Psi <- exp(-matrix(colSums(theta * AX), n, n))
if (useLambda) {
Psi <- Psi + diag(lambda, n)
}
if (any(is.infinite(Psi))) {
# this is required especially if distance matrices are forced to be CNSD/NSD and hence have zero distances
penalty <- penval
return(list(
NegLnLike = penalty,
Psi = NA,
Psinv = NA,
mu = NA,
ssq = NA
)) #todo: add something like smallest eigenvalue to penalty?
}
## Check whether Psi is ill-conditioned
kap <- rcond(Psi)
if (is.na(kap)) {
kap <- 0
}
if (kap < 1e-10) {
#warning("Correlation matrix is ill-conditioned (During Maximum Likelihood Estimation in buildKriging). Returning penalty.")
return(
list(
NegLnLike = penval - log10(kap),
Psi = NA,
Psinv = NA,
mu = NA,
ssq = NA,
a = NA,
U = NA,
isIndefinite = TRUE
)
)
}
## (L-3) cholesky decomposition
cholPsi <- try(chol(Psi), TRUE)
## give penalty if fail
if (inherits(cholPsi, "try-error")) {
#warning("Correlation matrix is not positive semi-definite (During Maximum Likelihood Estimation in buildKriging). Returning penalty.")
penalty <-
penval - log10(min(eigen(
Psi, symmetric = TRUE, only.values = TRUE
)$values))
return(list(
NegLnLike = penalty,
Psi = NA,
Psinv = NA,
mu = NA,
ssq = NA
))
}
## (L-4) Determininant
## calculate natural log of the determinant of Psi
# (numerically more reliable and also faster than using det or determinant)
LnDetPsi <- 2 * sum(log(abs(diag(cholPsi))))
## (L-5.1) Psi Inverted
#inverse with cholesky decomposed Psi
Psinv <- try(chol2inv(cholPsi), TRUE)
## give penalty if failed
if (inherits(Psinv, "try-error")) {
#warning("Correlation matrix is not positive semi-definite (During Maximum Likelihood Estimation in buildKriging). Returning penalty.")
penalty <-
penval - log10(min(eigen(
Psi, symmetric = TRUE, only.values = TRUE
)$values))
return(list(
NegLnLike = penalty,
Psi = NA,
Psinv = NA,
mu = NA,
ssq = NA
))
}
psisum <-
sum(Psinv)
# this sum of all matrix elements may sometimes become zero,
# which may be caused by inaccuracies. then, the following may help
if (psisum == 0) {
psisum <- as.numeric(rep(1, n) %*% Psinv %*% rep(1, n))
if (psisum == 0) {
#if it is still zero, return penalty
#warning("Sum of elements in inverse correlation matrix is zero (During Maximum Likelihood Estimation in buildKriging). Returning penalty.")
return(list(
NegLnLike = penval,
Psi = NA,
Psinv = NA,
mu = NA,
ssq = NA
))
}
}
## (L-5.2) Mean
mu <- sum(Psinv %*% Ay) / psisum
if (is.infinite(mu) | is.na(mu)) {
#warning("MLE estimate of mu is infinite or NaN (During Maximum Likelihood Estimation in buildKriging). Returning penalty.")
return(list(
NegLnLike = penval,
Psi = NA,
Psinv = NA,
mu = NA,
ssq = NA
))
}
## (L-5.3) yoneMu
yonemu <- Ay - mu
SigmaSqr <- (t(yonemu) %*% Psinv %*% yonemu) / n
if (SigmaSqr < 0) {
#warning("Maximum Likelihood Estimate of model parameter sigma^2 is negative (During Maximum Likelihood Estimation in buildKriging). Returning penalty. ")
return(list(
NegLnLike = penval - SigmaSqr,
Psi = NA,
Psinv = NA,
mu = NA,
ssq = NA
))
}
## (L-5.4) Log Likelihood
NegLnLike <- n * log(SigmaSqr) + LnDetPsi
if (is.na(NegLnLike) | is.infinite(NegLnLike)) {
return(list(
NegLnLike = 1e4,
Psi = NA,
Psinv = NA,
mu = NA,
ssq = NA
))
}
## (L-5.5) Compile Result
list(
NegLnLike = NegLnLike,
Psi = Psi,
Psinv = Psinv,
mu = mu,
yonemu = yonemu,
ssq = SigmaSqr
)
}
#' Predict Kriging Model
#'
#' Predict with Kriging model produced by \code{\link{buildKriging}}.
#'
#' @param object Kriging model (settings and parameters) of class \code{kriging}.
#' @param newdata design matrix to be predicted
#' @param ... not used
#'
#' @return list with predicted mean \code{y}, uncertainty / standard deviation \code{s} (optional)
#' and expected improvement \code{ei} (optional).
#' Whether \code{s} and \code{ei} are returned is specified by the vector of strings
#' \code{object$target}, which then contains \code{"s"} and \code{"ei"}.
#'
#'
#' @examples
#' \donttest{
#' ## Create design points
#' x <- cbind(runif(20)*15-5,runif(20)*15)
#' ## Compute observations at design points (for Branin function)
#' y <- funBranin(x)
#' ## Create model
#' fit <- buildKriging(x,y)
#' fit$target <- c("y","s","ei")
#' ## first estimate error with regressive predictor
#' predict(fit,x)
#' }
#' @seealso \code{\link{buildKriging}}
#' @export
#' @keywords internal
predict.kriging <- function(object, newdata, ...) {
x <- newdata
if (object$reinterpolate) {
return(predictKrigingReinterpolation(object, x, ...))
}
x <- repairNonNumeric(x, object$types) #round to nearest integer.
#normalize input x
x <-
normalizeMatrix2(data.matrix(x),
0,
1,
object$normalizexmin,
object$normalizexmax)
AX <- object$scaledx
theta <- object$dmodeltheta
Psinv <-
object$Psinv #fixed: does not need to be computed, is already done in likelihood function
n <- dim(AX)[1]
mu <- object$mu
yonemu <- object$yonemu
SigmaSqr <- object$ssq
k <- nrow(x)
nvar <- ncol(x)
if (object$optimizeP)
p <- object$P
else
p <- rep(2, nvar)
psi <- matrix(0, k, n)
for (i in 1:nvar) {
#todo nnn number variables
#psi[,i] <- colSums(theta*(abs(AX[i,]-t(x))^p))
tmp <- expand.grid(AX[, i], x[, i])
if (object$types[i] == "factor") {
tmp <- as.numeric(tmp[, 1] != tmp[, 2]) ^ p[i]
} else{
tmp <- abs(tmp[, 1] - tmp[, 2]) ^ p[i]
}
psi <- psi + theta[i] * matrix(tmp, k, n, byrow = T)
}
psi <- exp(-psi)
f <-
as.numeric(psi %*% (Psinv %*% yonemu)) + mu #vectorised #TODO: Psinv %*% yonemu can be computed ahead of time
res <- list(y = f)
if (any(object$target %in% c("s", "ei"))) {
lambda <- object$dmodellambda
SSqr <-
SigmaSqr * (1 + lambda - diag(psi %*% (Psinv %*% t(psi))))
s <- sqrt(abs(SSqr))
res$s <- s
if (any(object$target == "ei")) {
ei <- expectedImprovement(f, s, object$min)
# return y, s and ei
res <- list(y = f, s = s, ei = ei)
}
}
if (object$returnCrossCor)
res$psi <- psi
return(res)
}
#' Predict Kriging Model (Re-interpolating)
#'
#' Kriging predictor with re-interpolation to avoid stalling the optimization process
#' which employs this model as a surrogate.
#' This is supposed to be used with deterministic experiments,
#' which do need a non-interpolating model that avoids predicting non-zero error at sample locations.
#' This can be useful when the model is deterministic
#' (i.e. repeated evaluations of one parameter vector do not yield different values)
#' but does have a "noisy" structure (e.g. due to computational inaccuracies, systematical error).
#'
#' Please note that this re-interpolation implementation will not necessarily yield values of
#' exactly zero at the sample locations used for model building. Slight deviations can occur.
#'
#' @param object Kriging model (settings and parameters) of class \code{kriging}.
#' @param newdata design matrix to be predicted
#' @param ... not used
#'
#' @importFrom MASS ginv
#'
#' @return list with predicted mean \code{y}, uncertainty \code{s} (optional) and expected improvement \code{ei} (optional).
#' Whether \code{s} and \code{ei} are returned is specified by the vector of strings \code{object$target},
#' which then contains "s" and "ei.
#'
#' @examples
#' \donttest{
#' ## Create design points
#' x <- cbind(runif(20)*15-5,runif(20)*15)
#' ## Compute observations at design points (for Branin function)
#' y <- funBranin(x)
#' ## Create model
#' fit <- buildKriging(x,y,control=list(reinterpolate=FALSE))
#' fit$target <- c("y","s")
#' ## first estimate error with regressive predictor
#' sreg <- predict(fit,x)$s
#' ## now estimate error with re-interpolating predictor
#' sreint <- predictKrigingReinterpolation(fit,x)$s
#' ## equivalent:
#' fit$reinterpolate <- TRUE
#' sreint2 <- predict(fit,x)$s
#' print(sreg)
#' print(sreint)
#' print(sreint2)
#' ## sreint should be close to zero, significantly smaller than sreg
#' }
#' @seealso \code{\link{buildKriging}}, \code{\link{predict.kriging}}
#' @export
#' @keywords internal
predictKrigingReinterpolation <- function(object, newdata, ...) {
x <- newdata
x <- repairNonNumeric(x, object$types) #round to nearest integer.
#normalize input x
x <-
normalizeMatrix2(data.matrix(x),
0,
1,
object$normalizexmin,
object$normalizexmax)
AX <- object$scaledx
theta <- object$dmodeltheta
lambda <- object$dmodellambda
Psi <- object$Psi
Psinv <-
object$Psinv #fixed: does not need to be computed, is already done in likelihood function
n <- dim(AX)[1]
mu <- object$mu
yonemu <- object$yonemu
#
PsiB <-
Psi - diag(lambda, n) + diag(.Machine$double.eps, n) #TODO this and the following can be precomputed during model building
SigmaSqr <-
as.numeric(t(yonemu) %*% Psinv %*% PsiB %*% Psinv %*% yonemu) /
n
#
k <- nrow(x)
nvar <- ncol(x)
if (object$optimizeP)
p <- object$P
else
p <- rep(2, nvar)
psi <- matrix(0, k, n)
for (i in 1:nvar) {
#todo nnn number variables
#psi[,i] <- colSums(theta*(abs(AX[i,]-t(x))^p))
tmp <- expand.grid(AX[, i], x[, i])
if (object$types[i] == "factor") {
tmp <- as.numeric(tmp[, 1] != tmp[, 2]) ^ p[i]
} else{
tmp <- abs(tmp[, 1] - tmp[, 2]) ^ p[i]
}
psi <- psi + theta[i] * matrix(tmp, k, n, byrow = T)
}
psi <- exp(-psi)
f <- as.numeric(psi %*% (Psinv %*% yonemu)) + mu #vectorised
res <- list(y = f)
if (any(object$target %in% c("s", "ei"))) {
## Check whether PsiB is ill-conditioned
kap <- rcond(PsiB)
if (is.na(kap)) {
kap <- 0
}
if (kap < 1e-10) {
vmessage(0, "Warning. PsiB is ill conditioned",kap)
}
Psinv <- try(solve.default(PsiB), TRUE)
if (inherits(Psinv, "try-error")) {
#message("Trying ginv for PsiB")
Psinv <- try(ginv(PsiB), TRUE)
if (inherits(Psinv, "try-error")) {
message("Warning: Inverting PsiB failed. Using uncorrected Psinv.")
Psinv <- object$Psinv
}
}
SSqr <- SigmaSqr * (1 - diag(psi %*% (Psinv %*% t(psi))))
s <- sqrt(abs(SSqr))
res$s <- s
if (any(object$target == "ei")) {
ei <- expectedImprovement(f, s, object$min)
# return y, s and ei
res <- list(y = f, s = s, ei = ei)
}
}
if (object$returnCrossCor)
res$psi <- psi
res
}
#' Print Function Forrester Kriging
#'
#' Print information about a Forrester Kriging fit, as produced by \code{\link{buildKriging}}.
#'
#' @rdname print
#' @method print kriging
# @S3method print kriging
#' @param x fit returned by \code{\link{buildKriging}}.
#' @param ... additional parameters
#' @export
#' @keywords internal
print.kriging <- function(x, ...) {
cat("------------------------\n")
cat("Forrester Kriging model.\n")
cat("------------------------\n")
cat("Estimated activity parameters (theta) sorted \n")
cat("from most to least important variable \n")
cat(paste(
"x",
order(x$dmodeltheta, decreasing = TRUE),
sep = "",
collaps = " "
))
cat("\n")
cat(sort(x$dmodeltheta, decreasing = TRUE))
cat("\n \n")
cat("exponent(s) p:\n")
if (x$optimizeP)
cat(x$P)
else
cat(2)
cat("\n \n")
cat("Estimated regularization constant (or nugget) lambda:\n")
cat(x$dmodellambda)
cat("\n \n")
cat("Number of Likelihood evaluations during MLE:\n")
cat(x$nevals)
cat("\n")
cat("------------------------\n")
}
Try the SPOT package in your browser
Any scripts or data that you put into this service are public.
SPOT documentation built on June 26, 2022, 1:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions print.kriging predictKrigingReinterpolation predict.kriging krigingLikelihood normalizeMatrix buildKriging
Documented in buildKriging krigingLikelihood normalizeMatrix predict.kriging predictKrigingReinterpolation print.kriging
#' @title Build Kriging Model
#'
#' @description This function builds a Kriging model based on code by Forrester et al..
#' By default exponents (p) are fixed at a value of two, and a nugget (or regularization constant) is used.
#' To correct the uncertainty estimates in case of nugget, re-interpolation is also by default turned on.
#'
#' @details The model uses a Gaussian kernel: \code{k(x,z)=exp(-sum(theta_i * |x_i-z_i|^p_i))}. By default, \code{p_i = 2}.
#' Note that if dimension \code{x_i} is a factor variable (see parameter \code{types}), Hamming distance will be used
#' instead of \code{|x_i-z_i|}.
#'
#' @param x design matrix (sample locations)
#' @param y vector of observations at \code{x}
#' @param control (list), with the options for the model building procedure. Note: This can also be changed after the model has been built, by manipulating the respective \code{object$target} value.
#' \describe{
#' \item{\code{types}}{a character vector giving the data type of each variable.
#' All but "factor" will be handled as numeric, "factor" (categorical) variables
#' will be subject to the hamming distance.}
#' \item{\code{thetaLower}}{lower boundary for theta, default is \code{1e-4}}
#' \item{\code{thetaUpper}}{upper boundary for theta, default is \code{1e2}}
#' \item{\code{algTheta}}{algorithm used to find theta, default is \code{optimDE}.}
#' \item{\code{budgetAlgTheta}}{budget for the above mentioned algorithm,
#' default is \code{200}. The value will be multiplied with the length of
#' the model parameter vector to be optimized.}
#' \item{\code{optimizeP}}{boolean that specifies whether the exponents (\code{p}) should be optimized. Else they will be set to two. Default is \code{FALSE}.}
#' \item{\code{useLambda}}{whether or not to use the regularization constant lambda
#' (nugget effect). Default is \code{TRUE}.}
#' \item{\code{lambdaLower}}{lower boundary for log10{lambda}, default is \code{-6}}
#' \item{\code{lambdaUpper}}{upper boundary for log10{lambda}, default is \code{0}}
#' \item{\code{startTheta}}{optional start value for theta optimization, default is \code{NULL}}
#' \item{\code{reinterpolate}}{whether (\code{TRUE},default) or not (\code{FALSE}) reinterpolation
#' should be performed.}
#' \item{\code{target}}{values of the prediction, a vector of strings.
#' Each string specifies a value to be predicted, e.g., "y" for mean, "s" for standard deviation, "ei" for expected improvement.
#' See also \code{\link{predict.kriging}}.}
#' }
#'
#' @return an object of class \code{kriging}. Basically a list, with the options
#' and found parameters for the model which has to be passed to the predictor function:\cr
#' \code{x} sample locations (scaled to values between 0 and 1)\cr
#' \code{y} observations at sample locations (see parameters)\cr
#' \code{thetaLower} lower boundary for theta (see parameters)\cr
#' \code{thetaUpper} upper boundary for theta (see parameters)\cr
#' \code{algTheta} algorithm to find theta (see parameters)\cr
#' \code{budgetAlgTheta} budget for the above mentioned algorithm (see parameters)\cr
#' \code{optimizeP} boolean that specifies whether the exponents (\code{p}) were optimized (see parameters)\cr
#' \code{normalizeymin} minimum in normalized space\cr
#' \code{normalizeymax} maximum in normalized space\cr
#' \code{normalizexmin} minimum in input space\cr
#' \code{normalizexmax} maximum in input space\cr
#' \code{dmodeltheta} vector of activity parameters\cr
#' \code{Theta} log_10 vector of activity parameters (i.e. \code{log10(dmodeltheta)})\cr
#' \code{dmodellambda} regularization constant (nugget) \cr
#' \code{Lambda} log_10 of regularization constant (nugget) (i.e. \code{log10(dmodellambda)})\cr
#' \code{yonemu} \code{Ay-ones*mu} \cr
#' \code{ssq} sigma square\cr
#' \code{mu} mean mu\cr
#' \code{Psi} matrix large Psi\cr
#' \code{Psinv} inverse of Psi\cr
#' \code{nevals} number of Likelihood evaluations during MLE
#'
#' @importFrom stats predict
#' @importFrom stats dist
#' @importFrom stats var
#'
#' @export
#' @seealso \code{\link{predict.kriging}}
#' @references Forrester, Alexander I.J.; Sobester, Andras; Keane, Andy J. (2008). Engineering Design via Surrogate Modelling - A Practical Guide. John Wiley & Sons.
#'
#' @examples
#' ## Create design points
#' set.seed(1)
#' x <- cbind(runif(20)*15-5,runif(20)*15)
#' y <- funBranin(x)
#' ## Create model with default settings
#' fit <- buildKriging(x,y,control = list(algTheta=optimLHD))
#' ## Print model parameters
#' print(fit)
#' ## Predict at new location
#' predict(fit,cbind(1,2))
#' ## True value at location
#' funBranin(matrix(c(1,2), 1))
#' ##
#' ## Next Example: Handling factor variables
#' \donttest{
#' ## create a test function:
#' braninFunctionFactor <- function (x) {
#' y <- (x[2] - 5.1/(4 * pi^2) * (x[1] ^2) + 5/pi * x[1] - 6)^2 +
#' 10 * (1 - 1/(8 * pi)) * cos(x[1] ) + 10
#' if(x[3]==1)
#' y <- y +1
#' else if(x[3]==2)
#' y <- y -1
#' y
#' }
#' ## create training data
#' set.seed(1)
#' x <- cbind(runif(50)*15-5,runif(50)*15,sample(1:3,50,replace=TRUE))
#' y <- as.matrix(apply(x,1,braninFunctionFactor))
#' ## fit the model (default: assume all variables are numeric)
#' fitDefault <- buildKriging(x,y,control = list(algTheta=optimDE))
#' ## fit the model (give information about the factor variable)
#' fitFactor <- buildKriging(x,y,control =
#' list(algTheta=optimDE,types=c("numeric","numeric","factor")))
#' ## create test data
#' xtest <- cbind(runif(200)*15-5,runif(200)*15,sample(1:3,200,replace=TRUE))
#' ytest <- as.matrix(apply(xtest,1,braninFunctionFactor))
#' ## Predict test data with both models, and compute error
#' ypredDef <- predict(fitDefault,xtest)$y
#' ypredFact <- predict(fitFactor,xtest)$y
#' mean((ypredDef-ytest)^2)
#' mean((ypredFact-ytest)^2)
#' }
buildKriging <- function(x, y, control = list()) {
## (K-1) Set Parameters
x <-
data.matrix(x)
# data.matrix is better than as.matrix,
# because it always converts to numeric, not character! (in case of mixed data types.)
k <- ncol(x) # dimension
n <- nrow(x) # number of observations
con <- list(
thetaLower = 1e-4,
thetaUpper = 1e2,
types = rep("numeric", k),
algTheta = optimDE,
budgetAlgTheta = 200,
optimizeP = FALSE,
useLambda = TRUE,
lambdaLower = -6,
lambdaUpper = 0,
startTheta = NULL,
reinterpolate = TRUE,
target = "y"
)
con[names(control)] <- control
control <- con
fit <- control
fit$x <- x
fit$y <- y
LowerTheta <- rep(1, k) * log10(fit$thetaLower)
UpperTheta <- rep(1, k) * log10(fit$thetaUpper)
## (K-2) Normalize input data
ymin <- 0
ymax <- 1
fit$normalizeymin <- ymin
fit$normalizeymax <- ymax
# normalization in each column:
res <- normalizeMatrix(fit$x, ymin, ymax)
fit$scaledx <- res$y
fit$normalizexmin <- res$xmin
fit$normalizexmax <- res$xmax
## (K-3) Prepare distance/correlation matrix
A <- matrix(0, k, n * n)
for (i in 1:k) {
if (control$types[i] != "factor") {
A[i, ] <-
as.numeric(as.matrix(dist(fit$scaledx[, i]))) #euclidean distance
} else {
tmp <-
outer(fit$scaledx[, i], fit$scaledx[, i], '!=') #hamming distance
class(tmp) <- "numeric"
A[i, ] <- tmp
}
}
## (K-4) Prepare starting points, search bounds and penalty value for MLE optimization
## 4.1 theta
if (is.null(fit$startTheta))
x1 <- rep(n / (100 * k), k) # start point for theta
else
x1 <- fit$startTheta
## 4.2 p
if (fit$optimizeP) {
LowerTheta <- c(LowerTheta, rep(1, k) * 0.01)
UpperTheta <- c(UpperTheta, rep(1, k) * 2)
x3 <- rep(1, k) * 1.9 #start values for p
x0 <- c(x1, x3)
} else{
# p is fixed to 2 and the array A is completely precalculated
A <- A ^ 2
x0 <- c(x1)
}
## 4.3 lambda
if (fit$useLambda) {
# start value for lambda:
x2 <- (fit$lambdaUpper + fit$lambdaLower) / 2
x0 <- c(x0, x2)
#append regression constant lambda (nugget)
LowerTheta <- c(LowerTheta, fit$lambdaLower)
UpperTheta <- c(UpperTheta, fit$lambdaUpper)
}
#force x0 into bounds
x0 <- pmin(x0, UpperTheta)
x0 <- pmax(x0, LowerTheta)
x0 <- matrix(x0, 1) #matrix with one row
opts <- list(funEvals = fit$budgetAlgTheta * ncol(x0))
## 4.4 penalty
#determine a good penalty value (based on number of samples and variance of y)
penval <- n * log(var(y)) + 1e4
## (K-5) MLE objective function
# Wrapper for optimizing theta based on krigingLikelihood:
fitFun <-
function (x, fX, fy, optimizeP, useLambda, penval) {
#todo vectorize, at least for cma_es with active vectorize?
as.numeric(krigingLikelihood(x, fX, fy, optimizeP, useLambda, penval)$NegLnLike)
}
## (K-7) MLE optimization
res <- fit$algTheta(
x = x0,
fun =
function(x, fX, fy, optimizeP, useLambda, penval) {
if (!is.matrix(x)) {
fitFun(x, fX, fy, optimizeP, useLambda, penval)
}
else{
apply(x, 1, fitFun, fX, fy, optimizeP, useLambda, penval)
}
},
lower = LowerTheta,
upper = UpperTheta,
control = opts,
fX = A,
fy = fit$y,
optimizeP = fit$optimizeP,
useLambda = fit$useLambda,
penval = penval
)
if (is.null(res$xbest))
res$xbest <- x0
## (K-8): compile results from MLE optimization (to fit object)
Params <- res$xbest
nevals <- as.numeric(res$count[[1]])
fit$Theta <- Params[1:k]
fit$dmodeltheta <- 10 ^ Params[1:k]
if (fit$optimizeP) {
fit$P <- Params[(k + 1):(2 * k)]
}
if (fit$useLambda) {
fit$Lambda <- Params[length(Params)]
fit$dmodellambda <- 10 ^ Params[length(Params)]
} else{
fit$Lambda <- -Inf
fit$dmodellambda <- 0
}
## (K-9) Evaluate with optimized parameters
res <-
krigingLikelihood(c(fit$Theta, fit$P, fit$Lambda),
A,
fit$y,
fit$optimizeP,
fit$useLambda)
if (is.na(res$Psinv[1]))
stop(
"buildKriging failed to produce a valid correlation matrix. Consider activating the nugget/regularization constant via useLambda=TRUE in the control list."
)
## (K-10) Add results from MLE evaluation to fit object
fit$yonemu <- res$yonemu
fit$ssq <- as.numeric(res$ssq)
fit$mu <- res$mu
fit$Psi <- res$Psi
fit$Psinv <- res$Psinv
fit$nevals <- nevals
fit$like <- res$NegLnLike
fit$returnCrossCor <- FALSE
## (K-11) Calculate observed minimum
xlist <- split(x, 1:nrow(x))
uniquex <- unique(xlist)
ymean <- NULL
for (xi in uniquex) {
ind <- xlist %in% list(xi)
ymean <- c(ymean, mean(y[ind]))
}
fit$min <- min(ymean)
class(fit) <- "kriging"
return(fit)
}
#' Normalize design 2
#'
#' Normalize design with given maximum and minimum in input space.
#' Supportive function for Kriging model, not to be used directly.
#'
#' @param x design matrix in input space (n rows for each point, k columns for each parameter)
#' @param ymin minimum vector of normalized space
#' @param ymax maximum vector of normalized space
#' @param xmin minimum vector of input space
#' @param xmax maximum vector of input space
#'
#' @return normalized design matrix
#' @seealso \code{\link{buildKriging}}
#' @export
normalizeMatrix2 <- function (x, ymin, ymax, xmin, xmax) {
rangex <- xmax - xmin
rangey <- ymax - ymin
s <- dim(x)[1]
rangey * (x - matrix(rep(xmin, s), nrow = s, byrow = TRUE)) / matrix(rep(rangex, s), nrow =
s, byrow = TRUE) + ymin
}
#' @title Normalize design matrix
#'
#' @description Normalize design by using minimum and maximum of the design values
#' for input space. Each column has entries in the range from \code{ymin} to \code{ymax}.
#'
#' @param x design matrix in input space
#' @param ymin minimum vector of normalized space
#' @param ymax maximum vector of normalized space
#' @param MARGIN a vector giving the subscripts which the function will be applied over.
#' 1 indicates rows, 2 indicates columns. Default: \code{2}.
#'
#' @return \code{list} with the following entries:
#' \describe{
#' \item{\code{y}}{normalized design matrix in the range [ymin, ymax]}
#' \item{\code{xmin}}{min in each column}
#' \item{\code{xmax}}{max in each column}
#' }
#'
#' @seealso \code{\link{buildKriging}}
#' @examples
#' set.seed(1)
#' x <- matrix(c(rep(1,3), rep(2,3),rep(3,3), rep(4,3)),3,4)
#' ## columnwise:
#' normalizeMatrix(x, ymin=0, ymax=1)
#' ## rowwise
#' normalizeMatrix(x, ymin=0, ymax=1, MARGIN=1)
#' # rows with identical values are mapped to the mean:
#' x <- matrix(rep(0,4),2,2)
#' normalizeMatrix(x, ymin=0, ymax=1)
#'
#' @export
normalizeMatrix <- function(x, ymin, ymax, MARGIN = 2) {
if (MARGIN == 1)
x <- t(x)
# Return the maximum from each row:
xmax <- apply(x, 2, max)
# Return the minimum from each row:
xmin <- apply(x, 2, min)
s <- dim(x)[1]
rangex <- xmax - xmin
rangey <- ymax - ymin
xmin[rangex == 0] <- xmin[rangex == 0] - 0.5
xmax[rangex == 0] <- xmax[rangex == 0] + 0.5
rangex[rangex == 0] <- 1
y <-
rangey * (x - matrix(rep(xmin, s), nrow = s, byrow = TRUE)) / matrix(rep(rangex, s), nrow =
s, byrow = TRUE) + ymin
if (MARGIN == 1)
y <- t(y)
list(y = y, xmin = xmin, xmax = xmax)
}
#' Calculate negative log-likelihood
#'
#' Used to determine theta/lambda values for the Kriging model in \code{\link{buildKriging}}.
#' Supportive function for Kriging model, not to be used directly.
#'
#' @param x vector, containing parameters log10(theta), log10(lambda) and p.
#' @param AX 3 dimensional array, constructed by buildKriging from the sample locations
#' @param Ay vector of observations at sample locations
#' @param optimizeP boolean, whether to optimize parameter p (exponents) or fix at two.
#' @param useLambda boolean, whether to use nugget
#' @param penval a penalty value which affects the value returned for invalid correlation matrices / configurations
#'
#' @return list with elements\cr
#' \code{NegLnLike} concentrated log-likelihood *-1 for minimising \cr
#' \code{Psi} correlation matrix\cr
#' \code{Psinv} inverse of correlation matrix (to save computation time during prediction)\cr
#' \code{mu} \cr
#' \code{ssq}
#' @seealso \code{\link{buildKriging}}
#' @export
#' @keywords internal
krigingLikelihood <-
function(x,
AX,
Ay,
optimizeP = FALSE,
useLambda = TRUE,
penval = 1e8) {
## note: this penalty value should not be a hard constant.
## the scale of the likelihood (n*log(SigmaSqr) + LnDetPsi)
## at least depends on log(var) and the number of samples
## hence, large number of samples may lead to cases where the
## penality is lower than the likelihood of most valid parameterizations
## suggested penalty:
#penval <- n*log(var(y)) + 1e4
## currently, this better penalty is set in the buildKriging function,
## when calling krigingLikelihood
## (L-1) Starting Points
nx <- nrow(AX)
theta <- 10 ^ x[1:nx]
if (optimizeP) {
AX <- abs(AX) ^ (x[(nx + 1):(2 * nx)])
}
lambda <- 0
if (useLambda) {
lambda <- 10 ^ x[length(x)]
}
if (any(theta == 0) ||
any(is.infinite(c(theta, lambda)))) {
#for instance caused by bound violation
return(list(
NegLnLike = 1e4,
Psi = NA,
Psinv = NA,
mu = NA,
ssq = NA
))
}
n <- dim(Ay)[1]
## (L-2) Correlation Matrix Psi
Psi <- exp(-matrix(colSums(theta * AX), n, n))
if (useLambda) {
Psi <- Psi + diag(lambda, n)
}
if (any(is.infinite(Psi))) {
# this is required especially if distance matrices are forced to be CNSD/NSD and hence have zero distances
penalty <- penval
return(list(
NegLnLike = penalty,
Psi = NA,
Psinv = NA,
mu = NA,
ssq = NA
)) #todo: add something like smallest eigenvalue to penalty?
}
## Check whether Psi is ill-conditioned
kap <- rcond(Psi)
if (is.na(kap)) {
kap <- 0
}
if (kap < 1e-10) {
#warning("Correlation matrix is ill-conditioned (During Maximum Likelihood Estimation in buildKriging). Returning penalty.")
return(
list(
NegLnLike = penval - log10(kap),
Psi = NA,
Psinv = NA,
mu = NA,
ssq = NA,
a = NA,
U = NA,
isIndefinite = TRUE
)
)
}
## (L-3) cholesky decomposition
cholPsi <- try(chol(Psi), TRUE)
## give penalty if fail
if (inherits(cholPsi, "try-error")) {
#warning("Correlation matrix is not positive semi-definite (During Maximum Likelihood Estimation in buildKriging). Returning penalty.")
penalty <-
penval - log10(min(eigen(
Psi, symmetric = TRUE, only.values = TRUE
)$values))
return(list(
NegLnLike = penalty,
Psi = NA,
Psinv = NA,
mu = NA,
ssq = NA
))
}
## (L-4) Determininant
## calculate natural log of the determinant of Psi
# (numerically more reliable and also faster than using det or determinant)
LnDetPsi <- 2 * sum(log(abs(diag(cholPsi))))
## (L-5.1) Psi Inverted
#inverse with cholesky decomposed Psi
Psinv <- try(chol2inv(cholPsi), TRUE)
## give penalty if failed
if (inherits(Psinv, "try-error")) {
#warning("Correlation matrix is not positive semi-definite (During Maximum Likelihood Estimation in buildKriging). Returning penalty.")
penalty <-
penval - log10(min(eigen(
Psi, symmetric = TRUE, only.values = TRUE
)$values))
return(list(
NegLnLike = penalty,
Psi = NA,
Psinv = NA,
mu = NA,
ssq = NA
))
}
psisum <-
sum(Psinv)
# this sum of all matrix elements may sometimes become zero,
# which may be caused by inaccuracies. then, the following may help
if (psisum == 0) {
psisum <- as.numeric(rep(1, n) %*% Psinv %*% rep(1, n))
if (psisum == 0) {
#if it is still zero, return penalty
#warning("Sum of elements in inverse correlation matrix is zero (During Maximum Likelihood Estimation in buildKriging). Returning penalty.")
return(list(
NegLnLike = penval,
Psi = NA,
Psinv = NA,
mu = NA,
ssq = NA
))
}
}
## (L-5.2) Mean
mu <- sum(Psinv %*% Ay) / psisum
if (is.infinite(mu) | is.na(mu)) {
#warning("MLE estimate of mu is infinite or NaN (During Maximum Likelihood Estimation in buildKriging). Returning penalty.")
return(list(
NegLnLike = penval,
Psi = NA,
Psinv = NA,
mu = NA,
ssq = NA
))
}
## (L-5.3) yoneMu
yonemu <- Ay - mu
SigmaSqr <- (t(yonemu) %*% Psinv %*% yonemu) / n
if (SigmaSqr < 0) {
#warning("Maximum Likelihood Estimate of model parameter sigma^2 is negative (During Maximum Likelihood Estimation in buildKriging). Returning penalty. ")
return(list(
NegLnLike = penval - SigmaSqr,
Psi = NA,
Psinv = NA,
mu = NA,
ssq = NA
))
}
## (L-5.4) Log Likelihood
NegLnLike <- n * log(SigmaSqr) + LnDetPsi
if (is.na(NegLnLike) | is.infinite(NegLnLike)) {
return(list(
NegLnLike = 1e4,
Psi = NA,
Psinv = NA,
mu = NA,
ssq = NA
))
}
## (L-5.5) Compile Result
list(
NegLnLike = NegLnLike,
Psi = Psi,
Psinv = Psinv,
mu = mu,
yonemu = yonemu,
ssq = SigmaSqr
)
}
#' Predict Kriging Model
#'
#' Predict with Kriging model produced by \code{\link{buildKriging}}.
#'
#' @param object Kriging model (settings and parameters) of class \code{kriging}.
#' @param newdata design matrix to be predicted
#' @param ... not used
#'
#' @return list with predicted mean \code{y}, uncertainty / standard deviation \code{s} (optional)
#' and expected improvement \code{ei} (optional).
#' Whether \code{s} and \code{ei} are returned is specified by the vector of strings
#' \code{object$target}, which then contains \code{"s"} and \code{"ei"}.
#'
#'
#' @examples
#' \donttest{
#' ## Create design points
#' x <- cbind(runif(20)*15-5,runif(20)*15)
#' ## Compute observations at design points (for Branin function)
#' y <- funBranin(x)
#' ## Create model
#' fit <- buildKriging(x,y)
#' fit$target <- c("y","s","ei")
#' ## first estimate error with regressive predictor
#' predict(fit,x)
#' }
#' @seealso \code{\link{buildKriging}}
#' @export
#' @keywords internal
predict.kriging <- function(object, newdata, ...) {
x <- newdata
if (object$reinterpolate) {
return(predictKrigingReinterpolation(object, x, ...))
}
x <- repairNonNumeric(x, object$types) #round to nearest integer.
#normalize input x
x <-
normalizeMatrix2(data.matrix(x),
0,
1,
object$normalizexmin,
object$normalizexmax)
AX <- object$scaledx
theta <- object$dmodeltheta
Psinv <-
object$Psinv #fixed: does not need to be computed, is already done in likelihood function
n <- dim(AX)[1]
mu <- object$mu
yonemu <- object$yonemu
SigmaSqr <- object$ssq
k <- nrow(x)
nvar <- ncol(x)
if (object$optimizeP)
p <- object$P
else
p <- rep(2, nvar)
psi <- matrix(0, k, n)
for (i in 1:nvar) {
#todo nnn number variables
#psi[,i] <- colSums(theta*(abs(AX[i,]-t(x))^p))
tmp <- expand.grid(AX[, i], x[, i])
if (object$types[i] == "factor") {
tmp <- as.numeric(tmp[, 1] != tmp[, 2]) ^ p[i]
} else{
tmp <- abs(tmp[, 1] - tmp[, 2]) ^ p[i]
}
psi <- psi + theta[i] * matrix(tmp, k, n, byrow = T)
}
psi <- exp(-psi)
f <-
as.numeric(psi %*% (Psinv %*% yonemu)) + mu #vectorised #TODO: Psinv %*% yonemu can be computed ahead of time
res <- list(y = f)
if (any(object$target %in% c("s", "ei"))) {
lambda <- object$dmodellambda
SSqr <-
SigmaSqr * (1 + lambda - diag(psi %*% (Psinv %*% t(psi))))
s <- sqrt(abs(SSqr))
res$s <- s
if (any(object$target == "ei")) {
ei <- expectedImprovement(f, s, object$min)
# return y, s and ei
res <- list(y = f, s = s, ei = ei)
}
}
if (object$returnCrossCor)
res$psi <- psi
return(res)
}
#' Predict Kriging Model (Re-interpolating)
#'
#' Kriging predictor with re-interpolation to avoid stalling the optimization process
#' which employs this model as a surrogate.
#' This is supposed to be used with deterministic experiments,
#' which do need a non-interpolating model that avoids predicting non-zero error at sample locations.
#' This can be useful when the model is deterministic
#' (i.e. repeated evaluations of one parameter vector do not yield different values)
#' but does have a "noisy" structure (e.g. due to computational inaccuracies, systematical error).
#'
#' Please note that this re-interpolation implementation will not necessarily yield values of
#' exactly zero at the sample locations used for model building. Slight deviations can occur.
#'
#' @param object Kriging model (settings and parameters) of class \code{kriging}.
#' @param newdata design matrix to be predicted
#' @param ... not used
#'
#' @importFrom MASS ginv
#'
#' @return list with predicted mean \code{y}, uncertainty \code{s} (optional) and expected improvement \code{ei} (optional).
#' Whether \code{s} and \code{ei} are returned is specified by the vector of strings \code{object$target},
#' which then contains "s" and "ei.
#'
#' @examples
#' \donttest{
#' ## Create design points
#' x <- cbind(runif(20)*15-5,runif(20)*15)
#' ## Compute observations at design points (for Branin function)
#' y <- funBranin(x)
#' ## Create model
#' fit <- buildKriging(x,y,control=list(reinterpolate=FALSE))
#' fit$target <- c("y","s")
#' ## first estimate error with regressive predictor
#' sreg <- predict(fit,x)$s
#' ## now estimate error with re-interpolating predictor
#' sreint <- predictKrigingReinterpolation(fit,x)$s
#' ## equivalent:
#' fit$reinterpolate <- TRUE
#' sreint2 <- predict(fit,x)$s
#' print(sreg)
#' print(sreint)
#' print(sreint2)
#' ## sreint should be close to zero, significantly smaller than sreg
#' }
#' @seealso \code{\link{buildKriging}}, \code{\link{predict.kriging}}
#' @export
#' @keywords internal
predictKrigingReinterpolation <- function(object, newdata, ...) {
x <- newdata
x <- repairNonNumeric(x, object$types) #round to nearest integer.
#normalize input x
x <-
normalizeMatrix2(data.matrix(x),
0,
1,
object$normalizexmin,
object$normalizexmax)
AX <- object$scaledx
theta <- object$dmodeltheta
lambda <- object$dmodellambda
Psi <- object$Psi
Psinv <-
object$Psinv #fixed: does not need to be computed, is already done in likelihood function
n <- dim(AX)[1]
mu <- object$mu
yonemu <- object$yonemu
#
PsiB <-
Psi - diag(lambda, n) + diag(.Machine$double.eps, n) #TODO this and the following can be precomputed during model building
SigmaSqr <-
as.numeric(t(yonemu) %*% Psinv %*% PsiB %*% Psinv %*% yonemu) /
n
#
k <- nrow(x)
nvar <- ncol(x)
if (object$optimizeP)
p <- object$P
else
p <- rep(2, nvar)
psi <- matrix(0, k, n)
for (i in 1:nvar) {
#todo nnn number variables
#psi[,i] <- colSums(theta*(abs(AX[i,]-t(x))^p))
tmp <- expand.grid(AX[, i], x[, i])
if (object$types[i] == "factor") {
tmp <- as.numeric(tmp[, 1] != tmp[, 2]) ^ p[i]
} else{
tmp <- abs(tmp[, 1] - tmp[, 2]) ^ p[i]
}
psi <- psi + theta[i] * matrix(tmp, k, n, byrow = T)
}
psi <- exp(-psi)
f <- as.numeric(psi %*% (Psinv %*% yonemu)) + mu #vectorised
res <- list(y = f)
if (any(object$target %in% c("s", "ei"))) {
## Check whether PsiB is ill-conditioned
kap <- rcond(PsiB)
if (is.na(kap)) {
kap <- 0
}
if (kap < 1e-10) {
vmessage(0, "Warning. PsiB is ill conditioned",kap)
}
Psinv <- try(solve.default(PsiB), TRUE)
if (inherits(Psinv, "try-error")) {
#message("Trying ginv for PsiB")
Psinv <- try(ginv(PsiB), TRUE)
if (inherits(Psinv, "try-error")) {
message("Warning: Inverting PsiB failed. Using uncorrected Psinv.")
Psinv <- object$Psinv
}
}
SSqr <- SigmaSqr * (1 - diag(psi %*% (Psinv %*% t(psi))))
s <- sqrt(abs(SSqr))
res$s <- s
if (any(object$target == "ei")) {
ei <- expectedImprovement(f, s, object$min)
# return y, s and ei
res <- list(y = f, s = s, ei = ei)
}
}
if (object$returnCrossCor)
res$psi <- psi
res
}
#' Print Function Forrester Kriging
#'
#' Print information about a Forrester Kriging fit, as produced by \code{\link{buildKriging}}.
#'
#' @rdname print
#' @method print kriging
# @S3method print kriging
#' @param x fit returned by \code{\link{buildKriging}}.
#' @param ... additional parameters
#' @export
#' @keywords internal
print.kriging <- function(x, ...) {
cat("------------------------\n")
cat("Forrester Kriging model.\n")
cat("------------------------\n")
cat("Estimated activity parameters (theta) sorted \n")
cat("from most to least important variable \n")
cat(paste(
"x",
order(x$dmodeltheta, decreasing = TRUE),
sep = "",
collaps = " "
))
cat("\n")
cat(sort(x$dmodeltheta, decreasing = TRUE))
cat("\n \n")
cat("exponent(s) p:\n")
if (x$optimizeP)
cat(x$P)
else
cat(2)
cat("\n \n")
cat("Estimated regularization constant (or nugget) lambda:\n")
cat(x$dmodellambda)
cat("\n \n")
cat("Number of Likelihood evaluations during MLE:\n")
cat(x$nevals)
cat("\n")
cat("------------------------\n")
}
Try the SPOT package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.