R Pierre

Bertrand

Mehdi Fhima

Arnaud Guillin

Off-line detection of multiple change points by the Filtered Derivative with p-Value method

Keywords: Change points detection, Filtered Derivative, Strong approximation

This paper deals with off-line detection of change points, for time series of independent observations, when the number of change points is unknown. We propose a sequential analysis method with linear time and memory complexity. Our method is based at first step, on Filtered Derivative method which detects the right change points as well as the false ones. We improve the Filtered Derivative method by adding a second step in which we compute the p-values associated to every single potential change point. Then we eliminate false alarms, i.e. the change points which have p-value smaller than a given critical level. Next, we apply our method and Penalized Least Squares Criterion procedure to detect change points on simulated data sets and then we compare them. Eventually, we apply the Filtered Derivative with p-Value method to the segmentation of heartbeat time series, and the detection of change points in the average daily volume of financial time series.

INTRODUCTION

In a wide variety of applications including health and medicine, finance, civil engineering, one models time dependent systems by a sequence of random variables described by a finite number of structural parameters. These structural parameters can change abruptly and it is relevant to detect the unknown change points. Both on-line and off-line change point detection have their own relevance, but in this work we are concerned with off-line detection.

Statisticians have studied change point detections since the 1950's and there is a huge literature on this subject, see e.g. the textbooks [START_REF] Basseville | Detection of abrupt changes: theory and application[END_REF]; [START_REF] Brodsky | Nonparametric methods in change-point problems[END_REF]; [START_REF] Csörgo | Limit Theorem in Change-Point Analysis[END_REF]; [START_REF] Montgomery | Introduction to Statistical Quality Control, 3rd Edition[END_REF], the interesting papers of [START_REF] Chen | Inference in a simple change-point model[END_REF]; [START_REF] Steinebach | Heart rate variability. Standards of measurement, physiological interpretation, and clinical use[END_REF], and let us also refer to [START_REF] Hušková | Change point analysis based on the empirical characteristic functions of ranks[END_REF]; [START_REF] Kirch | Bootstrapping sequential change-point tests[END_REF]; [START_REF] Gombay | Monitoring parameter change in ar(p) time series models[END_REF] for update reviews, [START_REF] Fliess | Algebraic change-point detection[END_REF] who propose an algebraic approach for change point detection problem, and [START_REF] Birgé | Minimal penalties for Gaussian model selection[END_REF] for a good summary of the model selection approach.

Among the popular methods, we find the Penalized Least Squares Criterion (PLSC). This algorithm is based on the minimization of the contrast function when the number of change points is known, see [START_REF] Bai | Estimating and testing linear models with multiple structural changes[END_REF], [START_REF] Lavielle | Least-squares estimation of an unknown number of shifts in a time series[END_REF]. When the number of changes is unknown, many authors use the penalized version of the contrast function, see e.g. [START_REF] Lavielle | Detection of multiple change points in multivariate time series[END_REF] or [START_REF] Lebarbier | Detecting multiple change-points in the mean of gaussian process by model selection[END_REF]. From a numerical point of view, the least squares methods are based on dynamic programming algorithm which needs to compute a matrix. Therefore, the time and memory complexity are of order O(n 2 ) where n is the size of data sets. So, complexity becomes an important limitation with technological progress.

Indeed, recent measurement methods allow us to record and to stock large data sets. For example, in Section 5, we present change point analysis of heartbeat time series: It is presently possible to record the duration of each single heartbeat during a marathon race or for healthy people during 24 hours. This leads to data sets of size n ≥ 40, 000 or n ≥ 100, 000, respectively.

Actually, this phenomenon is general: time dependent data are often recorded at very high frequency (VHF), which combined with size reduction of memory capacities allows recording of millions of data.

This technological progress leads us to revisit change point detection methods in the particular case of large or huge data sets. This framework constitutes the main novelty of this work: we have to develop embeddable algorithms with low time and memory complexity. Moreover, we can adopt an asymptotic point of view. The non asymptotic case will be treated elsewhere.

For the off-line detection problem, many papers propose procedures based on the dynamic programming algorithm. Among them, as previously said, we find the Penalized Least Squares Criterion (PLSC) which relies on the minimization of contrast function. It performs with time and space complexities of order O n 2 , see e.g. [START_REF] Bai | Estimating and testing linear models with multiple structural changes[END_REF]; [START_REF] Lavielle | Detection of multiple change points in multivariate time series[END_REF] or [START_REF] Lebarbier | Detecting multiple change-points in the mean of gaussian process by model selection[END_REF]. Then, [START_REF] Guédon | Exploring the state sequence space for hidden markov and semi markov chains[END_REF] has developed a method based on the maximization of the log-likelihood. Using a forward-backward dynamic programming algorithm, he obtains a time complexity of order O (K + 1)n 2 and a memory complexity of order O ((K + 1)n), where K represents the number of change points. These methods are then time consuming when the number of observations n becomes large. For this reason, many other procedures have been proposed for a faster segmentation. In [START_REF] Chopin | Dynamic detection of change points in long time series[END_REF] or [START_REF] Fearnhead | On-line inference for multiple change points problems[END_REF], the problem of multiple change-points is reformulated as a Bayesian model. The segmentation is computed by using a particle filtering algorithm which requires complexity of order O (M × n) both in time and memory, where M is the number of particles. For instance, [START_REF] Chopin | Dynamic detection of change points in long time series[END_REF] obtains reasonable estimations of change points with M = 5000 particles. Recently, a segmentation method based on wavelet decomposition in Haar basis and thresholding has been proposed by [START_REF] Ben-Yaacov | A fast and flexible method for the segmentation of acgh data[END_REF], with time complexity of order O (n log(n)) and memory complexity of order O (n). Though, no mathematical proofs has been given. Moreover, no results has been proposed to detect change points in the variance or in the slope and the intercept of linear regression model.

In this paper, we investigate the properties of a new off-line detection methods for multiple change points, so-called Filtered Derivative with p-value method (FDp-V). Filtered Derivative has been introduced by [START_REF] Benveniste | Detection of abrupt changes in signals and dynamical systems: some statistical aspects[END_REF]. [START_REF] Basseville | Detection of abrupt changes: theory and application[END_REF], next [START_REF] Antoch | Procedures for the detection of multiple changes in series of independent observations[END_REF] propose an asymptotic study and [START_REF] Bertrand | A local method for estimating change points: the "hat-function[END_REF] gives some non asymptotic results. On the one hand, the advantage of Filtered Derivative method is its time and memory complexity, both of order O(n). On the other hand, the drawback of Filtered Derivative method is that if it detects the right change points it also gives many false alarms. To avoid this drawback, we introduce a second step in order to disentangle right change points and false alarms. In this second step we calculate the p-value associated to each potential change point detected in the first step. Stress that the second step has still time and memory complexity of order O(n).

Our belief is that FDp-V method is quite general for large datasets. However, in this work, we restrict ourselves to detection of change points on mean and variance for a sequence of independent random variables and change point on slope and intercept for linear model. The rest of this paper is organized as follows: In Section 2, we describe Filtered Derivative with p-value Algorithm. Then, Section 3 is concerned with theoretical results and FDp-V method for detecting changes on mean and variance. In Section 4, we present theoretical results of FDp-V method for detecting changes on slope and intercept for linear regression model. Finally, in Section 5, we give numerical simulations with a comparison with PLSC algorithm, and we present some results on real data showing the robustness of our method (as no independance is guaranteed). All the proofs are postponed to Section 6.

DESCRIPTION OF THE FILTERED DERIVATIVE WITH p-VALUE METHOD

In this section, we describe the Filtered Derivative with p-value method (FDp-V). First, we describe precisely the statistical model we will use throughout our work. Next, we describe the two steps of FDp-V method:

Step 1 is based on Filtered Derivative and select the potential change points, whereas Step 2 calculate the p-value associated to each potential change point, for disentangling right change points and false alarms.

Our model:

Let (X t ) t=1,...,n be a sequence of independent r.v. with distribution M θ(t) , where θ ∈ R d is a finite dimensional parameter. We assume that the maps t → θ(t) is piecewise constant, i.e. there exists a configuration of change points 0

