ProbaMin {anMC}

Probability of exceedance of minimum of Gaussian vector

Description

Computes P(min X \le 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

ProbaMin(
  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 [0,\infty]) 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<-ProbaMin(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 ProbaMin with user defined function for active dimension probability estimate
if(!requireNamespace("TruncatedNormal", quietly = TRUE)) {
stop("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<-ProbaMin(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 ProbaMin 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<-ProbaMin(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)