Gilles Chabriel
email: chabriel@univ-tln.fr-barrere@univ-tln.fr

Jean Barrère

SECOND-ORDER BLIND SOURCE SEPARATION: A NEW EXPRESSION OF INSTANTANEOUS SEPARATING MATRIX FOR MIXTURES OF DELAYED SOURCES

In this paper, we study the blind separation of mixtures of propagating waves (delayed sources) encountered for example in underwater telephone (UWT) systems . We suggest a new second-order statistics method using as many observations as sources. First, we show that each of the N delayed sources can be developed as a particular linear combination of the different temporal-derivatives of the N observations. Under some assumptions, an instantaneous rectangular separating matrix is then identified by the joint diagonalization of a set of covariance matrices estimated from the observations and its derivatives. The algorithm used takes into account the particular structure of the spectral mixing matrix encountered. A numerical simulation is provided in a 3-sources/3-observations case for propagating audio signals.

INTRODUCTION

Consider a set of N propagating waves s j (t) in an echo-free and noise-free medium. These are recorded on a set of M identical sensors. Let x j (t) denote the contribution of the source s j (t) on a sensor arbitrary indexed by 1, i.e. x j (t) is a filtered version of the source s j (t). Assuming that the set of sensors is sufficiently compact, then the contribution of the source s j (t) recorded on a sensor i is the same one as that recorded on the first sensor, except for an attenuation factor and a propagation delay. The observation y i (t) recorded on the sensor i, is a linear combination of the delayed contribu-tions x j (t):

y 1 (t) = N j=1
x j (t),

y i (t) = N j=1 c ij x j (t -τ ij ), i ∈ [2, . . . , M ], (1)
where c ij is the relative attenuation coefficient of source j on sensor i, and τ ij a relative propagation delay of source j on sensor i.

The different contributions are assumed to be band-limited, differently colored, statistically independent and zero-mean. Thereafter in the paper, the contributions will be called sources.

The problem is to provide an estimation of the N different sources x j (t) from the observations y i (t).

A first way to solve the problem using truncated Taylor series of each delayed source x j (t -τ ij ) has been successfully treated in [START_REF] Chabriel | Blind Identification of Slightly Delayed Mixtures[END_REF]- [START_REF] Yeredor | Blind Source Separation with Pure Delay Mixtures[END_REF]. It consists of the estimate of an instantaneous square mixing matrix M M ×M between the observations and the different derivatives of the sources:

[y 1 , y 2 , . . .] T = M [x 1 , x 2 , . . . , ẋ1 , ẋ2 , . . .] T .
For low delays, second-order statistics are sufficient to extract the different contributions from the observations (see [START_REF] Barrère | A Compact Sensor Array for Blind Separation of Sources[END_REF], [START_REF] Cauwenberghs | Blind Broadband Source Localization and Separation in Miniature Sensor Arrays[END_REF], [START_REF] Stanacevic | Micropower Gradient Flow VLSI Acoustic Localizer[END_REF]). For higher delays, a combination of two-order and highorder statistics methods is proposed in [START_REF] Chabriel | An Instantaneous Formulation of Mixtures for Blind Separation of Propagating Waves[END_REF] to achieve the separation.

The main limitation of these methods is the high necessary number of observations (up to several times the number of sources).

We propose here a N -sources N -sensors original approach based on the direct estimation of an instantaneous rectangular separating matrix S between the different derivatives of the observations and filtered sources:

x h 1 , x h 2 , . . . x h N , T = S [y 1 , ẏ1 , . . . , y 2 , ẏ2 , . . . , y N , ẏN , . . .] T ,
where x h j is the source x j filtered by an invertible known filter h.

It is important to note that the two approaches are not equivalent: the matrix M is not a simple pseudoinverse of S.

The paper is organized as follows: in the next section we explain how to construct the separating matrix S from the propagation model described by the system of equations [START_REF] Chabriel | Blind Identification of Slightly Delayed Mixtures[END_REF]. Then the third section explains the second-order statistics method implemented to estimate this separating matrix. The estimation uses a new version of the joint-diagonalization algorithm proposed by S.Dégerine in [START_REF] Dégerine | Sur la diagonalisation conjointe approchée par un critère des moindres carrés[END_REF]. A description of this algorithm has to appear in english, a preprint can already be read in [START_REF] Dégerine | A new non-orthogonal approximate joint diagonalization algorithm for blind source separation[END_REF]. A numerical simulation in the last section illustrates the effectiveness of our approach.

EXPRESSION OF THE SEPARATING MATRIX

In the frequency domain, the system (1) becomes

Y 1 (ν) = N j=1 X j (ν), Y i (ν) = N j=1 c ij e -2πντij X j (ν), i ∈ [2, . . . , N ].
(2) X j (ν) and Y i (ν) are respectively the Fourier Transforms (FT) of the j th source and of the i th observation.

The system (2) can be rewritten in matrix notation as:

Y(ν) = M f (ν)X(ν), (3)
where

Y(ν) = [Y 1 (ν), . . . , Y N (ν)] T , X(ν) = [X 1 (ν), . . . , X N (ν)] T .
The N × N spectral mixing matrix M f (ν) is defined as:

M f (ν) =      1 . . . 1 c 21 e -2πντ21 . . . c 2N e -2πντ2N . . . . . . . . . c n1 e -2πντN1 . . . c N N e -2πντNN      . (4)
We assume that M f (ν) is regular for any frequency. The inverse matrix of M f (ν) is

M f (ν) -1 = 1 det M f (ν) (adj M f (ν)) T ,
where adj M f (ν) and det M f (ν) are respectively is the adjoint matrix and the determinant of the matrix M f (ν).

From equation ( 3), one has:

det M f (ν) X(ν) = (adj M f (ν)) T Y(ν). (5)
According to the particular structure of the spectral mixing matrix M f (ν), its determinant (denoted by H d (ν)) is a sum of weighted complex exponentials:

H d (ν) = det M f (ν) = i β i exp(-2πντ di ),
where the weights β i are particular products of coefficients c ij , and τ di are particular sums of delays τ ij .

Let H(ν) denote the transposed adjoint of M(ν). Each entries H kl (ν) of H(ν) can also be expressed as a sum of weighted exponentials:

H kl (ν) = (adj M f (ν)) T kl = i α kli exp(-2πντ kli ),
where the weights α kli are also particular products of coefficients c ij , and τ kli are particular sums of delays τ ij .

Back to the temporal representation, the system (5) becomes:

{h d * x k }(t) = N l=1 {h kl * y l }(t),
where Each filtered source {h d * x k }(t) is then a finite sum of delayed observations:

h kl (t) = FT -1 (H kl (ν)), h d (t) = FT -1 (H d (ν))
{h d * x k }(t) = l,i α kli y l (t -τ kli ). (6)
Assuming that the τ kli are "small" for all indices, y l (t -τ kli ) can be approximated by its truncated up to P-order Taylor series expansion. Each filtered source can be then approximated by a linear combination of the observations and its derivatives:

{h d * x k }(t) ≈ l,i α kli P p=0 (-τ kli ) p p! y (p) l (t).
For k = 1, . . . , N the system is rewritten as:

x h (t) ≈ Sỹ(t), (7)
where N ] T is the observations and their derivatives vector. S is the N × (N P + N ) separating matrix.

x h (t) = [x h 1 (t), . . . , x h N (t)]
The analytical expression of each entry of S with respect to the parameters c ij and τ ij is heavy but does not present theoretical difficulties. This expression is detailed in the Appendix A for the 2-sources/2-sensors case.

The joint estimation of the parameters c ij ,τ ij and of the separating matrix S forms the basis of the iterative algorithm proposed in the following section.

PRESENTATION OF THE ALGORITHM

From [START_REF] Dégerine | Sur la diagonalisation conjointe approchée par un critère des moindres carrés[END_REF] and because the sources x h k (t) are mutually uncorrelated, it follows that S has to diagonalize the set of covariance matrices obtained at different lags R ỹỹ (τ k ) = E{ỹ(t)ỹ T (t + τ k )}:

SR ỹỹ (τ k )S T = Λ(τ k ), ∀τ k Λ(τ k ) = E{x h (t)x h T (t+τ k )} being a diagonal matrix what- ever the lag τ k .
S. Dégerine propose in [START_REF] Dégerine | Sur la diagonalisation conjointe approchée par un critère des moindres carrés[END_REF], [START_REF] Dégerine | A new non-orthogonal approximate joint diagonalization algorithm for blind source separation[END_REF] a new algorithm (called Least Square on B or LSB) for approximate non-orthogonal joint diagonalization of a set of matrices. In its original form, this algorithm iteratively searches an optimal joint diagonalizer S of a set of K matrices to the mean square sense. The secondorder criterion to be optimized is then

C( S) = k Off( SR ỹỹ (τ k ) ST ), (8)
where Off(.) is the sum of the square non diagonal elements of the considered matrix. To perform the mean square optimization, S. Dégerine uses a relaxation on the lines of S conducting to solve N eigenvalue problems for each iteration.

Here, the expected diagonalizer S being non square, the optimization [START_REF] Dégerine | A new non-orthogonal approximate joint diagonalization algorithm for blind source separation[END_REF] does not conduct to an unique solution for S. The main idea is then to force the algorithm to hold the theoretical structure of the separation matrix at each iteration. This theoretical expression depends on the parameters c ij and τ ij . At each step of relaxation, we estimate the parameters c ij and τ ij from their formal expression wrt the entries of the current estimated S. The next step is processed using an new S built on these estimated parameters.

An illustration of such theoretical expressions can be found in Appendices for the 2-sources case. For the N -sources case, we implemented formal calculus subroutine in order to provide automatically the formal expressions we need. The final algorithm is summed up as follows:

Numerical Algorithm 1) From the observations y(t), build the observation vector ỹ(t) at the order P .

2) Compute K covariance matrix R ỹỹ (τ k ).

3) Intitialize S(0) with arbitrary values.

4) Proceed to an iteration k with LSB algorithm.

5) From S(k) obtained at step 4), estimate the attenuation coefficients c k ij and the propagation delays τ k ij .

6) From c k ij and τ k ij , compute a new matrix S(k + 1).

7) Repeat steps 4), 5) and 6) until the evolution of the estimated matrix S becomes sufficiently weak.

RESULTS

We present now a numerical simulation in the 3-sources/3observations case. The sources are 24s long music extracts. The sampling frequency is F s = 22.05kHz.

The matrix C of relative attenuations is:

C =  
1.0000 1.0000 1.0000 2.0000 1.3000 1.2000 3.0000 1.6000 2.2000

  .
The matrix D of relative propagation delays is:

D =   0
0 0 30.00µs 18.00µs 5.00µs 25.00µs 12.50µs 10.00µs   .

For the 3 observations, we take a Taylor series expansion up to 3-order (P = 3). 151 covariance matrices R ỹỹ (τ k ) are used (K = 151), with τ k = k/F s, k ∈ {-75, 75}. Taking the matrices C and D into account, the filter H d (ν) can be identified and its effect reversed after separation. We obtain the following source estimations (Fig. 2) that have to be compared to the original ones (Fig. 3).

CONCLUSION

This paper present a novel iterative method for blind source separation from delayed mixtures. The delayed mixtures are approached by instantaneous mixtures of the observation derivatives. This approach uses only second-order statistics and needs no more mixtures than sources. The algorithm is based on the joint diagonalization of a set of spatial covariance matrices but the non square separating matrix is constrained to hold a theoretical structure at each iteration. A numerical simulation shows the efficiency of the approach.

APPENDIX A FROM M f (ν) TO S IN THE 2-OBSERVATIONS CASE. Now, we have the following spectral mixing matrix

M f (ν) = 1 1 c 21 e -2πντ21 c 22 e -2πντ22 . (9)
In this case (6) leads to the 2 following filtered sources: The inverse matrix of S f (ν) is also a matrix of polynomials:

x h 1 (t) = c 22 y 1 (t -τ 22 ) -y 2 (t), x h 2 (t) = -c 21 y 1 (t -τ 21 ) + y 2 (t), where x h k (t) = x k (t -τ 22 ) -x k (t -τ 21 ), k ∈ [1, 2]. With x h (t) = [x h 1 , x h 2 ] T ,
S f (ν) -1 = 1 det S f (ν)
p 11 (ν) p 12 (ν) p 21 (ν) p 22 (ν) .

An approximation of the spectral mixing matrix M f (ν) is obtained normalizing each column of the previous matrix by p 1k (ν)

M f (ν) ≈ 1 1 p 21 (ν)/p 11 (ν) p 22 (ν)/p 12 (ν) .

The entries of the second line of this matrix can be developed in power series expansion in ν:

p 21 (ν)/p 11 (ν) = q 0 + q 1 ν + O 2 (ν), p 22 (ν)/p 12 (ν) = r 0 + r 1 ν + O 2 (ν).

Identifying these expansions with the Taylor-series expansion of the entries of M f (ν) (see 9) we find:

c 21 = q 0 , τ 21 = -1 j2π q 1 q 0 , c 22 = r 0 , τ 22 = -1 j2π r 1 r 0 .

  and where {h . * x . }(t) = h . (t -τ )x . (τ )dτ denotes the convolution of h . (t) and x . (t).

Fig. 1 (Fig. 1 .

 11 Fig. 1 (a) shows the evolution of the optimization criterion C( S) at each iteration number i and part (b) presents the Frobenius norm of the matrix S(k + 1) -S(k)

Fig. 2 .Fig. 3 .

 23 Fig. 2. Estimated sources

Table 1 .

 1

		source #1 source #2 source #3
	MSE	0.0141	0.0155	0.0164
	SIR	41 dB	44 dB	36 dB

Table 1 .

 1 Performances

  and ỹ(t) = [y 1 , ẏ1 , ÿ1 , . . . , y TO c ij , τ ij IN THE 2-OBSERVATIONS CASE. Here, we have the following separating matrix S = s 11 s 12 . . . s 1P s 1P +1 s 21 s 22 . . . s 2P s 2P +1 , one has the (2 × 2) spectral separating matrix

					(P ) 1 , y 2 ] T ,
	the (2 × P + 2) unmixing matrix of equation (7) becomes:
	S =	c 22 -c 22 τ 22 . . . -c 21 c 22 τ 21 . . . -c 21 c 22 (-τ22) P P ! (-τ21) P P !	-1 1	.
			APPENDIX B
	FROM S S f (ν) =	s f 11 (ν) s f 12 (ν) s f 21 (ν) s f 22 (ν)	,
	where the entries of S f (ν) are polynomials in ν:
		s f 11 (ν) =	P k=1 (2jπν) k s 1k ,
		s f 12 (ν) =	s 1P +1 ,
		s f 21 (ν) =	P k=1 (2jπν) k s 2k ,
		s f 22 (ν) =	s 2P +1 .


M Baidada

K Mansouri

F Poirier

Vers un modèle d'utilisation de la traçabilité dans les plateformes e-learning, en vue de la prise de décision dans le recrutement ou la promotion des ressources humaines

Résumé

Le domaine de la gestion des ressources humaines (GRH) est un contexte qui accentue de plus en plus l'utilisation des nouvelles technologies de l'information pour assister les responsables dans la prise de décision lors d'un recrutement d'embauche ou de promotion. Un CV d'un candidat ou un dossier d'un employé peut prendre un format électronique, mais peuvent-ils apporter plus de précision par rapport aux compétences des candidats, et épargner les responsables les efforts supplémentaires des tests techniques ? Il est alors fortement envisageable d'utiliser ces outils technologiques, notamment les systèmes d'apprentissage en ligne, pour avoir une traçabilité des compétences acquises, et donner une visibilité rapide aux responsables. Cet article sera présenté comme suit : nous commencerons par poser le problème, nous poursuivrons par faire une présentation de l'état de l'art dans le domaine, nous présenterons ensuite notre approche proposée, et nous terminerons avec une conclusion et perspectives.

Mots clé

Position du problème

Le module de gestion des ressources humaines est certainement un des modules essentiels composant un système d'information [1]. Les responsables des ressources humaines utilisent toutes les possibilités offertes par les systèmes informatisés pour être plus efficaces :

-

Etat de l'art dans le domaine

Les solutions e-learning peuvent s'intégrer à d'autres systèmes d'information : scolarité, ressources humaines, etc., et peuvent être représentées par le schéma global de la figure 1: -Les outils de reporting : pour générer toute information traduisant une traçabilité de tout ce qui se passe dans le système [START_REF] Alowayr | Review of monitoring tools for e-learning platforms[END_REF] ; -Les e-portfolio : pour permettre aux apprenants de mettre à jour, sous forme de CV résumés, leurs savoirs acquis. Plusieurs solutions existent telles que MAHARA [START_REF]Mahara user manual[END_REF]

Approche proposée

Notre solution consiste à proposer un modèle, permettant d'intégrer un système de gestion des ressources humaines avec une plateforme e-learning, tout en gardant les deux composantes : outil de reporting et système e-portfolio. On peut envisager alors d'exploiter des données analysées dans l'outil de reporting pour automatiser la mise à jour des systèmes e-portfolio ; ainsi les responsables trouveront à leur portée des données pertinentes sur tout candidat qui est passé par la plateforme d'apprentissage de e-learning. Le système e-learning, comme il a été présenté ci haut, génére des logs contenant toutes les informations traduisant le déroulement de l'apprentissage dans le système, et des évaluations qui les accompagnent. Certaines recherches se sont intéréssées à l'extraction et au traitement de ces informations [START_REF] Delgado Calvo-Flores | Predicting students' marks from Moodle logs using neural network models[END_REF]. Comme taces dans les logs, on peut trouver : le nombre d'accès à une ressource en ligne, le nombre d'essais pour chaque réponse lors d'un test, le temps passé dans chaque activité, le taux de réussite [START_REF] Settouti | Vers des Systèmes à Base de Traces modélisées pour les EIAH[END_REF].

Nous proposons, dans notre modèle, d'automatiser la mise à jour du e-portfolio d'un apprenant en exploitant les indicateurs receuillis dans les outils de reporting, et offrir ainsi une traçabilité pour les responsables ressources humaines, qui les aidera à la prise de décision lors d'un recrutement ou d'une promotion.

On peut donner comme exemple, un poste pour lequel un certain niveau d'anglais est requis, les responsables auront alors devant eux une traçabilité de l'évolution du niveau d'un candidat en langue anglaise, à travers ce qui aété recceuilli dans les logs lors des formations suivi par le candidats, et ce qui a été automatiquement inscrit dans son e-protfolio.

On peut envisager d'adapter le schéma proposé avec une multitude de plateformes e-learning, si une interopérabilité est assurée entre celles-ci. Ainsi, les informations sur un apprenant qui est passé par une plateforme précédente, se trouveront facilement exploitables par la plateforme utilisée par les responsables en ressources humaines, et assurer ainsi une continuité et une transparence vis-à-vis de l'usage de différentes plateformes.

Conclusion

La

Figure 1 : 3 ]

 13 Figure 1 : Schéma global d'une solution e-learning [3] Chaque composante présente un intérêt particulier, nous nous intéresserons surtout aux outils d'analyse et de reporting et aux systèmes e-portfolio.

