picksimulationpoints_weightedPI: Choosing simulation points to maximize an exceedance...

In IRSN/RobustInv: "Robust inversion of expensive black-box functions"

Description

Selects (without replacement) p points among some candidates in order to maximize the exceedance probability over a threshold T of a Gaussian process simulated in these p points. The algorithm consist in choosing sequentially the points with highest individual exceedance probability and to use a penalty which depends on the distance to the points already selected.

Usage

picksimulationpoints_weightedPI(model, p, fullpoints,
  opt.simulation.points, topt.simulation.points = NULL, opt.index, T,
  type = "UK", unscale = FALSE)

Arguments

model	The current kriging model. km object.
p	Number of simulation points to select, without replacement.
fullpoints	A matrix with d columns containing all the candidate points.
opt.simulation.points	A matrix with d.opt columns, where d.opt is the number of nuisance parameters, containing the projection of fullpoints on the space Dopt of nuisance parameters.
topt.simulation.points	The transpose of opt.simulation.points. Used to improve speed.
opt.index	Array with integers corresponding to the indices of the nuisance parameters.
T	Target threshold.
type	String. The type of kriging used here. Recommended value is "UK".
unscale	Boolean. Used to apply scaling on the simulation points. This does have an influence on the computed distances between the points.

Value

An array of size p containing the indices of the selected points. If we call this output result, then the final simulation points are fullpoints[result,].

Author(s)

Clement Chevalier [e-mail]

Examples

library(KrigInv)
library(randtoolbox)
myfun <- branin_robinv
d <- 3

set.seed(8)

n0 <- 30
T <- 10
opt.index <- c(3)
inv.index <- c(1,2)
lower <- rep(0,times=d)
upper <- rep(1,times=d)
d.inv <- length(inv.index);d.opt <- length(opt.index)
lower.inv <- lower[inv.index];upper.inv <- upper[inv.index]
lower.opt <- lower[opt.index];upper.opt <- upper[opt.index]

design <- matrix(runif(d*n0),nrow=n0)
response <- myfun(design)
model <- km(formula = ~1,design = design,response = response,covtype = "matern3_2")

p <- n.optpoints <- 40

n <- n.points <- 50 # number of integration points
n.optpoints.candidates <- 500
inv.integration.points <- t(lower.inv + t(sobol(n=n.points,dim=d.inv))*(upper.inv-lower.inv))
opt.simulation.points <- t(lower.opt + t(sobol(n=n.optpoints.candidates,dim=d.opt))*(upper.opt-lower.opt))

# build the simulation points associated to the first integration point only
fullpoints <- matrix(c(0), nrow= n.optpoints.candidates,ncol=d)
fullpoints[,inv.index[1]] <-  inv.integration.points[1,1]
fullpoints[,inv.index[2]] <-  inv.integration.points[1,2]
fullpoints[,opt.index] <- opt.simulation.points

result <- picksimulationpoints_weightedPI(model = model,p = p,fullpoints = fullpoints,
                                         opt.simulation.points = opt.simulation.points,
                                         opt.index=opt.index,T = T,unscale=FALSE)


############################################################
# an example with scaling

library(KrigInv)
myfun <- function(x){ return(-1*branin_robinv(x) - 50*sin(min(100,1/x[3]))  ) }
d <- 3

set.seed(8)

n0 <- 60
T <- -60
opt.index <- c(3)
inv.index <- c(1,2)
lower <- rep(0,times=d)
upper <- rep(1,times=d)
d.inv <- length(inv.index);d.opt <- length(opt.index)
lower.inv <- lower[inv.index];upper.inv <- upper[inv.index]
lower.opt <- lower[opt.index];upper.opt <- upper[opt.index]

design <- matrix(runif(d*n0),nrow=n0)
response <- apply(X = design,FUN = myfun,MARGIN = 1) 
knots.number <- c(0,0,3)
knots <- generate_knots(knots.number = knots.number , d = d)

model <- km(formula = ~1,design = design,response = response,covtype = "matern3_2",
scaling = TRUE,knots=knots)
model@covariance@eta

p <- n.optpoints <- 40
n <- n.points <- 50 # number of integration points
n.optpoints.candidates <- 500
inv.integration.points <- t(lower.inv + t(sobol(n=n.points,dim=d.inv))*(upper.inv-lower.inv))
opt.simulation.points <- t(lower.opt + t(sobol(n=n.optpoints.candidates,dim=d.opt))*(upper.opt-lower.opt))

# more candidates in regions where the gp moves
opt.simulation.points <- unscalingFun(mat = opt.simulation.points,model = model,indices = opt.index,standardize = TRUE,lower=lower,upper=upper)

fullpoints <- matrix(c(0), nrow= n.optpoints.candidates,ncol=d)
fullpoints[,inv.index[1]] <-  inv.integration.points[1,1]
fullpoints[,inv.index[2]] <-  inv.integration.points[1,2]
fullpoints[,opt.index] <- opt.simulation.points

result <- picksimulationpoints_weightedPI(model = model,p = p,fullpoints = fullpoints,opt.index = opt.index,
                                          opt.simulation.points = opt.simulation.points,
                                          T = T,unscale=TRUE)

plot(opt.simulation.points[result,])

IRSN/RobustInv documentation built on Nov. 20, 2019, 10:46 p.m.