Alexandru Amȃrioarei
email: alexandru.amarioarei@inria.fr

Cristian Preda
email: cristian.preda@polytech-lille.fr

Approximations for two-dimensional discrete scan statistics in some block-factor type dependent models

Approximations for two-dimensional discrete scan statistics in some block-factor type dependent models

Alexandru Amarioarei, Cristian Preda

Introduction

Let N 1 , N 2 be positive integers, R = [0, N 1 ] × [0, N 2 ] be a rectangular region and {X i,j | 1 ≤ i ≤ N 1 , 1 ≤ j ≤ N 2 } be a family of random variables from a specified distribution. When the random variables X i,j take nonnegative integer values it is common to interpret them as the number of occurrences of some events observed in the elementary square sub-region

r i,j = [i -1, i] × [j -1, j]. Let m 1 , m 2 be positive integers such that 1 ≤ m 1 ≤ N 1 , 1 ≤ m 2 ≤ N 2 . For 1 ≤ i 1 ≤ N 1 -m 1 + 1, 1 ≤ i 2 ≤ N 2 -m 2 + 1 define Y i1,i2 = Y i1,i2 (m 1 , m 2 ) = i1+m1-1 i=i1 i2+m2-1 j=i2 X i,j
(1.1)

as the number of events in the rectangular region R(i

1 , i 2 ) = [i 1 -1, i 1 + m 1 -1] × [i 2 -1, i 2 + m 2 -1],
comprised of m 1 × m 2 adjacent elementary squares r i,j . The two-dimensional discrete scan statistic is defined as the largest number of events in any rectangular scanning window R(i 1 , i 2 ), within the rectangular region R, i.e.

S = S m1,m2 (N 1 , N 2 ) = max 1≤i1≤N1-m1+1 1≤i2≤N2-m2+1 Y i1,i2 . (1.2)
Most of research devoted to the two-dimensional discrete scan statistic considers the i.i.d. model for the random variables X i,j . Then, the statistic S is used for testing the null hypothesis of randomness (H 0 ), that assumes that X i,j 's are independent and identically distributed according to some specified probability law, in general Bernoulli, binomial or Poisson [START_REF] Chen | Two-dimensional discrete scan statistics[END_REF], [START_REF] Glaz | Scan statistic[END_REF]), against an alternative (H 1 ) of clustering. Under H 1 , one suppose that there is a change, with respect to H 0 , in the distribution of the random field within a rectangular sub-region R(i * , j * ) ⊂ R, with 1 ≤ i * ≤ N 1m 1 + 1, 1 ≤ j * ≤ N 2m 2 + 1, while outside this region X i,j 's are distributed according to the null hypothesis distribution. As an example, consider that under H 0 , X i,j 's are i.i.d.

Poisson random variables with mean λ 0 . In this setting, the alternative hypothesis assumes that there exists a rectangular sub-region R(i * , j * ) such that for i * ≤ i ≤ i * + m 1 -1 and j * ≤ j ≤ j * + m 2 -1 the distribution of X i,j is given by a Poisson distribution of mean λ 1 > λ 0 whereas in R \ R(i * , j * ) the events occur according to the distribution specified by the null hypothesis.

The distribution of the two-dimensional scan statistic,

P (S m1,m2 (N 1 , N 2 ) ≤ n)
is successfully applied in brain imaging [START_REF] Naiman | Computing Scan Statistic p Values Using Importance Sampling, with Applications to Genetics and Medical Image Analysis[END_REF]), astronomy [START_REF] Darling | Extreme value distribution for the largest cube in a random lattice[END_REF], [START_REF] Marcos | From star complexes to the field: open cluster families[END_REF]), target detection in sensors networks (Goerriero, Willett and Glaz [2009]), reliability theory [START_REF] Boutsikas | Reliability approximations for Markov chain imbeddable systems[END_REF]) among many other domains. For an overview of the potential applications of the scan statistics one can refer to the monographs of [START_REF] Glaz | Scan statistic[END_REF] and more recently the one of Glaz, Pozdnyakov and Wallenstein [2009, Chapter 6].

Since there are no exact formulas for P(S ≤ n), various methods of approximation and bounds have been proposed by several authors. An overview of these methods as well as a complete bibliography on the subject can be found in [START_REF] Chen | Two-dimensional discrete scan statistics[END_REF], Glaz, Naus and Wallenstein [2001, Chapter 16], [START_REF] Boutsikas | Bounds for the distribution of two dimensional binary scan statistics[END_REF], [START_REF] Haiman | Estimation for the distribution of two-dimensional scan statistics[END_REF] and the references therein.

In this paper we introduce a dependence structure for the underlying random field

({X i,j | 1 ≤ i ≤ N 1 , 1 ≤ j ≤ N 2 })
based on a block-factor model and approximate the distribution of the two dimensional discrete scan statistics in this setting. Writing the scan statistics random variable S as the maximum of a 1-dependent stationary sequence, we approximate its distribution employing a result obtained by [START_REF] Haiman | First passage time for some stationary processes[END_REF] and later improved by Amȃrioarei [2012]. This approach was successfully used to evaluate the distribution of scan statistics, both in discrete and continuous cases, in a series of articles: for one-dimensional case in [START_REF] Haiman | Estimating the distributions of scan statistics with high precision[END_REF] and [START_REF] Haiman | Estimating the distribution of one-dimensional discrete scan statistics viewed as extremes of 1-dependent stationary sequences[END_REF], for two-dimensional case in [START_REF] Haiman | Estimation for the distribution of two-dimensional scan statistics[END_REF] and [START_REF] Haiman | A new method for estimating the distribution of scan statistics for a two dimensional Poisson process[END_REF] and for three-dimensional case in Amȃrioarei and Preda [2013]. The advantage of our approach is that it can be applied under very general conditions and provides accurate approximations and sharp bounds for the errors. The paper is organized as follows. In Section 2, we introduce the block-factor type model that will generate the random field to be scanned. The methodology for approximating the distribution of the scan statistics generated by the block-factor model as well as the associated error bounds are presented in Section 3. Section 4 includes numerical results based on simulations for a particular block-factor model.

Block-factor type model

In this section we introduce a particular dependence structure for the random field

{X i,j | 1 ≤ i ≤ N 1 , 1 ≤ j ≤ N 2 } based on a block-factor type model.
Recall that (see [START_REF] Burton | On one-dependent processes and k-block factors[END_REF]) a sequence (W l ) l≥1 of random variables with state space S W is said to be a k block-factor of the sequence ( Wl ) l≥1 with state space S W , if there is a measurable function f :

S k W → S W such that W l = f Wl , Wl+1 , . . . , Wl+k-1
for all l.

Our block-factor type model is defined in the following way. Let Ñ1 , Ñ2 be positive integers and Xi,j | 1 ≤ i ≤ Ñ1 , 1 ≤ j ≤ Ñ2 be a family of independent and identically distributed real valued random variables. Notice that if the region R = [0, Ñ1 ] × [0, Ñ2 ] is divided in a grid with step 1, then we can locate the random variables Xi,j as being at the intersection of the j-th row with the i-th column. Let x 1 , x 2 , y 1 , y 2 be nonnegative integers such that x 1 + x 2 ≤ Ñ1 -1 and y 1 + y 2 ≤ Ñ2 -1. Define c 1 = x 1 + x 2 + 1, c 2 = y 1 + y 2 + 1 and take N s = Ñsc s + 1 for s ∈ {1, 2}. To each pair (i, j) ∈ {x 1 + 1, . . . , Ñ1x 2 } × {y 1 + 1, . . . , Ñ2y 2 } we associate the random matrix of size c 2 × c 1 , C (i,j) ∈ M c2,c1 (R), with entries

C (i,j) (k, l) = Xi-x1-1+l,j+y2+1-k , 1 ≤ k ≤ c 2 , 1 ≤ l ≤ c 1 . (2.1)
If T : M c2,c1 (R) → R is a measurable function then the block-factor type model is given by

