Elena Di Bernardino

Thomas Laloë
email: laloe@unice.fr

Rémi Servien

Elena Di Bernardino

Estimating covariate functions associated to multivariate risks: a level set approach

Keywords: Multidimensional distribution function, plug-in estimation, regression function Mathematics Subject Classification (2000) 62G05, 62G20, 60E05, 91B30

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Traditionally, risk measures are thought of as mappings from a set of realvalued random variables to the real numbers. However, it is often insufficient to consider a single real measure to quantify risks, especially when the riskproblem is affected by other external risk factors. Note that the evaluation of an individual risk may strongly be affected by the degree of dependence amongst all risks. Modeling the dependency structure of multivariate data is helpful to obtain meaningful and accurate inference and prediction results in risk analysis. Several important challenges exist in this sense: the ability to incorporate important known covariates in the model; the dependence between variables is neither linear nor constant; they exhibit tail dependence, and so on. For instance, several hydrological phenomena are described by two or more correlated characteristics. These dependent characteristics should be considered jointly to be more representative of the multivariate nature of the phenomenon. Consequently, probabilities of occurrence of risks cannot be estimated on the basis of univariate analysis. The multivariate hydrological risks literature mainly treated one or more of the following three elements:

(1) showing the importance and explaining the usefulness of the multivariate framework, (2) fitting the appropriate multivariate distribution in order to model risks and (3) defining and studying multivariate return periods (see [START_REF] Chebana | Multivariate quantiles in hydrological frequency analysis[END_REF]). One of the most popular measures in hydrology and climate is undoubtedly the return period. It is closely related to the notion of quantile which has therefore been extensively studied in dimension one. For a random variable X that represents the magnitude of an event that occurs at a given time and at a given location, the quantile of order 1 -1 T expresses the magnitude of the event which is exceeded with a probability equal to 1 T . T is then called the return period. In univariate risk theory the quantile is known as the Value-at-Risk (VaR) and is defined by Q X (c) = inf{x ∈ R : F X (x) ≥ c}, for c ∈ (0, 1), with F X the univariate distribution of random variable X. A second important univariate risk measure, based on the quantile notion, is the Conditional-Tail-Expectation (CTE) defined by CTE c (X) = E[ X | X > Q X (c) ], for c ∈ (0, 1).

From the year 2000 onward, much research has been devoted to risk measures, and many extensions to multidimensional settings have been suggested (see, e.g., [START_REF] Jouini | Vector-valued coherent risk measures[END_REF]; [START_REF] Bentahar | Tail conditional expectation for vector-valued risks[END_REF]; [START_REF] Embrechts | Bounds for functions of multivariate risks[END_REF]; [START_REF] Nappo | Kendall distributions and level sets in bivariate exchangeable survival models[END_REF]; [START_REF] Ekeland | Comonotonic measures of multivariate risks[END_REF]). As a starting point, in the following, we consider the multivariate version of the CTE measure, proposed by Di Bernardino et al (2013) and [START_REF] Cousin | On multivariate extensions of Value-at-Risk[END_REF]. It is constructed as the conditional expectation of a multivariate random vector given that the latter is located in the c-upper level set of the associated multivariate distribution function. In this sense this measure is essentially based on a "multivariate distributional approach". More precisely they define, for i = 1, . . . , d and for c ∈ (0, 1),

CTE i c (X) = E[ X i | X ∈ L(c) ],
(1.1)

where X = (X 1 , . . . , X d ) is a non-negative multivariate risk portfolio with distribution function F . In particular, [START_REF] Cousin | On multivariate extensions of Value-at-Risk[END_REF] proved that properties of the multivariate Conditional-Tail-Expectation in (1.1) turn to be consistent with existing properties on univariate risk measures (positive homogeneity, translation invariance, increasing in risk-level c, . . . ). In the financial econometrics literature, we are often interested in analyzing the behavior of a univariate return measure (average return, skewness, . . . ) with respect to a set of d risk factors X (volatility or variance, kurtosis, . . .). In other words, we consider a dependent multivariate vector of risk-factors X and a univariate covariate Y (i.e. a dependent variable on X). Furthermore, in climatology, one may be interested in how climate change over years might affect high temperatures. Multivariate examples include the study of rainfall as a covariate function represented by the geographical location. In this sense [START_REF] Daouia | Kernel estimators of extreme level curves[END_REF] deal with the problem of estimating quantiles when covariate information is available.

So, the goal of this paper is the study of the behavior of a covariate Y on the level sets of a d-dimensional vector of risk-factors X. More precisely, adapting the multivariate risk measure in (1.1), we deal with the multivariate Covariate-Conditional-Tail-Expectation (CCTE) defined by:

CCTE c (X, Y ) := E[Y | X ∈ L(c)], (1.2)
where c ∈ (0, 1). In order to estimate this risk measure, we first need to estimate the level sets L(c) associated to the d-dimensional distribution function F of X. For a non-negative d-dimensional risk portfolio with distribution function F , the c-upper level set of F (i.e., L(c) = {x ∈ R d + : F (x) ≥ c}) and its associated c-level curve (i.e., ∂L(c) = {x ∈ R d + : F (x) = c}) have recently been proposed as risk measures in multivariate hydrological models. Among their many advantages, it appears that they are simple, intuitive, interpretable and probability-based (see [START_REF] Chebana | Multivariate quantiles in hydrological frequency analysis[END_REF]. A risk-problem of flood in the bivariate setting using an estimator of level curves ∂L(c) of the bivariate distribution function is proposed by de [START_REF] De Haan | Large quantile estimation in a multivariate setting[END_REF]. Furthermore, as noticed by [START_REF] Embrechts | Bounds for functions of multivariate risks[END_REF], ∂L(c) can be viewed as a natural multivariate version of the univariate quantile. The interested reader is also referred to [START_REF] Tibiletti | On a new notion of multidimensional quantile[END_REF], [START_REF] Belzunce | Quantile curves and dependence structure for bivariate distributions[END_REF], [START_REF] Nappo | Kendall distributions and level sets in bivariate exchangeable survival models[END_REF].

The problem of estimating level sets of an unknown function (for instance of a density function and more recently a regression function) has received attention. However, most of the existing literature has focused on the density or regression function [START_REF] Baíllo | Convergence rates in nonparametric estimation of level sets[END_REF]; [START_REF] Cavalier | Nonparametric estimation of regression level sets[END_REF]; [START_REF] Cuevas | Plug-in estimation of general level sets[END_REF]; [START_REF] Laloë | Nonparametric estimation of regression level sets using kernel plug-in estimator[END_REF]; [START_REF] Rigollet | Optimal rates for plug-in estimators of density level sets[END_REF]). [START_REF] Mason | Asymptotic normality of plug-in level set estimates[END_REF] obtained the asymptotic normality of plug-in level set estimates in the density case. As we consider the level sets of a multivariate distribution function, i.e. L(c) for some fixed c ∈ (0, 1), the commonly assumed property of compactness for these sets (required both in the density and in the regression cases) is no more reasonable. Then, differently from the literature cited above, in the present work, a special attention is given to this non-compact setting.

Considering a consistent estimator F n of the distribution function F , we propose a plug-in approach to estimate the c-upper level set L(c) by

L n (c) = {x ∈ R d + : F n (x) ≥ c}, for c ∈ (0, 1).
The regularity properties of F and F n as well as the consistency properties of F n will be specified in the statements of our theorems. Our consistency result for L(c) is stated with respect to a criterion of "physical proximity" between sets: the volume of the symmetric difference. Obviously, the convergence rate suffers from the well-known curse of dimensionality (see Theorem 1). Using Theorem 1, we state L p -consistency with a convergence rate for the estimation of the regression function

r(x) = E[Y | X = x],
on these level sets L(c), i.e., for x ∈ L(c), and c ∈ (0, 1) (see Theorems 2 and 3). The motivation behind the point-wise estimation of r(x) for x ∈ L(c) is an interesting problem, for different practical problems. Indeed, this represents the expected value of a covariate Y given that a dependent multivariate vector of risk-factors X takes value in a specific risk area L(c) (for instance L(c) can represent a risk-scenario or critical-layer and so on). For the importance of L(c) in the risk-management in the environmental or hydrological fields see for [START_REF] Chebana | Multivariate quantiles in hydrological frequency analysis[END_REF], [START_REF] Salvadori | On the return period and design in a multivariate framework[END_REF]. Finally, we provide a consistency result for the estimation of the CCTE risk measure, if Y is completely available (see Theorem 4) or not (see Theorem 5). Furthermore, we investigate the impact of a change in the scale of data on our results. In particular this property is related to the suitable positive homogeneity property of risk measures (e.g. see [START_REF] Artzner | Coherent measures of risk[END_REF]).

The paper is organized as follows. We introduce some notation, tools and technical assumptions in Section 2. Consistency and asymptotic properties of our estimator of L(c) are given in Section 3. Section 4 is devoted to the L pconsistency of the estimation of the regression function r on the level set L(c) and Section 5 to the consistency of the CCTE's estimation. The effects of scaling data are analyzed in Section 6. Illustrations with simulated data are presented in Section 7, and a comparison with parametric and semi-parametric approaches is detailed. A real example is studied in Section 8. Section 9 summarizes and briefly mentions directions for future research. Finally, proofs are postponed to Section 10.

Notation and preliminaries

In this section we introduce some notations and tools which will be useful later.

Let

N * = N \ {0}, R * + = R + \ {0} and R d + * = R d + \ {0}.
Let F be the set of continuous distribution functions R d + → [0, 1] and X := (X 1 , X 2 , . . . , X d ) a random vector with distribution function F ∈ F. Given an i.i.d. sample {X i } n i=1 in R d + with distribution function F ∈ F, we denote by F n an estimator of F based on this finite sample. Let Y be a random variable with values in J ⊂ R + , where J is supposed to be bounded. We denote by {(X i , Y i )} i=1,...,n the associated i.i.d. sample. Define, for c ∈ (0, 1), the c-upper level set of F ∈ F and its plug-in estimator

L(c) = {x ∈ R d + : F (x) ≥ c}, L n (c) = {x ∈ R d + : F n (x) ≥ c}, and
∂L(c) = {x ∈ R d + : F (x) = c}.
In addition, given T > 0, we set

L(c) T = {x ∈ [0, T ] d : F (x) ≥ c}, L n (c) T = {x ∈ [0, T ] d : F n (x) ≥ c}, ∂L(c) T = {x ∈ [0, T ] d : F (x) = c}.
Note that, in the presence of a plateau at level c, ∂L(c) can be a portion of quadrant R d + instead of a set of Lebesgue measure null in R d + . In the statement of our results we will require suitable conditions in order to avoid this situation. where (∇F ) x is the gradient vector of F evaluated at x and (∇F ) x its Euclidean norm, (HF ) x the Hessian matrix evaluated in x and (HF ) x M its matrix norm induced by the Euclidean norm.

We define:

g p := R d | g(x) | p f (x) dx 1/p
, where f denotes the density function associated to the probability measure µ,

g p,λ := R d | g(x) | p dx 1/p and g ∞ := sup | g(x) | : x ∈ R d .
For the sake of simplicity, we use from now on the notation . instead of . 2 .

Let r(x) be the regression function such that

r(x) := E[Y | X = x]
and let r n be a consistent estimate of r.

Finally, given two functions f and g, f (n) = o(g(n)) as n → ∞ means that for every positive constant there exists a constant N such that |f (n)| ≤ |g(n)|, for all n ≥ N. The o(•)-notation is called in the literature the Landau symbol.

