ProbaMax: Probability of exceedance of maximum of Gaussian vector

In anMC: Compute High Dimensional Orthant Probabilities

Probability of exceedance of maximum of Gaussian vector

Description

Computes P(max X > threshold) with choice of algorithm between ANMC_Gauss and MC_Gauss. By default, the computationally expensive sampling parts are computed with the Rcpp functions.

Usage

ProbaMax(
  cBdg,
  threshold,
  mu,
  Sigma,
  E = NULL,
  q = NULL,
  pn = NULL,
  lightReturn = T,
  method = 4,
  verb = 0,
  Algo = "ANMC",
  trmvrnorm = trmvrnorm_rej_cpp,
  pmvnorm_usr = pmvnorm
)

Arguments

cBdg	computational budget.
threshold	threshold.
mu	mean vector.
Sigma	covariance matrix.
E	discretization design for the field. If NULL, a simplex-lattice design n,n is used, with n=length(mu). In this case the choice of method=4,5 are not advised.
q	number of active dimensions, it can be either an integer: in this case the optimal q active dimension are chosen; a numeric vector of length 2: this is the range where to search for the best number of active dimensions; NULL: q is selected as the best number of active dimensions in the feasible range.
pn	coverage probability function evaluated with mu, Sigma. If NULL it is computed automatically.
lightReturn	boolean, if TRUE light return.
method	method chosen to select the active dimensions. See selectActiveDims for details.
verb	level of verbosity (0-5), selects verbosity also for ANMC_Gauss (verb-1) and MC_Gauss (verb-1).
Algo	choice of algorithm to compute the remainder Rq ("ANMC" or "MC").
trmvrnorm	function to generate truncated multivariate normal samples, it must have the following signature trmvrnorm(n,mu,sigma,upper,lower,verb), where n: number of simulations; mu: mean vector of the Normal variable of dimension d; sigma: covariance matrix of dimension d x d; upper: vector of upper limits of length d; lower: vector of lower limits of length d; verb: the level of verbosity 3 basic, 4 extended. It must return a matrix d x n of realizations. If not specified, the rejection sampler trmvrnorm_rej_cpp is used.
pmvnorm_usr	function to compute core probability on active dimensions. Inputs: lower: the vector of lower limits of length d. upper: the vector of upper limits of length d. mean: the mean vector of length d. sigma: the covariance matrix of dimension d. returns a the probability value with attribute "error", the absolute error. Default is the function pmvnorm from the package mvtnorm.

Value

A list containing

• probability: The probability estimate

• variance: the variance of the probability estimate

• q:the number of selected active dimensions

If lightReturn=F then the list also contains:

• aux_probabilities: a list with the probability estimates: probability the actual probability, pq the biased estimator p_q, Rq the conditional probability R_q

• Eq: the points of the design E selected for p_q

• indQ: the indices of the active dimensions chosen for p_q

• resRq: The list returned by the MC method used for R_q

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.

Chevalier, C. (2013). Fast uncertainty reduction strategies relying on Gaussian process models. PhD thesis, University of Bern.

Dickmann, F. and Schweizer, N. (2014). Faster comparison of stopping times by nested conditional Monte Carlo. arXiv preprint arXiv:1402.0243.

Genz, A. (1992). Numerical computation of multivariate normal probabilities. Journal of Computational and Graphical Statistics, 1(2):141–149.

Examples

## Not run: 
# Compute probability P(X \in (-\infty,0]) with X~N(0,Sigma)
d<-200     # example dimension
mu<-rep(0,d)    # mean of the normal vector
# correlation structure (Miwa et al. 2003, Craig 2008, Botev 2016)
Sigma<-0.5*diag(d)+ 0.5*rep(1,d)%*%t(rep(1,d))

pANMC<-ProbaMax(cBdg=20, q=min(50,d/2), E=seq(0,1,,d), threshold=0, mu=mu, Sigma=Sigma,
 pn = NULL, lightReturn = TRUE, method = 3, verb = 2, Algo = "ANMC")
proba<-1-pANMC$probability

# Percentage error
abs(1-pANMC$probability-1/(d+1))/(1/(d+1))


# Implement ProbaMax with user defined function for active dimension probability estimate
if(!requireNamespace("TruncatedNormal", quietly = TRUE)) {
stop("Package TruncatedNormal needed for this example to work. Please install it.",
     call. = FALSE)
}

# define pmvnorm_usr with the function mvNcdf from the package TruncatedNormal
pmvnorm_usr<-function(lower,upper,mean,sigma){
    pMET<-TruncatedNormal::mvNcdf(l = lower-mean,u = upper-mean,Sig = sigma,n = 5e4)
    res<-pMET$prob
    attr(res,"error")<-pMET$relErr
    return(res)
}

pANMC<-ProbaMax(cBdg=20, q=min(50,d/2), E=seq(0,1,,d), threshold=0, mu=mu, Sigma=Sigma,
 pn = NULL, lightReturn = TRUE, method = 3, verb = 2, Algo = "ANMC",pmvnorm_usr=pmvnorm_usr)
proba<-1-pANMC$probability

# Percentage error
abs(1-pANMC$probability-1/(d+1))/(1/(d+1))


# Implement ProbaMax with user defined function for truncated normal sampling

if(!requireNamespace("tmg", quietly = TRUE)) {
stop("Package tmg needed for this example to work. Please install it.",
     call. = FALSE)
}
trmvrnorm_usr<-function(n,mu,sigma,upper,lower,verb){
  M<-chol2inv(chol(sigma))
 r=as.vector(M%*%mu)

 if(all(lower==-Inf) && all(upper==Inf)){
   f<- NULL
   g<- NULL
 }else{
   if(all(lower==-Inf)){
     f<--diag(length(mu))
     g<-upper
     initial<-(upper-1)/2
   }else if(all(upper==Inf)){
     f<-diag(length(mu))
     g<- -lower
     initial<-2*(lower+1)
   }else{
     f<-rbind(-diag(length(mu)),diag(length(mu)))
     g<-c(upper,-lower)
     initial<-(upper-lower)/2
   }
 }
 reals_tmg<-tmg::rtmg(n=n,M=M,r=r,initial = initial,f=f,g=g)

 return(t(reals_tmg))
}

pANMC<-ProbaMax(cBdg=20, q=min(50,d/2), E=seq(0,1,,d), threshold=0, mu=mu, Sigma=Sigma,
 pn = NULL, lightReturn = TRUE, method = 3, verb = 2, Algo = "ANMC",trmvrnorm=trmvrnorm_usr)
proba<-1-pANMC$probability

# Percentage error
abs(1-pANMC$probability-1/(d+1))/(1/(d+1))

## End(Not run)

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