X i,j = T C (i+x1,j+y1) with 1 ≤ i ≤ N 1 , 1 ≤ j ≤ N 2 . ( 2
T X i-x 1 ,j-y 1 (a) 1 . . . x 1 + 1 Ñ1 -x 2 . . . Ñ1 1 . . . y 1 + 1 Ñ2 -y 2 . . . Ñ2 T 1 N1 1 N2 Xi,j Xi-x 1 ,j-y 1 (b) Figure 1. Illustration of the block-factor type model Obviously, {X i,j | 1 ≤ i ≤ N 1 , 1 ≤ j ≤ N 2 } forms a dependent family of random variables (see Fig 2).
Recall that a sequence (W k ) k≥1 is m-dependent with m ≥ 1 (see [START_REF] Burton | On one-dependent processes and k-block factors[END_REF]), if for any h ≥ 1 the σ-fields generated by {W 1 , . . . , W h } and {W h+m+1 , . . . } are independent. From the definition of the random variables X i,j given by Eq.(2.2), we observe that for each 1 ≤ i ≤ N 1 the sequence (X i,j ) 1≤j≤N2 is (c 2 -1)-dependent and for each 1

≤ j ≤ N 2 the sequence (X i,j ) 1≤i≤N1 is (c 1 -1)- dependent (see also Fig 2).
Remark 2.1. Notice that if c 1 = c 2 = 1 (x 1 = x 2 = 0 and y 1 = y 2 = 0) then the sequence X i,j = Xi,j and we are in the i.i.d. situation. In this case the distribution of the two-dimensional scan statistics can be approximated using the known methods

Xi,j Xi-c 1 +1,j-c 2 +1 Xi-c 1 +1,j+c 2 -1 Xi+c 1 -1,j-c 2 +1 Xi+c 1 -1,j+c 2 -1 c2 -1 c2 -1 c1 -1 c1 -1 Figure 2. The dependence of X i,j
(see Glaz, Naus and Wallenstein [2001, Chapter 16] and [START_REF] Glaz | Scan statistic: Methods and Applications[END_REF]).

If we take Ñ2 = 1, which automatically implies that c 2 = 1, we obtain an one dimensional block-factor model W i = X i,1 and the two-dimensional scan statistic becomes the usual discrete scan statistics in one dimension over a (c 1 -1)-dependent sequence. The distribution of one dimensional scan statistics over this type of dependence was studied by [START_REF] Haiman | One dimensional scan statistics generated by some dependent stationary sequences[END_REF] in the particular case of Gaussian stationary 1-dependent (x 1 = 0, x 2 = 1 and c 1 = 2) sequences W i ∼ N (0, 1) of random variables generated by a two block-factor of the form

W i = aU i + bU i+1 , i ≥ 1,
where a 2 + b 2 = 1 and (U i ) i≥1 is an i.i.d. sequence of N (0, 1) random variables.

An application of the one dimensional scan statistics over a sequence of moving average of order q (c 1 = q + 1) is presented in Section 4.2.

Based on the model presented in this section, in Section 3 we give an approximation for the distribution of two-dimensional scan statistic over the random field generated by the family X i,j and the corresponding error bounds.

Approximation and error bounds

In this section we present the methodology used to obtain the approximation of the two-dimensional discrete scan statistics distribution over the field generated by the block-factor model described in Section 2. Let's consider the scanning window of size m 1 × m 2 with m 1 ≥ 2, m 2 ≥ 2 and assume that for s ∈ {1, 2}, Ñs = (L s + 1)(m s + c s -2) where L 1 , L 2 are positive integers. Observe that

N s = L s (m s + c s -2) + m s -1, s ∈ {1, 2}
and define the sequence

Z k = max 1≤i1≤L1(m1+c1-2) (k-1)(m2+c2-2)+1≤i2≤k(m2+c2-2) Y i1,i2 , k ∈ {1, 2, . . . , L 2 }.
(3.1)

The random variables Z k represent the scan statistics on the overlapping

N 1 × 2(m 2 + c 2 -2) -(c 2 -1) rectangular regions R k = [1, N 1 ] × [(k -1)(m 2 + c 2 -2) + 1, (k + 1)(m 2 + c 2 -2) -(c 2 -1)]. m 2 m 1 Z 1 Z 2 Z 3 c 2 m 2 + c 2 -1 2m 2 + c 2 -3 2m 2 + 2c 2 -3 3m 2 + 2c 2 -5 4m 2 + 3c 2 -7 N 2 N 1 1 Figure 3. Illustration of Z k emphasizing the 1-dependence Remark 3.1. If we consider the extreme situation when c 2 = m 2 = 1 (or c 1 = m 1 = 1)
, that is when we have row (column) independence and we are scanning only on rows (columns), then the sequence described by Eq.(3.1) is no longer well defined. In this case we define

Z k = max 1≤i1≤L1(m1+c1-2) Y i1,k , k ∈ {1, 2, . . . , L 2 }, (3.2) where Y i1,k = i1+m1-1 i=i1 X i,k .
We observe that from Eq.(3.1) the set of random variables {Z 1 , . . . , Z L2 } is 1dependent (see also Figure 3). Indeed, we have

Z k-1 ∈ σ {X i,j |1 ≤ i ≤ N 1 , (k -2)(m 2 + c 2 -2) + 1 ≤ j ≤ k(m 2 + c 2 -2) -(c 2 -1)} ∈ σ Xi,j |1 ≤ i ≤ Ñ1 , (k -2)(m 2 + c -2 -2) + 1 ≤ j ≤ k(m 2 + c 2 -2)
and similarly,

Z k ∈ σ Xi,j |1 ≤ i ≤ Ñ1 , (k -1)(m 2 + c 2 -2) + 1 ≤ j ≤ (k + 1)(m 2 + c 2 -2) , Z k+1 ∈ σ Xi,j |1 ≤ i ≤ Ñ1 , k(m 2 + c 2 -2) + 1 ≤ j ≤ (k + 2)(m 2 + c 2 -2) .
From the above relations, the measurability of T from the definition of the depen-

dent model {X i,j |1 ≤ i ≤ N 1 , 1 ≤ j ≤ N 2 } and the independence of the sequence Xi,j |1 ≤ i ≤ Ñ1 , 1 ≤ j ≤ Ñ2 we conclude that the sequence (Z k ) 1≤k≤L2 is 1- dependent.
Since Xi,j are identically distributed we deduce stationarity of the random variables Z k .

Notice that from Eq.(3.1) and the definition of the two-dimensional scan statistics in Eq.( 2.1) we have the following relation

S = max 1≤k≤L2 Z k . (3.3)
The relation described by Eq.(3.3) is the key idea behind our approximation, i.e. the scan statistic random variable can be expressed as a maximum of 1-dependent stationary sequence of random variables. The approximation methodology that we use is based on the following result developed in Haiman [1999, Theorem 4] and improved in Amȃrioarei [2012, Theorem 2.6]: Let (W k ) k≥1 be a stationary 1-dependent sequence of random variables and for x < sup{u|P(W 1 ≤ u) < 1}, consider

q m = q m (x) = P(max(W 1 , . . . , W m ) ≤ x). (3.4)
Theorem 3.2. Assume that x is such that q 1 (x) ≥ 1-α ≥ 0.9 and define η = 1+lα with l = l(α) > t 3 2 (α) and t 2 (α) the second root in magnitude of the equation αt 3t + 1 = 0. Then the following relation holds

q m - 2q 1 -q 2 [1 + q 1 -q 2 + 2(q 1 -q 2 ) 2 ] m ≤ mF (α, m)(1 -q 1 ) 2 , (3.5) with F (α, m) = 1 + 3 m + Γ(α) m + K(α) (1 -q 1 ), (3.6)
and where Γ(α) = L(α) + E(α) and

K(α) = 11-3α (1-α) 2 + 2l(1 + 3α) 2+3lα-α(2-lα)(1+lα) 2 [1-α(1+lα) 2 ] 3 1 -2α(1+lα) [1-α(1+lα) 2 ] 2 , (3.7) L(α) = 3K(α)(1 + α + 3α 2 )[1 + α + 3α 2 + K(α)α 3 ] + α 6 K 3 (α), + 9α(4 + 3α + 3α 2 ) + 55.1 (3.8) E(α) = η 5 [1 + (1 -2α)η] 4 [1 + α(η -2)] 1 + η + (1 -3α)η 2 2(1 -αη 2 ) 4 [(1 -αη 2 ) 2 -αη 2 (1 + η -2αη) 2 ] .
(3.9)

Following the approach in Amȃrioarei and Preda [2013] for three dimensional scan statistics, we obtain an approximation formula for the distribution of two-dimensional scan statistic S along with the corresponding error bounds, in two steps as follows.

Define for r ∈ {2, 3},

Q r = Q r (n) = P r-1 k=1 {Z k ≤ n} = P    max 1≤i1≤L1(m1+c1-2) 1≤i2≤(r-1)(m2+c2-2) Y i1,i2 ≤ n    . (3.10)
For n such that Q 2 (n) ≥ 1α 1 ≥ 0.9, we apply the result in Theorem 3.2 to obtain the first step approximation

P (S ≤ n) - 2Q 2 -Q 3 [1 + Q 2 -Q 3 + 2(Q 2 -Q 3 ) 2 ] L2 ≤ L 2 F (α 1 , L 2 )(1 -Q 2 ) 2 . (3.11)
In order to evaluate the approximation in Eq.(3.11) one has to find approximations for the quantities Q 2 and Q 3 . To achieve this, we apply again the result of Theorem 3.2. We define, as in Eq.(3.1), for each r ∈ {2, 3} and l ∈ {1, 2, . . . , L 1 } the random variables

Z (r) l = max (l-1)(m1+c1-2)+1≤i1≤l(m1+c1-2) 1≤i2≤(r-1)(m2+c2-2)
Y i1,i2 .

(3.12)

As described in the case of the sequence Z k , we deduce that the random variables Z (r) l defined by Eq.(3.12) are stationary, 1-dependent and the following relation holds:

Q r = P max 1≤l≤L1 Z (r) l ≤ n , r ∈ {2, 3}.
(3.13)

Denoting, for u, v ∈ {2, 3}

Q uv = Q uv (n) = P u-1 l=1 {Z (v) l ≤ n} = P    max 1≤i1≤(u-1)(m1+c1-2) 1≤i2≤(v-1)(m2+c2-2) Y i1,i2 ≤ n   
(3.14) then, under the supplementary condition that n is such that Q 23 (n) ≥ 1α 2 ≥ 0.9, we apply Theorem 3.2 to obtain

Q r - 2Q 2r -Q 3r [1 + Q 2r -Q 3r + 2(Q 2r -Q 3r ) 2 ] L1 ≤ L 1 F (α 2 , L 1 )(1 -Q 2r ) 2 . (3.15)
Combining Eq.(3.11) and Eq.(3.15) we find an approximation formula for the distribution of the two-dimensional scan statistic depending on the values of Q 22 , Q 23 , Q 32 and Q 33 . There are no exact formulas for Q uv , u, v ∈ {2, 3}, thus these quantities will be evaluated using Monte Carlo simulation. The approximation process is summarized by the diagram in Figure 4:

m 2 m 1 N 2 N 1 R 2m 2 + c 2 -3 m 1 N 2 N 1 Q2 2m 2 + c 2 -3 2m 1 + c 1 -3 N 2 N 1 Q22 2m 2 + c 2 -3 3m 1 + 2c 1 -5 N 2 N 1 Q32 3m 2 + 2c 2 -5 m 1 N 2 N 1 Q3 3m 2 + 2c 2 -5 2m 1 + c 1 -3 N 2 N 1 Q23 3m 2 + 2c 2 -5 3m 1 + 2c 1 -5 N 2 N 1 Q33 Figure
P (S m1,m2 (M 1 , M 2 ) ≤ n) ≤ P (S m1,m2 (N 1 , N 2 ) ≤ n) ≤ P (S m1,m2 (T 1 , T 2 ) ≤ n) ,
(3.16) where for j ∈ {1, 2} we consider M j = (L j + 2)(m j + c j -2) -(c j -1) and T j = (L j + 1)(m j + c j -2) -(c j -1), we can approximate the distribution of the scan statistics by linear interpolation.

3.1.

Computing the error bounds. For the error computation we have to notice that there are three expressions involved: the first one is the theoretical error (E app ) obtained from the substitution of Eq.(3.15) in Eq.(3.11) whereas the other two are simulations errors, one corresponding to the approximation formula (E sf ) and the other to the error formula (E sapp ). In what follows we will deal with each of them separately. To simplify the presentation it will be convenient to introduce the following notations:

H(x, y, m) = 2x -y [1 + x -y + 2(x -y) 2 ] m , α 1 = 1 -Q 3 , α 2 = 1 -Q 23 , F 1 = F (α 2 , L 1 ), F 2 = F (α 1 , L 2 ), R s = H (Q 2s , Q 3s , L 1 ) , s ∈ {2, 3}.
Notice that the choice for the thresholds α 1 and α 2 is natural since we have the inequalities

Q 3 ≤ Q 2 and Q 23 ≤ Q 22 .
Based on mean value theorem in two dimensions, one can easily verify that if y i ≤ x i , i ∈ {1, 2} then we have the relation

|H(x 1 , y 1 , m) -H(x 2 , y 2 , m)| ≤ m [|x 1 -x 2 | + |y 1 -y 2 |] .
(3.17)

Rewriting Eq.(3.11) using the above notations and applying the inequality in Eq.(3.17) we can write

|P(S ≤ n) -H (R 2 , R 3 , L 2 )| ≤ |P(S ≤ n) -H (Q 2 , Q 3 , L 2 )| + |H (Q 2 , Q 3 , L 2 ) -H (R 2 , R 3 , L 2 )| ≤ L 2 F 2 (1 -Q 2 ) 2 + L 2 [|Q 2 -R 2 | + |Q 3 -R 3 |] . (3.18)
If we substitute Eq.(3.15) in Eq.(3.18) and take

B 2 = 1 -R 2 + L 1 F 1 (1 -Q 22 ) 2 ,
then the theoretical approximation error is given by

E app = L 2 F 2 B 2 2 + L 1 L 2 F 1 (1 -Q 22 ) 2 + (1 -Q 23 ) 2 . (3.19)
To compute the simulation error corresponding to the approximation formula let us denote with Quv the simulated values corresponding to Q uv for each u, v ∈ {2, 3}.

Usually between the true and the estimated values we have a relation of the type

Q uv -Quv ≤ β uv . (3.20)
Indeed, if IT ER is the number of iterations used in the Monte Carlo simulation algorithm for the estimation of Q uv then, one can consider, for example, the bound

β uv = 1.96 Quv(1-Quv) IT ER
with a 95% confidence level. Taking for r ∈ {2, 3}, Qr = H Q2r , Q3r , L 1 to be the simulated values that corresponds to Q r and applying Eq.(3.17) whenever Q3 ≤ Q2 we get

H (R 2 , R 3 , L 2 ) -H Q2 , Q3 , L 2 ≤ L 2 R 2 -Q2 + R 3 -Q3 ≤ L 1 L 2 Q 22 -Q22 + Q 23 -Q23 + Q 32 -Q32 + Q 33 -Q33 . (3.21)
Combining Eq.(3.21) and Eq.(3.20) we obtain the simulation error associated with the approximation formula

E sf = L 1 L 2 (β 22 + β 23 + β 32 + β 33 ). (3.22)
Finally, introducing

C 2v = 1 -Q2v + β 2v , v ∈ {2, 3}, C 2 = 1 -Q2 + L 1 (β 22 + β 32 ) + L 1 F 1 C 2 22 ,
and substituting them in the theoretical approximation error formula in Eq.(3.19), we obtain the simulation error corresponding to the approximation error formula

E sapp = L 2 F 2 C 2 2 + L 1 L 2 F 1 C 2 22 + C 2 23 . ( 3

Examples and numerical results

In order to illustrate the efficiency of the approximation and the error bounds obtained in Section 3, we consider the following examples: a minesweeper game presented in Section 4.1 and an one dimensional scan statistics over a moving average model described in Section 4.2.

4.1. Example 1: minesweeper game. Let Ñ1 , Ñ2 be positive integers and Xi,j | 1 ≤ i ≤ Ñ1 , 1 ≤ j ≤ Ñ2 be a family of i.i.d. Bernoulli random variables of parameter p. We interpret the random variable Xi,j as representing the presence ( Xi,j = 1) or absence ( Xi,j = 0) of a mine in the elementary square region ri,j =

[i -1, i] × [j -1, j].
In this example we consider x 1 = x 2 = 1 and y 1 = y 2 = 1. Based on the notations introduced in Section 2, we observe that

c 1 = c 2 = 3, N 1 = Ñ1 -2 and N 2 = Ñ2 -2.
For each (i, j) ∈ {2, . . . , Ñ1 -1} × {2, . . . , Ñ2 -1} the configuration matrix is given by

C (i,j) = C (i,j) (k, l) 1≤k≤3 1≤l≤3
, where and define for 1 ≤ i ≤ N 1 and 1 ≤ j ≤ N 2 , the block-factor model

C (i,j) (k, l) = Xi+l-2,j+2-k . (4.1) Let T : M 3,3 (R) → R T  
X i,j = T C (i+1,j+1) = (s,t)∈{0,1,2} 2 (s,t) =(1,1) Xi+s,j+t . (4.3)
The random variable X i,j can be interpreted as the number of neighboring mines associated with the location (i, j). In Figure 5 we present a realization of the introduced model. On the left, we have the realization of the initial set of random variables where the gray squares represent the presence of mines while the white squares signifies the absence of mines. On the right side we have the realization of the X i,j random variables, that is the corresponding number of neighboring mines associated to each site. We present numerical results (Table 1-Table 8) for the described block-factor model with Ñ1 = Ñ2 = 44 (that is N 1 = N 2 = 42), m 1 = m 2 = 3 and the underlying random field generated by i.i.d. Bernoulli random variables of parameter p ( Xi,j ∼ B(p)) in the range {0.1, 0.3, 0.5, 0.7}. We also include numerical values for the corresponding i.i.d. model:

N 1 = N 2 = 42, m 1 = m 2 = 3 and X i,j ∼ B(8, p).
For all our results presented in the tables we used Monte Carlo simulations with 10 8 iterations for the block-factor model and with 10 5 replicas for the i.i.d. model. Notice that the contribution of the approximation error (E app ) to the total error is almost negligible in most of the cases with respect to the simulation error (E sim ). Thus, the precision of the method will depend mostly on the number of iterations Xi,j : X i,j : In this example we consider the particular situation of an one dimensional scan statistics over a M A(q) model.

1 2 . . . Ñ2 -1 Ñ2 1 2 . . . Ñ1 -1 Ñ1 1 2 . . . N2 1 2 . . . N1 2 2 1 2 2 2 3 2 4 5 5 3 2 2 2 2 2 2 3 2 1 1 2 1 3 4 4 4 3 1 2 1 1 2 2 1 2 0 1 0 3 2 1 1 2 1 1 0 2 0 0 0 0 0 1 1 4 4 1 2 1 1 1 0 T Figure 5. A realization of the minesweeper related model Table 1. Block-factor: m 1 = m 2 = 3, Ñ1 = Ñ2 = 44, N 1 = N 2 = 42, p = 0.1, IT ER = 10 8 n Sim Approx Eapp Esim E total Table 2. Independent: m 1 = m 2 = 3, N 1 = N 2 = 42, B(n = 8, p = 0.
m 1 = m 2 = 3, Ñ1 = Ñ2 = 44, N 1 = N 2 = 42, p = 0.3, IT ER = 10 8 n Sim Approx Eapp Esim E total Table 5. Block-factor: m 1 = m 2 = 3, Ñ1 = Ñ2 = 44, N 1 = N 2 = 42, p = 0.
In the two dimensional block-factor model introduced in Section 2 we consider Ñ2 = 1, which in particular implies that c 2 = 1 and m 2 = 1, x 1 = 0 and x 2 = q for q ≥ 1 a positive integer. Let m 1 ≥ 2, Ñ1 ≥ m 1 + q + 1 be positive integers and Xi = Xi,1 | 1 ≤ i ≤ Ñ1 be a sequence of i.i.d. Gaussian random variables with known mean µ and variance σ 2 . We observe that N 1 = Ñ1q and that for each i ∈ {1, . . . , N 1 } the configuration matrix becomes

C (i) = Xi , Xi+1 , . . . , Xi+q . (4.4)
Let the transformation T : M 1,q+1 (R) → R be defined by

T (x 1 , . . . , x q+1 ) = a 1 x 1 + a 2 x 2 + • • • + a q+1 x q+1 , (4.5)
where a = (a 1 , . . . , a q+1 ) ∈ R q+1 a not null vector and consider the block-factor model

X i = T C (i) = a 1 Xi + a 2 Xi+1 + • • • + a q+1 Xi+q , 1 ≤ i ≤ N 1 . (4.6)
Clearly, the sequence X 1 , . . . , X N1 forms a M A(q) model. Notice that the moving sums

Y t = Y t,1 , 1 ≤ t ≤ N 1 -m 1 + 1, can be expressed as Y t = t+m1-1 i=t X i = b 1 Xt + b 2 Xt+1 + • • • + b m1+q Xt+m1-1+q . (4.7)
If, for example, m 1 ≥ q then the coefficients b 1 , . . . , b m1+q are given by

b k =                        k j=1 a j , k ∈ {1, . . . , q -1} q+1 j=1 a j , k ∈ {q + 1, . . . , m 1 } k j=k-m1+1
a j , k ∈ {m 1 + 1, . . . , m 1 + q}.

(4.8)

Therefore, for each t ∈ {1, . . . , N 1m 1 + 1}, the random variable Y t follows a normal distribution with mean

E [Y t ] = (b 1 + • • • + b m+q )µ and variance V ar [Y t ] = b 2 1 + • • • + b 2 m+q σ 2 . The covariance matrix Σ = {Cov [Y t , Y s ]} has the entries Cov [Y t , Y s ] =          m1+q-|t-s| j=1 b j b |t-s|+j   σ 2 , |t -s| ≤ m 1 + q -1 0 , otherwise.
(4.9)

Given the mean and the covariance matrix of the vector (Y 1 , . . . , Y N1-m1+1 ), one can use the importance sampling algorithm developed by [START_REF] Naiman | Computing Scan Statistic p Values Using Importance Sampling, with Applications to Genetics and Medical Image Analysis[END_REF] (see also [START_REF] Malley | A Comprehensive Method for Genome Scans[END_REF] and Shi, Siegmund and Yakir [2001]) to estimate the distribution of the one dimensional scan statistics S = S m1 (N 1 ). Another way is to use the algorithm developed by [START_REF] Genz | Computation of Multivariate Normal and T Probabilities[END_REF] to approximate the multivariate normal distribution. In this paper we adopt the importance sampling procedure.

In order to evaluate the accuracy of the approximation developed in Section 3, we consider q = 2, N 1 = 1000, m 1 = 20, Xi ∼ N (0, 1) and the coefficients of the moving average model (a 1 , a 2 , a 3 ) = (0.3, 0.1, 0.5). In Table 9 we present numerical results for the setting described above. In our algorithms we used IT ER app = 10 6 iterations for the approximation and IT ER sim = 10 5 replicas for the simulation. In Figure 8 we illustrate the cumulative distribution functions obtained by approximation and simulation. For the approximation we present also the corresponding lower and upper bounds (computed from the total error of the approximation process (E total column in Table 9)). In this article we derived an approximation for the two dimensional discrete scan statistic generated by a block-factor type model obtained from an i.i.d. sequence. Our method provides a sharp approximation for the high order quantiles of the distribution of the scan statistics along with the corresponding error bounds. A simulation study was included to show the accuracy of our method.

  Figure 4. Illustration of the approximation process Remark 3.3. If Ñ1 and Ñ2 are not multiples of m 1 + c 1 -2 and m 2 + c 2 -2, respectively, then we take L j + 1 = Ñj mj +cj -2 for j ∈ {1, 2}. Based on the inequalities

Figure 6 .Figure 7 .

 67 Figure 6. Cumulative distribution function for block-factor and i.i.d. models

Figure 8 .

 8 Figure 8. Cumulative distribution function for approximation and simulation along with the corresponding error under M A model

Table 3 .

 3 1), IT ER = 10 5 Block-factor:

	n	Sim	Approx	Eapp	Esim	E total
	17 0.789376 0.788934 0.005813 0.011393 0.017206
	18 0.925456 0.925186 0.000529 0.002095 0.002625
	19 0.976889 0.976763 0.000045 0.000455 0.000500
	20 0.993444 0.993447 0.000003 0.000105 0.000108
	21 0.998288 0.998287 0.000000 0.000023 0.000024
	22 0.999584 0.999583 0.000000 0.000005 0.000005
	23 0.999905 0.999905 0.000000 0.000001 0.000001
	24 0.999980 0.999980 0.000000 0.000000 0.000000

Table 6 .

 6 5, IT ER = 10 8 Independent:m 1 = m 2 = 3, N 1 = N 2 = 42, B(n = 8, p = 0.5), IT ER = 10 5

	n	Sim	Approx	Eapp	Esim	E total
	61 0.725109 0.701781 0.085110 0.093544 0.178654
	62 0.828019 0.888902 0.004453 0.008665 0.013118
	63 0.899560 0.888902 0.004453 0.008665 0.013118
	64 0.945304 0.939436 0.001049 0.004054 0.005103
	65 0.972203 0.969026 0.000235 0.002334 0.002569
	66 0.986999 0.985439 0.000047 0.001460 0.001507
	67 0.994506 0.993814 0.000008 0.000927 0.000935
	68 0.997851 0.997605 0.000001 0.000572 0.000573
	69 0.999326 0.999230 0.000000 0.000320 0.000320
	70 0.999826 0.999786 0.000000 0.000171 0.000171
	71 0.999968 0.999952 0.000000 0.000083 0.000083
	72 1.000000 1.000000 0.000000 0.000000 0.000000
	n	Sim	Approx	Eapp	Esim	E total
	50 0.741089 0.735210 0.010514 0.018002 0.028516
	51 0.882209 0.880827 0.001499 0.004196 0.005695
	52 0.952545 0.952389 0.000200 0.001098 0.001299
	53 0.982842 0.982891 0.000024 0.000307 0.000331
	54 0.994328 0.994337 0.000002 0.000084 0.000087
	55 0.998282 0.998278 0.000000 0.000022 0.000022
	56 0.999517 0.999518 0.000000 0.000005 0.000005
	57 0.999876 0.999876 0.000000 0.000001 0.000001
	58 0.999971 0.999971 0.000000 0.000000 0.000000
	59 0.999994 0.999994 0.000000 0.000000 0.000000
	60 0.999999 0.999999 0.000000 0.000000 0.000000

Table 7 .

 7 Block-factor: m 1 = m 2 = 3, Ñ1 = Ñ2 = 44, N 1 = N 2 = 42, p = 0.7, IT ER = 10 8

	n	Sim	Approx	Eapp	Esim	E total
	70 0.729239 0.705944 0.074290 0.082392 0.156682
	71 0.876484 0.864370 0.006976 0.011623 0.018600
	72 1.000000 1.000000 0.000000 0.000000 0.000000

Table 8 .

 8 Independent: m 1 = m 2 = 3, N 1 = N 2 = 42, B(n = 8, p = 0.7), IT ER = 10 5

	n	Sim	Approx	Eapp	Esim	E total
	62.0 0.620295 0.611819 0.030328 0.042319 0.072646
	63.0 0.847421 0.846730 0.002591 0.005851 0.008442
	64.0 0.952524 0.952588 0.000194 0.000978 0.001172
	65.0 0.987854 0.987887 0.000011 0.000168 0.000179
	66.0 0.997472 0.997460 0.000000 0.000026 0.000027
	67.0 0.999568 0.999568 0.000000 0.000003 0.000003
	68.0 0.999943 0.999943 0.000000 0.000000 0.000000
	69.0 0.999994 0.999994 0.000000 0.000000 0.000000

Table 9 .

 9 MA model: m 1 = 20, N 1 = 1000, X i = 0.3 Xi + 0.1 Xi+1 + 0.5 Xi+2 , IT ER app = 10 6 , IT ER sim = 10 5

	n	Sim	Approx	Eapp	Esim	E total
	11 0.582252 0.584355 0.011503 0.003653 0.015156
	12 0.770971 0.771446 0.002319 0.001691 0.004010
	13 0.889986 0.889431 0.000434 0.000733 0.001167
	14 0.951529 0.951723 0.000073 0.000297 0.000370
	15 0.980653 0.980675 0.000011 0.000113 0.000124
	16 0.992827 0.992791 0.000001 0.000040 0.000042
	17 0.997486 0.997499 0.000000 0.000013 0.000014
	18 0.999186 0.999188 0.000000 0.000004 0.000004
	19 0.999754 0.999754 0.000000 0.000001 0.000001
	20 0.999930 0.999930 0.000000 0.000000 0.000000


Ting Wang
email: ting.wang@inria.fr

Jean-Pierre Merlet

Guillaume Sacco
email: sacco.g@chu-nice.fr

Philippe Robert
email: robert.ph@chu-nice.fr

Jean-Michel Turpin
email: turpin.jm@chu-nice.fr

Bernard Teboul
email: bernard.teboul@libertysurf.fr

Audrey Marteu
email: marteu.a@chu-nice.fr

Olivier Guerin
email: guerin.o@chu-nice.fr

Walking analysis of young and elderly people by using an intelligent walker ANG

Keywords: Intelligent walker, Walking analysis, Walking quality index, Elderly

This paper proposes a new method to analyze human walking by using a 3-wheels rollator walker instrumented with encoders and a 3D accelerometer/gyrometer. In order to develop walking quality index and monitor the health state of elderly people at home, the walking of 23 young adults and 25 elderly people (> 69 years) with the help of the walker, are compared. The results show that many general walking indicators such as walking speed, stride length do not present obvious difference between the two groups, but that new indicators obtained by using the walker measurements and not

Introduction

The elderly population is growing fast all over the world. Population ageing will cause significant challenges of care giving. One of the problems that affect the most of the elderly population is the reduction of mobility.

Therefore, personal assistance mobility devices are strongly desired to keep the elderly independent. Among the possible assistance devices, the walkers have large number of users because of its simplicity while using the person's remaining locomotion capability in order to move. Besides of the physical benefits of maintaining the standing position, there are also other important psychological benefits, such as maintaining self-esteem and social relationship.

There are many studies and projects regarding advanced versions of walkers. According to the user's needs, the functions of the proposed walkers are not restricted to their primary task, i.e. physical support and mobility assistance. There are other functions such as sensorial assistance, cognitive assistance and health monitoring [START_REF] Frizera | The smart walkers as geriatric assistive device. the simbiosis purpose[END_REF]. These walkers focus on mobility assistance [START_REF] Alwan | Stability margin monitoring in steering-controlled in-telligent walkers for the elderly[END_REF], [START_REF] Bostelman | Robotic patient transfer and rehabilitation device for patient care facilities or the home[END_REF], [START_REF] Glover | A robotically-augmented walker for older adults[END_REF], sit-to-stand transfer [START_REF] Chugo | A moving control of a robotic walker for standing, walking and seating assistance[END_REF], [START_REF] Pasqui | Elderly people sit to stand transfer experimental analysis[END_REF], navigation help [START_REF] Grasse | Assisted navigation for persons with reduced mobility: path recognition through particle filtering (condensation algorithm)[END_REF], [START_REF] Kulyukin | iWalker: toward a rollator mounted wayfinding system for the elderly[END_REF], obstacle avoidance and fall detection [START_REF] Hirata | Fall prevention of passive intelligent walker based on human model[END_REF]. Besides of these, there are other multifunctional walkers such as PAMM SmartWalker [START_REF] Spenko | Robotic personal aids for mobility and monitoring for the elderly[END_REF], which was designed to offer extra support for walking, guidance, scheduling (reminding the time of medicines, for an example) and health monitoring for elderly users. The Medical Automation Research Center (MARC) smart walker [START_REF] Wasson | Effective shared control in cooperative mobility aids[END_REF], which was installed a pair of tridimensional force/torque sensors on it's handles, can be used to determine gait characteristics such as the heel strike, toe-off, double support, and single support [START_REF] Alwan | Passive derivation of basic walker-assisted gait characteristics from measured forces and moments[END_REF], [START_REF] Alwan | Basic walkerassisted gait characteristics derived from forces and moments exerted on the walker's handles: results on normal subjects[END_REF].

To study the extension of the functions of walkers we have developed our own family of walking aids, the ANG family [START_REF] Merlet | ANG, a family of multi-mode, low cost walking aid[END_REF], i.e., ANG-light and ANG-2.

It has been tested that the walker ANG-2 can perform multifunctions such as navigation, street mapping, fall detection/prevention and autonomous object recovery [START_REF] Merlet | ANG, a family of multi-mode, low cost walking aid[END_REF]. Our walker is also a communication device that can receive and emit information (e.g., fall detection for emission, change of user profile for reception). In addition, the walker is able to download and execute new services according to the end-user's trajectory of life (e.g., specific rehabilitation program). It is therefore much more flexible in its application than other walkers. In this work we use the simplest version, ANG-light (Fig. 1) which is based on a commercially available 3-wheels Rollator. It has been instrumented with encoders at the two fixed rear wheels and a 3D accelerometer/gyrometer at the front with the purpose of determining the walker's trajectory on a 24/24 basis. A small, low energy consumption fit-pc computer manages the measurements and records all the data. Additionally the walker has GPS, GSM and infrared network connection that are not used in this paper. Compared with the walkers proposed above, our walker is low cost, simple to be used at home and its functions can be easily extended. This paper will present how it can be used for medical monitoring of walking patterns and what kind of medical information may be obtained. Many studies have examined the effect of age on the walking by comparing younger with older adults [START_REF] Hollman | Age-related differences in spatiotemporal markers of gait stability during dual task walking[END_REF], [START_REF] Menz | Age-related differences in walking stability[END_REF], [START_REF] Woledge | The variable component of lateral body sway during walking in young and older humans[END_REF], [START_REF] Owings | Variability of step kinematics in young and older adults[END_REF], [START_REF] Grabiner | Age-related changes in spatial and temporal gait variables[END_REF], [START_REF] Hausdorff | Increased gait unsteadiness in community-dwelling elderly fallers[END_REF], [START_REF] Barbara | Gait characteristics of younger-old and older-old adults walking overground and on a compliant surface[END_REF], [START_REF] Su | Age-associated differences in the gait pattern changes of older adults during fast-speed and fatigue conditions: results from the baltimore longitudinal study of ageing[END_REF]. These studies calculated some gait parameters, such as step length, gait cycle, step width, cadence and gait speed [START_REF] Su | Age-associated differences in the gait pattern changes of older adults during fast-speed and fatigue conditions: results from the baltimore longitudinal study of ageing[END_REF], [START_REF] Callisaya | Sex modifies the relationship between age and gait: a population-based study of older adults[END_REF], [START_REF] Zijlstra | Assessment of spatio-temporal gait parameters from trunk accelerations during human walking[END_REF], [START_REF] Studenski | Gait speed and survival in older adults[END_REF]. Usually gait speed or walking velocity is regarded as a very important indicator of health. Some studies claimed that older subjects exhibited significantly reduced gait speed compared to younger adults [START_REF] Hollman | Age-related differences in spatiotemporal markers of gait stability during dual task walking[END_REF], [START_REF] Menz | Age-related differences in walking stability[END_REF], [START_REF] Su | Age-associated differences in the gait pattern changes of older adults during fast-speed and fatigue conditions: results from the baltimore longitudinal study of ageing[END_REF]. Other studies showed that there were little or no differences in the gait speed between the healthy younger and elderly people [START_REF] Hausdorff | Increased gait unsteadiness in community-dwelling elderly fallers[END_REF], [START_REF] Grabiner | Age-related changes in spatial and temporal gait variables[END_REF]. In fact, [START_REF] Studenski | Gait speed and survival in older adults[END_REF] has showed that the measured gait speed is based on age, sex [START_REF] Callisaya | Sex modifies the relationship between age and gait: a population-based study of older adults[END_REF], use of mobility aids, chronic conditions, smoking history, blood pressure, body mass index, and hospitalization. Therefore, the classical gait parameters are not sufficient to monitor the user's health. Some studies have considered also gait variability [START_REF] Owings | Variability of step kinematics in young and older adults[END_REF], [START_REF] Callisaya | Ageing and gait variability-a population-based study of older people[END_REF], [START_REF] Grabiner | Age-related changes in spatial and temporal gait variables[END_REF], [START_REF] Hausdorff | Gait variability and fall risk in community-living older adults: a 1-year prospective study[END_REF]. The variability of gait parameters can be characterized by the coefficient of variation (CV) of kinematic gait parameters [START_REF] Dubost | Relationship between dual task related changes in stride velocity and stride time variability in healthy older adults[END_REF], [START_REF] Hollman | Age-related differences in spatiotemporal markers of gait stability during dual task walking[END_REF], which is an index of gait stability and complexity. The increased variability of gait parameters corresponds to decreased gait stability, complexity and increased risk of falling.

However, gait instability is multifactorial and the results of previous studies are often inconsistent according to the conditions of experiment. Therefore, we need to do more tests and find more pertinent indicators of walking. As proposed in [START_REF] Thomas | Evaluation of gait disorders in the elderly[END_REF], at least the following components of a person's gait should be considered: initiation of gait, step length, height, and symmetry, step continuity, step path, trunk motion, walking stance, turning, and heel-to-toe walking. Presently, although some studies began to analyze other gait characteristics such as medial-lateral displacement, center of mass [START_REF] Hahn | Can motion of individual body segments identify dynamic instability in the elderly[END_REF] and foot placement [START_REF] Barbara | Gait characteristics of younger-old and older-old adults walking overground and on a compliant surface[END_REF], [START_REF] Rebula | Measurement of foot placement and its variability with inertial sensors[END_REF], they are still not sufficient to describe comprehensively walking motion. This paper proposes a new method of walking analysis by using our instrumented smart walker. A 10 meters straight line walking test has been performed by two groups of younger and elderly people. Thanks to encoders and a 3D accelerometer/gyrometer, we can calculate accurately the classical gait parameters such as gait cadence, walking speed and stride cycle, stride length and their variability without the subjectivity of an human examination. We also can obtain the trajectory of the walker and therefore compare it with the reference trajectory. This comparison will allow to establish several original gait parameters which have not been considered in previous studies and the walking accuracy of two groups of people can be compared. We define a mobile frame attached to the walker with X 1 being the forward direction and Y 1 being the lateral direction (see Fig. 3). Because of the elasticity of the wheels, when leaning forward to push the walker there will be a clock-wise rotation around the Y 1 axis and when leaning on the right (left) handle a rotation around the X 1 axis should be observed. Therefore, by using the measured angular velocity of the walker, the proportion of forward/backward support force and left/right support force on the handles during the whole walking can be estimated. Overall, using our walker the gait characteristics can be described more comprehensively. Secondly, a drawback of most studies is that these measures are obtained on a reduced space with specialized laboratory equipment such as motion capture systems and instrumented walkways, which may not be available in many clinics and certainly not during daily activities. In contrast, our walker can be easily used at home and outdoors, so it is possible to develop it for individual medical monitoring of walking patterns at any time of the day and in any context. Moreover, motion capture systems don't have large workspace and are not very accurate. Conversely, we have measured the error on the final pose of the walker after a 10 meters walk: the positioning error was less than 1 cm and 0.1 degree in orientation, which is much better than a motion capture system on such large workspace. Besides of this, a very small abnormal limb motion during one step, which characterizes an emerging pathology, cannot be detected by motion capture but can be detected on the walking trajectory.

Description of experiments

Physical functioning tests have showed significant aged-related differences for older adults [START_REF] Isles | Normal values of balance tests in women aged 20-80[END_REF]. Several classical tests used to assess the mobility of elderly people are the 10m walk test (10mWT) (measure: time duration) [START_REF] Salbach | Responsiveness and predictability of gait speed and other disability measures in acute stroke[END_REF], Timed Up and Go test (TUG) (measure: time duration) [START_REF] Podsiadlo | The timed 'up & go': a test of basic functional mobility for frail elderly persons[END_REF], Tinetti Test (TT) (analysis of gait parameters through a video) [START_REF] Tinetti | Performance-oriented assessment of mobility problems in elderly patients[END_REF]. Such tests are easy to implement but are basically global (the time for the 10mWT and the TUG may be identical for two subjects which have however very different walking patterns) or are subjective (for the TT [START_REF] Panella | Towards objective quantification of the tinetti test[END_REF]). Furthermore these tests are performed only during medical visits and consequently are not appropriate to detect abnormal events in the walking patterns. Hence we have decided to examine if the measurements of our walking aid allow one to refine the output of the above walking tests.

For that purpose we have led a large scale experiment that was approved by the regional ethical committee (Comité de Protection des Personnes).

In this paper only the results of a 10mWT will be presented. Each subject was asked to walk along a 10m straight line trajectory with the help of the walker. The experiment takes place at INRIA and at Nice hospital. The subjects were 23 INRIA members (with age between 25 and 65 years, mean value 32) and 25 elderly people (age over 65 years) recruited by Nice hospital (see Fig. 2). No subject has pathological walking deceases. All the subjects were asked to perform twice the trajectory with the walking aid.

Methods

This section will explain how the walker can obtain the walking trajectory and determine the stride. During all the measurements, the calculation of walking trajectory and the detection of stride are the two most important issues as the measurement of all the gait parameters and their variability are based on the detection of stride.

Calculation of the trajectory

As shown in Fig. 3, the origin of the walker frame O 1 X 1 Y 1 is supposed to be the position of the middle point between two rear wheels. The position of the walker in a reference frame OXY is described by [x, y, θ], where θ describes the walking direction of the rollator and represents the angle between the horizontal axis of two rear wheels and the X axis. In our experiment

dL = 2πr 4C • 360 ∆ L (1)
and

dR = 2πr 4C • 360 ∆ R ( 2
)
where r is the radius of the rear wheel and C is a constant parameter of the transformation between the value of encoder and the wheel radius. The change of the direction angle θ during the (k + 1) th sampling time can be estimated as:

dθ = dL -dR D , (3)
where D is the distance between the two rear wheels.

According to the kinematic model shown in Fig. 3, the changes of the walker's position can be obtained as follows [START_REF] Lee | Accurate calibration of kinematic parameters for two wheel differential mobile robots[END_REF]:

dx = dL + dR 2 sin(θ k + dθ 2 ) ( 4
)
dy = dL + dR 2 cos(θ k + dθ 2 ) (5)
Finally, the new position of the walker can be calculated by using:

         x k+1 = x k + dx y k+1 = y k + dy θ k+1 = θ k + dθ (6)
Using the above equations, the trajectory of the walker can be determined by using the encoders. The experiments have shown that after a straight line walking of 10 meters the estimated positioning has an absolute accuracy better than 1cm.

Detection of the stride

The instruments generally used to evaluate human's gait are pedometers, accelerometers or gyrometers. To be appropriate for long-term measurements in everyday environments, these devices should be practical and not interfere with normal movement behaviour. Pedometers are small, easy to use and count the number of steps. The Yamax Digi-Walker SW-200 is considered the most accurate electronic pedometer, but its precision decreases at slower walking speeds, making it less suitable for seniors with low physical fitness or gait abnormalities [START_REF] Dijkstra | Detection of walking periods and number of steps in older adults and patients with parkinson's disease: accuracy of a pedometer and an accelerometry-based method[END_REF]. We have also used a pedometer (Omron) during the experiment and we have noticed large errors on the number of steps.

Accelerometers are utilized to detect the walking stride in many studies [START_REF] Zijlstra | Assessment of spatio-temporal parameters during unconstrained walking[END_REF], [START_REF] Huang | An orientation free adaptive step detection algorithm using a smart phone in physical activity monitoring[END_REF]. Most of methods use the peak value of forward acceleration to detect the walking cycle. However, some steps often does not lead to a high-peak forward acceleration, and hence they are not counted although there is displacement during these periods. A recent study [START_REF] Rebula | Measurement of foot placement and its variability with inertial sensors[END_REF] used thresholds on the magnitude of the gyroscope and accelerometer signals to identify the zero velocity instant and regarded it as the end of a step.

Our walker ANG also uses the gyrometer data to detect the walking stride. An interesting contribution of ANG is that it allows one to differentiate the right and left steps. Indeed when the subject is on the left (right) support phase the walking aid rotates on the left (right). Hence the rotational velocity of the walker around the vertical axis, which can be easily obtained by the gyrometer, is used to detect the walking stride. Its zero value instant is regarded as the end of a step. An example of rotational velocity for an elderly people is shown in Fig. 4.

Since the position of the walker at every moment has been calculated by using the method presented in Subsection 3.1, the displacement of the walker during every step, which is regarded as the step length of the subject, can be easily calculated as soon as all the steps are detected, as shown in Fig. 5.

Accordingly, all gait speed characteristics (such as mean value, minimum and maximum value) can be obtained for each step. Moreover, other spatialtemporal gait parameters such as minimum and maximum acceleration can be analyzed. With a sampling time of 1ms for the encoders and 4.8ms for the gyrometer, we may obtain a quiet reasonable accuracy on these parameters.

Results

In order to analyze the result comprehensively, the walking accuracy, ability and stability of the younger and elderly people are compared respectively in the following subsections.

Comparison of the walking accuracy

Using the calculation method proposed in section 3.1, the trajectories of all the subjects for the 10mWT are presented in Fig. 6. Here and in the following figures the younger adults' trajectories are presented in red while the trajectories of the elderly are presented in blue. The reference trajectory is the horizontal axis and the vertical axis is scaled to illustrate the lateral deviations between the real and reference trajectories. Fig. 6 clearly shows that the elderly subjects have larger deviations than the younger.

Several indicators about the gait can be calculated from a trajectory, such as the maximum and mean value of the lateral deviations between the real and reference trajectory, the domain of the later deviation, the area between the real and reference trajectory, and the relative Standard Deviation (SD)

values. Detailed results are given in the Appendix. Fig. 7 and Fig. 8 show the maximum lateral deviation and the area between the real and reference trajectory respectively, where the results of every group of subjects are sorted into ascending order. These figures illustrate that the results of the elderly subjects have a significantly larger deviation from the reference trajectory than that of the young subjects. For example, the mean value (with a standard deviation) of the maximum lateral deviation for the elderly people is 11.048 ± 5.99 cm while that for the young people is only 3.963 ± 2.301 cm. Similarly, the mean value of the traveled area relative to the reference trajectory for the elderly people is 5930.639 ± 3218.4 cm 2 while that for the young people is 2085.702 ± 1708.313 cm 2 . [START_REF] Hahn | Can motion of individual body segments identify dynamic instability in the elderly[END_REF] has proposed that there exists significant group difference in the medio-lateral displacement of center of mass between healthy elderly adults and elderly patients. It is consistent with our result and they reveal that walking accuracy can be regarded as a pertinent walking quality index.

In addition, several other indicators such as the traveled Manhattan distance and the absolute mean orientation angle presented in Appendix also can be used to measure the walking accuracy of the subjects, and we found that their values for the elderly people are significantly larger than that of the younger people. For the elderly people the mean value of the Manhattan distance is 1070.003 ± 37.077 cm while for the young people this value is 1020.863 ± 14.885 cm. Similarly, the mean value of the absolute mean orientation angle for the elderly people is 1.538 ± 0.553 • while that for the young people is 0.951 ± 0.311 • . In summary, the lateral motion of the elderly is larger than the younger, and the indicators that can be used to characterize this deviation are:

• the relative values of the lateral deviations between the real and reference trajectory,

• the area between the real and reference trajectory,

• the travelled Manhattan distance,

• the relative values of the orientation angle of the walking aid.

Comparison of the walking ability

By using the walker, most of gait parameters presented in the usual walking test can be calculated or estimated, such as gait cycle, gait or walking speed, step length, cadence and forward acceleration. Detailed results are given in Table A.19 in the Appendix. Although the step width cannot be calculated precisely, the analysis of the walker' lateral motion in the previous section can reflect the characteristic of the subjects' step width.

In addition, the instantaneous walking velocity can be derived from the encoder measurements. The instantaneous walking velocity is given in Fig. 9 which shows that there is no obvious difference between the elderly and young subjects. [START_REF] Grabiner | Age-related changes in spatial and temporal gait variables[END_REF] and [START_REF] Hausdorff | Increased gait unsteadiness in community-dwelling elderly fallers[END_REF] also mentioned that the age-related differences in walking velocity were not significant. It has been estimated that a comfortable walking speed for young adult lies in the range 130 cm/s-160 cm/s while for elderly people this speed is given by the formula 117 -0.4 × age [START_REF] Salbach | Responsiveness and predictability of gait speed and other disability measures in acute stroke[END_REF]. The mean speed value for elderly people is coherent with this formula while the result of the younger adults is lower than the normal walking speed. Experiences without the walking aid have shown that the younger subjects were presenting a mean velocity that was close to the normal walking speed. Our interpretation is that elderly people are more familiar with walking aids and have walking patterns that benefit from such an aid, while younger people have a more dynamic pattern that is jeopardized by the aid. This can explain why the maximum velocities of the younger are higher, as shown in Fig. 10.

The results of step period and step length for the two groups are given in Fig. 11 and Fig. 12. Our previous work [START_REF] Merlet | Using a robotized aid for walking analysis: experiments and preliminary results[END_REF] has shown that there is no difference between the left steps and the right steps so here and in the following contents the results of two steps are put together. The results illustrate that there is almost no difference between the two groups which Next we have considered the mean value of forward acceleration shown in Fig. 13. It illustrates that the forward accelerations of the elderly are larger than that of the younger. In addition, almost 70% of the younger's (16 of 23) mean forward accelerations are less than zero while for the elderly this number is only 40% (10 of 25). Therefore, we can deduce that the minimum velocity of the younger is less than that of the elderly although their mean speed is almost the same. As a result, the elderly subjects can use less time to arrive to the goal, as shown in Fig. 14. This contradicts the usual assumption (which is based on the 10mWT) that the walking time may be used to evaluate the walking ability. In summary, the elderly people have similar walking speed, step length, step period as the younger people, but there are three indicators for which the difference exists between the younger and elderly adults:

• maximum instantaneous walking velocity,

• mean value of the forward acceleration,

• time used for the total test.

Comparison of the walking stability

Gait variability is an index of gait stability and complexity. The increased variability of gait parameters corresponds to decreased gait stability, complexity and increased risk of falling. Gait variability is defined as changes in gait parameters from one stride to the next. It can be characterized by the coefficient of variation (CV) of kinematic gait parameters [START_REF] Hollman | Age-related differences in spatiotemporal markers of gait stability during dual task walking[END_REF], [START_REF] Dubost | Relationship between dual task related changes in stride velocity and stride time variability in healthy older adults[END_REF]. The coefficient of variation (CV) is defined as the ratio of the standard deviation (SD) to the mean, i.e., for a set of gait parameter A, it's CV is:

CV (A) = SD(A) mean(A) . (7)
Hence CV shows the extent of variability in relation to mean of the population. Here we have compared the CV of step length, step period and walking speed between the two groups. The results are given in Fig. 15 and detailed information are given in Table A.19 in the Appendix. Fig. 15 shows that the CV of step length for the two groups are similar. The comparison of the CV of stride time illustrates the same characteristic thus it was not presented here. Our results are similar to those of some previous studies [START_REF] Grabiner | Age-related changes in spatial and temporal gait variables[END_REF], [START_REF] Hausdorff | Increased gait unsteadiness in community-dwelling elderly fallers[END_REF], [START_REF] Gabell | The effect of age on variability in gait[END_REF],

which presented that there were no significant differences between young and elderly healthy people in CV of step length and stride time. For the elderly people the mean value of CV of step length is 0.545 ± 0.194 cm and for the young people this value is 0.513 ± 0.143 cm. Next, the CV of walking speed is shown in Fig. 15, as expected the walking speed of elderly exhibits less variability than that of younger adults but the difference is small. Exactly, the mean values of CV of walking speed of the elderly and young people are respectively 0.206 ± 0.09 cm/s and 0.238 ± 0.07 cm/s.

As indicated in the introduction, other information about the pressure on the handles can be used to analyze the walking stability by measuring the angular velocity around the forward direction X 1 and the lateral direction Y 1 (see Fig. 3). We were surprised that the angular velocity measurements were sensitive enough to estimate changes in the forward/backward support force (change on the angular velocity around Y 1 ) and on the left/right support force (change on the angular velocity around X 1 ).

Hence we have been able to determine the respective percentage of forward/backward and left/right support with a reasonable accuracy without any force sensors in the handles [START_REF] Bachschmidt | Development of an instrumented walker for measurement of unilateral hand loads[END_REF]. Fig. 17 shows the percentage of forward support. It is interesting to note that the results of the younger people are much farther away from 50% than that of the elderly people. Exactly, the mean value of the respective percentage of forward support of the young people is 61.902% ± 5.093% while that of the elderly people is 58.160% ± 4.378%.

That means for younger people the difference between forward and backward support is larger. Surprisingly the younger adults are leaning significantly more on the aid than the elderly people.

Based on the above analysis, it appears that the following three indicators should be investigated in the future: • variability of walking speed,

• percentage of forward/backward support,

• percentage of right/left support.

Other applications of the walker ANG

Besides of gait analysis, our walker ANG [START_REF] Merlet | ANG, a family of multi-mode, low cost walking aid[END_REF] can perform multifunctions such as navigation, fall detection/prevention and street mapping. Here we will give an example of its application for updating collaborative maps with sideway slopes and location of lowered kerbs.

The walker is instrumented with an accelerometer that measures both the gravity and the walker's acceleration, which is usually much lower than the gravity. Hence an appropriate processing of the measurement allows to determine the direction of gravity in the walker's frame. In turn this directly allows to determine the roll and tilt of the walker. According to these information, we can measure the slope of sideways, detect lowered kerbs and measure the quality of the sidewalk surfaces. For example, Fig. 17 shows a qualified map of INRIA Sophia Antipolis, where the red lines denote the roads with high slope and the green lines denote the almost horizontal roads that we walked on with the walker. Post-processing of this information allows to report the results in OpenStreetMap to produce an ad-hoc map. Hence it is significant for the end-users of walker or wheelchair because they can share more information of a city while planners can take into account sideway slopes and location of lowered kerbs to determine an optimal itinerary for them.

Conclusions

This paper has proposed a gait analysis method based on the use of an instrumented walker. The results of a 10 m straight line test for 23 younger people and 25 elderly people were compared comprehensively. The comparison includes the relative information about walking accuracy, ability and stability. Several important indicators that exhibit significant difference between the two groups were obtained, such as the maximum lateral deviations between the real and reference trajectory, the area and the Manhattan distance between the real and reference trajectory. For the elderly people these indicators are much larger than that of the young people, and exhibit also significant difference within the same group.

As for the gait parameters describing the walking ability, it appears that there is no obvious difference in step length, step period and walking speed between two groups. The instantaneous walking velocity has been obtained and we found out the maximum instantaneous walking velocity of the younger people is a little larger than that of the elderly people. On the other hand, surprisingly, the time required to perform the trajectory is usually lower for the elderly than for young people. Our explanation is that the normal gait of younger adults was more affected by the presence of the walker than for elderly that are more accustomed with the use of the walker. In addition, when we tried to use the variability of gait parameters to analyze the walking stability, there are similar results for the two groups in the variability of step length and step period while for the younger subjects the variability of the walking speed is larger.

Moreover, we also found that the younger adults are leaning significantly more on the aid than the elderly people. Hence the influence of the walker on the walk pattern of the younger people should be investigated. Besides of this, another walking test with a returning trajectory for two groups people will be studied. We also want to examine if a learning process may be implemented in order to characterize the walking pattern at a given time and customize the walking analysis software in order to better determine future trends. Finally, it is worth mentioning that the present analysis method depends on the measurement of sensors with noise that do not satisfy usual statistical hypothesis (e.g., they are not Gaussian). In the future we plan to take measurement uncertainties into account by using interval analysis for obtaining the indicators as interval values that will be guaranteed to include the real values, with the advantage that the width of the interval will be a measure of the quality of the result.

Results of the elderly subjects

Figure 1 :

 1 Figure 1: The walking aid ANG-light

Figure 2 :

 2 Figure 2: The walking aid ANG-light at Nice Hospital

1 Figure 3 :

 13 Figure 3: Simple kinematic model of the walker

Figure 4 :

 4 Figure 4: The rotational velocity of the walker around the vertical axis as function of time ( • /s). One step is finished when it passes zero.

Figure 5 :

 5 Figure 5: Displacement (cm) of the walker during every step as function of time (s). The results of left steps and right steps were put together and they appeared alternately.

Figure 6 :

 6 Figure 6: Trajectory of the subjects in the xy plane, where the blue color denotes the elderly subjects and the red one denotes the young subjects. The reference trajectory is the horizontal line: y = 0.

Figure 7 :Figure 8 :

 78 Figure 7: The maximum lateral deviation (cm) between the real and reference trajectory, where the blue color denotes the elderly subjects and the red one denotes the young subjects. The results of every group of subjects are sorted into ascending order.

Fig. 10 givesFigure 9 :

 109 Figure 9: Instantaneous walking velocity (cm/s) of elderly and young adults. Blue color denotes the elderly subjects and the red one denotes the young subjects. In order to estimate it more precisely, only the middle part of the trajectory is used to do the derivation.

Figure 10 :

 10 Figure 10: Maximum instantaneous walking velocity (cm/s) of elderly and young adults.Blue color denotes the elderly subjects and the red one denotes the young subjects. The results of every group of subjects are sorted into ascending order.

Figure 11 :

 11 Figure 11: Mean value of step period (s), where the blue color denotes the elderly subjects and the red one denotes the young subjects. The results of every group of subjects are sorted into ascending order.

Figure 12 : 2 )Figure 13 :

 12213 Figure 12: Mean value of step length (cm), where the blue color denotes the elderly subjects and the red one denotes the young subjects. The results of every group of subjects are sorted into ascending order.

Figure 14 :

 14 Figure 14: Time used for 10mWT (s) of elderly and young adults. Blue color denotes the elderly subjects and the red one denotes the young subjects. The results of every group of subjects are sorted into ascending order.

Figure 15 :

 15 Figure 15: Coefficient of variation for the step length and walking speed, where the blue color denotes the elderly subjects and the red one denotes the young subjects. The results of every group of subjects are sorted into ascending order.

Figure 16 :

 16 Figure 16: Percentage of the forward support, where the blue color denotes the elderly subjects and the red one denotes the young subjects. The results of every group of subjects are sorted into ascending order.

Figure 17 :

 17 Figure 17: Qualified map of INRIA Sophia Antipolis by using the walker ANG, where the red lines denote the roads with high slope and the green lines denote the almost horizontal roads that we walked on with the walker.

Figure A. 18 :

 18 Figure A.18: Result of trajectory for the elderly (1 -25) and younger subjects (26 -48).

Figure A. 19 :

 19 Figure A.19: Result of gait parameters for the elderly (1 -25) and younger subjects(26 -48).

  :///auto/sop-nas2a/u/sop-nas2a/vol/home_coprin/tiwan...

					Analyse of two legs Results of the elderly (1-25) and young subjects (26-48)
	Subject	Step number	Mean step length (cm)	Minimum step length (cm)	Maximum step length (cm)	Variation of step length	Mean walking cycle (s)	Minimun walking cycle(s)	Maximum walking cycle(s)	Variation of walking cycle	Mean walking speed(cm/s)	Variation of walking speed	Mean foward acceleration (m/s^2)	Standard deviation of foward acceleration (m/s^2)
	1	12	82.921 52.920	202.720	0.551	0.734	0.482	1.692	0.501	112.608	0.129	-0.102	0.208
	2	25	42.404 7.022	136.153	0.658	0.389	0.078	1.073	0.553	109.091	0.268	0.031	0.233
	3	18	56.293 1.427	101.454	0.555	0.461	0.001	0.818	0.507	217.939	1.780	-0.053	0.385
	4	23	45.676 6.320	69.053	0.319	0.610	0.081	1.045	0.349	75.833	0.106	0.219	0.144
	5	17	72.040 1.126	235.058	0.872	0.661	0.012	1.974	0.800	107.812	0.135	0.107	0.213
	6	20	51.100 20.202	99.100	0.459	0.577	0.200	0.960	0.406	89.776	0.193	0.075	0.120
	7	14	74.004 17.926	130.485	0.360	0.581	0.328	1.233	0.400	127.207	0.207	-0.231	0.261
	8	19	58.561 27.199	86.262	0.276	0.523	0.339	0.689	0.201	110.614	0.131	0.035	0.198
	9	17	59.859 3.074	207.145	0.871	0.504	0.016	1.513	0.740	126.354	0.288	-0.034	0.259
	10	21	47.998 2.916	152.318	0.728	0.480	0.028	1.577	0.721	100.605	0.168	0.014	0.189
	11	20	51.994 1.416	100.955	0.473	0.465	0.005	0.866	0.447	121.855	0.411	-0.132	0.494
	12	31	37.701 0.669	110.098	0.687	0.607	0.049	2.673	0.833	68.111	0.374	0.165	0.199
	13	27	40.483 1.137	99.439	0.678	0.451	0.006	1.016	0.658	94.144	0.272	0.108	0.240
	14	29	33.312 5.858	97.523	0.588	0.379	0.055	1.006	0.568	90.361	0.167	0.131	0.266
	15	18	56.441 9.332	91.475	0.455	0.529	0.079	0.915	0.476	109.992	0.136	0.087	0.224
	16	19	48.610 8.972	97.786	0.438	0.401	0.255	0.706	0.311	118.447	0.267	-0.119	0.459
	17	20	52.996 1.223	106.687	0.529	0.592	0.006	1.748	0.679	99.928	0.262	0.173	0.173
	18	17	62.476 6.570	89.165	0.351	0.545	0.107	0.789	0.336	112.329	0.179	-0.077	0.291
	19	24	43.558 15.888	87.297	0.437	0.511	0.174	1.443	0.551	87.792	0.080	0.198	0.191
	20	18	60.378 18.027	75.697	0.247	0.548	0.438	0.657	0.117	110.144	0.223	-0.099	0.240
	21	17	57.511 6.455	85.754	0.512	0.375	0.041	0.625	0.502	151.234	0.104	0.018	0.198
	22	20	57.823 19.110	89.430	0.315	0.432	0.220	0.606	0.261	134.145	0.186	-0.066	0.408
	23	17	69.029 20.903	198.043	0.659	0.739	0.234	2.075	0.640	93.700	0.122	0.114	0.195
	24	22	51.904 15.072	145.279	0.597	0.487	0.130	1.211	0.491	103.673	0.181	0.082	0.482
	25	18	56.484 1.473	221.339	0.999	0.577	0.007	2.124	0.946	103.840	0.356	-0.032	0.197
	26	17	63.562 4.135	128.368	0.454	0.786	0.315	2.275	0.630	83.235	0.288	0.022	0.099
	27	23	42.397 2.444	128.813	0.861	0.404	0.017	1.104	0.798	108.786	0.277	-0.101	0.285
	28	21	44.727 13.856	79.080	0.441	0.500	0.242	0.806	0.326	88.666	0.252	-0.086	0.193
	29	24	48.638 6.452	99.169	0.465	0.450	0.045	0.852	0.488	116.145	0.240	0.104	0.317
	30	19	65.546 17.305	169.675	0.646	0.668	0.211	1.453	0.592	99.821	0.261	0.026	0.293
	31	19	60.559 18.135	110.274	0.421	0.528	0.137	0.947	0.425	117.764	0.173	-0.040	0.269
	32	22	48.086 9.139	93.776	0.473	0.472	0.083	1.084	0.559	106.607	0.147	-0.000	0.159
	33	21	54.244 7.390	85.969	0.420	0.578	0.073	0.853	0.365	93.539	0.214	-0.038	0.295
	34	21	52.800 1.019	78.531	0.405	0.554	0.007	0.764	0.334	97.361	0.244	-0.016	0.244
	35	16	71.941 32.556	184.565	0.509	0.655	0.280	1.779	0.540	112.527	0.186	-0.043	0.265
	36	16	68.243 46.257	90.528	0.169	0.591	0.375	1.289	0.374	120.280	0.133	-0.092	0.201
	37	19	62.723 3.583	109.991	0.413	0.551	0.026	0.935	0.350	114.700	0.201	-0.109	0.336
	38	17	56.930 9.549	137.246	0.588	0.610	0.248	1.428	0.486	91.555	0.321	-0.088	0.314
	39	20	58.976 13.937	106.636	0.475	0.662	0.173	1.430	0.476	90.444	0.244	0.021	0.271
	40	25	40.416 1.155	79.878	0.637	0.482	0.004	1.038	0.588	94.237	0.478	0.048	0.207
	41	21	49.478 5.605	126.980	0.733	0.514	0.083	1.248	0.646	95.306	0.266	-0.041	0.227
	42	18	58.275 2.376	129.605	0.635	0.524	0.038	1.143	0.622	109.730	0.261	-0.131	0.331
	43	19	54.777 7.646	111.694	0.555	0.463	0.056	1.253	0.609	123.488	0.252	-0.067	0.276
	44	22	51.677 7.011	79.897	0.397	0.454	0.135	0.747	0.349	112.852	0.241	-0.091	0.363
	45	20	54.936 2.875	79.032	0.372	0.478	0.105	0.662	0.309	110.902	0.240	-0.148	0.395
	46	22	49.220 2.809	101.165	0.651	0.381	0.017	0.895	0.702	137.518	0.180	0.066	0.384
	47	16	64.334 11.461	116.511	0.547	0.479	0.075	1.001	0.595	141.711	0.162	-0.130	0.359
	48	27	43.658 9.531	89.937	0.529	0.578	0.141	1.135	0.488	76.389	0.216	0.066	0.208
	1 of 1												01/31/2014 02:12 PM

file


Sylvain Conchon

Mohamed Iguernelala

Tuning the Alt-Ergo SMT Solver for B Proof Obligations

Keywords: SMT solvers, B Proof Obligations, B Method

In this paper, we present recent developments on the Alt-Ergo SMT-solver to efficiently discharge proof obligations (PO) generated by Atelier B. This includes a new plugin architecture to facilitate experiments with different SAT engines, new heuristics to handle quantified formulas, and important modifications in its internal data structures to boost performances of core decision procedures. Benchmarks realized on more than 10,000 PO generated from industrial B projects show significant improvements.

1 The Alt-Ergo SMT Solver Alt-Ergo is an open-source SMT solver capable of reasoning in a combination of several built-in theories such as uninterpreted equality, integer and rational arithmetic, arrays, records, enumerated data types and AC symbols. It is the unique SMT solver that natively handles polymorphic first-order quantified formulas, which makes it particularly suitable for program verification. For instance, Alt-Ergo is used as a back-end of SPARK and Frama-C to discharge proof obligations generated from Ada and C programs, respectively.

Recently, we have started using Alt-Ergo in the context of the ANR project BWare [START_REF]The BWare Project[END_REF] which aims at integrating SMT solvers as back-ends of Atelier B. The proof obligations sent to Alt-Ergo are extracted from Atelier B as logical formulas that are combined with a (polymorphic) model of B's set theory [START_REF] Mentré | Discharging proof obligations from atelier b using multiple automated provers[END_REF]. This process relies on the Why3 plateform [START_REF] Filliâtre | Why3 -where programs meet provers[END_REF] which can target a wide range of SMT solvers. However, we show (Section 2) on a large benchmark of industrial B projects that it is not immediate to obtain a substantial gain of performances by using only SMT solvers. Without a specific tunning for B, Alt-Ergo together with other SMT solvers compete just equally with Atelier B's solver on those industrial benchmarks.

In this paper, we report on recent developments in Alt-Ergo that significantly improve its capacities to handle PO of Atelier B. Our improvements are: (1) better heuristics for instantiating polymorphic quantified formulas from B model;

(2) new efficient internal data structures; (3) a plugin architecture to facilitate experiments with different SAT engines; and (4) the implementation of a new CDCL-based SAT solver.

Benchmarks

The test-suite of BWare contains 10572 formulas obtained from three industrial B projects provided by ClearSy and Mitsubishi Electric. The first one, called DAB, is an automated teller machine (ATM). The last two ones, called P4 and P9, are obfuscated and unsourced programs.

The formulas generated from these projects are composed of two parts. The first one is the context, a large set of axioms (universally and polymorphic quantified formulas). The second one is the goal, a ground, unique and large formula. A more thorough investigation of the shape of these formulas shows that they mainly contain equalities over uninterpreted function symbols and atoms involving enumerated data types. Only a small portion of atoms contains linear integer arithmetic and polymorphic records. In comparison with other benchmarks coming from deductive program verification plateforms, the average number of axioms, as well as the size of these PO are much larger in this test suite, as shown in the following

Conclusion and Future Works

Originally developed at LRI, it is now maintained and distributed by the OCaml-Pro company.

As can be expected, the BWare project consists of several tasks, which cannot be exhaustively described in this paper due to space restrictions. We therefore choose to focus on two major current lines of work of the project.

The first current line of work is upstream and consists in completing the axiomatization of the B set theory in Why3 in order to be able to consider all the provided proof obligations. This is mainly carried out according to what is described in [START_REF] Mentré | Discharging Proof Obligations from Atelier B Using Multiple Automated Provers[END_REF], i.e. by adding set constructs to the axiomatization and modifying the translator of proof obligations from Atelier B to Why3 accordingly. This line of work is quite important in the project as it will allow us to consider a broad scope of proof obligations related to different application domains and test the scalability of our platform as well.

A second current line of work focuses on the first order provers to make them able to reason modulo the B set theory. To do so, we rely on deduction modulo. The theory of deduction modulo is an extension of predicate calculus, which allows us to rewrite terms as well as propositions, and which is well suited for proof search in axiomatic theories, as it turns axioms into rewrite rules. Both first order provers considered in the project have been extended to deduction modulo to obtain Zenon Modulo [START_REF] Delahaye | Proof Certification in Zenon Modulo: When Achilles Uses Deduction Modulo to Outrun the Tortoise with Shorter Steps[END_REF][START_REF] Delahaye | Zenon Modulo: When Achilles Outruns the Tortoise using Deduction Modulo[END_REF] and iProver Modulo [START_REF] Burel | Experimenting with Deduction Modulo[END_REF]. Both tools have also been extended to produce Dedukti proofs (see [START_REF] Delahaye | Proof Certification in Zenon Modulo: When Achilles Uses Deduction Modulo to Outrun the Tortoise with Shorter Steps[END_REF][START_REF] Delahaye | Zenon Modulo: When Achilles Outruns the Tortoise using Deduction Modulo[END_REF] and [START_REF] Burel | A Shallow Embedding of Resolution and Superposition Proofs into the λΠ-Calculus Modulo[END_REF]), which is natural as Dedukti relies on deduction modulo as well. Currently, most of the efforts in this line of work consist in building a B set theory modulo, which is appropriate for automated deduction and keeps some properties such as cut-free completeness. [START_REF] Déharbe | SMT Solvers for Rodin[END_REF]