3 Estimating level sets of a multidimensional distribution function using a plug-in method

We consider the consistency in terms of the volume (in the Lebesgue measure sense) of the symmetric difference between L(c) Tn and L n (c) Tn . This means that we define the distance between two subsets A 1 and A 2 of R d + by

d λ (A 1 , A 2 ) = λ(A 1 A 2 ),
where λ stands for the Lebesgue measure on R d and for the symmetric difference.

Let us introduce the following assumption:

A1 There exist positive increasing sequences

(v n ) n∈N * and (T n ) n∈N * such that v n → ∞, T n → ∞, and
v n [0,Tn] d | F (x) -F n (x) | p λ(dx) P → n→∞ 0, for some 1 ≤ p < ∞.
We now establish our consistency result with convergence rate.

Theorem 1 Let c be in (0, 1). Let F ∈ F be a twice differentiable distribution function on R d * + . Assume that there exist κ > 0, ζ > 0 such that m > 0 and M H < ∞. Let T 1 > 0 such that for all t : | t -c | ≤ κ, ∂L(t) T1 = ∅. Assume that for each n, F n is measurable. Let (v n ) n∈N * and (T n ) n∈N * positive increasing sequences such that Assumption A1 is satisfied. Then, it holds that

p n d λ (L(c) Tn , L n (c) Tn ) P → n→∞ 0, with p n an increasing positive sequence such that p n = o v 1 p+1 n /T (d-1) p p+1 n .
Remark 1 Note that in the univariate case (d = 1) the convergence rate of Theorem 1 does not depend on the truncation sequence T n . The interested reader is referred to the proof of this result in Section 10 for further details.

Theorem 1 provides a convergence rate, which is closely related to the choice of the sequence T n , for d > 1. Note that, as in Theorem 3 in [START_REF] Cuevas | Plug-in estimation of general level sets[END_REF], Theorem 1 above does not require any continuity assumption on F n . Furthermore, we remark that a sequence T n whose divergence rate is large, implies a convergence rate p n quite slow. Moreover, this phenomenon is emphasized by the dimension d of the data, and we face here the well-known curse of dimensionality. In the following we will illustrate this aspect by giving convergence rate in the case of the empirical distribution function (see Example 1). Firstly, from Theorem 1 we can derive the following result.

Corollary 1 Let c ∈ (0, 1). Let F ∈ F be a twice differentiable distribution function on R d * + . Assume that there exist κ > 0, ζ > 0 such that m > 0 and M H < ∞. Let T 1 > 0 such that for all t : | t -c | ≤ κ, ∂L(t) T1 = ∅.
Assume that for each n, F n is measurable. Assume that there exists a positive increasing sequence

(w n ) n∈N * such that w n F -F n ∞ P → n→∞ 0. Then, it holds that p n d λ (L(c) Tn , L n (c) Tn ) P → n→∞ 0,
with p n an increasing positive sequence such that

p n = o w n p p+1 /T n d+(d-1) p p+1
.

This result comes trivially from Theorem 1 and the fact that

w n F -F n ∞ P → n→∞ 0 implies ∀ p ≥ 1, v n [0,Tn] d | F -F n | p λ(dx) P → n→∞ 0, with v n = w p n T d n .
Let us now present a more practical example in the case of a d-variate empirical distribution function.

Example 1 Let F n the d-variate empirical distribution function. Then, it holds that

w n F -F n ∞ P → n→∞ 0, with w n = o( √ n).
From Corollary 1, with p = 2, we obtain for instance:

p n = o n 1/3
Tn 7/3 , for d = 3;

p n = o n 1/3
Tn 10/3 , for d = 4.

L p -consistency of r n on the level sets

In this section we study the L p -consistency of an estimator r n of the regression function r(x), for x ∈ L(c). More precisely, in the first part we provide a consistency result in terms of the L p -distance of the absolute error between r n 1 Ln(c) Tn and r 1 L(c) (see Theorem 2). In the second part we analyze the problem of a convergence rate (see Theorem 3). Remind that the level set L(c) is not compact. The point-wise estimation of r(x) for x ∈ L(c) is an interesting aim, for different practical problem, as previously discussed in the introduction of this paper. Then, since the final goal is to study the behavior of the regression function r(x), knowing that x is in the "risk area" L(c), it is natural to check if r n (x), for x in the estimated level set L n (c) is globally close to r(x), for x in the theoretical level set L(c).

Let us introduce the following assumptions:

A2 F n -F p,λ P → 0 and r n -r p P → 0, for some 1 ≤ p < ∞. A3 The density function f of X is such that f 1+ ,λ < ∞, with > 0.
Theorem 2 Let c be in (0, 1). Let F ∈ F be a twice differentiable distribution function on R d * + with an associated density f such that Assumption A3 is satisfied. Assume that there exist κ > 0, ζ > 0 such that m > 0 and M H < ∞. Assume that for each n, F n is measurable. Assume that condition A2 is satisfied. Let (T n ) n∈N * be an increasing sequence of positive values with

T 1 > 0 such that for all t : | t -c | ≤ κ, ∂L(t) T1 = ∅. Then it holds that r n 1 Ln(c) Tn -r 1 L(c) p P → 0, for n → ∞. Remark 2 If r n -r ∞ → 0, a.s. we have r n 1 Ln(c) -r 1 L(c) ∞ = sup R d |r n 1 Ln(c) -r 1 L(c) | 0, a.s. Indeed, we have sup R d |r n 1 Ln(c) -r 1 L(c) | = max x∈Ln(c) L(c) (r(x), r n (x)) which is greater than zero if L n (c) L(c) = ∅.

Let us introduce the following assumption:

A4 There exist positive increasing sequences

(v 1,n ) n∈N * , (v 2,n ) n∈N * and (T n ) n∈N * such that v 1,n → ∞, v 2,n → ∞, T n → ∞, and
v 1,n r n -r p P → 0, v 2,n [0,Tn] d | F (x) -F n (x) | p λ(dx) P → 0,
for some 1 ≤ p < ∞.

Note that the control of the convergence rate of the distribution function estimate which was previously handled by Assumption A1 is now included in Assumption A4.

Theorem 3 Let c be in (0, 1). Let F ∈ F be a twice differentiable distribution function on R d * + with an associated density f such that Assumption A3 is satisfied. Assume that there exist κ > 0

, ζ > 0 such that m > 0 and M H < ∞. Let T 1 > 0 such that for all t : | t -c | ≤ κ, ∂L(t) T1 = ∅. Assume that for each n, F n is measurable. Let (v 1,n ) n∈N * , (v 2,n ) n∈N * and (T n ) n∈N * positive increasing sequences such that Assumption A4 is satisfied. Then it holds that w n r n 1 Ln(c) Tn -r 1 L(c) Tn p P → 0, for n → ∞.
where

w n = min {v 1,n , a n } with a n = o v ( +1)p(p+1) 2,n /T (d-1) ( +1)(p+1) n .
Remark 3 Note that we do not provide a convergence rate for r n 1 Ln(c) Tnr 1 L(c) p because it implies the knowledge of the rate of convergence of

r 1 L(c) Tn -r 1 L(c) p = L(c)\[0,Tn] d |r(x)| p f (x)dx 1/p
. Indeed the rate of

convergence of P[X 1 > T n , . . . , X d > T n ] → 0, for n → ∞, is unknown.
In order to overcome this problem we should have to assume that the vector X belongs to a particular distribution class in order to know the rate of decay of its multivariate tails.

5 Covariate-Conditional-Tail-Expectation consistency

A risk measure has recently received growing attention in risk theory literature: the CTE measure. According to [START_REF] Artzner | Coherent measures of risk[END_REF], [START_REF] Dedu | Restricted optimal retention in stop-loss reinsurance under VaR and CTE risk measures[END_REF], [START_REF] Denuit | Plug-in estimation of level sets in a non-compact setting with applications in multivariable risk theory[END_REF], for a continuous univariate loss distribution function F X the CTE at level c ∈ (0, 1) is defined by

CTE c (X) = E[ X | X ≥ F -1 X (c)].
Several multivariate generalizations of the classical univariate CTE have been proposed (see for instance [START_REF] Cai | Conditional tail expectations for multivariate phase-type distributions[END_REF]). Using the same approach as Di Bernardino et al ( 2013), we define the multivariate CTE in such a way as to preserve both the complete information about dependence structure between X and Y , and the marginal behavior of each component.

Definition 1 Consider a random vector X with distribution function F and a random variable Y . For c ∈ (0, 1), we define the theoretical multivariate c-Covariate-Conditional-Tail-Expectation as

CCTE c (X, Y ) = E [Y | X ∈ L(c)] .
Using the truncated version of the c-upper level set defined in Section 2 we also define

CCTE Tn c (X, Y ) = E Y | X ∈ L(c) Tn .
In the following we define consistent estimates for the CCTE in Definition 1.

Covariable Y is measured

In this subsection, we assume to have an i.i.d. sample {(X i , Y i )} i=1,...,n (see Section 2) and we introduce the following estimate for the CCTE.

Definition 2 Consider a random vector X with distribution function F and a random variable Y . For c ∈ (0, 1), using the truncated theoretical multivariate c-Covariate-Conditional-Tail-Expectation introduced in Definition 1, we define the associated estimate as Tn , where E n denotes the empirical version of the expected value.

CCTE Tn c,n (X, Y ) = E n Y |X ∈ L n (c)
Using Definition 2, we now establish the consistency of this estimate with a convergence rate.

Theorem 4 Let c be in (0, 1). Let F ∈ F be a twice differentiable distribution function on R d * + with an associated density f such that Assumption A3 is satisfied. Assume that there exist κ > 0, ζ > 0 such that m > 0 and M H < ∞.

Let T 1 > 0 such that for all t : | t -c | ≤ κ, ∂L(t) T1 = ∅. Assume that for each n, F n is measurable. Let (v 1,n ) n∈N * , (v 2,n ) n∈N * and (T n ) n∈N * positive increasing sequences such that Assumption A4 is satisfied. We have β n CCTE Tn c,n (X, Y ) -CCTE Tn c (X, Y ) P ----→ n→∞ 0, with β n = min p 2(1+ ) n , d n , where p n = o v 1 p+1 2,n /T (d-1)p p+1 n and d n = o ( √ n).
Remark 4 Obviously, it could also be very interesting to consider the con-

vergence CCTE Tn c,n (X, Y ) -CCTE c (X, Y ) .
We remark that in this case the speed of convergence will also depend on the convergence rate to zero of

CCTE Tn c (X, Y ) -CCTE c (X, Y ) , then, in particular of P[X ∈ L(c) \ L(c) Tn ],
for n → ∞. This could not be done without adding strong assumptions and careful developments and, by consequence, remains out of the scope of this paper.

Example 2 Let F n the d-variate empirical distribution function. Then, as in Example 1 and using Theorem 4 with p = 2, we obtain for instance:

p n = o n 1/3 T 7/3 n , for d = 3.
This gives us β n = min o n 1/3 T 7/3 n 2(1+ ) , o n 1/2 . In the case of a bounded density function f we obtain

β n = min o n 1/3 T 7/3 n 1 2 , o n 1/2 = o n 1/6
T 7/6 n .

Covariable Y is partially unknown

In this subsection, we deal with a more difficult case. We suppose that the covariable Y cannot be measured for all the individuals. It could happen if a measure of Y is very expensive or invasive (in some medical treatment, for example). So we have two different i.i.d. samples:

S 1 N = {(X i , Y i )} N i=1 and S 2 n = {X j } n j=1
, with n potentially much bigger than N .

In this case we use S1 N to get an estimate r N of the regression function r. Then, we apply this estimate on the sample S 2 n in order to estimate the CCTE measure. To this aim we define: Definition 3 Consider a random vector X with distribution function F , a random variable Y and two i.i.d. samples S 1 N and S 2 n . Using the truncated version of the c-upper level set defined in Section 2, for c ∈ (0, 1), we define

CCTE Tn c,n,N (X, Y ) = E n r N (X)|X ∈ L n (c) Tn .
The following result proves the consistency of the estimate introduced in Definition 3.

Theorem 5 Let c be in (0, 1). Let F ∈ F be a twice differentiable distribution function on R d * + with an associated density f such that Assumption A3 is satisfied. Assume that there exist κ > 0,

ζ > 0 such that m > 0 and M H < ∞. Let T 1 > 0 such that for all t : | t -c | ≤ κ, ∂L(t) T1 = ∅. Assume that for each n, F n is measurable. Let (v 1,N ) N ∈N * , (v 2,n ) n∈N * and (T n ) n∈N *
positive increasing sequences such that Assumption A4 is satisfied. Let r be a continuous positive regression function such that E r(X) 2 < ∞. We have

β n,N CCTE Tn c,n,N (X, Y ) -CCTE Tn c (X, Y ) P → 0, for N, n → ∞, with β n,N = min p 2(1+ ) n , c N , d n , where p n = o v p p+1 2,n /T d+(d-1)p p+1 n , c N = o (E |r N (X) -r(X)|) and d n = o ( √ n).
Note that we have a supplementary term (c N ) comparing to Theorem 4. This term controls the rate of convergence of r N toward r.

Example 3 Using the same context as in Example 2, taking r N the kernel regression estimator of [START_REF] Kohler | Optimal global rates of convergence for nonparametric regression with unbounded data[END_REF], assuming that r is a (1 + k, C)smooth with 0 < k ≤ 1 function 1 and that there exist ψ > 2 (1 + k) such that E X ψ < ∞, we have c N = o N 2/5 . Then we obtain

β n,N = min o n 1/3 T 7/3 n 2(1+ ) , o N 2/5 , o n 1/2 .
In the case of a bounded density function f and as long as N >> n 5/12 T 35/12 n , we obtain

β N,n = min o n 1/3 T 7/3 n 1 2 , o N 2/5 , o n 1/2 = o n 1/6 T 7/6 n ,
which is the same rate as in Example 2.

6 About the effects of scaling data

In this section we study the impact of a change in the scale of data with respect to the estimate CCTE Tn c,n (X, Y ). In particular, this property is related to the suitable positive homogeneity property of risk measures (e.g. see [START_REF] Artzner | Coherent measures of risk[END_REF]). In the literature, the homogeneity property of a risk measure is often motivated by a change of currency argument: the amount of required capital in order to manage risks should be independent of the currency in which it is expressed. Suppose that we scale our data using a scale parameter a ∈ R * + . In our case, the scaled random vector will be (a X 1 , a X 2 , . . . , a X d ) =: a X. From now on we denote F a X (resp. F X ) the distribution function associated to a X (resp. to X). Using notation of Section 2, let

L a (c) = {x ∈ R d + : F a X (x) ≥ c}.
We can now consider the effects of scaling data with respect to the volume of the symmetric difference. Corollary 2 is a straightforward consequence of Theorem 1 using the fact that

L Tn a (c) L Tn n,a (c) = a L Tn (c) L Tn n (c) , fora ∈ R * + .
Corollary 2 Assume that F n,aX (x) = F n,X (x/a).

Under same notations and assumptions of Theorem 1 it holds that Using Theorem 6 we obtain that the rates of convergence of Theorem 3 and Theorem 4 are not affected by any scaling on the data.

