Nicolas Gast

A REFINED MEAN FIELD APPROXIMATION OF SYNCHRONOUS DISCRETE-TIME POPULATION MODELS

Mean field approximation is a popular method to study the behaviour of stochastic models composed of a large number of interacting objects. When the objects are asynchronous, the mean field approximation of a population model can be expressed as an ordinary differential equation. When the objects are (clock-) synchronous the mean field approximation is a discrete time dynamical system. We focus on the latter.

We study the accuracy of mean field approximation when this approximation is a discrete-time dynamical system. We extend a result that was shown for the continuous time case and we prove that expected performance indicators estimated by mean field approximation are O(1/N )-accurate. We provide simple expressions to effectively compute the asymptotic error of mean field approximation, for finite time-horizon and steady-state, and we use this computed error to propose what we call a refined mean field approximation. We show, by using a few numerical examples, that this technique improves the quality of approximation compared to the classical mean field approximation, especially for relatively small population sizes.

Introduction

Stochastic models are often used to model and analyse the performance of computer (and many other) systems. A particularly rich and popular class of models is given by stochastic population models. These have been used, for instance, to model biological systems [START_REF] Wilkinson | Stochastic Modelling for Systems Biology[END_REF], epidemic spreading [START_REF] Andersson | Stochastic Epidemic Models and Their Statistical Analysis[END_REF] or queuing networks [START_REF] Vvedenskaya | Queueing system with selection of the shortest of two queues: An asymptotic approach[END_REF]. These systems are composed of a set of homogeneous objects interacting with one another. These models have a high expressive power, but an exact analysis of any such a model is often computationally prohibitive when the number of objects of the system grows. This results in the need for approximation techniques.

A popular technique is to use mean field approximation. The idea behind mean field approximation is to replace the study of the original stochastic system by the one of a, much simpler, deterministic dynamical system. The success of mean field approximation can be explained by multiple factors : (a) it is fast -many models can be solved in closed form [START_REF] Vvedenskaya | Queueing system with selection of the shortest of two queues: An asymptotic approach[END_REF][START_REF] Mitzenmacher | The power of two choices in randomized load balancing[END_REF][START_REF] Tsitsiklis | On the power of (even a little) centralization in distributed processing[END_REF][START_REF] Minnebo | A fair comparison of pull and push strategies in large distributed networks[END_REF] or easily solved numerically [START_REF] Massoulié | Coupon replication systems[END_REF][START_REF] Gast | A mean field model of work stealing in large-scale systems[END_REF][START_REF] Van Houdt | A mean field model for a class of garbage collection algorithms in flashbased solid state drives[END_REF] -(b) it is proven to be asymptotically optimal as the number of objects in the system goes to infinity [START_REF] Thomas | Solutions of Ordinary Differential Equations as Limits of Pure Jump Markov Processes[END_REF][START_REF] Boudec | A generic mean field convergence result for systems of interacting objects[END_REF][START_REF] Benaïm | A class of mean field interaction models for computer and communication systems[END_REF][START_REF] Gast | Markov chains with discontinuous drifts have differential inclusion limits[END_REF][START_REF] Bortolussi | Continuous approximation of collective system behaviour: A tutorial[END_REF]; and (c) it is often very accurate also for systems of moderate size, composed of N ≈ 100 objects.

The mean field approximation of a given model is constructed by considering the limit of the original stochastic model as the number of objects N goes to infinity. There can be two types of limits. The first type arises when the dynamics of the objects are asynchronous. In this case the mean field approximation is given by a continuous time dynamical system (often a system of ordinary differential equations) -this is the most studied case e.g. [START_REF] Thomas | Solutions of Ordinary Differential Equations as Limits of Pure Jump Markov Processes[END_REF][START_REF] Benaïm | A class of mean field interaction models for computer and communication systems[END_REF][START_REF] Bortolussi | Continuous approximation of collective system behaviour: A tutorial[END_REF]. The second type arises when the objects are synchronous. In this case the mean field approximation is a discrete time dynamical system [START_REF] Boudec | A generic mean field convergence result for systems of interacting objects[END_REF][START_REF] Gast | A mean field approach for optimization in discrete time[END_REF][START_REF] Tinnakornsrisuphap | Limit behavior of ECN/RED gateways under a large number of TCP flows[END_REF]. We focus on the latter. Contributions. Our main contribution is an extension to (synchronous) DTMC population models of the results proposed in [START_REF] Gast | A refined mean field approximation[END_REF] for (asynchronous) CTMC population models, thus providing a new approximation technique that is significantly more accurate than classical mean field approximation, especially for relatively small systems. Our results apply to the classical model of [START_REF] Boudec | A generic mean field convergence result for systems of interacting objects[END_REF][START_REF] Gast | A mean field approach for optimization in discrete time[END_REF][START_REF] Latella | On-the-fly PCTL fast mean-field approximated model-checking for self-organising coordination[END_REF]. We prove our result for the transient and the steady state dynamics. Moreover, it retains an interesting feature of mean field approximation by being computationally non-intensive.

More precisely, if M (N ) i (t) denotes the proportion of objects in a state i at time t, then the classical result of [START_REF] Boudec | A generic mean field convergence result for systems of interacting objects[END_REF] states that, as N grows large, if the vector M (N ) (0) converges almost surely to m, for some vector m, then the vector M (N ) (t) converges almost surely to a deterministic quantity µ(t) that satisfies a recurrence equation of the form µ(t + 1) = µ(t)K(µ(t)) with µ(0) = m. We show that, for any twice differentiable function h, there exists a constant V t,h such that lim

N →∞ N (E h(M (N ) t ) -h(µ(t))) = V t,h . (1)
We provide an algorithm to compute the constant V t,h by a linear dynamical system that involves the first and second derivative of the functions m → mK(m) and h. We also show that if the function m → mK(m) has a unique fixed point µ(∞) that is globally exponentially stable, then the same result holds for the steady-state : in this case, V ∞,h = lim t→∞ V t,h exists and can be expressed as the solution of a discrete-time Lyapunov equation that involves the first and second derivative of m → mK(m) and h evaluated at the point µ(∞).

By using these results, we define a quantity h(µ(t)) + V t,h /N that we call the refined mean field approximation. As opposed to the classical mean field approximation, this approximation depends on the system size N . We illustrate our theoretical results with four different examples. While these examples all show that our refined model is clearly more accurate than the classical approximations, they illustrate different characteristics. The first two examples are cases where the dynamical system has a unique exponentially stable attractor. In these examples, refined mean field provides performance estimates that are extremely accurate (the typical error between M (N ) and µ is less than 1% for N = 10). The third example is different as it is a case when the stochastic system has two absorbing states. In this case, the refined mean field is still more accurate than the classical mean field approximation but remains far from the exact values for N = 10. It is only for larger values of N that the refined mean field provides a very accurate estimate. Finally, the fourth example is a case where the mean field approximation has a unique attractor that is not exponentially stable. We observe that in this case the refined approximation provides an accurate approximation of E[M (N ) (t)] for small values of t but fails to predict correctly what happens when t is large compared to N . In fact in this case, one cannot refine the steady-state expectation by a term in O(1/N ) because the convergence is only in O(1/ √ N ) in this case. This suggests that, when using a mean field or refined mean field approximation, one has to be careful : the approximations of a system with more than one stable equilibrium or a unique but non-exponentially stable equilibrium is likely to be inaccurate for small values of N , even when one focuses on the transient behaviour. Related work. Our results extend the recent results of [START_REF] Gast | A refined mean field approximation[END_REF]. The authors of [START_REF] Gast | A refined mean field approximation[END_REF] study the steady-state of stochastic models that have a continuous-time mean field approximation. They show that Equation ( 1) is true in this case and provide a numerical algorithm to compute the constant. Our paper has two theoretical contributions with respect to [START_REF] Gast | A refined mean field approximation[END_REF] : First we show that the results also hold for models that have a discrete-time mean field approximation, and second we show how to derive these equations for the transient and steady state regimes of such systems. This means that our results remain in the realm of discrete-time models whereas in [START_REF] Gast | A refined mean field approximation[END_REF] it is shown how some discrete-time models can be transformed into density-dependent (continuous time) population models by replacing time steps with steps that last for a random time that is exponentially distributed with mean 1/N , where N is the population size. The resulting continuous time model can then be analysed using the approximation techniques for CTMC population models discussed in [START_REF] Gast | A refined mean field approximation[END_REF].

The results of [START_REF] Gast | A refined mean field approximation[END_REF] and the one of the current paper follow from a series of recent results concerning the rate of convergence of stochastic models to their mean field approximation [START_REF] Gast | Expected values estimated via mean-field approximation are 1/n-accurate: Extended abstract[END_REF][START_REF] Ying | On the approximation error of mean-field models[END_REF][START_REF] Ying | Stein's method for mean field approximations in light and heavy traffic regimes[END_REF][START_REF] Kolokoltsov | Mean Field Games and Nonlinear Markov Processes[END_REF]. The key idea behind these works is to study the convergence of the generator of the stochastic processes to the one of its mean field approximation and to use this convergence rate to obtain a bound on Equation [START_REF] Andersson | Stochastic Epidemic Models and Their Statistical Analysis[END_REF]. For the steady-state regime, this is made possible by using Stein's methods [START_REF] Stein | Approximate computation of expectations[END_REF][START_REF] Braverman | Steins method for steady-state diffusion approximations of M/P h/n + M systems[END_REF][START_REF] Braverman | Steins method for steady-state diffusion approximations: An introduction through the erlang-a and erlang-c models[END_REF]. Note that the approach taken in the current paper is fundamentally different from the one that is usually used to obtain convergence rates, like [START_REF] Gast | Markov chains with discontinuous drifts have differential inclusion limits[END_REF][START_REF] Bortolussi | Bounds on the deviation of discrete-time markov chains from their mean-field model[END_REF][START_REF] Gast | A mean field approach for optimization in discrete time[END_REF] in which the authors focus on sample path convergence and obtain bounds on the convergence of the expected distance E[ M (N ) -µ ] between the stochastic system and its mean field approximation. When focusing on sample path convergence, the refinement of the mean field approximation would be to consider an additive term of 1/ √ N times a Gaussian noise as for example in [START_REF] Gast | Markov chains with discontinuous drifts have differential inclusion limits[END_REF] and not a 1/N term as in this paper. Outline. The rest of the paper is organised as follows. In Section 2, we introduce the model that we study. In Section 3 we provide the main results and in particular Theorem 1. In Sections 4, 5, 6 and 7, we provide a few numerical examples that demonstrate the accuracy of the refined mean field approximation and its limits. Finally, we conclude in Section 8.

Preliminaries

In this section we introduce some terminology and notation as well as some preliminary definitions, setting the context for the rest of the paper.

2.1. Notations. We let IN denote the set of natural numbers and IR n ≥0 the set of n-tuples of non-negative real numbers; we conventionally see any such an n-tuple m = (m 1 , . . . , m n ) as a row-vector, i.e. a 1 × n matrix. We let

