max_sur: Minimizer of the SUR criterion

Description

Minimization, based on the package rgenoud (or on exhaustive search on a discrete set), of the SUR criterion. If not provided, the integration points in X space are generated before optimization by latin hypercube sampling. The integration points in Y space are fixed. The default number of integration points in X is 100 times the dimension of the problem.

Usage

max_sur(lower, upper, parinit=NULL, control=NULL, discrete.X=NULL, 
n.int.y=10, integration.points=NULL, T, model, new.noise.var=0, type="UK")

Arguments

lower

vector containing the lower bounds of the variables to be optimized over

upper

vector containing the upper bounds of the variables to be optimized over

parinit

optional vector containing the initial values for the variables to be optimized over

control

optional list of control parameters for optimization. One can control "pop.size" (default : [4+3*log(nb of variables)]), "max.generations" (5), "wait.generations" (2) and "BFGSburnin" (0) of fun

target value (a real number)

model

An object of class km (Kriging model)

new.noise.var

optional scalar value of the noise variance for the new observations

type

Kriging type (string): "SK" or "UK" (default)

n.int.y

optional scalar value of the the number of integration points in the y space.

discrete.X

optional matrix of candidate points. If provided, the search for new observations is made on this discrete set instead of running the continuous optimisation, and it is also used as integration points

integration.points

optional scalar or matrix. If it is scalar, it defines the number of integration points, then generated as an LHS design; if it is a matrix, it defines the integration points directly. If discrete.X is provided, thi

Value

A list with components:
parthe best set of parameters found.
valuethe value of the SUR criterion at par.

References

Bect J., Ginsbourger D., Li L., Picheny V., Vazquez E. (2010), Sequential design of computer experiments for the estimation of a probability of failure, accepted with minor revisions to the Journal of Statistics and Computing, http://arxiv.org/abs/1009.5177 Vazquez, E., Bect, J.: A sequential Bayesian algorithm to estimate a probability of failure. In: Proceedings of the 15th IFAC Symposium on System Identification, (SYSID 2009), Saint-Malo, France (2009)

Examples

Run this code

##################################################################
#a 9-point full factorial initial design
design.fact <- expand.grid(seq(0,1,length=3), seq(0,1,length=3))

design.fact <- data.frame(design.fact)
names(design.fact) <- c ( "x1","x2")
testfun <- camelback2			#our test function

#the response
response <- testfun(design.fact)

#the initial km model
model <- km(formula=~., design = design.fact, response = response, 
covtype="matern5_2")

#the integration points
n.grid <- 20
x.grid <- y.grid <- seq(0,1,length=n.grid)
design.grid <- expand.grid(x.grid, y.grid)

#evaluate criterion on the grid
T <- 0
		
sur.grid <- apply(design.grid, 1, sur_optim, T=T, model=model, 
integration.points=design.grid)
z.grid <- matrix(sur.grid, n.grid, n.grid)

#plots: contour of the criterion, doe points and new point
contour(x.grid,y.grid,z.grid,25)
points(design.fact, col="black", pch=20, lwd=4)

#plots: contour of the actual function at threshold
testfun.grid <- testfun(design.grid)
z.grid.2 <- matrix(testfun.grid, n.grid, n.grid)
contour(x.grid,y.grid,z.grid.2,levels=T,col="blue",add=TRUE)
title("Contour lines of sur criterion (black) and of f(x)=T (blue)")

#search best point with Genoud
opt <- max_sur(lower=c(0,0), upper=c(1,1), T=T, model=model, 
integration.points=design.grid)
points(opt$par, col="blue", pch=20, lwd=4)
              
###################################################################