Aida Toma
email: aida-toma@yahoo.com

Michel Broniatowski
email: michel.broniatowski@upmc.fr

Dual divergence estimators and tests: robustness results

Keywords: Location model, minimum divergence estimators, robust estimation, robust test, scale model

The class of dual φ-divergence estimators (introduced in Broniatowski and Keziou (2009) [6]) is explored with respect to robustness through the influence function approach. For scale and location models, this class is investigated in terms of robustness and asymptotic relative efficiency. Some hypothesis tests based on dual divergence criterions are proposed and their robustness properties are studied. The empirical performances of these estimators and tests are illustrated by Monte Carlo simulation for both noncontaminated and contaminated data.

Introduction

Minimum divergence estimators and related methods have received considerable attention in statistical inference because of their ability to reconcile efficiency and robustness. Among others, Beran [START_REF] Beran | Minimum Hellinger distance estimates for parametric models[END_REF], Tamura and Boos [START_REF] Tamura | Minimum Hellinger distance estimation for multivariate location and covariance[END_REF], Simpson [START_REF] Simpson | Minimum Hellinger distance estimation for the analysis of count data[END_REF][START_REF] Simpson | Hellinger deviance tests: efficiency, breakdown points, and examples[END_REF] and Toma [START_REF] Toma | Minimum Hellinger distance estimators for multivariate distributions from Johnson system[END_REF] proposed families of parametric estimators minimizing the Hellinger distance between a nonparametric estimator of the observations density and the model. They showed that those estimators are both asymptotically efficient and robust. Generalizing earlier work based on the Hellinger distance, Lindsay [START_REF] Lindsay | Efficiency versus robustness: the case of minimum Hellinger distance and related methods[END_REF], Basu and Lindsay [START_REF] Basu | Minimum disparity estimation for continuous models[END_REF], Morales et al. [START_REF] Morales | Asymptotic divergence of estimates of discrete distribution[END_REF] have investigated minimum divergence estimators, for both discrete and continuous models. Some families of estimators based on approximate divergence criterions have also been considered; see Basu et al. [START_REF] Basu | Robust and efficient estimation by minimising a density power divergence[END_REF].

Broniatowski and Keziou [START_REF] Broniatowski | Parametric estimation and tests through divergences and duality technique[END_REF] have introduced a new minimum divergence estimation method based on a dual representation of the divergence between probability measures. Their estimators are defined in an unified way for both continuous and discrete models. They do not require any prior smoothing and include the classical maximum likelihood estimators as a benchmark. A special case for the Kullback-Leibler divergence is presented in Broniatowski [START_REF] Broniatowski | Estimation through Kullback-Leibler divergence[END_REF]. The present paper presents robustness studies for the classes of estimators generated by the minimum dual φ-divergence method, as well as for some tests based on corresponding estimators of the divergence criterion.

We give general results that allow to identify robust estimators in the class of dual φ-divergence estimators. We apply this study for the Cressie-Read divergences and state explicit robustness results for scale models and location models. Gain in robustness is often paid by some loss in efficiency. This is discussed for some scale and location models. Our main remarks are as follows. All the relevant information pertaining to the model and the true value of the parameter to be estimated should be used in order to define, when possible, robust and nearly efficient procedures. Some models allow for such procedures. The example provided by the scale normal model shows that the choice of a good estimation criterion is heavily dependent on the acceptable loss in efficiency in order to achieve a compromise with the robustness requirement. When sampling under the model is overspread (typically for Cauchy and logistic models), non surprisingly the maximum likelihood estimator is both efficient and robust and therefore should be prefered (see subsection 3.2).

On the other hand, these estimation results constitute the premises to construct some robust tests. The purpose of robust testing is twofold. First, the level of a test should be stable under small arbitrary departures from the null hypothesis (i.e. robustness of validity). Second, the test should have a good power under small arbitrary departures from specified alternatives (i.e. robustness of efficiency). To control the test stability against outliers in the aforementioned senses, we compute the asymptotic level of the test under a sequence of contaminated null distributions, as well as the asymptotic power of the test under a sequence of contaminated alternatives. These quantities are seen to be controlled by the influence function of the test statistic. In this way, the robustness of the test is a consequence of the robustness of the test statistic based on a dual φ-divergence estimator. In many cases, this requirement is met when the dual φ-divergence estimator itself is robust.

The paper is organized as follows: in Section 2 we present the classes of es-timators generated by the minimum dual φ-divergence method. In Section 3, for these estimators, we compute the influence functions and give the Fisher consistency. We particularize this study for the Cressie-Read divergences and state robustness results for scale models and location models. Section 4 is devoted to hypothesis testing. We give general convergence results for contaminated observations and use it to compute the asymptotic level and the asymptotic power for the tests that we propose. In Section 5, the performances of the estimators and tests are illustrated by Monte Carlo simulation studies. In Section 6 we shortly presents a proposal for the adaptive choice of tuning parameters.

2 Minimum divergence estimators

Minimum divergence estimators

Let ϕ be a non-negative convex function defined from (0, ∞) onto [0, ∞] and satisfying ϕ(1) = 0. Also extend ϕ at 0 defining ϕ(0) = lim x↓0 ϕ(x). Let (X , B) be a measurable space and P be a probability measure (p.m.) defined on (X , B). Following Rüschendorf [START_REF] Rüschendorf | On the minimum discrimination information theorem[END_REF], for any p.m. Q absolutely continuous (a.c.) w.r.t. P , the φ-divergence between Q and P is defined by

φ(Q, P ) := ϕ dQ dP dP. ( 1
)
When Q is not a.c. w.r.t. P, we set φ(Q, P ) = ∞. We refer to Liese and Vajda [START_REF] Liese | On divergences and informations in statistics and information theory[END_REF] for an overview on the origin of the concept of divergence in Statistics.

A commonly used family of divergences is the so-called "power divergences", introduced by Cressie and Read [START_REF] Cressie | Multinomial goodness of fit tests[END_REF] and defined by the class of functions

x ∈ R * + → ϕ γ (x) := x γ -γx + γ -1 γ(γ -1) (2)
for γ ∈ R \ {0, 1} and ϕ 0 (x) :

= -log x + x -1, ϕ 1 (x) := x log x -x + 1 with ϕ γ (0) = lim x↓0 ϕ γ (x), ϕ γ (∞) = lim x→∞ ϕ γ (x)
, for any γ ∈ R. The Kullback-Leibler divergence (KL) is associated with ϕ 1 , the modified Kullback-Leibler (KL m ) to ϕ 0 , the χ 2 divergence to ϕ 2 , the modified χ 2 divergence (χ 2 m ) to ϕ -1 and the Hellinger distance to ϕ 1/2 .

Let {P θ : θ ∈ Θ} be some identifiable parametric model with Θ a subset of R d . Consider the problem of estimation of the unknown true value of the parameter θ 0 on the basis of an i.i.d. sample X 1 , . . . , X n with p.m. P θ 0 . When all p.m. P θ share the same finite support S which is independent upon the parameter θ, the φ-divergence between P θ and P θ 0 has the form φ(P θ , P θ 0 ) = j∈S ϕ P θ (j) P θ 0 (j) P θ 0 (j).

In this case, Liese and Vajda [START_REF] Liese | Convex Statistical Distances[END_REF], Lindsay [START_REF] Lindsay | Efficiency versus robustness: the case of minimum Hellinger distance and related methods[END_REF] and Morales et al. [START_REF] Morales | Asymptotic divergence of estimates of discrete distribution[END_REF] investigated the so-called "minimum φ-divergence estimators" (minimum disparity estimators in Lindsay [START_REF] Lindsay | Efficiency versus robustness: the case of minimum Hellinger distance and related methods[END_REF]) of the parameter θ 0 defined by

θ n := arg inf θ∈Θ φ(P θ , P n ), (3)
where φ(P θ , P n ) is the plug-in estimator of φ(P θ , P θ 0 )

φ(P θ , P n ) = j∈S ϕ P θ (j) P n (j) P n (j),
P n being the empirical measure associated to the sample. The interest on these estimators is motivated by the fact that a suitable choice of the divergence may leads to an estimator more robust than the maximum likelihood one (see also Jiménez and Shao [14]). For continuous models, the estimators in (3) are not defined. Basu and Lindsay [START_REF] Basu | Minimum disparity estimation for continuous models[END_REF], among others, proposed smoothed versions of (3) in this case.

In the following, for notational clearness we write φ(α, θ) for φ(P α , P θ ) for α and θ in Θ. We assume that for any θ ∈ Θ, P θ has density p θ with respect to some dominating σ-finite measure λ.

The divergence φ(α, θ 0 ) can be represented as resulting from an optimization procedure. This result has been obtained independently by Liese and Vajda [START_REF] Liese | On divergences and informations in statistics and information theory[END_REF] and Broniatowski and Keziou [START_REF] Broniatowski | Minimization of φ-divergences on sets of signed measures[END_REF] who called it the dual form of a divergence, due to its connection with convex analysis. Assuming the strict convexity and the differentiability of the function ϕ, it holds

ϕ(t) ≥ ϕ(s) + ϕ ′ (s)(t -s) (4)
with equality only for s = t. Let α and θ 0 be fixed and put t = p α (x)/p θ 0 (x) and s = p α (x)/p θ (x) in ( 4) and then integrate with respect to P θ 0 . This gives