U n ⊂ IR n ≥0 be the unit simplex of IR n ≥0 , that is U n = {m ∈ [0, 1] n | m 1 + . . . + m n = 1}.
For function f : IR n → IR p continuous and twice differentiable, with f (m) = (f 1 (m), . . . , f p (m)), we denote by Df and D 2 f its first and second derivatives, respectively. Df (m) is the p × n (function) matrix such that (Df (m)

) ij = ∂fi(m) ∂mj . D 2 f (m) is the p × n × n tensor such that (D 2 f (m)) ijk = ∂ 2 fi(m) ∂mj ∂m k .
Moreover, for a p × n × n tensor P and a n × n matrix Q, we let P • Q be the (row) vector in IR p such that (P • Q) i = n j,k=1 (P ) ijk (Q) jk , for i = 1, . . . , p. In addition, for p × n matrix A and n × m matrix B we use the notation AB for standard matrix product: (AB) ij = n k=1 (A) ik (B) kj ; obviously this includes the case of vector inner product u T v for u and v being n × 1 (column) vectors, where A T is the transpose of A, i.e. (A T ) ij = (A) ji ; the vector outer product u ⊗ v, sometimes also denoted by u v T , is the matrix such that (u ⊗ v) i,j = u i v j .

Finally, for any vector v, we let v denote the norm of v and v 2 denote the square of v . Note that the choice of the specific norm is left unspecified because the results presented in the present paper hold for any norm 2.2. Synchronous Mean Field Model. We consider a system of 0 < N ∈ IN identical 2 interacting objects; N is called the size of the system. We assume that the number of local states of each object is finite 3 , say n; for the sake of simplicity, in the sequel we let the set of states of local objects be {1, . . . , n}, when not specified otherwise. Time is discrete and the behaviour of the system is characterised by a (time homogeneous) discrete time Markov chain (DTMC)

X (N ) (t) = (X (N ) 1 (t), . . . , X (N )
N (t)), where X (N ) i (t) is the state of object i at time t, for i = 1, . . . , N .

We define the occupancy measure at time t as the row-vector

M (N ) (t) = (M (N ) 1 (t), . . . , M (N )
n (t)) where, for j = 1, . . . , n, M (N ) j (t) is the fraction of objects in state j at time t, over the total population of N objects:

M (N ) j (t) = 1 N N i=1 1 {X (N ) i (t)=j}
where 1 {x=j} is equal to 1 if x = j and 0 otherwise. At each time step t ∈ IN each object performs a local transition, that may also be a self-loop to its current state. The transition probabilities of an object state depend on the current local state of the object and may depend also on M (t).

We let K(m) denote the one-step transition probability n×n matrix of an object in the system: K ij (m) is the probability for the object to jump from state i to state j in the system when the occupancy measure vector is m. We assume that, given the occupancy measure, the transitions made by two objects are independent. Our model is identical to the one of [START_REF] Boudec | A generic mean field convergence result for systems of interacting objects[END_REF] up to the fact that the authors of [START_REF] Boudec | A generic mean field convergence result for systems of interacting objects[END_REF] add a continuous resource to the model and allow object transition matrix K to depend also on the size N of the system (in which case they assume that the sequence of transition matrices K N converges to a function K as N goes to infinity). To 1 Of course, technical details in the proofs may depend on the specific norm (see e.g. footnote 4 in the proof of Lemma 4)

2 It is worth pointing out here that the requirement that the N objects be identical can be relaxed, since a system with different classes of identical objects can easily be modelled by considering an equivalent system with instances of an object whose set of states is the union of those of the original objects and similarly for the set of its transitions, as shown in the example in Section 5. 3 In fact, the same theoretical results could be derived for infinite dimensional models with two additional assumptions to cope with the fact that the simplex U ∞ is not compact : (i) imposing all functions to be uniformly continuous and (ii) Imposing tightness assumptions.

simplify the exposition, in this paper we consider a case without resource and we assume K(m) is a continuous function of m that does not depend on N . The results presented in this paper could be extended to the more general case where K(m) is also a function of N and could be modified to incorporate a resource, essentially by replacing our equation for the variance Γ by the one presented in [START_REF] Gast | A mean field approach for optimization in discrete time[END_REF]Equation (7)].

Below we recall Theorem 4.1 of [START_REF] Boudec | A generic mean field convergence result for systems of interacting objects[END_REF] on classical mean field approximation, under the simplifying assumptions mentioned above: Theorem 4.1 of [START_REF] Boudec | A generic mean field convergence result for systems of interacting objects[END_REF] (Convergence to Mean Field) Assume that the initial occupancy measure M (N ) (0) converges almost surely to the deterministic limit µ(0). Define µ(t) iteratively by (for t ≥ 0):

µ(t + 1) = µ(t) K(µ(t)). ( 2
)
Then for any fixed time t, almost surely:

lim N →∞ M (N ) (t) = µ(t).
In the sequel, we will write M (t) or simply M instead of M (N ) (t), leaving N and t implicit, when this does not cause confusion.

Refined Deterministic Approximation Theorem

In this section, we present our main results (Theorem 1 and Theorem 2) and provide their proofs.

3.1. First Main Result : Transient Behaviour. The iterative procedure of Theorem 4.1 of [START_REF] Boudec | A generic mean field convergence result for systems of interacting objects[END_REF] can be formalised as a t-indexed family of functions Φ t from U n to U n , where the functions Φ t : U n → U n are defined as follows:

Φ 0 (m) = m; (Φ 1 (m)) j = n i=1 m i K ij (m); Φ t+1 (m) = Φ 1 (Φ t (m))
This implies that for all m ∈ U n and t ∈ IN, we have Φ 1 (Φ t (m)) = Φ t (Φ 1 (m)).

In the following we assume that Φ 1 is continuous and twice differentiable with respect to m and that its second derivative is continuous (note that as U n is compact, this implies that Φ 1 and its first two derivative are uniformly continuous). Moreover, in what follows, we assume that M (N ) (0) converges to µ(0) (a deterministic value) as N goes to infinity and we let µ(t) be defined as in Equation ( 2), or, equivalently, µ(t + 1) = Φ 1 (µ(t)) = Φ t+1 (µ(0)).

Our main theorem can be stated as follows.

Theorem 1. Assume that the function Φ 1 is twice differentiable with continuous second derivative and that M (N ) (0) converges weakly to µ(0). Let A t and B t be respectively the n × n matrix A t = (DΦ 1 )(µ(t)) and the n × n × n tensor B t = (D 2 Φ 1 )(µ(t)). Then for any continuous and twice differentiable function with continuous second derivative h : U n → IR p ≥0 we have:

lim N →∞ N E h(M (N ) (t)) -h(Φ t (M (N ) (0))) = Dh(µ(t))V t + 1 2 D 2 h(µ(t)) • W t ,
where V t is an n × 1 vector and W t is an n × n matrix, defined as follows:

V t+1 = A t V t + 1 2 B t • W t W t+1 = Γ(µ(t)) + A t W t A T t , with V 0 = 0, W 0 = 0 and Γ(m) is the following n × n matrix: Γ jj (m) = n i=1 m i K ij (m)(1 -K ij (m)) Γ jk (m) = - n i=1 m i K ij (m)K ik (m)
The key idea of the proof is to use a Taylor expansion of E[h(Φ 1 (m))] around Φ 1 (m). We postpone the proof to Section 3.3.1.

One of the main consequences of Theorem 1 is that it allows us to compute precisely a development in O(1/N ) of the mean and the covariance of the vector M (N ) (t). This first order development is what we call the refined mean field approximation. In our numerical simulations, we will show that this refined approximation can greatly improve the accuracy of the original mean field approximation when the number of entities N is relatively small. Corollary 1. Let t ∈ IN. Under the assumptions of Theorem 1, and denoting µ(t) = Φ t (m), it holds that (i) For any coordinate i and any time-step

E M (N ) i (t) = µ i (t) + (V t ) i N + o 1 N .
(ii) For any pair of coordinates i, j, the co-variance satisfies

cov M (N ) i (t), M (N ) j (t) = 1 N (W t ) i,j + o 1 N .

3.2.

Second Main Result : Steady-State. Mean field approximation can also be used to characterise the steady-state behaviour of a population model. It has been shown that, in the case of continuous time or discrete-time mean field approximation, if this approximation has a unique attractor, then the stationary distribution of the system of size N concentrates on this attractor. In this section, we refine this result by computing the rate of convergence and by defining the refined approximation in this case. We say that a point µ(∞) is an exponentially stable attractor if

• For any m ∈ U n : lim t→∞ Φ t (m) = µ(∞) (i.e.

it is a global attractor).

• There exists an open neighbourhood

V of µ(∞) and two constants a, b such that all m ∈ V : Φ t (m) -µ(∞) ≤ ae -bt m -µ(∞) (i.e. it is exponentially stable).
Theorem 2. Assume that M (N ) has a unique stationary distribution (for each N ), that the function Φ 1 is twice differentiable and that the flow has a unique exponentially stable attractor µ(∞). Then there exists a n × 1 vector V ∞ and a n × n matrix W ∞ such that the constants V t and W t defined in Theorem 1 satisfy:

lim t→∞ V t = V ∞ and lim t→∞ W t = W ∞ Moreover (i) W ∞
is the unique solution of the discrete-time Lyapunov equation:

A ∞ W A T ∞ -W + Γ(µ(∞)) = 0 and V ∞ is uniquely determined by V ∞ = 1 2 (I -A ∞ ) -1 B ∞ W ∞ ,
where

A ∞ = DΦ 1 (µ(∞)), B ∞ = D 2 Φ 1 (µ(∞)
) and I is the identity matrix. (ii) for any twice differentiable function h, we can exchange the limits :

lim N →∞ lim t→∞ N E[h(M (t))] -h(Φ t (M (N ) (0))) = lim t→∞ lim N →∞ N E[h(M (t))] -h(Φ t (M (N ) (0))) = Dh(µ(∞))V ∞ + 1 2 D 2 h(µ(∞)) • W ∞ ,
This result is interesting for at least two reasons. First, it generalises Theorem 1 to the case of stationary distribution. Second, it is also the first result that provides a rate of convergence for the steady-state distribution of a model that has a discretetime mean field approximation (to the best of our knowledge, the rate of convergence had only been obtained for the finite-time horizon in [START_REF] Gast | A mean field approach for optimization in discrete time[END_REF]). We postpone the proof to Section 3.3.2.

Proofs.

To ease notation, in all the proofs we denote M N (0) by m, unless specified otherwise. In particular, when we write

E[h(M (t))] we formally mean E[h(M (N ) (t))|M (N ) (0) = m].

3.3.1.

Proof of Theorem 1. One of the key ingredients to prove our result is to study what happens for t = 1. This is what we do in Lemma 1. Then Theorem 1 will follow by using an induction on t. Note that one of the main technicalities in all these lemmas is to prove that the convergence is uniform in the initial condition. This is why we make use of various functions ε(N ) or ε t,g (N ) that control the uniform convergence to 0. Lemma 1. Let h : U n → IR p ≥0 (for p ≥ 1), be a twice differentiable function such that D 2 h is continuous. Then, there exists a function ε(N ) such that lim N →∞ ε(N ) = 0 and for all M (N ) (0) = m ∈ U n , the following holds:

