conservativeEstimate: Computationally efficient conservative estimate

In anMC: Compute High Dimensional Orthant Probabilities

Computationally efficient conservative estimate

Description

Computes conservative estimates with two step GANMC procedure for a Gaussian vector. The probability is approximated with a biased low dimensional estimator and the bias is corrected with a MC estimator.

Usage

conservativeEstimate(
  alpha = 0.95,
  pred,
  design,
  threshold,
  pn = NULL,
  type = ">",
  verb = 1,
  lightReturn = T,
  algo = "GANMC"
)

Arguments

alpha	probability of conservative estimate.
pred	list containing mean vector (pred$mean) and covariance matrix (pred$cov).
design	a matrix of size length(pred$mean)x(input space dimension) that contains the design where pred$mean was computed.
threshold	threshold, real number.
pn	coverage probability function, vector of the same length as pred$mean (if not specified it is computed).
type	type of excursion: ">" for excursion above threshold or "<" for below.
verb	level of verbosity, integer from 1–7.
lightReturn	boolean for light return.
algo	choice of algorithm for computing probabilities ("GANMC", "GMC").

Value

A list containing the conservative estimate (set), the Vorob'ev level (lvs). If lightReturn=FALSE, it also returns the actual probability of the set (proba) and the variance of this estimate (vars).

References

Azzimonti, D. and Ginsbourger, D. (2018). Estimating orthant probabilities of high dimensional Gaussian vectors with an application to set estimation. Journal of Computational and Graphical Statistics, 27(2), 255-267. Preprint at hal-01289126

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

Bolin, D. and Lindgren, F. (2015). Excursion and contour uncertainty regions for latent Gaussian models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 77(1):85–106.

Examples

if (!requireNamespace("DiceKriging", quietly = TRUE)) {
stop("DiceKriging needed for this example to work. Please install it.",
     call. = FALSE)
}
# Compute conservative estimate of excursion set of testfun below threshold
# Initialize
testfun<-function(x){return(((3*x^2+7*x-3)*exp(-1*(x)^2)*cos(5*pi*x^2)-1.2*x^2))}
mDet<- 1500

# Uniform design points
set.seed(42)
doe<-runif(n = 8)
res<-testfun(doe)
threshold<-0
# create km
smallKm <- DiceKriging::km(design = matrix(doe,ncol=1),
response = res,covtype = "matern5_2",coef.trend = -1,coef.cov = c(0.05),coef.var = 1.1)
# prediction at newdata
newdata<-data.frame(matrix(seq(0,1,,mDet),ncol=1)); colnames(newdata)<-colnames(smallKm@X)
pred<-DiceKriging::predict.km(object = smallKm,newdata = newdata,type = "UK",cov.compute = TRUE)

## Not run: 
# Plot (optional)
plot(seq(0,1,,mDet),pred$mean,type='l')
points(doe,res,pch=10)
abline(h = threshold)
lines(seq(0,1,,mDet),pred$mean+pred$sd,lty=2,col=1)
lines(seq(0,1,,mDet),pred$mean-pred$sd,lty=2,col=1)

## End(Not run)
# Compute the coverage function
pn<-pnorm((threshold-pred$mean)/pred$sd)

## Not run: 
pred$cov <- pred$cov + 1e-7*diag(nrow = nrow(pred$cov),ncol = ncol(pred$cov))
CE<-conservativeEstimate(alpha = 0.95,pred = pred,design = as.matrix(newdata),
threshold = threshold,type = "<",verb=1, pn=pn,algo = "ANMC")
points(newdata[CE$set,],rep(-0.1,mDet)[CE$set],col=4,pch="-",cex=2)

## End(Not run)

anMC documentation built on Aug. 22, 2023, 5:08 p.m.