integration_design_robinv: Constructing a set of integration and simulation points in...

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

Description

Constructing a set of integration and simulation points in robust inversion

Usage

integration_design_robinv(integcontrol = NULL, d = NULL, lower, upper,
  opt.index, inv.index, model, T, min.prob = 0.001, seed = NULL)

Arguments

integcontrol	An important list with the following possible fields. distrib: This field needs to be set to either "sur" or "surnew" depending on the sampling criterion that is used. (I) If integcontrol$distrib = "surnew" (default) the user can set the field n.points (default: 20) which is the number of integration points. The points live in the space of controlled parameters, called Dinv here. The distribution of these n.points points is controlled through the finaldistrib field. Possible values are "spec" to set up manually the integration points through the inv.integration.points field, "MC" for uniformly distributed points over Dinv, "sobol" (default) to use the sobol sequence and "surnew" to use importance sampling based on a density proportional to pn (1-pn) with pn being the excursion probability of the integration points. In that last case, extra fields are required: n.candidates is the number of points which are candidate to become integration points (default: n.points*3). The important sampling scheme will choose (with replacement) n.points integration points among the candidates. The distribution of the candidates is set with the field init.distrib which can be either "MC" for uniform points, "sobol" (default) to use the sobol sequence or "spec" to set up manually the matrix of candidates with the inv.integration.points field. When scaling is used in the kriging model, the candidate points can be unscaled by setting the unscale.inv.integration.points to TRUE (default: FALSE). Regarding the nuisance parameters, the number of simulation points used, e.g., to compute pn is set with the field n.optpoints (default: 100). There is the option to choose these points among candidates by setting the field choose_optpoints to TRUE (default: FALSE). In that case, the number of candidates is set with the field n.optpoints.candidates (default: 10*n.optpoints). The field opt.simulation.points.distrib refers to the distribution of the simulation points if integcontrol$choose_optpoints=FALSE and to the candidates otherwise. Possible values are "spec" for points set up manually with the opt.simulation.points field and "sobol" (default) for the sobol sequence. When the sobol sequence is used, there is the option to add the vertices of the input domain (of the nuisance parameters) with the field addvertices (default: TRUE). When scaling is used in the kriging model, the simulation points (or the candidates) can be unscaled by setting the unscale.opt.simulation.points field to TRUE (default: TRUE). The number of simulation of Gaussian process sample paths per integration points is controlled with the field nsimu (default: 1000). Finally, when multiple calls are performed, the seed argument can be overwritten by setting the seed field. (II) The option integcontrol$distrib = "sur" needs to be used only when the sur criterion (and not surnew) is used. In that case, the fields n.points and finaldistrib work similarly except that the importance sampling option is set with integcontrol$finaldistrib="sur" and not "surnew". In that case, the field init.distrib works as in (I). Since there is no Gaussian process simulation here, the excursion probability computation is set with the field pnmethod. Possible values are "fixed" and "optimMC" (recommended). The number of MC trials to get pn (through a multivariate normal cdf call) is set with the field pop.size. Finally, the discretization in the space of nuisance parameters is also obtained with the field n.optpoints (maximum value: 10).
d	Dimension (i.e. number of scalar input variables) of the objective function.
lower	Array of size d. Lower bound of the input domain.
upper	Array of size d. Upper bound of the input domain.
opt.index	Array with integers corresponding to the indices of the nuisance parameters.
inv.index	Array with integers corresponding to the indices of the controlled parameters. We must have d = length(opt.index) + length(inv.index).
model	The current kriging model. km object.
T	Target threshold.
min.prob	Minimum probability for a candidate integration point to be picked.
seed	The random seed. For repeatability.

Details

(I) If the "surnew" criterion is used, then the output is a list with the following fields:
(i) inv.integration.points: n x dinv matrix containing the integration points. n is the number of integration points and dinv the dimension of the space of controlled parameters.
(ii) integration.weights: array of size n containing the weights given to each integration points.
(iii) pn: array of size n containing the current excursion probabilities for each integration points. These excursion probabilities are based on nsimu conditional Gaussian process simulations.
(iv) 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.
(v) allsimucentered: p x (nsimu*n) matrix containing the later simulations centered by substracting the kriging mean of the points where the simulations are performed.
(vi) 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.
(vii) allprecomp: a 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.
(viii) 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.
(ix) allmn: n*p array with the kriging means of all n*p simulation points stored in allsimupoints.
(II) If the "sur" criterion is used, then the output is still a list, which however has different fields than in (I):
(i) integration.points : (p*n) x d matrix containing all the simulation points. The first p points are associated to the first integration point and so on.
(ii) cov: see description of the field allKn.inv (viii) in (I).
(iii) mean: see description of the field allmn (ix) in (I).
(iv) pn: array of size n containing the current excursion probabilities for each integration points.

Value

Depeding on the fields in the integcontrol list, the return is different. See the description on the integcontrol argument and the details.

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)

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

integcontrol <- list(distrib = "surnew",n.points = 20,finaldistrib="surnew",
                     n.candidates=50,nsimu=1000,n.optpoints = 50,
                     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)                                  

## End(Not run)
#####################################
# 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 <- 30
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 <- 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)

integcontrol <- list(distrib = "surnew",n.points = 20,finaldistrib="surnew",
                     n.candidates=200,nsimu=100,n.optpoints = 50,
                     choose_optpoints=TRUE,n.optpoints.candidates=500,
                     unscale.opt.simulation.points=TRUE)
## 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)

## End(Not run)                                

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