N E[h(M (1))] -h(Φ 1 (m)) - 1 2 (D 2 h)(Φ 1 (m)) • Γ(m) ≤ ε(N ).
where Γ(m) is defined in Theorem 1.

Proof. The key idea of this proof is to consider the Taylor expansion of h for M (1) in the neighbourhood of Φ 1 (m). Let E = M (1) -Φ 1 (m). We get the following :

h(M (1)) -h(Φ 1 (m)) = (Dh)(Φ 1 (m))E + 1 2 (D 2 h)(Φ 1 (m)) • (E ⊗ E) + o(||E|| 2 )
Taking the expectation on both sides, we get :

E[h(M (1))] -h(Φ 1 (m)) = (Dh)(Φ 1 (m))E[E] + 1 2 (D 2 h)(Φ 1 (m)) • E[(E ⊗ E)] + E[o(||E|| 2 )]
The result follows by using Lemma 3, that establishes that

E[E] = 0 and N E[(E ⊗ E)] = Γ(m) and Lemma 4 that shows ||N E[o(||E|| 2 )]|| ≤ ε(N ).
The following lemma is a direct generalisation of Lemma 1 and it will be used in the proof of Theorem 1. The proof is essentially the same as that of Lemma 1 and exploits the time-homogeneity of the Markov chain.

Lemma 2. Let h : U n → IR p ≥0 be a twice differentiable function whose second derivative is continuous. Then there exists a function ε(N ) such that lim N →∞ ε(N ) = 0 and for all t ∈ IN, N ∈ IN >0 and M (N ) (t) = m ∈ U n , the following holds:

N E[h(M (t + 1))|M (t) = m ] -h(Φ 1 (m )) - 1 2 (D 2 h)(Φ 1 (m )) • Γ(m ) ≤ ε(N ) Lemma 3.
Under the assumptions of Lemma 1, we obtain that

E = E[M (N ) (1) -Φ 1 (m) | M (N ) (0) = m] satisfies E[E] = 0 and E[E ⊗ E] = Γ(m)/N .
Proof. We observe that by definition of our model, M j (1) is the following random variable:

(3)

M j (1) = 1 N n i=1 B ij
where ( B i,. ) is a random vector with multinomial distribution, with parameters N m i and (K i,• (m)). The variables are independent for different values of i (in particular, if i = i we have cov B ij , B i k = 0). Moreover, for all i and all j = k :

E B ij = N m i K ij (m) var B ij = N m i K ij (m)(1 -K ij (m)) cov B ij , B ik = N m i K ij (m)K ik (m).
This implies that

E[M j (1)] = E 1 N n i=1 B ij = 1 N n i=1 E B ij = Φ 1 (m) j . (4)
The case of E ⊗ E makes again use of Equation [START_REF] Bortolussi | Bounds on the deviation of discrete-time markov chains from their mean-field model[END_REF]. Note that by ( 4),

E j = M j (1)-Φ 1 (m) j = M j (1)-E[M j (1)]. This shows that E[(E ⊗ E) jk ] = cov (M j (1), M k (1)),
i.e. the covariance of M j (1) and M k (1). We consider the case k = j and the case k = j separately.

Case

k = j. N E [(E ⊗ E) jj ] = N var (M j (1)) = N var 1 N n i=1 B ij = 1 N n i=1 var B ij = n i=1 m i K ij (m)(1 -K ij (m)),
where the one but last equality comes from the independence of the variables B ij for i ∈ {1 . . . n}.

Case k = j. This case is similar.

N E [(E ⊗ E) jk ] = N cov (M j (1), M k (1)) = N n i=1 n i =1 1 N 2 cov B ij , B i k = - n i=1 m i K ij (m)K ik (m),
where in the double sum, only the terms i = i are non-zero because B ij and B i k are independent when i = i .

Lemma 4. Under the assumptions of Lemma 1 and using the notations of Lemma 1, there exists a function ε(N

) such that lim N →∞ ε(N ) = 0 and that ||N E[o(||E|| 2 )]|| ≤ ε(N ).
Proof. First of all we note that, as D 2 h is continuous, it is uniformly continuous (because U is compact). Hence, the term o(||E|| 2 ) ∈ IR p is uniform in m and ||E|| 2 , i.e., there exist a function δ : IR n → IR and a constant γ > 0 such that lim ||e||→0 δ(e) = 0, δ(e) ≤ γ for all e ∈ IR n , and

||o(||E|| 2 )|| ≤ ||E|| 2 δ(E). This implies E[||o(||E|| 2 )||] ≤ E[||E|| 2 δ(E)].
The proof proceeds with the following derivation: First, note that lim ||e||→0 δ(e) = 0 implies that for all > 0, there exists a > 0 such that δ(e) ≤ for all e such that ||e|| ≤ a . Therefore

E[||E|| 2 δ(E)] ≤ E[||E|| 2 δ(E)1 ||E||≥a ] + E[||E|| 2 1 ||E||<a ] ≤ η γ E[1 ||E||≥a ] + E[||E|| 2 ],
for some constant4 η < ∞.

As indicated by Equation (3),

E j = i ( B ij /N -m i K ij (m)) is the sum of the n independent random variables ( B ij /N -m i K ij (m)) and B ij has a binomial distribution of parameters (N m i , K ij (m)).
Hence, E j can be expressed as the sum of N independent Bernoulli random variables. By Hoeffding's inequality,

P [ E j ≥ t] ≤ e -2N t 2 .
This implies that

P[ E ≥ a ] ≤ n j=1 P[ E j ≥ a /n] ≤ ne -2N a /n 2 . Moreover, by Lemma 3, E[||E|| 2 ] ≤ 1/N . This shows that E[||E|| 2 δ(E)] ≤ η γ P(||E|| ≥ a ) + E[||E|| 2 ] ≤ η γ ne -2N a /n 2 + N
The assert follows by using ε(N ) = inf >0 (N η γ ne -2N a /n 2 + ).

Lemma 5. Under the assumptions of Lemma 1 and that M (N ) (0) converges weakly to µ(0) as N goes to infinity, for any continuous function g : U n → IR p ≥0 and all t, there exists a function ε t,g (N ) such that lim N →∞ ε t,g (N ) = 0 and

E[g(M (t))] -g(µ(t)) ≤ ε t,g (N ).
Proof. We proceed by induction on t. The lemma holds for t = 0 because M (N ) (0) converges weakly to µ(0). As g is continuous, there exists a function δ : IR + → IR + such that g(m) -g(m ) ≤ δ( m -m ) and lim r→0 δ(r) = 0. Moreover, as g and Φ t are continuous, g •Φ t is also uniformly continuous (and the continuity is uniform since U n is compact). Hence

E [ g(M (t + 1)) -g(µ(t + 1)) ] = E [ g(M (t + 1)) -g(Φ 1 (M (t))) + g(Φ 1 (M (t))) -g(µ(t + 1)) ] ≤ E [δ( M (t + 1) -Φ 1 (M (t)) )] + E [g • Φ 1 (M (t)) -g • Φ 1 (µ(t))] ≤ E [E [δ( E ) | M (t)]] + ε t,g•Φ1 ,
where E = M (t+1)-Φ 1 (M (t)) converges to 0 (uniformly in M(t)) by Lemma 3.

Of the main theorem. We proceed by induction on t. The theorem clearly holds for t = 0 because Φ 0 (M (N ) (0)) = M (N ) (0) by definition of Φ 0 . Assume that the theorem now holds for some t ≥ 0. We have :

N (E[h(M (t + 1))] -h(µ(t + 1))) =N E[h(M (t + 1)) -h(Φ 1 (M (t)))] + N (E[h(Φ 1 (M (t)))] -h(µ(t + 1))).
We will analyse the two lines separately. For the first line, the idea is to use Lemma 1. Indeed this line is equal to

E [N E [h(M (t + 1)) -h(Φ 1 (M (t))) | M (t)]] = E 1 2 (D 2 h)(Φ 1 (M (t))) • Γ(M (t)) + Θ(N ),
where by Lemma 1, Θ(N ) is such that Θ(N ) ≤ ε(N ). By Lemma 5 with g = (D 2 h)(Φ 1 ) • Γ, as N goes to infinity, this quantity converges to

1 2 (D 2 h)(Φ 1 (µ(t))) • Γ(µ(t)) = 1 2 (D 2 h)(µ(t + 1)) • Γ(µ(t)). (5)
For the second line, the idea is to apply the induction hypothesis to h • Φ 1 which can be done because the h • Φ 1 is twice differentiable (because both h and Φ 1 are). This shows that

N (E[h(Φ 1 (M (t)))] -h(µ(t + 1))) = N (E[h(Φ 1 (M (t)))] -h(Φ 1 (µ(t)))) = D(h • Φ 1 )(µ(t))V t + 1 2 D 2 (h • Φ 1 )(µ(t)) • W t + ε t,h•Φ1 (N )
The first term can be dealt with by applying the chain rule D(h•Φ 1 ) = (Dh)(Φ 1 )(DΦ 1 ) which shows that:

D(h • Φ 1 )(µ(t))V t = (Dh)(Φ 1 (µ(t)))(DΦ 1 )(µ(t))V t = (Dh)(µ(t + 1))(DΦ 1 )(µ(t))V t = (Dh)(µ(t + 1))A t V t . (6)
For the second term, we apply the product rule and again the chain rule:

1 2 D 2 (h • Φ 1 ) • W t = 1 2 D (Dh)(Φ 1 )(DΦ 1 ) • W t = D (Dh)(Φ 1 ) • (DΦ 1 ) + (Dh)(Φ 1 ) • D(DΦ 1 ) • 1 2 W t = (D 2 h)(Φ 1 ) • (DΦ 1 )(DΦ 1 ) T + (Dh)(Φ 1 )(D 2 Φ 1 ) • 1 2 W t
By applying the last function at the point µ(t) we get that:

1 2 D 2 (h • Φ 1 )(µ(t)) • W t =(D 2 h)(Φ 1 (µ(t))) • (DΦ 1 )(µ(t)) 1 2 W t (DΦ 1 ) T (µ(t)) + (Dh)(Φ 1 (µ(t)))(D 2 Φ 1 (µ(t))) • 1 2 W t ,
which, using the definition of µ(t + 1) = Φ 1 (µ(t)) and the assumptions A t = (DΦ 1 )(µ(t)) and B t = (D 2 Φ 1 (µ(t))), is the same as:

1 2 (D 2 h)(µ(t + 1)) • A t W t A T t + 1 2 (Dh)(µ(t + 1))(B t • W t ). ( 7
)
The theorem holds by combining Equations ( 5), ( 6) and [START_REF] Braverman | Steins method for steady-state diffusion approximations of M/P h/n + M systems[END_REF].

3.3.2. Proof of Theorem 2.
Proof. The proof is inspired by the proof of [13, Theorem 3.1] and uses ideas of Stein's method. Because many details are similar to the proof of Theorem 1, we skip some details of computation in this proof.

Let h be a twice-differentiable function and let G h be the function defined for all m by:

G h (m) = ∞ t=0 [h(Φ t (m)) -h(µ(∞))]. G h (m) is well defined because µ(∞) is exponentially stable attractor.
By construction, for any m we have