p n, a d λ (L a (c) Tn , L n, a (
Remark 5 Note that if r n is the classical kernel estimator, assumption r n,aX (x) = r n,X

x a is not automatically satisfied. However, it can be satisfied if the scaling is also applied to the bandwidth.

Illustrations

In the following we consider some different simulated cases for which we illustrate the finite sample properties of our estimation of r and CCTE c (X, Y ). In particular, we will consider an independent copula (Section 7.1), and Ali-Mikhail-Haq copula (Section 7.2). To compare the estimated results with the theoretical ones we consider cases for which we can calculate (using Maple) the explicit value of the theoretical CCTE c (X, Y ). However, our estimator can be applied to much more general cases.

In this section we consider the kernel regression estimate proposed by [START_REF] Kohler | Optimal global rates of convergence for nonparametric regression with unbounded data[END_REF] in order to estimate r. Furthermore, the plug-in estimation of level sets, i.e., L n (c), is constructed using the empirical estimator F n of the distribution function. Remark that these two estimates both satisfy Assumptions A1, A2 and A4, and that the considered random vectors X in this section satisfy the conditions required for the kernel regression estimate (see Example 3 page 11).

Independent Copula

We consider here a bivariate independent copula with exponentially distributed marginals with parameter 1. Furthermore we choose the bounded regression function y = r(x 1 , x 2 ) = 1/(e x1 e x2 ), for all x 1 , x 2 ≥ 0 (Figure 1).

L p -consistency of r n

In this section we provide an illustration of Theorem 3. Remark that the assumptions of Theorem 3 are satisfied. In particular, we denote

E L2 := r n 1 Ln(c) Tn -r 1 L(c) Tn 2 . x1 x 2 r(x1, x2) Figure 1 Regression function: y = r(x 1 , x 2 ) = 1/(e x 1 e x 2 ) for all x 1 , x 2 ≥ 0.
According to Remark 3, we take T n = ln(n) (the interested reader is also referred to Di Bernardino et al ( 2013)). In Table 1 we illustrate the mean of E L2 (i.e., L 2 -consistency for the estimator r n ) on M = 50 simulated samples, for different level of risk c and different sample size n. Table 1 E L 2 for independent copula and exponentially distributed marginals with parameter 1. The bounded regression function is y = r(x 1 , x 2 ) = 1/(e x 1 e x 2 ) for all x 1 , x 2 ≥ 0.

As expected, the performance of the estimation increases with the size of the sample. However, the fact that the quality for the E L2 error seems better for large risk level c is quite surprising. This behavior of the estimation comes from the fact that the regression r has a flat plateau for x 1 , x 2 → ∞ (see Figure 1).

CCTE c (X, Y ) estimation

In the following, we compare CCTE Tn c,n (X, Y ) with the theoretical CCTE c (X, Y ). Following Remark 4, we obtain that CCTE c (X, Y ) -CCTE Tn c (X, Y ) decays to zero at least with a convergence rate o n 1/6 ln(n) 4/6 , with a choice of sequence T n = ln(n). This kind of compromise provides an illustration on how to choose T n , apart from satisfying the assumptions of Theorem 4.

In the following tables we denote Mean= CCTE Tn c (X, Y ), i.e., the mean of CCTE Tn c (X, Y ) on M = 50 simulated samples. We denote σ the empirical standard deviation

σ = 1 M -1 M j=1 CCTE Tn c (X, Y ) j -CCTE Tn c (X, Y ) 2 .
Finally, we denote RMAE the Relative Mean Absolute Error, i.e.,

RMAE= 1 M M j=1 CCTE Tn c (X,Y )j -CCTEc(X,Y ) CCTEc(X,Y ) .
The results are gathered in Table 2. Furthermore, in Table 3, we provide an illustration of the convergence rate of Theorem 4. In this case, we remark that As expected, the greater n is, the better the estimations are. Furthermore, results in Table 3 set out how β n = n 1/6 ln(n) 4/6 is at least the convergence rate of CCTE Tn c,n (X, Y ) -CCTE c (X, Y ) , in this particular case.

β n = n 1/6 ln(n) 4/

Ali-Mikhail-Haq Copula

We consider here a bivariate Ali-Mikhail-Haq copula with parameter 0.5 and exponentially distributed marginals with parameter 1 (see [START_REF] Nelsen | An introduction to copulas[END_REF]). We take the bounded regression function y = r(x 1 , x 2 ) = -(-2 + e -x1-x2 ) 3 , for all x 1 , x 2 ≥ 0. Conversely to the previous section we now deal with dependent variables X 1 , . . . , X n (i.e., only the coordinates of the X i are dependent) and an increasing regression function r (see Figure 2).

= r(x 1 , x 2 ) = -(-2 + e -x 1 -x 2 ) 3 for all x 1 , x 2 ≥ 0.

L p -consistency of r n

In this section we provide an illustration of Theorem 3. We can verify that the assumptions of Theorem 3 are satisfied. As explained in Remark 4, we choose T n = n 0.2 . In Table 4

= r(x 1 , x 2 ) = -(-2 + e -x 1 -x 2 ) 3 , for all x 1 , x 2 ≥ 0.
As before, we see that the quality of the estimation is better when the size of the sample growths. Again, the presence of a plateau for r(x 1 , x 2 ) when x 1 , x 2 → ∞ explains the good results for large risk levels c (see Figure 2).

CCTE c (X, Y ) estimation

In the following, we compare CCTE Tn c,n (X, Y ) with the theoretical CCTE c (X, Y ) in the case of Ali-Mikhail-Haq copula with parameter 0.5 (see Table 5). Furthermore, following Remark 4, we obtain that CCTE c (X, Y ) -CCTE Tn c (X, Y ) decays to zero at least with a convergence rate o n 1/10 , with a choice of sequence T n = n 1/5 . Again, our theoretical results are confirmed by these simulations (see Table 6).

A comparative study

In this section we provide a detailed comparative study to illustrate the advantages of our non-parametric estimation procedure. In particular, we compare the performances of the estimation of CCTE c using a parametric (resp.

semi-parametric) estimation procedure based on Maximum Likelihood Estimation (MLE), by fitting a parametric copula model with parametric (resp. non-parametric) marginals (see Tables 7, 8 and9).

We consider the following parametric estimator of the Covariate-Conditional-Tail-Expectation, i.e.,

CCTE P c = E n [Y |X ∈ L-→ θ ,c
], for c ∈ (0, 1), where

L-→ θ ,c = {x ∈ R d + : F-→ θ ,c (x) ≥ c}, and F-→ θ ,c (x) = C θ ( F θ X 1 (x 1 ), . . . , F θ X d (x d )),
with C θ a copula with parameter θ fitted by MLE, and F θ X i the parametric marginal distribution of X i , for i = 1, . . . , d with parameter (or eventually a vector of parameters) θ Xi fitted by MLE.

We also introduce the semi-parametric estimator of the Covariate-Conditional-Tail-Expectation, i.e.,

CCTE SP c = E n [Y |X ∈ L θ,n,c ], for c ∈ (0, 1), where L θ,n,c = {x ∈ R d + : F θ,n,c (x) ≥ c}, F θ,n,c (x) = C θ (F X1,n (x 1 ), . . . , F X d ,n (x d )),
with C θ a copula with parameter θ fitted by MLE, and F Xi,n the empirical marginal distribution of X i , for i = 1, . . . , d.

The idea is to compare the performance of these three different estimators (non-parametric, semi-parametric and parametric). We aim to consider two different situations: when the copula model used in the parametric and semiparametric estimator is the good one or not. For the second case we consider the Clayton and Gumbel copulas. The Clayton copula is closer to the true one (Ali-Mikhail-Haq) since they are both non-heavy tail copulas (upper-tail dependence coefficient λ U = 0, see [START_REF] Nelsen | An introduction to copulas[END_REF]). Conversely the Gumbel copula presents an upper-tail dependence structure. In this sense the Gumbel copula is substantially different from the others two. The results are gathered in Tables 7, 8 and9.

Not surprisingly, we see that our estimator performs better than the parametric and the semi-parametric one in the case when the chosen class of copulas does not contain the true model (i.e., here Clayton and Gumbel copulas). And as expected the worst results are for the Gumbel copula as it is the farthest from the true model. We also note that in this case the semi-parametric approach gives no better results than the parametric one. Remark also that our non-parametric estimation is comparable to the parametric and semiparametric ones in case the chosen class of copulas contains the true model. 8 Real data study: waves and water levels in coastal engineering design On coasts with high tidal ranges, or subject to high surges, both still water levels and waves can be important in assessing flood risk; their relative importance depends on location and on the type of sea defence. The simultaneous occurrence of large waves and a high still water level is, therefore, important in estimating their combined effect on sea defences. In design of a sea defence, a key step is the estimation of the probability of failure to protect against sea conditions. It is important in engineering design to identify the combinations of sea condition variables which cause each failure. The interested reader is referred for instance to [START_REF] Hawkes | The joint probability ofwaves andwater levels in coastal engineering design[END_REF].

For any particular mode of failure (structural failure, excessive overtopping, . . . ), the regression function r is dependent on the sea condition variables. Then, in particular, at any particular time t, the overtopping covariate Y will be related to the sea condition vector (X). In the literature the sea condition variables are often represented by the significant Wave height Hm0 (X 1 ), the Still Water level SWL(X 2 ), and the Wave period Tpb(X 3 ), then X = (X 1 , X 2 , X 3 ) (see Figure 3).

The regression function r(x) := E[Y | X = x] represents the relationship between the sea conditions and the overtopping at a given time t. This relationship could be complex and in some real analysis, can be represented by equations. The most advantage of the use of this regression function is to reduce a joint probability risk problem to a single covariate problem.

In this section we analyze the Wave height (Hm0), Still Water level (SWL), Wave period (Tpb) data, recorded during 828 storm events spread over 13 years in front of the Dutch coast near the town of Petten (Figure 4). These data has been recently studied in the literature (for details see for instance [START_REF] Draisma | Bivariate tail estimation: dependence in asymptotic independence[END_REF]). Following [START_REF] Tau | Preliminary design study. Project group Flood Defence[END_REF], at a given time t, the principal equation used for overtopping discharge (l/m/s) Y is given by:

Y = a e -b (h-SWL) Hm0 g (Hm0) 3 , (8.1) with a = 0.04 tan(β) L Hm0 , b = 7.05 √ Hm0 L tan(β) , L = g (Tpb) 2 2 π
, where -Hm0 (m) is the wave height at the toe of the structure at time t; -Tpb(s) is the wave period at time t. In particular the number of waves in a storm (N ) can easily be computed from information about the wave period and the duration of the storm, i.e., N = (duration(h)/Tpb(s))•3600; -SWL(m) the level of the sea if it is flat, without any waves at time t; -h (m) is the height of the costal design above SWL= 0 (see Figure 3); -g (m/s 2 ) the gravitational acceleration (i.e., 9.8 m/s 2 ); -β (rad) is the seaward slope steepness. In the following we consider the case tan(β) = 0.3 and tan(β) = 0.6 (see Figure 3). Note that Theorem 6 allows us to make any changes on the units before or after the estimation. We now estimate the mean overtopping rate (i.e., CCTE c (X, Y )) using Equation (8.1) for Y and conditionally to the fact that the sea variable conditions X belong to the joint risk area L(c). In particular, we consider the jointly large realisations of sea variables Hm0, SWLand Tpbfor different values of risk level c. Obviously, the dangerous effects of the sea conditions will be increasing with respect to c. The results are gathered in Figure 5.

As expected, we obtain a decreasing behavior of mean overtopping discharge according to the height of the costal design h for a fixed risk level c. Furthermore, for a fixed height of the costal design, the higher the risk c, the higher the mean overtopping discharge. In our study we consider both the case tan(β) = 0.3 (Figure 5, left) and tan(β) = 0.6 (Figure 5, right). As expected, we remark that for fixed risk level c and height of costal design h, a lower seaward slope steepness β generates a smaller mean overtopping discharge.

Conclusion

We propose in this paper a generalization to the estimation of the level sets of a d-variate distribution function. The non-compactness of the level sets requires special attention in the statement of the problem. The consistency results with a convergence rate are stated in terms of the volume of the symmetric difference. In a second part, we analyze the problem of the estimation of a regression function on the level sets of a d-variate distribution function and we obtain the consistency with a convergence rate in terms of the L p -distance. Then, we study a new multivariate risk measure: the Covariate-Conditional-Tail-Expectation, i.e. the Conditional-Tail-Expectation of the regression function.

A consistent estimator and a rate of convergence are provided. Moreover, we analyze the impact of scaling data on our results. Our theoretical results are illustrated on a complete simulation study. We discuss a real application in the evaluation of the mean overtopping discharge conditionally to the fact that the sea variable conditions belong to some joint risk area. It highlights the importance of the parameter T n (which solved the problem of the non-compactness of the level sets) as well as the curse of dimensionality. An interesting future work could be a deep investigation about these points, with a focus on the optimal choice for this parameter. Furthermore, the proposed methods are based on an i.i.d. samples framework. We remark that in real applications such as seasonal pattern in the temperature and water level rise series, data can have different types of serial correlations like nonlinear or non-stationary correlations [START_REF] Fan | Nonlinear time series[END_REF]. It would be interesting in a future work to include other more interesting and complex types of serial correlation structures and to analyze how this affects the performance of the proposed procedure. Finally, a study on the lower bounds of our estimation problem could be an interesting development of the present work.

Proofs

Proof of Theorem 1: Under assumptions of Theorem 1, we can always take T 1 > 0 such that for all t : | t -c | ≤ κ, ∂L(t) T1 = ∅. Then for each n, for all t : | t -c | ≤ κ, ∂L(t) Tn is a non-empty (and compact) set on R d + .

We consider a positive sequence

ε n such that ε n → n→∞ 0. For each n ≥ 1 the random sets L(c) Tn L n (c) Tn , Q εn = {x ∈ [0, T n ] d : | F -F n |≤ ε n } and Q εn = {x ∈ [0, T n ] d : | F -F n |> ε n } are measurable and λ(L(c) Tn L n (c) Tn ) = λ(L(c) Tn L n (c) Tn ∩ Q εn )+λ(L(c) Tn L n (c) Tn ∩ Q εn ). Since L(c) Tn L n (c) Tn ∩ Q εn ⊂ {x ∈ [0, T n ] d : c -ε n ≤ F < c + ε n } we obtain λ(L(c) Tn L n (c) Tn ) ≤ λ({x ∈ [0, T n ] d : c -ε n ≤ F < c + ε n }) + λ( Q εn ).
From assumptions of Theorem 1 and Proposition 2.1 in Di Bernardino et al (2013), it follows that there exists a γ > 0 such that, if 2 ε n ≤ γ then

d H (∂L(c + ε n ) Tn , ∂L(c -ε n ) Tn ) ≤ 2 ε n A
where A = 2 m ∇ and d H is the Hausdorff distance. From assumptions on first derivatives of F and Property 1 in [START_REF] Imlahi | Estimación de la curva mediana de una cópula C(x 1 , . . . , xn)[END_REF], we can write

λ({x ∈ [0, T n ] d : c -ε n ≤ F < c + ε n }) ≤ (2 ε n A) d T d-1 n .
Interestingly we remark that in the univariate case (d = 1) the Hausdorff distance between the two points ∂L(c -ε n ) Tn and ∂L(c + ε n ) Tn is also the Lebesgue measure (in dimension 1) for this interval. Then λ({x

∈ [0, T n ] : c -ε n ≤ F < c + ε n }) ≤ 2 ε n A.
This means that in this case, the result does not depend on the truncation sequence T n .

If we now choose

ε n = o 1 p n T d-1 n , (10.1)
we obtain that, for n large enough, 2 ε n ≤ γ and

p n λ({x ∈ [0, T n ] d : c -ε n ≤ F < c + ε n }) → n→∞ 0.
Let us now prove that p n λ( Q εn ) P → n→∞ 0. To this end, we write

p n λ( Q εn ) = p n 1 {[0,Tn] d : |F -Fn|>εn} λ(dx) ≤ p n ε p n [0,Tn] d | F -F n | p λ(dx).
Take ε n such that n we can choose ε n that satisfies (10.1) and (10.2). Hence the result.

ε n = p n v n 1 p . ( 10
In the following proofs, K denotes a constant whose value may change from line to line.

Proof of Theorem 2: We have

r 1 Ln(c) Tn -r 1 L(c) Tn p = R d r(x) 1 Ln(c) Tn -r(x) 1 L(c) Tn p f (x) dx 1/p ≤ L(c) Tn Ln(c) Tn |r(x)| p f (x) dx 1/p ≤ K µ L(c) Tn L n (c) Tn 1/p ≤ K f 1+ ,λ λ L(c) Tn L n (c) Tn 1+ 1/p which gives us r n 1 Ln(c) Tn -r 1 L(c) Tn p ≤ r n -r p + K f 1+ ,λ λ L(c) Tn L n (c) Tn 1+ 1/p . Since r 1 L(c) Tn -r 1 L(c) p = L(c)\[0,Tn] d |r(x)| p f (x)dx 1/p , we finally get r n 1 Ln(c) Tn -r 1 L(c) p = r n 1 Ln(c) Tn -r 1 L(c) Tn + r 1 L(c) Tn -r 1 L(c) p ≤ r n 1 Ln(c) Tn -r 1 L(c) Tn p + r 1 L(c) Tn -r 1 L(c) p ≤ r n -r p + K f 1+ ,λ λ L(c) Tn L n (c) Tn 1+ 1/p + L(c)\[0,Tn] d |r(x)| p f (x)dx 1/p .
From Theorem 1 and assumptions of Theorem 2, the last inequality concludes the proof.

Proof of Theorem 3: Note that the proofs of Theorem 3 and Theorem 1 are strongly related. We have Hence the result.

  We denote by B(x, ρ) the closed ball centered on x ∈ R d + and with positive radius ρ. Let B(S, ρ) = x∈S B(x, ρ), with S a closed set of R d + . For κ > 0 and ζ > 0, define E = B({x ∈ R d + : | F (x) -c |≤ κ}, ζ), and, for a twice differentiable function F ,

  means that the scale factor a ∈ R * + impacts the volume in R d with an exponent d. Let us now investigate the effects of scaling data on Theorem 3 and Theorem 4. Let denote r aX,bY (x) = E[b Y | a X = x] (r n,aX,bY its estimate version) and CCTE Tn c (aX, bY ) = E b Y | a X ∈ L a (c) Tn ( CCTE Tn c,n (aX, bY ) its estimate version). We obtain the following result. Theorem 6 Let a and b in R * + and assume that r n,aX,bY (x) = b r n,X,Y x a . 1. Under same hypotheses and notations of Theorem 3 we have b r n,X,Y 1 Ln(c) Tn -r X,Y 1 L(c) Tn p = r n,aX,bY 1 La,n(c) Tn -r aX,bY 1 La(c) Tn p . 2. Under same hypotheses and notations of Theorem 4 we have b CCTE Tn c,n (X, Y ) -CCTE Tn c (X, Y ) = CCTE Tn c,n (aX, bY ) -CCTE Tn c (aX, bY ) .

Figure 2

 2 Figure 2 Regression function y = r(x 1 , x 2 ) = -(-2 + e -x 1 -x 2 ) 3 for all x 1 , x 2 ≥ 0.

  of the estimations of the CCTE using our estimator, CCTE SP c (semiparametric), CCTE P c (parametric) for Ali-Mikhail-Haq, Clayton and Gumbel copulas for different levels of c. (X 1 , X 2 ) is a vector with Ali-Mikhail-Haq copula with parameter 0.5 and exponentially distributed marginals with parameter 1. Sample size n = 1000, Monte-Carlo iterations M = 300. of the estimations of the CCTE using our estimator, CCTE SP c (semiparametric), CCTE P c (parametric) for Ali-Mikhail-Haq, Clayton and Gumbel copulas for different levels of c. (X 1 , X 2 ) is a vector with Ali-Mikhail-Haq copula with parameter 0.5 and exponentially distributed marginals with parameter 1. Sample size n = 5000, Monte-Carlo iterations M = 300. of the estimations of the CCTE using our estimator, CCTE SP c (semiparametric), CCTE P c (parametric) for Ali-Mikhail-Haq, Clayton and Gumbel copulas for different levels of c. (X 1 , X 2 ) is a vector with Ali-Mikhail-Haq copula with parameter 0.5 and exponentially distributed marginals with parameter 1. Sample size n = 10000, Monte-Carlo iterations M = 300.

Figure 3

 3 Figure 3 Definition of some parameters for the calculation of overtopping.

Figure 4

 4 Figure 4 Representation of the data in three dimensions.

Figure 5

 5 Figure 5 Mean overtopping discharge CCTEc(X, Y ) in function of the height of the costal design h for different condition of risk-sea variables. We take here tan(β) = 0.3 (left) and tan(β) = 0.6 (right).

  .2) So, from Assumption A1, we obtain p n λ( Q εn )

w n r n 1 Lemma 3

 13 Ln(c) Tn -r 1 L(c) Tn p ≤ w n r n 1 Ln(c) Tn -r 1 Ln(c) Tn p + w n r 1 Ln(c)Tn -r 1 L(c) Tn p ≤ w n r n -r p + w n K f 1+ ,λ λ L(c) Tn L n (c) Tn 1+ 1/pand Theorem 1 concludes the proof.Proof of Theorem 4: The proof is a straightforward application of Lemma 1 and Lemma 2.Lemma 1 Under assumptions of Theorem 4, we havep 2(1+ ) n | E Y | X ∈ L(c) Tn -E Y | X ∈ L n (c) Tn | such that f is 1 + integrable.Proof of Lemma 1: Using Theorem 1, we obtainp 2(1+ ) n |P X ∈ L(c) Tn L n (c) Tn | ≤ p 2(1+ ) n d λ L(c) Tn , L n (c) Tn 1+ f 1+ P L(c) Tn -P X ∈ L n (c)g the density of the pair (X, Y ). Under assumptions, and since Y is bounded, we also obtain p 2(1+ ) n y1 L(c) Tn g(x, y)λ(dxdy) -y1 Ln(c) Tn g(x, y)λ(dxdy)≤ p 2(1+ ) n y1 L(c) Tn Ln(c) Tn g(x, y)λ(dxdy) y)1 L(c) Tn Ln(c) Tn λ(dxdy) 1 L(c) Tn Ln(c) Tn λ(dx) L(c) Tn , L n (c) Tn 2(1+ ) f |X ∈ L(c) Tn ] -E[Y |X ∈ L n (c) Tn ] = p 2(1+ ) n y 1 L(c) Tn g(x, y)λ(dxdy) P[X ∈ L(c) Tn ] -1 -y 1 Ln(c) Tn g(x, y)λ(dxdy) P[X ∈ L n (c) Tn ] ∈ L(c) Tn ] P[X ∈ L n (c) Tn ] P[X ∈ L(c) Tn ] y 1 L(c) Tn g(x, y)λ(dxdy)-y 1 Ln(c) Tn g(x, y)λ(dxdy)+ y 1 Ln(c) Tn g(x, y)λ(dxdy) • P[X ∈ L(c) Tn ] -P[X ∈ L n (c) Tn ] .Using (10.3) and (10.4) we obtain the result.Lemma 2 Under assumptions of Theorem 4 and withd n = o ( √ n), we have d n E Y |X ∈ L n (c) Tn -E n Y |X ∈ L n (c) We have d n E Y |X ∈ L n (c) Tn -E n Y |X ∈ L n (c) Tn = d n y1 Ln(c) Tn g(x, y)λ(dxdy) P[X ∈ L n (c) Tn ] -n i=1 Y i 1 {Xi∈Ln(c) Tn } n i=1 1 {Xi∈Ln(c) Tn } .Under assumptions of the Lemma 2 and using Theorem 27.2 in[START_REF] Billingsley | Probability and measure, 3rd edn[END_REF], we obtain thatd n P[X ∈ L n (c) Tn ]-Under assumptions of Theorem 5, we havec N E[Y | X ∈ L n (c) Tn ] -E[r N (X) | X ∈ L n (c) Tn ] c N = o(E |r N (X) -r(X)|).Proof of Lemma 3:We haveE[Y | X ∈ L n (c) Tn ] -E[r N (X) | X ∈ L n (c) Tn ] = Ln(c) Tn (r(x) -r N (x))f (x)dλ(x) P[X ∈ L n (c) Tn ] ≤ K Ln(c) Tn |r(x) -r N (x)| f (x)dλ(x) ≤ R d |r(x) -r N (x)| f (x)dλ(x) = E |r N (X) -r(X)|.Hence the result.Lemma 4 Under assumptions of Theorem 5, we haved n E r N (X)|X ∈ L n (c) Tn -E n r N (X)|X ∈ L n (c) Tn P → n→∞ 0 where d n = o ( √ n).Proof of Lemma 4: We haved n E r N (X)|L n (c) Tn -E n r N (X)|L n (c) Tn = d n Ln(c) Tn r N (x)f (x)λ(dx) P[X ∈ L n (c) Tn ] -n i=1 r N (X i )1 {Xi∈Ln(c) Tn } n i=1 1 {Xi∈Ln(c) Tn } .Under the assumptions of the lemma and using Theorem 27.2 in[START_REF] Billingsley | Probability and measure, 3rd edn[END_REF], we obtain thatd n P[X ∈ L n (c) Tn ] -Proof of1. in Theorem 6: We have r n,aX,bY 1 La,n(c) Tn -r aX,bY 1 La(c) Tn p = R d r n,aX,bY (x)1 La,n(c) Tn -r aX,bY (x)1 La(c) Tn p f aX (x)dx 1/p = R d b p r n,X,Y x a 1 {x/a∈Ln(c) Tn } -r X,Y x a 1 {x/a∈L(c) Tn } p and taking t = x/a we obtain r n,aX,bY 1 La,n(c) Tn -r aX,bY 1 La(c) Tnp = R d b p r n,X,Y (t) 1 Ln(c) Tn -r X,Y (t) 1 L(c) Tn p f X (t) dt 1/p = b r n,X,Y 1 Ln(c) Tn -r X,Y 1 L(c) Tn pHence the result.Proof of 2. inTheorem 6: We have CCTE Tn c,n (aX, bY ) -CCTE Tn c (aX, bY ) = E n bY |aX ∈ L n,a (c) Tn -E bY |aX ∈ L a (c) Tn . Using L a (c) Tn = aL(c) Tn and the assumptions, we obtain CCTE Tn c,n (aX, bY ) -CCTE Tn c (aX, bY ) = 1 b E n Y |aX ∈ aL n (c) Tn -E Y |aX ∈ aL(c) Tn .

