Sylvain Coly

Pierre Druilhet

Nourddine Azzaoui

Efficient spatial designs for targeted regions or level set detection

Keywords: Optimal designs, computer experiment, kriging, Gaussian processes

Acquiring information on spatial phenomena can be costly and time-consuming. In this context, to obtain reliable global knowledge, the choice of measurement location is a crucial issue. Space-filling designs are often used to control variability uniformly across the whole space. However, in a monitoring context, it is more relevant to focus on crucial regions, especially when dealing with sensitive areas such as the environment, climate or public health. It is therefore important to choose a relevant optimality criterion to build models adapted to the purpose of the experiment. In this article, we propose two new optimality criteria: the first aims to focus on areas where the response exceeds a given threshold, while the second is suitable for estimating sets of levels. We introduce several algorithms for constructing optimal designs. We also focus on cost-effective algorithms that produce non-optimal but efficient designs. For both sequential and non-sequential contexts, we compare our designs with existing ones through extensive simulation studies.

Introduction

In many situations, collecting spatial data can be costly or time-consuming. This is the case in epidemic propagation studies, in the spatial control of pollutants, in climatic phenomena, and so on. When the number of observations is limited, it is crucial to optimize the positions of sampling locations, taking into account prior information, generally obtained from previous or indirect observations or diffusion models. In the case of sequential monitoring, the next sampling locations can be selected on the basis of up-to-date information. For single-stage sampling, we can only rely on prior information to select sampling positions.

Space filling designs are the most popular non-sequential designs when no prior information is available. They aim to ensure uniformly accurate estimation of the phenomenon over the whole space, see, e.g. [START_REF] Niederreiter | Point sets and sequences with small discrepancy[END_REF][START_REF] Johnson | Minimax and maximin distance designs[END_REF]. For example, distance-based approaches such as Minimax, Maximin or other discrepancy criteria have been proposed for the construction of optimal space filling designs. To reduce computation costs, low discrepancy sequences such as Halton, Hammersley, Sobol and Faure sequences have provided a first attempt to build easy to compute designs. For high dimensional problems, Latin Hypercube Designs and orthogonal arrays have been introduced to ensure equilibrium on axes [START_REF] Stein | Large sample properties of simulations using latin hypercube sampling[END_REF][START_REF] Owen | Randomly orthogonal arrays for computer experiments, integration and visualization[END_REF].

In the context of Gaussian fields or Kriging methods, entropy [START_REF] Shewry | Maximum entropy sampling[END_REF][START_REF] Currin | Bayesian prediction of deterministic functions, with applications to the design and analysis of computer experiments[END_REF][START_REF] Bueso | A state-space-model approach to optimal spatial sampling design based on entropy[END_REF][START_REF] Angulo | Space-time adaptive sampling and data transformations[END_REF]] and Integrated Mean Square Error (IMSE) [START_REF] Sacks | Design and analysis of computer experiments[END_REF] are the two main criteria proposed for optimal sensor deployment. Numerous alternatives have been proposed in recent decades: Conditional Minimizer Entropy (CME) [START_REF] Villemonteix | An informational approach to the global optimization of expensive-to-evaluate functions[END_REF], Generalized Relative Complexity [START_REF] Alonso | Dependence assessment based on generalized relative complexity: application to sampling network design[END_REF], Expected Improvement (EI) [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF] and Expected Improvement Gain (EIG) [START_REF] Ryan | Estimating expected information gains for experimental designs with application to the random fatigue-limit model[END_REF]. All these methods aim to determine the design carrying the greatest expected amount of information [START_REF] Mateu | Spatio-temporal Design: Advances in Efficient Data Acquisition[END_REF]. Recently, many variants of Expected Improvement (EI) [START_REF] Bichon | Efficient global reliability analysis for nonlinear implicit performance functions[END_REF][START_REF] Ranjan | Sequential experiment design for contour estimation from complex computer codes[END_REF], such as the quasi-Expected Improvement (q-EI) [START_REF] Chevalier | Fast computation of the multi-points expected improvement with applications in batch selection[END_REF] have been introduced and widely used in industrial fields.