G h (m) = h(m) -h(µ(∞)) + ∞ t=1 [h(Φ t (m)) -h(µ(∞))] = h(m) -h(µ(∞)) + G h (Φ 1 (m))
The above equation is a discrete time Poisson equation and implies that for any m:

h(m) -h(µ(∞)) = G h (m) -G h (Φ 1 (m)) (8)
Assume that at time 0, the initial state M (0) is distributed according to the stationary distribution of the system of size N . By the definition of stationarity, at time 1, M (1) is also distributed according to the same stationary distribution and we have

E[G h (M (0))] = E[G h (M (1))].
By using (8) and then the above equation, we get:

E [M (0)] -h(µ(∞)) = E [M (0) -h(µ(∞))] = E [G h (M (0)) -G h (Φ 1 (M (0)))] = E [G h (M (1)) -G h (Φ 1 (M (0)))]
N (E [M (∞)] -h(µ(∞))) = E [G h (M (1)) -G h (Φ 1 (M (0)))] = E 1 2 D 2 (G h )(Φ 1 (M (0))) • Γ(M (0)) + o(1) = 1 2 D 2 (G h )(Φ 1 (µ(∞))) • Γ(µ(∞)) + o(1),
where the last equality comes from the fact that the stationary distribution of the system of size N converges weakly to a Dirac measure in µ(∞) as N goes to infinity (see [START_REF] Gast | A mean field approach for optimization in discrete time[END_REF]Corollary 14]).

To conclude the proof, the only remaining step is to compute the second differential of G h which can be expressed as the infinite sum:

D 2 (G h )(µ(∞)) = ∞ t=0 D 2 (h • Φ t )(µ(∞)).
The expressions for V ∞ and W ∞ come from plugging the above equations into Equation ( 5), ( 6) and [START_REF] Braverman | Steins method for steady-state diffusion approximations of M/P h/n + M systems[END_REF]. The uniqueness of the solution of the Lyapunov equation is due to the fact that the fixed point µ(∞) is exponentially stable and therefore also linearly stable.

Refined Mean Field Model for SEIR

In this section we provide a simple example that illustrates the results for the refined mean field model of the simple computer epidemic SEIR example presented in [START_REF] Bortolussi | Continuous approximation of collective system behaviour: A tutorial[END_REF]. Each object in the model consists of four local states: Susceptible (S), Exposed (E), Infected (I) (and active) and Recovered (R). The four-state SEIR model of an individual object is shown in Figure 1.

Its discrete time evolution is given by the following probability transition matrix K in which m S , m E , m I and m R denote the fraction of objects in the system that are in local state S, E, I and R, respectively:

K(m S , m E , m I , m R ) =     1 -(α e + α i m I ) α e + α i m I 0 0 0 1 -α a α a 0 0 0 1 -α r α r α l 0 0 1 -α l    
In other words, a susceptible becomes exposed with probability (α e + α i m I ) -i.e., α e denotes the external and α i the internal infection probability -; An exposed node activates his infection with probability α a ; An infected recovers with probability α r ; and α l is the probability to loose the protection against infection.

4.1. Computation of A, B and Γ. We illustrate how to apply Theorem 1 in its simplified form, when h is the identity function, as in Corollary 1-(i). The first step is to compute the Jacobian and the Hessian of the function Φ 1 for a generic occupancy measure vector m at time step t. Written as a column vector, the function Φ 1 (m) = mK(m) is given by

Φ 1 (m) =     m S (1 -α e -α i m I ) + α l m R m S (α e + α i m I ) + (1 -α a )m E m E α a + (1 -α r )m I α r m I + (1 -α l )m R    
Hence, the Jacobian is the following 4 × 4 matrix:

D(Φ 1 )(m S , m E , m I , m R ) =     1 -(α e + α i m I ) 0 -α i m S α l α e + α i m I 1 -α a α i m S 0 0 α a 1 -α r 0 0 0 α r 1 -α l    
The Hessian is a 4 × 4 × 4 tensor. We provide them as 4 matrices of 4 × 4, one for each function (Φ 1 ) j , where j ∈ {S, E, I, R}:

D 2 ((Φ 1 ) S )(m S , m E , m I , m R ) =     0 0 -α i 0 0 0 0 0 -α i 0 0 0 0 0 0 0     D 2 ((Φ 1 ) E )(m S , m E , m I , m R ) =     0 0 α i 0 0 0 0 0 α i 0 0 0 0 0 0 0    
The matrices for I and R are two 4×4 zero-matrices. The 4×4 matrix Γ depends on K and the occupancy measure m as defined in Lemma 1. The refined mean field approximation of the occupancy measure vector is thus given by E[M (N ) (t)] ≈ µ(t) + V t /N , where V t is computed recursively, according to Theorem 1.

4.2.

Dynamics of the SEIR Model and its Approximations. We consider a model with the following parameter values for the local transition probabilities: α e = 0.01, α i = 0.08, α r = 0.02, α l = 0.01 and α a = 0.04. Initially, M (0) = (0.2, 0.2, 0.2, 0.4). Figure 2 shows the results for the classical mean field approximation, the refined mean field approximation and the average of 100,000 runs of a stochastic simulation of the model obtained. The results are given for population size N = 10, N = 20, N = 50 and N = 100, respectively; time t ranges from 0 to 500 time units.

We observe that, as exepcted, the gap between the classical mean field approximation and the simulation is relatively small and decreases with N . Still, for N = 10, we observe a clear difference between the classical mean field approximation and the simulation, whereas the refined mean field provides a much closer approximation (in this case, the graphs overlap almost everywhere). With the increase of the population size N both approximations converge to the same value, as well as the value obtained by simulation: for N ≥ 50, the curves are almost indistinguishable.

N = 50 N = 100
To highlight the differences, we plot in Figure 3 the difference between the two approximations with respect to the simulation : On the left panel, we plot a function of time for the quantities E[M (t)]-µ(t); On the right we plot E[M (t)]+V t /N -µ(t) (in both cases for N = 10). We observe that the refined mean field approximation (right panel) is an order of magnitude closer to the value obtained by simulation: while the error of the classical mean field approximation can be larger than 0.05, the error of the refined mean field approximation remains always smaller than 0.01.

These two figures illustrate that Theorem 1 is not just valid asymptotically, but actually it refines the classical mean field approximation for relatively small values of N . To go further, we study the steady-state distribution in Table 1 in which we display the average proportion of objects in states S, E, I or R estimated by simulation, refined mean field approximation and classical mean field approximation. This illustrates the approximation accuracy for the steady state of the SEIR example for each local state of an object. As for the two previous figures, this table illustrates that the refined mean field approximation provides very accurate estimates of the true stationary distribution even for very small values of N , which shows that the asymptotic results presented in Theorem 2 are also useful for small values of N . These results are in line with the results presented in [START_REF] Gast | A refined mean field approximation[END_REF] for continuous time mean field models. Table 1. SEIR model: Comparison of the accuracy of the mean field and refined mean field approximation. The columns show the average proportion of objects in the states susceptible (S), exposed (E), infected (I) and recovered (R), respectively. Each item in the table was computed by measuring the occupancy measure for times t = 1000, i.e. when the systems' occupancy measure has reached a sufficiently stable value. Simulation values are averages over 100, 000 simulations.

Refined Mean Field Model for WSN

The next example concerns a simple model of a wireless sensor network [START_REF] Bortolussi | Bounds on the deviation of discrete-time markov chains from their mean-field model[END_REF]. Such networks are composed of wireless sensor nodes and gateways. This example serves two purposes. First it shows that the assumption of homogeneous objects is not restrictive: In this example, there are two classes of objects, which is represented by having a block-diagonal matrix K. Second, we use it to consider a function h that is not just the projection on one coordinate.

Wireless sensor nodes have three local states. In the initial state, e, a sensor node waits for detecting an event of interest and collects data for that event. After that, the node moves to state c to communicate its data to an available gateway. The communication attempt may timeout if no gateway is available. In that case Gateway nodes have two states. Initially they are in state a and available to receive data from a sensor node. Upon connection to a sensor node they move to state b during which they are busy processing the data. When in state b they are temporarily unavailable for communication with other sensor nodes. After processing the batch of data they move back to state a.

We consider a model where objects have five local states {a, b, c, d, e}, where a and b are states of a gateway node and c, d and e are states of a sensor node, i.e. each object in the model can behave either as a gateway or as a sensor, but it cannot change its behaviour from that of a gateway to that of a sensor, or viceversa. A system is then composed of N (syntactical) homogeneous objects with a fixed fraction m G = m a +m b of gateway nodes and a fraction m W = m c +m d +m e of wireless sensor nodes, such that m W = 1 -m G . To keep the model simple for the purpose of illustrating the refined mean field approach, we do not consider interference due to collision in the communication between nodes and gateways.

The probability transition matrix is given below:

K(m) =       1 -βm c βm c 0 0 0 α 1 -α 0 0 0 0 0 1 -γ -βm a γ βm a 0 0 η 1 -η 0 0 0 λ 0 1 -λ      
, where α denotes the probability of the gateway to get again available, β the probability of data communication between the gateway and a sensor node, λ the probability that a sensor node is ready to send data, γ the probability that a sensor node performs a time-out and η the probability that a delayed sensor node tries to communicate again.

In the example we will use the following values for the above parameters: α = 0.09, β = 0.9, λ = 0.09, γ = 0.01, and η = 0.01, and let M (N ) (t) denote, as usual, the occupancy measure process of the WSN model (leaving N and t implicit for the sake of notation simplicity). We are interested in the average response time of a sensor node, i.e. the time a sensor node needs to wait to be able to communicate its data to the gateway. This expected response time can be defined as the fraction between the sensor nodes that are already waiting to communicate their data, i.e. the sensor nodes in state c and state d, and the new sensor nodes that became ready in the current time step, i.e. λ times the nodes in local state e:

E[R] = E (M c + M d ) λM e ,
With reference to Theorem 1, we define h(x 1 , x 2 , x 3 , x 4 , x 5 ) = (x3+x4) λx5 . 5.1. Computation of A, B and Γ for the WSN model. In the sequel, we make reference to m = (m a , m b , m c , m d , m e ) ∈ U 5 . The Jacobian of function Φ 1 :

D(Φ 1 )(m) =       1 -βm c α -βm a 0 0 βm c 1 -α βm a 0 0 -βm c 0 1 -γ -βm a η λ 0 0 γ 1 -η 0 βm c 0 βm a 0 1 -λ      
The Hessian of the function Φ 1 satisfies D 2 ((Φ 1 ) d )(m) = 0:

D 2 ((Φ 1 ) a )(m) = D 2 ((Φ 1 ) c )(m) =       0 0 -β 0 0 0 0 0 0 0 -β 0 0 0 0 0 0 0 0 0 0 0 0 0 0       D 2 ((Φ 1 ) b )(m) = D 2 ((Φ 1 ) e )(m) =       0 0 β 0 0 0 0 0 0 0 β 0 0 0 0 0 0 0 0 0 0 0 0 0 0       The Jacobian of function h is D(h)(m) = (0, 0, 1 λme , 1 λme , -mc +m d λm 2 e
) and its Hessian is