Table 2

 2 6 is at least the convergence rate of this CCTE estimation. Estimation of CCTEc(X, Y ) in the case of independent copula and exponentially distributed marginals with parameter 1.

	n		c = 0.1	c = 0.25	c = 0.5	c = 0.7	c = 0.9
		CCTE=0.142 CCTE=0.0792 CCTE=0.0279 CCTE=0.0087 CCTE=0.0009
		RMAE	0.0357	0.0388	0.0725	0.1187	0.2898
	1000	Mean	0.1432	0.0798	0.0281	0.0090	0.0009
		σ	0.0062	0.0038	0.0026	0.0012	0.0003
		RMAE	0.0151	0.0194	0.0317	0.0484	0.1623
	5000	Mean	0.1425	0.0795	0.0278	0.0087	0.0008
		σ	0.0027	0.0021	0.0011	0.0004	0.0002
		RMAE	0.0110	0.0137	0.0211	0.0453	0.0875
	10000	Mean	0.1428	0.0792	0.0279	0.0088	0.0008
		σ	0.0019	0.0013	0.0007	0.0004	9.43 10 -5
			c = 0.1 c = 0.25 c = 0.5	c = 0.7	c = 0.9
		n = 1000	0.0311	0.0339	0.0632	0.1035	0.2526
		n = 5000	0.0151	0.0191	0.0314	0.0479	0.1609
		n = 10000	0.0116	0.0145	0.0223	0.0478	0.0925

Table 3

 3 Approximated βn •RMAE, with βn = n 1/6 ln(n) 4/6 , in the case of independent copula and exponentially distributed marginals with parameter 1.

Table 4 E

 4 L 2 for Ali-Mikhail-Haq copula with parameter 0.5 and exponentially distributed marginals with parameter 1. The bounded regression function is y

	n	c = 0.1	c = 0.25	c = 0.5	c = 0.7	c = 0.9
	n = 1000	0.7736	0.7980	0.7599	0.5817	0.2928
	n = 5000	0.4846	0.5801	0.5047	0.3975	0.1781
	n = 10000	0.4167	0.4883	0.4459	0.3612	0.1679

we illustrate E L2 (i.e. the L 2 -consistency for the estimator r n ) for different level of risk c and different sample size n.

Table 5

 5 Estimation of CCTEc(X, Y ) in the case of Ali-Mikhail-Haq copula with parameter 0.5 and exponentially distributed marginals with parameter 1.

		c = 0.1	c = 0.25 c = 0.5 c = 0.7 c = 0.9
	n = 500	0.0236	0.0169	0.0068	0.0043	0.0009
	n = 1000	0.0161	0.0121	0.0076	0.0030	0.0008
	n = 5000	0.0075	0.0066	0.0035	0.0014	0.0002
	n = 10000	0.0053	0.0057	0.0020	0.0013	0.0002

Table 6

 6 Approximated βn • RMAE, with βn = n 1/10 in the case of Ali-Mikhail-Haq copula with parameter 0.5 and exponentially distributed marginals with parameter 1.

This requires that the partial derivatives of the regression function r are k-Hölderian with a constant C (for further details see Definition 1 in[START_REF] Kohler | Optimal global rates of convergence for nonparametric regression with unbounded data[END_REF]).

Acknowledgements The authors thank the two anonymous reviewers and an associated editor for their useful comments and suggestions. This work has been partially supported by the French research national agency (ANR) under the references ANR-08BLAN-0314-01 and ANR 2011 BS01 010 01 projet Calibration. The authors thank Yannick Baraud, Roland Diel, Christine Tuleau-Malo and Patricia Reynaud-Bouret for fruitful discussions.


est destinée au dépôt

Reducing Manipulability