The approaches cited above aim to control the global variance without considering the expected values of the response. The aim of the paper is to propose optimal designs based on criteria targeting an area of interest. The approaches cited above aim to control the global variance without considering the expected values of the response. The aim of the paper is to propose optimal designs based on criteria targeting an area of interest. We mainly consider two cases: the first aims to target regions where the response values is significantly high. The second, more suited to sequential designs, aims to estimate a given set of levels.

In section 2, we motivate and define the optimality criteria used to target regions of interest. In section 3, we propose some algorithms for obtaining efficient designs. In section 4, we carry out simulation studies for sequential designs when the target area is a level set. We propose several performance indicators to compare our designs with those obtained in [START_REF] Picheny | Adaptive designs of experiments for accurate approximation of a target region[END_REF]. We then discuss the relative merits of each method. In section 5, we computationally evaluate the effectiveness of our methods for single-stage designs on a few examples.

Optimality criteria focusing on areas of interest

We consider a grid E of size N × N , where N is an integer. We denote by y(x) ≥ 0 the variable of interest at x ∈ E, and by y the N 2 -vector with entries y(x). We consider a grid E of size N × N , where N is an integer. We denote by y(x) the true response at point x in E, and by y the N 2 -vector with entries y(x).We assume that knowledge of y can be modeled by a Gaussian field:

y ∼ N (µ, Σ), (1)
which corresponds to prior knowledge from a Bayesian perspective, or, in the field of computer experiments, to a meta-model. For x ∈ E, we denote µ(x) = E(y(x)) and σ 2 x = Σ x,x = Var(y(x)). Let d = {x i } 1≤i≤n be a n-point design and y d = (y(x i )) 1≤i≤n be the nvector y restricted to d. Knowing the response y d on the design points, the updated knowledge on the field is given by the linear updating formula:

µ d|y d = E(y d|y d ) = µ d + Cov (y d, y d ) Var -1 (y d ) (y d -µ d ), (2)
and for the variance: In most applications, measurement errors are negligible compared to spatial variability. Therefore, we assume that Var (y(x i )|d) = 0 for x i belonging to d and we will confound y(x i ) with its measurement.

Var (y d|y d ) = Var (y d) -Cov (y d, y d ) Var -1 (y d ) Cov (y d , y d) , (3)
The aim of the paper is to construct a design d that provides accurate knowledge on the area of interest. In order to define optimality criteria, we first define for any x in E the weighted variance as in [START_REF] Picheny | Adaptive designs of experiments for accurate approximation of a target region[END_REF] :

c(x; d) = w(x) × Var (y(x)|d) , (4)
where w(x) is a weight function that depends on the area of interest. Two types of global criteria can be derived from this weighted variance function: a max-criterion

MC(d) = max x∈E c(x; d), (5)
and an integrated criterion

IC(d) = x∈E c(x; d). (6)
Then, we seek the design d * that minimizes either MC(c) or IC(d). In the next sections, we propose several weight functions, depending on the goal of the experiment and the target zone. Note that a design d * that minimizes MC(d) also minimizes max x∈E h(c(x; d)) for any increasing transformation h. Such monotonic invariance property does not hold for IC(d).

In the next sections, we propose two different weight functions, depending on the goal of the experiment and the target area.

Level set detection

Here, we aim to estimate the level set L = {x ∈ E ; y(x) = T } associated to a given threshold T > 0. We propose the weight function defined by