D 2 (h)(m) =        0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 λm 2 e 0 0 0 0 -1 λm 2 e 0 0 -1 λm 2 e -1 λm 2 e 2 * (mc +m d ) λ * m 3 e        .
The 5 × 1 vector V t and the 5 × 5 matrix W t are computed recursively according to Theorem 1, using the new Jacobian and Hessian for function Φ; thus the refined mean field approximation of the measure of interest is given by

E[h(M (N ) (t))] ≈ h(Φ t (m)) + (D(h) • V t + 1 2 * (D 2 (h) • W t ))/N . 5.2.
Results. In Figure 5 various approximations of the expected response time for the WSN model are shown, for time values t ranging from 0 to 400 time units. We consider a relatively small system with 15 nodes (10 sensor and 5 gateway nodes). Recall that the function h is defined by : h(x 1 , x 2 , x 3 , x 4 , x 5 ) = (x3+x4) λx5 . We compare five curves :

(1) The blue curve labelled Classic Mean Field (1) (obtained with Octave) shows the expected response time when this is approximated by defining E[R] as in [START_REF] Bortolussi | Bounds on the deviation of discrete-time markov chains from their mean-field model[END_REF]:

E[R] = (m c + m d ) λm e = h(m),
where m c , m d and m e denote the classical mean field approximation values for the fractions of sensor nodes being in state c, d and e, respectively. (2) The red curve (2) is the expectation of h(M ), computed by stochastic simulation:

E [h(M )] = min(E (M c + M d ) λM e , 100)
In the case of individual simulation runs it may of course happen that M e = 0 occasionally. This is why we defined E[h(M )] as the minimum between the actual value and 100. The latter is the value one obtains when one but all nodes are waiting and the last node is getting ready for communication too: M c + M d = 9 and M e = 1, i.e. 9/(0.09*1) = 100. (3) The orange curve (3) shows the expected response time approximated using the refined mean field approximation of Theorem 1 with the function h. For comparison, we also compute two other quantities :

4. The purple curve (4) shows the response time approximated as follows, using the refined mean field approximation for the fraction of sensor nodes in each state (i.e. we use Theorem 1 with the identity function and then apply h):

h(rmf ) = (rmf c + rmf d ) λrmf e
where rmf c = µ c + V c /N denotes the refined mean field approximation of the fraction of sensor nodes in state c, and similarly for rmf d and rmf e . 5. Finally, the green curve [START_REF] Boudec | A generic mean field convergence result for systems of interacting objects[END_REF] shows the expected response time defined as:

h(E [M ]) = (E[M c ] + E[M d ]) λE[M e ]
where E[M c ] is the average fraction of sensor nodes in state c obtained via the average of 100,000 individual simulation runs of the model. E[M d ] and E[M e ] are obtained in a similar way. In Figure 5(a), we plot these various curves for 15 nodes in total. We make two observations. First, in this case the value obtained by simulation (red curve (2)) is almost 50% larger than the classic mean field approximation (1) whereas the refined approximation (3) is much closer. Second, the purple curve (4) is close to the green curve (5) but quite far away from the red (2) and orange (3) curves. This shows that when applying Theorem 1, computing a refined model for E[h(M )] and for h(E[M ]) might lead to very different results.

Of course, the larger N gets (in an otherwise equal model), the closer the orange (3) and red (2) curves will get to the blue curve (1), i.e. the classic mean field approximation, as illustrated in Figure 5(b). In both cases, all curves collapse into to a single curve.

Refined Mean Field Model for Majority Rule Decision-making

The example in this section concerns a model for collective decision-making. The model is inspired by the work of Montes de Oca et al. (see [START_REF] Antonio Montes De Oca | Majority-rule opinion dynamics with differential latency: a mechanism for self-organized collective decision-making[END_REF][START_REF] Scheidler | Dynamics of majority rule with differential latencies[END_REF] and references therein). Collective decision-making is a process whereby the members of a group decide on a course of action by consensus. Such collective decision-making processes have also been applied in swarm robotics. In particular, in that context the robots where asked to choose between two actions that have the same effect but differ in their execution times [START_REF] Antonio Montes De Oca | Majority-rule opinion dynamics with differential latency: a mechanism for self-organized collective decision-making[END_REF].

One strategy of collective decision making is the use of the majority rule. In this strategy the agents in a population are initially divided into two groups. One in which all members have opinion A and one where all members have opinion B. In every step three agents are selected randomly from the total population to form a temporary team. The team applies the majority rule such that all its members adopt the opinion held by the majority (i.e. at least two) of the members, after which they return to the total population until the population has reached a consensus on one of the two opinions.

In the majority rule strategy extended with differential latency the agents in the population are not all the time available for team formation. Both types of agents are assumed to perform an action with a certain duration during which they cannot participate in team formation. For example, agents with opinion B perform such actions taking (on average) relatively more time than those with opinion A. In [START_REF] Antonio Montes De Oca | Majority-rule opinion dynamics with differential latency: a mechanism for self-organized collective decision-making[END_REF] such latency periods for agents with opinion A and B are modelled by random variables with exponential distributions with rate λ A and λ B respectively. For simplicity, it can also be assumed that the A-type actions take 1 time unit on average (i.e. λ a = 1) and that B-type actions take 1/λ time units on average, where λ takes a value in (0, 1]. This variant of self-organised collective decisionmaking is known as majority rule with differential latency (MRDL).

In the following we develop a probabilistic, discrete time variant of the MRDL strategy which we call MRDL-DT. In this variant agents can have either opinion A or opinion B, and in both cases they can be either be latent or not, leading to a partition of the population into exactly four classes: LA (latent A), NA (non-latent A), LB (latent B), NB (non-latent B). It is assumed that if an agent is latent it cannot be selected for team formation. The four state MRDL-DT model of an individual object of the population is shown in Figure 6. The name of the states indicate in which class the object is. The behaviour of an individual object is as follows. Initially the object is latent, and, assuming it has opinion A (state LA), it finishes its job and becomes available for team-formation (transition actA) moving to state NA with probability 1/q, for appropriate q. When in NA it gets selected in a team with two other members. If the two other members have opinion B it changes its opinion into B and moves to state LB . This can happen with a probability 3m 2 NB where factor 3 models the fact that we abstract from the exact order in which the members of the team are selected, which can happen in 3 different ways. In alternative the two other members can have both opinion A, or one opinion A and the other opinion B. In that case the opinion of the object does not change and the object moves back to LA with probability 3 q (m 2 NA + m NA m NB ). If the object is in state LB it becomes available for team formation, moving to state NB with a probability λ/q, where λ is a value in (0, 1]. The latter models the relative longer duration of activity B with respect to A. The behaviour in state NB is similar to that in state NA, except that now the opinion may change from B to A.

The discrete time evolution of the model is given by probability transition matrix K in which m LA , m NA , m LB and m NB denote the fraction of objects in the system that are in local state LA, NA, LB and NB , respectively and m = (m LA , m NA , m LB , m NB ):

K(m) =      1 -1 q 1 q 0 0 3 q (m 2 NA + m NA m NB ) NA(m) 3 q m 2 NB 0 0 0 1 -λ q λ q 3 q m 2 NA 0 3 q (m 2 NB + m NA m NB ) NB (m)     
where

NA(m) = 1 - 3 q (m 2 NB + m 2 NA + m NA m NB )
and

NB (m) = 1 - 3 q (m 2 NA + m 2 NB + m NA m NB ).
Since we are dealing with clock-synchronous discrete systems, we also introduced the discretisation factor q = 10 so that only a fraction of the population is moving from the latent to the non-latent state at any time.

In the example we use the values q = 10 and λ taking values 1.0, 0.5 and 0.25 in the various analyses, modelling that task B takes the same time as task A, or twice as much time or four times as much time as task A, on average, respectively. We are interested in the evolution of the consensus, C A , on opinion A as a function of the initial values and the differential latency λ. Let

E[C A ] = E[M LA + M N A ] = E[h(M LA , M N A , M LB , M N B )]
where, with reference to Theorem 1, h(x 1 , x 2 , x 3 , x 4 ) = x 1 + x 2 .

6.1. Computation of A, B and Γ for the MRDL-DT model. As before, we first need to compute the Jacobian and the Hessian of the function Φ 1 for a generic occupancy measure vector m at time step t. The Jacobian of function Φ 1 is:

D(Φ 1 )(m) =      1 -1/q 9 q m 2 NA + 12 q m NA m NB 0 6 q m 2 NA 1/q JNA(m) 0 -3 q m 2 NA -6 q m NA m NB 0 6 q m 2 NB 1 -λ q 12 q m NB m NA + 9 q m 2 NB 0 -3 q m 2 NB -6 q m NA m NB λ q JNB (m)     
where

JNA(m) = 1 -( 9 q m 2 NA + 6 q m NA m NB + 3 q m 2 NB )
and

JNB (m) = 1 -( 3 q m 2 NA + 9 q m 2 NB + 6 q m NA m NB ).
The Hessian of function Φ 1 is:

D 2 ((Φ 1 ) LA )(m) =     0 0 0 0 0 18 q m NA + 12 q m NB 0 12 q m NA 0 0 0 0 0 12 q m NA 0 0     D 2 ((Φ 1 ) N A )(m) =     0 0 0 0 0 -18 q m NA -6 q m NB 0 -6 q m NA -6 q m NB 0 0 0 0 0 -6 q m NA -6 q m NB 0 -6 q m NA     D 2 ((Φ 1 ) LB )(m) =     0 0 0 0 0 0 0 12 q m NB 0 0 0 0 0 12 q m NB 0 18 q m NB + 12 q m NA     D 2 ((Φ 1 ) N B )(m) =     0 0 0 0 0 -6 q m NB 0 -6 q m NA -6 q m NB 0 0 0 0 0 -6 q m NA -6 q m NB 0 -18 q m NB -6 q m NA     6.2.
Results for the MRDL-DT example. We first show some results for a medium size population of N = 160. Figure 7 shows the dynamics of the fractions of the population having opinion A and B, for both the latent and non-latent objects, for the first 200 time units. Similarly to the previous examples, a good correspondence can be observed between the results of the mean of 1000 simulation runs and the mean field approximations. Also in this case the refined mean field provides a better approximation than the classical mean field approximation for the period ranging approximately from 50 to 150 time units. However, if a much smaller population is considered, e.g. N = 32, the mean field approximation differs considerably from the simulation results and the refined approximation does not really improve the accuracy of the approximation. This is what can be observed in Figure 8, where the results for N = 32 are shown for a model without differential latency (i.e. λ = 1.0), for t ∈ [0, 500]. This can be explained as follows : For a large population, a system that is initially biased towards one opinion will reach consensus on this opinion. For a small population, however, a system that is initially biased towards one opinion still can reach consensus on the other opinion due to intrinsic stochastic fluctuations. This cannot be taken into account in a mean field population model. However, both analytical models, derived from a master equation approach, and simulation show that the probability to reach consensus on a given opinion rapidly converges to a step function for a growing population size N [START_REF] Scheidler | Dynamics of majority rule with differential latencies[END_REF], where the critical density is given by m A = λ (1+λ) , where m A is the initial fraction of the population with opinion A. For large N , if m A > λ (1+λ) the system almost surely reaches consensus on A, whereas for m A < λ (1+λ) it almost surely reaches consensus on B. This explains why for larger populations the mean field approximations become increasingly accurate as shown in Figure 7. This example shows a limit of the refined mean field approximation: when a system has multiple equilibrium point (and in particular when there are multiple absorbing states as in this example), the dynamics of the mean field approximation depends on the initial state of the system: For a given initial state the mean field will always follow the same trajectory (it is a deterministic system). When the system is large, the random fluctuations will remain small and the corresponding stochastic system will stay in the same basin of attraction. In the case of a small population, however, the dynamics will be greatly affected by the random fluctuations. These fluctuations can lead the system to another basin of attraction than the original one.

