conservativeEstimate: Computationally efficient conservative estimate
In dazzimonti/ConservativeEstimates: Conservative estimates for excursion sets

Description Usage Arguments Value References Examples

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.

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

`alpha`	probability of conservative estimate.
`pred`	list containing mean vector (pred$mean) and covariance matrix (pred$cov).
`design`	the discretization design for the field [to be removed?]
`Thresh`	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 Thresh 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").

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

Azzimonti, D. and Ginsbourger, D. (2016). Estimating orthant probabilities of high dimensional Gaussian vectors with an application to set estimation. Preprint at hal-01289126

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.

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 Thresh
# 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)
Thresh<-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 = Thresh)
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((Thresh-pred$mean)/pred$sd)

## Not run: CE<-conservativeEstimate(alpha = 0.95,pred = pred,design = as.matrix(newdata),
Thresh = Thresh,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)