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

Description Usage Arguments Value Author(s) Examples

View source: R/picksimulationpoints_weightedPI.R

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

1
2
3
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 [email protected]

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
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)
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 Dec. 8, 2018, 2:17 a.m.