Figure 2 :

 2 Figure 2 : Schéma de l'apprcohe proposée Des archictures de médiation telle que l'architecture SOA, pour intéropérer le système e-learning avec le système d'information des ressources humaines [8], et assurer ainsi l'échange d'informations telles que, les sessions de formations planifiées, les tests en ligne réalisés ou les résultats obtenus.

  Cette traçabilité peut être assurée par une exploitation des outils de reporting, et une mise à jour automatique des e-portfolio des candidats.Comme perspective, nous pouvons conceptualiser l'approche, et proposer un modèle précisant les choix des technologies à utiliser. Information des Ressources Humaines (SIRH) : un atout dans l'optimisation de la GRH au service de l'entreprise", p51 à p55, Mémoire MBA Management des Ressources Humaines, Univ. Paris Dauphine, Octobre 2010 [2] : S.Gottrand, V.Queant, "Le e-learning comme innovation en Ressources Humaines", dossier sous la direction de P.Louart, 2003

	prise de décision rapide et efficace d'un responsable
	en ressources humaines pour un recrutement ou une
	promotion est un défi, qui peut se concrétiser par une
	meilleure utilisation des nouvelles technologies,
	notamment les plateformes e-learning avec tout ce
	qu'elles peuvent offrir comme traçabilité des acquis
	des différents candidats, en donnant avec précision les
	niveaux d'aptitude et d'habilité pour tout candidat.

Références

[1] : G.Exbrayat, N.Fisteberg, R.Fouesnant, "Le Système d'


Lluís Antoni

Jiménez Rugama

¨laurent Gilquin

Laurent Gilquin
email: gilquin.laurent@inria.fr

Reliable error estimation for Sobol' indices

Keywords:

In the field of sensitivity analysis, Sobol' indices are sensitivity measures widely used to assess the importance of inputs of a model to its output. The estimation of these indices is often performed trough Monte Carlo or quasi-Monte Carlo methods. A notable method is the replication procedure that estimates first-order indices at a reduced cost in terms of number of model evaluations.

An inherent practical problem of this estimation is how to quantify the number of model evaluations needed to ensure that estimates satisfy a desired error tolerance. This paper addresses this challenge by proposing a reliable error bound for first-order and total effect Sobol' indices. Starting from the integral formula of the indices, the error bound is defined in terms of the discrete Walsh coefficients of the different integrands.

We propose a sequential estimation procedure of Sobol' indices using the error bound as a stopping criterion. The sequential procedure combines Sobol' sequences with either Saltelli's strategy to estimate both first-order and total effect indices, or the replication procedure to estimate only firstorder indices.

Mathematics Subject Classification (2000) 49Q12

62L12 ¨65R10

Introduction

Let f represent a deterministic numerical model in r0, 1s d , d ě 1. Sensitivity measures, also known as Sobol' indices, are used to assess which inputs of f are influential for the output. The normalized indices are scalars between 0 and 1 whose values are interpreted as follows: the closer to 1 the more influential the index. Alternatively, they can be interpreted as the percentage of the variance explained by the inputs. Among all Sobol' indices one can distinguish firstorder and total effect indices. The first measure the effect of a single input, while the second measure the effect of a single input plus all its interactions with other inputs. When dealing with complex numerical models, analytical expressions of Sobol' indices are often inaccessible. In such cases, one must rely on an estimation of these indices. The original estimation procedure is due to Sobol' [START_REF] Sobol | Sensitivity indices for nonlinear mathematical models[END_REF]. However, this procedure requires several model evaluations which are usually expensive. Later on, Saltelli [START_REF] Saltelli | Making best use of models evaluations to compute sensitivity indices[END_REF] proposed strategies to estimate sets of Sobol' indices at once through the use of a combinatorial formalisms. While elegant, these strategies still require a large number of model evaluations. A cost efficient alternative to estimate first-order indices was introduced in [START_REF] Mara | Comparison of some efficient methods to evaluate the main effect of computer model factors[END_REF]. This alternative, called the replication procedure, has recently been further studied in [START_REF] Tissot | A randomized Orthogonal Array-based procedure for the estimation of first-and second-order Sobol' indices[END_REF] and generalized to the estimation of closed second-order indices.

A practical problem concerning the use of these methods is how to quantify the number of model evaluations required to ensure that Sobol' estimates are accurate enough. This article addresses this challenge by proposing a reliable error bound for Sobol' indices based on digital sequences. The error bound is defined in terms of the discrete Walsh coefficients of the integrands involved in the Sobol' indices definition. We propose a sequential estimation procedure of Sobol' indices using the error bound as our stopping criterion. The procedure operates under the assumption that all integrands lie inside a particular cone of functions (see [START_REF] Hickernell | Reliable Adaptative Cubature Using Digital Sequences: Monte Carlo and Quasi-Monte Carlo Methods[END_REF]).

Firstly, Section 2 introduces Sobol' indices and reviews both Saltelli's strategy to estimate first-order and total effect Sobol' indices, and the replication procedure. Our main contribution is detailed in Section 3. There, we review the construction of the error bound proposed in [START_REF] Hickernell | Reliable Adaptative Cubature Using Digital Sequences: Monte Carlo and Quasi-Monte Carlo Methods[END_REF] for the estimation of integrals based on digital sequences, and then we generalize it for Sobol' indices. Section 4 is devoted to analyze the cost in terms of model evaluations of our sequential estimation algorithm. It combines the error bound in Section 3 and either one of the two estimation procedures of Section 2. We also discuss a potential improvement to estimate small first-order indices according to [START_REF] Owen | Better estimation of small Sobol' sensitivity indices[END_REF]. Finally, examples and illustrations of our procedure are presented in Section 5.