w ls T (x) = 2 1 2 -F µ x -T σ x . ( 7
)
where F is the cumulative distribution function of the standard normal distribution. There are two possible interpretations of the weight function w ls T (x). The first one is frequentist: consider y(x) as the unknown fixed quantity to be evaluated and µ(x) as the realization of a normally distributed random variable with mean y(x) and variance Var(y(x)). In that case, w ls T (x) corresponds to the p-value of the two-tailed test : H 0 : "y(x) = T " vs H 1 : "y(x) ̸ = T ". The weight w ls T (x) is close to 1 when the hypothesis "x belongs to the level set" is rejected and is close to 0 otherwise. The second interpretation is Bayesian: put a flat prior on y(x) that reflects the ignorance on y(x). Consider the HPD-credible set C β (x) of y(x) with credible level β, then w ls T (x) is the maximal credible level β such that T does not belongs to C β (x). This corresponds to the construction of Bayesian two-sided hypothesis testing based on confidence intervals (see [START_REF] Lindley | Introduction to probability and statistics from a Bayesian viewpoint. Part II: Inference[END_REF][START_REF] Thulin | Decision-theoretic justifications for Bayesian hypothesis testing using credible sets[END_REF])

We denote by MC ls (d) and IC ls (d) the max and integrated criteria related to w ls T (x).

In [START_REF] Picheny | Adaptive designs of experiments for accurate approximation of a target region[END_REF], another weight function has been proposed to estimate level sets. It is denoted by w σ 2 ϵ (x) in this paper and defined by

w σ 2 ϵ (x) = 1 (2π(σ 2 ϵ + σ 2 x )) 1/2 exp{-(µ x -T ) 2 / 2 σ 2 ϵ + σ 2 x }, (8)
where σ 2 ϵ is a smoothing parameter that needs to be calibrated. We denote by MC W (d) and IC W (d) the max and integrated criteria derived from w σ 2 ϵ (x). The authors in [START_REF] Picheny | Adaptive designs of experiments for accurate approximation of a target region[END_REF] consider only IC W . We will show on examples that the performance of designs based on MC W (d) are less efficient than that based on IC W (d).

When

σ 2 ϵ = 0, w σ 2 ϵ (x) is unbounded for σ 2
x close to 0 (see Fig. 1a). This occurs for points located in the close neighborhood of an observation point and therefore, the future design points will be concentrated around already observed points. At the opposite, when σ 2 ϵ is large, weights tend to be uniform over E, resulting in a space-filling design. The choice of σ 2 ϵ will therefore influence the related optimal design.

As shown in Fig. 1b, the main difference between the weight functions w ls T (x) and w σ 2 ϵ (x) lies in their behaviors when uncertainty is large. For a given value of |µ x -T | and a large value of σ 2

x , w σ 2 ϵ (x) ≈ 0 whereas w ls T (x) ≈ 1. In section 4, we compare sequential designs build w.r.t. w ls T (x) and w σ 2 ϵ (x) by simulation studies. Note that the weight function w ls T (x) has no parameter to calibrate.

Areas exceeding a given threshold

Here, we aim to build designs that give an accurate knowledge of the area {x ; y(x) > T } where the response y(x) exceeds a given threshold T . In that case, we propose the weight function defined by

w exc T (x) = P(y(x) > T ) = F µ(x) -T σ x . (9)
We give two interpretations of the weight function w exc T (x) as above: w exc T (x) is the p-value of the one-tailed hypothesis test H 0 : "y(x) > T " vs H 1 : "y(x) ≤ T ". In the Bayesian interpretation with a flat prior on y(x),

w exc T (x) = P(y(x) > T )
is the posterior distribution of the event "y(x) > T ", given µ x . Therefore, w exc T (x) corresponds to the faith that x belongs to the target area.

We denote by MC exc (d) and IC exc (d) the max-criterion ( 5) and the integrated criterion (6) related to w exc T (x).

Algorithms for optimal designs

We review some algorithms for obtaining optimal designs according to a generic criterion C(d) that is to be minimized. In our applications, C(d) will be either MC(d) = max x∈E c(x; d) or IC(d) = x∈E c(x; d), where c(x; d) is defined by (4). We consider two types of design: non-sequential and sequential.

