expDistMeasure: Compute expected distance in measure of approximate excursion...

In pGPx: Pseudo-Realizations for Gaussian Process Excursions

Compute expected distance in measure of approximate excursion set

Description

Computes expected distance in measure between the excursion set of the approximated process and the true excursion set.

Usage

expDistMeasure(
  simupoints,
  model,
  threshold,
  batchsize,
  integration.param = NULL
)

Arguments

simupoints	a numeric array of size batchsize*d containing the simulation points.
model	a km model
threshold	threshold value
batchsize	number of simulations points
integration.param	a list containing parameters for the integration of the criterion A, see max_sur_parallel for more details.

Value

A positive value indicating the expected distance in measure.

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))


threshold <- -10

# Obtain simulation point sampling from maximin LHS design
batchsize <- 50
set.seed(1)
mmLHS_simu_points <-  DiceDesign::maximinSA_LHS(DiceDesign::lhsDesign(n=batchsize,
                                                                     dimension = d,
                                                                     seed=1)$design)$design


# Compute expected distance in measure for approximation obtain from random simulation points
EDM_mmLHS <- rep(NA,batchsize)
integcontrol <- list(distrib="sobol",n.points=1000)
integration.param <- KrigInv::integration_design(integcontrol,d=d,
                                        lower=c(0,0),upper=c(1,1),
                                        model=kmModel,T=threshold)
integration.param$alpha <- 0.5
for(i in seq(1,batchsize)){
EDM_mmLHS[i]<-expDistMeasure( mmLHS_simu_points[1:i,],model = kmModel,
                             threshold = threshold,batchsize = i,
                             integration.param = integration.param  )
}

plot(EDM_mmLHS,type='l',main="Expected distance in measure",xlab="batchsize")


## Not run: 
# Get optimized simulation points with algorithm B
simu_points <- optim_dist_measure(model=kmModel,threshold = threshold,
                                  lower = c(0,0),upper = c(1,1),
                                  batchsize = batchsize,algorithm = "B")
# plot the criterion value
plot(1:batchsize,simu_points$value,type='l',main="Criterion value")

# Compute expected distance in measure for approximation obtained from optimized simulation points
EDM_optB <- rep(NA,batchsize)
for(i in seq(1,batchsize)){
  EDM_optB[i]<-expDistMeasure( simu_points$par[1:i,],model = kmModel,threshold = threshold,
                                 batchsize = i,integration.param = integration.param  )
}
plot(EDM_mmLHS,type='l',main="Expected distance in measure",
     xlab="batchsize",ylab="EDM",
     ylim=range(EDM_mmLHS,EDM_optB))
lines(EDM_optB,col=2,lty=2)
legend("topright",c("Maximin LHS","B"),lty=c(1,2),col=c(1,2))

# Get optimized simulation points with algorithm A
simu_pointsA <- optim_dist_measure(model=kmModel,threshold = threshold,
                                   lower = c(0,0),upper = c(1,1),
                                   batchsize = batchsize,algorithm = "A")
# plot the criterion value
plot(1:batchsize,simu_pointsA$value,type='l',main="Criterion value")

# Compute expected distance in measure for approximation obtained from optimized simulation points
EDM_optA <- rep(NA,batchsize)
for(i in seq(1,batchsize)){
  EDM_optA[i]<-expDistMeasure( simu_pointsA$par[1:i,],model = kmModel,threshold = threshold,
                                 batchsize = i,integration.param = integration.param  )
}
plot(EDM_mmLHS,type='l',main="Expected distance in measure",
     xlab="batchsize",ylab="EDM",
     ylim=range(EDM_mmLHS,EDM_optB,EDM_optA))
lines(EDM_optB,col=2,lty=2)
lines(EDM_optA,col=3,lty=3)
legend("topright",c("Maximin LHS","A","B"),lty=c(1,3,2),col=c(1,3,2))


## End(Not run)

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