Dario Azzimonti

Estimating orthant probabilities of high dimensional Gaussian vectors with an application to set estimation

Keywords: Conservative set estimation, Gaussian probabilities, Gaussian random fields, Monte Carlo

The computation of Gaussian orthant probabilities has been extensively studied for low-dimensional vectors. Here, we focus on the high-dimensional case and we present a two-step procedure relying on both deterministic and stochastic techniques. The proposed estimator relies indeed on splitting the probability into a lowdimensional term and a remainder. While the low-dimensional probability can be estimated by fast and accurate quadrature, the remainder requires Monte Carlo sampling. We further refine the estimation by using a novel asymmetric nested Monte Carlo (anMC) algorithm for the remainder and we highlight cases where this approximation brings substantial efficiency gains. The proposed methods are compared against state-of-the-art techniques in a numerical study, which also calls attention to the advantages and drawbacks of the procedure. Finally, the proposed method is applied to derive conservative estimates of excursion sets of expensive to evaluate deterministic functions under a Gaussian random field prior, without requiring a Markov assumption. Supplementary material for this article is available online.

Introduction

Assume that X = (X 1 , . . . , X d ) is a random vector with Gaussian distribution N d (µ, Σ).

We are interested in estimating, for any fixed t ∈ R, the following probability π(t) = P (X ≤ (t, . . . , t)).

(1)

The general problem of evaluating π(t), which, for a full rank matrix Σ, is the integral of the multivariate normal density φ(•; µ, Σ) over the one-sided d-dimensional rectangle (-∞, t] d , has been extensively studied in moderate dimensions with many different methods. In low dimensions tables are available (see, e.g., [START_REF] Owen | Tables for computing bivariate normal probabilities[END_REF] for d = 2). Furthermore, when the dimension is smaller than 20, there exist methods (see, e.g., [START_REF] Abrahamson | Orthant probabilities for the quadrivariate normal distribution[END_REF], [START_REF] Moran | The monte carlo evaluation of orthant probabilities for multivariate normal distributions[END_REF], [START_REF] Miwa | The evaluation of general non-centred orthant probabilities[END_REF] and [START_REF] Craig | A new reconstruction of multivariate normal orthant probabilities[END_REF]) exploiting the specific orthant structure of the probability in (1). Currently, however, most of the literature uses numerical integration techniques to approximate the quantity. In moderate dimensions fast reliable methods are established to approximate π(t) (see, e.g. [START_REF] Cox | A simple approximation for bivariate and trivariate normal integrals[END_REF]) and more recently the methods introduced in [START_REF] Schervish | Algorithm AS 195: Multivariate normal probabilities with error bound[END_REF]; [START_REF] Genz | Numerical computation of multivariate normal probabilities[END_REF] and [START_REF] Hajivassiliou | Simulation of multivariate normal rectangle probabilities and their derivatives theoretical and computational results[END_REF] (see also [START_REF] Genz | Comparison of methods for the computation of multivariate t probabilities[END_REF], [START_REF] Ridgway | Computation of gaussian orthant probabilities in high dimension[END_REF] and the book [START_REF] Genz | Computation of Multivariate Normal and t Probabilities[END_REF] for a broader overview) provide state-of-the-art algorithms when d < 100. The method introduced by Genz (1992) has been recently revised in [START_REF] Botev | The normal law under linear restrictions: simulation and estimation via minimax tilting[END_REF] where a more efficient tilted estimator is proposed. Those techniques rely on fast quasi Monte Carlo (qMC) methods and are very accurate for moderate dimensions. Here we focus on problems where d is larger than 1000 and π(t) is not a rare event probability. Such estimation problems occur, for example, if π(t) comes from a discretization of a Gaussian random field and t is a fixed finite threshold. In such cases, existing techniques are not computationally efficient or become intractable. Commonly used alternative methods are standard Monte Carlo (MC)

techniques (see [START_REF] Tong | The multivariate normal distribution[END_REF], Chapter 8 for an extensive review), for which getting accurate estimates can be computationally prohibitive.

We propose here a two step method that exploits the power of qMC quadratures and the flexibility of stochastic simulation for the specific problem of estimating π(t). We rely on the following equivalent formulation.

π(t) = 1 -P (max X > t),
where max X denotes max i=1,...,d X i . In the following we fix t and denote p = P (max X > t).

The central idea is using a moderate dimensional subvector of X to approximate p and then correcting bias by MC. Let us fix q d and define the active dimensions as E q = {i 1 , . . . , i q } ⊂ {1, . . . , d}. Let us further denote with X q the q dimensional vector X q = (X i 1 , . . . , X iq ) and with X -q the (d -q) dimensional vector X -q = (X j ) j∈E\Eq . Then, p = P (max X > t) = p q + (1 -p q )R q , (2) p q = P (max X q > t), R q = P (max X -q > t | max X q ≤ t).

The quantity p q is always smaller or equal to p as E q ⊂ {1, . . . , d}. Selecting a nondegenerate vector X q , we propose to estimate p q with the QRSVN algorithm [START_REF] Genz | mvtnorm: Multivariate Normal and t Distributions[END_REF] which is efficient as we choose a number of active dimensions q much smaller than d.

In [START_REF] Chevalier | Fast uncertainty reduction strategies relying on Gaussian process models[END_REF], Chapter 6, the similar problem of approximating the non-exceedance probability of the maximum of a Gaussian random field (GRF) ξ based on a few wellselected points is presented. Each component of X stands for the value of ξ at one point of a given discretization of the field's domain. Active dimensions (i.e. the well-selected points)

were chosen by numerically maximizing p q , and the remainder was not accounted for. Our proposed method, instead, does not require a full optimization of the active dimensions as we exploit the decomposition in (2) to correct the error introduced by p q . For this task, we propose two techniques to estimate the reminder R q : a standard MC technique and a novel asymmetric nested Monte Carlo (anMC) algorithm. The anMC technique draws samples by taking into account the computational cost, resulting in a more efficient estimator.

The anMC method presented is quite general, however its overall performance is depends on the techniques chosen to estimate p q and R q . The choices described in this paper are implemented as default in the R programming language (R Core Team, 2017) in the package anMC, however numerical experiments presented in Appendix C and in supplementary material show that, for some specific problems, alternative choices might be better suited.

In the remainder of the paper, we propose an unbiased estimator for p and we compute its variance in Section 2. In Section 3 we introduce the anMC algorithm in the more general setting of estimating expectations depending on two vectors with different simulation costs.

It is then explicitly applied to efficiently estimate R q . In Section 4 the results of two numerical studies are reported. The first one studies the efficiency of the anMC algorithm compared to standard MC. The second one is a benchmark study where the efficiency of the proposed methods is compared with a selection of state-of-the-art techniques. This study is extended to the case of small and very high probabilities in Appendix C. In Section 5, we present an implementation of this method to compute conservative estimates of excursion sets for expensive to evaluate functions under non-necessarily Markovian Gaussian random field priors. More details on the choice of active dimensions are presented in Appendix A.

All proofs are in Appendix B. Computer code for partially replicating the experiments presented here is attached in supplementary material, where we also report the results of an additional numerical experiment and a study on the computational times for the application of Section 5. The figures summarizing the benchmark results were produced with the package ggplot2 [START_REF] Wickham | ggplot2: Elegant Graphics for Data Analysis[END_REF].

2 The estimator properties 2.1 An unbiased estimator for p Equation (2) gives us a decomposition that can be exploited to obtain an unbiased estimator for p. In the following proposition we define the estimator and we compute its variance.

Proposition 1. Consider p q and R q , independent unbiased estimators of p q and R q respectively, then p = p q + (1 -p q ) R q is an unbiased estimator for p. Moreover its variance is

var( p) = (1 -R q ) 2 var( p q ) + (1 -p q ) 2 var( R q ) + var( p q ) var( R q ). ( 3
)
This property is independent from the choice of estimators p q and R q . In what follows we consider different efficient computational strategies for p q and R q .

2.2 Quasi Monte Carlo estimator for p q

The quantity p q can also be computed as

p q = 1 -P (X q ≤ t q ) ,
where t q denotes the q dimensional vector (t, . . . , t). The approximation of p q thus requires only an evaluation of the c.d.f. of X q . We denote with p q a generic estimator for p q and, since we assume that q d, we propose to estimate p q with the estimator p q G that uses the randomized quasi Monte Carlo integration method QRSVN introduced in Genz (1992), [START_REF] Hajivassiliou | Simulation of multivariate normal rectangle probabilities and their derivatives theoretical and computational results[END_REF] and refined in [START_REF] Genz | Computation of Multivariate Normal and t Probabilities[END_REF]. In particular we consider here the implementation of QRSVN with the variable reordering described in Genz and Bretz (2009, Section 4.1.3). The estimate's error is approximated with the variance of the randomized integration. The quantity p q G obtained with this procedure is an unbiased estimator of p q , see [START_REF] Genz | Computation of Multivariate Normal and t Probabilities[END_REF]. While this choice is implemented by default in the R package anMC, it is not the only possible choice. The package allows for user-defined functions to estimate p q . In supplementary material, Section C we present an numerical study where the MET method introduced in Botev ( 2017) is used in place of QRSVN.