Background on Sobol' indices

Definition of Sobol' indices

Denote by x " px 1 , . . . , x d q the vector of inputs of f . We assume that f is in some subset of P L 2 pr0, 1s d q for which f pxq is defined for all x P r0, 1s d , and D " t1, . . . , du the set of dimension indexes. Let u be a subset of D, ´u its complement and |u| its cardinality. Then, x u represents a point in r0, 1s |u| with components x j , j P u. Given two points x and x 1 , the hybrid point w " px u : x 1 ´uq is defined as w j " x j if j P u and w j " x 1 j if j R u. The uncertainty on x is modeled by a uniform random vector, namely x iid " U pr0, 1s d q. The Hoeffding decomposition [START_REF] Hoeffding | A class of statistics with asymptotically normal distribution[END_REF][START_REF] Sobol | Sensitivity indices for nonlinear mathematical models[END_REF] of f is:

f pxq " f ∅ `ÿ uĎD,u‰∅ f u pxq, (1)
where:

f ∅ " Er f pxqs " µ, f u pxq " ż r0,1s |u| f pxqdx ´u ´ÿ vĂu f v pxq.
Due to orthogonality, the variance of equation ( 1) leads to the variance decomposition of f :

σ 2 " Varr f pxqs " ÿ uĎD,u‰∅ σ 2 u , with σ 2 u " ż r0,1s |u| f u pxq 2 dx u .
From this decomposition, one can define the following two quantities:

τ 2 u " ÿ vĎu σ 2 v , τ 2 u " ÿ vXu‰∅ σ 2 v , u Ĺ D.
These two quantities τ 2 u and τ 2 u measure the importance of variables x u : τ 2 u quantifies the main effect of x u , that is the effect of all interactions between variables in x u , and τ 2 u quantifies the main effect of x u plus the effect of all interactions between variables in x u and variables in x ´u.

Both τ 2 u and τ 2 u satisfy the following relations: 0 ď τ 2 u ď τ 2 u and τ 2 u " σ 2 ´τ2 ´u. These two measures are commonly found in the literature in their normalized form: S u " τ 2 u {σ 2 is the closed |u|-order Sobol' index for inputs u, while S u " τ 2 u {σ 2 is the total effect Sobol' index of order |u|.

In our framework, we are only interested in single input indices, namely |u| " 1. The computation of the normalized indices is performed based on the following integral formulas for their numerators:

τ 2 u " ż r0,1s 2d `f pxq ´f px u : x 1 ´uq ˘f px 1 qdxdx 1 , (2)
τ 2 u " 1 2 ż r0,1s d`1 p f px 1 q ´f px u : x 1 ´uqq 2 dx u dx 1 , u P D, (3)
while variance and mean of f are evaluated as:

σ 2 " ż r0,1s d f pxq 2 dx ´µ2 , µ " ż r0,1s d f pxqdx, (4)
Usually the complexity of f causes the solution of integrals (2), ( 3) and ( 4) to be intractable. In such cases, one can instead estimate these quantities.

Estimation of Sobol' indices

In this section we review two Monte Carlo procedures for the estimation of Sobol' indices. A design is a point set P " tx i u n´1 i"0 where each point is obtained by sampling each variable x j n times. Each row of the design is a point x i in r0, 1s d and each column of the design refers to samples of a variable x j . Consider P " tx i u n´1

i"0 and P 1 " tx 1 i u n´1 i"0 two designs where px i , x 1 i q

iid " U pr0, 1s 2d q. One way to estimate the two quantities (2) and ( 3) is via:

p τ 2 u " 1 n n´1 ÿ i"0 `f px i q ´f px i,u : x 1 i,´u q ˘f px 1 i q, (5)
p τ 2 u " 1 2n n´1 ÿ i"0 p f px 1 i q ´f px i,u : x 1 i,´u qq 2 , u P D, (6)
using for σ 2 :

p σ 2 " 1 n n´1 ÿ i"0 f px i q 2 ´p µ 2 , with p µ " 1 n n´1 ÿ i"0 f px i q. (7)
Then, the Sobol' indices estimators are:

p S u " p τ 2 u { p σ 2 , p S u " p τ 2 u { p σ 2 . ( 8
)
Using this method, the estimation of a single pair ( p S u , p S u ) requires 3n evaluations of the model f . Thus, for all first order and total effect indices, one would need 3nd model evaluations. Using a combinatorial formalism, in [START_REF] Saltelli | Making best use of models evaluations to compute sensitivity indices[END_REF] Saltelli proposes the following estimation strategy:

Theorem 1 The d `2 designs tx i,u : x 1 i,´u u n´1 i"0 constructed for u P t∅, t1u, . . . , tdu, Du allow to estimate all first-order and total effect Sobol' indices at a cost of npd `2q evaluations of the model.

The main idea of the theorem is that there is no need to reevaluate f px i q and f px 1 i q for each u. Hence, we can simply evaluate f px i q n times, f px 1 i q n times, and f px i,u : x 1 i,´u q nd times, which accounts for the npd `2q evaluations.

The d `2 designs of Theorem 1 are obtained by substituting columns of P for columns of P 1 accordingly to u. While elegant, this approach requires a number of model evaluations that grows linearly with respect to the input space dimension.

A more efficient alternative to evaluate all first-order indices was proposed in [START_REF] Mara | Comparison of some efficient methods to evaluate the main effect of computer model factors[END_REF] and only requires 2n model evaluations. This alternative relies on the construction of two replicated designs. The notion of replicated designs was first introduced by McKay through his replicated Latin Hypercubes in [START_REF] Mckay | Evaluating prediction uncertainty[END_REF]. In order to apply this definition to other types of points, in [START_REF] Gilquin | Iterative construction of replicated designs based on Sobol' sequences[END_REF] we generalized this notion to the following. Definition 1 Let P " tx i u n´1

i"0 and P 1 " tx 1 i u n´1 i"0 be two point sets in r0, 1s d . Let P u " tx i,u u n´1 i"0 (resp. P 1 u ), u Ĺ D, denote the subset of dimensions of P (resp. P 1 ) indexed by u. We say that P and P 1 are two replicated designs of order a P t1, . . . , d ´1u if for any u Ĺ D such that |u| " a, P u and P 1 u are the same point set in r0, 1s a . We define by π u the permutation that rearranges the rows of P 1 u into P u .

The method introduced in [START_REF] Mara | Comparison of some efficient methods to evaluate the main effect of computer model factors[END_REF] allows to estimate all firstorder Sobol' indices with only two replicated designs of order 1. The key point of this method is to use the permutations resulting from the structure of the two replicated designs to mimic the hybrid points in equation [START_REF] Hong | Algorithm 823: Implementing scrambled digital nets[END_REF].

More precisely, let P " tx i u n´1 i"0 and P 1 " tx 1 i u n´1 i"0 be two replicated designs of order 1. Denote by ty i u n´1 i"0 " t f px i qu n´1 i"0 and ty 1 i u n´1 i"0 " t f px 1 i qu n´1 i"0 the two sets of model evaluations obtained with P and P 1 . From Definition 1, we know that x 1 π u piq,u " x i,u . Then, y 1 π u piq " f px 1 π u piq,u : x 1 π u piq,´u q " f px i,u : x 1 π u piq,´u q.

Hence, each τ 2 u can be estimated via formula (5) by using y 1 π u piq instead of f px i,u : x 1 i,´u q without requiring further model evaluations for each u. This estimation method has been studied deeply and generalized in Tissot et al. [START_REF] Tissot | A randomized Orthogonal Array-based procedure for the estimation of first-and second-order Sobol' indices[END_REF] to the case of closed second-order indices. In the following we will refer to this method as the replication procedure.

Towards a reliable estimation

The aim of this article is to propose a sequential procedure to estimate first-order and total effect Sobol' indices. A practical problem concerning the estimation of these indices is how large to choose the number of evaluations to ensure that Sobol' estimates are accurate enough. Asymptotical results show that Sobol' estimates are normally distributed ([7, Proposition 2.2], [START_REF] Tissot | A randomized Orthogonal Array-based procedure for the estimation of first-and second-order Sobol' indices[END_REF]Proposition 3.5]). As a consequence, errors can be estimated through confidence intervals. However, these error estimates are only guaranteed asymptotically as the number of model evaluations goes to infinity. Additional sequential procedures are the replicated procedure and McKay's procedure respectively proposed in [START_REF] Gilquin | Recursive estimation procedure of Sobol' indices based on replicated designs[END_REF] and [START_REF] Tong | Self-validated variance-based methods for sensitivity analysis of model outputs[END_REF]. Nevertheless, in those two cases, the stopping criterion is a purely empirical quantity of interest, built directly upon the estimates. Such stopping criteria often involve hyper-parameters that are difficult to tweak but more importantly, fail to guarantee any error bound on the estimates.

