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.

Usage

max_sur(lower, upper, optimcontrol = NULL, 
integration.param = NULL, T, model, 
new.noise.var = 0, real.volume.variance = FALSE, real.volume.constant = FALSE)

Arguments

lower

Vector containing the lower bounds of the design space.

upper

Vector containing the upper bounds of the design space.

optimcontrol

Optional list of control parameters for the optimization of the sampling criterion. The field method defines which optimization method is used: it can be either "genoud" (default) for an optimisation using the genoud algorithm, or "discrete" for an optimisation over a specified discrete set. If the field method is set to "genoud", one can set some parameters of this algorithm: pop.size (default : 50*d), max.generations (10*d), wait.generations (2), BFGSburnin (2) and the mutations P1, P2, up to P9 (see genoud). Numbers into brackets are the default values. If the field method is set to "discrete", one can set the field optim.points: p * d matrix corresponding to the p points where the criterion will be evaluated. If nothing is specified, 100*d points are chosen randomly.

integration.param

Optional list of control parameter for the computation of integrals, containing the fields integration.points: a p*d matrix corresponding to p integrations points and integration.weights: a vector of size p corresponding to the weights of these integration points. If nothing is specified, default values are used (see: function integration_design for more details).

Target value (scalar).

model

A Kriging model of km class.

new.noise.var

Optional scalar value of the noise variance at the new observation.

real.volume.variance

Boolean to decide if the "jn" sampling criterion should be used instead of the "sur" criterion. When it is equal to FALSE (default), the "sur" criterion is used (faster computation).

real.volume.constant

When the"jn" criterion is chosen, this argument decides whether a constant part of the "jn" criterion should be computed or not. Computing this constant does NOT change the optimum of the "jn" criterion. Default value: FALSE. This argument is ignored if the argument real.volume.variance is set to FALSE.

Value

A list with components:

par

Best set of parameters found.

value

Value of the criterion at par.

allvalues

If an optimization on a discrete set of points is chosen, the value of the criterion at all these points

variance.volume

If the arguments real.volume.variance and real.volume.constant are both set to TRUE,the value of the computed constant

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, Statistics and Computing, pp.1-21, 2011, 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)

Chevalier C., Bect J., Ginsbourger D., Vazquez E., Picheny V., Richet Y. (2011), Fast parallel kriging-based stepwise uncertainty reduction with application to the identification of an excursion set ,http://hal.archives-ouvertes.fr/hal-00641108/

Examples

Run this code

# NOT RUN {
#max_sur

set.seed(8)
N <- 9 #number of observations
T <- 80 #threshold
testfun <- branin
lower <- c(0,0)
upper <- c(1,1)

#a 9 points initial design
design <- data.frame( matrix(runif(2*N),ncol=2) )
response <- testfun(design)

#km object with matern3_2 covariance
#params estimated by ML from the observations
model <- km(formula=~., design = design, 
	response = response,covtype="matern3_2")

optimcontrol <- list(method="genoud",pop.size=50)
integcontrol <- list(distrib="sur",n.points=50,init.distrib="MC")
integration.param <- integration_design(integcontrol=integcontrol,d=2,
                                        lower=lower,upper=upper,model=model,
                                        T=T)

# }
# NOT RUN {
obj <- max_sur(lower=lower,upper=upper,optimcontrol=optimcontrol,T=T,
                model=model,integration.param=integration.param)

obj$par;obj$value
new.model <- update_km(model=model,NewX=obj$par,NewY=testfun(obj$par),
                       CovReEstimate=TRUE)

par(mfrow=c(1,2))
print_uncertainty(model=model,T=T,type="pn",lower=lower,upper=upper,
cex.points=2.5,main="probability of excursion")

print_uncertainty(model=new.model,T=T,type="pn",lower=lower,upper=upper,
new.points=1,col.points.end="red",cex.points=2.5,
main="updated probability of excursion")
# }