computeVolumes: Compute Excursion Volume Distribution

In pGPx: Pseudo-Realizations for Gaussian Process Excursions

Compute Excursion Volume Distribution

Description

Compute the volume of excursion for each realization, includes a bias.correction for the mean. If the input is the actual GP values, compute also the random sets.

Usage

computeVolumes(
  rand.set,
  threshold,
  nsim,
  n.int.points,
  bias.corr = F,
  model = NULL,
  bias.corr.points = NULL
)

Arguments

rand.set	a matrix of size n.int.pointsxnsim containing the excursion set realizations stored as long vectors. For example the excursion set obtained from the result of simulate_and_interpolate.
threshold	threshold value
nsim	number of simulations of the excursion set
n.int.points	total length of the excursion set discretization. The size of the image is sqrt(n.int.points).
bias.corr	a boolean, if TRUE it computes the bias correction for the volume distribution.
model	the km model for computing the bias correction. If NULL the bias correction is not computed.
bias.corr.points	a matrix with d columns with the input points where to compute the bias correction. If NULL it is initialized as the first n.int.points of the Sobol' sequence.

Value

A vector of size nsim containing the excursion volumes for each realization.

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

### Simulate and interpolate for 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=50,
                                                                  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))
# Get simulation points
# Here they are not optimized, you can use optim_dist_measure to find optimized points
simu_points <- DiceDesign::maximinSA_LHS(DiceDesign::lhsDesign(n=100,
                                                               dimension = d,
                                                               seed=1)$design)$design
# obtain nsims posterior realization at simu_points
nsims <- 30
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 = nsims, simupoints = simu_points,
                                        interpolatepoints = as.matrix(nn_data),
                                        nugget.sim = 0, type = "UK")
exVol<- computeVolumes(rand.set = approx.simu,threshold = -10,
                             nsim = nsims,n.int.points = 50^2,bias.corr=TRUE, model=kmModel)

hist(exVol, main="Excursion Volume")

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