In general, the estimator p q requires two choices: q, the number of active dimensions, and the dimensions themselves. The decomposition of p in Equation (2) leads to computational savings if we can approximate most of p with p q for a small q. On the other hand a large number of active dimensions allows to intercept most of the probability mass in p. Here we adopt a heuristic approach to select both q and E q sequentially by increasing the number of active dimensions until we meet an appropriate stopping condition. This approach, detailed in Algorithm 3, Appendix A, was chosen in the current implementation because it represents a good trade-off between speed and accuracy.

For a fixed q, the choice of E q plays an important role in the approximation of p because it determines the error p q -p, which is always negative. Selecting E q such that P (max X q > t) is numerically maximized, as in [START_REF] Chevalier | Fast uncertainty reduction strategies relying on Gaussian process models[END_REF], optimally reduces the bias of p q as an estimator for p. Here we are not interested in a fully fledged optimization of this quantity as the residual bias is removed with the subsequent estimation of R q , therefore, we exploit fast heuristics methods. The main tool used here is the excursion probability

p t (i) = P (X i > t) = Φ µ i -t Σ i,i ,
where Φ is the standard normal c.d.f. The function p t is widely used in spatial statistics (see, e.g. [START_REF] Bolin | Excursion and contour uncertainty regions for latent Gaussian models[END_REF] and Bayesian optimization (see, e.g. [START_REF] Kushner | A new method of locating the maximum point of an arbitrary multipeak curve in the presence of noise[END_REF][START_REF] Bect | Sequential design of computer experiments for the estimation of a probability of failure[END_REF]. In our setting it can be used to quickly identify the active dimensions.

In fact this function takes high values at dimensions with high probability of exceeding the threshold and thus contribute the most to p. We propose the following methods.

Method A: sample q indices with probability given by p t .

Method B: sample q indices with probability given by p t (1 -p t ).

These methods require only one evaluation of the normal c.d.f. at each element of the

vector µ i -t √ Σ i,i i=1,...,d
, and are thus very fast. Both methods were already introduced for sequential evaluations of expensive to evaluate functions, see, e.g., [START_REF] Chevalier | The KrigInv package: An efficient and user-friendly R implementation of kriging-based inversion algorithms[END_REF].

Figure 1 shows a comparison of the estimates p q obtained with different methods to select E q . We consider 30 replications of an experiment where p q G is used to approximate p. The dimension of the vector X is d = 1000, the threshold is fixed at t = 11. The vector X is obtained from a discretization of a six dimensional GRF, defined on [0, 1] 6 , over the first 1000 points of the Sobol' sequence [START_REF] Bratley | Algorithm 659: Implementing Sobol's quasirandom sequence generator[END_REF]. The GRF has a tensor product Matérn (ν = 5/2) covariance kernel. We generate a non constant mean function m by imposing the interpolation of a deterministic function at 70 points. The covariance kernel's hyperparameters are fixed as θ = [0.5, 0.5, 1, 1, 0.5, 0.5] T and σ 2 = 8, see [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF], Chapter 4, for details on the parametrization. In this example, the two methods clearly outperform a random choice of active dimensions.

Methods A and B work well for selecting active dimensions when the mean vector µ and the covariance matrix diagonal are anisotropic. In such cases both methods select dimensions that are a good trade-off between high variance and mean close to t.

The choices of q and of the active dimensions influence the behaviour of the estimator for R q . This aspect is discussed in more details in the next section.

Monte Carlo estimator for R q

Debiasing p q as an estimator of p can be done at the price of estimating

R q = P max X -q > t | max X q ≤ t .
There is no closed formula for R q , so it is approximated here via MC. Since X is Gaussian then so are X q , X -q and X -q | X q = x q , for any deterministic vector x q ∈ R q .

In order to estimate R q = P max X -q > t | X i 1 ≤ t, . . . , X iq ≤ t , we first generate n realizations x q 1 , . . . , x q n of X q such that X q ≤ t q . Second, we compute the mean and covariance matrix of X -q conditional on each realization x q l , l = 1, . . . , n with the following formulas

µ -q|x q l = µ -q + Σ -q,q (Σ q ) -1 (x q l -µ q ), Σ -q|q = Σ -q -Σ -q,q (Σ q ) -1 Σ q,-q , ( 4
)
where µ q , Σ q and µ -q , Σ -q are the mean vector and covariance matrix of X q and X -q respectively, Σ -q,q is the cross-covariance between the dimensions E \ E q and E q , Σ q,-q is the transpose of Σ -q,q . Note that the conditional covariance Σ -q|q does not depend on the realization x q l , therefore it can be computed before the sampling procedure. Given the mean and covariance matrix conditional on each sample x q l , we can easily draw a realization y -q|q l from X -q | X q = x q l . Once n couples (x q l , y -q|q l

), l = 1, . . . , n are drawn from the respective distributions, an estimator for R q is finally obtained as follows

R q MC = 1 n n l=1 1 max y -q|q l >t .
There exists many technique to draw realizations from X q conditional on X q ≤ t q . Here we use a crude multivariate rejection sampling algorithm [START_REF] Robert | Simulation of truncated normal variables[END_REF][START_REF] Horrace | Some results on the multivariate truncated normal distribution[END_REF], however this is not the only method possible. In numerical examples in Appendix C and in supplementary material we show that replacing crude rejection sampling with another sampler might be beneficial in some situations. In any case the cost of this step can be very high, in particular if we use rejection sampling then that cost is driven by the acceptance probability. The accepted samples satisfy the condition X q ≤ t q thus we have that the acceptance probability is P (X q ≤ t q ) = 1 -p q . This shows that the choice of q and of the active dimensions play an important role. If p q is much smaller than p, then the rejection sampler will have a high acceptance probability, however the overall method will be less efficient as most of the probability is in the remainder. On the other hand, if q and the active dimensions are well chosen, the value of p q could be very close to p. This will also lead to a slower rejection sampler as the acceptance probability would be small.

The second part of the procedure for R q

MC

, drawing samples from the distribution of X -q | X q = x q l , is instead less dependent on q and generally less expensive than the first step. The mean vector and covariance matrix computations requires only linear algebra operations as described in Equation ( 4) and realizations of X -q | X q = x q l can be generated by sampling from a multivariate normal distribution.

The difference in computational cost between the first step and the second step of the MC procedure can be exploited to reduce the variance at a fixed computational cost. This idea is exploited by the asymmetric nested MC procedure presented in Section 3.

We denote with p GMC the unbiased estimator of p defined as

p GMC = p q G + (1 -p q G ) R q MC ,
where GMC denotes the use of Genz's method for p q and MC for R q .

Figure 2 shows the box plots of 30 replications of an experiment where p is approximated with p GMC . The set-up is the same as in Fig. 1. The core of the probability is approximated with p q G and the active dimensions are chosen with Method 1. The residual . The remainder allows to correct the bias of p q G even with a small number of active dimensions. As comparison the results of the same experiment with a full MC estimator for p are also shown. For all experiments and for each method the number of samples was chosen in order to have approximately the same computational cost. The estimator p GMC exploits an almost exact method to estimate the largest part of the probability p, therefore the MC estimator R q MC has less variance than a full MC procedure for a fixed computational cost.

3 Estimation of the residual with asymmetric nested

Monte Carlo

In section 2, R q was estimated by R q

MC

. There exists many methods to reduce the variance of such estimators, including antithetic variables [START_REF] Hammersley | A new monte carlo technique: antithetic variates[END_REF], importance sampling [START_REF] Kahn | Random sampling (Monte Carlo) techniques in neutron attenuation problems-I[END_REF][START_REF] Kahn | Methods of reducing sample size in monte carlo computations[END_REF] or conditional Monte Carlo [START_REF] Hammersley | Conditional Monte Carlo[END_REF] among many others; see, e.g. Robert and Casella (2013, Chapter 4), for a broader overview. Here we focus on reducing the variance at a fixed computational cost, i.e. we are interested in increasing the estimator efficiency (Lemieux, 2009, Section 4.2).

We propose a so-called asymmetric nested Monte Carlo (anMC) estimator for R q that increases the efficiency with a parsimonious multiple use of conditioning data. In this section we develop some useful theoretical properties of anMC estimators.

The idea is to use an asymmetric sampling scheme that assigns computational resources by taking into account the actual cost of simulating each component. A similar asymmetric sampling scheme was introduced in the particular case of comparing the performance of stopping times for a real-valued stochastic process in discrete times in [START_REF] Dickmann | Faster comparison of stopping times by nested conditional Monte Carlo[END_REF]. Here we introduce this procedure in a general fashion and, in the next section, we detail it to R q

MC

. For two measurable spaces W, Z, consider two random elements W ∈ W and Z ∈ Z, defined on the same probability space and not independent.

We are interested in estimating the quantity

G = E [g(W, Z)] , (5)
where g : W ×Z → R is a measurable function, assumed integrable with respect to (W, Z)'s probability measure. Let us also assume that it is possible to draw realizations from the marginal distribution of W , Z and from the conditional distribution of Z | W = w i , for each w i sample of W . In the spirit of a Gibbs sampler, we can then obtain realizations (w i , z i ), i = 1, . . . , n of (W, Z) by simulating w i from the distribution of W and then z i from the conditional distribution Z | W = w i , leading to:

G = 1 n n i=1 g(w i , z i ). ( 6
)
This MC estimator can actually be seen as the result of a two step nested MC procedure where, for each realization w i , one inner sample z i is drawn from Z | W = w i . Note that the estimator R q MC used in Section 2 is a particular case of Equation ( 6) with W = X q | X q ≤ t q , Z = X -q and g(x, y) = 1 max y>t . As noted in Section 2, drawing realizations of X q | X q ≤ t q has a higher computational cost than simulating X -q because rejection sampling is required in the first case. More generally, let us denote with C W (n) the cost of n realizations of W and with C Z|W (m; w i ) the cost of drawing m conditional simulations from

Z | W = w i . If C W (1
) is much higher than C Z|W (1; w i ) then sampling several conditional realizations for a given w i might bring computational savings.

In the proposed asymmetric sampling scheme for each realization w i we sample m realizations z i,1 , . . . , z i,m from Z | W = w i . Assume that we use this sampling scheme for the couples (w i , z i,j ), i = 1, . . . , n, j = 1, . . . , m, then an estimator for G is

G = 1 nm n i=1 m j=1 g(w i , z i,j ). ( 7
)
For a fixed number of samples, the estimator G may have a higher variance than G due to the dependency between pairs sharing the same replicate of W . However, in many cases, the estimator G may be relevant to reduce the variance at a fixed computational time. In fact, let us fix the computational budget instead of the number of samples. If

C Z|W (1; w i ) < C W (1)
, then anMC may lead to an overall variance reduction thanks to an increased number of simulated pairs. In the remainder of the section, we show that, in the case of an affine cost functions, there exists an optimal number of inner simulations m such that var( G) < var( G). Assume If W = X q | X q ≤ t q , Z = X -q as in Section 2, then Z | W is Gaussian with mean and covariance matrix described in (4). In this case, the cost for sampling Z | W is affine, with α describing preliminary computations and β random number generation and algebraic operations. Denote with W 1 , . . . , W n replications of W . For each W i we consider the conditional distribution Z | W i and m replications Z 1,i , . . . , Z m,i . Under these assumption the total simulation budget is

C tot (n, m) = c 0 + n(c + α + βm).
If the total budget is fixed, C tot (n, m) = C fix ∈ R + , then the number of replications of W as a function of m is

N C fix (m) = C fix -c 0 c + α + βm .
The following proposition shows a decomposition of var( G) that is useful to find the optimal number of simulations m * under a fixed simulation budget C tot (n, m) = C fix .

Proposition 2. Consider n independent copies W 1 , . . . , W n of W and, for each W i , m copies Z i,j = Z j | W i j = 1, . . . , m, independent conditionally on W i . Then,

var( G) = 1 n var(g(W 1 , Z 1,1 )) - m -1 nm E var(g(W 1 , Z 1,1 ) | W 1 ) . ( 8
)
Corollary 1. Under the same assumptions, G has minimal variance when

m = m = (α + c)B β(A -B) ,
where

A = var(g(W 1 , Z 1,1 )) and B = E var(g(W 1 , Z 1,1 ) | W 1 ) . Moreover denote with ε = m -m , then the optimal integer is m * = m if ε < (2 m + 1) -4( m) 2 + 1 2 (9)
or m * = m otherwise.

Algorithmic considerations

In order to compute m * , we need the quantities A = var(g(W 1 , Z 1,1 )) and B = E var(g(W 1 , Z 1,1 ) | W 1 ) and the constants c 0 , c, α and β. A and B depend on the specific problem at hand and are usually not known in advance. Part of the total computational budget is then needed to estimate A and B. This preliminary phase is also used to estimate the system dependent constants c and β. Algorithm 1 reports the pseudo-code for anMC.

Estimate p with p GanMC

The anMC algorithm can be used to reduce the variance compared to R q 's MC estimate proposed in Section 2.3. In fact, let us consider W = X q | X q ≤ t q and Z = X -q . We have that W is expensive to simulate as it requires sampling from a truncated normal Algorithm 1: Asymmetric nested Monte Carlo.

Input : µ W , µ Z , Σ W , Σ Z , Σ W Z , g, C tot
Output: G Part 0: estimate c 0 , c, β, α ;

initialize compute the conditional covariance Σ Z|W and initialize n 0 , m 0 ;

Part 1: for i ← 1 to n 0 do estimate A, B simulate w i from the distribution of W and compute µ Z|W =w i ; draw m 0 simulations z i,1 , . . . , z i,m 0 from the conditional distribution

Z | W = w i ; estimate E [g(W, Z) | W = w i ] with Ẽi = 1 m 0 m 0 j=1 g(w i , z i,j ); estimate var (g(W, Z) | W = w i ) with Ṽi = 1 m 0 -1 m 0 j=1 (g(w i , z i,j ) -Ẽi ) 2 ; end compute m = (α+c) 1 n 0 n 0 i=1 Ṽi β 1 n 0 -1 n 0 i=1 ( Ẽi -1 n 0 n 0 i=1 Ẽi ) 2 , m * as in Corollary 1 and n * = N C fix (m * ); Part 2: for i ← 1 to n * do compute G if i ≤ n 0 then for j ← 1 to m * do if j ≤ m 0 then use previously calculated E i and V i ; else simulate z i,j from the distribution Z | W = w i ; compute E i = 1 m * m *
j=1 g(w i , z i,j ); end end else simulate w i from the distribution of W and compute µ Z|W =w i ; We estimate R q via

for j ← 1 to m * do simulate z i,j from the conditional distribution Z | W = w i ; end compute E i = 1 m * m * j=1 g(w i , z i,j ); end end estimate E [g(W, Z)] with G = 1 n * n * i=1 Ẽi ;
R q anMC = 1 n * m * n * i=1 m * j=1 1 max z i,j >t .
Finally plugging in R q anMC and p q G in Equation ( 2), we obtain

p GanMC = p q G + (1 -p q G ) R q anMC .
Figure 3a shows a comparison of results using 30 replications of the experiment presented in Section 2.3. Results obtained with a MC estimator are shown for comparison.

While the simulations of all experiments were obtained under the constraint of a fixed computational cost, the actual time to obtain the simulations was not exactly the same. In order to be able compare the methods in more general settings we further rely on the notion of efficiency. For an estimator p, we define the efficiency (Lemieux, 2009, Section 4.2) as

Eff[ p] = 1 var( p) time[ p] , (10)
where time[ p] denotes the computational time of the estimator p.

Figure 3b shows a comparison of the efficiency of p GMC and p GanMC with a full Monte Carlo estimator. With as few as q = 50 active dimensions we obtain an increase in efficiency of around 10 times on average over the 30 replications of the experiment with the estimator p GMC . The estimator p GanMC shows a higher median efficiency than the others for all q ≥ 20.

4 Numerical studies

Choice of the number of inner samples

In this section we study the efficiency of the anMC method compared with a standard MC method for different choices of m. Here we do not select the optimal m * defined in

Corollary 1, but we study the efficiency as a function of m. In many practical situations even if part 1 of Algorithm 1 does not render the optimal m * the anMC algorithm is still more efficient than a standard MC if the chosen m is close to m * .

We consider a similar setup to the experiment presented in Section 2.2. Here we start from a GRF with tensor product Matérn (ν = 5/2) and a non constant mean function m different from the example in Section 2.2, initialized as conditional mean on 60 randomly generated values at a fixed design on [0, 1] 6 . The hyperparameters are fixed as θ = [0.5, 0.5, 1, 1, 0.5, 0.5] T and σ 2 = 8. The GRF is then discretized over the first d = 1000

points of the Sobol sequence to obtain the vector X. We are interested in 1-p = P (X < t), with t = 5. We proceed by estimating p with p GMC and p GanMC for different choices of m to compare their efficiency. The initial part p q is computed once with estimator p q G with q and the active dimensions chosen with Algorithm 3, Method B. The number of outer simulations in the anMC algorithm is kept fixed to n = 10, 000 and we only vary m. For each m, the anMC estimation is replicated 20 times.

The median estimated value for p is p = 0.9644. Most of the probability is estimated with p q , in fact p q G = 0.9636. Figure 4a shows Eff[ p] computed with the overall variance of p. A choice of m = 10 leads to a median increase in efficiency of 73% compared to the MC case. In this example, both the probability to be estimated and p q are close to 1, thus the acceptance probability for R q is low. In this situation the anMC method is able to exploit the difference in computational costs to provide a more efficient estimator for R q .

In order to study the effect of the acceptance probability on the method's efficiency we change the threshold in the previous example to t = 7.5 by keeping the remaining parameters fixed. The value of p is smaller, p = 0.1178. The number of active dimensions q, chosen with Algorithm 3, is smaller (q = 90) as the probability mass is smaller. The value of p q ( p q G = 0.1172) is much smaller than in the previous case and this leads to a higher acceptance probability for R q . Figure 4b shows efficiency of the method as a function of m. Here the anMC method does not bring significant gains over the MC method as the the ratio between the cost of rejection sampling and the conditional simulations in R q is close to one. The estimated m * is equal to 1.91, thus it is smaller than the minimum threshold of Proposition 3 that guarantees a more efficient anMC algorithm. Appendix C shows a comparison of the efficiency for t = 7.5 and for t = 3. By changing the level t we obtain radically different situations: in the first case the acceptance probability of the rejection sampler becomes quite high, thus limiting the benefit of the anMC procedure. In the second case π(t) becomes very small as the dimension increases, thus making the problem of estimating R q out of reach with rejection sampling. In this example we present an alternative to the default choice implemented in anMC.

Application: efficient computation of conservative estimates

A problem where the anMC method leads to substantial increases in efficiency is conservative excursion set estimation relying on Gaussian field models. We consider an expensive to evaluate system described by a continuous function f : D ⊂ R → R, ≥ 1, where D is a compact domain, and we focus on estimating, for some fixed threshold t ∈ R, the set

Γ * = {x ∈ D : f (x) ≤ t}.
Such problems arise in many applications such as reliability engineering (see, e.g., [START_REF] Picheny | Quantile-based optimization of noisy computer experiments with tunable precision[END_REF], Chevalier et al. (2014a)) climatological studies [START_REF] Bolin | Excursion and contour uncertainty regions for latent Gaussian models[END_REF][START_REF] French | Spatio-temporal exceedance locations and confidence regions[END_REF] or in natural sciences [START_REF] Bayarri | Using statistical and computer models to quantify volcanic hazards[END_REF]. Often f is seen as expensive to evaluate black-box [START_REF] Sacks | Design and analysis of computer experiments[END_REF] and can only be evaluated with computer simulations. We assume here that f was only evaluated at points

χ k = {x 1 , . . . , x k } ⊂ D
and the associated responses are denoted with f (χ k ) = (f (x 1 ), . . . , f (x k )) ∈ R k and we are interested in giving an estimate of Γ * starting from these k evaluations.

In a Bayesian framework we consider f as a realization of a GRF (ξ x ) x∈D with prior mean function m and covariance kernel K. A prior distribution of the excursion set is hence obtained by thresholding ξ, thus obtaining the following random closed set

Γ = {x ∈ D : ξ x ≤ t}.
Denoting with ξ χ k the random vector (ξ x 1 , . . . , ξ x k ), we can then condition ξ on the observations f (χ k ) and obtain a posterior distribution for the field ξ x | ξ χ k = f (χ k ). This gives rise to a posterior distribution for Γ. Different definitions of random closed set expectation [START_REF] Molchanov | Theory of Random Sets[END_REF], Chapter 2) can be used to summarize this posterior distribution and to provide estimates for Γ * . In [START_REF] Chevalier | Estimating and quantifying uncertainties on level sets using the Vorob'ev expectation and deviation with Gaussian process models[END_REF], for example, the Vorob'ev expectation was introduced in this setting. Let us briefly recall this definition. We denote with

p Γ,k : D → [0, 1] the coverage function of the posterior set Γ | ξ χ k = f (χ k ), defined as p Γ,k (x) = P k (x ∈ Γ), x ∈ D, where P k (•) = P (• | ξ χ k = f (χ k ))
. This function associates to each point in D its probability of being inside the posterior excursion set. The function p Γ,k gives rise to a family of excursion set estimates: for each ρ ∈ [0, 1] we can define the posterior ρ-level Vorob'ev quantile of Γ

Q ρ = {x ∈ D : p Γ,k (x) ≥ ρ}.
The Vorob'ev expectation of Γ [START_REF] Molchanov | Theory of Random Sets[END_REF] is the quantile

Q ρ V that satisfies |Q ρ | ≤ E k [|Γ|] ≤ |Q ρ V | for all ρ ≥ ρ V ,
where |A| denotes the volume of a set A ⊂ R l . The set Q ρ V consists of the points that have high enough marginal probability of being inside the excursion set. In some applications, however, it is important to provide confidence statements on the whole set estimate. Conservative estimates introduced in [START_REF] Bolin | Excursion and contour uncertainty regions for latent Gaussian models[END_REF] for Gaussian Markov random fields address this issue. A conservative estimate of Γ * is

C Γ,k = arg max C⊂D {|C| : P k (C ⊂ {ξ x ≤ t}) ≥ α}, (11)
where |C| denotes the volume of C.

The estimation of the object in Equation ( 11), however, leads to major computational issues. First of all we need to select a family of sets to use for the optimization procedure in Equation ( 11). Here we follow [START_REF] Bolin | Excursion and contour uncertainty regions for latent Gaussian models[END_REF] and select the Vorob'ev quantiles as family of sets. This family has the advantage that it is parametrized by one real number ρ and thus it renders the optimization straightforward. Algorithm 2 details the optimization procedure.

Second, for each candidate Q we need to evaluate

P next = P k (Q ⊂ {ξ x ≤ t}), the
probability that Q is inside the excursion. In fact, this quantity is a high dimensional orthant probability. For a Vorob'ev quantile Q ρ , discretized over the points c 1 , . . . , c r ,

P k (Q ρ ⊂ {ξ x ≤ t}) = P k (ξ c 1 ≤ t, . . . , ξ cr ≤ t) = 1 -P k ( max i=1,...,r ξ c i > t).
Thus we use the estimator p GanMC to approximate 1-P k (Q ρ ⊂ {ξ x ≤ t}). The use of anMC allows resolutions for the discretized Vorob'ev quantiles that seem out of reach otherwise.

We apply Algorithm 2 to a two dimensional artificial test case. We consider as function standard Monte Carlo and asymmetric nested Monte Carlo (anMC). Both methods showed higher efficiency than other state-of-the-art methods for dimensions higher than 1000 when the orthant probability π(t) is not a rare event in high dimensions.

The version of the anMC method proposed here relies on the QRSVN algorithm to estimate p q and a rejection sampling to the estimate of R q . These choices, implemented as default in the package anMC, can be easily changed to improve the method. In fact those particular choices do not prove to be the most efficient when π(t) is very small in high dimensions, as shown in Appendix C. For example, in supplementary material, we

show that it is possible to implement the anMC method using the MET algorithm for the estimation of p q and a Hamiltonian Monte Carlo technique as truncated normal sampler.

Such choices provide efficiencies close to the best state-of-the-art method in all cases and allow for efficient estimates in higher dimensions.

Within its computational limits, the efficiency of p GanMC with default choices is mainly driven by the acceptance probability of the rejection sampler in R q anMC , which in turn depends on p q . This highlights the existence of a trade-off between p q G and R q . If the choice of q and active dimensions is not optimal, then the acceptance probability of the rejection sampler becomes larger, making the estimation of R q easier. An estimator p q closer to p makes the quantity R q harder to estimate, however, in this case, R q anMC becomes more efficient than R q

MC

as the ratio between computational costs becomes more favourable.

In general, anMC relies on an initial step where several constants and probabilistic quantities are empirically estimated to choose the optimal m * , the number of inner samples.

In particular the cost parameters c, β, the slopes of the linear costs, might be hard to estimate if the constants c 0 , α are comparatively large. In this case Algorithm 1 might not choose the optimal m * . However, a numerical study of the algorithm behaviour for different choices of m showed that, on the considered examples, even if the chosen m is not optimal but it is close to optimal, the efficiency gain is very close to the optimal efficiency gain.

The estimator p GanMC made the computation of conservative estimates possible for excursion sets of expensive to evaluate functions under general GRF priors. The R implementation of the algorithm is included in the package anMC currently available on CRAN and on GitHub.

q by sequentially increasing the number of active dimensions until the relative change of p q G is less than the estimate's error.

Algorithm 3: Select q, active dimensions and compute p q G .

Input : q 0 , small initial q, e.g. q 0 = d 1/3 , and q step the increment of q, γ > 0

Output: q, p G q Compute p G q 0 and save err( p G q 0 ) := 3 var( p G q 0 ) ; initialize k = 0;

repeat increment k = k + 1 ; q k := q 0 + kq step ; choose q k active dimensions, compute p G q k and err( p G q k ) ; compute ∆( p G q k ) = p G q k -p G q k-1 1+ p G q k ; until ∆( p G q k ) < γ err( p G q k ) or q k > 300; q = q k and p q G = p G q k ;
The constant γ > 0 is chosen equal to 1 in our implementation. Moreover the algorithm stops if q k > 300 to avoid using Genz's algorithm in high dimensions.

A.1 Add spatial information

If the random vector X comes from a GRF discretized over a set of points E spat = {e 1 , . . . , e d } ⊂ R l , then we can exploit this information to choose E q . Let us consider the sequence of vectors (δ j ) j=1,...,q , defined for each j as

δ j = j k=1 dist(e i k , E spat ) (j = 1, . . . , q)
where dist(e i k , E spat ) denotes the d-dimensional vector of Euclidean distances between e i k and each point in E spat and {e i 1 , . . . , e iq } are the points corresponding to the selected active dimensions E q . We then adjust Methods A, B by sampling the jth active dimension with probabilities given by the component-wise products p t δ j δ j and p t (1 -p t ) δ j δ j respectively.

B Proofs

Proof of Proposition 1

Proof. We have that E[ p q ] = p q and E[ R q ] = R q . Then we have

var( p) = var( p q ) + var((1 -p q ) R q ) = +2 cov( p q , (1 -p q ) R q ) = . (12)
We can write the variance and the covariance as

= var((1 -p q ) R q ) = (1 -p q ) 2 var( R q ) + R 2 q var( p q ) + var( p q ) var( R q ), = cov[ p q , (1 -p q ) R q ] = -var( p q )R q ,
respectively, by exploiting the independence of p q and R q . By plugging in those expressions in Equation ( 12) we obtain the result in Equation (3).

Proof of Proposition 2

Proof.

var( G) = 1 n 2 m 2 var n i=1 m j=1 g(W i , Z i,j ) = 1 nm 2 var m j=1 g(W 1 , Z 1,j ) = 1 nm 2 m j=1 m j =1 cov g(W 1 , Z 1,j ), g(W 1 , Z 1,j ) = 1 nm 2 m var(g(W 1 , Z 1,1 )) + m(m -1) cov(g(W 1 , Z 1,1 ), g(W 1 , Z 1,2 )) = 1 nm 2 [m var(g(W 1 , Z 1,1 )) + m(m -1) ] . (13)
where the first equality is a consequence of the independence of W 1 , . . . , W n and the third equality is a consequence of the independence of Z i,j and Z i,j conditionally on W i . Moreover the covariance denoted by in ( 13) can be written as follows.

= E cov(g(W 1 , Z 1,1 ), g(W 1 , Z 1,2 ) | W 1 ) =0 Z 1,1 ,Z 1,2 independent conditionally on W 1 + cov E[g(W 1 , Z 1,1 ) | W 1 ], E[g(W 1 , Z 1,2 ) | W 1 ] =var(E[g(W 1 ,Z 1,1 )|W 1 ]) = var E[g(W 1 , Z 1,1 ) | W 1 ] = var g(W 1 , Z 1,1 ) -E var g(W 1 , Z 1,1 ) | W 1 . (14)
Equations ( 13) and ( 14) give the result (8).

Proof of Corollary 1

Proof. Denote with e

= β(A -B), f = (α + c)(A -B) + βB, g = (c + α)B, h = C tot -c 0 , then var( G)(m) = em 2 + f m + g hm . (15)
Observe that the first and second derivatives of var( G) with respect to m are respectively changing the threshold t to t = 7.5 for the high probability values and t = 3 for the small values.

∂ var( G) ∂m = 1 h e - g m 2 , ∂ 2 var( G) ∂m 2 = 2g hm
In Figure 7 we show an efficiency comparison for estimating π(t), with t = 3 and t = 7.5, with the algorithms GanMC, GMC, GHK, MET and QRSVN, see Section 4.2 for details.

The threshold t = 7.5 leads to a high probability setup as the median estimate for π(t) ranges between 0.99685, for d = 100, and 0.67436, for d = 7000. Figure 7a compares the estimators' efficiency. In this case the QRSVN is the most efficient algorithm in low dimensions. The GMC and the GanMC algorithms however are the most efficient for all dimensions higher than 2000. The GanMC algorithm is 3 times more efficient than the MET for d = 3000 and 9 times more efficient than GHK for d = 5000. In this setup the computational cost of the rejection sampler in R q anMC is not much higher than the conditional sampler. In fact, the acceptance probability of the rejection sampler is always higher than 0.6. In most replications this leads to a choice of m * very small or even equal to 1. Thus GanMC is slower than GMC because of Part 1 in Algorithm 1 while achieving the same variance. This is the main reason why the GMC algorithm proves to be more efficient in dimensions 2000, 3000, 5000, 6000, in fact for d = 5000 GMC is 3.6 times more efficient than GanMC and for d = 6000 the ratio is 1.9. Computations for this experiment were carried on the cluster of the University of Bern on machines with Intel Xeon CPU 2.40GHz and 16 GB RAM.

For both the case t = 7.5 and t = 5 the probability π(t) does not qualify as small.

Figure 7b shows an efficiency comparison for t = 3, where the value of π(t) is equal to 9.02 × 10 -3 for d = 100, decreases to 1.33 × 10 -13 for d = 2000 and becomes too small to be estimated reliably for larger dimensions. The GHK method while performing well in low dimensions is not able to estimate reliably the probability for d = 1000 or higher. Notice that the GanMC method with the default rejection sampler for the residual part does not work when d > 500 because the acceptance probability is too low. If the rejection sampler is replaced by the Hamiltonian Monte Carlo method for truncated normal vectors described in [START_REF] Pakman | Exact hamiltonian monte carlo for truncated multivariate gaussians[END_REF] and implemented in the R package tmg, then the GanMC method performs better than GHK and MET for dimensions lower that 2000, however it is still 3 times less efficient than MET for d = 2000. Higher dimensional comparison were not

C

  W (n) = c 0 + cn and, for each sample w i C Z|W (m; w i ) = C Z|W (m) = α + βm, with c 0 , c, α, β ∈ R + dependent on the simulators of W and Z | W . The second equation entails that the cost of conditional simulations does not depend on the conditioning value.

Proposition 3 .

 3 Under the same assumptions, if m * > 2(α+c)B (c+α)B+β(A-B) then var( G) = var( G) [1 -η], where η ∈ (0, 1).

  Efficiency. Values in logarithmic scale.

Figure 3 :

 3 Figure 3: Comparison of results with p q G , p GMC , p GanMC and standard MC on 30 replications of the example introduced in Fig. 1.

  High probability state, t = 5, p GanMC = Low probability state, t = 7.5, p GanMC = 0.1178.

Figure 4 :

 4 Figure 4: Efficiency of p GanMC estimator versus the number of inner simulations m. For each m the experiment is reproduced 30 times.

Figure 5 :

 5 Figure 5: Efficiency of the probability estimator for π(t), with t = 5, versus the dimension d. For each d the experiment is reproduced 15 times. Values in logarithmic scale. The median estimated value for p = 1 -π(t) ranges from 0.33845 to 0.99876.

f

  a realization of a GRF (ξ x ) x∈D , where D ⊂ R 2 is the unit square. We consider two parametrizations for the prior covariance kernel: a tensor product Matérn covariance kernel with ν = 5/2, variance σ 2 = 0.5 and range parameters θ = [0.4, 0.2] and a Gaussian covariance kernel with variance σ 2 = 0.5 and range parameters θ = [0.2, 0.4]. In both cases we assume a prior constant mean function. We are interested in the set Γ * with t = 1. For both cases we consider k = 15 evaluations of f at the same points chosen by Latin hypercube sampling. Figures6a and 6bshow the conservative estimate at level 95% compared with the true excursion, the Vorob'ev expectation and the 0.95-quantile for the Matérn and the Gaussian kernel. The 0.95-quantile does not guarantee that the estimate is included in the true excursion with probability 0.95 in both examples. The conservative sition of the covariance matrix. The remainder R q is instead estimated with two methods:

  3 . The second derivative is positive for all m > 0 then var( G) is a convex function for m > 0 and the point of minimum is equal to the zero of ∂ var( G)/∂m, which is m = g/e = m.Since var( G) is convex in m, the integer that realizes the minimal variance is either m or m . By plugging in m = m -ε = g/e -ε and m = m -ε + 1 = g/e -ε + 1 in Equation (15), we obtain the condition in (9).Proof of Proposition 3Proof. First of all notice that the total cost of sampling G isC tot = c 0 + n(c + C Z|W ) = c 0 n(c + α + β). By isolating n in the previous equation we obtain n = Ctot c+α+β , where C tot for the sake of brevity and, by computations similar to those in Proposition 2 we obtainvar( G) = c + α + β C tot var(g(W 1 , Z 1,1 )) = c + α + β C tot A,where A = var(g(W 1 , Z 1,1 )). In the following we will also denote B = E var(g(W 1 , Z 1,1 ) | W 1 ) as in Corollary 1. Let us now substitute N C fix (m * ) in equation (8), thus obtainingvar( G) = (c + α + βm * )Am * -(m * -1)(c + α + βm * )B C tot m * = var( G) (m * ) 2 β(A -B) + m * [(c + α)(A -B) + βB] + (c + α)B A(c + α + β)m * = var( G) 2(α + c)B + m * [(c + α)(A -B) + βB] A(c + α + β)m * , (16)where in (16) we substituted (m * ) 2 from Corollary 1. By rearranging the terms, we obtainvar( G) = var( G) 1 -(m * -2)(c + α)B + m * β(B -A) A(c + α + β)m * = var( G) [1 -η] .Since A -B, B, c, β, α are always positive, then η < 1 for all m * > 0. Moreover η > 0 if m * > 2(α + c)B (α + c)B + β(A -B) .

The first author gratefully acknowledges the support of the Swiss National Science Foundation, grant number 146354 and of the Hasler Foundation, grant number 16065.

Comparison with state of the art

In this section we compare the GMC and GanMC methods, with the default options as implemented in the R package anMC, with available state-of-the-art algorithms to estimate π(t). In particular, we compare this implementation with:

QRSVN an implementation of Genz method [START_REF] Genz | Computation of Multivariate Normal and t Probabilities[END_REF] in the R package mvtnorm [START_REF] Genz | mvtnorm: Multivariate Normal and t Distributions[END_REF], function pmvnorm;

GHK an implementation of GHK method [START_REF] Geweke | Efficient simulation from the multivariate normal and student-t distributions subject to linear constraints and the evaluation of constraint probabilities[END_REF][START_REF] Hajivassiliou | The method of simulated scores for the estimation of LDV models[END_REF] in the R package bayesm [START_REF] Rossi | bayesm: Bayesian Inference for Marketing/Micro-Econometrics. R package version 3[END_REF], function ghkvec;

MET R implementation of the minimax-exponentially-tilted (MET) method [START_REF] Botev | The normal law under linear restrictions: simulation and estimation via minimax tilting[END_REF] in the package TruncatedNormal [START_REF] Botev | The normal law under linear restrictions: simulation and estimation via minimax tilting[END_REF], function mvNcdf;.

We consider the example introduced in Section 4.1 and we increase the dimension of the problem d by considering finer discretizations of the underlying GRF. For example, the vector X of dimension d = 100 is obtained from the GRF discretized on the first 100 points of the 6-dimensional Sobol' sequence. As the dimension d increases the probability π(t) changes, thus providing different setups. Each experiment is replicated 15 times.

Figure 5 presents a comparison of the estimator's efficiency for the problem of computing π(t), with t = 5. In this setup the value of π(t) varies between 0.66155 for d = 100 and 0.00124 for d = 7000. The most efficient algorithm is the QRSVN Genz method, however this implementation does not scale to dimensions higher than 1000. The GMC algorithm is the second most efficient in all dimensions except d = 2000 where it is the most efficient. The GanMC algorithm is instead the most efficient when d is greater than 2000. This effect is explained by the efficiency gains brought by R q anMC when the rejection sampler is expensive. If d > 2000, the probability P (X q ≤ t q ) is always smaller than 0.01, thus the rejection sampler becomes much more expensive than the conditional sampler in the estimation of the remainder R q . Algorithms GHK and MET allowed estimates until dimension d = 5000 and d = 4000 respectively before running in memory overflows. The GanMC algorithm is 45 times more efficient than the GHK algorithm for d = 5000 and 3.8 times more efficient than MET for d = 4000. It is also 5 times more efficient than GMC for Algorithm 2: Conservative estimates algorithm.

Input : m k , K k , conditional mean and covariance of ξ | ξ χ k = f (χ k ), and G, fine discretization design; Output: Conservative estimate for Γ * at level α.

Part 0: sort the points in G in decreasing order of p Γ,k , with indices

end end estimates instead are guaranteed to be inside the true excursion with probability α = 0.95. They correspond to Vorob'ev quantiles at levels 0.998 (Matérn) and 0.993 (Gaussian). The conservative estimates were obtained with a 100×100 discretization of the unit square. Such high resolution grids lead to very high dimensional probability calculations. In fact, the dichotomy algorithm required 11 computations of the probability 1

for each case. The discretization's size for Q ρ varied between 1213 and 3201 points in the Matérn kernel case and between 1692 and 2462 points in the Gaussian case. Such high dimensional probabilities cannot be computed with the current implementation of the algorithm by Genz, however they could be computed with other Monte Carlo methods at higher computational costs. Instead, with the proposed method, the total computational time on a laptop with Intel Core i7 1.7GHz CPU and 8GB of RAM was equal to 365 and 390 seconds respectively for Matérn and Gaussian kernel. In Supplementary Materials, Section D, we compare the time required for conservative estimates when the core orthant probability estimate is computed with the GanMC, MET, GHK or QRSVN algorithms.

Discussion

In this paper we introduced a new method to approximate high dimensional orthant Gaussian probabilities based on a decomposion of the probability in a low dimensional part p q and a remainder R q . The number of active dimensions q and the dimensions themselves are chosen with two heuristic algorithms which provide good results in case of dense covariance matrix with anisotropic diagonal and anisotropic mean vector. An alternative proposal is choosing the first q dimensions ordered according to the inverse Genz variable reordering proposed in Genz and Bretz (2009, Section 4.1.3). While similar to the heuristics proposed here, this method is not efficient in high dimensions as it requires a full Cholesky decompo-

SUPPLEMENTARY MATERIAL

SupplementaryMaterial: short description of the supplementary material provided, additional numerical results and computational times for Section 5. (PDF file) R-package anMC: R-package implementing the GanMC, GMC procedures and the con-

A Choice of active dimensions

The estimator p q G , introduced in Section 2.2, requires the choice of q, the number of active dimensions and the choice of the dimensions themselves. Algorithm 3 describes the heuristic procedure implemented in anMC to select q and obtain the active dimensions. Here we select C Numerical study for small and large probabilities

The algorithm GanMC was originally developed to estimate probabilities in the form π(t)

for high dimensional problems where this probability does not decrease rapidly with the dimension, as, for example, in the numerical study shown in Section 4. In this section we study two limit cases where this assumption is challenged: first we consider very high probability values for π(t), second we consider the very small probability case. For very small probabilities, the anMC method with the default choices implemented in anMC might not be the correct choice. In fact, if π(t) becomes too small then p = 1 -π(t) is very close to 1 and p q is very close to p. As suggested in Section 2.3, the residual term R q becomes increasingly hard to estimate with rejection sampling as p q becomes small. An alternative sampler for truncated normal vectors might improve performance. This can be achieved in the function ProbaMax of the package anMC by choosing a user defined truncated normal sampler.

We construct two benchmark studies with the problem defined in Section 4.1 and by possible as MET and GHK method required more memory than allowed. Computations for the experiment t = 3 were carried out on the Idiap Research Institute computing grid, on machines with Intel Xeon E312xx (Sandy Bridge) CPU 3.00GHz and 8 GB RAM.


Baptiste Magnier
email: baptiste.magnier@mines-ales.fr

Alexandre Le
email: alexandre.le@mines-ales.org

Albert Zogo
email: albert.zogo@mines-ales.org

A Quantitative Error Measure for the Evaluation of Roof Edge Detectors

Crest line extraction remains a hard task in image processing. Indeed, these roof edges represent narrow edges on the image surface and whatever undesirable pixel close to or on the crest line may disturb the detection. This communication presents a new crest line detection overall evaluation. Comparing the ground truth contour image and the candidate crest line image, the proposed algorithm is based upon a new criterion that take into account the list of ground truth, the recall and their associated spacial nearness. Doubtlessly, an efficient evaluation penalizes a misplaced edge point proportionally to the distance to the true contour.

This quantitative performance evaluation proves its efficiency on several crest line images of different types, bringing a favorable indicator for tow closest edge images or a poor in the presence of a degraded/distorted candidate edge image.

I. INTRODUCTION

In computer science, all the systems, especially automated information processing structures, must be evaluated before being developed in the industry or the medical data for principal examples. Image processing is no exception to this rule. In image segmentation, all these methods have to be tested and their quality must be measured, whether it is edge detection, point matching, region segmentation or image restoration (using the well known measures: PSNR or Similarity [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF]).

In image processing applications, edge detection remains a key point as it can provide geometrical information [START_REF] Torre | On edge detection[END_REF] [START_REF] Ziou | Edge detection techniques: an overview[END_REF]. Moreover, effective and precise edge detection [START_REF] Canny | A computational approach to edge detection[END_REF] can prove highly useful. In most cases, due to the undesirable noise, textures, some other disturbances or parameters of the edge extractor themselves, most of edge detectors do not lead directly to object edges; then contours must be searched among numerous contour points. Furthermore, as shown in [START_REF] Montesinos | A new perceptual edge detector in color images[END_REF], crossing edges and corners which are cornerstone in image processing are not well detected with most of edge extractors. Nevertheless, evaluation of all image segmentation method type is an important labor. To perform this task, error measures must be quantified between the ground truth solution and a candidate edge image. These measures could be either quantitative or qualitative.

Qualitative assessments require interaction with human observers, and the performance indicators are based on subjective evaluation of the quality of the detected contours [START_REF] Martin | Learning to detect natural image boundaries using local brightness, color, and texture cues[END_REF]. However, evaluations made by different human beings judging the quality of contour detector outputs might be very diverse.

Considering the error quantified as the difference between a truth solution and a candidate edge image, quantitative measures are assigned numerical values. As described in [START_REF] Papari | Edge and line oriented contour detection: State of the art[END_REF], quantitative assessments are generally missing in proposed methods on edge detection. Nevertheless, in literature, some authors suggest [START_REF] Baddeley | An error metric for binary images[END_REF] [START_REF] Román-Roldán | A measure of quality for evaluating methods of segmentation and edge detection[END_REF] or compare different performance evaluations of contour detector measures [START_REF] Chabrier | A comparative study of supervised evaluation criteria for image segmentation[END_REF] [START_REF] Lopez-Molina | Quantitative error measures for edge detection[END_REF].

Digital images embody different type of edges. Crest lines (or called roof edges) represent a special type of contours, as shown in Fig. 1. Classical edge detectors are optimized to extract step or ramp edges [START_REF] Canny | A computational approach to edge detection[END_REF]; however, they fail to detect crest lines in images. Indeed, a step/ramp edge extractor will return two edges of both sides of the crest line because crest lines in an image represent narrow ridges or valleys on the image surface. Roof edges are defined as thin nets inside the image; describing roads [START_REF] Laptev | Automatic extraction of roads from aerial images based on scale space and snakes[END_REF] or rivers in satellite images, blood vessels in medical images [START_REF] Kirbas | A review of vessel extraction techniques and algorithms[END_REF] or roots in underground images. Therefore, finding these dense and thin structures is a significant task in image processing [START_REF] Ziou | Line detection using an optimal IIR filter[END_REF][15] [16][17][18] [START_REF] Magnier | Ridges and valleys detection in images using difference of rotating half smoothing filters[END_REF].

Concerning a curve in a grey level image I, ridges correspond maxima, valleys to minima [START_REF] Haralick | Ridges and valleys on digital images[END_REF] [START_REF] Tremblais | Algorithms for the extraction of various diameter vessels[END_REF]. Indeed, considering the image surface S, in the Cartesian space, S defines all pixel coordinates: S(x, y) = (x y I(x, y)) T . Let S T (x, y) = { S x (x, y), S y (x, y)} be the tangent plan of the surface S(x, y) in all points where:

S x (x, y) = ∂ S(x,y) ∂x = (1 0 I x (x, y)) T S y (x, y) = ∂ S(x,y) ∂y = (0 1 I y (x, y)) T ,
with I x and I y the partial derivatives of I respectively, along the x and the y axis. Ridges and valleys are given by the points where the values of S(x, y) are maxima or, respectively, minima in the orthogonal direction of the curve at (x, y).

Even though the evaluation of crest lines lingers an open problem as it is for edge detection, some formulas suggest efficient ways to evaluate edge detection. Note that the main difference between step/ramp edge detection and roof

Original image Image surface

Crest line edge/crest line extraction is that, in the first case, edges are chosen between two pixels, whereas in the second crest lines site in the middle of the ridge/valley (single pixel: maxima -top of ridges-or minima -bottom of valleys-). A non exhaustive review is proposed in the following Section, leading to a contour/crest line detectors assessment.

Contour

Homogeneous region

II. A REVIEW OF DISTANCE-BASED ERROR MEASURES

To provide an overall evaluation of the error assessment of edge detectors, the measure has to take into account both the amount of false positive points, false negative, but also of the displacement or distortion between the candidate edge image and the ground truth.

A. Image Pixel Positions Assessment

Position is the first criterion to be assessed over an edge image. Its basic evaluation is compounded of statistics. To that effect, ground truth edge image or ground true skeleton (T s ) and detected contours/crest lines (DC) are merged. Afterward, all points are partitioned into three sets:

• True Positive points (TPs), common points of T s and DC, • False Positive points (FPs), spurious detected edges of DC, • False Negative points (FNs), missing skeleton points of DC. First, FPs appear in the presence of noise, texture or other contours influencing the detection in the image. Then, FNs represent a hole in a contour of DC (i.e. a discontinuous edge), generally caused by a blur or a poor pronounced edge/crest line in the original image I. Finally, a failing threshold of the segmentation could create both FPs and FNs. Computing only FPs and FNs [START_REF] Magnier | Ridges and valleys detection in images using difference of rotating half smoothing filters[END_REF] [START_REF] Magnier | Multi-scale crest line extraction based on half gaussian kernels[END_REF] or combining these two statistics allows to display evaluations like Receiver Operating Characteristic (ROC) [START_REF] Bowyer | Edge detector evaluation using empirical roc curves[END_REF] or Precision-Recall (PR) [START_REF] Martin | Learning to detect natural image boundaries using local brightness, color, and texture cues[END_REF]. As illustrated in Fig. 7 (top left), these procedures do not use spatial interpretation because they consider only a one-to-one matching process.

In the results, presented in Section IV, TPs, FPs and FNs (respectively represented in green, red and blue, see illustrations in Fig. 2 (b),(c) and (d) for more details) are outlined the proportions of points by color set within the resultant image. Regardless the image dimension and non spatial tolerance, an indicator reflecting the spacial precision of points called KP I p can be created:

KP I p = 1 - 1 P b 1+ F N + 1-P b 1+ F P , with P b ∈ [0, 1], (1)
where In literature, a widely used similarity measure [START_REF] Abdou | Quantitative design and evaluation of enhancement/thresholding edge detectors[END_REF] between ground truth contours and detected edges, in respect with these criteria, is Pratt's Figure of Merit (F oM ):

F
F oM = 1 max(card(T s ), card(DC)) • p∈DC 1 1 + k • d 2 Ts (p) , (2)
with d Ts the distance between the p th pixel of DC and the nearest pixel of T s and k ∈ R * + a positive real. The edge detection is considered as good if F oM is close to 1 and poor when this value tends to 0. This measure, however, behaves abnormally, as shown in [START_REF] Peli | A study of edge detection algorithms[END_REF] and recalled in [START_REF] Lopez-Molina | Quantitative error measures for edge detection[END_REF], two different edge candidate images produce the same F oM , whereas one is visually closer to the ground truth. Moreover, as pointed out in [START_REF] Lopez-Molina | Quantitative error measures for edge detection[END_REF] and illustrated in Fig. 3, F oM is particularly sensitive to displacements, i.e. the errors grow up very fast when a small translation of the detected edges appears.

C. Baddeley's Delta Metric (BDM)

Inspired by Haudorff distance, Baddeleys Delta Metric (BDM) [START_REF] Baddeley | An error metric for binary images[END_REF] is a measure which intends to estimate the dissimilarity between each element of two edge images: The main drawback of the BDM is its hypersensitivity to false positive points, i.e. this measure tends to over penalize images with false detections. Indeed, when a false positive pixel is far from the true edge, the |w(d Ts (p)) -w(d DC (p))| value creates a high impact on the evaluation, thus penalizing the measure.

∆ k w (T s , DC) = 1 |I| •   p∈I |w(d Ts (p)) -w(d DC (p))| k   1/k , (3)

III. IMAGE SHAPE RESPECT ASSESSMENT

Even though the evaluation of crest lines lingers an open problem as it is for edge detection, some formulas suggest efficient ways to evaluate edge detection (see above). An evaluation of spatial position is far from being enough, since a dispersion (FP/FN points) around well-detected points (TP points) might be observed in the resultant image (see Fig. 8 for superposition of ground truth and segmentation/candidate images). To evaluate this dispersion on respect of image shape, a measure must be regarded. To that effect, we considered the Yasnoff et al. [START_REF] Yasnoff | Error measures for objective assessment of scene segmentation algorithms[END_REF] measure as a start:

M YAS = 100 card(F P ∪ T P ) • p∈F P d 2 Ts (p). (4)
This measure has been used for binary images in pattern recognition, compared to other methods and thus proving its efficiency. The first noticeable fact is the possible values of this measure are between zero and infinity. Because it is not plain to draw out conclusions about this measure within its value window, we perform a transformation to only observe values between zero and one and build a KP I (see Section II-A) so that one would correspond to infinity in the distance and error:

KP I : [0; ∞[ → [0; 1[ u → 1 -1 1+k•u h , (5)
where k ∈ R * + and h represents a power for observation. The parameter out of this formula corresponds to the distance d 2 Ts represented by u. In this work, a displaced edge expects to be penalized in function of the false pixels number (i.e. card(F P ∪ F N )) and in function of the distance to the position it should be located at. Indeed, since F N ∩F P = ∅, using the Euclidian distance, the d Ts is defined as:

∀ p ∈ F P ∪ F N, d Ts (p) = Inf { ||p -t||, t ∈ T s }, if p ∈ F P Inf { ||p -t||, t ∈ T P }, if p ∈ F N.
Inspired by Eq. 4, the first KP I formulation for shape respect is the following:

KP I 1 = 1 - 1 1 + √ p∈F P ∪F N d 2 Ts (p) card(F P ∪F N ) h . (6)
The latter by considering the whole image area, return non observable values by virtue of the image area importance. For example, two different images I 1 and I 2 both containing only one FP point at a same distance of T s (T s = 10 in I 1 and T s = 100 in I 2 ) obtain a common KP I 1 value (according to the formula) through the evaluation interpretation is different, because of the unlike ratio values. This led the optimum KP I research to a second formulation. That is the reason why we approached a formula to quantitatively observe dispersion regardless the image nature. This new term measures a pure dispersion of misplaced points compared to the cardinality of the ground truth T s :

KP I 2 = 1 - 1 1 + √ p∈F P ∪F N d 2 Ts (p) card(T S ) h . (7)
Thereafter, this consideration partially resolves the problem stated above. Indeed, considering two images with identical ground truth points T s and dispersion distance, but different ratio, so that in one image a single FP point behaves as n FP points (n>1) on another image; these images have the same interpretation in terms of KP I 2 which is not the case with ratios 1/100 and 90/100. The table I illustrates this drawback.

A same reasoning is applicable to FN points. That means the distance divided by the cardinality of T s must be balanced with the proportion of misplaced points, i.e. card(F P ∪ F N ).

As demonstrated above and illustrated in Fig. 4, KP I 1 and KP I 2 suffer from different drawbacks. In edge/crest line detection assessment, a performing KP I has to take into account all the following input parameters T s , F P , F N and d 2

Ts (p) (∀ p ∈ F P ∪ F N ). Ultimately, this new KP I is proposed by the following formula:

KP I 3 = 1 - 1 1 + √ p∈F P ∪F N d 2 Ts (p) card(T S ) • card(F P ∪F N ) card(T S ) h . (8)
This last formulation takes into consideration all observable cases and theoretically observable. Since distance evaluation remains local so that there is no difference in the evaluation of these two images, it is to be noticed that image dimensions are not considered. Moreover, a huge image may contain a skeleton with a small number of pixels. On the graphs presented in Fig. 4, we can notice how KP I 3 evolves with the image characteristics (number of ground truth and misplaced points), compared to KP I 1 and KP I 2 .

An undeniable parameter of KP I 3 formula is the power of the denominator term. It may be called a power of observation. Inasmuch as KP I 3 depends on its value, it evolves more or less quickly around 0.5 and embodies a range of observable cases. We determined average values for the distance term d Ts . The advice to choose values between 1 and 2 can be easily checked. Otherwise, the more KP I 3 evolution will be abrupt, the less the transition between 0.5 and 1 will be marked (i.e. the slope of the KP I 3 curve, for example h = 1 in Fig. 6). Moreover, fixing h = 1, KP I 3 stagnates far from 1 when d Ts becomes high. Additionally, when h = 2, KP I 3 starts to increase slowly and the slope becomes sharp around 0.5 to converge quickly towards 1. Finally, to fix the power at the golden ratio φ 1.6180339887 in order to ensure an evolution of KP I 3 that would not be too abrupt from 0 to 1 and also not penalize when d Ts is not elevated (contrary to KP I 3 with h = 1, see Fig. 6). Fig. 5 illustrates KP I 3 values in function of the distance of the mistake points and card(F P ∪ F N ). For a small skeleton T s = 200, it shows that KP I 3 draws near 1 quickly when card(F P ∪ F N ) and d Ts are growing. Alternatively, for a bigger skeleton, KP I 3 rises slower. These three evolution surfaces show the coherence and robustness of KP I 3 regarding our starting conditions: a displaced edge is penalized in function of the false pixel number and also of the distance to the position it should be located at.

IV. EXPERIMENTAL RESULTS

In natural images, obtaining true contours or crest lines is not trivial, for this reason, some authors create their own benchmark and then evaluate the detection. New crest line images are proposed in the rest of this paper to evaluate roof edge extractors 1 .

A. Construction of Test Images

Images are first created by building elementary shapes like lines or circles, since a skeleton can be decomposed into lines or arcs. Then, to create more complex images, we sum, then concatenate elementary images. Here is an example Fig. 7, a vertical line, a set of circular arcs and a tree.

B. Results of Quality Measures

In the proposed experiments, to create candidate images, ground truth images are progressively degraded with FNs (i.e. substitution of pixels in T s ), adding FPs (i.e. randomly noise pixels) or a displacement of each edge point of T s and then compared with the initial image. Indeed, Fig. 7 shows KP I {p,2,3} values in function of these 3 different deteriorations. Note that shape evaluation measurements (KP I {2,3} ) depend on the image nature, i.e. edge geometry, especially for displacements and the edge shapes can be interpreted as decent. Curves presented in Fig. 7 and values available in Table II indicate that KP I p (precision) is accurate and grows along the importance of the translations. It is the same for for KP I {2,3} (shape). However, on the one hand, KP I 2 is sensitive to displacements and FPs but, on the other hand, does not penalize enough FNs addition (particularly visible for the vertical line case and in Fig. 9 (Arch vs Circle).

Fig. 8 and9 present other images with their evaluations in Table II. Circle vs line and Comb images show how KP I 3 behaves in the presence of two spatially distant edges and also contours images having DC and T s close enough. Building Circle arc and Tree images with multiple sub images where each sub image possesses a valley which is different, thus a specific morphological operation is applied to the binary image 1 These images and the edge detection quality assessment are created using M atlab software and available on: http://www.mathworks.com/matlabcentral/fileexchange/58415-evaluationof-roof-edge-detectors-with-a-quantitative-error-measure until obtaining an image skeleton (involving a 8-connected neighborhood center element). Then, this skeleton is compared to ground truth valley images (presented in Fig. 7 left) and evaluated through the proposed KP I {p,2,3} . Table II states also when Circle arc and Tree images endure a translation before the morphological operation, obviously, KP I 2 penalizes sharply edge displacements (contrary to KP I 3 ). Results presented in these experiments illustrate the importance of a coherent evaluation. Indeed, penalizing a detected edge/crest line having either FPs and FNs or a small displacement is necessary proportionally to both the number of the mistake points and distance from ground truth. However, this penalty must not have a great influence in the presence of weak errors. Furthermore, the edge detector evaluation process must evolve from a good criterion to a poor criterion without steep translations. KP I 3 formula portrays better this state, and the use of this criterion enables to say whether a crest line extractor algorithm is efficient. Eventually, this crest line detector evaluation adapts to step/ramp edge detection evaluation, taking into account a spatial distance of the candidate edge from the true contours, contrary to the crest line where there is no distance (more details at the end of Section I).

V. CONCLUSION AND FUTURE WORKS

This work aimed at a roof edge extraction evaluation using ground truth contour image and candidate edge image. Combining the number of true pixels compared to the total number of false pixels and their distances to the ground truth enable to build a criterion which predicts if an image is correctly segmented or not. The proposed algorithm allows to evaluate precisely the efficiency of a crest line detection method for different types of images of different sizes. Additionally, the method proves its efficiency to measure under-and over-edge detection. This approach could be employed to determine the best fitting crest line detectors parameters, like thresholds in the last stage of the detection.

Future works suggest an evaluation of several methods used in the crest lines detections. Another work consists of a multiscale roof edge detection evaluation [START_REF] Magnier | Multi-scale crest line extraction based on half gaussian kernels[END_REF][21] [START_REF] Lindeberg | Edge detection and ridge detection with automatic scale selection[END_REF] [START_REF] Armande | Thin nets extraction using a multi-scale approach[END_REF]. Indeed, adapting the proposed evaluation would measure the accuracy of crest line scale detection of multi-scale images.

Fig. 1 .

 1 Fig. 1. Crest lines in an image.

Fig. 2 .

 2 Fig. 2. Crest lines extracted in red on the image in Fig. 1, and illustration of TP (green), FP (red) and FN (blue) points. In (c), DC is contaminated with 6 FPs ans 4 FNs, illustrated in (d). Note that F N ∪ T P = Ts.

  N and F P represent, respectively the proportion of FNs and FPs in the union of T s and DC (denoted by the usual set-theoretic operation T s ∪ DC) and P b represents the weight of FNs (in our experiments, P b =0.5). Finally, this type of formula is called Key Performance Indicator (KP I), with KP I ∈ [0, 1]. For a given experiment, a KP I value close to 1 means a poor segmentation. Alternatively, a KP I value close to 0 translates a good segmentation. B. Pratt's Figure of Merit criterion

  where k ∈ R * + and w represents a weighting concave function (in general w(x) = min( √ n 2 + m 2 , x) for an image of size m×n, with m and n ∈ N *[START_REF] Lopez-Molina | A gravitational approach to edge detection based on triangular norms[END_REF]).

FoMFig. 3 .Fig. 4 .

 34 Fig. 3. FoM is very sensitive to the translations, here, card(DC) = 200.

Fig. 5 .

 5 Fig. 5. Evolution of the KP I 3 in function of the mistake points distance and card(F P ∪ F N ), for different card(Ts) and fixing h = φ.

2 TsFig. 6 .

 26 Fig. 6. Evolution of the KP I 3 in function of the mistake points distance, for different powers h, with card(F P ∪ F N ) = 4000 and card(T S ) = 2200.

1 %Fig. 7 .

 17 Fig. 7. Evolution of the KP I {p,2,3} in function of the FPs, FNs and displacements. card(Ts) = 100, for the top image (100×101), card(Ts) = 2176, for the middle image (300×300) and card(Ts) = 2200, for the bottom image (580×300). Note that KP I 2 and KP I 3 curves are the same for pixel translation evaluation of the first image because card(Ts) = card(DC). FNs are created by adding horizontal lines of null pixel values and FPs by randomly adding isolated pixels of value 1. For a better visualization, black and white colors are inverted and the reader is referred to the web version of this article.

Fig. 9 .

 9 Fig. 9. Representation and ratio of TPs, FPs and FNs in: Circle vs line 100×100, Comb 100×100 and Arch vs line 201×201 pixels.

Fig. 8 .

 8 Fig. 8. Representation and ratio of TPs, FPs and FNs in: Circle arcs image 300×300 and Tree image 580×300. For a better visualization and an interpretation of the reference to color in legend, the reader is referred to the web version of this article.

TABLE II DIFFERENT

 II KP I {p,2,3} VALUES AND FOR THE IMAGES IN FIG. 8 AND 9

	Image	x translation	y translation	%TP	%FP	%FN	card(Ts)	card(F P ∪F N )	d 2	KP I P	KP I 2	KP I 3
		0	0	51.20	24.15	24.64	2176	1400	76	0.244	0.002	0.004
	Circle	+1	+2	10.79	44.35	44.86	2176	3488	519	0.446	0.174	0.089
	arcs	-2	+2	5.40	47.04	47.55	2176	3887	816	0.473	0.343	0.169
		+30	+20	1.22	47.10	51.69	2176	4063	1286	0.495	0.540	0.299
		0	0	47.22	24.16	28.61	2200	1531	153	0.265	0.007	0.013
	Tree	-5	+10	1.03	46.96	52.00	2200	4105	1811	0.496	0.667	0.422
		+5	+10	1.59	46.92	51.48	2200	4079	1927	0.493	0.687	0.447
		+50	+100	0.49	42.95	56.56	2200	3837	3175	0.507	0.817	0.644
	circle/line	0	0	0.26	74.02	25.71	100	828	384	0.620	0.996	0.968
	Comb	0	0	64.90	19.50	15.60	124	54	23	0.176	0.016	0.060
	arch/circle	0	0	50.10	0	49.90	1012	505	1531	0.5	0.388	0.661