Our sequential procedure stands apart from others since it proposes a robust stopping criterion, not costly to compute. This criterion is an error bound based on the Walsh series decomposition of the integrands in (2), ( 3) and ( 4), and exploits the group properties of digital nets. As such, our procedure relies on an iterative construction of Sobol' sequences. This construction is performed accordingly to the multiplicative approach presented in [START_REF] Gilquin | Iterative construction of replicated designs based on Sobol' sequences[END_REF].

The description of the error bound is introduced in the following section and our sequential procedure is detailed in Section 4.

Reliable error bound for Sobol' indices

We start by reviewing the construction of the error bound proposed in [START_REF] Hickernell | Reliable Adaptative Cubature Using Digital Sequences: Monte Carlo and Quasi-Monte Carlo Methods[END_REF] for the estimation of d-dimensional integrals. Then, we present an extension of this error bound for normalized Sobol' indices. This extension is built upon the integral formula of a Sobol' index.

Reliable integral estimation using digital sequences

We assume we have an embedded sequence of digital nets in base b as in [4, Sec. 2-3],

P 0 " t0u Ă ¨¨¨Ă P m " tx i u b m ´1
i"0 Ă ¨¨¨Ă P 8 " tx i u 8 i"0 .

Each P m has a group structure under the digitwise addition:

x ' t " ˜8 ÿ "1 rpx j `t j q mod bsb ´ pmod 1q ¸dj"1 ,
where x j and t j are the b-adic decompositions of the j th component of points x and t.

To relate the group structure of P m with the integration error, we introduce the dual net which establishes the relationship between any digital net and the wavenumber space of non-negative integers N d 0 . A dual net is

P K m " tk P N d 0 : xk, xy " 0, x P P m u, xk, xy " d ÿ j"1 8 ÿ "0 k j x j, `1 pmod bq,
and inherits the same embedded structure as for digital nets,

P K 0 " N d 0 Ą ¨¨¨Ą P K 8 " t0u. (9)
As shown in [START_REF] Hickernell | Reliable Adaptative Cubature Using Digital Sequences: Monte Carlo and Quasi-Monte Carlo Methods[END_REF]Sec. 3], the group structure of digital nets guarantees the property below affecting any Walsh basis

ϕ k pxq " e 2π ? ´1xk,xy{b , 1 b m ÿ xPP m ϕ k pxq " # 1, k P P K m , 0, k R P K m . (10)
Therefore, considering the Walsh decomposition of any f P L 2 pr0, 1s d q:

f pxq " ÿ kPN d 0 fk ϕ k pxq,
I ´1 b m ÿ xPP m f pxq " f0 ´1 b m ÿ xPP m ÿ kPN d 0 fk ϕ k pxq , " ÿ kPP K m zt0u fk ď ÿ kPP K m zt0u fk . (11)
Based on the size of fk and the structure of dual nets [START_REF] Lemieux | Monte Carlo and Quasi-Monte Carlo Sampling[END_REF], in [START_REF] Hickernell | Reliable Adaptative Cubature Using Digital Sequences: Monte Carlo and Quasi-Monte Carlo Methods[END_REF]Sec. 4.1] we proposed an ordering of the wavenumbers kp¨q : N 0 Ñ N d 0 . Defining fκ " fkpκq , the error bound (11) becomes,

I ´1 b m ÿ xPP m f pxq ď 8 ÿ λ "1 fλb m . ( 12
)
Nonetheless, because the knowledge of the Walsh coefficients fκ is not assumed, we will estimate them using the fast transform obtained with the precomputed function values, and refer to them as r f m,κ ,

r f m,κ " r f m,kpκq " 1 b m ÿ xPP m f pxqϕ kpκq pxq.
Note that for function values evaluated at P m , one only generates b m discrete Walsh coefficients r f m,κ . For , m P N 0 and ď m we introduce the following notation,

S m p f q " b m ´1 ÿ κ"tb m´1 u fκ , p S ,m p f q " b ´1 ÿ κ"tb ´1u 8 ÿ λ "1 fκ`λb m , q S m p f q " p S 0,m p f q `¨¨¨`p S m,m p f q " 8 ÿ κ"b m fκ , r S ,m p f q " b ´1 ÿ κ"tb ´1u r f m,κ .
Finally, we define the set of functions C ,

C :" t f P L 2 pr0, 1s d q : p S ,m p f q ď p ωpm ´ q q S m p f q, ď m, q S m p f q ď ωpm ´ qS p f q, ˚ď ď mu. ( 13
)
for ˚P N, p ω and ω two non-negative valued functions with lim mÑ8 ωpmq " 0, and such that p ωprq ωprq ă 1 for some fixed r P N.

In [4, Sec. 4.2] it is shown that for any f P C , [START_REF] Owen | Better estimation of small Sobol' sensitivity indices[END_REF] where one may increase m until the error bound ε I is small enough.

I ´1 b m ÿ xPP m f pxq ď r S ,m p f q p ωpmq ωpm ´ q 1 ´p ωpm ´ q ωpm ´ q looooooooooooomooooooooooooon Cpmq " ε p I ,
Details concerning the algorithm, the mapping of the wavenumber space, or the meaning and properties of C , are provided in [START_REF] Hickernell | Reliable Adaptative Cubature Using Digital Sequences: Monte Carlo and Quasi-Monte Carlo Methods[END_REF].

For our problem, only Sobol' sequences [START_REF] Sobol | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF] have been considered. Their major interest comes from their fast and easy implementation. Further details concerning Sobol' sequences can be found in [START_REF] Lemieux | Monte Carlo and Quasi-Monte Carlo Sampling[END_REF][START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF].

Extension to Sobol' indices

The idea here is to extend the definition of the error bound [START_REF] Owen | Better estimation of small Sobol' sensitivity indices[END_REF] to Sobol' indices. To achieve this goal, we consider the two integral formulas of the first-order and total effect Sobol' indices with |u| " 1:

S u pIq " ş r0,1q 2d p f pxq ´f px u : x 1 ´uqq f px 1 qdxdx 1 ş r0,1q d f pxq 2 dx ´´ş r0,1q d f pxqdx ¯2 , S u pIq " I 1 I 3 ´pI 4 q 2 , S u pIq " 1 2 ş r0,1q d`1 p f px 1 q ´f px u : x 1 ´uqq 2 dx u dx 1 ş r0,1q d f pxq 2 dx ´´ş r0,1q d f pxqdx ¯2 , S u pIq " I 2 I 3 ´pI 4 q 2 ,
where I " pI 1 , I 2 , I 3 , I 4 q is a vector of integral values. Indices S u pIq and S u pIq are defined as functions over the vector I.

Because numerator and denominator are both known to be positive, maximizing S u pIq (resp. S u pIq) is done through maximizing the numerator I 1 (resp. I 2 ) and minimizing the denominator I 3 ´I2 4 . Analogously, to minimize S u pIq (resp. S u pIq) one minimizes the numerator I 1 (resp. I 2 ) and maximizes the denominator I 3 ´I2 4 . As an example, Figure 1 illustrates the region of possible values of S 2 pIq given the true value of I 1 for the test function of Bratley et al. [START_REF] Bratley | Implementation and tests of low-discrepancy sequences[END_REF]. This function is further described in Section 5.