= τ 0 < τ 1 < • • • < τ K < τ K+1 = n such that θ(t) = θ k for τ k ≤ t < τ k+1 .
The integer K corresponds to the number of change times and (K + 1) to the number of segments. In summary, if j ∈ [τ k , τ k+1 [, the r.v. X j are independent and identically distributed with distribution M θ k . We stress that the number of abrupt changes K is unknown, leading to a problem of model selection. There is a huge literature on change point analysis and model selection see e.g. the monographs [START_REF] Basseville | Detection of abrupt changes: theory and application[END_REF]; [START_REF] Brodsky | Nonparametric methods in change-point problems[END_REF] or [START_REF] Birgé | Minimal penalties for Gaussian model selection[END_REF]. As pointed out in the introduction, large and huge datasets lead to revisit change point analysis by taking into account time and memory complexity of the different methods.

Filtered Derivative:

Filtered Derivative is defined as the difference between the estimators of the parameter θ computed on two sliding windows respectively at the right and at the left of the index k, both of size A, that is specified by the following function:

D(k, A) = θ(k, A) -θ(k -A, A), (2.1)
where θ(k, A) is an estimator of θ on the sliding box [k + 1, k + A]. Eventually, this method consists on filtering data by computing the estimators of the parameter θ before applying a discrete derivation. So, this construction explains the name of the algorithm, so-called Filtered Derivative method.

The norm of the filtered derivative function, say D , exhibits "hats" at the vicinity of parameter change points. To simplify the presentation, without losing generality, we will assume now that θ is a one dimensional parameter, say for example the mean. Thus, the change points can be estimated as arguments of local maxima of |D|, see Figure 2.1 below. Moreover, the size of changes in θ is equal to the height of positive or negative hats of D.

Nevertheless, we remark through the graph of the function |D| that there are not only the "right hats" (surrounded in blue in Figure 2.2) which gives the right change points, but also false alarms (surrounded in green in Figure 2.2). Consequently, we have introduced another idea in order to keep just the right change points. This objective is reached by splitting the detection procedure into two successive steps: In Step 1, we detect potential change points as local maxima of the filtered derivative function. In Step 2, we test wether a potential change point is a false alarm or not. Both steps use estimation of the p-value of existence of change point.

The construction of two different statistical tests and the computation of p-values is detailed below. Step 1: Detection of the potential change points

In order to detect the potential change points, we test the null hypothesis (H 0 ) of no change in the parameter θ (H 0 ) :

θ 1 = θ 2 = • • • = θ n-1 = θ n
against the alternative hypothesis (H 1 ) indicating the existence of one or more than one changes (H 1 ) : There is an integer K ∈ N * and 0 = τ

0 < τ 1 < • • • < τ K < τ K+1 = n such that θ 1 = • • • = θ τ 1 = θ τ 1 +1 = • • • = θ τ 2 • • • = θ τ K +1 = • • • = θ τ K+1 .
where θ i is the value of the parameter θ for 1 ≤ i ≤ n.

In [START_REF] Antoch | Procedures for the detection of multiple changes in series of independent observations[END_REF] or [START_REF] Bertrand | A local method for estimating change points: the "hat-function[END_REF], potential change points are selected as times τ k corresponding to local maxima of the absolute value of the filtered derivative |D(τ k , A)|, when moreover this last quantity exceeds a given threshold λ. However, the efficiency of the approach is strongly linked to the choice of the threshold λ. Therefore, in this work, we have a slightly different approach: we fix a probability of type I error at level p * 1 , and we determine the corresponding critical value C 1 given by

P max k∈[A:n-A] |D(k, A)| > C 1 | H 0 is true = p * 1 .
Of course, such a probability is usually not available, so that we only have the asymptotic distribution of the maximum of |D|, which will be the main part of Section 3 and Section 4. Then, roughly speaking, we select as potential change points, local maxima for which |D(τ k , A)| > C 1 . By doing so, the p-value p * 1 appears as an hyper parameter. The second hyper parameter is the window size A. As pointed out in [START_REF] Antoch | Procedures for the detection of multiple changes in series of independent observations[END_REF]; [START_REF] Bertrand | A local method for estimating change points: the "hat-function[END_REF], Filtered Derivative method works under the implicit assumption that the minimal distance between two successive change points is greater than 2A. An adaptive and automatic choice of A is an open and particularly relevant question, but there is no such result in the litterature. Thus, we need some a priori knowledge on the minimal length between two successive changes.

More formally, we have the following algorithm:

Step 1 of the algorithm 1. Choice of the hyper parameters

• Choice of the window size A from information of the practitioners.

• Choice of p * 1 : First we fix the significance level of type I error at p * 1 . Then, the theoretic expression of type I error, given in Section 3 and Section 4, fixes the value of the threshold C 1 . The list of potential change points, which contains right change points as well as false ones, can be too large. So the level of significance for Step 1 can be chosen large, i.e. p * 1 = 5% or 10%.

2. Computation of the filtered derivative function 3. Determination of the potential change points

• Initialization: Set counter of potential change point k = 0 and τ k = arg max k∈[A,n-A] |D(k, A)|. • While (|D( τ k , A)| > C 1 ) do -k = k + 1 -D(k, A) = 0 for all k ∈ ( τ k -A, τ k + A). -τ k = arg max k∈[A,n-A] |D(k, A)|
We increment the change point counter and we set the values of the function D to zero because the width of the hat is equal to 2A.

• Finally, we sort the vector ( τ 1 , . . . , τ Kmax ) in increasing order. The integer K max represents the number of potential change points detected at the end of the Step 1.

By construction, we have K max ≤ ⌊n/A⌋ -2, where ⌊x⌋ corresponds to the integer part of x.

End of the

Step 1 of the algorithm Step 2: Removing false detection

The list of potential change points ( τ 1 , . . . , τ Kmax ) obtained at step 1 contains right change points but also false detections. In the second step a test is carried out to remove the false detection from the list of change points found at step 1. By doing so, we obtain a subset τ 1 , . . . , τ K of the first list.

More precisely, for all potential change point τ k , we test wether the parameter is the same on the two successive intervals ( τ k-1 + 1, τ k ) and ( τ k + 1, τ k+1 ), or not. Formally, for all 1 ≤ k ≤ K max , we apply the following hypothesis testing where θ k is the value of θ supposed to be constant on the segment ( τ k-1 + 1, τ k ). By using this second test, we calculate new p-values ( p 1 , . . . , p Kmax ) associated respectively to each potential change points ( τ 1 , . . . , τ Kmax ). Then, we only keep the change points which have a p-value smaller than a critical level denoted p * 2 .

(H 0,k ) : θ k = θ k+1 versus (H 1,k ) : θ k = θ k+1
Step 2 must be much more selective with a significance level p * 2 = 10 -3 or 10 -4 . Consequently, Step 2 permits us to remove most of false detections, and so to deduce an estimator of the piecewise constant map t → θ(t), see Figure 2.3 below.

Remark 2.1

We note that τ k , for 1 ≤ k ≤ K max , is a random variable. So that at first sight, the average value of the parameter θ on the interval ( τ k-1 +1, τ k ) may also be a random variable. However, [START_REF] Bertrand | A local method for estimating change points: the "hat-function[END_REF] has showed that, with probability converging to 1, the right change point is in

( τ k + ǫ k , τ k -ǫ k )
, where ǫ k is a given bounded real number. For the detection of change in the mean in a series of independent Gaussian random variables, Bertrand (2000, Theorem 3.1, p.255) has proved that ǫ k = 5(σ/δ k ) 2 + 1 where δ k represents the size of change in the mean at τ k . This result can be easily generalized to a series of independent non Gaussian random variables with finite second moment. Therefore, on a set Ω n with P(Ω n ) → 1 as n → ∞, the average value of θ on the segment ( τ k + ǫ k , τ kǫ k ) becomes not random and actually is equal to θ k .

The time and memory complexity of this second step of the algorithm is still O(n) because we only need to compute and to stock (K max + 1) estimators of parameter θ which are successively compared two by two in order to eliminate false alarms.

To sum up, FDp-V method is a two step procedure, in which the filtered derivative method is applied first and then a second test is carried out to remove the false detection from the list of change points found in Step 1. The first step has both time and memory complexity of order O(n). At the second step, the number of selected potential change points is smaller. As a consequence, both time and memory complexity of Step 2 are still of order O(n).

Let us finish the presentation of the FDpV algorithm by the following remark: in this paper we will restrict ourselves to the detection of change points of a one dimensional parameter, such as mean, variance, slope and intercept of linear regression. However, our algorithm still works in any finite dimension parameter space. Indeed, for each parameter, we can compute the corresponding filtered derivative sequence and its asymptotic distribution and then compare it to the corresponding type I errors p * 1 and p * 2 .

THEORETICAL RESULTS

In this section, we present theoretical results of the FDp-V method for detecting changes on the mean and on the variance. The two following subsections give the asymptotic distribution

Change in the mean

Let (X i ) i=1,...,n be a sequence of independent r.v. with mean µ i and a known variance σ 2 . We assume that the map i → µ i is piecewise constant, i.e. there exists a configuration 0

= τ 0 < τ 1 < • • • < τ K < τ K+1 = n such that E(X i ) = µ k for τ k ≤ i < τ k+1 .
The integer K corresponds to the number of changes. However, in any real life situation, the number of abrupt changes K is unknown, leading to a problem of model selection, see e.g. [START_REF] Birgé | Minimal penalties for Gaussian model selection[END_REF]. Filtered Derivative method applied to the mean is based on the difference between the empirical mean computed on two sliding windows respectively at the right and at the left of the index k, both of size A, see [START_REF] Antoch | Procedures for the detection of multiple changes in series of independent observations[END_REF]; [START_REF] Basseville | Detection of abrupt changes: theory and application[END_REF]. This difference corresponds to a sequence (D 1 (k, A)) A≤k≤n-A defined by

D 1 (k, A) = μ(k, A) -μ(k -A, A) (3.1) where μ(k, A) = 1 A k+A j=k+1 X j is the empirical mean of X on the (sliding) box [k + 1, k + A].

These quantities can easily be calculated by recurrence with complexity O(n). It suffices to remark that

A × D 1 (k + 1, A) = A × D 1 (k, A) + X k+A+1 -2X k + X k-A+1 . (3.2)
First we give in Theorem 3.1 the asymptotic behaviour of the maximum of |D 1 | under null hypothesis of no change in the mean and with size of the sliding windows tending to infinity at a certain rate. In the sequel, we denote by A n the size of the sliding windows and we will always suppose that

lim n→+∞ A n n = 0 and lim n→+∞ (log n) 2 A n = 0. (3.3)
Theorem 3.1 (Change point in the mean with known variance) Let (X i ) i=1,...,n be a sequence of independent identically distributed random variables with mean µ, variance σ 2 and assume that one of the following assumptions is satisfied

(A 1 ) X 1 ∼ N (µ, σ 2 ). (A 2 ) ∃ t > 0 such as E[exp(tX 1 )] < +∞, and lim n→+∞ (log n) 3 A n = 0. (A 3 ) ∃ p > 2 such as E[|X 1 | p ] < +∞, and lim n→+∞ n 2/p log n A n = 0.
Let D 1 be defined by (3.1). Then under the null hypothesis

lim n→+∞ P max k∈[An:n-An] |D 1 (k, A n )| ≤ σ √ A n c n (x) = exp(-2e -x ), (3.4) c n (x) = c n A n -1, x (3.5) c(y, x) = 1 √ 2 log y x + 2 log y + 1 2 log log y - 1 2 log π . (3.6)
A similar version of Theorem 3.1 was first proved by [START_REF] Chen | Inference in a simple change-point model[END_REF], via the strong invariance principle, in order to detect changes in the mean. Later, [START_REF] Steinebach | Heart rate variability. Standards of measurement, physiological interpretation, and clinical use[END_REF] have studied the same results for the detection of changes in the intensity of the renewal counting processes, and have improved them by giving another rate for the size of the sliding window under assumptions (A 2 ) and (A 3 ). The reader is also referred to the book of [START_REF] Révész | Random Walks in Random and non-random enviroments[END_REF] for a good summary about the increments on partial sums.

Furthermore, for a known change point τ k , the law of the random variable D 1 (τ k , A n ) defined by (3.1) is known. Therefore, the probability to detect a known change point in the mean is described in the following remark.

Remark 3.1 (Probability to detect known change point in the mean) Let τ k be a known change point in the mean of size

δ k := |µ τ k -µ τ k +1 |.
Then, under assumption (A 1 ), the probability to detect it is given by

P (|D 1 (τ k , A n )| ≥ C 1 ) = 1 -φ (C 1 -δ k ) √ A n √ 2σ + φ (-C 1 -δ k ) √ A n √ 2σ , (3.7)
and under assumptions (A 2 ) and (A 3 ), we have the following asymptotic probability

lim n→+∞ P A n |D 1 (τ k , A n ) -δ k | ≥ C 1 = 2φ -C 1 √ 2σ (3.8)
where C 1 is the threshold fixed in Step 1 and φ is the standard normal distribution. The proof can be deduced from

D 1 (τ k , A n ) = δ k + √ 2σ √ A n Z where Z is the standard normal r.v.
In applications, the variance σ 2 is unknown. For this reason we may replace it by its empirical estimator, σ 2 n . But, in order to keep the same result as in Theorem 3.1, the estimator σ 2 n has to verify a certain condition given by the following theorem.

Theorem 3.2 (Change point in the mean with unknown variance)

We apply to the same notations and the same assumptions as in Theorem 3.1. Moreover, we assume that σ n is an estimator of σ satisfying

lim n→+∞ |σ -σ n | log n P = 0, (3.9)
where the sign P = means convergence in probability. Then, under the null hypothesis,

lim n→+∞ P max k∈[An:n-An] |D 1 (k, A n )| ≤ σ n √ A n c n (x) = exp(-2e -x ).
(3.10)

Remark that condition (3.9) is not really restrictive, indeed as soon as σ 2 n satisfies a CLT, the condition is verified. For example, with the usual empirical variance estimator, a fourth order moment (of X) is sufficient.

Change in the variance

Now, we consider the case where we have a set of observations (X i ) i=1,...,n and we wish to know whether their variance has changed at an unknown time. If µ is known, then the problem is very simple. Testing (H 0 ) against (H 1 ) means that we are looking for a change in the mean of the sequence (X iµ) 2 i=1,...,n . Filtered Derivative method applied to the variance is based on the difference between the empirical variance computed on two sliding windows respectively at the right and at the left of the index k, both of size A n which satisfy condition (3.3). This difference is in fact a sequence of random variables denoted by (D 2 (k, A n )) An≤k≤n-An and defined as follows

D 2 (k, A n ) = σ2 (k, A n ) -σ2 (k -A n , A n ) (3.11) where σ2 (k, A n ) = 1 A n k+An j=k+1 (X j -µ) 2
denotes the empirical variance of X on the box [k + 1, k + A n ]. By using Theorem 3.1, we can deduce the asymptotic distribution of the maximum of |D 2 | under null hypothesis (H 0 ). This gives straightforwardly the following corollary:

Corollary 3.1 (Change point in the variance with known mean) Let (X i ) i=1,...,n be a sequence of independent identically distributed random variables with mean µ and assume that one of the following assumptions is satisfied

(A 4 ) ∃ t > 0 such as E[exp(t(X 1 -µ) 2 )] < +∞, and lim n→+∞ (log n) 3 A n = 0. (A 5 ) ∃ p > 2 such as E |X 1 -µ| 2p < +∞, and lim n→+∞ n 2/p log n A n = 0.
Let D 2 be defined by (3.11) and ν 2 = Var (X 1µ) 2 . Then, under the null hypothesis,

lim n→+∞ P max k∈[An:n-An] |D 2 (k, A n )| ≤ ν √ A n c n (x) = exp(-2e -x ), (3.12)
where c n (.) is defined by (3.5).

In practical situations we rarely know the value of the constant mean. However, μ(k, A n ) is a consistent estimator on the box [k + 1, k + A n ] for the mean µ under both null hypothesis (H 0 ) and alternative hypothesis (H 1 ). Thus, in the definition of the sequence D 2 , µ can be replaced by its estimator μ(k,

A n ) on the box [k + 1, k + A n ]. To be precise, let D 2 (k, A n ) = σ(k, A n ) -σ(k -A n , A n ), (3.13) where σ(k, A n ) = 1 A n k+An j=k+1 (X j -μ(k, A n )) 2
denotes the empirical variance of X on the sliding box [k + 1, k + A n ] with unknown mean. So, in order to obtain the same asymptotic distribution as the one obtained in Corollary 3.1, we must add extra conditions on the estimators μ(k, A n ). These new conditions are given in the next corollary.

Corollary 3.2 (Change point in the variance with unknown mean) With the notations and assumptions of Corollary 3.1, we suppose moreover that

lim n→+∞ max 0≤k≤n-An |µ -μ(k, A n )| (A n log n) 1 4 a.s = 0 (3.14)
where the sign a.s = means almost surely convergence. Then under the null hypothesis

lim n→+∞ P max k∈[An:n-An] | D 2 (k, A n )| ≤ ν √ A n c n (x) = exp(-2e -x ). (3.15)
Let us remark that (3.14) is not a very stringent condition.

Step 2: calculus of p-values for changes on the mean and changes on the variance

In this subsection, we recall p-value formula associated with the second test in order to remove false alarms. Let us stress that the only novelty of this subsection is the idea to divide the detection of abrupt change into two steps, see section 3. Since the r.v. X i are independent, the calculus of p-value relies on well known results that can be found in any statistical textbook. First, consider the Gaussian case and let us introduce some notations: For 1 ≤ k ≤ K max , let (X τ k-1 +1 , . . . , X τ k ) and (X τ k +1 , . . . , X τ k+1 ) two successive samples of i.i.d. Gaussian random variables such that

X τ k-1 +1 ∼ N (µ k , σ 2 k ) and X τ k +1 ∼ N (µ k+1 , σ 2 k+1 ).
We can use Fisher's F-statistic to determine the p-value of the existence of a change on the variance at time τ k under the null assumption (H0):

σ 2 k = σ 2 k+1 .
Then, we can use Student statistic to determine the p-value of a change on the mean at time τ k , that is under the null assumption (H0):

µ k = µ k+1 .
Second, consider the general case. Since A n → ∞ and by construction τ k+1τ k ≥ A n , we can apply CLT as soon as the r.v. X i satisfy Lindeberg condition. Thus, the empirical mean and the empirical variance converge to the corresponding ones in the Gaussian case. Eventually, if one of the assumptions (A i ) for i = 1, 2 or 3, then we can still apply Fisher and Student statistics to compute the p-value.

LINEAR REGRESSION

In this Section, we consider linear regression model:

Y i = aX i + b + e i , for 1 ≤ i ≤ n, (4.1)
where the terms (e i ) i=1,...,n are independent and identically distributed Gaussian random errors with zero-mean and variance σ 2 and, to simplify, (X i ) i=1,...,n are equidistant time points given by

X i = i∆ with ∆ > 0. (4.2)
Our aim is to detect change points on the parameters (a, b) of the linear model. As the Filtered Derivative is a local method, and to simplify notations, we will restrict ourselves here to one change point.

The two following subsections give the result for detection of potential change points on the slope and on the intercept. Subsection 4.3 provides formulas for calculating the p-values during Step 2.

Change in the slope

In this subsection, we are concerned with detection of change points in the slope a when the intercept remains constant.

Filtered Derivative method applied to the slope is based on the differences between estimated values of the slope a computed on two sliding windows at the right and at the left of the index k, both of size A n . These differences, for k ∈ [A n , n -A n ], form a sequence of random variables, given by

D 3 (k, A n ) = â(k, A n ) -â(k -A n , A n ) (4.3) where â(k, A n ) =   A k+A j=k+1 X j Y j - k+A j=k+1 X j k+A j=k+1 Y j     A k+A j=k+1 X 2 j -   k+A j=k+1 X j   2   -1 , (4.4)
is the estimator of the slope a on the (sliding

) box [k + 1, k + A n ].
Let us stress that these quantities can be calculated by recurrence with complexity O(n).

Our first result gives the asymptotic distribution of the maximum of |D 3 | under the null hypothesis of no change on the linear regression.

Theorem 4.1 (Change point in the slope) Let (X i ) i=1,...,n and (Y i ) i=1,...,n be given by (4.2) and (4.1) where e i is a family of i.i.d. mean zero Gaussian r.v. with variance σ 2 . Let D 3 be defined by (4.3) and assume that A n satisfies condition (3.3). Then under the null hypothesis

lim n→+∞ P max k∈[An:n-An] |D 3 (k, A n )| ≤ 2 √ 6σ ∆ A n (A 2 n -1) d n (x) = exp(-2e -x ), (4.5) with d n (x) = c (A n , x
) and c(., .) is defined by (3.6).

Moreover, for a given change point τ k , the law of the random variable D 3 (τ k , A n ) defined by (4.3) is known. Therefore, the probability to detect a known change point in the slope is specified in the following remark.

P (|D 3 (τ k , A n )| ≥ C 1 ) = 1 -φ (C 1 -δ a,k )∆ A n (A 2 n -1) 2 √ 6σ +φ (-C 1 -δ a,k )∆ A n (A 2 n -1) 2 √ 6σ (4.6)
where C 1 is the threshold fixed in Step 1 and φ is the standard normal distribution. The proof

is directly deduced from D 3 (τ k , A n ) = δ a,k + 2 √ 6σ ∆ A n (A 2 n -1)
Z where Z is a standard normal random variable.

Change in the intercept

In this subsection, we focus on the detection of change points in the intercept with known slope a. To do this, we calculate the differences between estimators of the intercept b computed on two sliding windows respectively at the right and at the left of the index k, both of size A n .

For A n ≤ k ≤ n -A n , these differences form a sequence of random variables given by

D 4 (k, A n ) = b(k, A n ) -b(k -A n , A n ) (4.7) where b(k, A n ) = 1 A k+A j=k+1 Y j -a × 1 A k+A j=k+1 X j , ( 4
|D 4 (k, A n )| ≤ σ √ A n c n (x) = exp(-2e -x ), (4.9)
where c n (.) is defined by (3.5).

In addition, for a known change point τ k , the distribution of the random variable D 4 (τ k , A n ) defined by (4.7) is known. Therefore, the probability to detect a known change point in the intercept is given in the following remark.

P (|D 4 (τ k , A n )| ≥ C 1 ) = 1 -φ (C 1 -δ b,k ) √ A n √ 2σ + φ (-C 1 -δ b,k ) √ A n √ 2σ (4.10)
where C 1 is the threshold fixed in Step 1 and φ is the standard normal distribution. The proof may be easily deduced from

D 4 (τ k , A n ) = δ b,k + √ 2σ √ A n Z where Z is a standard normal random variable.
In real applications the slope is often unknown. In this case, we replace it in (4.7) by its empirical estimator a n . This leads to the definition

D 4 (k, A n ) = b(k, A n ) -b(k -A n , A n ) (4.11) where b(k, A n ) = 1 A k+A j=k+1 Y j -a n × 1 A k+A j=k+1 X j ,
is the estimator of the intercept b on the (sliding) box [k + 1, k + A n ] with unknown slope. Naturally, we must assume that the estimator of the slope satisfy a certain convergence condition which is given in the following corollary.

| D 4 (k, A n )| ≤ σ √ A n c n (x) = exp(-2e -x ). (4.13)

Step 2: calculus of p-values

Let us give here p-value formulae associated to Step 2 in linear regression model (4.2) and (4.1). More precisely, we are concerned with detection of right change points on the slope or intercept. We recall that Step 2 of FDp-V method has been introduced in order to eliminate false alarms. Indeed, at Step 1, a time X τ k has been selected as a potential change point.

At

Step 2, we test whether the slopes and intercepts of two data sets at left and right of the potential change point X τ k are significantly different or not, and we measure this by the corresponding p-value. Before going further, let us introduce some notations. For 1

≤ k ≤ K max , let R k = (X τ k-1 +1 , Y τ k-1 +1 ), . . . , (X τ k , Y τ k ) and R k+1 = (X τ k +1 , Y τ k +1 ), . . . , (X τ k+1 , Y τ k+1
) two successive samples of observations such that the relationship between variables X and Y is given by (4.2) and (4.1), or more explicitly by

Y j = a k X j + b k + e j for τ k-1 + 1 ≤ j ≤ τ k .
By using definition of the error terms,we can easily see that slope estimator a k and intercept estimator b k have Gaussian distribution given by

a k ∼ N a k , σ 2 a k with σ 2 a k σ 2 =   τ k j= τ k-1 +1 (X j -X k ) 2   -1 , b k ∼ N b k , σ 2 b k with σ 2 b k σ 2 = 1 n k + X k   τ k j= τ k-1 +1 (X j -X k ) 2   -1
, respectively, where X k is the empirical mean of the sequence X τ k-1 +1 , . . . , X τ k . In the sequel, we denote by σ 2 a k and σ 2 b k the empirical variance of respectively the random variables a k and b k .

Comparing slope

We want to test if the samples R k and R k+1 present a change in slope or not.

H a 0,k : a k = a k+1 against H a 1,k : a k = a k+1 .
Then, the p-value p k,a associated to the potential change point τ k in order to eliminate false alarms for changes in slope is given by

p k,a = 1 -φ a     | a k -a k+1 | σ 2 a k n k + σ 2 a k+1 n k+1    
where φ a is a Student T-distribution with υ a degrees of freedom

υ a = N σ 2 a k , σ 2 a k+1 , n k , n k+1 =        σ 2 a k n k + σ 2 a k+1 n k+1 2 σ 2 a k n k √ n k -1 2 + σ 2 a k+1 n k+1 √ n k+1 -1 2       
and ⌊x⌋ is the integer part of x.

Comparing intercept

We desire to study if the intercept before and after a detected change points τk are equal or not, i.e. if it is a real change point, resulting in testing

(H b 0,k ) b k = b k+1 against (H b 1,k ) b k = b k+1 .
Then, the p-value p k,b associated to the potential change point τ k in order to remove false alarms from the list of changes in intercept found in step 1, is given by

p k,b = 1 -φ b     | b k -b k+1 | σ 2 b k n k + σ 2 b k+1 n k+1     where φ b is a Student T-distribution with υ b = N σ 2 b k , σ 2 b k+1 , n k , n k+1 degrees of freedom.

NUMERICAL RESULTS

In this section, we apply FDpV and PLSC methods to detect abrupt changes in the mean of simulated Gaussian r.v and we compare the results obtained by both methods. Next, we use the FDpV for the detection of change points in the slope of linear regression model. Finally, we run the FDpV algorithm to the segmentation of heartbeat time series, and average daily volume drawn from the financial market.

A toy model: off-line detection of abrupt changes in the mean of independent Gaussian random variables with known variance

In the following subsection, we consider the elementary problem, namely: The off-line detection of multiple change points in the mean of simulated independent Gaussian r.v with known variance, and we numerically compare the efficiency of the different estimators given by the FDpV and the PLSC procedures.

Numerical simulation

To begin, for n = 5000 we have simulated a sequence of Gaussian r.v (X 1 , . . . , X n ) with known variance σ 2 = 1 and mean µ i = g(i/n) where g is a piecewise-constant function with five change points, i.e. we have chosen a configuration 1 = τ

0 < τ 1 < • • • < τ K < τ K+1 = n with K = 5

Monte-Carlo simulation

In this paragraph, we have made M = 1000 simulations of independent copies of sequences of Gaussian r.v.

X (k) 1 , . . . , X (k)
n with variance σ 2 = 1 and mean µ(i) = g(i/n), for k = 1, . . . , M . On each sample, we apply the FDp-V algorithm and the PLSC method. We find the right number of changes in 98.1% of all cases for the first method and in 97.9% for the second one, see Figure 5.1. Then, we compute the mean square errors. There are two kinds of mean square errors:

• Mean Integrate Square Error: E 1 n n i=1 | g(i/n) -g(i/n)| 2 = E g -g 2 L 2 ([0,1])
.The estimated function is obtained in two steps: first we estimate the configuration of change points (τ k ) k=1,..., K , then we estimate the value of g between two successive change points as the empirical mean.

• Square Error on Change Points: Table 1 gives the result of Monte Carlo simulation mean errors, and also the comparison between the mean time complexity and the mean memory complexity. We have written the two programs in Matlab and have runned it with computer system which has the following characteristics: 1.8GHz processor and 512MB memory.

E K k=1 | τ k -τ k | 2 ,

Square Error on Change Points

Mean Integrated Squared Error FDp-V method 1.1840 × 10 -4 0.0107 PLSC method 1.2947 × 10 -4 0.0114 Memory allocation (in Megabytes) CPU time (in second) FDp-V method 0.04 MB 0.005 s PLSC method 200 MB 240 s Table 5.1: Errors and complexities given by FDp-V method and PLSC method.

Numerical conclusion

On the one hand, both methods have the same accuracy in terms of percentage of the right number of changes, and in terms of Square Error or Mean Integrate Square Error. On the other hand, the Filtered Derivative with p-Value is less expensive in terms of time complexity and memory complexity, see

Off-line detection of changes in the slope of simple linear regression

In this subsection, we consider the problem of multiple change points detection in the slope of linear model corrupted by an additive Gaussian noise. At first, for n = 1400 we have simulated the sequences X = (X 1 , . . . , X n ) and Y = (Y 1 , . . . , Y n ) defined by (4.2) and (4.1) with ∆ = 1, σ = 30 and a i = h(i) where h is a piecewise-constant function with four change points, i.e. we have chosen a configuration 1 = τ

0 < τ 1 < • • • < τ K < τ K+1 = n with K = 4

Application to real data

In this subsection, we apply our algorithm to detect change points in the mean of two real samples: the first one is concerned with health and wellbeing, and the second one with finance. Our main purpose here is to show that our estimation procedure is sufficiently robust to consider non-independent time series, that we will study theoretically in a subsequent paper. We moreover provide an analysis of the obtained results.

An Application to Change Point Detection of Heartbeat Time Series

In this paragraph, we give examples of application of FDp-V method to heartbeat time series analysis. Electrocardiogram (ECG) has been processed from a long time since the implementation of monitoring by Holter in the fifties. We consider here the RR interval, which provides an accurate measure of the length of each single heartbeat and corresponds to the instantaneous speed of the heart engine, see Task force of the European Soc. Cardiology and the North American Society of Pacing and Electrophysiology (1996). From the beginning of 21st century, the size reduction of the measurement devices allows to record heartbeat time series for healthy people in ecological situations throughout long periods of time: Marathon runners, individuals daily (24 hours) records, etc. We then obtain large data sets of more than 40.000 observations, resp. 100,000 observations, that address change detection of heart rate. In

$) **!%$ %$ (* ('+$-% (!) (* %$ !# ()%"+*!%$ #!$+*) )%$) **!%$ * () %" *) -#!$+* $+#( % $ &%!$*) !) & ,"* $(- & ,"* $(- "$ %(# (* Figure 5
.9: Segmentation of heartbeat time series of a marathon runner. Blue: Heart rate segmented by FDpV. Red: Low Frequency wavelet energy segmented by FDpV. Green: High Frequency wavelet energy segmented by FDpV ondly, in Figure 5.9, we notice beginning and end of the marathon race, but also training before the race, and small breaks during the race. The same FDpV compression technology has been applied to wavelet energy corresponding to High Frequency, resp. Low Frequency band. Energy into these two frequencies bands is interpreted by cardiologists as corresponding to heart rate regulation by sympathetic and ortho-sympathetic systems. We refer to [START_REF] Ayache | Discretization error of wavelet coefficient for fractal like process[END_REF] and [START_REF] Khalfa | Heart rate regulation processed through wavelet analysis and change detection[END_REF] for detailed explanations. In this paper, we just want to show that FDpV detects changes on HF and LF energy at the beginning of the race. Fig. 5.10 is concerned with heart beats of a shift worker Y1 through out a day its the life. The shift worker Y1 has manually reported changes of activity on a diary, as shown in plotted the segmentation resulting from manually recorded diary. In blue, we have plotted the automatic segmentation resulting from FDpV method. Note that the computation time is 15 seconds for 120,000 data with a code written in Matlab in a 2.8 GHz processor. Moreover, the FDpV segmentation is more accurate than the manual one. For instance, Y1 has reported "Football training" from 20h to 21h38, which is supported by our analysis. But with FDpV method, we can see more details of the training namely the warming-up, the time for coach's recommendations, and the football game with two small breaks.

Let us remark once again that both heartbeat time series and wavelet coefficient series do not fulfill the assumption of independency. However, FDpV segmentation provides accurate information on the heart rate, letting suppose a certain robustness of the method with respect to the stochastic model.

In this case, the FD-pV method has the advantage of being a fast and automatic procedure of segmentation of a large dataset, on the mean in this example, but possibly on the slope or the variance. Combined with the "tachogram cleaning" algorithm, we then have an entirely automatic procedure to obtain apparently homogeneous segment. The next step will be to detect change on hidden structural parameters. First attempts in this direction are exposed in [START_REF] Ayache | Discretization error of wavelet coefficient for fractal like process[END_REF] and [START_REF] Khalfa | Heart rate regulation processed through wavelet analysis and change detection[END_REF], but should be completed by forthcoming studies.

Changes in the average daily volume

Trading volume represents number of shares or contracts traded in a financial market during a specific period. Average traded volume is an important indicator in technical analysis as it is used to measure the worth of a market move. If the markets move significantly up or down, the perceived strength of that move depends on the volume of trading in that period. The higher the volume during that price move, the more significant the is move. Therefore, the detection of abrupt changes in the average daily volume provide relevant information for financial engineer, trader, etc. Then, we consider here a daily volume of Carbone Lorraine compagny observed during 02 January 2009. These data have been kindly given by Charles-Albert Lehalle from Crédit Agricole Cheuvreux, Groupe CALYON (Paris). The results obtained with our algorithm for A = 300, p * 1 = 0.05 and p * 2 = 10 -5 are illustrated in Figure 5.11. It appears that FDp-V procedure detects major changes observed after each huge variations. In future works, we will investigate sequential detection of change points in the average daily volume in connection with the worth of a market move.

Remark 5.1

We note that in our procedure, the observations are assumed to be independent. However, by Durbin and Watson test, we show that this condition is not checked, and we obtain a global correlation ρ Heartbeat = 0.9982 for the heartbeat data, and ρ volume = 0.4663 for the financial time series. But this is not a restriction, on the contrary, this shows the robustness of our algorithm. In future works, this algorithm may be adapted to dependent data.

CONCLUSION

It appears that both methods, namely FDp-V and PLSC, give right results with practically the same precision. But, when we compare the complexity, we remark that the FDp-V method is less expensive in terms of time and memory complexity. Consequently, FDp-V method is faster (time) and cheaper (memory), and so it is more adapted to segment random signals with large or huge datasets.

In future works, we will develop the Filtered Derivative with p-value method in order to detect abrupt changes in parameters of weakly or strongly dependent time series. In particularly, we will consider the detection problem on the Hurst parameter of multifractional Brownian motion and apply it to physiological data as in [START_REF] Billat | Detection of changes in the fractal scaling of heart rate and speed in a marathon race[END_REF]. Let us also mention that the FDp-V method is based on sliding window and could be adapted to sequential detection, see for instance [START_REF] Bertrand | Detecting small shift on the mean by finite moving average[END_REF]; [START_REF] Bertrand | Filtered derivative with p-value method for multiple change-points detection[END_REF].

PROOF

Proof of Theorem 3.1. Under the null hypothesis (H 0 ) the filtered derivative can be rewritten as follows

√ A n σ D 1 (k, A n ) = S k-An -2S k + S k+An √ A n ,
where

S k = k j=1 ξ j , with (ξ j ) j∈[1,n] a sequence of i.i.d r.v such as E[ξ 1 ] = 0, E[ξ 2 1 ] = 1, S 0 = 0.
To achieve our goal, we state three lemmas. First, we show in Lemma 6.1 that if a positive sequence (a n ) has the following asymptotic distribution

lim n→+∞ P (a n ≤ c n (x)) = exp(-2e -x ),
where c n (.) is defined by (3.5) and if there is a second positive sequence (b n ) which converges almost surely (a.s) to (a n ) with rate of convergence of order O √ log n , then (b n ) has the same asymptotic distribution as (a n ). We then prove in Lemma 6.2 that under one of the assumptions (A i ) with i ∈ {1, 2, 3}, the maximum of the increment

√ A n σ |D 1 (k, A n )| converges
a.s. to the maximum of discrete Wiener process' increment with rate O √ log n . We show in Lemma 6.3 that the maximum of discrete Wiener process' increment converges a.s. to the maximum of continuous Wiener process' increment with rate O √ log n . Then, by applying Qualls and Watanabe (1972, Theorem 5.2, p. 594), we deduce the asymptotic distribution of the maximum of continuous Wiener process' increment. Finally, by combining these results, we get directly (3.4). To begin with, let us state the first lemma. We denote by |∆x| an infinitesimally small change in x. Then, for n large enough and |∆x| small enough, c n (x + |∆x|) and c n (x -|∆x|) satisfy the following inequalities

c n (x + |∆x n |) ≥ c n (x) + |∆x n | √ 2 log n , (6.2) c n (x -|∆x n |) ≤ c n (x) - |∆x n | √ 2 log n . (6.3)
Following [START_REF] Chen | Inference in a simple change-point model[END_REF], we supply a lower and an upper bounds of P (b n ≤ c n (x)). On the one hand, by using the inequality (6.2), the upper bound results from the following calculations

P (b n ≤ c n (x)) ≤ P b n ≤ c n (x + |∆x|) - |∆x| √ 2 log n ≤ P a n -η n ≤ c n (x + |∆x|) - |∆x| √ 2 log n ≤ P (a n ≤ c n (x + |∆x|)) + P -η n ≤ - |∆x| √ 2 log n = P (a n ≤ c n (x + |∆x|)) + P η n ≥ |∆x| √ 2 log n .
On the other hand, by using the inequality (6.3), the lower bound results from analogous calculations

P (a n ≤ c n (x -|∆x n |)) ≤ P (b n ≤ c n (x)) + P η n ≥ |∆x n | √ 2 log n .
Then, by putting together the two previous bounds, we obtain

P (a n ≤ c n (x -|∆x n |)) -P η n ≥ |∆x n | √ 2 log n ≤ P (b n ≤ c n (x)) ≤ P (a n ≤ c n (x + |∆x n |)) + P η n ≥ |∆x n | √ 2 log n .
Finally, by taking the limit in n and as |∆x| is arbitrary small, we deduce (6.1). This finishes the proof of Lemma 6.1. Now, we show that the maximum of

√ An σ |D 1 (k, A n )| converges a.
s to the maximum of discrete Wiener process' increment with rate of convergence of order O √ log n . This is stated in Lemma 6.2 below (which gives also corrections to previous results of Chen (1988): Lemma 6.2 Let (W t , t ≥ 0) be a standard Wiener process and Z An (q), 0 ≤ q ≤ n A n -1 be the discrete sequence defined by

Z An (q) = W q-1 -2W q + W q+1 if 1 ≤ q ≤ n An -1, 0 else . (6.4)
Let (X i ) i=1,...,n be a sequence of independent identically distributed random variables with mean µ, variance σ 2 , D 1 be defined by (3.1), and denote

η 1,n = max 0≤k≤n-An √ A n σ |D 1 (k, A n )| -max 0≤q≤ n An -1 |Z An (q)| .
Moreover, we suppose that one of the assumptions (A i ), with i ∈ {1, 2, 3} is in force. Then there exists a Wiener process (W t ) t such that

lim n→+∞ η 1,n log n a.s = 0. (6.5) ♦ Proof of Lemma 6.2. We consider a new discrete sequence, (B(k, A n ), 0 ≤ k ≤ n -A n ),
obtained by scaling from the sequence Z An (q), 0 ≤ q ≤ n An -1 . It is defined as follows

B(k, A n ) =    W k-An -2W k + W k+An √ A n if A n ≤ k ≤ n -A n 0 else Then η 1,n D = max 0≤k≤n-An √ A n σ |D(k, A n )| -max 0≤k≤n-An |B(k, A n )| ,
where the sign D = means equality in law. Depending on which assumption (A i ) is in force, we have three different proofs:

1. Assuming (A 1 ). This is the simplest case. We can choose a standard Wiener process, (W t , t ≥ 0), such that

S k = W k at all the integers k ∈ [0, n -A n ]. Hence, √ An σ D(k, A n ) = B(k, A n ).
Then, we can deduce (6.5).

Assuming (A 2 ).

We have

0 ≤ η 1,n ≤ max 0≤k≤n-An |D(k, A n ) -B(k, A n )| ≤ 4 √ A n max 0≤k≤n-An |S k -W k | and after 0 ≤ η 1,n log n ≤ 4(log n) 3 2 √ A n × max 0≤k≤n-An |S k -W k | log n
However, lim n→+∞ (log n) 3 An = 0 and according to Komlós et al. (1975, Theorem 3, p.34), there is a Wiener process, (W t , t ≥ 0), such as lim

n→+∞ max 0≤k≤n |S k -W k | log n
< +∞. Then, we can deduce (6.5).

3. Assuming (A 3 ). By using the same tricks than in the proof of the case (A 2 ), we can show that 0 ≤ η 1,n log n ≤ 4n

1 p √ log n √ A × max 0≤k≤n-An |S k -W k | n 1 p
However, lim n→+∞ n 2 p log n A n = 0 and according to Komlós et al. (1975, Theorem 3, p.34), there is a Brownian motion, (W t , t ≥ 0), such as lim

n→+∞ max 0≤k≤n |S k -W k | n 1 p < +∞. Then,
we can deduce (6.5).

This finishes the proof of Lemma 6.2.

In order to apply Qualls and Watanabe (1972, Theorem 5.2, p. 594) theorem, we have to consider a continuous version of the process Z An . For this reason, we define the continuous

process Z An (t), t ∈ 0, n A n -1 such as Z An (t) = W t-1 -2W t + W t+1 . (6.6)
Then, in Lemma 6.3, we show that the maximum of |Z An (q)| converges a.s to the maximum of |Z An (t)| with rate of convergence of order O √ log n .

Lemma 6.3

Let Z An be defined by (6.4) and set η 2,n = sup t∈ 0, n An -1

|Z An (t)| -max 0≤q≤ n An -1 |Z An (q)| .
Then lim n→+∞ η 2,n log n a.s = 0. (6.7)

♦

Proof of Lemma 6.3. We have

0 ≤ η 2,n ≤ 4 sup t∈ 0, n An -1 |W t+1 -W t | ≤ 4 sup t∈ 0, n An -1 An sup s∈ 0, 1 An |W t+s -W t | . This implies 0 ≤ η 2,n log n ≤ 4 log n sup t∈ 0, n An -1 An sup s∈ 0, 1 An |W t+s -W t |
and after

P η 2,n log n ≥ δ ≤ P    sup t∈ 0, n An -1 An sup s∈ 0, 1 An |W t+s -W t | ≥ δ √ A n 4 √ log n 1 √ A n   
According to Csörgö and Révész (1981, Lemma 1.2.1, p. 29), and by taking T

= n/A n , h = 1/A n , ε = 1, ν = δ √ A n 4 √ log n
, and C a non negative real, we deduce that

P η 2,n log n ≥ δ ≤ C n A n exp - δ 2 A n 48 log n .
Next, by using lim This result can be proved by applying Theorem 5.2 of [START_REF] Qualls | Asymptotic properties of gaussian processes[END_REF] to the centered stationary Gaussian process (Z An (t), t ≥ 0). The covariance function of Z An (t) is given by

ρ(τ ) = Cov (Z An (t), Z An (t + τ )) V ar (Z An (t)) V ar (Z An (t + τ )) =      1 -|τ | if 0 ≤ |τ | ≤ 1, -1 + |τ | 2 if 1 < |τ | ≤ 2, 0 if 2 < |τ | < +∞.
So, the conditions of Theorem 5.2 of [START_REF] Qualls | Asymptotic properties of gaussian processes[END_REF] are satisfied in the following way: α = 1, H α = 1 (according to Pickands (1969, p. 77), and σ -1 (x) = 2x -2 . This finishes the proof of (6.8). Eventually, by combining Lemma 6.1, Lemma 6.3 and result (6.8), we get the asymptotic distribution of the maximum of the sequence |Z An (q)|, 0 ≤ q ≤ n A n -1 . Then, by using Lemma 6.1 and Lemma 6.2, we immediately get the distribution of the maximum of the filtered derivative sequence (|D

1 (A n , k)|, 0 ≤ k ≤ n -A n ).
End of the proof of Theorem 3.1.

Proof of Theorem 3.2. Fix ε > 0, the key argument is to divide Ω into two complemen- Eventually, combined with (6.9) this implies (3.10). To finish the proof, it just remains to verify that (6.10) is satisfied. Indeed, Then, by using (6.13), we obtain (4.5). This finishes the proof of Theorem 4.1.

Proof of Corollary 4.1. First we note that, under the null hypothesis (H 0 ), the Filtered Derivative applied to the intercept satisfy

D 4 (k, A n ) = A -1 n k+An j=k+1 e j -A -1 n k j=k-An+1 e j .
Therefore, (D 4 (k, A n )) An≤k≤n-An corresponds to a sequence of Filtered Derivative of the mean, applied to the particular case of i.i.d centered Gaussian r.v with known variance σ 2 . Then, by applying Theorem 3.1 under assumption A 1 , we obtain (4.9). This finishes the proof of Corollary 4.1.

Proof of Corollary 4.2. Let D 4 and D 4 be defined respectively by (4.7) and (4.11), and set

η 5,n = max 0≤k≤n-An √ A n σ |D 4 (k, A n )| -max 0≤k≤n-An √ A n σ D 4 (k, A n ) .
The key argument is to prove that lim n→+∞ η 5,n log n a.s = 0. (6.15) By using assumption (4.12). Then, we apply Lemma 6.1 which combined with Corollary 4.1 implies (4.13). So, to finish the proof we must verify (6.15). Next, we have

Figure 2

 2 Figure 2.1: Above: Gaussian random variables simulated with constant variance and presenting changes in the mean. Below: Filtered derivative function D.

Figure 2 .

 2 Figure 2.2 below provides an example: The family of potential change points contains the right change points (surrounded in blue in Figure 2.2) as well as false alarms (surrounded in green in Figure 2.2).

Figure 2 . 2 :

 22 Figure 2.2: Detection of potential change points

Figure 2 . 3 :

 23 Figure 2.3: Above: Detection of right change points. Below: Theoretical value of the piecewiseconstant map t → µ t (black), and its estimators given by PLSC method (blue) and FDp-V method (red).

  Remark 4.1 (Probability to detect known change point in the slope) Let τ k be a known change point in the slope of size δ a,k := |a τ ka τ k +1 |. Then the probability to detect it is given by

  .8) is the estimator of the intercept b on the (sliding) box [k+1, k+A n ]. By applying Theorem 3.1, we get the asymptotic distribution of the maximum of |D 4 | under the null hypothesis of no change on the linear regression.Corollary 4.1 (Change point in the intercept) With the notations and assumptions of Theorem 4.1, we have under the null hypothesis lim n→+∞ P max k∈[An:n-An]

  Remark 4.2 (Probability to detect known change point in the intercept) Let τ k be a known change point in the intercept of size δ b,k := |b τ kb τ k +1 |. Then the probability to detect it is given by

Corollary 4. 2 (

 2 Change point in the intercept with unknown slope) Under the same notations and the same assumptions than in Corollary 4.1. Moreover, we assume that the estimator ân of a satisfies the following condition

  and such as µ i = g(k) for τ k ≤ i < τ k+1 . The size of the change point at τ k specified by δ k := |µ τ kµ τ k +1 | are to be chosen in the intervalle [0.5, 1.25]. Then, we have computed the function k → |D 1 (A, k)| with A = 300, see Figure 2.1. Both methods, namely Filtered Derivative with p-values, p * 1 = 0.05 and p * 2 = 10 -4 , and Penalized Least Squares Criterion provide right results, see Figure 2.3. This example is plainly confirmed by Monte-Carlo simulations.

Figure 5 . 1 :

 51 Figure 5.1: Distribution of the estimated number of change points K for M = 1000 realizations. Left: Using PLSC method. Right: Using Filtered Derivative method.

  in the case where we have found the right number of change points.

Figure 5 . 3 :Figure 5 . 4 :

 5354 Figure 5.2: Case 1: Scatter plot of the simulated data (X i , Y i ) for 1 ≤ i ≤ n.

Figure 5 . 6 :Figure 5 . 7 :

 5657 Figure 5.5: Case 2: Scatter plot of the simulated data (X i , Y i ) for 1 ≤ i ≤ n.

Figure 5 . 8 :

 58 Figure 5.8: Raw heart rate time series of a marathon runner

Figure 5 . 10 :

 510 Figure 5.10: Segmentation of heartbeat time series of shift worker Y1. Yellow: cleaned heart rate. Blue: heart rate segmented by FDpV. Red: Heart rate segmented following the dairy.

  Figure 5.11: Segmentation of the daily volume of Carbone Lorraine compagny observed during 02 January 2009.

  Let (a n ) and (b n ) two sequences of positive random variables and we denote η n = |a nb n |. We assume that 1. lim n→+∞ P (a n ≤ c n (x)) = exp(-2e -x ).where c n (.) is defined by (3.5). Thenlim n→+∞ P (b n ≤ c n (x)) = exp(-2e -x ).(6.1)♦Proof of Lemma 6.1. Without any restriction, we can consider the case where η n > 0.

  n = {|σσ n | log n ≤ σε} and Ω 2,n = {|σσ n | log n > σε} . k∈[An:n-An] |D 1 (k, A n )| ≤ σ n √ A n c n (x) and Ω 2,nOn the one hand, we remark thatP max k∈[An:n-An] |D 1 (k, A n )| ≤ σ n √ A n c n (x) and Ω 2,n ≤ P (Ω 2,n ) , hand, for all ω ∈ Ω 1,n , we have σ n = σ(1 + λ n (ω)) with |λ n (ωk, A n )|, we get P max k∈[An:n-An] |D 1 (k, A n )| ≤ σ n √ A n c n (x) and Ω 1,n = P (a n ≤ c n (x)(1 + λ n (ω))) = P (a nη 3,n ≤ c n (x))with η 3,n = c n (x)λ n . Therefore, after having checked Lemma 6.1 which combined with Theorem 3.1 implies thatP max k∈[An:n-An] |D 1 (k, A n )| ≤ σ n √ A n c n (x)and Ω 1,n = exp(-2e -x )

  Proof of Corollary 3.2. Let D 2 and D 2 be defined respectively by(3.11) and (3(3.14). Then, we apply Lemma 6.1 which combined with Corollary 3.1 implies (3.15). So, to finish the proof we must verify (6|µµ k | 2 -|µµ k-An | 2 and after 0 ≤ η 4,n log n ≤ 2 ν max 0≤k≤n-An |µµ k | 2 A n log n.Therefore, by using condition (3.14), we get lim n→+∞ η 4,n log n a.s = 0. This finishes the proof of Corollary 3.2. Hence, by applying Csáki and Gonchigdanzan (2002, Theorem 2.1, p. 3) to the sequence D Std 3 (k, A n ) An≤k≤n-An , we can deduce that lim n→+∞ P max k∈[An:n-An] |D Std 3 (k, A n )| ≤ d n (x) = exp(-2e -x ). (6.14)

  k, A n ) -D 4 (k, A n ) n log n ≤ |aân |A 3 2 n ∆ n log n.Therefore, by using condition (4.12), we get lim n→+∞ η 5,n log n a.s = 0. This finishes the proof of Corollary 4.2

Table 5

 5

.1. Indeed, Penalized Least Squares Criterion algorithm needs 200 Megabytes of computer memory, while Filtered derivative method only needs 0.008%. This plainly confirms the difference of time and memory complexity, i.e. O(n 2 ) versus O(n).

Table 5 .

 5 2. In yellow, we have plotted the cleaned heart rate time series. In red, we have

	Time	7h24-8h23 8h23-11h05 11h05-12h53	12h53-20h	20h-21h38	22h30-4h30
	Activity	Task 1	Task 2	Picking	Free afternoon playing football	Sleeping

Table 5

 5

.2: A shift worker's activities record

  Cantelli lemma, we can deduce (6.7). This finishes the proof of Lemma 6.3. Next, we applyQualls and Watanabe (1972, Theorem 5.2, p. 594) to the continuous process |Z An (t)| . We obtain the asymptotic distribution of its maximum given by

	n→+∞	(log n) 2 An	= 0, we can deduce
				n≥1	n A n	exp -	δ 2 A n 48 log n	< +∞,
	and after
				n≥1	P η 2,n log n ≥ δ < +∞.
	Then, according to Borel-lim n→+∞	P	  sup t∈[0, n A -1]	 |Z An (t)| ≤ c n (x)  = exp(-2e -x ).	(6.8)

Acknowledgments. We are grateful to the editors and the referees for their helpful comments. A first version of this work was presented during the "International Workshop in Sequential Methodologies" at UTT (Troyes, France, June 15-17, 2009). We would like to thank both organizers and participants of this conference for stimulating discussions and suggestions.

PROOF FOR LINEAR REGRESSION

Proof of Theorem 4.1. First we note that, under the null hypothesis (H 0 ), the sequence

are respectively the empirical mean and the empirical variance of X on the (sliding

. By using the definition (4.2), we see that

Therefore, we can deduce

(6.12)

Remark that the mean and the variance of the Gaussian sequence (D 3 (k, A n )) An≤k≤n-An verify

.

Moreover the variance does not depend on k. Then, (D 3 (k, A n )) An≤k≤n-An is a centered stationary Gaussian sequence. Next, theorem 4.1 becomes an application of Csáki and Gonchigdanzan (2002, Theorem 2.1, p. 3) which gives the asymptotic distribution of the maximum of standardized stationary Gaussian sequences (Z k ) k≥1 with covariance r n = cov (Z 1 , Z n+1 ) under the condition r n log(n) → 0. Let us define the standardized version of D 3 as (6.13) and its covariance

The following lemma provides the value of the covariance:

An≤k≤n-An the standardized stationary Gaussian sequence defined by (6.13). Then, its covariance matrix denoted (Γ(k 1 , k 2 )) 1≤k 1 ,k 2 ≤n is given by

♦

Proof of Lemma 6.4. First we note that, by symmetry property of covariance matrix, we can restrict ourselves to the case k 2k 1 ≥ 0. Next, let us distinguish three different expressions of the covariance according to the value of k 2k 1 .

By replacing jk 1 with u, we get

But, in order to use formula (6.12) we must find the case where the sign of u and u + k 2k 1 remain constant where u ∈

That is why we must distinguish the following subcases

This finishes the proof of Lemma 6.4.


Qiao Chen
email: qiao.chen@femto-st.fr

Yassine Haddab
email: yassine.haddab@femto-st.fr

Philippe Lutz
email: philippe.lutz@femto-st.fr

Microfabricated Bistable Module For Digital Microrobotics

High precision microrobots are needed more and more to perform micro/nanomanipulation and microassembly tasks in various environments like microrobotic stations, electronic microscopes (SEM, TEM), etc. Current microrobots are based on the use of smart materials to perform proportional or incremental actuation. To avoid the main drawbacks of these microrobots (non linearities, integration of sensors, robust control, energy consumption, sensitivity to noise), we propose a new type of microrobots, called digital microrobots, based on microfabricated bistable modules. The study presented in this paper is dedicated to the microfabricated bistable modules, notably the structure and the actuators design and characterization. The results open a new paradigm in the field of microrobotics leading to open loop control and the design of various kinematics adapted to the microworld. Moreover, no external energy is required to maintain the microrobot in its position.

Microrobotics, Bistable mechanism, Thermal actuator, Open-loop control of robot.

INTRODUCTION

During the last decade, significant research activities have been performed in the field of microrobotics. Microrobotics deals with the design, the fabrication and the control of microrobots.

These microrobots are intended to perform various tasks in the so-called Microworld (i.e. the world of submillimetric objects), in particular micromanipulation of single objects (artificial or biological) for positioning, characterizing or sorting as well as for industrial microassembly [START_REF] Lopez-Walle | Principle of a submerged freeze gripper for microassembly[END_REF].

Research already done has shown that the use of active materials to actuate microrobots gives better performances than the use of traditional actuators. However, despite their intrinsic high resolution, these materials have some drawbacks, making the design of efficient controllers a hard task. Their behavior is often complex, nonlinear and sometimes non stationary (environmental conditions as temperature and humidity have a great influence on static and dynamic behaviors [START_REF] Rakotondrabe | Hysteresis and vibration compensation in a nonlinear unimorph piezocantilever[END_REF]). Closed-loop control of the microrobots requires the integration of very small sensors and the use of bulky and expensive instruments for signal processing and real-time operating. Packaging and integration of sensors and actuators are also hard problems to resolve. This is why building multidegrees of freedom microrobots able to perform complex tasks is difficult [START_REF] Bamberger | A novel six degrees-of-freedom parallel Robot for MEMS fabrication[END_REF], [START_REF] Rakotondrabe | Development, modelling and control of micro/nano positionning 2 dof stick-slip device[END_REF].

Generally, two actuation modes are used in the control of microrobots to generate movements: the proportional mode and the incremental mode. The proportional actuation is based on the controlled deformation of an active material that directly produces a proportional displacement of the microrobot [START_REF] Rakotondrabe | Quadrilateral modelling and robust control of a nonlinear piezoelectric cantilever[END_REF] [START_REF] Rakotondrabe | Voltage/frequency proportional control of stickslip microsystems[END_REF]. The use of such a continuous actuation leads to a high cost, the need of a closed-loop control which takes into account the non-linearities [START_REF] Rakotondrabe | Voltage/frequency proportional control of stickslip microsystems[END_REF] and also the use of sensors which are difficult to integrate [START_REF] Haddab | Improvement of strain gauges micro-forces measurement using Kalman optimal filtering[END_REF]. The incremental actuation mode is brought into play in the stick-slip [START_REF] Rakotondrabe | Voltage/frequency proportional control of stickslip microsystems[END_REF] and the inchworm systems [START_REF] Chen | Locomotive gait generation for Inchworm-Like robots using finite state approach[END_REF]. The displacement of the microrobot is then discrete.

The need of sensors and closed-loop control is still necessary.

In this paper, we propose a new approach for the design and the actuation of microrobots based on the use of bistability. This new concept based on mechanical robustness allows open-loop control of the microrobots, low cost, and no external energy consumption to maintain a position.

In literature, several conventional robots based on the use of bistable mechanisms are described.

However, in the microworld, there is no microrobot based on the bistability concept although bistable structures are widely used. In the macroworld, the most typical discretely-actuated robot is the Variable Geometry Truss (VGT) manipulator (see Fig. 1) [10][12]. This planar binary manipulator consists of several modules. Each one includes three binary actuators. As a result, every module has eight states (2 3 ). Fig. 1. A 3-bit binary VGT, "0" or "1" represents one of the states. The δ is the displacement from state 0 to state 1. G.S. Chirikjian firstly presented the concept of the binary paradigm for robotic manipulators in [START_REF] Chirikjian | A binary paradigm for robotic manipulators[END_REF]. There have been several improvements in this binary hyper-redundant manipulator concept.

A lot of work about the calculation of the forward and inverse kinematics is done in [START_REF] Lees | A combinatorial approach to trajectory planning for binary manipulators[END_REF]- [START_REF] Uphoff | Inverse kinematics of discretely actuated hyper-redundant manipulators using workspace densities[END_REF].

Another robot called Robotic Articulated Intelligent Device (BRAID) is described in [START_REF] Sujan | Design of a lightweight hyper-redundant deployable binary manipulator[END_REF]; it is a 3-D configuration of binary actuated manipulators. There are 2 15 states by cascading 5 modules in the 3D workspace. Nevertheless, the VGT and the BRAID present some drawbacks: trajectory following is costly in term of computation because the number of logical state is huge [START_REF] Lees | A combinatorial approach to trajectory planning for binary manipulators[END_REF], [START_REF] Lees | An efficient method for computing the forward kinematics of binary manipulators[END_REF] and the inverse kinematics is hard to obtain as indicated in [START_REF] Lees | Inverse kinematics of binary manipulators with applications to service robotics[END_REF] [START_REF] Uphoff | Inverse kinematics of discretely actuated hyper-redundant manipulators using workspace densities[END_REF]. Moreover, the use of joints leads to a repeatability loss due to the joint clearances or gaps. So this kind of structure can not be used as an accurate positioning device. In the microworld, submicrometric resolutions are required. This is why the use of assembled manipulators and joints is avoided. Compliant joints are preferred in order to achieve high performances adapted to the microworld [START_REF] Gomm | Development of in-plane compliant bistable microrelays[END_REF], [START_REF] Howell | Compliant mechanisms[END_REF].

Miniaturization of the mechanical structure of VGT or BRAID is not suitable.

In this paper, we propose a new kind of microrobot to get over these problems. An approach to perform monolithic microfabrication of microrobots takes use of a modular concept and an open-loop control strategy. The microrobots, named "digital microrobots" are built from several elementary modules. Each one offers a very good repeatability and stable positions. An elementary module gives two stable states. The position of the whole microrobot is controlled by a digital word representing the state of the modules. This opens a new paradigm in the field of microrobotics, allowing the design of various kinematics adapted to the microworld.

The paper is organized as follows: first, the concept of microrobots based on bistable modules is presented. After that, the design of the bistable element is presented, and finally the conclusion and future work are given.

General concepts for digital microrobotics 2.1 Bistable module

The bistable module is the fundamental element for designing a digital microrobot. The bistable mechanism is a structure or mechanism that has two stable states, Fig. 2 shows the principle of bistable modules. Fig. 2. Principle of bistable modules, the δ or θ stands for motion between the two stable states.

We define state 0 or 1 respectively as one of the two stable states. Between these two stable states there are translations δ or rotation θ. The bistable mechanism is switched from one state to another by applying an external force, and no external energy is needed to maintain the stable states. Many applications of this structure can be found, for example, RF switches, relays, microvalves, etc.

Robot axis based on bistable modules

A robot axis can be built using several cascaded bistable modules. The number of logical states is thus increased. The final position of the axis is obtained by the accumulation of the discrete displacement (δ) or discrete rotation (θ). The design of each module and consequently their δ or θ defines the characteristics of the reachable workspace. Although only a discrete area can be reached, the number of logical states is an exponent function of the number of bistable modules. A robust mechanical design of the bistable module allows a high positioning resolution to be obtained. In this paper, we only study the case of linear displacements.

The number of the logical states of the robot axis is 2 n , where n is the number of bistable modules. Fig. 3 shows an example where 5 modules are used. There are 32 logical states. If the displacements of the modules are defined as: d, 2d, 4d, 8d and 16d (as shown in Fig. 3), where d is the minimum displacement needed, the range of the robot axis is:

( ) 5 1 1 1 1 . 1 2 4 8 16 .
. 2 31. 1 1

1 i i D d d d - =         = = =         ∑ (1)
The minimum displacement d represents the resolution of the axis position. Fig. 3. The reachable workspace of a robot axis with five cascaded bistable modules.

The position of the tip of the robot axis is represented by a binary word that could be used as a control instruction.

Discrete reachable workspace and Microrobot applications

The use of bistable modules to design microrobots leads to a discrete reachable workspace. If robust bistable modules with small displacements can be built, a high resolution is obtained and the position of the robot can be controlled in open-loop. Moreover, this approach offers good immunity to noise thanks to the mechanical performances.

However, a lot of modules must be built in order to obtain high resolutions. Microfabrication technology is well adapted to the fabrication of monolithic high performance digital microrobots which do not require any assembly process. From discretely actuated axis, robots with various kinematics can be built. Examples are given in Fig. 4.

Comparison between current microrobots and digital microrobots

In Table 1, we give a detailed comparison between the characteristics of current microrobots and the proposed digital microrobots.

Control

Closed loop [START_REF] Rakotondrabe | Quadrilateral modelling and robust control of a nonlinear piezoelectric cantilever[END_REF] Non-linear control [START_REF] Rakotondrabe | Quadrilateral modelling and robust control of a nonlinear piezoelectric cantilever[END_REF] Closed loop [START_REF] Rakotondrabe | Voltage/frequency proportional control of stickslip microsystems[END_REF] Control [START_REF] Rakotondrabe | Voltage/frequency proportional control of stickslip microsystems[END_REF] Open loop The current microrobots are mainly built by assembling the actuators and the mechanical structure.

(
Digital microrobots can be fabricated at low cost in a monolithic manner using microfabrication technology. Moreover, this novel concept has many advantages. Good repeatability and accuracy are obtained thanks to the mechanical performances of the bi-stable modules, that is why neither proprioceptive sensors nor bulky and expensive instruments for signal processing and control are needed. Low power consumption can be obtained because external energy is not needed to maintain the modules in a given state but only during the transition phases. The immunity to noise and to changes of the environment is improved. Using the parallel control of the modules, fast displacements can be obtained. Using bistable structures gives an approach that turns the difficulties of non-linear control into mechanical structure design.

The inchworm or stepper motors can offer also a discrete position. The principle of these motors depends on the friction between the actuator and the contact surface. Each step has an uncertain linked to friction surface. When it goes through the open loops of great number of steps, the final error could be great [START_REF] Rakotondrabe | Development, modelling and control of micro/nano positionning 2 dof stick-slip device[END_REF], [START_REF] Rakotondrabe | Voltage/frequency proportional control of stickslip microsystems[END_REF].

Digital microrobots present the advantages discussed before. However, they have some drawbacks: discrete reachable workspace and limited load. Let us consider an example of a piezoelectric unimorph beam, commonly used in microrobotics, with 14mm length that enables a maximum displacement of 60µm at its tip under 100V. Taking into account the creep effect, a resolution of 500nm can be reached. According to our proposal, to obtain similar performances, the smallest module must have a displacement of 500nm and 7 cascaded modules are necessary according to equation 1. Microfabrication of such an axis is difficult. Cascading several modules requires that each module carries the weight of the following ones which limits the load.

Design of bistable compliant structures

An elastic deformed structure usually goes back to its original shape after unloading the applied force, but some structures go to another position. A structure with the latter behavior is called a bistable structure. This type of structure shows similar characteristics between force and displacement or between elastic deformation energy and displacement. In 2002, M.S. Baker [START_REF] Sujan | Design of a lightweight hyper-redundant deployable binary manipulator[END_REF] presented on-chip actuation of a compliant bistable mechanism. In which the V-shape thermal actuators were used. In 2004, Jinqiu [START_REF] Qiu | A curved-beam bistable mechanism[END_REF] presented a curved-beam bistable mechanism for designing a relay, in which the actuator and the bistable structure were made in the wafer by micro fabrication. Many bistable structures were fabricated through MEMS technologies [START_REF] Qiu | A curved-beam bistable mechanism[END_REF][23] [START_REF] Dong | Mechanical design and modeling of MEMS thermal actuators for RF applications[END_REF].

One of them has been widely studied: the curved-beam bistable mechanism. This compliant mechanism does not need any assembly. Moreover, it can be fabricated by batch processing and it perfectly suits our objective.

To define the dimensions of the bistable mechanism, two factors are considered: the fabrication limits that define the minimal width of the beams of the structure and the size of the thermal actuators required to switch the bistable mechanism. According to the pseudo-rigid-body model [START_REF] Howell | Compliant mechanisms[END_REF], [START_REF] Hwang | Modeling and experimental characterization of the chevron-type bistable microactuator[END_REF], a calculated model is presented in Fig. 5. The bistability is determined by the torsional spring stiffness ks, the elastic beam stiffness kb, and the initial height of the curved beam h.

× - + - + = Y h L h L x (2) ) ( tan ) ( tan 1 1 L Y h L h - - = - - θ (3)
L, h and Y are shown in Fig. 5.

The strain energy U can be derived from the torsional spring stiffness ks and the elastic beam stiffness kb:

2 2 2 1 4 2 1 8 x kb ks U ⋅ ⋅ ⋅ + ⋅ ⋅ ⋅ = θ (4) Y x x kb Y ks Y U F ∂ ∂ ⋅ × + ∂ ∂ ⋅ ⋅ = ∂ ∂ = 4 8 θ θ (5) 2 2 2 2 2 2 2 1 1 ) ( ) ( ) ) ( ( 4 ) ) ( 1 ( 1 )) ( tan ) ( (tan 8 Y h L Y h Y h L h L kb L L Y h L Y h L h ks Y U F - - - ⋅ - + - + ⋅ ⋅ + - + ⋅ - - ⋅ ⋅ = ∂ ∂ = - - (6)
Ignoring the high order about Y, The force F can be written as:

L Y h L Y h L Y hY kb L ksY F 4 ) ( 3
) ( 2 ) 2 ( 4 8 2 2 2 - + - ⋅ - + ≈ (7)
The two stable states Y s1 and Y s2 and the unstable state Y u are obtained for F=0.

L kb L ks h L kb h L kb h ks Y L kb ks ks L kb L ks h L kb h L kb h ks Y Y u s s ⋅ ⋅ - ⋅ ⋅ ⋅ ⋅ - ⋅ ⋅ + ⋅ ⋅ - ⋅ ⋅ ⋅ - ⋅ ⋅ = ⋅ ⋅ - ⋅ ⋅ ⋅ ⋅ - ⋅ ⋅ + ⋅ ⋅ + ⋅ ⋅ ⋅ - ⋅ ⋅ = = (8)
In our case, ks and kb can be calculated according to the pseudo-rigid-body theory. Where E is the Young's modulus, I is the inertial moment and A is the section's area.

According to this model, the required actuation force and displacement stroke can be found. The dimensions adapted to the microfabrication limits in our clean room facility are shown in Fig. 6. These dimensions are defined by making the trade-off between the fabrication limits and the actuation force. In our case, a force of 3mN has been considered according to the actuators defined in the next section. This choice allows the fabricated bistable mechanism to be actuated using the designed force. Fig. 6. Curved beams dimensions TABLE 2 The process summary to fabricate the bistable module. This bistable mechanism is fabricated using silicon on insulator (SOI) wafers with 100 m thick device layer, 300 m thick handle and 1 m sacrificial oxide layer. The handle layer is used to support the entire structure, and the bistable structure is fabricated in the device layer.

Firstly, the photoresist is deposited, and then the device layer is etched using DRIE (Deep Reactive Ionic Etching). After that, the aluminum mask is made for the etching of the handle layer. Finally, the oxide layer is etched using hydrofluoric (HF) acid. Since the last step is done in the liquid, some of structures are broken. Table 2 summarizes the fabrication process. The fabricated structure is shown in Fig. 7. The etching of the handle layer permits to free the structure by an aluminum mask. A finite element analysis (FEA) is made by using the commercial FEA software ANSYS. The results from the model, the FEA and the static experimentation are presented in Fig. 9.

100µm Fig. 9. Force (F) vs displacement (Y) characteristic. This static experimentation is made by using a commercial force sensing probe (ST-S270 from FEMTO TOOLS [START_REF] Tools | FT-S Force Sensing Probes[END_REF] with a sensitivity of 899.2µN/V).

This model gives a greater stiffness than the stiffness from FEA because of the simplification of the torsional spring model but it can be considered as a sufficient approximation. Therefore, for switching from a stable state to another one, a microactuator able to produce a force of 2.15mN with a displacement of 65µm is needed. The next section is dedicated to the design of the actuator.

Actuator design

To switch the bistable mechanism, the potential candidates could be piezoelectric, thermal, electrostatic and shape memory alloy (SMA) actuators. Due to the use of the bulk micromachine technology, piezoelectric actuators and SMA actuators are difficult to integrate. The calculation of the electrostatic actuator shows that it takes a lot of space to reach the required force. So we decided to use thermal actuators.

A U-shape thermal actuator has been built from two beams. One beam is called "hot beam" and the other is called "cold beam" [START_REF] Dong | Mechanical design and modeling of MEMS thermal actuators for RF applications[END_REF]. This actuator takes use of the asymmetric deformation to create in-plane displacement. For a given temperature, the tip's displacement magnitude mainly depends on the length of the thermal actuator [START_REF] Dong | Mechanical design and modeling of MEMS thermal actuators for RF applications[END_REF], [START_REF] Cragun | Thermal microactuators for microelectromechancial Syst. (Mems)[END_REF]. Fig. 10 shows the variation of the length Vs tip's displacement at 300 o C. We use commercial software ANSYS which offers a structural thermoelectric tetrahedron element (SOLID227) permitting to couple electro-thermal-mechanical analysis. We took into account all the thermal effects: the conduction, the convection, and the radiation (see Fig. 10). The thermal actuator is fabricated in the device layer of the SOI wafer, in the same way as the bistable mechanism. To obtain a control voltage lower than 24V which is supplied by common voltage generators, we have chosen a wafer with electric resistivity of 0.02 ohm.cm. The pads of the actuator's supply are electrically insulated from the other parts of the structure in order to avoid short circuits but the handle part is kept for the integration of the whole structure. The contact pads are made of aluminum to ensure a good ohmic contact through the diffusion of Al atoms in the silicon. Fig. 12 shows two microfabricated thermal actuators. According to [START_REF] Mankame | The effect of thermal boundary conditions and scaling on Electro-thermalcompliant micro devices[END_REF][30], the thermal convection and radiation can be neglected in the microworld, so the heat is accumulated in the structure. Energy of 12mJ is needed for every switching operation. This energy is quickly transferred to the substrate thanks to the good thermal conductivity of the silicon. Two pairs of thermal actuators are used to switch the bistable structure. This energy is needed only for switching. Since the positions are given by the stop blocks, no energy is required for maintaining these positions.

Stop block design

In order to ensure the bistability and a blocking force to the bistable structure in each state, we need to limit the displacement of the structure before it reaches its maximum position. The use of two stop blocks is thus necessary. The previous realizations [START_REF] Qiu | A curved-beam bistable mechanism[END_REF], [START_REF] Hwang | Modeling and experimental characterization of the chevron-type bistable microactuator[END_REF] of bistable mechanism can not meet this requirement. This bistable mechanism has to move between the blocked positions in order to produce a blocking force for each of the two positions. It is difficult to obtain naturally two pre-stressed positions with a high blocking force because of the monolithic microfabrication.

To solve this problem, we designed dedicated compliant mechanisms. The first stop block is a compliant mechanism (see Fig. 15). The beam deforms under an external activation force, and finally is forced into the stop block. The second one is located at the other side of the unstable position (see Fig. 16). The distance between these two stable positions can be chosen according to the design and the targeted displacement of the module. The position of the two stop blocks defines the displacement of the module and the blocking forces (the maximum forces that can be applied to the shuttle without losing the contact between the shuttle and the stop blocks).

Due to the monolithic microfabrication of the overall module, a first activation is necessary in order to insert the shuttle between the two stop blocks. The first activation gives a pre-energy in the structure which allows reducing the actuation energy to switch the bistable structure. The whole bistable module which contains two bistable beams, four thermal actuators, a frame and a shuttle is finished (see Fig. 17). Compared to the previous works on bistable structures, our design allows getting a blocked force in two positions. Fig. 18 The experimental set-up. In this test, the bistable modules are prepared for an in-plane measurement. For this measurement, we use a high resolution (10 nm) interferometer from SIOS Technology to measure the displacement of the bistable mechanism. All the test devices are mounted on an antivibration table (see Fig . 18). In the module shown in Fig. 18, the measured displacement of the shuttle is 10.1µm. The stop blocks give 1.54mN blockage force to limit the displacement of the shuttle.

Overall designs

Because the overall structure is composed of a large number of similar bistable modules, it is very important to consider the relationship between the overall structure and the bistable module.

Firstly, two cascaded bistable modules are considered (see Fig. 20). Since the two curved-beams will support all the structure weight in the first bistable module, it is necessary to study the gravity effect in the structure. Fig. 20. n cascaded bistable modules.

In order to find the out-of plan stiffness caused by gravity, we performed a calculation model presented in Fig. 21. In this calculation model, the two bistable curved beams are considered as two springs, and the actuators and the frame are considered as a rigid beam. Only the spring causes a flexion. The gravity center of every bistable module is located in the middle of the two springs.

For n bistable modules, we simplify the model as in Fig. 22. The total moment M T = (M (n 2n)/2+MG) and total force F T = (nF+G) are the equivalent moment and force of the n bistable modules which are applied on the first bistable module. The out of plane flexion can be obtained:

k L L L L L n L G L L n n F k G nF L L k L M k F T T n 2 1 1 1 2 1 2 5 . 0 ) ( ) 2
)) )( 1 The total displacement can be obtained:

∑ ∑ = = + + - + + + - + + = + + = n n n n T T T k L L L L L n L G L L n n F k G nF L L k L M k F 1 2 1 1 1 2 1 1 2 5 . 0 ) ( ) 2 )) )( 1 ( 5 . 0 ( 2 ) )( ( ( ) ( 5 . 0 δ (12)
The spring is a variable section cantilever beam, of which stiffness k can be determined by a classic cantilever model:

k = × + + × + + × × = ( 13
)
Where E is the Young's modulus, I is the inertial moment and l is the length of cantilever beam.

Each bistable module has a mass of 2.128mg. We simulate the deflection of cascaded structures carrying at the tip a silicon sphere with a diameter of 300µm (mass: 0.0325mg). The comparison of the tip's flexion between the model and the FEA is presented in table 3 for cascaded bistable modules with 1 to 4 modules. The tip's flexion is accumulated with the increase of the bistable modules. The analysis shows that there isn't much stress in the first bistable module. Nevertheless, the tip of the structure has an important flexion. According to the equations 12 and 13, we can obtain:

      ∝       ∝       ∝ 3 1 1 1 h I k n δ (14)
According to the equation 14, the tip's flexion could be reduced with the increase of the thickness of the structure. In our case, the tip's flexion is 2µm with a thickness of 100µm.This flexion is neglected.

Conclusion and future work

In this paper, we have presented a new concept for microrobotics: the digital microrobot. More precisely, we have designed the elementary component of this kind of microrobot: the microfabricated bistable module. This module is fabricated in a monolithic manner on a SOI wafer.

We have described and characterised each part of this module: the bistable compliant structures, the thermal actuators and the stop blocks which enable to give robust positions and blockage forces. The batch microfabrication permits to fabricate monolithic structures containing several modules without any assembly. For example, a microrobot axis is built by cascading several bistable modules to reach a discrete workspace. This digital microrobot takes advantage of bistable mechanism, such as repeatability and reliability. The proposed design of bistable module allows getting two stables positions with a blocking force useful for micromanipulation operations.

The study shown in this paper opens the possibility to cascade several modules to get robust positions in a discrete workspace. The main advantages of the digital microbotics approach are: repeatability and robustness obtained with mechanically blocked positions, and ease of control because of the open loop condition. Moreover, because each module can be triggered independently, parallel control allows reducing the time to get a position. Future work will include the dynamic characterization of the dynamic of bistable mechanism, and the fabrication of several cascaded bistable modules. Fig. 1. A 3-bit binary VGT, "0" or "1" represents one of the states. The δ is the displacement from state 0 to state 1. Fig. 2. Principle of bistable modules, the δ or θ stands for motion between the two stable states.

Fig. 4 .

 4 Fig. 4. (a) Parallel structured microrobot. (b)Serial structured microrobot

Fig. 5 .

 5 Fig. 5. An equivalent model for bistable beam The angle θ of torsional spring and the deformation x of the elastic beam can be expressed as a function of the vertical displacement Y, which is caused by an external force F.8 . 0 ) ) ( ( 2 2 2 2

Fig. 7 .

 7 Fig. 7. Microfabricated module We use the experimental set-up to measure the characteristics of the bistable module. The fabrication result is shown in Fig. 8.

Fig. 8 .

 8 Fig. 8. The experimentation set-up.

Fig. 10 .

 10 Fig. 10. Tip's displacement Vs length of actuator without any load.

Fig. 11 .

 11 Fig. 11. Designed thermal actuator.

Fig. 12 .

 12 Fig. 12. Two microfabricated thermal actuators. The static displacement of the tip of the thermal actuator is shown in Fig. 13.

13

 13

  .

Fig. 13 .

 13 Fig. 13. The static behavior of the thermal actuator.

Fig. 14

 14 Fig.14shows the response of the actuator to a square input voltage with 8 volts amplitude. The settling time is about 0.37s for a rising edge and 1.38s for a falling edge. The difference is due to the dissipation of the heat in the silicon. The characteristics required to switch the bistable structure are obtained with this actuator. Indeed, at a displacement of 140µm without any load (which is equivalent to a force of 2.2 mN at a displacement of 65µm) is obtained under a driving voltage of 20V.

Fig. 14 .

 14 Fig. 14. Dynamic response of the thermal actuator.

Fig. 15 .

 15 Fig. 15. (a) The stop block after microfabrication, (b) the stop block after insertion of the shuttle

Fig. 16 .Fig. 17 .

 1617 Fig. 16. The locations of the stop blocks in the force/displacement and energy/displacement characteristics.

Fig. 19

 19 Fig. 19 shows photos of the stop blocks during operation.

Fig. 19 .

 19 Fig. 19. (a1) The position just after fabrication (a2) The activated position (b1) The activated position (b2) The switched position.

Fig. 21 .

 21 Fig. 21. The equivalent out-of plane model of bistable module.

Fig. 22 .

 22 Fig. 22. Model of n cascaded bistable modules.

Fig. 3 .

 3 Fig.3. The reachable workspace of a robot axis with five cascaded bistable modules.

Fig. 4 .Fig. 8 .

 48 Fig. 4. (a) Parallel structured microrobot. (b)Serial structured microrobot Fig. 5. An equivalent model for bistable beam Fig. 6. Curved beams dimensions Fig.7. Microfabricated module Fig. 8. The experimentation set-up.

Fig. 9 .Fig. 10 .

 910 Fig.[START_REF] Chen | Locomotive gait generation for Inchworm-Like robots using finite state approach[END_REF]. Force (F) vs displacement (Y) characteristic. This static experimentation is made by using a commercial force sensing probe (ST-S270 from FEMTO TOOLS[START_REF] Gomm | Development of in-plane compliant bistable microrelays[END_REF] with a sensitivity of 899.2µN/V) Fig.10. Tip's displacement Vs length of actuator.

Fig. 11 .

 11 Fig. 11. Designed thermal actuator.

Fig. 12 .

 12 Fig. 12. Two microfabricated thermal actuators.

Fig. 13 .

 13 Fig. 13. The static behavior of the thermal actuator.

Fig. 14 .

 14 Fig.[START_REF] Lees | A combinatorial approach to trajectory planning for binary manipulators[END_REF]. Dynamic response of the thermal actuator. The settling time is about 0.37s for a rising edge and 1.38s for a falling edge. The difference is due to the dissipation of the heat in the silicon.

Fig. 15 .Fig. 16 .

 1516 Fig. 15. (a) The stop block after microfabrication, (b) the stop block after insertion of the shuttle Fig. 16. The locations of the stop blocks in the force/displacement and energy/displacement characteristics.

Fig. 17 .

 17 Fig. 17. One bistable module: four thermal actuators, two bistable beams, stop blocks, and a frame.

Fig. 18 .

 18 Fig. 18. The experimental set-up.

Fig. 19 .

 19 Fig.19. (a1) The position just after fabrication (a2) The activated position (b1)The activated position (b2) The switched position.

Fig. 20 .

 20 Fig. 20. n cascaded bistable modules.

Fig. 21 .

 21 Fig. 21. The equivalent out-of plane model of bistable module.

Fig. 22 .

 22 Fig.22. Model of n cascaded bistable modules.

TABLE 1 A

 1 Comparison between current microrobots and digital microrobots

	Characteristic	CURRENT MICROROBOTS	DIGITAL MICROROBOTS
	Actuation approach	ACTUATION PROPORTIONAL	INCREMENTAL ACTUATION	DISCRETE ACTUATION
	Size	Medium	Medium	Small (microfabrication)
	Cost	High	Medium	Low

  no sensors needed)

	Sensitivity to noise	High [8]	Low	Low
	Use of sensors	Yes	Yes	No
	Displacement	Continuous	Discrete	Discrete
	Range	Small	Theoretically unlimited [6]	Depends on the number of bistable module

TABLE 3 tip

 3 's flexion comparison between the model and the FEA

	n	Model (µm)	FEA (µm)
	1	0.24	0.22
	2	2.11	1.9.1
	3	16.8	14.4
	4	41.6	38.2

TABLE 1 A

 1 Comparison between current microrobots and digital microrobots TABLE 2 The process summary to fabricate the bistable module. TABLE 3 tip's flexion comparison between the model and the FEA

Acknowledgment

The fabrication of the bistable modules has been realized in the clean room of Femto-st institute.

The authors would like to thank the engineers of the clean room and Joel Agnus for sharing the microfabrication know-how.