Non-sequential designs rely entirely on prior information since all observation points are chosen before the experiment. They are suitable for parallelized computer experiments or sampling campaigns where responses are analyzed after the fact.

In adaptive sequential designs, on the other hand, responses on design points are observed after each stage. Therefore, the optimal designs points at a given stage can be based on the meta-model updated by the observations from the previous stages.

Non-sequential design

The standard strategy for computationally obtaining an optimal design is to iteratively improve a fictitious starting design d (0) using an exchange algorithm (see [START_REF] Jin | An efficient algorithm for constructing optimal design of computer experiments[END_REF] for a review). Usually, the starting design d (0) is a standard space-filling or random design. However, these designs can be far from optimal and therefore require a large number of iterations to improve. We propose here a starting design, denoted d † , based on a cost-effective algorithm. Our simulation studies show that it is highly efficient (see figures 9 and 10).

Construction of the starting design d †

The starting design d (0) = d † = {x 1 , x 2 , . . . , x n } is constructed from an empty design by sequentially adding one point at a time. Each new point maximizes (4) with updated variance and unchanged mean functions. Note that this algorithm only requires the calculation of n -1 updated covariance matrices.

Start with the empty design d 0 = ∅ and a given meta-model M 0 : N (µ, Σ). Define the weight function w 0 (x), as in ( 9), ( 7) or ( 8) depending on the goal of the experiment. Choose the point x 1 that maximizes h 0 (x) = w 0 (x) × Var (y(x)|d 0 ) and set d 1 = {x 1 }. Define the updated metamodel M 1 , with the same mean function µ(x) as the initial meta-model M 0 and with the updated variance Σ 1 given by ( 3) with d = d 1 . Note that, since the point x 1 is not observed, the mean function cannot be updated. Based on the new meta-model M 1 , we can define a new weight function w 1 (x).

At the second step, choose x 2 that maximizes w 1 (x) × Var (y(x)|d 1 ). Put d 2 = {x 1 , x 2 }, define a new meta-model M 2 with the same mean function and the updated variance matrix given d 2 . Then, define the new weight function w 2 (x), and so on... After n iteration, we get the n points {x 1 , ..., x n } of d † . This is detailed in Algorithm 1.

Note that the weight function w i (x) is based only on the mean function and the updated variance function var(y(x)|d i-1 ). Therefore, it is not necessary to compute the entire updated variance matrix, only the diagonal terms.

Exchange algorithm

From a starting design d (0) , the exchange algorithm, see [START_REF] Kennard | Computer aided design of experiments[END_REF][START_REF] Mitchell | An algorithm for the construction of "d-optimal" experimental designs[END_REF], displayed in Algorithm 2, consists of iterative permutations of the design's inner and outer points. A permutation is accepted if the corresponding criterion is improved. We denote d (k) the design obtained after k iterations.

In addition, simulated annealing [START_REF] Morris | Exploratory designs for computational experiments[END_REF] can be used within the exchange algorithm to avoid local extrema issues. However, we did not observe any

(x) from M 0 . foreach i from 1 to n do foreach x in E do Compute h i (x) = w i-1 (x) × Var (y(x)|d i-1 ) ; end Choose x i = ArgMax x∈E (h i (x)) ; Put d i = {x i } ∪ d i-1 ; Put M i : N (µ, Σ i ), where Σ i is the updated covariance w.r.t. d i (µ is unchanged). Define the function w i (x) from M i . end Output: design d † = d n .
improvement over the exchange algorithm in our simulation studies.

Sequential designs

A sequential design is divided in N stages. At Stage i, we have to build a n i -point design d i based on the updated meta-model, which consists of an initial meta-model updated with the observations obtained in the previous stages (see [START_REF] Williams | Sequential design of computer experiments to minimize integrated response functions[END_REF][START_REF] Bect | Sequential design of computer experiments for the estimation of a probability of failure[END_REF][START_REF] Villanueva | Sequential experimental design and model calibration for targeted events[END_REF]). The new design d i can be constructed by using Algorithm 1 and 2 with the updated meta-model obtained with Formulas (2) and (3). Examples of sequential designs are given in Section 4.