François Durand, Fabien Mathieu and Ludovic Noirie

CRITERIA FOR VOTING SYSTEMS

Set of voting systems Veto. . .

Informed majority coalition criterion (InfMC)

A majority may choose the outcome when they know the other votes. Borda. . .

Ignorant majority coalition criterion (IgnMC)

A majority may choose the outcome.

Range voting, Approval, Coombs. . .

Majority favorite criterion (MF)

Elects a candidate when she is preferred by a majority.

Plurality, Two-round, IRV, Bucklin. . .

Condorcet criterion (Cond)

Elects the Condorcet winner when there is one.

Baldwin, Dodgson, Kemeny, Maximin, Nanson, Schulze, Tideman. . .

GOAL: MINIMIZE THE MANIPULABILITY RATE

ρ(f ) = P(voting system f is manipulable)

CONDORCIFICATION

A new voting system based on f . Elects the Condorcet winner when there is one.

Otherwise, same outcome as f .

CONDORCIFICATION THM.

If f meets InfMC:

.

SLICE

Yet another voting system based on f . Depends only on orders of preference. Meets the Condorcet criterion.

SLICING THEOREM

If voters are independent:

f cs

Initial voting system f

Condorcification of f

Elects Condorcet winner when she exists. Otherwise, same outcome as f .

Best slice of f c

Depends only on orders of preference. Meets the Condorcet criterion.

CONSEQUENCES

To minimize manipulability while keeping InfMC, one may restrict to voting systems that: Depend only on orders of preference, And meet the Condorcet criterion.


Clement Papadacci
email: papadacci.clement@gmail.com

Mickael Tanter

Mathieu Pernot

Mathias Fink

UltrasoundBackscatter Tensor Imaging (BTI): Analysis of the spatial coherence of ultrasonic speckle in anisotropic soft tissues

Keywords: anisotropy, fiber, myocardium, spatial coherence, plane wave imaging, coherent compounding

The assessment of fiber architecture is of major interest in the progression of myocardial disease. Recent techniques such as Magnetic Resonance (MR) Diffusion Tensor Imaging or Ultrasound Elastic Tensor Imaging (ETI) can derive the fiber directions by measuring the anisotropy of water diffusion or tissue elasticity, but these techniques present severe limitations in clinical setting. In this study, wepropose a new technique, the Backscatter Tensor Imaging (BTI) whichenables determining the fibers directions in skeletal muscles and myocardial tissues, by measuring the spatial coherence of ultrasonic speckle. We compare the results to ultrasound ETI.Acquisitions were performed using a linear transducer array connected to an ultrasonic scanner mounted on a motorized rotation device with angles from 0° to 355° by 5° increments to image ex vivo bovine skeletal muscle and porcine left ventricular myocardial samples. At each angle, multiple plane waves were transmitted and the backscattered echoes recorded. The coherence factor was measured as the ratio of coherent intensity over incoherent intensity of backscattered echoes. In skeletal muscle, maximal/minimal coherence factor was found for the probe parallel/perpendicular to the fibers. In myocardium, the coherence was assessed across the entire myocardial thickness, and the position of maxima and minima varied transmurally due to the complex fibers distribution. In ETI, the shear wave speed variation with the probe angle was found to follow the coherence variation. Spatial coherence can thus reveal the anisotropy of the ultrasonic speckle in skeletal muscle and myocardium. BTI could be used on any type of ultrasonic scanner with rotative phased-array probes or 2-D matrix probes for non-invasive evaluation of myocardial fibers.

I. INTRODUCTION

The structure of the heart wall iscomplex.The orientation of the myofibersvaries smoothly and continuously through the wall thickness [START_REF] Streeter | An engineering analysis of myocardial fiber orientation in pig's left ventricle in systole[END_REF]. This specific fiber architecture is relatedto mechanical [START_REF] Arts | Relating myocardial laminar architecture to shear strain and muscle fiber orientation[END_REF][START_REF] Costa | Laminar fiber architecture and threedimensional systolic mechanics in canine ventricular myocardium[END_REF][START_REF] Waldman | Relation between transmural deformation and local myofiber direction in canine left ventricle[END_REF]and electrical [START_REF] Hooks | Laminar Arrangement of Ventricular Myocytes Influences Electrical Behavior of the Heart[END_REF][START_REF] Kadish | Interaction of fiber orientation and direction of impulse propagation with anatomic barriers in anisotropic canine myocardium[END_REF][START_REF] Roberts | Influence of cardiac fiber orientation on wavefront voltage, conduction velocity, and tissue resistivity in the dog[END_REF][START_REF] Taccardi | Effect of myocardial fiber direction on epicardial potentials[END_REF] properties of the myocardium. Therefore, imaging the fiber architecture in vivois of major interest in the understanding of cardiac function and in the progression of myocardial diseases.

Opticalimaging modalities such as optical coherent tomography [START_REF] Fleming | Quantification of cardiac fiber orientation using optical coherence tomography[END_REF] or two-photon microtomy [START_REF] Huang | Threedimensional cardiac architecture determined by two-photon microtomy[END_REF]can image the microscopic structure but with low penetration capabilities and are thus limited to the mapping of superficial structures. Other imaging tools such as Magnetic Resonance Diffusion Tensor Imaging (MR-DTI) [START_REF] Hsu | Magnetic resonance myocardial fiber-orientation mapping with direct histological correlation[END_REF], [START_REF] Scollan | Histological validation of myocardial microstructure obtained from diffusion tensor magnetic resonance imaging[END_REF]have also been used to map the myocardial fiber structure, but mostly in ex vivo tissues. DTI can also measure the fiber orientationin vivo [START_REF] Reese | Imaging myocardial fiber architecture in vivo with magnetic resonance[END_REF][START_REF] Tseng | Imaging myocardial fiber disarray and intramural strain hypokinesis in hypertrophic cardiomyopathy with MRI[END_REF][START_REF] Wu | Diffusion Tensor Magnetic Resonance Imaging Mapping the Fiber Architecture Remodeling in Human Myocardium After Infarction: Correlation With Viability and Wall Motion[END_REF]but this remains challenging to perform in human patients because of long acquisition times.

In contrast, ultrasonic imaging modalitiesenable real-time visualization of the heart. Echocardiography is routinely employed in clinical practice to examine the heart motionandto provide estimations of the global cardiac function such as cardiac output, derived from the measure of diastolic and systolic ventricular volume [START_REF] Rimington | Echocardiography: A Practical Guide for Reporting[END_REF]. Ultrasound can also be used to characterize the myocardial structure, through the analysis of backscattered ultrasound. The dependence of backscattered intensity and ultrasound attenuation with fiber orientation was investigated extensivelyover the past decades [START_REF] Mottley | Anisotropy of the ultrasonic attenuation in soft tissues: Measurements in vitro[END_REF][START_REF] Baldwin | Measurements of the anisotropy of ultrasonic attenuation in freshly excised myocardium[END_REF][START_REF] Wickline | Three-dimensional characterization of human ventricular myofiber architecture by ultrasonic backscatter[END_REF][START_REF] Madaras | Anisotropy of the ultrasonic backscatter of myocardial tissue: II. Measurements in vivo[END_REF]. Nevertheless,backscattered intensity depends on many parameters including heterogeneity of the material, presence of bright specular echoes,angle view and period in the cardiac cycle.Therefore, clinical implementation of these techniques remains challenging.

Ultrasound and MR elastography techniques have also been proposed for mapping the anisotropy of the elastic properties in fibrous soft tissues [START_REF] Sinkus | Imaging anisotropic and viscous properties of breast tissue by magnetic resonance-elastography[END_REF], [START_REF] Gennisson | Viscoelastic and Anisotropic Mechanical Properties of in vivo Muscle Tissue Assessed by Supersonic Shear Imaging[END_REF]. Because elastic properties of biological tissues are linked to the microstructure organization, elastic anisotropy can provide information on the fiber organization. Shear Wave Imaging (SWI), an ultrasound based technique that can map quantitatively the elastic properties of soft tissues in real-time [START_REF] Bercoff | Supersonic shear imaging: a new technique for soft tissue elasticity mapping[END_REF], was used to determine the fiber directions in skeletal muscle [START_REF] Gennisson | Viscoelastic and Anisotropic Mechanical Properties of in vivo Muscle Tissue Assessed by Supersonic Shear Imaging[END_REF]and myocardium [START_REF] Lee | Mapping myocardial fiber orientation using echocardiography-based shear wave imaging[END_REF][START_REF] Lee | Noninvasive assessment of myocardial anisotropy in vitro and in vivo using Supersonic Shear Wave Imaging[END_REF][START_REF] Lee | Ultrasound elastic tensor imaging: comparison with MR diffusion tensor imaging in the myocardium[END_REF].The propagation of shear waves was measured along several directions to determine the anisotropy of the shear modulus and the fibers direction. SWI has been employed in vivoon the beating heart [START_REF] Couade | In vivo quantitative mapping of myocardial stiffening and transmural anisotropy during the cardiac cycle[END_REF][START_REF] Pernot | Real-time assessment of myocardial contractility using shear wave imaging[END_REF][START_REF] Papadacci | Shear Wave Imaging of the heart using a cardiac phased array with coherent spatial compound[END_REF][START_REF] Song | Robust shear wave motion tracking using ultrasound harmonic imaging[END_REF]to derive myocardial stiffness and contractility estimation [START_REF] Couade | In vivo quantitative mapping of myocardial stiffening and transmural anisotropy during the cardiac cycle[END_REF], [START_REF] Pernot | Real-time assessment of myocardial contractility using shear wave imaging[END_REF] and was used to determine the fiber orientation in vitro in porcine myocardial samples [START_REF] Lee | Mapping myocardial fiber orientation using echocardiography-based shear wave imaging[END_REF], [START_REF] Lee | Noninvasive assessment of myocardial anisotropy in vitro and in vivo using Supersonic Shear Wave Imaging[END_REF] and in vivo in open-chest animals [START_REF] Lee | Mapping myocardial fiber orientation using echocardiography-based shear wave imaging[END_REF], [START_REF] Lee | Ultrasound elastic tensor imaging: comparison with MR diffusion tensor imaging in the myocardium[END_REF]. This technique, named Elastic Tensor Imaging (ETI) was shown to provide comparable information on the fiber direction as DTI in a small region of an ex vivo heart tissue [START_REF] Lee | Ultrasound elastic tensor imaging: comparison with MR diffusion tensor imaging in the myocardium[END_REF].Despite the potential interest of anisotropic elastic estimates for cardiac pathologies screening or diagnosis, an important limitation of this technique for future fiber tracking clinical applications is the need of generating a shear waveat every locations of interest, which require long acquisitions.