Under the assumption that each integrand in S u and S u lies inside C , these new estimators satisfy:

S u P " p S u ´εp S u , p S u `εp S u ı , S u P " p S u ´εp S u , p S u `εp S u ı .
For the rest of the paper, we will consider p S u and p S u as defined in formulas ( 15) and ( 16). We start by detailing our procedure under the form of an algorithm. Then, we discuss a possible improvement by considering a new estimator introduced in [START_REF] Owen | Better estimation of small Sobol' sensitivity indices[END_REF] for the estimation of first-order indices.

Sequential algorithm and cost

Algorithm 1 summarizes the main steps of our sequential procedure. First, one must fix the tolerance ε ą 0 at which Sobol' indices must be estimated, and set m " ˚`r, where ˚and r are the two hyper-parameters defined in [START_REF] Owen | Randomly permuted pt, m, sq-nets and pt, sqsequences: Monte Carlo and Quasi-Monte Carlo Methods in[END_REF]. Then, one must construct the two designs P m " tx i u 2 m ´1

i"0 and P 1 m " tx 1 i u 2 m ´1 i"0 at Step 5 according to the multiplicative approach detailed in [START_REF] Gilquin | Iterative construction of replicated designs based on Sobol' sequences[END_REF]. The choice of ˚must be large enough to ensure that each integrand of S u and S u (Section 3.2) is in C (see [START_REF] Owen | Randomly permuted pt, m, sq-nets and pt, sqsequences: Monte Carlo and Quasi-Monte Carlo Methods in[END_REF]).

Under this construction, P m and P 1 m correspond to the first 2 m points of two independent Sobol' sequences. Furthermore, they possess the structure of two replicated designs of order 1, as shown in [START_REF] Gilquin | Iterative construction of replicated designs based on Sobol' sequences[END_REF]. In addition, P m and P 1 m can also be scrambled as long as both point sets share the same scrambling dimensionwise.

The sets P m and P 1 m can be used with Saltelli's strategy to estimate all first-order indices and total effect Sobol' indices. This option is referred as Variant A in Algorithm 1.

Alternatively, P m and P 1 m can be used with the replication procedure to estimate all first-order Sobol' indices. This option is referred as Variant B in Algorithm 1.

In both cases we always check if the respective error bounds ε The cost of our algorithm, in terms of model evaluations, varies whether Variant A or Variant B is selected. To discuss this cost we note by m ‹ the ending iteration. If Variant A is selected, the cost of our algorithm is:

ÿ uPD 2 m u `2 ˆ2m ‹ , m ‹ " max uPD m u ,
where:

' 2 m u is the number of evaluations f px i,u : x 1 i,´u q used to estimate both the first-order index S u and the total effect index S u , ' 2 ˆ2m ‹ is the number of evaluations f px i q and f px 1 i q used in the estimation of each first-order and total effect indices.

If all m u are equal, the cost of Variant A becomes 2 m ‹ pd `2q and we recover the cost specified in Theorem 1 with n " 2 m ‹ .

If Variant B is selected, the cost of our algorithm equals 2 ˆ2m ‹ . This cost corresponds to the one of the replication procedure introduced in Section 2.2, where 2n " 2 ˆ2m ‹ independent of d.

Improvement

When S u is small, it usually becomes harder to estimate. For this reason, we consider the use of a new estimator to evaluate small first-order Sobol' indices in Variant A. This estimator called "Correlation 2" has been introduced by Owen in [START_REF] Owen | Better estimation of small Sobol' sensitivity indices[END_REF]. In this article, he discussed and highlighted the efficiency of "Correlation 2" when estimating small first-order indices. Our aim is to show that the use of "Correlation 2" in Variant A may reduce the total number of model evaluations. This new estimator is,

p τ 2 u " 1 n n´1 ÿ i"0 p f px i q ´f pz i,u : x i,´u qqp f px i,u : x 1 i,´u q ´f px 1 i qq,
where px i , x 1 i , z i q iid " U pr0, 1sq 3d . It requires an additional set of n model evaluations to estimate τ 2 u .

We discuss below the potential improvement brought by using "Correlation 2" in Variant A. The idea is to replace the current estimator (5) by "Correlation 2" for each small first-order index.

Assume that the number of small first-order indices is known and equals γ. We denote by u 1 , . . . , u γ the indices of the corresponding inputs and Γ " t1, . . . , γu. The cost of Variant A including "Correlation 2" becomes,

ÿ jPΓ 2 m 2 u j `ÿ jPΓ 2 m 1 u j `ÿ jPDzΓ 2 m u j `2 ˆ2m ‹ , (17)
where:

' for j P Γ , 2 m 2
u j is the number of evaluations f `zi,u j : x i,´u j ˘to estimate S u j , ' for j P Γ , 2 m 1 u j is the number of evaluations f `xi,u j : x 1 i,´u j ˘to estimate both S u j and S u j , ' likewise, for j P DzΓ , 2 m u j is the number of evaluations f `xi,u j : x 1 i,´u j ˘to estimate both S u j and S u j , ' 2 ˆ2m ‹ is the number of evaluations f px i q and f px 1 i q used in the estimation of each first-order and total effect index.

We recall the cost of Variant

A without "Correlation 2", ÿ jPΓ 2 m u j `ÿ jPDzΓ 2 m u j `2 ˆ2m ‹ . ( 18
)
The difference [START_REF] Sobol | Sensitivity indices for nonlinear mathematical models[END_REF] ´( 18) equals:

ÿ jPΓ 2 m u j ´2m 2 u j ´mu j `2m 1 u j ´mu j ´1¯" ÿ jPΓ c j . (19)
Hence, the sign of this difference indicates whether or not using "Correlation 2" brings an improvement to Variant A.

We distinguish two cases :

1) for j P Γ , if the total effect index S u j requires as much or more evaluations than the first-order index S u j . Since the total effect estimator is the same, as a consequence we have m 1 u j " m u j and c j ą 0. 2) for j P Γ , if the total effect index S u j requires less evaluations than the first-order index S u j . In this case, if both m 2 u j ă m u j and m 1 u j ă m u j then c j ď 0.

