computeVorobTerms {ARCHISSUR}R Documentation

Compute conditional probabilities

Description

computeVorobTerms compute the future uncertainty.

Usage

computeVorobTerms(
  i,
  object,
  intpoints.oldmean,
  krig2,
  sk.new,
  alpha,
  gpc,
  X.new,
  seed = NULL
)

Arguments

i

the index for which conditional probabilities are computed. It ranges from 1 to 2^{nrow(X.new)} which represents all possible combinations.

object

an object of class gpcm.

intpoints.oldmean

vectors of size ncol(integration.points) corresponding to the mean at the integration points before adding the batchsize points x to the design of experiments.

krig2

a list containing the output of the function predict_update_gpc_parallel.

sk.new

conditional standard deviations vector at the points intpoints.

alpha

a scalar representing the Vorob'ev threshold.

gpc

a list containing the output of the predict function predict at X.new.

X.new

input vector of size batchsize*d at which one wants to evaluate the criterion.

seed

to fix the seed.

Value

a list with:

term2

Vector equal to p_{n+q}p(y) where p_{n+q} is the updated feasibility probability at integration.points and p(y) is the probability of outcomes y at new points X.new

term1

Vector equal to the product of term2 with the indicator of p_{n+q}> \alpha where \alpha is the Vorob'ev threshold

term3

Vector equal to the product of p_{y} with the indicator of p_{n+q}< \alpha

Author(s)

Morgane MENZ, Delphine SINOQUET, Miguel MUNOZ-ZUNIGA. Contributors: Naoual SERRAJI.

References

Menz, M., Munoz-Zuniga, M., Sinoquet, D. Estimation of simulation failure set with active learning based on Gaussian Process classifiers and random set theory (2023). https://hal.science/hal-03848238.

Examples

#-------------------------------------------------------------------
#--------------------------- computeVorobTerms -----------------------------
#-------------------------------------------------------------------

## 20-points DoE, and the corresponding response
d <- 2
nb_PX <- 20
x <- matrix(c(0.205293785978832, 0.0159983370750337,
              0.684774733109666, 0.125251417595962,
              0.787208786290006, 0.700475706055049,
              0.480507717105934, 0.359730889653793,
              0.543665267336735, 0.565974761807069,
              0.303412043992361, 0.471502352650857,
              0.839505250127309, 0.504914690245002,
              0.573294917143728, 0.784444726564573,
              0.291681289223421, 0.255053812451938,
              0.87233450888786, 0.947168337730927,
              0.648262257638515, 0.973264712407035,
              0.421877310273815, 0.0686662506387988,
              0.190976166753807, 0.810964668176754,
              0.918527262507395, 0.161973686467513,
              0.0188128700859558, 0.43522031347403,
              0.99902788789426, 0.655561821513544,
              0.741113863862512, 0.321050086076934,
              0.112003007565305, 0.616551317575545,
              0.383511473487687, 0.886611679106771,
              0.0749211435982952, 0.205805968972305),
            byrow = TRUE, ncol = d)
require(DiceKriging)
fx <- apply(x, 1, branin)
f <- ifelse(fx < 14, -1, 1)
Xf <- as.matrix(x)

require(future)  # load future package for parallelization 
plan(multisession) # activate parallel calculations (with available cores automatic detection)

## gpcm object
require(GPCsign) 
model <- gpcm(f, Xf, coef.m = -1.25, coef.cov = c(1.17,0.89))

##
n.grid <- 20
x.grid <- seq(0,1,length=n.grid)
newdata <- expand.grid(x.grid,x.grid)
newdata <- as.matrix(newdata)
pred1 <- predict(object=model,newdata=newdata)
precalc.data <- list()
precalc.data$c.K <- crossprod(pred1$c, model@invK)
newdata.oldsd <- sqrt(pred1$Zsimu_var)

# new points added
new.x <- matrix(c(0.1,0.2),ncol=2,byrow=TRUE)

# predicion at new points
pred2 <- predict(object=model,newdata=new.x)
Sigma.r <- pred2$cov

newdata <- scale(x = newdata, center = model@X.mean, scale = model@X.std)
new.x <- scale(x = new.x, center = model@X.mean, scale = model@X.std)
kn <- computeQuickgpccov(object = model,
                         integration.points = newdata,
                         X.new = new.x,
                         precalc.data = precalc.data,
                         c.newdata = pred2$c)

updated.predictions <- predict_update_gpc_parallel(Sigma.r = Sigma.r,
                                                   newdata.oldsd = newdata.oldsd,
                                                   kn = kn)
# parameters for comp_term function
i <- 1
intpoints.oldmean <- pred1$Zsimu_mean
sk.new <- updated.predictions$sd
pn <- pred1$prob
alpha <- KrigInv::vorob_threshold(pn)

result <- computeVorobTerms(i = i, object = model,
                    intpoints.oldmean = intpoints.oldmean,
                    krig2 = updated.predictions,
                    sk.new = sk.new,
                    alpha = alpha,
                    gpc = pred2,
                    X.new = new.x)
plan(sequential) # deactivate parallel calculations: back to sequential mode 

[Package ARCHISSUR version 0.0.1 Index]