Efficiency factor

For any n-point design d, the theoretical efficiency factor of d is defined by

eff(d) = C(d opt ) C(d) .
where d opt is the theoretical n-point design that minimizes the criterion C(d). Since d opt is not attainable, we approximate eff(d) by

eff(d) = C(d * * ) C(d) . ( 10
)
Algorithm 2: Exchange algorithm.

Input: Choose a starting design d (0) , maximal number of iterations

M ; foreach k from 1 to M do Randomly draw x ∈ d (k-1) ; Randomly draw x ′ ∈ d(k-1) ; Permute x and x ′ considering d temp = d (k-1) ∪ {x ′ }\{x} ; if C (d temp ) < C d (k-1) then d (k) = d temp ; else d (k) = d (k-1) ; end end Output: design d * = d (M ) .
where d * * is computed as follow: we draw 1000 n-point designs (d

Simulation studies for sequential designs

In this section, the aim of the experiment is to provide an accurate estimation of the level set L = {x : y(x) = T } for a given T on a 150 × 150 grid. We consider sequential adaptive designs with one additional point at each stage. Once the design point is observed, the level set is estimated by L = {x ∈ E ; y(x) = T } where y(x) is the mean of the updated meta-model.

Before starting the experiments, we assume no prior information on the mean. In this case, we recommend to choose, for all x ∈ E, µ(x) = T as a default choice for the mean function of the meta-model (1). This choice implies that any point x potentially belongs to the level set L. For the covariance matrix of the meta-model, we assume that Σ the Matérn covariance function

cov(x i , x j ) = 2 1-ν σ 2 (d √ 2ν/κ) ν K ν (d √ 2ν/κ)/Γ(ν) (11)
with σ = 0.7, ν = 0.7 and κ = 0.2, where d = d(x i , x j ) is the Euclidean distance between x i and x j , Γ is the Gamma function and K the modified

Bessel function.

As we have no prior information on the mean, we start with a space-filling design with few points. Once the design points are observed, we update the mean and the variance of the meta-model by applying Formula 2 and 3. Then, we start the sequential design with one of our criteria, IC ls or MC ls , based on the updated meta-model. At each stage, we construct the 1-point design that minimizes the updated criterion.

We compare our designs with those obtained by minimizing the criteria IC W or MC W . To calibrate the parameter σ 2 ε , we follows the recommendation given in [START_REF] Picheny | Adaptive designs of experiments for accurate approximation of a target region[END_REF], i.e. we choose choose σ 2 ε = (max x µ(x) -min x µ(x))/20 which is updated at each stage.

To compare the relative performance of the models in terms of level set estimation accuracy, we propose three quality scores below:

Score Q dist : distance average
This score evaluates the symmetric distance between the estimated and actual level sets. It is defined by :

Q dist = 1 2 Q actual dist + Q est dist ,
where

Q actual dist = 1 #{L} x∈L min x ′ ∈ L d(x, x ′ ),
is the average distance between each point of the actual level set and the nearest point in the estimated level set (see Fig. 3a). Symmetrically,

Q est dist = 1 #{ L} x∈ L min x ′ ∈L d(x, x ′ ),
is the average distance between each point in the estimated level set and the nearest point in the actual level set.

Score Q value : values average

This score evaluates the average discrepancy between the estimated and actual values on both the actual and estimated level sets:

Q value = 1 2 Q actual value + Q est value ,
where

Q est value = 1 #{ L} x∈ L |y(x) -T |,
and

Q actual value = 1 #{L} x∈L | y(x) -T |.

Score Q area : area between the two level sets

This score corresponds to the proportion of area between the actual level set L and the estimated one L, (see Fig. 3b):

Q area = 1 N 2 #
x ∈ E : y(x) < T and y(x) > T or y(x) > T and y(x) < T .

It can also be seen as the proportion of points misclassified in the sets {x : y(x) > T } and {x : y(x) < T }.

Other scores

Both Q dist and Q value scores are defined as symmetrized average distances. The individual components Q est value , Q actual value , Q est dist and Q actual dist and their maximal alternatives have also been tested as quality scores. For smooth one-piece level sets, they give similar results to those obtained with Q dist and Q value . When the target level set is disconnected, they are difficult to interpret. So, they are not considered in this paper.

Example 1: smooth level set

We consider a smooth field displayed in Fig. 2a with the smooth level set L = {x : y(x) = T }. In this example, we choose T = 0.85. First, we perform a 4-point space-filling minimax design to acquire information. Then, for each criterion MC ls , IC ls , MC W and IC W , we construct the sequential design. The estimated level sets L = {x : y(x) = T } are displayed in Fig. 4 for 8-point designs, which correspond to four stages after the starting space-filling design.

In Fig. 5, we compare the performance of the four strategies against the three quality scores proposed above. There is no significant difference between the criteria, except with MC W which seems to be globally less efficient. After 12 stages, all the 16-point designs have the same performance.

Example 2: complex level set

In this example, we consider a more complex field displayed in 2b. The level set to be estimated has three irregular disconnect components. First, we perform a 3-point minimax design. Then, for each criterion MC ls , IC ls , MC W and IC W , we construct a sequential design with one additional point at each stage. After 3 stages, all the 6-point designs identify the upper left level set. The 6-point design based on MC ls also identifies a second level set in the righthand part of the area (see Fig. 6). The designs based on the two integrated criteria need four additional stages to identify a second component (see Fig. 6) whereas the design based on MC W requires eight additional stages.

After just 7 stages, the 10-point design based on MC ls is the only one capable of finding the three components of the level set (see Fig. 6). The three components are correctly identified by the designs based on the integrated criteria for the first time with a 14-point design (not displayed in the paper) and by the MC W criterion with a 17-point design.

The four 20-point designs approximately identify the shape of the three components of the level set. For the three quality scores (Fig. 7), the two integrated criteria give slightly better results.

At any stage, the three quality scores consistently show that MC W provide poor designs. The integrated criteria IC W and IC ls provide the best designs w.r.t. the quality scores Q dist and Q area for more than 14 points. Up to 13 points, our criterion MC ls provides the best designs w.r.t. Q value and Q dist .

Comments

The examples presented above show that all four criteria provide highly relevant models for finding smooth level sets. For irregular and disconnected level sets, the models behave very differently, depending on the criterion.

Sequential designs based on the max-criterion MC ls appear to be more efficient in detecting disconnected components of the target level set more quickly, especially when they are located close to the domain boundary. Con- sequently, MC ls seems to be the best option in the early stages of sequential design. Conversely, the maximum criterion MC W usually gives poor results. Designs based on the integrated criteria IC ls and IC W give, in average, similar results. Overall, the integrated criteria lead to designs that are more space filling designs than those based on maximal criteria: they efficiently control uncertainty over the whole area.

6-point designs

For both examples several initial designs were tested, including 3/4/5/6point space-filling designs and random designs. When all the space-filling design points are located on the same side of the level set, results similar to those presented here were obtained. When the points are located on either sides of the level set, each criterion has resulted in equivalent designs in terms of quality scores.

Simulation studies for non-sequential designs

In this section, we propose some simulation studies that show how our criterion influence the location of the designs points. We consider here a Gaussian random field y ∼ N (µ, Σ) on a 50 × 50 grid E over [0, 1] 2 , where Σ is defined, as in Section 4 by

cov(x i , x j ) = 2 1-ν σ 2 (d √ 2ν/κ) ν K ν (d √ 2ν/κ)/Γ(ν) (12)
with σ = 0.7, ν = 0.7 and κ = 0.2, where d = d(x i , x j ) is the Euclidean distance between x i and x j , Γ is the Gamma function and K the modified Bessel function. The mean function is given by the following formula µ(x) = 2 × exp -{(x 1 -1) 2 + 3(x 2 -0.5) 2 } 1/2 /3 , with x = (x 1 , x 2 ), see Fig. 8. Mean and variance functions can be interpreted as a prior knowledge on y(x) based, for example, on previous observations or on forecasting models.

Efficient design for area exceeding a given threshold

We aim to construct a 10-point design that focuses on the area where the expected values exceed a given threshold T . In our example, we choose T = 0.85.

We consider two criteria, MC exc T and IC exc T , obtained from ( 5) and ( 6) by choosing the weight function w exc T (x). For each criterion, we start with the design d † obtained by Algorithm 1. Then, we apply the exchange algorithm with 10 000 iterations to obtain the design d * . The resulting designs are displayed in Fig. 8a.

In Fig. 9a, we compare the efficiency factors [START_REF] Villemonteix | An informational approach to the global optimization of expensive-to-evaluate functions[END_REF] of the designs d † , d * and usual space filling designs. Space filling designs are obtained from the R packages randtoolbox, minimaxdesign and maximin. As these packages use optimization algorithm which include a degree of randomness, we ran times each package several. For each type of space filling designs, we display the box plots of the efficiency factors.

We can see that the starting design d † is highly efficient w.r.t. both MC exc T and IC exc T . We can also observe on our example that space filling designs are less efficient against MC exc T than w.r.t. IC exc T . As it will be seen in Section 4, this illustrates the fact that integrated criteria lead to more space filling designs than max-criteria.

Efficient design to detect level sets

We aim to estimate the level set y(x) = T with T = 0.85. We use the same meta-model as in Section 5.1 with the same parameters. So, we seek the optimal 10-point design w.r.t. the MC ls T and IC ls T criterion, where

IC ls = x∈E w ls T (x) × Var (y(x)|d) or MC ls = max x∈E {w ls T (x) × Var (y(x)|d)}.
Clearly, a sequential design is far more appropriate to this situation, since observations are needed to estimate the level set (see Section 4). However, a non-sequential design can be seen as a stage of a group-sequential design based on the meta-model updated with the previous observation.

The resulting designs are displayed in Fig. 8b. In Fig. 10, we compare their efficiency factors with usual space filling designs. Similarly to Section 5.1, designs d † obtained by Algorithm 1 are highly efficient with an efficiency factor 0.993 for MC ls and 0.94 for IC ls . In this example, space-filling designs are highly efficient for IC ls , but not for the MC ls criterion.

Further experiments were carried out by alternately considering 3-point designs and 20-point designs. As with 10-point designs, d † is more relevant than space-filling designs (in particular for the MC ls criterion) and is slightly improved by the exchange algorithm. However, as the number of design points increases, the differences of efficiencies between the designs are reduced.

Conclusion

We proposed new criteria to build optimal designs that aim to accurately estimate the response over a given target area when the response is modelled by a Gaussian field. From a Bayesian point of view, these criteria are based on the faith the a point belong or not to the are of interest. From a frequentist point of view, they are based on the test of whether the point belongs to the target area or not.

When the aim of the experiment is to estimate a level set, we have proposed three quality scores to evaluate and compare the performance of the designs. In the case of sequential designs, the max-criterion MC ls T appears to be more effective than the other criteria in exploring areas of high uncertainty and therefore to detect disconnected areas. Integrated criteria are more space filling and attempt to control overall uncertainty, with larger weights on the already explored areas. This suggests a hybrid strategy for investigating a target area using sequential designs: first, use the max-criterion MC ls to quickly identify areas of interests. Then, after several stages, use the integrated criterion IC ls or to accurately control the overall uncertainty on that area. The evaluation of the performance of this approach is left for future research.

For non-sequential designs, optimal designs are highly dependent on prior information, so there is no natural way of comparing designs with each other. For this situation, we proposed a non-optimal but efficient design, d † , based on a computationally-cheap algorithm.

  where d = E\d is the complement of d, y d the vector of responses restricted to d and µ d = E(y d ). Since Var (y d|y d ) does not depend on the actual value of y d but only on the position of the observation points, we write Var (y d|d) instead of Var (y d|y d ). Similarly we write Var (y(x)|d) instead of Var (y(x)|y d ).

1 Figure 1 :

 11 Figure 1: Comparison of weight functions w ls T (blue, straight) and w σ 2 ϵ for σ 2 ϵ = 0.005 (red, dashed), σ 2 ϵ = 0.2 (green, dot-dashed) and σ 2 ϵ = 0.8 (orange, dotted)

Algorithm 1 :

 1 Construction of the n-point design d † . Input: Put d 0 = ∅ the empty design ; Initialize the prior meta-model as M 0 : N (µ; Σ) Define the function w 0

  where the n points are chosen at random over the grid. For each design d (0) i , we apply the exchange algorithm and obtain a design d * i . The design d * * is the design that minimizes C(d) among (d * 1 , ...d * 1000 ).

Figure 2 :

 2 Figure 2: Map of y(x) and level set L = {x : y(x) = 0.85} (black line).

Figure 3 :

 3 Figure 3: Quality scores for estimated level set

  (a) Map for MC ls . (b) Map for IC ls . (c) Map for MC W . (d) Map for IC W .

Figure 4 :

 4 Figure 4: comparison of y(x) maps for 8-point designs (black dots) for Example 1. The actual (white) and estimated (black) level sets are displayed.

  (a) Q value score (b) Q dist score (c) Qarea score

Figure 5 :

 5 Figure 5: quality scores against the number of design points w.r.t criteria MC ls (black), IC ls (blue), IC W (red) and MC W (green) for Example 1.

  (a) Map for MC ls . (b) Map for IC ls . (c) Map for MC W . (d) Map for IC W . 10-point designs (e) Map for MC ls . (f) Map for IC ls . (g) Map for MC W . (h) Map for IC W . 13-point designs (i) Map for MC ls . (j) Map for IC ls . (k) Map for MC W . (l) Map for IC W .17-point designs(m) Map for MC ls . (n) Map for IC ls . (o) Map for MC W . (p) Map for IC W . 20-point designs (q) Map for MC ls . (r) Map for IC ls . (s) Map for MC W . (t) Map for IC W .

Figure 6 :

 6 Figure 6: comparison of y(x) maps for several n-point designs (black dots) for Example 2. The actual (white) and estimated (black) level sets are displayed.

  (a) Q value score (b) Q dist score (c) Qarea score

Figure 7 :

 7 Figure 7: quality scores against the number of design points w.r.t criteria MC ls (black), IC ls (blue), IC W (red) and MC W (green) for Example 2.

(a) 10 -

 10 point designs w.r.t. MC exc (black) and IC exc (white) (b) 10-point designs w.r.t. MC ls (black) and IC ls (white)

Figure 8 :

 8 Figure 8: Map of µ(x), level set {x : µ(x) = 0.85} (black line) and 10-point designs.

  (a) for MC exc T . (b) for IC exc T .

Figure 9 :

 9 Figure 9: Efficiency factors of usual designs w.r.t. MC exc T and IC exc T . The purple line corresponds to the efficiency of d † .

Figure 10 :

 10 Figure 10: Efficiency factors of usual designs w.r.t. MC LS T and IC LS T . The purple line corresponds to the efficiency of d † .

Preprint submitted to

Acknowledgment

This research was financed by the French government IDEX-ISITE initiative 16-IDEX-0001 (CAP 20-25).


Pranevicius

Influence of surface preparation and ion flux on the nitriding efficiency of austenitic stainless steel