Overall we expect to observe case 2q more often than case 1q. Indeed, the numerator of S u requires to estimate 2d dimensional integrals against only d `1 dimensional integrals for the numerator of S u . Hence, it seems reasonable to expect that it will take less points to estimate S u than S u . Furthermore, in case 2q, we expect the two conditions m 2 u j ă m u j and m 1 u j ă m u j to usually hold as in [START_REF] Owen | Better estimation of small Sobol' sensitivity indices[END_REF], "Correlation 2" is shown to perform better for small first-order indices. To support this discussion, numerical examples are presented in Section 5.

In practice, one does not know which are the small Sobol' indices (u 1 , . . . , u γ q. To overcome this issue, we propose the following alternative for Variant A. If at the end of the first iteration p S u is smaller than a threshold (in our case 0.1, as suggested by Owen), then the estimator ( 5) is switched to "Correlation 2" for this particular u and a third Sobol' sequence P 2 m " tz i u 2 m ´1 i"0 is constructed to obtain the corresponding evaluations f pz i,u : x i,´u q.

Applications

We illustrate our sequential estimation procedure with two classical test functions and one real example. In each case, Sobol' indices are estimated with the following three variants:

' Variant A.a as Variant A without using "Correlation 2". ' Variant A.b as Variant A using "Correlation 2". ' Variant B.

The threshold to decide whether or not a first-order index is small is set to 0.1.

For the two test functions, results of the three variants are compared based on the true estimation errors. These errors correspond to the absolute difference between the true values and their estimates:

δ S u " S u ´p S u , δ S u " S u ´p S u .
We also compare the total number of evaluations of each variant. Results are averaged over 100 repetitions using Owen's scrambled Sobol' sequences [START_REF] Hong | Algorithm 823: Implementing scrambled digital nets[END_REF][START_REF] Owen | Randomly permuted pt, m, sq-nets and pt, sqsequences: Monte Carlo and Quasi-Monte Carlo Methods in[END_REF]. For each repetition, we fix the tolerance ε " 5.10 ´3, set ˚" 5, and r " 4 (the algorithm starts with 512 Sobol' points). The algorithm error bound factor described in ( 14) is set to Cpmq " 10 ˆ2´m . The failure rates r δ S u and r δ Su of each estimate are evaluated over the 100 repetitions. These rates are the proportion of the repetitions in which the true estimation errors failed to satisfy δ S u ď ε and δ S u ď ε. When δ S u ą ε and δ S u ą ε, we can infer that our algorithm parameters that define the cone are not conservative enough and the integrands fall outside the cone (13).

Classical test functions

The two classical test functions considered in this article are the g-function introduced by Sobol' [START_REF] Sobol | Sensitivity indices for nonlinear mathematical models[END_REF], and the function introduced by Bratley et al. [START_REF] Bratley | Implementation and tests of low-discrepancy sequences[END_REF]. The idea is to test our method over two categories of functions: additive (Bratley et al. function) and multiplicative (g-function).

Sobol' g-function

The g-function is defined as follows:

f pxq " d ź j"1
g j px j q, g j px j q " |4x j ´2| `a j 1 `a j , a j ě 0.

Each value a j determines the relative importance of the x j .

When the value of a j gets closer to zero the variable x j becomes more influent. For this example, we chose d " 6 and a 1 " 0, a 2 " 0.5, a 3 " 3, a 4 " 9, a 5 " 99, and a 6 " 99. Table 1 shows the averaged estimation errors obtained with Variant A.a as well as the averaged total number of evaluations performed. Table 2 shows the same results obtained with Variant A.b. The main observation is that both approaches give similar results both in terms of estimation errors, failure rates and total number of evaluations. The use of "Correlation 2" in Variant A.b to estimate the four small first-order Sobol' indices S 3 , S 4 , S 5 , S 6 does not seem to bring much improvement. We propose to draw boxplots of the estimation errors to further investigate the results.

Figure 2 shows boxplots of the 100 estimation errors δ S u obtained with both Variant A.a and Variant A.b for the four inputs x 3 , x 4 , x 5 , x 6 . The dashed horizontal line marks the tolerance ε " 5.10 ´3. As expected, we observe that the use of "Correlation 2" in Variant A.b results in lower estimation errors for these four inputs. The discrepancy observed does not stand out in Table 2 due to its low magnitude (10 ´4 " 10 ´5).

Table 3 shows averaged estimation errors δ S u and total number of evaluations for Variant B. The results show that the use of Variant B leads to slightly higher estimation errors and higher failure rates than those obtained with Variant A.a

p´1q i i ź j"1 x j .
The importance of each variable x j depends on their own rank. More explicitly, x 1 is more influent than x 2 that is respectively more influent than x 3 and so on.

As for the g-function, Tables 4 and5 show averaged estimation errors and total number of evaluations for Variant A.a and Variant A.b. Variant A.b gives lower estimation errors and fewer failure rates than Variant A.a on the four small first-indices S 3 , S 4 , S 5 , S 6 which highlights the performance of "Correlation 2". The discrepancy is particularly notable for input x 3 where Variant A.a reaches a failure rate of 71% against only 4% for Variant A.b. Furthermore, Variant A.b requires less evaluations than Variant A.a. The boxplots represented in Figure 3 emphasize the latter observations. Table 6 shows averaged estimation errors δ S u and total number of evaluations for Variant B. The estimation errors are lower than those of Variant A.a but higher than those of Variant A.b. Failure rates are similar to those of Variant A.b. Since the total number of evaluations is close to the one of Variant A.b, the conclusion is that Variant B does not bring much improvement for this example.

From the results on these two test functions, the main conclusion is that Variant A.b performs the best. This highlights the efficiency of "Correlation 2" to estimate small first-order indices with quasi-Monte Carlo methods. Even if "Correlation 2" originally requires more model evaluations, this drawback is completely overridden when this estimator is included into a sequential procedure such as ours. For the case where we are only interested in estimating first-order indices, Variant B should be the best choice as illustrated with the g-function.

Real case model

As a real case example, we will study the payoff function of an arithmetic mean Asian call option. For a discretized Brownian motion B " pB t 1 , B t 2 , . . . , B t d q at times t i " T ˆi{d, the payoff of the option is:

f payoff pBq " e ´ρT max ˜1 d d ÿ i"1
S 0 e pρ´σ 2 {2qt i `σ B t i ´K, 0 where T is the maturity of the option, S 0 is the initial stock price, σ the volatility of the stock, ρ the interest rate, and K the strike price.

The discretized Brownian motion B follows a multivariate normal distribution with mean 0 and covariance matrix Σ, with Σ i j " minpt i ,t j q. If Y " Np0, I d q, then for any matrix M, MY " Np0, Σq where Σ " MM t . Hence, we can generate the discretized Brownian motion using the inverse normal cumulative distribution function φ ´1: B " Mφ ´1pXq with φ ´1pXq " pφ ´1pX 1 q, . . . , φ ´1pX d qq t and X " U pr0, 1s d q. Thus, we can estimate the Sobol' indices over f payoff pMφ ´1pXqq. In this particular application, the choice of M is pivotal to reduce the effective dimensionality of the model. In this regard, Sobol' indices can help us assess whether we reduce the effective dimensionality. We will compare the case with either the Cholesky decomposition or the PCA construction.

We propose to analyze the results with averaged Sobol' estimates and averaged total number of evaluations over 100 repetitions. The tolerance ε and parameters ˚, r and Cpmq are set as previously.

Tables 7,8, and 9 show the results when M corresponds to the Cholesky decomposition of Σ. It turns out that Variant A.b does not perform better than Variant A.a. Averaged Sobol' estimates are similar but the average total number of model evaluations is slightly worse. The explanation is that total effect Sobol' indices require more evaluations than their corresponding first-order indices. As such, "Correlation 2" falls into case 1) explained in Section 4.2. Variant B shows a drastic improvement on the total number of evaluations with the drawback of a false estimation for input x 4 . Its discrepancy is higher than ε when comparing to the other two variants.

Tables 10,11, and 12 show the results when M corresponds to the PCA decomposition of Σ. Variant A.a performs as well as Variant A.b in terms of average total number of evaluations. However, Variant A.a is most probably giving worse estimates for the small first- order indices. Table 10 shows that averaged p S 4 is greater than p S 4 , p S 5 greater than p S 5 , and p S 6 greater than p S 6 , which is inconsistent with the property S u ď S u . Once again, variant B shows an improvement on the total number of evaluations with the drawback of slightly overestimating the small firstorder indices. Overall, each variant captures well how the PCA Brownian motion construction reduces the high effective dimensionality of the payoff function, in contrast to the initial Cholesky decomposition.

Conclusion

When estimating Sobol' indices, the question of how many evaluations must be performed to reach a desired precision is often raised by practitioners. This question is difficult to address mostly because the number of evaluations needed depends on the complexity of the model studied. As such, it is hard to bring out a general rule of thumb.

The sequential estimation procedure proposed in this article offers a practical solution with the construction of an estimator and error bound for Sobol' indices. The number of points is progressively augmented until the error bound becomes lower than a user specified tolerance. The procedure presented combines Sobol' sequences with either Saltelli's strategy to estimate both first-order and total effect indices, or the replication procedure to estimate only first-order indices. Furthermore, we investigated the use of a recent es-timator well-suited to the estimation of small first-order indices using quasi-Monte Carlo methods. The efficiency of this estimator, called "Correlation 2", was assessed and highlighted on two test functions. Overall, the variant combining Saltelli's strategy and "Correlation 2" gave the best results, with low failure rates across all indices.

As a future project, the same estimators and algorithms can be designed for rank-1 lattices using the results in [START_REF] Jiménez Rugama | Adaptive Multidimensional Integration Based on Rank-1 Lattices: Monte Carlo and Quasi-Monte Carlo Methods[END_REF], and noticing that rank-1 lattices of the same size are also replicated designs of order 1.

  If we estimate I by p I with vector of error bounds ε p I " pε p I 1 , ε p I 2 , ε p I 3 , ε p I 4 q according to (14), we know that I P B ε p I p p Iq " r p I ´εp I , p I ὲp I s. Then, as an alternative to the common Sobol' indices estimators (8), we can define the following two estimators with their respective error bounds: p S u " Iq S u pIq, 1 ¯`max ´min IPB ε p I p p Iq S u pIq, 0 ¯' ε Iq S u pIq, 1 ¯´max ´min IPB ε p I p p Iq S u pIq, 0 ¯'