Non-Exponentially Stable Equilibrium : Accuracy versus Time

In the previous sections, the dynamical system m = Φ 1 (m) has either one exponentially stable attractor (Section 4) or multiple exponentially stable attractors (Section 5). When the attractor is unique and exponentially stable, the accuracy of the mean field or refined mean field approximation is uniform in time (Theorem 2(ii)). In this section, we study the case of a system that has a unique attractor but that is not exponentially stable. We show that, in this case, the accuracy of the (refined) mean field approximation is no longer uniform in time.

We consider a system with N objects in which each object is in state 0 or 1. An object in state 1 goes to state 0 with probability 1 and an object in state 0 goes to 1 with probability αm 0 , where α ∈ (0, 1) is a parameter. The transition matrix K is 7.1. Transient regime and accuracy for large t. In Figure 9 and Figure 10, we plot the first component of the mean field µ(t) and refined mean field approximation µ(t)+V (t)/N as well an exact value of E[M (t)] for N = 10 and N = 30. The initial value is m = (0.7, 0.3). The exact value of E[M (t)] was computed by a numerical method that uses the fact that the system with N objects can be described by a Markov chain with N + 1 states.

K(m) = 1 -αm 0 αm 0 1 0 The function Φ 1 : m → mK(m) has a unique fixed point whose first component is µ 0 (∞) = ( √ 1 + 4α -1)/(2α
These figures show that the refined approximation always improves the accuracy compared to the classical mean field approximation for small values of t, both for α = 0.6 and α = 0.75. The situation for large values of t is quite different. On the one hand, when the fixed point is exponentially stable (α = 0.6, Figure 9), the refined approximation is very accurate for all values of t. On the other hand, when the fixed point is not exponentially stable (α = 0.75, Figure 10), the refined approximation seems to be unstable and is not a good approximation of E[M (t)] for values of t that are too large compared to N (t > 7 for N = 10 or t > 12 for N = 30). 7.2. Steady-state convergence. To explore how the non-exponentially stable case affects the accuracy of mean field approximation, we now study in more details the steady-state convergence when α = 0.75. It is known (see for example [START_REF] Gast | A mean field approach for optimization in discrete time[END_REF]Corollary 14]) that when the dynamical system m = Φ 1 (m) has a unique attractor µ(∞), then the steady-state expectation E[M (N ) ] converges to µ(∞) as N goes to infinity. Theorem 2 shows that if in addition the attractor is exponentially stable then E[M (N ) ] ≈ µ(∞) + V /N .

In Figure 11, we show that the latter no longer holds when the mean field system has a unique attractor that is not exponentially stable. We consider the same model with α = 0.75 for which µ(∞) = 2/3. We plot in Figure 11 √ N (E[M (N ) ] -µ(∞)), where we computed E[M (N ) ] by inverting the transition matrix of the system of size N . This figure shows that √ N (E[M (N ) ] -µ(∞)) does not converge to 0 as N goes to infinity but seems to converge to approximately -0.0975 (as indicated by the fitted line in orange). This suggests that for this model, one has in steady-state :

E M (N ) ≈ µ(∞) - 0.0975 √ N + 0.14 N . (9)
Note that the constants -0.0975 and 0.14 were obtained by a purely numerical method that consist in finding the best curve of the form a + b/ √ N that fits √ N (E[M (N ) -µ(∞)]). For now, we do not known if there exists a systematic way to obtain these values for another model that would also have a non-exponentially equilibrium point. We left this question for future work.

We remark that the convergence of E[M (N ) ] to µ(∞) observed Equation ( 9) is in O(1/ √ N ) and not in O(1/N ). This model satisfies all the assumptions of Theorem 2 but one : The attractor µ(∞) is not exponentially stable. This suggests that having an exponentially stable attractor is needed to obtain a convergence in 1/N .

Conclusion

In this paper we studied population models composed of (clock-)synchronous objects. A classical method to study such systems is to consider the mean field approximation. By studying the accuracy of this deterministic approximation, we developed a new approximation, that we call the refined mean field approximation. We illustrated on a few examples that this approximation can greatly improve the accuracy of the classical mean field limit, also for systems with a relatively small size (10 -20 objects). Yet, this refined approximation has some limitations when the deterministic approximation has multiple basins of attraction or has a unique attractor that is not exponentially stable. The proposed refined approximation is given by a set of linear equations that scales as the square of the dimension of the model (but does not depend on the system size). For now, we limited our study to relatively small models, for which the Jacobian and Hessian can be computed in closed form. We are currently investing means to make this computation automatic which will allow us to study large-scale examples.

Figure 1 .

 1 Figure 1. SEIR model of individual object

Figure 2 .

 2 Figure 2. Evolution of the fraction of objects in state S for population sizes N = 10, N = 20, N = 50 and N = 100. The figures compare the classical mean field approximation (obtained with python-numpy) with the refined one and with the average of 10,000 simulation runs of the system.

  Error of the mean field approximation Error of the refined mean field approx.

Figure 3 .

 3 Figure 3. SEIR model: Quantification of the difference (error) between the simulation results and the classical mean field approximation (left) and the refined mean field result (right), respectively, for N = 10. The results show simulation value minus mean field value.

Figure 4 .

 4 Figure 4. WSN model of individual objects: Sensor Node (left) and Gateway (right)

Figure 5 .

 5 Figure 5. Expected response time E[R] for a sensor node to communicate its data to a GW for a WSN model with N nodes of which 2N/3 are sensor nodes and N/3 are gateway nodes. The red line (2) is an average over 20,000 simulations for N = 15 and 1000 for N = 1500.

Figure 6 .

 6 Figure 6. Majority rule differential latency model of an individual object. Latent states are red, non-latent ones blue.

Figure 7 .

 7 Figure 7. Dynamics of opinion A, latent and non-latent. Classical mean field (plain lines), refined mean field (dashed lines) and simulation results (dotted lines) of MRDL model with 160 objects, λ = 1.0, q=10 and initially 0.6 * 160 = 96 have opinion A and the population is initially latent.

Figure 8 .

 8 Figure 8. Classical mean field (plain lines), refined mean field (dashed lines) and simulation results (dotted lines) of MRDL model with 32 objects, λ = 1.0, q=10 and initially 0.6 * 32 = 9 have opinion A and the population is initially latent.

Figure 9 .

 9 Figure 9. Exponentially stable case (α = 0.6).

Figure 10 .

 10 Figure 10. Non-exponentially stable case (α = 0.75).

Figure 11 .

 11 Figure 11. Non-exponentially stable case (α = 0.75) : convergence of E[M (t)] to µ(∞).

  1 . Note that in the proofs, what we denote by E, o( E

2 ), E[o( E 2 )] are p-dimensional vectors while their norm E , E 2 , E[o( E 2 )

] ∈ IR are real non-negative numbers.

  ). This fixed point is exponentially stable if and only if α < 0.75.

		0.705							0.705
		0.700							0.700
	m(t)	0.695		Refined mean-field approx. Mean-field approx. Exact value of E[M(t)]	m(t)	0.695
		0.685 0.690							0.685 0.690	Refined mean-field approx. Mean-field approx. Exact value of E[M(t)]
		0	5	10	15 Time	20	25	30	0	5	10	15 Time	20	25	30

NICOLAS GAST (INRIA), DIEGO LATELLA (CNR-ISTI), AND MIEKE MASSINK (CNR-ISTI)

The specific value of η depends on the norm used; for instance η = 1 for the infinity norm.10NICOLAS GAST (INRIA), DIEGO LATELLA (CNR-ISTI), AND MIEKE MASSINK (CNR-ISTI)


Near shore: wave transformation due to interaction with complex bathymetries

W = η u , A = 0 1 1 0 , |A| = Id2
where we recall that (all . are removed for simplicity, and

c 2 0 = gd0) η = η dim d0 , u = u dim c0 , x = x dim λ , t = t dim λ/c0
1st order upwind method

∂tW + A W i+1 -W i-1 2∆x = ∆x 2 W i+1 -2W i + W i-1 ∆x 2 2nd order centered differencing ∂tW + A W i+1 -W i-1 2∆x = 0
1st order upwind and 2nd order centered differencing: modified equation/TE

∂tW smooth + A∂xW smooth = ∆x 2 ∂xxW smooth - A∆x 2 6 ∂xxxW smooth + O(∆x 3 ) ∂tW smooth + A∂xW smooth = - A∆x 2 6 ∂xxxW smooth + O(∆x 4 )
It looks like "if we could see the dispersive effects of the 1st order scheme, they would be the same as those of the second order".. but the numerical viscosity is too high...

Example: linearized shallow water

Turn down the viscosity: 3 rd upwind scheme

∂tW + A W R i+1/2 + W L i+1/2 2∆x - W R i+1/2 -W L i+1/2 2∆x -A W R i-1/2 + W L i-1/2 2∆x - W R i-1/2 -W L i-1/2 2∆x = 0
With "quadratically reconstructed " left and right values at

x i ± ∆x/2 W L i+1/2 =W i + W i -W i-1 6 + W i+1 -W i 3 W R i+1/2 =W i+1 - W i -W i-1 6 - W i+2 -W i+1 3
3 rd order upwind vs 4 th order central differencing: modified eq./TE

∂tW smooth + A∂xW smooth = - ∆x 3 12 ∂xxxxW smooth - A∆x 4 30 ∂xxxxxW smooth + O(∆x 5 ) ∂tW smooth + A∂xW smooth = - A∆x 4 30 ∂xxxxxW smooth + O(∆x 6 )
Low dissipation: O(∆x 3 ) viscosity (ok for hyperbolic + explicit time stepping) Same dispersion of the fourth order FD ! ht+qx = 0

qt+ (ũq)x + hηx = µ 2 h ψ ψ =α ∂x( h2 ∂x(ũt + ũũx)) + (α -1) ∂x( h2 ∂xx η) + Q ψ (ũ, h, d; ũx, hx, dx)
With the notation ∂th + ∇ • q = 0 ;

I + αT h ∂tq + ∇ • q ⊗ q h + gh∇η -T h gh∇η + hQ(u) = 0 .

With the notation

depth : h = d + η (vector) volume flux : q = hu with q the (dimensionless) depth averaged velocity vector T h the dispersive elliptic operator given by

T h (•) = hT ( • h )
where T is a self adjoint operator6 :

T (•) = S * (S(•)) , S(•) = h √ 3 ∇ • (•) - √ 3 -1 2 ∇b • (•)

