optim_dist_measure: Choose simulation points

In pGPx: Pseudo-Realizations for Gaussian Process Excursions

Choose simulation points

Description

Selects batchsize locations where to simulate the field by minimizing the distance in measure criterion or by maximizing the integrand of the distance in measure criterion. Currently it is only a wrapper for the functions max_distance_measure and max_integrand_edm.

Usage

optim_dist_measure(
  model,
  threshold,
  lower,
  upper,
  batchsize,
  algorithm = "B",
  alpha = 0.5,
  verb = 1,
  optimcontrol = NULL,
  integration.param = NULL
)

Arguments

model	a km model
threshold	threshold value
lower	a d dimensional vector containing the lower bounds for the optimization
upper	a d dimensional vector containing the upper bounds for the optimization
batchsize	number of simulations points to find
algorithm	type of algorithm used to obtain simulation points: "A" minimize the full integral criterion; "B" maximize the integrand of the criterion.
alpha	value of Vorob'ev threshold
verb	an integer to choose the level of verbosity
optimcontrol	a list containing optional parameters for the optimization, see max_sur_parallel for more details.
integration.param	a list containing parameters for the integration of the criterion A, see max_sur_parallel for more details.

Value

A list containing

• par a matrix batchsize*d containing the optimal points

• value a vector of length batchsize with the values of the criterion at each step

References

Azzimonti D. F., Bect J., Chevalier C. and Ginsbourger D. (2016). Quantifying uncertainties on excursion sets under a Gaussian random field prior. SIAM/ASA Journal on Uncertainty Quantification, 4(1):850–874.

Azzimonti, D. (2016). Contributions to Bayesian set estimation relying on random field priors. PhD thesis, University of Bern.

Examples

### Compute optimal simulation points in a 2d example
if (!requireNamespace("DiceKriging", quietly = TRUE)) {
stop("DiceKriging needed for this example to work. Please install it.",
     call. = FALSE)
}
if (!requireNamespace("DiceDesign", quietly = TRUE)) {
stop("DiceDesign needed for this example to work. Please install it.",
     call. = FALSE)
}
# Define the function
g=function(x){
  return(-DiceKriging::branin(x))
}
d=2
# Fit OK km model
design<-DiceDesign::maximinESE_LHS(design = DiceDesign::lhsDesign(n=20,
                                                                  dimension = 2,
                                                                  seed=42)$design)$design
colnames(design)<-c("x1","x2")
observations<-apply(X = design,MARGIN = 1,FUN = g)
kmModel<-DiceKriging::km(formula = ~1,design = design,response = observations,
                         covtype = "matern3_2",control=list(trace=FALSE))

# Run optim_dist_measure, algorithm B to obtain one simulation point
# NOTE: the approximating process resulting from 1 simulation point
# is very rough and it should not be used, see below for a more principled example.
simu_pointsB <- optim_dist_measure(model=kmModel,threshold = -10,
                                  lower = c(0,0),upper = c(1,1),
                                  batchsize = 1,algorithm = "B")

## Not run: 
# Get 75 simulation points with algorithm A
batchsize <- 50
simu_pointsA <- optim_dist_measure(model=kmModel,threshold = -10,
                                  lower = c(0,0),upper = c(1,1),
                                  batchsize = batchsize,algorithm = "A")

# Get 75 simulation points with algorithm B
batchsize <- 75
simu_pointsB <- optim_dist_measure(model=kmModel,threshold = -10,
                                  lower = c(0,0),upper = c(1,1),
                                  batchsize = batchsize,algorithm = "B")
# plot the criterion value
critValA <-c(simu_pointsA$value,rep(NA,25))
par(mar = c(5,5,2,5))
plot(1:batchsize,critValA,type='l',main="Criterion value",ylab="Algorithm A",xlab="batchsize")
par(new=T)
plot(1:batchsize,simu_pointsB$value, axes=F, xlab=NA, ylab=NA,col=2,lty=2,type='l')
axis(side = 4)
mtext(side = 4, line = 3, 'Algorithm B')
legend("topright",c("Algorithm A","Algorithm B"),lty=c(1,2),col=c(1,2))
par(mar= c(5, 4, 4, 2) + 0.1)

# obtain nsims posterior realization at simu_points
nsims <- 1
nn_data<-expand.grid(seq(0,1,,50),seq(0,1,,50))
nn_data<-data.frame(nn_data)
colnames(nn_data)<-colnames(kmModel@X)
approx.simu <- simulate_and_interpolate(object=kmModel, nsim = 1, simupoints = simu_pointsA$par,
                                        interpolatepoints = as.matrix(nn_data),
                                        nugget.sim = 0, type = "UK")

## Plot the approximate process realization
image(matrix(approx.simu[1,],ncol=50),
      col=grey.colors(20))
contour(matrix(approx.simu[1,],ncol=50),
        nlevels = 20,add=TRUE)
points(simu_pointsA$par,pch=17)
points(simu_pointsB$par,pch=1,col=2)


## End(Not run)

pGPx documentation built on Aug. 23, 2023, 5:09 p.m.