Fig. 1 2 " 2 " 2 `Smax 2 q{2. 4

 12224 Fig. 1 The delimited region on the figure represents the values S 2 in I P B p0,0,0.05,0.1q p p Iq for the Bratley et al. function. There, S min 2 "

Algorithm 1

 1 are lower than the tolerance ε. For Variant A the stopping criterion is ε p S u ď ε and ε p S u ď ε for all u P D. For Variant B it is ε p S u ď ε for all u P D. If the stopping criterion is satisfied, the algorithm stops and Sobol' estimates are returned. Otherwise, m is incremented by one to perform a new estimation. Sequential estimation of Sobol' indices

Fig. 2

 2 Fig. 2 Boxplots of estimation errors δ S u for the four inputs x 3 , x 4 , x 5 , x 6 obtained with both Variant A.a and Variant A.b. The dashed horizontal line marks the tolerance ε " 5.10 ´3.

Table 1

 1 Averaged estimation errors δ S u , δ Su and total number of evaluations for Variant A.a.

	input	S u	δ S u	r δ S u	S u	δ Su	r δ Su
	x 1	0.5868 0.0004	0.02	0.6901	0.0005	0.02
	x 2	0.2608 0.0006	0.02	0.3562	0.0004	0
	x 3	0.0367 0.0011	0.03	0.0563	0.0008	0
	x 4	0.0058 0.0012	0.02	0.0092	0.0003	0
	x 5	5.10 ´5 0.0003	0	9.10 ´5 ă 10 ´4	0
	x 6	5.10 ´5 0.0003	0	9.10 ´5 ă 10 ´4	0
	Total number of evaluations: 63 088

Table 2

 2 Averaged estimation errors δ S u , δ Su and total number of evaluations for Variant A.b.

	input	S u	δ S u	r δ S u	S u	δ Su	r δ Su
	x 1	0.5868	0.0005	0.03	0.6901	0.0006	0.02
	x 2	0.2608	0.0007	0.03	0.3562	0.0003	0
	x 3	0.0367	0.0015	0.05	0.0563	0.0009	0
	x 4	0.0058	0.0003	0	0.0092	0.0003	0
	x 5	5.10 ´5 ă 10 ´4	0	9.10 ´5 ă 10 ´4	0
	x 6	5.10 ´5 ă 10 ´4	0	9.10 ´5 ă 10 ´4	0
	Total number of evaluations: 62 703

Table 3

 3 Averaged estimation errors δ S u and total number of evaluations for Variant B.

	input	S u	δ S u	r δ S u
	x 1	0.5868 0.0011	0.07
	x 2	0.2608 0.0013	0.08
	x 3	0.0367 0.0018	0.09
	x 4	0.0058 0.0013	0.1
	x 5	5.10 ´5 0.0031	0.04
	x 6	5.10 ´5 0.0036	0.05
	Total number of evaluations: 32 768
	or Variant A.b. However, the total number of evaluations is
	twice as small. As such, this approach remains interesting
	when one wants to estimate only first-order indices.
	5.1.2 Bratley et al. function
	In this second example, we consider the Bratley et al. func-
	tion defined by,
	f px 1 , . . . , x d q "

d ÿ i"1

Table 4

 4 Averaged estimation errors δ S u , δ Su and total number of evaluations for Variant A.a.

	input	S u	δ S u	r δ S u	S u	δ Su	r δ Su
	x 1	0.6529 0.0004	0	0.7396 0.0003	0.
	x 2	0.1791 0.0006	0	0.2659 0.0004	0
	x 3	0.0370 0.0058	0.71	0.0764 0.0012 0.02
	x 4	0.0133 0.0022	0.12	0.0343 0.0020 0.07
	x 5	0.0015 0.0024	0.11	0.0062 0.0014	0
	x 6	0.0015 0.0021	0.06	0.0062 0.0010 0.03
	Total number of evaluations: 70 129

Table 5

 5 Averaged estimation errors δ S u , δ Su and total number of evaluations for Variant A.b. Boxplots of estimation errors δ S u for the four inputs x 3 , x 4 , x 5 , x 6 obtained with both Variant A.a and Variant A.b. The dashed horizontal line marks the tolerance ε " 5.10 ´3.

	input	S u	δ S u	r δ S u	S u	δ Su	r δ Su
	x 1	0.6529 0.0004	0	0.7396 0.0003	0
	x 2	0.1791 0.0006	0	0.2659 0.0005	0
	x 3	0.0370 0.0019	0.04	0.0764 0.0013 0.01
	x 4	0.0133 0.0016	0.01	0.0343 0.0017 0.02
	x 5	0.0015 0.0003	0	0.0062 0.0015 0.04
	x 6	0.0015 0.0004	0	0.0062 0.0011	0
	Total number of evaluations: 68 045
	δ S u Fig. 3 1e-05 5e-05 5e-04 5e-03	Variant A.a Variant A.b
		x 3	x 4	x 5	x 6
				inputs

Table 6

 6 Averaged estimation errors δ S u and total number of evaluations for Variant B.

	input	S u	δ S u	r δ S u
	x 1	0.6529 0.0003	0
	x 2	0.1791 0.0006	0
	x 3	0.0370 0.0010	0
	x 4	0.0133 0.0038	0.25
	x 5	0.0015 0.0024	0
	x 6	0.0015 0.0026	0
	Total number of evaluations: 65 536

Table 7

 7 Averaged values of p S u , p S u , and total number of evaluations of Variant A.a with Cholesky decomposition.

	input	p S u	p S u
	x 1	0.2120	0.5650
	x 2	0.1370	0.4420
	x 3	0.0820	0.3150
	x 4	0.0408	0.1950
	x 5	0.0142	0.0944
	x 6	0.0048	0.0257
	Total number of evaluations: 648 821

Table 8

 8 Averaged values of p S u , p S u and total number of evaluations of Variant A.b with Cholesky decomposition.

	input	p S u	p S u
	x 1	0.2120	0.5650
	x 2	0.1370	0.4420
	x 3	0.0822	0.3150
	x 4	0.0434	0.1950
	x 5	0.0169	0.0944
	x 6	0.0049	0.0259
	Total number of evaluations: 699 487

Table 9

 9 Averaged values of p S u and total number of evaluations of Variant B with Cholesky decomposition.

	input	p S u
	x 1	0.2120
	x 2	0.1350
	x 3	0.0774
	x 4	0.0338
	x 5	0.0149
	x 6	0.0028
	Total number of evaluations: 262 144

Table 10

 10 Averaged values of p S u , p S u and total number of evaluations of Variant A.a with PCA decomposition.

	input	p S u	p S u
	x 1	0.9800	0.9960
	x 2	0.0037	0.0172
	x 3	0.0016	0.0023
	x 4	0.0014	0.0005
	x 5	0.0005	0.0001
	x 6	0.0002	0.0000
	Total number of evaluations: 396 221

Table 11

 11 Averaged values of p S u , p S u and total number of evaluations of Variant A.b with PCA decomposition.

	input	p S u	p S u
	x 1	0.9800	0.9960
	x 2	0.0034	0.0172
	x 3	0.0004	0.0023
	x 4	0.0001	0.0005
	x 5	0.0000	0.0001
	x 6	0.0000	0.0000
	Total number of evaluations: 398 648

Table 12

 12 Averaged values of p S u and total number of evaluations of Variant B with PCA decomposition.

	input	p S u
	x 1	0.9800
	x 2	0.0037
	x 3	0.0028
	x 4	0.0029
	x 5	0.0017
	x 6	0.0029
	Total number of evaluations: 262 144

Acknowledgements The authors thank Fred J. Hickernell and Clémentine Prieur for initiating this collaborative work and Elise Arnaud for her proofreading. The authors are grateful to Stephen Joe, Frances Y. Kuo and Art B. Owen for their helpful answers and suggestions.

Carlo

This work is supported by the CITiES project funded by the Agence Nationale de la Recherche (grant ANR-12-MONU-0020) and by the United States National Science Foundation (grant DMS-1522687).