φ(α, θ 0 ) = ϕ p α p θ 0 dP θ 0 = sup θ∈Θ m(θ, α)dP θ 0 (5) with m(θ, α) : x → m(θ, α, x) and m(θ, α, x) := ϕ ′ p α p θ dP α -ϕ ′ p α p θ (x) p α p θ (x) -ϕ p α p θ (x) . (6)
The supremum in ( 5) is unique and is attained in θ = θ 0 , independently upon the value of α. Naturally, a class of estimators of θ 0 , called "dual φ-divergence estimators" (DφE's), is defined by

θ n (α) := arg sup θ∈Θ m(θ, α)dP n , α ∈ Θ. (7)
Formula ( 7) defines a family of M-estimators indexed by some instrumental value of the parameter α and by the function ϕ defining the divergence. The choice of α appears as a major feature in the estimation procedure. Its value is strongly dependent upon some a priori knowledge on the value of the parameter to be estimated. In some examples in subsection 3.2, it even appears that a sharp a priori knowledge on the order of θ 0 leads to nearly efficient and robust estimates. This plays in favor of using the available information pertaining to the model and the data. Section 6 shortly presents some proposal for the adaptive choice of α.

For each α ∈ Θ, the divergence φ(P α , P θ 0 ) between P α and P θ 0 is estimated by

φ n (α, θ 0 ) := m( θ n (α), α)dP n = sup θ∈Θ m(θ, α)dP n . (8)
Further, since inf

α∈Θ φ(α, θ 0 ) = φ(θ 0 , θ 0 ) = 0,
and since the infimum in the above display is unique due to the strict convexity of ϕ, a natural definition of estimators of θ 0 , called "minimum dual φ-divergence estimators" (MDφE's), is provided by

α n := arg inf α∈Θ φ n (α, θ 0 ) = arg inf α∈Θ sup θ∈Θ m(θ, α)dP n . (9)
The DφE's enjoy the same invariance property as the maximum likelihood estimator does. Invariance with respect to a reparametrization (one to one transformation of the parameter space) holds with direct substitution in [START_REF] Broniatowski | An estimation method for Neyman chisquare divergence with application to test of hypothesis[END_REF]. Also, consider a one to one differentiable transformation of the observations, say Y = T (X) and the Jacobian J(x) = d dx T (x). Let θ n (α) defined in (7), based on the X i 's. Let f θ (y) denote the density of the transformed variable Y and θ * n (α) be the DφE based on the Y i 's in the transformed model (with the same parameter θ). Specifically,

θ * n (α) = arg sup θ∈Θ ϕ ′ f α f θ (y) f α (y)dy - 1 n n i=1 ϕ ′ f α f θ (Y i ) f α f θ (Y i ) -ϕ f α f θ (Y i ) . Since f θ (y) = p θ (T -1 (y))|J(T -1 (y))| -1
for all θ ∈ Θ, it follows that θ * n (α) = θ n (α), which is to say that the DφE's are invariant estimators under any regular transformation of the observation space. The same invariance properties hold for MDφE's.

Broniatowski and Keziou [START_REF] Broniatowski | Parametric estimation and tests through divergences and duality technique[END_REF] have proved both the weak and the strong consistency, as well as the asymptotic normality for the estimators θ n (α) and α n . In the next sections, we study robustness properties for these classes of estimators and robustness of some tests based on dual φ-divergence estimators.

Some comments on robustness

The special form of divergence based estimators to be studied in this paper leads us to handle robustness characteristics through the influence function approach. An alternative and appealing robustness analysis in the minimum divergence methods is provided by the Residual Adjustment Function (RAF) (introduced in Lindsay [START_REF] Lindsay | Efficiency versus robustness: the case of minimum Hellinger distance and related methods[END_REF]), which explains the incidence of non typical Pearson residuals, corresponding to over or sub-sampling, in the stability of the estimates. This method is quite natural for finitely supported models. In the case when the densities in the model are continuous, the Pearson residuals are estimated non parametrically which appears to cause quite a number of difficulties when adapted to minimum dual divergence estimation. This motivates the present choice in favor of the influence function approach.

Let α be fixed. For the Cressie-Read divergences, the equation whose solution is θ n (α) defined by [START_REF] Broniatowski | An estimation method for Neyman chisquare divergence with application to test of hypothesis[END_REF] is

- p α p θ γ ṗθ dλ + 1 n n i=1 p α (X i ) p θ (X i ) γ ṗθ (X i ) p θ (X i ) = 0, ( 10
)
where ṗθ is the derivative with respect to θ of p θ . Starting from the definition given by [START_REF] Broniatowski | An estimation method for Neyman chisquare divergence with application to test of hypothesis[END_REF], this equation is obtained by equalizing with zero the derivative with respect to θ of m(θ, α)dP n .

Let x be some outlier. The role of x in [START_REF] Dell'aquilla | Robust tests of predictive accuracy[END_REF] is handled in the term

p α (x) p θ (x) γ ṗθ (x) p θ (x) . ( 11
)
The more stable this term, the more robust the estimate. In the classical case of the maximum likelihood estimator (which corresponds to θ n (α) with γ = 0 and independent on α), this term writes as ṗθ (x) p θ (x) which is the likelihood score function associated to x. It is well known that, for most models, this term is usually unbounded when x belongs to R, saying that the maximum likelihood estimator is not robust. In this respect, [START_REF] Hampel | Robust Statistics: the Approach Based on Influence Functions[END_REF] appears as a weighted likelihood score function. In our approach, for several models, such as the normal scale, [START_REF] Hampel | Robust Statistics: the Approach Based on Influence Functions[END_REF] is a bounded function of x, although ṗθ (x) p θ (x) itself is not. Thus, in estimating equation [START_REF] Dell'aquilla | Robust tests of predictive accuracy[END_REF], the score function is downweighted for large observations. The robustness of θ n (α) comes as a downweight effect of the quantity ṗθ (X i ) p θ (X i ) through the multiplicative term pα(X i ) p θ (X i ) γ which depends on the choice of the divergence. This choice is dictated by the form of pα(x)

p θ (x)
for large x and α fixed. For the models we'll consider as examples, for large x and α fixed, the quantity pα(x) p θ (x) can be large, close to zero, or close to one. Then we appropriately choose γ to be negative, respectively positive in order to obtain the downweight effect. In the next section we study in detail these robustness properties by the means of the influence function.

Some alternative choice has been proposed in literature. Basu et al. [START_REF] Basu | Robust and efficient estimation by minimising a density power divergence[END_REF] proposed to alter the likelihood score factor by the multiplicative term p β θ (x), where β > 0. This induces an estimating procedure which is connected to the minimization of a density power divergence. Both their approach and the present one are adaptive in the sense that the downweight likelihood score factor is calibrated on the data.

Robustness as handled in the present paper is against the bias due to the presence of very few outliers in the data set. Bias due to misspecification of the model is not considered. It has been observed that DφE's are biased under misspecification even in simple situations (for example when estimating the mean in a normal model with assumed variance 1, whereas the true variance is not 1); see Broniatowski and Vajda [START_REF] Broniatowski | Several applications of divergence criteria in continuous families[END_REF]; similar bias are unavoidable in parametric inference and can only be reduced through adaptive specific procedures, not studied here. For alternative robust M-estimation methods using divergences we refer to Toma [START_REF] Toma | Optimal robust M-estimators using divergences[END_REF].

3 Robustness of the estimators

Fisher consistency and influence functions

In order to measure the robustness of an estimator it is common to compute the influence function of the corresponding functional.

A map T which sends an arbitrary probability measure into the parameter space is a statistical functional corresponding to an estimator T n of the parameter θ whenever

T (P n ) = T n .
This functional is called Fisher consistent for the parametric model

{P θ : θ ∈ Θ} if T (P θ ) = θ, for all θ ∈ Θ.
The influence function of the functional T in P measures the effect on T of adding a small mass at x and is defined as

IF(x; T, P ) = lim ε→0 T ( P εx ) -T (P ) ε (12)
where P εx = (1ε)P + εδ x and δ x is the Dirac measure putting all its mass at x.

The gross error sensitivity measures approximately the maximum contribution to the estimation error that can be produced by a single outlier and is defined as sup

x IF(x; T, P ) .

Whenever the gross error sensitivity is finite, the estimator associated with the functional T is called B-robust.

Let X 1 , . . . , X n be an i.i.d. sample with p.m. P .

Let α be fixed and consider the dual φ-divergence estimators θ n (α) defined in [START_REF] Broniatowski | An estimation method for Neyman chisquare divergence with application to test of hypothesis[END_REF]. The functional associated to an estimator θ n (α) is

T α (P ) := arg sup θ∈Θ m(θ, α, y)dP (y). (13)
The functional T α is Fisher consistent. Indeed, the function θ → m(θ, α)dP θ 0 has a unique maximizer θ = θ 0 . Therefore T α (P θ ) = θ, for all θ ∈ Θ.

We denote m ′ (θ, α) = ∂ ∂θ m(θ, α) the d-dimensional column vector with entries ∂ ∂θ i m(θ, α) and m ′′ (θ, α) the d × d matrix with entries ∂ 2 ∂θ i ∂θ j m(θ, α).

In the rest of the paper, for each α, we suppose that the function θ → m(θ, α) is twice continuously differentiable and that the matrix m ′′ (θ 0 , α)dP θ 0 exists and is invertible. We also suppose that, for each α, all the partial derivatives of order 1 and 2 of the function θ → m(θ, α) are respectively dominated on some neighborhoods of θ 0 by P θ 0 -integrable functions. This justifies the subsequent interchanges of derivation with respect to θ and integration.

Proposition 1 The influence function of the functional T α corresponding to an estimator θ n (α) is given by

IF(x; T α , P θ 0 ) = m ′′ (θ 0 , α)dP θ 0 -1 ϕ ′′ p α p θ 0 p α p 2 θ 0 ṗθ 0 dP α - -ϕ ′′ p α p θ 0 (x) p 2 α (x) p 3 θ 0 (x) ṗθ 0 (x) . Particularizing α = θ 0 in Proposition 1 yields IF(x; T θ 0 , P θ 0 ) = - m ′′ (θ 0 , θ 0 )dP θ 0 -1 ϕ ′′ (1) ṗθ 0 (x) p θ 0 (x)
and taking into account that

- m ′′ (θ 0 , θ 0 )dP θ 0 -1 = 1 ϕ ′′ (1) I -1 θ 0 it holds IF(x; T θ 0 , P θ 0 ) = I -1 θ 0 ṗθ 0 (x) p θ 0 (x) (14)
where I θ 0 is the information matrix

I θ 0 = ṗθ 0 ṗt θ 0 p θ 0 dλ.
We now look at the corresponding estimators of the φ-divergence. For fixed α, the divergence φ(P α , P ) between the probability measures P α and P is estimated by [START_REF] Broniatowski | Several applications of divergence criteria in continuous families[END_REF]. The statistical functional associated to φ n (P α , P θ 0 ) is

U α (P ) := m(T α (P ), α, y)dP (y). ( 15
)
The functional U α has the property that U α (P θ ) = φ(α, θ), for any θ ∈ Θ. Indeed, using the fact that T α is a Fisher consistent functional,

U α (P θ ) = m (T α (P θ ) , α, y)dP θ (y) = m(θ, α, y)dP θ (y) = φ(α, θ)
for all θ ∈ Θ.

Proposition 2

The influence function of the functional U α corresponding to the estimator φ n (P α , P ) is given by

IF(x; U α , P θ 0 ) = -φ(α, θ 0 ) + m(θ 0 , α, x). (16)
For a minimum dual φ-divergence estimator α n defined in ( 9), the corresponding functional is

V (P ) := arg inf α∈Θ U α (P ) = arg inf α∈Θ m(T α (P ), α, y)dP (y). ( 17
)
The statistical functional V is Fisher consistent. Indeed,

V (P θ ) = arg inf α∈Θ U α (P θ ) = arg inf α∈Θ φ(α, θ) = θ for all θ ∈ Θ.
In the following proposition, we suppose that the function m(θ, α) admits partial derivatives of order 1 and 2 with respect to θ and α and also we suppose that conditions permitting to derivate m(θ, α) under the integral sign hold.

The following result states that, unlike θ n (α), an estimator α n is generally not robust. Indeed, it has the same robustness properties as the maximum likelihood estimator, since it has its influence function which in most cases is unbounded. Whatever the divergence, the estimators α n have the same influence function.

Proposition 3

The influence function of the functional V corresponding to an estimator α n is given by

IF(x; V, P θ 0 ) = I -1 θ 0 ṗθ 0 (x) p θ 0 (x) . ( 18
)

Robustness of the estimators for scale models and location models

In this subsection, examining the expressions of the influence functions, we give conditions for attaining the B-robustness of the dual φ-divergence estimators θ n (α), as well as of the corresponding divergence estimators. The case of interest in our B-robustness study is α = θ 0 since, as observed above, the choice α = θ 0 generally leads to unbounded influence functions. For the Cressie-Read family of divergences (2) it holds

IF(x; T α , P θ 0 ) = [ m ′′ (θ 0 , α)dP θ 0 ] -1 p α p θ 0 γ ṗθ 0 dλ - p α (x) p θ 0 (x) γ ṗθ 0 (x) p θ 0 (x) (19) and IF(x; U α , P θ 0 ) = -φ(α, θ 0 ) + m(θ 0 , α, x) = -φ(α, θ 0 ) + 1 γ -1    p α p θ 0 γ-1 dP α -1    - 1 γ p α (x) p θ 0 (x) γ -1 .

Scale models

For a given density p, it holds p θ (x) = 1 θ p( x θ ) and ṗθ (x) = -1 θ 2 p x θ + x θ ṗ x θ . Consider the following conditions:

(A.1) |u ṗ(u)|du < ∞. (A.2) sup x p(α -1 x) p(θ -1 0 x) < ∞. (A.3) sup x p(θ -1 0 x) p(α -1 x) < ∞. (A.4) sup x ∂ ∂θ [log p(θ -1 0 x)] p(α -1 x) p(θ -1 0 x) γ < ∞.
Proposition 4 For scale models, if the conditions (A.2) (for the case γ > 0) or (A.3) (for the case γ < 0) together with (A.1) and (A.4) are satisfied, then

θ n (α) is B-robust.
As a particular case, consider the problem of robust estimation of the parameter θ 0 = σ of the univariate normal model, when the mean m is known, intending to use an estimator θ n (σ) with σ = σ. We are interested on those divergences from the Cressie-Read family and those possible values of σ for which θ n (σ) is B-robust. We have

IF(x; T σ , P σ ) = [ m ′′ (σ, σ)dP σ ] -1 p σ p σ γ ṗσ p σ dP σ - p σ (x) p σ (x) γ ṗσ (x)
p σ (x) . It is easily seen that IF(x; T σ , P σ ) is bounded whenever the function p σ (x)

pσ(x) γ ṗσ(x) pσ(x)
is bounded. Since

p σ (x) p σ (x) γ ṗσ (x) p σ (x) = σ γ-1 σ γ x -m σ 2 -1 exp - 1 2 x -m σ 2 -
xm σ calculation shows that these choices of γ and σ assure that m ′′ (σ, σ)dP σ is finite and non zero. However, when using the modified Kullback-Leibler divergence (γ=0), none of the estimators θ n (σ) is B-robust, the function (20) being unbounded. These aspects can also be observed in Figure 1, which presents influence functions for different divergences when σ = 1 and σ = 1.9. The negative values of the influence function in a neighborhood of 0 is explained by the decrease of the variance estimate when oversampling close to the mean.

The asymptotic relative efficiency of an estimator is the ratio of the asymptotic variance of the maximum likelihood estimator to that of the estimator in question. For the scale normal model, the choice of σ close to σ assures a good efficiency of θ n (σ) and also the B-robustness property. Then, the bigger is the value of |γ|, the smaller is the gross error sensitivity of the estimator. For example, for σ = 1 and σ = 0.99, the efficiency of θ n (σ) is 0.9803 when γ = 0.5, 0.9615 when γ = 1, 0.9266 when γ = 2 and 0.8947 when γ = 3, the most B-robust estimator corresponding to γ = 3. As can be inferred from Figure 1, the curves IF 2 (x; T σ , P σ ) are ordered decreasingly with respect to |γ|. Therefore, large values of |γ| lead to small gross error sensitivities and low efficiencies, since the asymptotic variance of

θ n (σ) is [ IF 2 (x; T σ , P σ )dP σ ] -1
(see also Hampel et al. [11] for this formula).

For scale models, conditions of Proposition 4 assure that θ n (α) and the corresponding divergence estimator φ n (α, θ 0 ) are B-robust.

Location models

It holds p θ (x) = p(xθ).

Proposition 5 For location models, if the condition

sup x p(x -α) p(x -θ 0 ) γ ∂ ∂θ log p(x -θ 0 ) < ∞ (21) is satisfied, then θ n (α) is B-robust.
For the Cauchy density the maximum likelihood estimator exists, it is consistent, efficient and B-robust and all the estimators θ n (α) exist and are B-robust. Indeed, condition [START_REF] Simpson | Hellinger deviance tests: efficiency, breakdown points, and examples[END_REF] writes

sup x 2 1 + (x -θ 0 ) 2 1 + (x -α) 2 γ x -θ 0 1 + (x -θ 0 ) 2 < ∞
and is fulfilled for any γ and any α. Also, the integral m ′′ (θ 0 , α)dP θ 0 exists and is different to zero for any γ and any α. This is quite natural since sampling of the Cauchy law makes equivalent outliers and large sample points due to heavy tails. However it is known that the likelihood equation for Cauchy distribution has multiple roots. The number of solutions behaves asymptotically as two times a Poisson(1/π) variable plus 1 (see van der Vaart [START_REF] Van Der Vaart | Asymptotic Statistics[END_REF] p. 74). The possible selection rule for the estimate is to check the nearly common estimates for different α and φ-divergences. Figure 2 presents influence functions IF(x; T α , P θ 0 ), when γ ∈ {-1, 0, 1, 2, 3, }, θ 0 = 0.5 and α = 0.8. For these choices of θ 0 and α, the efficiency of θ n (α) is 0.9775 when γ = 1, 0.9208 when γ = 2, 0.8508 when γ = 3. Here, when γ increases, the decrease of the efficiency is worsened by a loss in B-robustness. In this respect, the maximum likelihood estimator appears as a good choice in terms of robustness and efficiency.

In the case of the logistic location model, a simple calculation shows that the condition ( 21) is fulfilled for any γ and any α. Also, the integral m ′′ (θ 0 , α)dP θ 0 exists and is different from zero for any γ and any α. These conditions entail the fact that all the estimators θ n (α) are B-robust. Figure 4 presents influence functions IF(x; T α , P θ 0 ), when γ ∈ {-1, 0, 0.5, 1, 2, 3, }, θ 0 = 1 and α = 1.5.

As in the case of the Cauchy model, when γ increases, the decrease of the efficiency is worsened by the increase of the gross error sensitivity, such that the maximum likelihood estimator appears again as a good choice in terms of robustness and efficiency.

On the other hand, for the mean of the normal law, none of the estimators θ n (α) is B-robust, their influence functions being always unbounded. In the case of the Cauchy model, as well as in the case of the logistic model, IF(x; U α , P θ 0 ) is bounded for any γ and any α. In Figure 3, respectively in Figure 5, we present such influence functions for different choices of γ. Thus, for these two location models, all the estimators φ n (α, θ 0 ) are B-robust. 4 Robust tests based on divergence estimators

Asymptotic results for contaminated observations

This subsection presents some asymptotic results that are necessary in order to analyze the robustness of some tests based on divergence estimators. These asymptotic results are obtained for contaminated observations, namely X 1 , . . . , X n are i.i.d. with

P P n,ε,x := 1 - ε √ n P θn + ε √ n δ x (22)
where θ n = θ 0 + ∆ √ n , ∆ being an arbitrary vector from R d .

For α fixed consider the following conditions:

(C.1) The function θ → m(θ, α) is C 3 for all x and all partial derivatives of order 3 of θ → m(θ, α) are dominated by some P θn -integrable function x → H(x) with the property H 2 dP θn is finite, for any n and any ∆.

(C.2) m(θ 0 , α)dP θn and m 2 (θ 0 , α)dP θn are finite, for any n and any ∆.

(C.3) m ′ (θ 0 , α)dP θn and m ′ (θ 0 , α)m ′ (θ 0 , α) t dP θn exist, for any n and any ∆.

(C.4) m ′′ (θ 0 , α)dP θn and m ′′ (θ 0 , α) 2 dP θn exist, for any n and any ∆.

The estimators θ n (α) have good properties with respect to contamination in terms of consistency. Also, φ n (α, θ 0 ) enjoys normal convergence under [START_REF] Tamura | Minimum Hellinger distance estimation for multivariate location and covariance[END_REF].

Proposition 7 If α = θ 0 and the conditions (C.1) -(C.4) are satisfied, then √ n( φ n (α, θ 0 ) -U α (P P n,ε,x )) [ IF 2 (y; U α , P P n,ε,x )dP P n,ε,x (y)] 1/2
converges in distribution to a normal standard variable.

Robust tests based on divergence estimators

In this subsection we propose tests based on dual φ-divergence estimators and study their robustness properties. We mention that the use of the dual form of a divergence to derive robust tests was discussed in a different context by Broniatowski and Leorato [START_REF] Broniatowski | An estimation method for Neyman chisquare divergence with application to test of hypothesis[END_REF] in the case of the Neyman χ 2 divergence.

For testing the hypothesis θ = θ 0 against the alternative θ = θ 0 , consider the test of level α 0 defined by the test statistic φ n := φ n (α, θ 0 ) with α = θ 0 and by the critical region

C := √ n( φ n -φ(α, θ 0 )) [ IF 2 (y; U α , P θ 0 )dP θ 0 (y)] 1/2 ≥ q 1-α 0 2 where q 1-α 0 2 is the (1 -α 0
2 )-quantile of the standard normal distribution.

Due to the asymptotic normality of φ n , for n large, the level writes as

α 0 ≃ P θ 0 √ n( φ n -φ(α, θ 0 )) [ IF 2 (y; U α , P θ 0 )dP θ 0 (y)] 1/2 ≥ q 1-α 0 2 (23) = P θ 0 (| φ n -φ(α, θ 0 )| ≥ ( √ n) -1 [ IF 2 (y; U α , P θ 0 )dP θ 0 (y)] 1/2 q 1-α 0 2 ) (24) = 2P θ 0 ( φ n ≥ k n (α 0 )) ( 25
)
where

k n (α 0 ) = ( √ n) -1 [ IF 2 (y; U α , P θ 0 )dP θ 0 (y)] 1/2 q 1-α 0 2 + φ(α, θ 0 ).
We work with the form (25) of the level and consequently of the probability to reject the null hypothesis, this being easier to handle in the proofs of the results that follows.

Consider the sequence of contiguous alternatives θ n = θ 0 + ∆n -1/2 , where ∆ is any vector from R d . When θ n tends to θ 0 , the contamination must converge to 0 at the same rate, to avoid the overlapping between the neighborhood of the hypothesis and that of the alternative (see Hampel et al. [START_REF] Hampel | Robust Statistics: the Approach Based on Influence Functions[END_REF], p.198 and Heritier and Ronchetti [START_REF] Heritier | Robust bounded-influence tests in general parametric models[END_REF]). Therefore we consider the contaminated distributions

P L n,ε,x = 1 - ε √ n P θ 0 + ε √ n δ x (26)
for the level and

P P n,ε,x = 1 - ε √ n P θn + ε √ n δ x (27)
for the power.

The asymptotic level (the asymptotic power) under (26) (under (27)) will be evaluated now.

Let β 0 = lim n→∞ 2P θn ( φ n ≥ k n (α 0 )) be the asymptotic power of the test under the family of alternatives P θn . The test is robust with respect to the power if the limit of the powers under the contaminated alternatives stays in a bounded neighborhood of β 0 , so that the role of the contamination is somehow controlled. Also, the test is robust with respect to the level if the limit of the level under the contaminated null distributions stays in a bounded neighborhood of α 0 .

Let P n,ε,x = 2P P n,ε,x ( φ n ≥ k n (α 0 )). In the same vein as in Dell'Aquilla and Ronchetti [START_REF] Dell'aquilla | Robust tests of predictive accuracy[END_REF] it holds: Proposition 8 If the conditions (C.1) -(C.4) are fulfilled, then the asymptotic power of the test under P P n,ε,x is given by

lim n→∞ P n,ε,x = 2 -2Φ Φ -1 1 - α 0 2 -∆ c [ IF 2 (y; U α , P θ 0 )dP θ 0 (y)] 1/2 - -ε IF(x; U α , P θ 0 ) [ IF 2 (y; U α , P θ 0 )dP θ 0 (y)] 1/2 ( 28
)
where c = m(θ 0 , α, y) ṗθ 0 (y) p θ 0 (y) dP θ 0 (y) and Φ is the cumulative distribution function of the standard normal.

A Taylor expansion with respect to ε yields

lim n→∞ P n,ε,x = 2 -2Φ Φ -1 1 - α 0 2 -∆ c [ IF 2 (y; U α , P θ 0 )dP θ 0 (y)] 1/2 + +2εf Φ -1 1- α 0 2 -∆ c [ IF 2 (y; U α , P θ 0 )dP θ 0 (y)] 1/2 IF(x; U α , P θ 0 ) [ IF 2 (y; U α , P θ 0 )dP θ 0 (y)] 1/2 +o(ε) = β 0 +2εf Φ -1 1- α 0 2 -∆ c [ IF 2 (y; U α , P θ 0 )dP θ 0 (y)] 1/2 IF(x; U α , P θ 0 ) [ IF 2 (y; U α , P θ 0 )dP θ 0 (y)] 1/2 + +o(ε)
where β 0 is the asymptotic power for the non contaminated model and f is the density of the standard normal distribution.

In order to limit the bias in the power of the test it is sufficient to bound the influence function IF(x; U α , P θ 0 ). Bounding the influence function is therefore enough to maintain the power in a pre-specified band around β 0 .

Let L n,ε,x = 2P L n,ε,x ( φ n ≥ k n (α 0 )
). Putting ∆ = 0 in (28) yields:

Proposition 9 If the conditions (C.1) -(C.4) are fulfilled, then the asymptotic level of the test under P L n,ε,x is given by

lim n→∞ L n,ε,x = 2 -2Φ Φ -1 1 - α 0 2 -ε IF(x; U α , P θ 0 ) [ IF 2 (y; U α , P θ 0 )dP θ 0 (y)] 1/2 = α 0 + εf Φ -1 1 - α 0 2 IF(x; U α , P θ 0 ) [ IF 2 (y; U α , P θ 0 )dP θ 0 (y)] 1/2 + o(ε).
Hence, when IF(x; U α , P θ 0 ) is bounded, L n,ε,x remains between pre-specified bounds of α 0 .

As the Proposition 8 and Proposition 9 show, both the asymptotic power of the test under P P n,ε,x and the asymptotic level of the test under P L n,ε,x are controlled by the influence function of the test statistic. Hence, the robustness of the test statistic φ n , as discussed in the previous section, assures the stability of the test under small arbitrary departures from the null hypothesis, as well as a good power under small arbitrary departures from specified alternatives. Figures 3 and5 provide some specific values of γ and α inducing robust tests for θ 0 corresponding to those models.

Simulation results

Simulation were run in order to examine empirically the performances of the robust dual φ-divergence estimators and tests. The considered parametric model was the scale normal model with known mean. We worked with data generated from the model, as well as with contaminated data.

To make some comparisons, beside dual φ-divergence estimators, we considered minimum density power divergence estimators of Basu et al. [START_REF] Basu | Robust and efficient estimation by minimising a density power divergence[END_REF] (MD-PDE's) and the maximum likelihood estimator (MLE). Recall that a MDPDE of a parameter θ is obtained as solution of the equation ṗθ (z)p β θ (z)dz -

1 n n i=1 ṗθ (X i )p β-1 θ (X i ) = 0 (29)
with respect to θ, where β > 0 and X 1 , . . . , X n is a sample from P θ . In the case of the scale normal model N (m, σ), equation ( 29) writes as

1 σ β+2 ( √ 2π) β+1 e -1 2 ( z-m σ ) 2 β+1 -1 + z -m σ 2 dz - - 1 n n i=1 1 σ β+1 ( √ 2π) β e -1 2 ( X i -m σ ) 2 β -1 + X i -m σ 2 = 0
and the MDPDE of the parameter σ is robust for any β > 0.

In a first Monte Carlo experiment the data were generated from the scale normal model N (0, 1) with mean m = 0 known, σ = 1 being the parameter of interest. We considered different choices for the tuning parameter α and for the Cressie-Read divergence to compute DφE's, and different choices for the tuning parameter β in order to compute MDPDE's. For each set of configurations considered, 5000 samples of size n = 100 were generated from the model, and for each sample DφE's, MDPDE's and MLE were obtained.

In Table 1 we present the results of the simulations, showing simulation based estimates of the bias and MSE given by

Bias = 1 n s ns i=1 ( σ i -σ), MSE = 1 n s ns i=1 ( σ i -σ) 2 ,
where n s denotes the number of samples (5000 in our case) and σ i denotes an estimate of σ for the ith sample. Examination of the table shows that DφE's give as good results as MDPDE's or MLE.

In a second Monte Carlo experiment, we first generated samples with 100 observations, namely 98 coming from N (0, 1) and 2 outliers x = 10 and then we generated samples with 100 observations, namely 96 from N (0, 1) and 4 outliers x = 10. The tuning parameters were the same as in the non contaminated case and also n s = 5000. The simulation results are given in Table 2. As can be seen, the results for DφE's and MDPDE's are comparable, they being better than the results for MLE in both cases.

In Figure 7 are represented relative errors of the robust tests applied to the scale normal model N (0, 1), for samples with n = 100 data, namely 98 data generated from N (0, 1) and 2 outliers x = 10. We considered σ = 1.9 and γ ∈ {-2, -1.5}. Again, the approximation of the level of the test is good for all the considered divergences.

In Figure 8 we present relative errors of the robust tests applied to the scale normal model N (0, 1), for samples with n = 100 data, namely 96 data generated from N (0, 1) and 4 outliers x = 10. We considered σ = 1.9 and γ ∈ {-2, -1.5}.

Observe that the tests give good results for values of γ close to zero when the data are not contaminated, respectively for large negative values of γ when the data are contaminated.

Thus, the numerical results show that dual φ-divergence estimates and corresponding tests are stable in the presence of some outliers in the sample.

Table 1.

Simulation results for

An adaptive choice of the tuning parameter

At the present stage we only present some heuristic and defer the formal treatment of this proposal, which lays beyond the scope of the present work.

According to the model and the parameter to be estimated, the choice of γ should be considered with respect to the expression ( 11) which has to be bounded. We refer to the examples given in subsection 3.2 for some scale and location model.

Given a set of observations X 1 , . . . , X n an adaptive choice for α would aim at reducing the estimated maximal bias caused by an extraneous data. Define θ n (α, γ) the DφE of θ 0 on the entire set of observation. For 1 ≤ i ≤ n, let θ i n-1 (α, γ) be the DφE of θ 0 built on the leave one out data set X 1 , X 2 , . . . , X i-1 , X i+1 , . . . , X n . Define

B n (α, γ) := max i | θ n (α, γ) -θ i n-1 (α, γ)|
which measures the maximal bias caused by a single outlier and

α * (γ) := arg inf α B n (α, γ).
7 Proofs

Proof of Proposition 1

For fixed α, θ n (α) are M-estimators. In accordance with the theory regarding the M-estimators (see for example van der Vaart [START_REF] Van Der Vaart | Asymptotic Statistics[END_REF]), the so called ψ-function

corresponding to θ n (α) is ψ α (x, θ) = m ′ (θ, α, x)
and the influence function of T α is

IF(x; T α , P θ 0 ) = [M(ψ α , P θ 0 )] -1 ψ α (x, T α (P θ 0 )) (30)
where

M (ψ α , P θ 0 ) = - ∂ ∂θ [ψ α (y, θ)] θ 0 dP θ 0 (y) = -m ′′ (θ 0 , α, y)dP θ 0 (y).
Using the Fisher consistency of the functional T α ,

ψ α (x, T α (P θ 0 )) = ψ α (x, θ 0 ) = -ϕ ′′ p α p θ 0 p α p 2 θ 0 ṗθ 0 dP α + ϕ ′′ p α p θ 0 (x) p 2 α (x) p 3 θ 0 (x) ṗθ 0 (x)
which substituted in (30) leads to the announced result.

Proof of Proposition 2

Let ε > 0 and P θ 0 εx = (1ε)P θ 0 + εδ x be the contaminated model. Then

U α P θ 0 εx = m T α P θ 0 εx , α, y d P θ 0 εx (y)
= (1ε) m T α P θ 0 εx , α, y dP θ 0 (y) + εm T α P θ 0 εx , α, x and derivation yields

IF(x; U α , P θ 0 ) = ∂ ∂ε U α P θ 0 εx ε=0 = = -m(θ 0 , α, y)dP θ 0 (y) + IF(x; T α , P θ 0 ) t m ′ (θ 0 , α, y)dP θ 0 (y) + m(θ 0 , α, x) = -φ(α, θ 0 ) + m(θ 0 , α, x).

Proof of Proposition 3

For notational clearness, define T : Θ × M → Θ, T (α, P ) := T α (P ).

For each α ∈ Θ, the definition of T (α, P ) leads to m ′ (T (α, P ), α, y) dP (y) = 0.

By the very definition of V (P ) and T (V (P ), P ), they both obey

    
m ′ (T (V (P ), P ), V (P ), y) dP (y) = 0 m ′ (T (V (P ), P ), V (P ), y) dP (y) = 0 n (T (V (P ), P ), V (P ), y) dP (y) = 0

.

For the contaminated model n T V P θ 0 εx , P θ 0 εx , V P θ 0 εx , y d P θ 0 εx (y) = 0 and so

(1ε) n T V P θ 0 εx , P θ 0 εx , V P θ 0 εx , y dP θ 0 (y)+ +εn T V P θ 0 εx , P θ 0 εx , V P θ 0 εx , x = 0.

Now derivation yields

n(θ 0 , θ 0 , y)dP θ 0 (y) + ∂ ∂θ n (T (θ 0 , P θ 0 ), θ 0 , y) dP θ 0 (y) ∂ ∂α T (θ 0 , P θ 0 )IF(x; V, P θ 0 )+ +IF(x; T θ 0 , P θ 0 )} + ∂ ∂α n(θ 0 , θ 0 , y)dP θ 0 (y)IF(x; V, P θ 0 ) + n(θ 0 , θ 0 , x) = 0. ( As a consequence ∂ ∂α n(θ 0 , θ 0 , y)dP θ 0 (y) = ϕ ′′ (1) ṗθ 0 ṗt

θ 0 p 2 θ 0 dP θ 0 = ϕ ′′ (1)I θ 0 . (34)
Also

n(θ 0 , θ 0 , x) = -ϕ ′′ (1) ṗθ 0 (x) p θ 0 (x) . ( 35
)
Using the Fisher consistency of the functional T α and substituting (33), ( 34) and ( 35) in (32) it holds

IF(x; V, P θ 0 ) = I -1 θ 0 ṗθ 0 (x) p θ 0 (x) (36)
and this completes the proof.

Proof of Proposition 4

By replacing ṗθ ,

p α (x) p θ 0 (x) γ ṗθ 0 (x) = - θ γ-2 0 α γ p(α -1 x) γ p(θ -1 0 x) γ-1 + x θ 0 p(α -1 x) p(θ -1 0 x) γ ṗ(θ -1 0 x) = - θ γ-2 0 α γ p(α -1 x) γ p(θ -1 0 x) γ-1 -θ 0 p(α -1 x) γ p(θ -1 0 x) γ-1 ∂ ∂θ log p(θ -1 0 x) (37)
and similarly

p α (x) p θ 0 (x) γ ṗθ 0 (x) p θ 0 (x) = - θ γ-1 0 α γ p(α -1 x) p(θ -1 0 x) γ + x θ 0 p(α -1 x) p(θ -1 0 x) γ ṗ(θ -1 0 x) p(θ -1 0 x) = - θ γ-1 0 α γ p(α -1 x) p(θ -1 0 x) γ -θ 0 p(α -1 x) p(θ -1 0 x) γ ∂ ∂θ log p(θ -1 0 x) .
The condition (A.1) together with one of the conditions (A.2) or (A.3) (depending on the choice of γ) entails that the function (37) is integrable. On the other hand (A.2) or (A.3) together with (A.4) assure that the function in the above display is bounded. Then IF(x; T α , P θ 0 ) as it is expressed by ( 19) is a bounded function.

Proof of Proposition 5

It holds p α (x)

p θ 0 (x) γ ṗθ 0 (x) = p(x -α) γ p(x -θ 0 ) γ-1 ∂ ∂θ log p(x -θ 0 )
and

p α (x) p θ 0 (x) γ ṗθ 0 (x) p θ 0 (x) = p(x -α) p(x -θ 0 ) γ ∂ ∂θ log p(x -θ 0 ).
Then the condition [START_REF] Simpson | Hellinger deviance tests: efficiency, breakdown points, and examples[END_REF] allows to conclude that IF(x; T α , P θ 0 ), as it is expressed by [START_REF] Rüschendorf | On the minimum discrimination information theorem[END_REF], is bounded.

Proof of Proposition 6

First prove that θ n (α) -T α (P P n,ε,x ) = o P (1).

It holds m ′ (θ 0 , α)dP θ 0 = 0 and

m ′′ (θ 0 , α)dP θ 0 = -ϕ ′′ p α p θ 0 p 2 α p 3 θ 0 ṗθ 0 ṗt θ 0 dλ = -S. (38)
The matrix S is symmetric and positive since ϕ ′′ is positive by the convexity of ϕ. Using (C.3) in connection with the Lindeberg-Feller theorem for triangular arrays we have √ n m ′ (θ 0 , α)dP n = O P (1). Using (C.4) in connection with the Lindeberg-Feller theorem for triangular arrays yields m ′′ (θ 0 , α)dP n +S = o P (1). Now, for any θ = θ 0 +un -1/3 with u ≤ 1, a Taylor expansion of m(θ, α)dP n around θ 0 under (C.1) yields n m(θ, α)dP nn m(θ 0 , α)dP n = = n 2/3 u t m ′ (θ 0 , α)dP n + 2 -1 n 1/3 u t m ′′ (θ 0 , α)dP n u + O P [START_REF] Basu | Robust and efficient estimation by minimising a density power divergence[END_REF] uniformly on u with u ≤ 1. Hence

n m(θ, α)dP n -n m(θ 0 , α)dP n = O P (n 1/6 ) -2 -1 n 1/3 u t Su + O P (1)
uniformly on u with u ≤ 1. Hence uniformly on u with u = 1,

n m(θ, α)dP n -n m(θ 0 , α)dP n ≤ O P (n 1/6 ) -2 -1 cn 1/3 + O P (1) (39)
where c is the smallest eigenvalue of the matrix S. Note that c is positive since S is positive definite. In view of (39), by the continuity of θ → m(θ, α)dP n , it holds that as n → ∞, with probability one, θ → m(θ, α)dP n attains its maximum at some point θ n (α) in the interior of the ball {θ : θ-θ 0 ≤ n -1/3 }, and therefore θ n (α)θ 0 = o P (1).

On the other hand,

T α (P P n,ε,x ) = θ 0 + ε √ n IF(x; T α , P θ 0 ) + ∆ √ n 1 + ρ ε √ n , ∆ √ n √ ε 2 + ∆ 2 n
where 1 := (1, . . . , 1) t above coincides with ∂ ∂ ∆

[T α (P θ 0 + ∆ )] ∆=0 by the Fisher consistency of the functional T α and the function ρ satisfies lim n→∞ ρ ε √ n , ∆ √ n = 0. Then T α (P P n,ε,x )θ 0 converges to zero in probability as n → ∞. Combining this with (40) we obtain that θ n (α)-T α (P P n,ε,x ) converges to zero in probability.

In the following, we prove that √ n( θ n (α) -T α (P P n,ε,x )) = O P (1).

By Taylor expansion, exists θ n inside the segment that links T α (P P n,ε,x ) and θ n (α) such that

0 = m ′ ( θ n (α), α)dP n = m ′ (T α (P P n,ε,x ), α)dP n + m ′′ (T α (P P n,ε,x ), α)dP n ( θ n (α) -T α (P P n,ε,x )) + + 1 2 ( θ n (α) -T α (P P n,ε,x )) t m ′′′ (θ n , α)dP n ( θ n (α) -T α (P P n,ε,x )). (41)
By condition (C.1), using the sup-norm

m ′′′ (θ n , α)dP n = 1 n n k=1 m ′′′ (θ n , α)(X k ) ≤ 1 n n k=1 H(X k ).
Applying the Lindeberg-Feller theorem for triangular arrays yields m ′′′ (θ n , α)dP n = O P (1). Then the last term in (41) writes o P (1)( θ n (α) -T α (P P n,ε,x )).

Under (C.1) and (C.4), by applying a Taylor expansion and repeatedly the Lindeberg-Feller theorem for triangular arrays,

m ′′ (T α (P P n,ε,x ), α)dP n = 1 n n k=1 m ′′ (T α (P P n,ε,x ), α, X k ) = 1 n n k=1 m ′′ (θ 0 , α, X k ) + ε n √ n n k=1 m ′′′ (θ 0 , α, X k )IF(x; T α , P θ 0 ) + + ∆ n √ n n k=1 m ′′′ (θ 0 , α, X k )1 + ρ ε √ n , ∆ √ n √ ε 2 + ∆ 2 n = P θ 0 m ′′ (θ 0 , α) + o P (1)
where 1 := (1, . . . , 1) t above coincides with ∂ ∂ ∆

[T α (P θ 0 + ∆ )] ∆=0 by the Fisher consistency of the functional T α and the function ρ satisfies lim n→∞ ρ

ε √ n , ∆ √ n = 0. Therefore (41) becomes -m ′ (T α (P P n,ε,x ), α)dP n = ( m ′′ (θ 0 , α)dP θ 0 + o P (1))( θ n (α) -T α (P P n,ε,x )). (42)
We prove that √ n m ′ (T α (P P n,ε,x ), α)dP n is O P [START_REF] Basu | Robust and efficient estimation by minimising a density power divergence[END_REF]. By Taylor expansion,

m ′ (T α (P P n,ε,x ), α)dP n = 1 n n k=1 m ′ (θ 0 , α, X k )+ ε n √ n n k=1 m ′′ (θ 0 , α, X k )IF(x; T α , P θ 0 ) + + ∆ n √ n n k=1 m ′′ (θ 0 , α, X k )1 + ρ ε √ n , ∆ √ n √ ε 2 + ∆ 2 n
and therefore

√ n m ′ (T α (P P n,ε,x ), α)dP n = 1 √ n n k=1 m ′ (θ 0 , α, X k ) + + ε n n k=1 m ′′ (θ 0 , α, X k )IF(x; T α , P θ 0 ) + ∆ n n k=1 m ′′ (θ 0 , α, X k )1 + ρ ε √ n , ∆ √ n ε 2 + ∆ 2 n .
Under (C.3) and (C.4), by applying the Lindeberg-Feller theorem for triangular arrays it holds √ n m ′ (T α (P P n,ε,x ), α)dP n = O P (1). Then from (42) √ n( θ n (α) -T α (P P n,ε,x )) = O P (1).

Proof of Proposition 7

By Taylor expansion, there exists θ n inside the segment that links T α (P P n,ε,x ) and θ n (α) such that {m(T α (P P n,ε,x ), α, X k )-P P n,ε,x m(T α (P P n,ε,x ), α)}.

φ n (α, θ 0 ) = m( θ n (α), α)dP n = m(T α (P P n,ε,x ), α)dP n + m ′ (T α (P P n,ε,x ), α) t dP n ( θ n (α) -T α (P P n,ε,x )) + + 1 2 ( θ n (α) -T α (P P n,ε,x )) t m ′′ (T α (P P n,ε,x ), α)dP n ( θ n (α) -T α (P P n,ε,x )) + + 1 3! 1≤i,j,k≤d ( θ n (α)-T α (P P n,ε,x )) i ( θ n (α)-T α (P P n,ε,x )) j ( θ n (α)-T α (P P n,ε,x )) k ∂ 3 m(θ n , α) ∂θ i ∂θ j ∂θ k dP n . Then √ n( φ n (α, θ 0 ) -U α (P P n,ε,x )) [ IF 2 (y; U α , P P n,ε,x )dP P n,ε,x (y)]
Apply the Lindeberg-Feller theorem for the triangular array Z n,k := m(T α (P P n,ε,x ), α, X k ) -P P n,ε,x m(T α (P P n,ε,x ), α).

For this, compute first Var(Z n,k ). Observe that Var(Z n,k ) = m 2 (T α (P P n,ε,x ), α, y)dP P n,ε,x (y)m(T α (P P n,ε,x ), α, y)dP P n,ε,x (y) where P P n,ε,x ty = (1t)P P n,ε,x + tδ y . Also U α ( P P n,ε,x ty ) = (1t) m(T α ( P P n,ε,x ty ), α, z)dP P n,ε,x (z) + tm(T α ( P P n,ε,x ty ), α, y) whence IF(y; U α , P P n,ε,x ) =m(T α (P P n,ε,x ), α, z)dP P n,ε,x (z)+ + m ′ (T α (P P n,ε,x ), α, z)dP P n,ε,x (z)IF(y; T α , P P n,ε,x ) + m(T α (P P n,ε,x ), α, y).

By the definition of T α (P P n,ε,x ), it holds m ′ (T α (P P n,ε,x ), α, z)dP P n,ε,x (z) = 0 and hence (44) holds. Here we observe that IF(y; T α , P P n,ε,x ) is finite for any y, for any n and any ∆, since IF(y; T α , P θ 0 ) is. Thus, by Lindeberg-Feller theorem for triangular arrays, the first term in the expansion (43) converges in distribution to a variable N (0, 1).

We have m ′ (T α (P P n,ε,x ), α)dP n = o P (1) since √ n m ′ (T α (P P n,ε,x ), α)dP n = O P (1) (see the proof of Proposition 6). Also, it holds m ′′ (T α (P P n,ε,x ), α)dP n = O P (1) and ∂ 3 m(θn,α) ∂θ i ∂θ j ∂θ k dP n = O P (1).

Consequently, using Proposition 6 we obtain the announced result.

Proof of Proposition 8

The level α 0 is given by 2 (y; U α , P θ 0 )dP θ 0 (y)] 1/2 ≥ √ n(k n (α 0 ) -U α (P θ 0 )) [ IF 2 (y; U α , P θ 0 )dP θ 0 (y)] 1/2 .

α 0 = 2P θ 0 ( φ n ≥ k n (α 0 )) = 2P θ 0 √ n( φ n -U α (P θ 0 )) [ IF
Using the asymptotic normality of φ n in the case of uncontaminated observa-tions (see Broniatowski and Keziou [6]), √ n(k n (α 0 ) -U α (P θ 0 )) [ IF 2 (y; U α , P θ 0 )dP θ 0 (y)]

Fig. 1 .

 1 Fig. 1. Influence functions IF(x; T σ , P σ ) for normal scale model, when m = 0, the true scale parameter is σ = 1 and σ = 1.9.

Fig. 2 .

 2 Fig. 2. Influence functions IF(x; T α , P θ 0 ) for the Cauchy location model, when the true location parameter is θ 0 = 0.5 and α = 0.8.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Influence functions IF(x; U α , P θ 0 ) for the Cauchy location model, when the true location parameter is θ 0 = 0.5 and α = 0.8.

Fig. 5 .

 5 Fig. 5. Influence functions IF(x; U α , P θ 0 ) for the logistic location model, when the true location parameter is θ 0 = 1 and α = 1.5.

Proposition 6

 6 If the conditions (C.1), (C.3) and (C.4) are satisfied, then √ n( θ n (α) -T α (P P n,ε,x )) = O P (1).

  DφE, MDPDE and MLE of the parameter σ = 1 when the data are generated from the model N (0, 1).

Fig. 6 .Fig. 7 .Fig. 8 .

 678 Fig.[START_REF] Broniatowski | Parametric estimation and tests through divergences and duality technique[END_REF]. Relative errors of the robust tests applied to the scale normal model N (0, 1), when σ = 1.9 and 100 data are generated from model.

∂

  ∂α [m (T (α, P ), α, y)] V (P ) dP (y) = 0 .(31)Denoting n(θ, α, y) = ∂ ∂α m(θ, α, y)

Fromm

  ′ (T (V (P ), P ), V (P ), y) dP (y) = 0 ∂ ∂α [T (α, P )] V (P ) m ′ (T (V (P ), P ), V (P ), y) dP (y)+ n (T (V (P ), P ), V (P ), y)dP (y)

Table 2 .

 2 Simulation results for DφE, MDPDE and MLE of the parameter σ = 1 when 98 data are generated from the model N (0, 1) and 2 outliers x = 10 are added, respectively when 96 data are generated from the model N (0, 1) and 4 outliers x = 10 are added.

			2 outliers			4 outliers
		σ	Bias	MSE	σ	Bias	MSE
	DφE
	σ=1.5 γ = -2	1.01186 0.01186	0.00914	1.02540 0.02540 0.00946
	σ=1.5 γ = -1.5	1.00850 0.00850	0.00816	1.01911 0.01911 0.00833
	σ=1.5 γ = -1	1.00499 0.00499	0.00697	1.01210 0.01210 0.00707
	σ=1.5 γ = -0.5	1.00171 0.00171	0.00572	1.00526 0.00526 0.00583
	σ=1.5 γ = -0.1	1.09661 0.09661	0.01641	0.99766 -0.00233 0.00088
	σ=1.9 γ = -2	1.01589 0.01589	0.01059	1.03547 0.03547 0.01182
	σ=1.9 γ = -1.5	1.01236 0.01236	0.00942	1.02840 0.02840 0.01027
	σ=1.9 γ = -1	1.00785 0.00785	0.00785	1.01912 0.01912 0.00838
	σ=1.9 γ = -0.5	1.00274 0.00274	0.00598	1.00842 0.00842 0.00637
	σ=1.9 γ = -0.1	1.06708 0.06708	0.02241	1.10531 0.10531 0.02083
	MDPDE
	β = 0.1	1.01117 0.01117	0.00646	1.02676 0.02676 0.00891
	β = 0.5	1.00700 0.00700	0.00712	1.01417 0.01417 0.00743
	β = 1	1.01406 0.01406	0.00975	1.02892 0.02892 0.01062
	β = 1.5	1.01916 0.01916	0.01148	1.03876 0.03876 0.01297
	β = 2	1.02233 0.02233	0.01254	1.04448 0.04447 0.01450
	β = 2.5	1.02450 0.02450	0.01342	1.04771 0.04771 0.01556
	MLE	1.72587 0.72587	0.52852	2.22720 1.22720 1.50701

  (α) -T α (P P n,ε,x )) t m ′′ (T α (PP n,ε,x ), α)dP n ( θ n (α) -T α (P P n,ε,x )) 2[ IF 2 (y; U α , P P n,ε,x )dP P n,ε,x (y)] 1/2 + √ n 1≤i,j,k≤d ( θ n (α)-T α (P P n,ε,x )) i ( θ n (α)-T α (P P n,ε,x )) j ( θ n (α)-T α (P P n,ε,x )) k ∂ 3 m(θn,α) ∂θ i ∂θ j ∂θ k dP n 3![ IF 2 (y; U α , P P n,ε,x )dP P n,ε,x (y)] 1/2In the following we analyze each term in the above display. It holds m(T α (P P n,ε,x ), α)dP n -U α (P P n,ε,x ) =

		1/2 =	√	n( m(T α (P P n,ε,x ), α)dP n -U α (P P n,ε,x )) [ IF 2 (y; U α , P P n,ε,x )dP P n,ε,x (y)] 1/2	+
	+	+ n( θ n . ( m ′ (T α (P P n,ε,x ), α)dP n ) t √ n( θ n (α) -T α (P P n,ε,x )) [ IF 2 (y; U α , P P n,ε,x )dP P + n,ε,x (y)] 1/2 √
					(43)
			1	n
			n	k=1

  (θ 0 , α, y) t IF(x; T α , P θ 0 ) + , α, y)m ′ (θ 0 , α, y) t IF(x; T α , P θ 0 ) + Hence the conditions (C.2) and (C.3) assure that Var(Z n,k ) is finite.We now prove the equality Var(Z n,k ) = IF 2 (y; U α , P P n,ε,x )dP P n,ε,x (y).

								(44)
	By definition,
	IF(y; U α , P P n,ε,x ) =	∂ ∂t	[U α ( P P n,ε,x ty )] t=0 ,
								2
								=
	= 1 -	ε √ n		m 2 (T α (P P n,ε,x ), α, y)dP θn (y) +	ε √ n	m 2 (T α (P P n,ε,x ), α, x) -
	-1 -	ε √ n	m(T α (P P n,ε,x ), α, y)dP θn (y) +	ε √ n	m(T α (P P n,ε,x ), α, x)
	By Taylor expansions
	m(T α (P P n,ε,x ), α, y) = m(θ 0 , α, y) + m ′ + ε √ n ∆ √ n m ′ (θ 0 , α, y) t 1 + ρ	ε √ n	,	∆ √ n	√	ε 2 + ∆ 2 n
	m 2 (T α (P P n,ε,x ), α, y) = m 2 (θ 0 , α, y) + 2 m(θ 0 2 ε √ n ∆ √ n m(θ 0 , α, y)m ′ (θ 0 , α, y) t 1 + ρ	ε √ n	,	∆ √ n	√	ε 2 + ∆ 2 n	.

2

.

  1/2 = Φ -1 1 -Now we are interested in the value of the asymptotic power, when the underlying distribution deviates slightly from the model. Using (45)P n,ε,x = 2P P n,ε,x ( φ n ≥ k n (α 0 )) -U α (P P n,ε,x )) [ IF 2 (y; U α , P P n,ε,x )dP P n,ε,x (y)] 1/2 ≥ √ n(k n (α 0 ) -U α (P P n,ε,x )) [ IF 2 (y; U α , P P n,ε,x )dP P n,ε,x (y)] 1/2 -U α (P P n,ε,x )) [ IF 2 (y; U α , P P n,ε,x )dP P n,ε,x (y)] 1/2 ≥ -√ n(U α (P P n,ε,x ) -U α (P θ 0 )) [ IF 2 (y; U α , P P n,ε,x )dP P n,ε,x (y)] 1/2 + + Φ -1 1 -α 0 2 [ IF 2 (y; U α , P θ 0 )dP θ 0 (y)] 1/2 [ IF 2 (y; U α , P P n,ε,x )dP P n,ε,x (y)] 1/2 +o(1) . Expand U α (P P n,ε,x ) around to U α (P θ 0 ) to obtain √ n(U α (P P n,ε,x ) -U α (P θ 0 )) = εIF(x; U α , P θ 0 ) + ∆ Using the asymptotic normality of the test statistic when the observations are i.i.d. with P P n,ε,x and taking into account that lim n→∞ [ IF 2 (y; U α , P P n,ε,x)dP P n,ε,x (y)] 1/2 = [ IF 2 (y; U α , P θ 0 )dP θ 0 (y)] 1/2 [U α (P θ 0 + ∆ )] ∆=0[ IF 2 (y; U α , P θ 0 )dP θ 0 (y)] 1/2 -ε IF(x; U α , P θ 0 ) [ IF 2 (y; U α , P θ 0 )dP θ 0 (y)] 1/2 . (46)

								α 0 2	+ o(1).
	Therefore
	k n (α 0 ) = U α (P θ 0 ) +	1 √ n	Φ -1 1 -	α 0 2	[ IF 2 (y; U α , P θ 0 )dP θ 0 (y)] 1/2 + o	1 √ n	.
								(45)
	= 2P P n,ε,x	 	√	n( φ n  
	= 2P P n,ε,x	 	√	n( φ n ∂ ∂ ∆	[U α (P θ 0 + ∆ )] ∆=0 +
							+ρ	ε √ n	,	∆ √ n	ε 2 + ∆ 2 n	.
	it holds
	lim n→∞	P n,ε,x = 2 -2Φ	  Φ -1 1 -	α 0 2	-∆	∂ ∂ ∆

γ[START_REF] Simpson | Minimum Hellinger distance estimation for the analysis of count data[END_REF] boundedness of IF(x; T σ , P σ ) holds when γ > 0 and σ < σ or when γ < 0 and σ > σ, cases in which the conditions of Proposition 4 are satisfied. A simple

A close look at the results of the simulations show the DφE performs well under the model, when no outliers are generated; indeed the best results are obtained when γ = -0.1, whatever σ = 1.5 or σ = 1.9. The performance of the estimator under the model is comparable to that of some MDPDE's in terms of empirical MSE ( MSE): indeed the MSE for DφE with γ = -0.1 parallels MDPDE's for small β. It is also slightly shorter than the one obtained through the MLE. Under contamination, the DφE with γ = -0.5 yields clearly the most robust estimate and the empirical MSE is very small, indicating a strong stability of the estimate. It compares favorably with MDPE for all β, whatever σ = 1.5 or σ = 1.9. The simulation with 4 outliers at x = 10 provide a clear evidence of the properties of the DφE with γ = -0.5. Also small values of β give similar results as large negative values of γ, whatever σ, under contamination. Although γ = -0.1 is a good alternative to MLE under the model, γ = -0.5 behaves quite well in terms of bias while keeping short empirical MSE under the model or under contamination. These results are in full accordance with Figure 1; indeed the influence function is constant close to 0 for large values of x.

Thus, the DφE is shown to be an attractive alternative to both the MLE and MDPDE in these settings.

In order to test the hypothesis σ = 1 with respect to the alternative σ = 1, we considered the test statistic

(here θ 0 = σ = 1). Under the null hypothesis, this test statistic is asymptotically N (0, 1). We worked with data generated from the model N (0, 1), as well as with contaminated data. In each case, we simulated 5000 samples and we computed the actual levels

corresponding to the nominal levels α 0 = 0.01, 0.02, . . . , 0.1. We reported the corresponding relative errors

In Figure 6 we present relative errors for the robust tests applied to the scale normal model N (0, 1), when the data are generated from the model. The sample size is n = 100, the tuning parameter is σ = 1.9 and the Cressie-Read divergences correspond to γ ∈ {-1.5, -1, -0.5, -0.1}. The approximation of the level is good for all the considered divergences.

A simple calculation shows that ∂ ∂ ∆

[U α (P θ 0 + ∆ )] ∆=0 = m(θ 0 , α, y) ṗθ 0 (y) p θ 0 (y) dP θ 0 (y). Hence (28) holds.


978-1-4244-3304-9/09/$25.00 ©2009 IEEE DiaSim: A Parameterized Simulator for Pervasive Computing Applications Wilfried Jouve, Julien Bruneau, Charles Consel INRIA / LaBRI / ENSEIRB, Talence, France jouve@labri.fr, bruneau@enseirb.fr, consel@labri.fr Abstract-Pervasive computing applications involve both software concerns, like any software system, and integration concerns, for the constituent networked devices of the pervasive computing environment. This situation is problematic for testing because it requires acquiring, testing and interfacing a variety of software and hardware entities. This process can rapidly become costly and time-consuming when the target environment involves many entities.

In this demonstration, we present DiaSim, a simulator for pervasive computing applications. To cope with widely heterogeneous entities, DiaSim is parameterized with respect to a description of a target pervasive computing environment. This description is used to generate both a programming framework to develop the simulation logic and an emulation layer to execute applications. Furthermore, a simulation renderer is coupled to DiaSim to allow a simulated pervasive system to be visually monitored and debugged.

I. INTRODUCTION

Pervasive computing applications coordinate a variety of networked entities collecting context data from sensors and reacting by triggering actuators. To collect context data, sensors process stimuli that are observable changes of the environment (e.g., fire and motion). Developing a pervasive computing application requires to address a number of issues such as entity heterogeneity, physical constraints, and types of stimuli present in the target environment. Also, such an application needs to implement strategies to manage a variety of scenarios e.g., fire situations, intrusions, and crowd emergency-escape plans. Consequently, in addition to the challenges of developing any software system, a pervasive computing system needs to validate the application building blocks both individually and globally, to identify potential conflicts. In practice, the many parameters to take into account for the development of a pervasive computing application can considerably lengthen this process and a fully-equipped pervasive computing environment is still required to run and test an application. As a result, an iteration process is needed, involving the physical layout of the target environment and the application code. This process is cumbersome and hinders testing applications against a wide range of scenarios. To ease this process, we propose DiaSim, a parameterized simulator for pervasive computing applications.

II. OUR APPROACH

DiaSim relies on the DiaGen approach [START_REF] Jouve | A SIP-based programming framework for advanced telephony applications[END_REF], [START_REF] Jouve | High-level programming support for robust pervasive computing applications[END_REF] where a high-level specification of a pervasive computing environment written in the DiaSpec language is passed to a compiler to produce a customized programming framework. This programming framework provides developers with high-level abstractions to discover services and to communicate with them. In the DiaSim approach, DiaGen automates the production of simulation environments, enables transparent simulation of applications and allows actual entities to be tested together with simulated ones. In the rest of this section, we describe the key features of our approach.

Parameterized simulator. Pervasive computing systems target a variety of application areas, including home automation, building surveillance and assisted living. A simulation tool for the pervasive computing domain is required to deal with different application areas, enabling new classes of entities and stimuli to be introduced easily. To adapt to various application areas, the DiaSim simulator is parameterized with respect to a high-level specification of a pervasive computing environment written in the DiaSpec language.

Transparent simulation. Our approach makes it possible for the same code to be simulated or executed in the actual environment. DiaSim emulate the execution of an application without requiring any change in the application code. As a result, when the testing phase is completed, the application code can be uploaded as is and its logic does not require further debugging. To do that, we ensure a functional correspondence between a simulated environment and an actual one by requiring both implementations to be in conformance with the same DiaSpec specification.

Generated simulation support.

A DiaSpec specification is used to generate both an emulation layer to execute applications and a simulation programming framework to develop simulated entities. Furthermore, simulation scenarios are defined in a scenario editor which is parameterized by a DiaSpec specification. Simulated services and stimuli can then be graphically defined using a wizard. To facilitate the creation of simulated services and stimulus producers, we provide libraries of predefined behaviors.

Hybrid simulated environments. DiaSim builds simulated environments as extensions of actual environments. Because of this inheritance strategy, an application can be executed in an hybrid environment, combining simulated and actual services. This feature is particularly useful to perform unit testing of actual entities.

Simulation renderer. We present a simulation renderer that enables the developer to visually monitor and debug a pervasive computing system. The simulation renderer takes into account various features of the pervasive computing domain. Specifically, it supports visual representations for an openended set of entities and stimuli, visual support for scenario monitoring, and debugging facilities to navigate in scenarios in terms of time and space.

The DiaSim simulator architecture is shown in Figure 1. Stimulus producers emit stimuli of various types according to a predefined scenario. In place of actual sensors, simulated ones process these stimuli and produce events. The unchanged application reacts to these events by invoking actuator commands. In turn, actuators change the simulated environment, triggering stimulus producers. To demonstrate the key features of DiaSim, we will simulate various applications in the building management area. DiaSim will be used to test these applications and evaluate the feasibility of their deployment in ENSEIRB, an engineering school to which the authors are affiliated.

A. Target environment

The ENSEIRB school is a three-floor building of 13,500 m 2 , consisting of several lecture halls, labs and recreation rooms for students. ENSEIRB hosts up to 900 occupants, including students and faculty members. Many devices compose the ENSEIRB school. Sensors such as motion detectors, light and temperature sensors provide context information from all over the building to various applications. The applications trigger several types of actuators. These include alarms, lights, LCD screens, PDAs, loudspeakers and air conditioners.

B. Demonstration Steps

Specifying simulated environments. The first step of our demonstration will show how the DiaSim simulator is parameterized with respect to the DiaSpec specification of the ENSEIRB environment. We will show how to model a pervasive computing environment in DiaSpec and how the resulting DiaSpec specification parameterizes the DiaSim simulator.

Defining simulation scenarios. The second step of our demonstration will be the definition of simulation scenarios. Scenarios will be defined using a Java GUI called the scenario editor (Figure 2). From a DiaSpec specification, simulated services are either graphically defined using a wizard, or developed using the generated simulation programming framework. The demonstration will show examples of both ways to define simulated services. The scenario editor also supports the definition of stimulus producers by allowing the user to define stimulus intensities in areas of the simulated space at specific moments in time. For example, a producer of motion stimuli simulates a user moving in a school hallway at a given time. Furthermore, the definition of simulation scenarios is supported by libraries of generic simulated services and stimulus producers which can be parameterized following the target simulation scenario.

To see how the applications behave in a realistic environment, the simulated scenario defines hundreds of simulated persons moving around in the school. The simulated persons are students and faculty members. Students belong to one of four departments at ENSEIRB. Their behavior depends on the agenda of their study field. When a class they have to attend begins, they go to the class room. Furthermore, they modify the context of the simulated environment. For instance, a simulated person is detected when passing in front of a motion detector. As another example, simulated persons increase the CO 2 density in the room where they are located. Testing applications. Using the programming framework generated from a DiaSpec specification, we developed various applications which will be shown in our demonstration. In the third step, these applications will be tested against the previously defined simulation scenarios. Each application requires several simulated services. To follow the evolution of the simulated environments, DiaSim extends an existing visualization tool: the Siafu open source context simulator 1 . Siafu provides a 2D rendering and time-control functionalities.

Figure 3 shows the resulting simulation renderer. Services and stimuli defined in the scenario are displayed on top of a picture of the simulated space. The simulation renderer gives the state of services by displaying a bubble of raw text above services and/or modifying the visual representation of the service. To complement these macroscopic views, we enriched Siafu's rendering functionalities with Web Interfaces and audio streams. As shown in Figure 3, in the ENSEIRB simulation, clicking on school LCD screens runs a Web interface showing its display. Similarly, loudspeakers are rendered using audio streams. Demonstrating hybrid simulated environments. One of the key feature of DiaSim is to enable hybrid simulated environments, i.e., combining simulated and actual devices in the same simulation scenario. As illustrated in Figure 3, in the final step, we will demonstrate this feature by making actual devices interact with simulated ones. These actual devices are an IP camera and an actual PDA. The camera is used as a motion detector for the intrusion manager. Messages are received on the actual PDA when an intrusion in the simulated environment is detected. These actual devices will be added incrementally in the simulation, demonstrating the flexibility of our approach.

Actual Devices Actual Web Interfaces

C. Demonstrated Applications

The first application is a newscast manager. Its purpose is to display news and schedules on screens scattered throughout the school. This manager reads the school RSS feed in order to get the latest school news. It also receives schedule events from the ENSEIRB agenda. This agenda describes the courses of the four ENSEIRB departments. The newscast manager alternatively displays the news and the students' schedule. A key feature of the newscast manager is that it adapts the contents displayed on the screens with respect to the department affiliation and the nationality of the people surrounding the screens. Thus, for instance, the messages displayed on a screen are in English if this screen is mainly surrounded by English native speakers. Moreover, when a class begins, the newscast manager plays a message on the loudspeakers of the school using the text-to-speech technology.

The second application is an intrusion manager. It is activated when ENSEIRB is closed. The intrusion manager receives intrusion events from the motion detectors. They are installed in every corridor and room of the school. When a motion is detected, the intrusion manager turns on the school alarms. It also displays a message on the PDA of the school supervisor to warn him and to indicate where the intrusion took place.

The last application is a building automation manager. This manager controls the lights depending on the outside luminosity and the area occupancy. It also controls the school air conditioners to regulate the temperature based on temperature values produced by multiple temperature sensors. When the school is closed, all lights are turned off and the target temperature is lower (in winter) or higher (in summer) than during the day to save energy. Between 7 am and 8 pm, its behavior depends on luminosity values of the outside light sensors. If the luminosity is lower than a predefined threshold, the lights of a given area are turned on, when some people enter this area, and turned off, when the area is empty. If the outside luminosity is greater than the predefined threshold, the previous behavior only occurs for areas without windows. The manager also receives information from CO 2 sensors. When the CO 2 density is too high in a given area, the corresponding air conditioner is activated to decrease its density.

D. Technical constraints

The demonstration does not have special requirements other than a demo booth or table and an Internet connection.

Fig. 1 .

 1 Fig. 1. The DiaSim simulator architecture

Fig. 2 .

 2 Fig. 2. The scenario editor

Fig. 3 .

 3 Fig. 3. Demonstration of a hybrid simulated environment

Authorized licensed use limited to: UR Futurs. Downloaded on December 15, 2009 at 11:29 from IEEE Xplore. Restrictions apply.

http://siafusimulator.sourceforge.net/ Authorized licensed use limited to: UR Futurs. Downloaded on December 15,

at 11:29 from IEEE Xplore. Restrictions apply.


Brahim Hamadicharef

Haihong Zhang
email: hhzhang@i2r.a-star.edu.sg

Cuntai Guan
email: ctguan@i2r.a-star.edu.sg

Wang Chuanchu

Soon Kok

Keng Peng Phua
email: ksphua@i2r.a-star.edu.sg

Kai Keng Tee
email: kptee@i2r.a-star.edu.sg

Ang

Learning EEG-based Spectral-Spatial Patterns for Attention Level Measurement

In our every day life, our brain is constantly processing information and paying attention, reacting accordingly, to all sorts of sensory inputs (auditory, visual, etc.). In some cases, there is a need to accurately measure a person's level of attention to monitor a sportsman performance, to detect Attention Deficit Hyperactivity Disorder (ADHD) in children, to evaluate the effectiveness of neuro-feedback treatment, etc.

In this paper we propose a novel approach to extract, select and learn spectral-spatial patterns from electroencephalogram (EEG) recordings. Our approach improves over prior-art methods that was, typically, only concerned with power of specific EEG rhythms from few individual channels. In this new approach, spectral-spatial features from multichannel EEG are extracted by a two filtering stages: a filter-bank (FB) and common spatial patterns (CSP) filters. The most important features are selected by a Mutual Information (MI) based feature selection procedure and then classified using Fisher linear discriminant (FLD). The outcome is a measure of the attention level.

An experimental study was conducted with 5 healthy young male subjects with their EEG recorded in various attention and non-attention conditions (opened eyes, closed eyes, reading, counting, relaxing, etc.). EEGs were used to train and evaluate the model using 4x4fold cross-validation procedure. Results indicate that the new proposed approach outperforms the prior-art methods and can achieve up to 89.4% classification accuracy rate (with an average improvement of up to 16%). We demonstrate its application with a two-players attention-based racing car computer game.

I. INTRODUCTION

In our every day life, our brain is constantly processing information and paying attention, reacting accordingly, to all sorts of sensory inputs (auditory, visual, etc.). There is a need to accurately measure someone's level of attention for detection, performance, monitoring, in medicine, sports, gaming, driving, etc.

A large number of studies have reported clear association between attention or relevant mental conditions and spectral features of the EEG. For example, in [START_REF] Jung | Estimating alertness from the EEG power spectrum[END_REF] the authors applied principal component analysis to full EEG log spectrum and used sub-space features to estimate a local error rate in a sustained attention task. Many other studies focused on specific EEG rhythms. In [START_REF] Kelly | EEG Alpha power and coherence time courses in a sustained attention task[END_REF], alpha band activities and coherence time courses were examined in an attention study. In [START_REF] Huang | Multi-scale EEG brain dynamics during sustained attention tasks[END_REF], experimental results also indicated correlates between alpha power dynamics and errors by subjects in sustained attention tasks. In [START_REF] Haufler | Neuro-cognitive activity during a self-paced visuospatial task: comparative EEG profiles in marksmen and novice shooters[END_REF], neuro-cognitive activity during a selfpaced visuospatial task were compared from EEG profiles of marksmen and novice shooters. Furthermore, many studies suggest that self-regulation of theta/beta reduces Attention Deficit Hyperactivity Disorder (ADHD) symptoms [START_REF] Barry | A review of electrophysiology in attention-deficit/hyperactivity disorder: I. qualitative and quantitative electroencephalography[END_REF] [6] and the basis of neuro-feedback therapy [START_REF] Lubar | Discourse on the development of EEG diagnostics and biofeedback for attention-deficit/hyperactivity disorders[END_REF] [8] [START_REF] Friel | EEG biofeedback in the treatment of attention deficit/hyperactivity disorder[END_REF].

Rather than limiting potential useful features to only specific EEG rhythms, we propose to identify and exploit joint spectral-spatial patterns. We recognized the importance and usefulness of spatial filters from other EEG research studies, for e.g. in Brain-Computer Interface (BCI) to improve spatial resolution or for discriminating motor intentions [START_REF] Müller-Gerking | Designing optimal spatial filters for single-trial EEG classification in a movement task[END_REF].

In this paper we propose a novel approach that combines advanced signal processing and machine learning (ML) to overcome the limitations of the current prior-art methods. Specifically, we use filter banks to cover a broad range of EEG rhythms, together with common spatial pattern (CSP) filtering, that makes it possible to extract (and tailor) patterns in the EEG for discriminating different EEG conditions (i.e. attention and non-attention). Because of the large number of possible features, we use a MI-based feature selection procedure to reduce them. Notice that the Filter Bank Common Spatial Pattern (FBCSP) was first proposed in [START_REF] Ang | Filter Bank Common Spatial Pattern (FBCSP) in Brain-computer Interface[END_REF], but to the best of our knowledge, this is the first time CSP is used in attention detection.

An experimental study was conducted with 5 young healthy male subjects in various attention and non-attention conditions. Results using the novel approach indicate significant improvement in accuracy and outperforming the prior-art methods based on alpha, beta, and theta EEG rhythms [START_REF] Monastra | Electroencephalographic Biofeedback in the Treatment of Attention-Deficit/Hyperactivity Disorder[END_REF] [START_REF] Leins | Neurofeedback for Children with ADHD: A Comparison of SCP and Theta/Beta Protocols[END_REF].

One important advantage of using machine learning in our approach is that it can be learn and be tuned for subjectspecific performance. The final model can be trained from a relatively small amount of EEG data, stored, reused at later stage, retrained, etc. Such flexibility is important when, for e.g., continuously monitoring an athlete's performance over the course of training program. Models can also be trained on subject-specific EEG recording but also from a pool of EEG recording of a group of subjects.

The remainder of the paper is organized as follows. In Section II, the methods are introduced. In Section III, we present an the EEG data from an experimental study carried out to evaluate the novel approach. Results are presented in Section IV. Finally, we conclude the paper in Section V.

CSP Array Filter Banks

II. METHODS

The overall flowchart of the novel approach is shown in Figure 1. It consists of four stages: a multiple bandpass filtering stage with Filter-Bank (FB), a spatial filtering stage with common spatial pattern (CSP) algorithm, a feature selection procedure and finally a classifier. The classifier outcome is a measure of attention level.

A. Filter-Bank (FB)

The filter-bank (FB), front-end on the system, consists of 9 bandpass IIR filters (Chebyshev Type II) with equally spaced central frequency (4Hz, 8Hz, 12Hz, 16Hz, 20Hz, 24Hz, 28Hz, 32Hz, 36Hz), each with an 8Hz bandwidth, a stop-band attenuation of 40db, and a pass-band ripple of 3db.

B. Common Spatial Pattern (CSP)

In the present context, the CSP algorithm [START_REF] Müller-Gerking | Designing optimal spatial filters for single-trial EEG classification in a movement task[END_REF] explores spatial patterns of brain rhythm modulations and will discriminate between two classes of EEG patterns, of class O 1 and O 2 , to obtain an optimal projection vector w onto EEG signal x that satisfies

argmax w E[(w T x i -E[w T x]) 2 |x ∈ O 1 ] s.t. E[(w T x i -E[w T x]) 2 ] = 1 (1)
This can be reformulated as

argmax w w T Σ O1 w s.t. w T (Σ O1 + Σ O2 )w = 1 (2)
where Σ Oi denotes the covariance matrix for class O i . This optimization problem can be solved using joint diagonalization, which consists of computing a projection matrix W and a diagonal matrix D:

W Σ (1) W T = D (3) W Σ (2) W T = I -D (4)
Matrices W and D can be obtained using the following procedure [START_REF] Müller-Gerking | Designing optimal spatial filters for single-trial EEG classification in a movement task[END_REF]. First we calculate the matrix P which whitens the matrix Σ O1 + Σ O2 : P (Σ O1 + Σ O2 )P T = I. This is typically carried out using singular value decomposition (SVD) and normalization. We then calculate the whitened matrix Σ O1 = P Σ O1 P T , and calculate the eigenvector matrix Q for Σ O1 = QDQ T . Finally, we calculate the projection matrix as

W = P T Q.
W clearly satisfies Eq. 4 with the columns of W being the optimal CSP projection vectors. Typically, only the vectors with the m highest eigenvalues (for the positive class, i.e. O 1 ) and lowest m eigenvalues (for the negative class i.e. O 2 ) are used. Given one w from the selected 2m projection vectors, the CSP feature, f p , is computed as the normalized log power (variance) of the filtered signal.

f p = log var w T x 2m i=1 var w T i x (5)

C. Feature Selection

Feature selection is an important aspect in machine learning and can be defined as selecting a subset of size k, from a given a set of d features, that leads to the smallest classification errors. One commonly used method, based on statistics from the data distribution, is the Principal Component Analysis (PCA) which employs covariance, the second-order statistics of the empirical data distribution [START_REF] Chow | Estimating optimal feature subsets using efficient estimation of high-dimensional mutual information[END_REF]. In this work we choose Mutual Information (MI), a method from higher order statistics [START_REF] Cover | Elements of information theory[END_REF], also considered as a good indicator, robust to noise and data transformation, to find relevance between variables. Let us consider a variable X and its corresponding class label variable Y , then the MI between X and Y is defined as:

I(X; Y ) = H(Y ) -H(Y |X) (6)
where H(X) denotes the entropy of the feature variable X and H(Y |X) represents the conditional entropy of class label variable Y given feature variable X.

H(X) = - x∈X p(x)log 2 p(x)dx (7) H(Y |X) = - x∈X y∈Y P (y|x)log 2 P (y|x)dx (8)
We opted for the Naive Bayesian Parzen Window (NBPW) approach described in [START_REF] Ang | Rough set-based neuro-fuzzy system[END_REF]. The procedure consists of the following steps:

III. EXPERIMENTAL STUDY

An experiment was conducted to record EEG for the evaluation of the novel approach to discriminate the two conditions: attention and non-attention. The experiment consists of runs with 4 different tasks as shown in Figure 2. Task 1 was a concentration task with eyes opened, in which the subject was asked to continuously read an recent IEEE Transactions paper. Task 2 was a non-concentration task with eyes opened, in which the subject stay in a relax state and while looking at surroundings without paying particular attention to any objects. Task 3 was a concentration task but this time with eyes closed, during which the subject is asked to perform a mental arithmetic calculation (400 minus 7, minus 7, minus 7, etc.). Finally, Task 4 was a non-concentration task with with eyes closed, where a subject is in relax state and not focusing on any particular thought. The ground truth used by the classifier is shown as the attention curve in Figure 2. Five healthy young male subjects participated in the experiment. Two subjects completed 3 sessions of data collection, while the remaining 3 subjects completed only 2 sessions. Each session consists of 5 runs and each run include all 4 tasks. The user was permitted a short pause of few minutes between runs to avoid fatigue. EEG data was check for apparent EMG artefacts. Subjects were seated in a comfortable armchair, and reminded to not to have large body movements in order to avoid too much EMG artefacts. During Task 3 and Task 4 of the experiment, subjects were monitored to avoid dozing off episodes.

EEG was recorded using a 40-channel EEG amplifier (Compumedics Inc., USA). The EEG cap is standard 32 channels (See Figure 3) but we limited our recordings only 15 channels

FP1 FP2 F7 F3 Fz F4 F8 FT7 FC3 FCz FC4 FT8 T7 C3 Cz C4 T8 TP 7 CP3 CPz CP4 TP8 P7 P3 Pz P4 P8 Oz O1 O2 PO1 PO2
A2 A1 GND Fig. 3. EEG cap with highlighted in grey the electrodes used in this study (F3, Fz, F4, FC3, FCz, FC4, C3, Cz, C4, CP3, CPz, CP4, P3, Pz, P4) using the right earlobe (A2) as reference. The sampling rate was 250Hz, and a hardware 50Hz notch filter was used to reduce the main interference. All electrodes had less than 5kΩ impedance.

IV. RESULTS

The objective of this study was to examine the capability of the novel approach to differentiate attention from nonattention conditions. Thus, we considered Task 1 and Task 3 as representing attention class, and Task 1 plus Task 4 as the non-attention class. Continuous data were segmented into EEG epochs, and our method was tested to classify each epoch into attention or non-attention class. It was of interest to see the system performance over difference lengths of EEG epochs. Thus, we examined EEG epochs of 2s, 4s, 8s, or 16s long.

To avoid any bias in the evaluation results, we used a typical 4x4fold cross-validation. Briefly, each subject's EEG data was divided into n parts, then where in each fold a particular part is used for testing while the other n -1 parts were used for training. In other words the test part is not used in the training. This process is repeated n times and accuracy is averaged.

The CSP was applied to find the most discriminative spatial patterns between the attention and non-attention classes from channels filtered by the filter-bank. After the feature extraction stage, we obtain 9 × 2 = 18 candidate feature set. An MIbased feature selection procedure is then applied for further reduction (typically 2 to 4 features) in combination with the FLD classification. We used a threshold δ of 0.1 for all experiments.

The baseline for comparison is the prior-art method spectralbased method [START_REF] Monastra | Electroencephalographic Biofeedback in the Treatment of Attention-Deficit/Hyperactivity Disorder[END_REF] [START_REF] Leins | Neurofeedback for Children with ADHD: A Comparison of SCP and Theta/Beta Protocols[END_REF] with alpha (8-12Hz), beta (12-30Hz), and theta (4-8Hz) band joint to form a feature input vector to a FLD classifier.

We wanted to investigate how robust was our approach when changing the window length as typically, near real-time situation (online system with continuous EEG) requires short window length, while longer windows can be used for offline situation.

Results for classification accuracy are presented in Table I. Our method outperforms the prior-art method with higher accuracies from 60.1% to 89.4% (mean from 65.25% to 85.58%) while prior-art only achieving 50.60% to 71.9% (mean 55.50% to 69.55%). With subject #2, for e.g., the improvement is significant with more than 16% (mean from 69.55% to 85.58%). Classification accuracies are also shown to be constant across all different window lengths. We also notice that performance do not drop for short window (e.g. 2s), which is an important issue when aiming at near realtime operation.

V. CONCLUSIONS

In this paper we proposed a novel approach to concentration taking advantages of signal processing and machine learning paradigms. Spectral-Spatial features are extract using filter bank and CSP filtering from EEG data. Best features are selected using a MI-based procedure and classified using FLD. Results from an 5 subjects experiment showed that the novel approach provides significant improvement over prior-art method which uses specific EEG band powers. The resulting EEG-based attention level measurement system is robust, run in near real-time, and requires only few electrodes. The system can be implemented as a low-cost portable EEG system. The novel approach can be demonstrated with a two-players attention-based racing car computer game.

Future work will focus on EEG recording conditions with a EEG recording protocol based on variants of the stroop test [START_REF] Macleod | Half a Century of Research on the Stroop Effect: An Integrative Review[END_REF]. We will also investigate the use of Relevance Vector Machine (RVM) [START_REF] Tipping | Sparse Bayesian Learning and the Relevance Vector Machine[END_REF] as classifier.

Fig. 1 .

 1 Fig. 1. Schematic flowchart of attention monitoring

  Initial candidate set of d features F = {f 1 , f 2 , . . . , f d }, and the initial select feature set as a null set F opt = ∅, I 0 = 0; 2) For each feature f k in the candidate set: a) Form a tentative feature vector F k = F opt {f k }; b) Use F t and NBPW to predict the class label Y k ; c) Compute the MI of the predicted class label and the true label I(Y k ; Y ); 3) Select the feature, say f k , which maximizes I(Y k ; Y ); 4) If F opt = ∅ and mutual information gains less than a preset threshold δ: I(Y k ; Y ) -I 0 < δ, exit; 5) I 0 = I(Y k ); 6) Update the candidate set by F → F \ {f k } and update the select feature set by F opt → F opt {f k }; 7) If the candidate set is empty, exit; otherwise go to Step 2;

	Prev Task	Eye-opened Reading (Attention State)	Eye-opened Relax (Non-attention State)	Eye-closed Arithmetic (Attention State)	Eye-closed Relax (Non-attention State)	Next Task
	Attention Curve
		Fig. 2. Sequence of attention / non-attention tasks during one run
	1) Initialization:

ACKNOWLEDGMENT

This work is supported by the Agency for Science, Technology and Research (A*STAR), Singapore.