Another way to characterize the tissue microstructure with ultrasound, is to analyzethe spatial coherence of the backscattered echoes. Spatial coherence characterizes the level of similarity of backscattered signals received by two different distant elements of an ultrasonic probe. The spatial coherence of backscattered signals provide information on the distribution of scatterers within a focal region. In a random distribution of scatterers, Van Cittert-Zernike showed that spatial coherence of light was dependent only on the size of the focal zone. This result was extended to the field of pulse echo ultrasound by Fink and Mallart [START_REF] Mallart | The van Cittert-Zernike theorem in pulse echo measurements[END_REF]. In the field of diagnostic ultrasound, spatial coherence has been extensively investigated in many applications including aberration corrections [START_REF] Walker | Speckle coherence and implications for adaptive imaging[END_REF][START_REF] Mallart | Adaptive focusing in scattering media through sound-speed inhomogeneities: The van Cittert Zernike approach and focusing criterion[END_REF][START_REF] Hollman | Coherence factor of speckle from a multi-row probe[END_REF][START_REF] Li | Adaptive imaging using the generalized coherence factor[END_REF]and reduction of the clutter signal in B-mode images [START_REF] Lediju | Short-lag spatial coherence of backscattered echoes: imaging characteristics[END_REF], [START_REF] Dahl | Harmonic spatial coherence imaging: an ultrasonic imaging method based on backscatter coherence[END_REF].Additionally, Derode et al. [START_REF] Derode | Spatial coherence of ultrasonic speckle in composites[END_REF] proposed to use the spatial coherence information to characterize the scatterers distribution in composite materials. They estimated the spatial coherence of the backscattered signals on a linear array along different directions and found a higher spatial coherence when the linear array was oriented along the fibers than across them. Therefore, we assumed that spatial coherencecould also reveal the fibers orientation in soft biological tissues, and we introduce here a new technique called Backscatter Tensor Imaging (BTI) to characterize the human tissue microstructure. This approach takes advantage of ultrafast ultrasound imaging [START_REF] Tanter | Ultrafast imaging in biomedical ultrasound[END_REF]based on plane wave transmissions in order to strongly decrease the number of transmissions required to assess the backscatter tensor information. This is mandatoryto ensure a clinical applicability of the technique in a dynamic imaging mode.

In this study, weinvestigated the anisotropy of spatial coherence in fibrous tissues. First, a technique was developed to map rapidly the spatial coherence in 2D using alimited number of plane wave transmits. The anisotropy of the spatial coherence was then evaluated by rotating the ultrasonic probe. Experiments were performed in ex vivobovine muscles and porcine myocardial samples. The results obtained with BTI werecomparedtothe anisotropy of elastic modulus provided by Elastic Tensor Imaging. This technique has the potential of imaging in real-time the three-dimensional (3D) microstructure of fibrous tissues.

II. MATERIAL AND METHODS

A. Spatial coherence in random media and in anisotropic materials 1) Spatial coherence in random media: The Van-Cittert Zernike theorem In conventional ultrasound imaging, an ultrasound pulseis focused in the region of interest (Figure 1(a)). The ultrasound field backscattered by the random distribution of scatterersis received on all the array elements (Figure 1(b)) and appear mostly coherent at large distance from the focus. The spatial coherencecharacterizes the similarity between the signalsreceived by two distant elements of the array.VanCittert-Zernike determined the degree of coherence by defining a coherence function as the averaged cross-correlation between two signals received at two points of space. As the distance between elements increases, the degree of coherence decreases. The so-called Van Cittert-Zernike theorem states that the coherence function is the spatial Fourier transform of the intensity distribution at the focus. Therefore a focus generated by a rectangular aperture provides a triangle coherence function [START_REF] Mallart | The van Cittert-Zernike theorem in pulse echo measurements[END_REF], (Figure 1(c)).

Fig.1. Principle ofspatial coherence assessment in a random media. An ultrasound pulse is focused in the biological tissue with an ultrasonic probe (a). The backscatter echoes are received on the array (b). The coherence function R(m) is computed using cross-correlations performed between pairs of transducer elements distant by m (c).

The coherence function𝑅is assessed as a function of distance in number of elements 𝑚 (or lag)by making auto-correlations between all pairs of receiver elements( [START_REF] Derode | Spatial coherence of ultrasonic speckle in composites[END_REF]):

𝑅 𝑚 = 𝑁 𝑁 -𝑚 𝑐 𝑖, 𝑖 + 𝑚 𝑁-𝑚 𝑖=1 𝑐 𝑖, 𝑖 𝑁 𝑖=1 , ( 1
)
where N is the number of elements of the array and 𝑐 𝑖, 𝑗 is defined as:

𝑐 𝑖, 𝑗 = (𝑆 𝑖 𝑡 -𝑆 𝑖 )(𝑆 𝑗 𝑡 -𝑆 𝑗 ) 𝑇 2 𝑇 1 , (2)
where[𝑇 1 𝑇 2 ] is the temporal window centered on the focal time, 𝑆 𝑖 is the time-delayed signal received on transducer i and 𝑆 𝑘 is defined as follow:

𝑆 𝑘 = 1 𝑇 2 -𝑇 1 𝑆 𝑘 𝑡 𝑇 2 𝑇 1 , (3)
The degree of spatial coherence can also be evaluated by another parameter,the ratio of coherent intensity over incoherent intensity, the so called Mallart-Fink focusing factor orcoherence factor C [START_REF] Mallart | Adaptive focusing in scattering media through sound-speed inhomogeneities: The van Cittert Zernike approach and focusing criterion[END_REF], [START_REF] Hollman | Coherence factor of speckle from a multi-row probe[END_REF], defined as:

𝐶 = 𝑆 𝑖 𝑡 𝑁 𝑖=1 2 𝑇 2 𝑇 1 𝑁 𝑆 𝑖 (𝑡) 2 𝑇 2 𝑇 1 𝑁 𝑖=1 , (4)
2) Coherence in anisotropic media Derode et al.in [START_REF] Derode | Spatial coherence of ultrasonic speckle in composites[END_REF] demonstrated that the coherence function𝑅 varied with the fiber direction in anisotropic composite solid materials (e.g. Figure 2(a)). An ultrasound linear transducer mounted on a rotation devicewas used to evaluate the coherence function along different directions.Across the fibers,the coherence function was found to be a triangle (e.g. Figure 2(b) blue curve) whereas the degree of coherence increased along the fibers (e.g. Figure 2(b) red curve). Indeed, across the fibers, the scatterers appears randomly distributed; but along the fibers, scatterers are distributed preferentially along them which backscatters the transmitted beam more coherently. B. Experimental setup and data acquisition 1) Mapping the spatial coherence using plane wave coherent compounding Spatial coherence evaluation as described above requires to focus an ultrasound pulse within a region of interest. However, soft tissues are highly inhomogeneous in terms of tissue composition which implies the presence of inhomogeneities in spatial coherence. In fact, spatial coherence can be affected by strong specular scatterers. Specular scatterers increase spatial coherence locally and can occur both along the fibers and across the fibers. In such tissues, the coherence function or coherence factor has to be averaged. In order to map the spatial coherence overa larger region, the ultrasonic beam must be focused successively at different locations.However, a large number of focused transmits may strongly affect the frame rate and therefore limit the clinical development of this technique.

In order to reduce the number of transmits for real-time implementation; we propose here to use plane wave coherent compounding. The principle relies on the transmission of several plane waves with the full aperture at different incident angles [START_REF] Montaldo | Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography[END_REF]. Then, by applying coherent summation of the backscattered waves,a focus issynthetized everywhere in the imaging plane (e.g. Figure 3).

By changing the delay applied to each backscattered echoes from (a), the resulting waves can bevirtually focused at different depth (c) and lateral positions (d).

By synthesizinga focus at each point of the imaging plane, the spatial coherence of corresponding backscattered echoes can be evaluated everywhere.Therefore, it allows mapping the spatial coherence of the entire 2D imaging plane using a small number of transmits. A mean coherence degree can be computed over a region of interest (e.g. Figure 4). 3) Fiber angle estimation: Backscatter Tensor Imaging (BTI)

In soft tissues,the coherence factor was averaged around the central axis of the image, at each depth, over a small area (the lateral dimension was 4 mm and the axial dimension was 1 mm). Thisparameter represents the degree of spatial coherence over a region with a given lateral and axial extension and was defined as theSpatially Averaged Coherence Factor (SACF). SACF was estimated at each depth as a function of probe angle. The angle of the fiber was assigned to the probe angle for which the SACF was at its maximum.

4) Fractional anisotropy

For an anisotropic medium the fractional anisotropy (FA)was calculated at each depth as:

𝐅𝐀 = 𝟐 (𝐂 // -𝐂 ) 𝟐 + (𝐂 ⊥ -𝐂 ) 𝟐 𝐂 // 𝟐 + 𝐂 ⊥ 𝟐 , (5)
where 𝐂 // is the value ofSACFalong fibers, 𝐂 ⊥ is the value ofSACFacross fibers and 𝐂 isthe mean SACF over probe angles. FA allows the evaluation of the degree of anisotropy of the backscatter properties.FA quantifies the degree of anisotropy of a medium with respect to a technique.This parameter varies between 0 (isotropic medium) to 1.

C. Elastic Tensor Imaging (ETI) 1) Principle

Elastic tensor imaging was used here as a validated technique to determine the fiber direction in anisotropic tissue. The principle relies on the fact that shear waves propagate the fastest along the fibers and the slowest across the fibers.Lee et al [START_REF] Lee | Mapping myocardial fiber orientation using echocardiography-based shear wave imaging[END_REF] has also demonstrated that shear wave imaging has the sensitivity to detect the complex fiber orientation in myocardium as a function of depth. This technique has been called ETI (elastic tensor imaging). ETI has also been correlated with DTI (diffusor tensor imaging) in the myocardium [START_REF] Lee | Ultrasound elastic tensor imaging: comparison with MR diffusion tensor imaging in the myocardium[END_REF]. We propose in this study to use ETI to validate results from BTI.

2) Acquisition

ETI acquisitions were performed during the BTI experiments using the same probe. At each probe angle, just after BTI acquisition, Shear Wave Imaging was performed.40 frames were acquired to image the propagation of a shear wave induced by the acoustic radiation force of a focused ultrasound burst (duration 300µs, focused at 20mm depth for the bovine muscles and 30 mmfor the myocardial samples with anFnumber of 1.5) at 8,000 frames per second.

3) Fiber angle estimation Tissue velocities were obtained using a per pixel frame to frame 1D cross-correlation on demodulated IQ images with an axial kernel size of 3 pixels (385 µm) to obtain images of tissue frame-to-frame axial displacements. At each depth, shear wave speedswere estimated by tracking the maxima of tissue velocities along time in a small region laterally around the focal spot (4mm). Shear wave speeds were evaluatedfor each probe angle.At each depth, the fiber orientation was derived by detecting the maximum shear wave speed with respect to probe angles.For each sample, Spearman's rank correlation coefficient (𝜌)on fiber orientation estimated by ETI and BTI was calculated to compare the two techniques.

III. RESULTS

A. Imaging phantom Spatial coherence using plane wave coherent compounding was evaluated in the imaging phantom as a function of the number of plane wave transmits.In the top row ofFigure 5,normalized B-mode images are shown for increasing number of plane waves (e.g. (a),(c),(e),(g)). Improvement of both contrast and resolution is noticed by increasing the number of plane wave transmits as shown by [START_REF] Montaldo | Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography[END_REF].Signal to noise ratio (SNR) is estimated as described above (part. B.2.a) and averaged on the small and the large inclusions (e.g. Figure 6(b)). The blue curve on Figure 6(a), shows an improvement of about 10dB when 41 plane wavesare used (e.g. Figure 5(g)) compared to 1 (e.g. Figure 5(a)).

Coherence functions are calculated across the fibers (c)-blue curve and along the fibers (c)-red curve. The distance (transducer elements) varies with depth z, and is equal to z/F-number. It was normalized with

respect to a F-number of 1.5.

The SACF was then calculated as described in part II-B-3. It was calculated on the central axis for each probe angle. Figure 7(a) shows an example of SACF variation at a 20.5mm depth as a function of probe angle. The maxima of the curve give the parallel direction to the fibers (blue arrow) whereas minima give the transverse direction to the fibers (green arrow).

