R/computeVolumes.R

In pGPx: Pseudo-Realizations for Gaussian Process Excursions

#' @title 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.
# Input:
#' @param rand.set a matrix of size \code{n.int.points}x\code{nsim} containing the excursion set realizations stored as long vectors. For example the excursion set obtained from the result of \code{\link{simulate_and_interpolate}}.
#' @param threshold threshold value
#' @param nsim number of simulations of the excursion set
#' @param n.int.points total length of the excursion set discretization. The size of the image is \code{sqrt(n.int.points)}.
#' @param bias.corr a boolean, if \code{TRUE} it computes the bias correction for the volume distribution.
#' @param model the \link[DiceKriging]{km} model for computing the bias correction. If \code{NULL} the bias correction is not computed.
#' @param bias.corr.points a matrix with \eqn{d} columns with the input points where to compute the bias correction. If \code{NULL} it is initialized as the first \code{n.int.points} of the Sobol' sequence.
#' @return A vector of size \code{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)
#' \donttest{
#' hist(exVol, main="Excursion Volume")
#' }
#' @export
computeVolumes = function(rand.set,threshold,nsim,n.int.points,bias.corr=F,model=NULL,bias.corr.points=NULL){


  if(!is.logical(rand.set))
    rand.set <- (rand.set > threshold)#(sign(PG-threshold)+1)/2

  if(ncol(rand.set)!=nsim & nrow(rand.set)==nsim)
    rand.set<-t(rand.set)

  # for each random set: compute the volume
  # return the realizations of the volume
  volumes<-colSums(rand.set)/n.int.points

  if(bias.corr & !is.null(model)){
    if(is.null(bias.corr.points)){
      bias.corr.points<-sobol(n=n.int.points,dim = model@d)
    }
    evalOIP<-predict_nobias_km(object = model, newdata = as.data.frame(bias.corr.points),type = "UK", se.compute = TRUE, cov.compute = FALSE,low.memory = TRUE)
    pn <- pnorm((evalOIP$mean - threshold)/evalOIP$sd)
    volumes<-volumes - rep(mean(volumes)-mean(pn),nsim)
  }



  return(volumes)
}