ProbaMin {anMC}  R Documentation 
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.
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
)
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 (05), 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 
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
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), 255267. Preprint at hal01289126
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.
## 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<1pANMC$probability
# Percentage error
abs(1pANMC$probability1/(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 = lowermean,u = uppermean,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<1pANMC$probability
# Percentage error
abs(1pANMC$probability1/(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<(upper1)/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<(upperlower)/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<1pANMC$probability
# Percentage error
abs(1pANMC$probability1/(d+1))/(1/(d+1))
## End(Not run)