Fig.7. SACF variation with probe angle in skeletal bovine muscle,at a depth of 20.5 mm (a). The maximal SACF is obtained when the probe is parallel to the fibers (blue arrow). The minimal SACF is obtained when the probe is perpendicular to the fibers (green arrow). Fiber direction can be assessed at each depth (b).

The SACF variation was assessed as a function of probe angle at each depth and the fiber direction was assigned to the probe angle for with maximal SACF. As expected the fiber direction was found to be nearly the same as a function of depth.

The fractional anisotropy was calculated, as described previously, for all samples and averaged: FA = 0.46 (+/-0.05).

2) Validation of BTI against ETI

The fiber direction wasalso measured using ETI. The procedure is described in part II-C. At each probe angle and each depth,the shear wave speed (SWS)was assessed in a small region laterally around the focal spot. It was found to varyin average from1.9 (+/-1) m/s across the fibers up to 4.1 (+/-1) m/s along the fibers with FA = 0.42 (+/-0.07). Figure 8 shows an example at a 20.5mm depth. A good agreement was found between estimations of fiber directions made with BTI and with ETI. The agreement was quantified by a Spearman correlation between ETI and BTI on the estimations of fiber directions over depth. In average, over the three samples, the correlation coefficient was 𝜌 = 0.92 ± 0.04.

C. Inex vivo myocardial samples 1) BTI measurement

The myocardium is a much more complex anisotropic soft tissue compared to skeletal muscle. The fibers are not visible on focused B-mode images (e.g. Figure 9(a)) and their orientations change transmurally [START_REF] Streeter | An engineering analysis of myocardial fiber orientation in pig's left ventricle in systole[END_REF].The coherence function was shown to vary with fiber angle through the myocardium. For each depth, the direction of fibers was assessed by tracking the maxima of the coherence factor. Then corresponding parallel and perpendicular coherence functions were averaged along depth. Figure 9(b) displays averagedcoherence function when the probe is set along the fibers (red curve) and across the fibers (blue curve).The coherence function is higher along the fibers (red curve) than across the fibers (blue curve). SACFwas also found to vary with probe angle (Figure 10(b)). The maxima and minima of SACFwere also shown to vary with depth. Figure 10(b) displays an example of SACF variation in one myocardial sample as a function of probe angle and at three different depths. The red curve represents the SACF at sub-epicardium (e.g. Figure 10(a) red rectangle), the green curve at midwall (e.g. Figure 10(a) green rectangle) and the blue curve at sub-endocardium (e.g. Figure 10(a) blue rectangle). The SACF variation can be imaged as a function of probe angle and depth through the myocardial wall (wall thickness) (e.g. Figure 10(c)). This maxima shift gives a different fiber orientation as a function of depthas expected. The probe angle that corresponds to the maximum SACFwas estimated at each depth. The fiber directions estimated with BTI, over the three samples, are represented in figure 11(a),(c),(e). Fibers were found to vary transmurally, from approximately -55° near the epicardium to +60° near the endocardium. The fractional anisotropy was also calculated over the 3 samples and found to be: FA = 0.42 (+/-0.05).

2) Validation of BTI against ETI

The fiber direction was also obtained by ETI. The shear wave speed was found to vary from 3.1m/s across the fibers to 5.9m/s along the fibers. The shear wave speed variation with the probe angle was found to follow mostly the spatial coherencevariation at each depth.

The estimation of fiber directionprovided by the two methods were also compared on the three myocardial samples (e.g. Figure 11(b),(d),(f)) and show a good agreement. Spearman correlation between ETI and BTI on the estimation of fiber directions over all porcine myocardial samples was calculated. We found in average: 𝜌 = 0.89 ± 0.03.

3) BTI measurement witha small number of transmitted plane waves.

With the same set of data, SACF was calculated for different number of plane waves used in the coherent compounding process. The results show a degradation of fiber angle estimation as the number of plane waves decrease but the anisotropy is observed in all acquisitions (𝜌 ≥ 0.6). A number as low as 5 transmits could be enough to obtain the estimation of the fiber directions (𝜌 ≈ 0.7).

This may enable us to increase the acquisition frame rate for real-time in vivo experiments and to reduce artifacts from coherent compound due to high tissue velocity.

IV. DISCUSSION

In this paper, spatial coherence of backscattered echoes wasanalyzed to investigate the anisotropy of soft tissues. As previously found by Derode et al. [START_REF] Derode | Spatial coherence of ultrasonic speckle in composites[END_REF]in solid composite materials,we observeda strong anisotropy of spatial coherence in skeletal muscles as well as in myocardium. Based on this anisotropy, the direction of fibers of the skeletal muscle samples was successfully derived and significant FA were found.

Moreover, the complex fiber distribution of myocardial samples was measured in good agreement with ETI, Another original aspect of the method developed in this paper was the use of plane wave transmits to image the spatial coherence. Plane wave coherent compounding was used to synthetize a focus everywhere in the imaging plane in order to assess the coherence functions and the coherence factor at each point of space. In contrast to conventional transmit focusing, plane wave compounding presents several advantages for mapping the spatial coherence. First, a lower number of transmits is required to map the spatial coherence in 2D over the entire field of view. In this study, 41 plane waves were used but the number of plane waves could be reduced to increase the frame rate (e.g. Figure 12). When decreasing the number of transmits, the overall spatial coherence decreases because a larger focus are generated but the contrast of the coherence image remains high.

The results presented in this paper,are limited to the 1D mapping of the fiber distributionalong the rotation axis of the probe. Mapping the fiber distribution in 3D would require to mechanically scan the probe over the sample. However, using 2D matrix arrays the spatial coherence could be assessed in different directions using the same backscattered data, removing the need of rotation axis. Moreover, the plane wave approach could be extended in 3D which would provide the backscatter information over an entire volume. Thus, it may be possible with BTI to in real-time 3D mapsof the microstructure of tissues. The fiber orientation detection is also limited to projected component of real fiber direction on the transducer surface plane. As long as the fibers are oriented parallel to the transducer surface plane, the real direction can be detected. However, to detect fibers with a certain angle regarding to the transducer surface plane, a

Beamforming with subaperture strategy could be used [START_REF] Tanter | Ultrafast compound imaging for 2-D motion vector estimation: application to transient elastography[END_REF].

In this study, a linear transducer with a central frequency of 6MHz was used. For cardiac applications with smaller aperture such as in transthoracic and transoesophagial echocardiography, plane wave coherent compounding could be replaced by diverging waves [START_REF] Karaman | Synthetic aperture imaging for small scale systems[END_REF], [START_REF] Papadacci | High-contrast ultrafast imaging of the heart[END_REF]. The application to transoesophagial echocardiography is particularly interesting because most of the 2D transoesophagial transducers are mounted on a motorized rotation axis.

BTI could be implemented on any type of ultrasonic scanner with conventional probes and adapted sequences.Moreover, beyond cardiac applications, BTI could be used to map the anisotropy of fibrous tissues such as brain in transfontanellar imaging, skeletal muscle or tendon.

V. CONCLUSION

In this paper, spatial coherence of ultrasonic backscatter echoes from plane wave transmits was analyzed in anisotropic tissues. Plane wave coherent compounding enabled mapping the spatial coherence in 2D using a small number of emissions. We investigated the spatial coherence dependence with the probe orientation and showed the possibility to map the fibers distribution in soft tissues. This new technique called Backscatter Tensor Imaging (BTI),was used tomeasure the fibers orientation in (N=3) ex vivo bovine skeletal muscles and in (N=3) ex vivo porcine myocardial samples and compared it to the anisotropy measured by shear wave elastography. We demonstrated that BTI has the sensitivity to reveal the complex transmural fiber distribution inex vivo myocardium. Finally, the spatial coherence dependence with the number of plane waves was investigated and a small number of transmitted plane waves was found sufficient to measure the tissue anisotropy and fiber direction. BTI has a strong potential to become a major tool to explore non-invasively the anisotropy in vivo in human hearts or in other anisotropic organs such including brain, tendon and skeletal muscles.

Fig. 2 .

 2 Fig.2. Principles of spatial coherence in anisotropic media. Along the fibers (red curve), backscattered

Fig. 3 .

 3 Fig.3. Three tilted plane waves are sent successively by an ultrasonic linear transducer. Each plane wave is

Fig. 4 .

 4 Fig.4. Plane wave coherent compounding allows spatial coherence estimationeverywhere in the imaging

Fig. 5 .

 5 Fig.5. Comparison between B-mode images (a),(c),(e),(g) and Coherence B-mode images (b),(d),(f),(h), for

Fig. 6 .

 6 Fig.6. 41 plane waves B-mode image (a) and Coherence B-mode image (b) of a bovine skeletal muscle.

Fig. 8 .

 8 Fig.8.Validation of BTI against ETI at a depth of 20.5 mm (black line on (b)). The SACF is represented for

Fig. 9 .

 9 Fig.9. Focused B-mode image of a porcine myocardial sample embedded in a gel (a). Coherence

Fig. 10 .

 10 Fig.10. Three locations on the focused B-mode image of a myocardial sample (a) and the associated

Fig. 11 .

 11 Fig.11. The comparison of the fiber orientation estimated by BTI and ETI. The BTI fiber orientation plots

Figure 12

 12 displays the results of SACF for 1,5,11,21 ((a),(b),(c),(d) respectively) compared to 41 plane waves (e) in a myocardial sample at each depth (22 mm to 36 mm).To compare the results, a normalization was applied to set the maximum of coherence to 1 and minimum to 0.5. Spearman correlation was performed to compare the fiber angles derived with ETI compared to BTI for the 5 different configurations. Spearman correlation's coefficient was plotted against the number of plane wave use in the coherent summation (e.g. Figure12(f)).

Fig. 12 .

 12 Fig.[START_REF] Scollan | Histological validation of myocardial microstructure obtained from diffusion tensor magnetic resonance imaging[END_REF]. Normalized SACF as a function of probe angle for different number of plane waves[START_REF] Streeter | An engineering analysis of myocardial fiber orientation in pig's left ventricle in systole[END_REF][START_REF] Hooks | Laminar Arrangement of Ventricular Myocytes Influences Electrical Behavior of the Heart[END_REF][START_REF] Hsu | Magnetic resonance myocardial fiber-orientation mapping with direct histological correlation[END_REF][START_REF] Sinkus | Imaging anisotropic and viscous properties of breast tissue by magnetic resonance-elastography[END_REF][START_REF] Tanter | Ultrafast compound imaging for 2-D motion vector estimation: application to transient elastography[END_REF]

  another technique previously developed to map the elastic anisotropy of myocardium. Although the two techniques rely on the anisotropy of different physical properties (i.e backscatter and elastic properties), this result is linked in both cases to the microstructure of the tissue, and provided comparable results in term of fiber orientation. ETI was already compared in ex vivo myocardial samples to magnetic resonance diffusion tensor imaging (MR-DTI), which measures the anisotropy of water diffusion, and similar fiber distribution was found by both techniques. Current clinical uses of MR DTI are strongly limited by long acquisition times, cost and complexity of MR scanners, so that ultrasound based techniques such as BTI and ETI may have strong clinical potential for real-time measurements in the myocardium. One advantage of BTI over ETI is the lower energy required to map an entire region because ETI relies on the generation of shear waves at multiple locations using the acoustic radiation force induced by long ultrasonic burstswhereas BTI relies only on pulse emissions.

Acknowledgements

The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement n°311025