R/EI.R

In DiceOptim: Kriging-Based Optimization for Computer Experiments

#' Analytical expression of the Expected Improvement criterion
#' 
#' Computes the Expected Improvement at current location. The current minimum
#' of the observations can be replaced by an arbitrary value (plugin), which
#' is usefull in particular in noisy frameworks.
#' 
#' 
#' @param x a vector representing the input for which one wishes to calculate
#' EI,
#' @param model an object of class \code{\link[DiceKriging]{km}},
#' @param plugin optional scalar: if provided, it replaces the minimum of the
#' current observations,
#' @param type "SK" or "UK" (by default), depending whether uncertainty
#' related to trend estimation has to be taken into account,
#' @param minimization logical specifying if EI is used in minimiziation or in
#' maximization,
#' @param envir an optional environment specifying where to assign
#' intermediate values for future gradient calculations. Default is NULL.
#' @param proxy an optional Boolean, if TRUE EI is replaced by the kriging mean (to minimize)
#' @return The expected improvement, defined as \deqn{EI(x) := E[( min(Y(X)) -
#' Y(x))^{+} | Y(X)=y(X)],} where X is the current design of experiments and Y
#' is the random process assumed to have generated the objective function y.
#' If a plugin is specified, it replaces \deqn{min(Y(X))} in the previous
#' formula.
#' @author David Ginsbourger 
#' 
#' Olivier Roustant 
#' 
#' Victor Picheny 
#' @seealso \code{\link{max_EI}}, \code{\link{EGO.nsteps}}, \code{\link{qEI}}
#' @references
#' 
#' D.R. Jones, M. Schonlau, and W.J. Welch (1998), Efficient global
#' optimization of expensive black-box functions, \emph{Journal of Global
#' Optimization}, 13, 455-492.
#' 
#' J. Mockus (1988), \emph{Bayesian Approach to Global Optimization}. Kluwer
#' academic publishers.
#' 
#' T.J. Santner, B.J. Williams, and W.J. Notz (2003), \emph{The design and
#' analysis of computer experiments}, Springer.
#' 
#' M. Schonlau (1997), \emph{Computer experiments and global optimization},
#' Ph.D. thesis, University of Waterloo.
#' @keywords models
#' @examples
#' 
#' set.seed(123)
#' ##########################################################################
#' ### 	  EI SURFACE ASSOCIATED WITH AN ORDINARY KRIGING MODEL        ####
#' ###    OF THE BRANIN FUNCTION KNOWN AT A 9-POINTS FACTORIAL DESIGN    ####
#' ##########################################################################
#' 
#' # a 9-points factorial design, and the corresponding response
#' d <- 2; n <- 9
#' design.fact <- expand.grid(seq(0,1,length=3), seq(0,1,length=3))
#' names(design.fact)<-c("x1", "x2")
#' design.fact <- data.frame(design.fact) 
#' names(design.fact)<-c("x1", "x2")
#' response.branin <- apply(design.fact, 1, branin)
#' response.branin <- data.frame(response.branin) 
#' names(response.branin) <- "y" 
#' 
#' # model identification
#' fitted.model1 <- km(~1, design=design.fact, response=response.branin, 
#' covtype="gauss", control=list(pop.size=50,trace=FALSE), parinit=c(0.5, 0.5))
#' 
#' # graphics
#' n.grid <- 12
#' x.grid <- y.grid <- seq(0,1,length=n.grid)
#' design.grid <- expand.grid(x.grid, y.grid)
#' #response.grid <- apply(design.grid, 1, branin)
#' EI.grid <- apply(design.grid, 1, EI,fitted.model1)
#' z.grid <- matrix(EI.grid, n.grid, n.grid)
#' contour(x.grid,y.grid,z.grid,25)
#' title("Expected Improvement for the Branin function known at 9 points")
#' points(design.fact[,1], design.fact[,2], pch=17, col="blue")
#' 
#' @export EI
EI <- function (x, model, plugin=NULL, type="UK", minimization = TRUE, envir=NULL, proxy=FALSE) {

  if (is.null(plugin)){ 
    if (minimization) {
      plugin <- min(model@y)
    } else {
      plugin <- -max(model@y)
    }
  }
  m <- plugin
  
  ########################################################################################
  # Convert x in proper format(s)
  if (is.data.frame(x)) {
    d <- length(x)
    if (d != model@d){ stop("x does not have the right size") }
    newdata.num <- as.numeric(x)
    newdata <- data.frame(t(newdata.num))
  } else {
    if (is.null(dim(x))) {
      d <- length(x)
      if (d != model@d){ stop("x does not have the right size") }
      newdata <- data.frame(t(as.numeric(x)))
    } else {
      d <- ncol(x)
      if (d != model@d){ stop("x does not have the right size") }
      newdata <- data.frame(x)
    }
  }
  colnames(newdata) <- colnames(model@X)
  
  ########################################################################################
  if (!is.null(envir)) {
    predx <- predict(object=model, newdata=newdata, type=type, checkNames = FALSE)
  } else {
    predx <- predict(object=model, newdata=newdata, type=type, checkNames = FALSE, light.return=TRUE)
  }
  
  kriging.mean <- predx$mean
  if(!minimization) {
    kriging.mean <- -kriging.mean
  }
  kriging.sd   <- predx$sd
  
  if (proxy) {
    xcr <- xcr.prob <- xcr.dens <- NULL
    res <- m - kriging.mean
  } else {
    xcr <- (m - kriging.mean) / kriging.sd
    
    xcr.prob <- pnorm(xcr)
    xcr.dens <- dnorm(xcr)	        
    res <- (m - kriging.mean) * xcr.prob + kriging.sd * xcr.dens
  }
  
  ########################################################################################  
  if (!is.null(envir)) 
  { assign("xcr", xcr, envir=envir)
    assign("xcr.prob", xcr.prob, envir=envir)
    assign("xcr.dens", xcr.dens, envir=envir)
    assign("kriging.sd", kriging.sd, envir=envir)
    assign("c", predx$c, envir=envir)
    assign("Tinv.c", predx$Tinv.c, envir=envir)
  }
  return(res)
}