surnew_optim_robinv_fast: The surnew criterion in robust inversion
In IRSN/RobustInv: "Robust inversion of expensive black-box functions"

Description Usage Arguments Details Value Author(s) Examples

View source: R/surnew_optim_robinv_fast.R

Compute the surnew criterion at a given point x

surnew_optim_robinv_fast(x, inv.integration.points, allsimu,
  allsimucentered, allsimupoints, allprecomp, allKn.inv,
  integration.weights = NULL, model, T, new.noise.var = NULL,
  current.sur, randmatrix, penalty_visited = 0.01)

`x`	Array of size d: the point where the criterion is evaluated.
`inv.integration.points`	n x d.inv matrix containing the integration points
`allsimu`	p x (nsimu*n) matrix where p is the number of simulation points. The first nsimu columns are the nsimu simulations in p points for the first integration point, the next nsimu columns are linked to the next integration point and so on.
`allsimucentered`	p x (nsimu*n) matrix containing the later simulations centered by substracting the kriging mean of the points where the simulations are performed.
`allsimupoints`	(n*p) x d matrix containing all the points where simulations are performed. The first p rows correspond to the simulation points linked to the first integration points, and so on.
`allprecomp`	List with 3 fields obtained from a call to the function `precomputeUpdateData` on the points `allsimupoints`. Help about `precomputeUpdateData` is given in the KrigInv package.
`allKn.inv`	(n*p) x p matrix containing n different p x p matrices. Matrix number i is the inverse of the p x p non-conditional covariance matrix of the p simulation points associated to the integration point i.
`integration.weights`	Array of size n containing the weights given to each integration points.
`model`	The current kriging model. km object.
`T`	Target threshold.
`new.noise.var`	Noise variance of the new observations. Leave to NULL for noiseless functions. For noisy functions, any non zero value is valid and will give the same result.
`current.sur`	The current integral of pn (1-pn) where pn is the excursion probability. This argument is used to floor the value of the criterion since, if a new point is evaluated, the expected future uncertainty is supposed to be lower than current.sur.
`randmatrix`	n x nsimu matrix containing independent realizations of a standard gaussian random variable.
`penalty_visited`	For points which are too close to the already visited points, the function outputs `current.sur + penalty_visited` (i.e. a 'bad' value for the criterion) in order to force exploration.

The arguments inv.integration.points, allsimu, allsimucentered, allsimupoints, allprecomp, allKn.inv, integration.weights can be generated with a single call to the integration_design_robinv function.

The code of this function is meant to be as fast as possible. A large part of the computation of the criterion is performed through a call to an optimized C++ function. Suggestions to further improve the speed of these computations are of course welcome.

A scalar: the value of the surnew criterion

Clement Chevalier clement.chevalier@unine.ch

library(KrigInv)
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)

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

n <- 20 ;p <- 50; nsimu <- 1000

integcontrol <- list(distrib = "surnew",n.points = n,finaldistrib="surnew",
                     n.candidates=50,nsimu=nsimu,n.optpoints = p,
                     choose_optpoints=TRUE,n.optpoints.candidates=500)
## Not run: 
obj <- integration_design_robinv(integcontrol = integcontrol,d=d,lower=lower,upper=upper,
                                 opt.index=opt.index,inv.index=inv.index,model=model,T=T)

randmatrix <- matrix(rnorm(nsimu*n),nrow=n)
current.sur <- sum(obj$integration.weights*obj$pn*(1-obj$pn))

x <- c(0.841295 , 0.0757517 , 0.7507468)

result <- surnew_optim_robinv_fast(x = x,allsimu = obj$allsimu,
                                   inv.integration.points = obj$inv.integration.points,
                                   allsimucentered = obj$allsimucentered,
                                   allsimupoints = obj$allsimupoints,
                                   allprecomp = obj$allprecomp,
                                   allKn.inv = obj$allKn.inv,
                                   integration.weights = obj$integration.weights,
                                   model = model,T = T,new.noise.var = NULL,
                                   current.sur = current.sur,randmatrix = randmatrix)

result
current.sur

## End(Not run)