a 2 .C 2 = C 2 0tanh(κ d 0 ) κ d 0 +O(µ 2 ) 2 )x

 222022 and inundation; 2004 Sumatra tsunami reaching the coast of Thailand Sea waves diffracting around a peninsula Rip current, Park beach (Coffs harbour, NSW Australia) a Closure model 1. Detect breaking regions 2. Remove dispersive terms 3. → shallow water shock 4. Total energy E = gh 2 /2 + hu 2 /2 5. Dissipation : D b = [F E -σE] ≈ γ b [H] 3 Need to handle shallow water: "Upwinding"/numerical dissipation Shock capturing in breaking regions etc All the std. artillery.... Principle of the talk Dispersive models for propagation Which dispersion dominates ? 1. Continuous dispersive models have a range of validity related to the model dispersion error Discrete dispersive models have a range of validity related to the scheme dispersion error How can we exploit this knowledge to construct efficient low dispersion schemes on unstructured grids ? Bonneton et al J.Geophys.Res. 2015) due to (G.B. Airy, Encyclopaedia Metropolitana, 1841) η = A cos(κx -ωt) travelling wave Φ = B(z) A sin(κx -ωt) potential κ wave number C = ω/κ phase velocity x (d 0 + z)) cosh(κ(d 0 ) Phase relation and phase dispersion C 2 0 = g d 0 shallow water (hyperbolic) celerity. C = C(κ)!! a wave amplitude λ = 2π/κ wave length T wave period λ = C(κ)T = C(2π/λ)T phase relation water equations (zero-th order model) ht+qx = 0 qt+ (ũq)x + hηx = 0 With the notation dimensionless depth : h = d + η dimensionless volume flux : q = hũ dimensionless depth averaged velocity: ũ Red : hyperbolic shallow water equations Blue terms are responsible for dispersion Nonlinear -non dispersive Example: Peregrine's equations 2 ht+qx = 0 qt+ (ũq)x + hηx = µ 2 With the notation dimensionless depth : h = d + η dimensionless volume flux : q = hũ dimensionless depth averaged velocity: ũ Red : hyperbolic shallow water equations Blue terms are responsible for dispersion Weakly nonlinear -weakly dispersive 2 D.H. Peregnine. J.Fluid.Mech, 1967 ht+qx = 0 qt+ (ũq)x + hηx = µ 2 β d2 qtxx + d dx 3 qtx + B d3 ηxxx + 2B d2 dx ηxx +O(µ 4 , µ 2 ) With the notation dimensionless depth : h = d + η dimensionless volume flux : q = hũ dimensionless depth averaged velocity: ũ Red : hyperbolic shallow water equations Blue terms are responsible for dispersion Weakly nonlinear -weakly dispersive Phase enhancement vie the tunable coeff. B and β = B + 1/3 (cf. later) 3 P.A. Madsen and O.R. Sorensen Coast.Eng., 1992 Dispersive surface waves: Boussinesq approach (8/8) Example: enhanced Serre-Green-Naghdi equations 4 ht+qx = 0 qt+ (ũq)x + hηx = µ 2 h ψ +O(µ 4 ) ψ =α ∂x( h2 ∂x(ũt + ũũx)) + (α -1) ∂x( h2 ∂xx η) + Q ψ (ũ, h, d; ũx, hx, dx) With the notation dimensionless depth : h = d + η dimensionless volume flux : q = hũ dimensionless depth averaged velocity: ũ Red : hyperbolic shallow water equations Blue terms are responsible for dispersion Fully nonlinear -weakly dispersive Phase enhancement vie the tunable coeff. α (cf. later)Shoaling test, weakly nonlinear models 5 variations for a given asymptotic accuracy (same linear limit), e.g. µ 2 d ≈ µ 2h as µ 2 η is negligible (weakly nonlinear) µ 2 ( dũ)xxt ≈ µ 2 qxxt as µ 2 (η ũ)xx is negligible (weakly nonlinear) etc. 5 S.T. Grilli et al J.Waterw.Port.C.-ASCE, 1994 -A.G. Filippini et al. Coast.Eng., 2015 Shoaling test, fully nonlinear models x

  dimensionless depth : h = d + η dimensionless volume flux : q = hũ with ũ the (dimensionless) depth averaged velocity Red terms provide the hyperbolic shallow water equations Blue terms are responsible for dispersion Fully nonlinear -weakly dispersive Phase enhancement vie the tunable coeff. α (cf. later)

Node -centered FV scheme ∂Ui ∂t + 1

 1 |Ci| ∂Ω Fn x + Gn y = 1 |Ci| Ω S b + hΨ ; Roe's Riemann solver + Harten-Hyman entropy fix High order reconstruction: weighted least squares 11 (quadratic or cubic) ; Well-balanced treatment of topography, wet/dry fronts, etc 12 11 Ollivier-Gooch et al. AIAA J. 2009; Wang et al JCP 2017 10 Bermudez&Vazquez, CAF 1994; Hubbard&Garcia-Navarro, JCP 2000; Brufau et al. IJNMF 2002; Castro, Math.&Computer Mod. 2005; etc.etc. Discretization of the (breaking) eGN on unstructured grids Dispersion analysis of the scheme Benchmarking: propagation test 1 Test Description: a0 = 0.2 [m], h0 = 1 [m] : a/h 0 = 0.0515, T = 1 [s] ; Energy transfer to higher harmonics; Experiments: (Berkhoff et al., 1982) Adapted mesh 88760 nodes ;Benchmarking: overtopping on a three dimensional reefBenchmarking: overtopping on a three dimensional reef t (g/h)

A.E. Green J.Fluid Mech., 1976 -F. Chazel et al. J.Sci.Comp., 2011

Alvarez-Samaniego and Lannes, Indiana Univ. Math J., 2008

shallow water solver + non-hydrostatic term ψ → (h n+1 , q n+1 ) .

P. Bonneton Ocean Eng., 2007

Example: enhanced linearized Madsen-Sorensen

Truncation error heuristics: dispersive models Example: enhanced linearized Madsen-Sorensen

The eGN equations

We recast the system in two independent steps:

This reformulation aims at exploiting the self-adjoint character of T

Linear system : The short of it ...

Interaction modelling/discretization error for near shore Boussinesq

Frst (hyperbolic) order system vs second order elliptic operator :

• High order on hyperbolic : third order (at least) for good dispersion • Elliptic component: can be treated with a second order method • Unstructured grid eGN: high(er) order FV + P1 FEM for elliptic part Wave breaking: revert to SW + shock capturing


G Bouselmi
email: bousselm@loria.fr

D Fohr
email: fohr@loria.fr

I Illina
email: illina@loria.fr

J.-P Haton

DISCRIMINATIVE PHONEME SEQUENCES EXTRACTION FOR NON-NATIVE SPEAKER'S ORIGIN CLASSIFICATION

In this paper we present an automated method for the classification of the origin of non-native speakers. The origin of non-native speakers could be identified by a human listener based on the detection of typical pronunciations for each nationality. Thus we suppose the existence of several phoneme sequences that might allow the classification of the origin of non-native speakers. Our new method is based on the extraction of discriminative sequences of phonemes from a non-native English speech database. These sequences are used to construct a probabilistic classifier for the speakers' origin. The existence of discriminative phone sequences in non-native speech is a significant result of this work. The system that we have developed achieved a significant correct classification rate of 96.3% and a significant error reduction compared to some other tested techniques.

INTRODUCTION

The problem of non-native speaker origin classification consists in detecting the mother tong of speakers uttering non-native speech. For example, the detection of the nationality of Spanish or French people uttering English words. This is different from simple origin detection as in the latter case, the decision is taken over native speech (ex. English people uttering English speech vs. French people uttering French speech). The issue we target here is closer to regional accent detection for the same native language.

With the recent advances in the speech recognition field, the automatic speech recognition (ASR) is more and more used especially in call centers. ASR is beneficial for both the callers and the call center companies as it allows an automated and fast processing with natural language and allows the reduction of the number of human operators in the repetitive task of phone replying. Let's consider the application of a car renting call center based on automatic speech recognition. In such case, the ASR system will interact with the customer and collect the information of his order such as the car type, the duration of the rent, the pick point etc. In this case, if the origin of non-native speakers is known, an adapted ASR system can be used Thanks to HIWIRE project for funding. in order to have better recognition accuracy. For the latter application as for a plane ticketing call center, there is a high probability of encountering non-native speakers.

The work presented here is part of the European project HIWIRE : Human Input that Works in Real Environments. It aims at the developments of means for helping human operators performing their duties in real environment conditions. The HIWIRE project consists in developing an automated system based on ASR that assists aircraft pilots in their tasks and communications. As communications between pilots and control operators must be in English, the system that is under construction within HI-WIRE will inherently be confronted with non-native English speech. Traditional ASR systems would be inefficient in such case as their performance drops drastically when confronted with non-native speech. This performance drop is a well known problem (see [START_REF] Bouselmi | Multilingual Non-Native Speech Recognition using Phonetic Confusion-Based Acoustic Model Modification and Graphemic Constraints[END_REF]).

Recent research works for non-native speech have already allowed a significant improvement in that filed. The approaches described in [START_REF] Bouselmi | Multilingual Non-Native Speech Recognition using Phonetic Confusion-Based Acoustic Model Modification and Graphemic Constraints[END_REF], [START_REF] Schaden | Generating Non-Native Pronunciation Lexicons by Phonological Rule[END_REF] and [START_REF] Morgan | Making a Speech Recognizer Tolerate Non-Native Speech Through Gaussian Mixture Merging[END_REF] allowed significative performance enhancement against non-native speech. Nevertheless, those approaches require the knowledge of the origin of the speakers uttering the speech they are applied to. Indeed, the modifications applied to the ASR system depend on both the native language and the spoken language.

A foreign accent classification procedure could be a great asset to any system based on speech recognition and confronted with non-native speakers. Only few articles have been published concerning non-native accent classification. For that matter, the approach developed by C. Teixeira et al. [START_REF] Teixeira | Accent identification[END_REF] was based on HMM phone models. It achieved 65.5% classification rate on isolated words database of Danish, German, British, Spanish and Italian speech. Arslan et al. [START_REF] Arslan | Language accent classification in american english[END_REF] used HMM phone models and HMM word models to identify Neutral, Chinese, Turkish and German accents. The method achieved a 68.3% classification rate on 5 isolated words. The approach of P. Angkitirakul et al. [START_REF] Angkitirakul | Advances in phone-based modeling for automatic accent classification[END_REF], based on Stochastic Trajectory Models (STM) and Parametric Trajectory Models (PTM), performed by 40.6% classification rate in supervised mode.

In the next section, we will describe the extraction of discriminative phone sequences for each foreign language and will sketch the decision process based on conditional probabilities. In section 3, we will describe the tests that we have carried and discuss their results. We will also discuss the future research work. Finally, we will end with a brief conclusion.

FOREIGN ACCENT DETECTION

It is well known that non-native speakers might produce pronunciation errors when uttering foreign speech ([1], [START_REF] Schaden | Generating Non-Native Pronunciation Lexicons by Phonological Rule[END_REF], [START_REF] Morgan | Making a Speech Recognizer Tolerate Non-Native Speech Through Gaussian Mixture Merging[END_REF], [START_REF] Ladefoged | The Sounds of the World's Languages[END_REF], [START_REF] Jeffers | Principles and Methods for Historical Linguistics[END_REF]). These errors are due to the phonological and articulatory properties of both the spoken and native languages.

For instance, some phonemes of the spoken language (SL) might not exist in the native language (NL) of the speaker. The speaker may replace these phonemes by some acoustically close phonemes of his NL. For instance, in French, diphthongs do not exist and some French speakers pronounce instead a sequence of two French phones.

These non-native pronunciation error depends on the pair of spoken and native languages. P. Ladefoged et al. [START_REF] Ladefoged | The Sounds of the World's Languages[END_REF] and R.J. Jeffers et al. [START_REF] Jeffers | Principles and Methods for Historical Linguistics[END_REF] have depicted in their work a list of common phone uttering errors made by various groups of non-native speakers for the English language (French, Italian, Greek etc.). Indeed, speakers from the same origin are very likely to commit the same pronunciation errors as they share the same native language, and thus the same articulatory and phonological mechanisms. Human listeners rely on those common errors as hints and tips to decide on the origin of non-native speakers. In the work of Arslan et al. [START_REF] Arslan | Language accent classification in american english[END_REF] and Angkitirakul et al. [START_REF] Teixeira | Accent identification[END_REF], human listener achieved 54% and 84% in foreign speech classification.

Our approach described here is based on that feature. We suppose the existence of discriminative uttering structures at the phonetic level that are shared among speakers from a particular origin when they speak a foreign language. In other words, we suppose that speakers from a particular origin X utter some discriminative sequences of phonemes when they speak a foreign language Y . We suppose that for a set of origins L = {L 1 ..L n } and a foreign language F , there exist sets of phoneme sequences S 1 ..S n corresponding to the origins L 1 ..L n (respectively) that might discriminate the native languages of L 1 ..L n speakers when they utter F speech (S i = {s i,1 ..s i,ki }, and s i,m are sequences of phones).

Discriminative phoneme sequences extraction

In order to better model the non-native speech, we have chosen to use the phone acoustic models (HMM) of all the native languages L 1 ..L n (the models are noted M 1 ..M n respectively). Let the non-native database B = n i=1 B i , where B i is the part of the database composed of L i speakers uttering F speech. First, the native phone models are adapted on the respective non-native database: i.e. the models M 1 are adapted on the B 1 resulting in M ′ 1 , and so on. Then, to extract these discriminative phone sequences, we perform a phonetic recognition with a phonetic pool M = n m=1 M ′ m on each of the non-native databases. For each native language L i , we count the occurrences of all the phone sequences having a maximum length of max p phones in the phonetic recognition results. The number of occurrences of the phone sequences is normalized against the number of sentences that compose each non-native parts of the database (B 1 ..B n ). This processing results in sets of preliminary phone sequences with their normalized number of occurrences for each language L i . Those sets are noted S ′ i = {s i,1 ..s i,k ′ i }, and the number of occurrences is noted n i (s) for a sequence s (for a language L i ). The next step consists in retaining the sets of discriminative sequences S i . For L i , a sequence s ∈ S ′ i is considered discriminative only if it verifies the equation 1.

n i (s) >= α * n k (s), ∀k = i ( 1
)
where α is a discriminant factor, α ≥ 1.

Knowing the S i sets and the counts of appearances of each of their sequences, some probabilities can be computed. All the following probabilities are conditional probabilities conditioned by the acoustic models M and the sets S 1 ..S n . For readability reasons, we will omit these conditions in the probabilities notations. The maximum likelihood (ML) probability P (L i ), P (s) and P (s/L i ) are computed as follows :

P (L i ) = Σ ki m=1 n i (s i,m ) Σ n l=1 Σ k l m=1 n l (s l,m ) = Σ x∈Si n i (x) Σ n l=1 Σ x∈S l n l (x) (2)
P (s) = Σ n l=1 n i (s) Σ n l=1 Σ k l m=1 n l (s l,m ) = Σ n l=1 n i (s) Σ n l=1 Σ x∈S l n l (x) (3)
P (s/L i ) = n i (s) Σ ki m=1 n i (s i,m ) = n i (s) Σ x∈Si n i (x) (4)
Using the bayes rule and the equations 2, 3 and 4, the conditional probability of a language L i knowing a sequence s can be computed as in the equation 5.

P (L i /s) = P (s/L i ) * P (L i ) P (s) = ni(s) Σx∈S i ni(x) * Σx∈S i ni(x) Σ n l=1 Σx∈S l n l (x) Σ n l=1 ni(s) Σ n l=1 Σx∈S l n l (x) = n i (s) Σ n l=1 n l (s) (5)
The conditional probability of a language L i knowing a list of sequences O = {s 1 ..s h } can be computed as in equation 6 using the bayes rule, the equations above and the hypothesis that sequences of O are independent. The hypothesis of independence of the sequences is not true. Nevertheless, this hypothesis must be assumed in order to compute this probability. Indeed, determining the interrelations between sequences of phones might prove to be impossible to compute with regards to the small size of our database.

P (L i /O) = P (O/L i )P (L i ) P (O) = P (s 1 ..s h /L i )P (L i ) P (s 1 ..s h ) = P (L i ) h m=1 P (s m /L i ) h m=1 P (s m ) = P (L i ) 1-h h m=1 P (s m /L i )P (L i ) P (s m ) = P (L i ) 1-h h m=1 P (L i /s m ) (6)

Classification of a speaker

In order to detect the origin of a speaker X, some of its recorded utterances must be analyzed. First, a phonetic recognition is performed on those sentences using a the models M described in 2.1. All the sequence of phones that appear in the sets S 1 ..S n are retained in a list O = {s 1 ..s h }. The speaker X is classified in the L i native language as in equation 7.

L i = argmax l=1..n {P (L l /O)} (7)
Another local decision approach can be adopted. Instead of collecting all the sequences of phones from the phonetic recognition (see last paragraph) that appear in S 1 ..S n in a single list O, separate lists O 1 ..O n corresponding to the sequences that appear in S 1 ..S n (respectively) could be made up. I.e., the list O l corresponds to all the sequences of phones observed in the phonetic recognition and that appear in the set S l (l = 1..n). The decision is then made over the probabilities of the languages knowing the lists O 1 ..O n , i.e. P (L l /O l ), ∀l = 1..n. In this classification approach, the probabilities of each language must be normalized over the number of sequences of each list O 1 ..O n in order to allow the comparison between them. Besides, any language decider that has a too small corresponding list of observations should be ignored, i.e., for a list O i , if exists k verifying card(O k ) ≥ β card(O i ), the classifier L i is ignored (β is a factor). If we note I the set of language indices that are not ignored, the speaker X is classified in the language L i if the equation 8 is verified.

L i = argmax l∈I {P (L l /O l ) 1 card(O l ) } (8)
3. EXPERIMENTS

Experimental conditions

Our tests have been carried out on the HIWIRE non-native speech database. This database have been tested in the approaches presented in [START_REF] Bouselmi | Multilingual Non-Native Speech Recognition using Phonetic Confusion-Based Acoustic Model Modification and Graphemic Constraints[END_REF] and [START_REF] Bouselmi | Fully Automated Non-Native Speech Recognition Using Confusion-Based Acoustic Model Integration and Graphemic Constraints[END_REF]. It is composed of 81 speakers: 31 French, 20 Greek, 20 Italian and 10 Spanish speakers. Each of those speakers reads 100 English sentences. The used grammar is a strict command language composed of 134 words. This grammar is used by aircraft pilots when communicating with airport control agents.

The speech was recorded in 16 bits and 16 kHz format.

We chose an MFCC parametrization with 13 coefficients and their first and second time derivatives. The acoustic models are 3 states HMMs (Hidden Markov Models) with 128 Gaussian mixtures and diagonal covariance matrices.

Tests and results

In ours tests, we have used 39 French, 33 Greek, 32 Spanish and 49 Italian monophone HMM models trained on native speech databases (respectively). As described in section 2.1, those models were adapted on HIWIRE nonnative database. I.e., the French models were adapted on all the French speakers, etc.

The extraction of the phone sequences was done following the "leave one out" scheme. For instance, when testing a French speaker X, the discriminative phone sequences of the French language are extracted using all the French speakers except X. And in that example, the significant sequences of the other languages are extracted using all the respective speakers.

In our preliminary tests, we have chose some threshold values as follows :

α = 4 : the significance factor (see section 2.1).

max p = 3 : the maximum length of a phone sequence (number of phones in the sequence).

-50 as the minimum occurrences count per speaker for a sequence to be eligible as discriminative.

-30 as the maximum discriminative sequences count per language.

β = 2.5 : see section 2.2.

Table 1 shows the confusion matrix in terms of speakers percentage for the global decision matrix. The accuracy achieved is 96.29% of correct speaker classification. Table 2 shows the confusion matrix in terms of speakers percentage for the local decision matrix. The accuracy achieved is 87.65% of correct speaker classification.

Discussion and future work

The tests presented here are only preliminary and we intend to further investigate and tune the parameters of this method. The effect of some parameters like the significance factor α, the maximum number of sequences per language and the minimum appearance count for sequences will be investigated in our future work. Besides, we will test the use of English acoustic models adapted on non-native speech instead of native models. This might avoid the inconvenience of collecting all the acoustic models of all non-native language that will be classified.

The potential of native language classification based on discriminative phone sequence might be great. Indeed, we have tested other mother tong detection (over nonnative speech) techniques inspired from the state of the art. We have tested a global GMM classification where a GMM was trained for each native language. We have also tested three HMM based approaches using TIMIT context independent phonemes. We have adapted the TIMIT phone models on the French, Greek, Spanish and Italian databases in a supervised fashion. The best result obtained with those methods is only 84% on the same HIWIRE database we have used. The global decision approach described above achieved a significantly better result of 96.3%, giving an error reduction of 76.9% (relative).

Another significant result of our work is the existence of discriminative phone sequences -or syllabic realizationsin non-native speech. Those phone sequences can be relied on to classify the origin of non-native speakers.

CONCLUSION

In this paper, we presented a novel approach for the detection of the mother tong of non-native speakers based on discriminative phone sequences. We have determined that there exists some discriminative phone sequences in non-native speech that could help in the mother tong detection. The preliminary results we obtained show a great potential for this technique: 96.3% correct classification rate. Our method will be further tested and tuned for the best classification results.

Table 1 .

 1 Confusion matrix for the global decision method. The classification rate is 96.29%.

		French Greek Italian Spanish
	French	100.0	0.0	0.0	0.0
	Greek	5.0	95.0	0.0	0.0
	Italian	0.0	5.0	95.0	0.0
	Spanish	0.0	0.0	10.0	90.0

Table 2 .

 2 Confusion matrix for the local decision method.

	The classification rate is 87.65%.
		French Greek Italian Spanish
	French	77.1	0.0	0.0	12.9
	Greek	5.0	85.0	5.0	5.0
	Italian	5.0	5.0	85.0	5.0
	Spanish	0.0	0.0	0.0	100.0

ACKNOWLEDGMENTS

This work was partially funded by the European project HIWIRE (Human Input that Works In Real Environments), contract number 507943, sixth framework program, information society technologies.