ProbaMax {anMC} | R Documentation |
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 |
q |
number of active dimensions, it can be either
|
pn |
coverage probability function evaluated with |
lightReturn |
boolean, if |
method |
method chosen to select the active dimensions. See |
verb |
level of verbosity (0-5), selects verbosity also for |
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
It must return a matrix |
pmvnorm_usr |
function to compute core probability on active dimensions. Inputs:
returns a the probability value with attribute "error", the absolute error. Default is the function |
Value
A list containing
probability
: The probability estimatevariance
: the variance of the probability estimateq
: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 estimatorp_q
,Rq
the conditional probabilityR_q
Eq
: the points of the designE
selected forp_q
indQ
: the indices of the active dimensions chosen forp_q
resRq
: The list returned by the MC method used forR_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)