Cindie Andrieu
email: cindie.andrieu@ifsttar.fr

Guillaume Saint
email: guillaume.saintpierre@ifsttar.fr

X Bressaud
email: bressaud@math.univ-toulouse.fr

A functional analysis of speed profiles: smoothing using derivative information, curve registration, and functional boxplot

Keywords: Functional Data Analysis, Smoothing spline, Landmarks alignment, Functional boxplots, Speed profiles

come

A functional analysis of speed profiles: smoothing using derivative information, curve registration, and functional boxplot

Introduction

The knowledge of the actual vehicle speeds on roads is essential from several points of view: to locate blackspot in the network, to improve the knowledge of travel time and to evaluate the effects of the modification of the infrastructure (addition of speed bumps, roundabouts, ...). The speed choice of drivers is one of the most important components of their behavior and also their road usage. This continuous information of road user's speed is available with the development of probe vehicles, that can be seen as mobile sensors exploring continuously the road network. More particularly, the development of smartphones equipped with a GPS (Global Positioning System) has increased the number of digital "traces" left by vehicles, and leads to the obtention of individual space-speed profiles that represent speed as a function of vehicle position.

The collection of individual space-speed profiles can leads to large volume of data that require the use of appropriate methods. Indeed, since in practice space-speed profiles are composed of time-stamped measurements of speed and position, most studies consider them as R n vectors where n is the number of measurements. However, advances in sensors technology enable to collect data with high sampling rate that leads to high dimensional vectors (n is very large), for which classical multivariate statistical methods become inadequate because of problems related to the so-called "curse of dimensionality" and the significant correlation between close observations. The originality of the approach presented in the current paper is to propose a functional analysis of space-speed profiles, i.e. to treat these objects as functions rather than vectors. This approach takes inspiration from Functional Data Analysis, a statistical domain that has developed considerably over the last twenty years and that appears in several domains such as meteorology, chemometrics or economics (e.g. [START_REF] Ramsay | Applied Functional Data Analysis: Methods and Case Studies[END_REF][START_REF] Febrero | A functional analysis of nox levels: location and scale estimation and outlier detection[END_REF], but that is not yet widely used in road transport. Yet the functional approach is particularly suitable for the analysis of speed profiles since it allows to preserve the physical consistency between speed and position (and implicitly time), and their functional characteristics: computation of derivatives (that leads to acceleration or jerk profiles), regularity, shape constraints... An overview of the theory of statistics with functional data can be found in the monographs by Ramsay andSilverman (2002, 2005) or [START_REF] Ferraty | Nonparametric Functional Data Analysis: Theory and Practice[END_REF], and in the reviews by [START_REF] Levitin | Introduction to functional data analysis[END_REF], [START_REF] Valderrama | An overview to modelling functional data[END_REF] or more recently [START_REF] Cuevas | A partial overview of the theory of statistics with functional data. J Stat Plan Inference Cuevas A[END_REF]. Thus, after a description of a real data set composed of 78 individual space-speed profiles, we begin our analysis by a functional modeling of space-speed profiles with a definition of the corresponding functional space and the study of some mathematical properties (continuity and differentiability).

Then, the first step in the functional analysis of space-speed profiles is to convert the raw data including speed and position measurements into functional objects. It is shown that this smoothing problem can be viewed as an interesting nonparametric regression problem that needs to take into account two constraints: the use of the derivative information, and a monotonicity constraint. In this paper, we propose a two-step estimator : smoothing step, and then monotonization step. We show that the smoothing step with the constraint of the use of the derivative information can be seen as a special case of the general spline smoothing problem (see [START_REF] Wahba | Spline models for observational data[END_REF][START_REF] Wang | Smoothing Splines: Methods and Applications[END_REF] and can be solved by using the theory of reproducing kernel Hilbert spaces. Thus, the estimator proposed in this paper can be written as a linear combination of basis functions and kernel functions. However, from a computational point of view, we show that the use of a semi-kernel in place of reproducing kernel is more appropriate, and we propose to use the theory of thin-plate spline [START_REF] Wahba | Spline models for observational data[END_REF][START_REF] Wahba | Some new mathematical methods for variational objective analysis using splines and cross-validation[END_REF] in order to obtain an estimator with a simpler form. Then, a monotonization step is proposed based on a method developed by [START_REF] Ramsay | Estimating smooth monotone functions[END_REF] which has the advantage of being relatively simple to implement. However, if this smoothing procedure leads to a set of individual spacespeed profiles, when the volume of data is large, it is necessary to summarize the information contained in this set. So, in a second time, we proposed a methodology of construction of an aggregated speed profile, such as the average profile. It is then necessary to use curve registration method in order to correct phase variation (especially at stops' location), and to obtain a representative speed profile with similar features of corresponding individual speed profiles. Finally, in a third time, we propose to apply the functional boxplot developed by Sun and Genton (2011) using an appropriate functional depth to the set of individual space-speed profiles. This graphical tool which is an extension of the classical boxplot used in the univariate setting, is very interesting to explore the variability of a functional data set. The application of this tool to speed profiles data set leads to the construction of speed corridors that reflect the variability between road users and are very informative about actual operating speed. This speed corridors are particularly adapted to driving assistance system and to enrich or update digital maps.

The remainder of the paper is structured as follows. In Sect. 2, the real data set of individual space-speed profiles used in this study is presented. In Sect. 3, we propose a definition of the functional space of space-speed profiles and we study some mathematical properties (continuity and differentiability) of these functions. In Sect. 4, we propose a two-step smoothing procedure : smooth using derivative information, and then monotonize. Performance and limitations of the estimator are discussed on simulation studies and the real data set. In Sect. 5, we propose a methodology of construction of an average speed profile using a curve registration procedure based on the method of landmarks alignment. In Sect. 6, we apply the functional boxplot to the real data set of individual speed-profiles and show the interest of these speed corridors to explore the variability between road users. Finally, Sect. 7 provides the main conclusions of the present study.

The data

We consider in this paper a data set extracted from an experiment conducted by the French laboratory IFSTTAR-LIVIC and that took place in 2012 in Versailles, France. Thirty-nine drivers participated to this experiment and performed twice a road section of urban and inter-urban type with a length of about 1100 m. This road section, illustrated at Fig. 1, corresponds to the path from A to B and is composed of a stop sign, two roundabouts and a traffic light. For logistical reasons, two vehicles were used for this experiment: a Renault Clio III equipped with a Garmin GPS 16x LVC (for 20 drivers), and a Renault Modus with a GPS GlobalSat BR-355 (for 19 drivers). Note that the use of two vehicles and the fact that each driver performed twice the studied section lead to conditions close to naturalistic driving studies where different drivers were observed in a natural setting, in particular during regular travels such as the commute to work. Thus, we do not take into account the correlation between the two paths of the same driver. The data-logger collects vehicle position (latitude and longitude) and speed via GPS at a 1Hz sampling frequency (i.e. 1mes/sec). Note that GPS receivers use the Doppler shifts of the satellite signals to calculate vehicle speed, which implies that vehicle speed is independent of vehicle position. In order to reduce the GPS position measurements to a one-dimensional framework, it is assumed that these are map-matched, so that the vehicle is positioned on the correct road segment. Many map-matching algorithms have been developed to identify the correct road segment on which the vehicle is travelling. Thus, in this paper, GPS measurements represent the curvilinear abscissa of the vehicle on the studied road segment (absolute location) from the initial position (point A).

The aim of the study is to focus on space-speed profiles from this data set, i.e. speed as a function of the distance traveled by the vehicle from the point A. Since the 39 drivers performed twice the studied section, the data set is composed of 78 individual space-speed profiles illustrated at Fig. 2. 3 Functional modeling of space-speed profiles 3.1 Definition of space-speed profiles Before beginning a functional analysis of space-speed profiles, it is necessary to define the functional space of such objects. Indeed, any function f : R + -→ R + is not a space-speed profile (e.g. a constant function equal to zero). In practice, a space-speed profile is a sequence of time-stamped measurements of position (from GPS or odometer) and speed, so it can be studied in the three following study areas: distance × time, speed × time and speed × distance (see Fig. 3). The functions defined in each of these three study areas are related mathematically: if we denote F (t) a function defined in the study area distance × time that represents the distance traveled as function of time, the derivative function F ′ (t) represents the speed as function of time and is defined in the study area speed × time. So, by definition, the function F must be increasing and at least of class C 1 . In order to define an acceleration profile, we propose to require that F is at least C 2 , and so we propose the following definition of the functional space of space-speed profiles :

Definition 1 Let x f ∈ R + .
Then the space of space-speed profiles, denoted E SSP , is defined as follows : E SSP = {v S : [0, x f ] -→ R + such that there exists a positive real T and an increasing function

F : [0, T ] -→ [0, x f ] of class C 2 with F (0) = 0 such that v S (x) = F ′ • F -1 (x)), x ∈ [0, x f ]}, where F -1 is the generalized inverse of F defined by F -1 (x) = inf {t ∈ [0, T ], F (t) = x}.
The positive real numbers x f and T represent respectively the length and the travel time of the studied section. Fig. 4 illustrates Definition 1 by showing the functional link between distance, speed and implicitly time.

Mathematical properties of space-speed profiles

We studied some properties of the space-speed profiles, i.e. functions in the space E SSP (as defined in Definition 1). The continuity property is given by the following theorem whose proof is deferred to Appendix A:

Theorem 1 All functions v S : [0, x f ] -→ R + belonging to the space of spacespeed profiles E SSP (as defined in Definition 1) are continuous on [0,

x f ].
If the continuity property of space-speed profiles is obvious, the differentiability property is less intuitive as shown in the following theorem whose proof is also deferred to Appendix A: Theorem 2 Assume that v S : [0, x f ] -→ R + belongs to the space of spacespeed profiles E SSP (as defined in Definition 1). Let H 0 = {x ∈ [0, x f ], v S (x) = 0}, all points for which the speed is zero. The two following assumptions are added:

(H 1 ) Assume that F is of class C 2 on [0, T ] and strictly increasing, and ∃t 0 ∈ ]0, T [ such that F ′ (t 0 ) = 0 and F ′′′ (t 0 ) exists with F ′′′ (t 0 ) = 0. (H 2 ) Assume that F is of class C 2 on [0, T ] and increasing, ∃t 0 , t 1 ∈]0, T [, t 0 = t 1 such that F ′ (t) = 0 on [t 0 , t 1 ], and the function G defined on [0, T -(t 1 -t 0 )] by: for t ≤ t 0 , G(t) = F (t), for t ≥ t 0 , G(t) = F (t + t 1 -t 0 ), satisfies the assumptions (H 1 ).

If F satisfies the assumptions (H 1 ) or (H 2 ), then v S = F ′ • F -1 is not differentiable on H 0 .

The assumptions (H 1 ) and (H 2 ) are not restrictive and are satisfied in most cases. So, this theorem shows that space-speed profiles are not differentiable at points for which the speed is zero, i.e. when the vehicle is stopped. From a geometrical point of view, if we assume that x 0 ∈ H 0 (i.e. v S (x 0 ) = 0), it is easily shown that the graph of a space-speed profile v S has a half-tangent parallel to the axis of y in x 0 , i.e. a cusp at the point (x 0 , 0). This property of non differentiability at points for which the speed is zero, implies some difficulties in the calculation of an average profile, particularly in the case of stops. Indeed, if a space-speed profile v 1 is equal to zero at a point x 0 , and a space-speed profile v 2 is strictly positive at x 0 , then the sum v 1 + v 2 is not a space-speed profile as defined in the Definition 1 since v 1 + v 2 is not differentiable at x 0 but (v 1 + v 2 )(x 0 ) > 0. Thus, the calculation of an average profile is meaningful only in this two cases:

when all space-speed profiles are strictly positive (no stops) ; -when all space-speed profiles are equal to zero at the same points (i.e. all vehicles stop at the same location).

This second case raises the issue of registration of speed profiles that will be discussed in Sect. 5.

4 Estimation of a space-speed profile from noisy data: A smoothing problem under constraints

The first step of a functional analysis is to convert the raw data into functional objects that leads to the use of an adapted smoothing procedure. However, the estimation of a space-speed profile from noisy measurements of position and speed is a complex nonparametric regression problem (see for example [START_REF] Andrieu | Estimation of space-speed profiles: A functional approach using smoothing splines[END_REF]. Indeed, on the one hand, both the response variable (corresponding to speed) and the explanatory variable (corresponding to vehicle position) are noisy. And on the other hand, the regression function must belong to the space E SSP defined in Definition 1 and then check its properties, in particular the non differentiability when the speed is zero. To overcome these difficulties, we propose to change to a more suitable study area and start by estimating the function F representing the distance traveled as function of time (study area distance × time in Fig. 3). Then the new nonparametric model is

y i = F (t i ) + ε x,i , i = 1, . . . , n, (1)
where y i are noisy observations of the distance traveled, ε x,i are uncorrelated errors with zero mean and σ 2 x variance, and F (t) is the regression function. This change of study area leads to take into account the two following constraints:

1. Use the derivative information, i.e. estimate the regression function F (t) of the model (1) from both noisy observations of F corresponding to position measurements, and noisy observations of its derivative F ′ corresponding to speed measurements. 2. A monotonicity constraint since the function F representing the distance traveled as function of time must be increasing.

The consideration of these two constraints in the smoothing step is the subject of the two following subsections. Then, once we have obtained an estimator F of the function F , it is easy to deduce by differentiation an estimator F ′ of the time-speed profile, and finally to deduce an estimator v S of the space-speed profile v S by the transformation v S = F ′ • F -1 .

Smoothing using derivative information

The first constraint in the nonparametric model ( 1) is to use the derivative information, i.e. to estimate the regression function F (t) from both noisy measurements of position and speed. Assume that the domain of

F (t) is X = [0, T ],
where T is a positive real, and F ∈ W m [0, T ] where W m [0, T ] is the Sobolev space of order m with m > 1. Then the use of the derivative information leads to consider the following data model:

y i = F (t i ) + ε x,i , i = 1, . . . , n v i = F ′ (t i ) + ε v,i , i = 1, . . . , n (2)
where y i and v i are noisy measurements of distance traveled and speed respectively at each sampling time t i , and ε x,i and ε v,i are independent zero mean errors with variance σ 2 x and σ 2 v respectively. We also assume that for all i = 1, . . . , n, ε x,i and ε v,i are independent. Note that we have assumed that the observations y i and v i are obtained at the same times t i . Otherwise, a data resampling will lead to this case.

The problem of smoothing with derivative information appears in various applications such as economy [START_REF] Hall | Nonparametric estimation when data on derivatives are available[END_REF], molecular biology [START_REF] Calderon | P-splines using derivative information[END_REF] or image analysis [START_REF] Mardia | Kriging and splines with derivative information[END_REF]. We propose to solve this problem by using smoothing splines, which have the advantage of requiring the estimation of a single smoothing parameter λ, contrary to penalized splines used by [START_REF] Calderon | P-splines using derivative information[END_REF], for which the choice of the number of knots can be difficult. Thus, this problem can be seen as a special case of the general spline smoothing problem (see [START_REF] Wahba | Spline models for observational data[END_REF][START_REF] Wang | Smoothing Splines: Methods and Applications[END_REF] and can be solved by using the theory of reproducing kernel Hilbert spaces (see [START_REF] Cox | Approximation of method of regularization estimators[END_REF]. Then, the estimator can be written as a linear combination of basis functions and kernel functions. However, to compute the estimator, it is necessary to choose a norm associated with the design space (in this study, the Sobolev space W m [0, T ]) that is suitable. Indeed, the expression of the reproducing kernels used for calculating the estimator depends on the choice of this norm, and can leads to difficulties in the numerical computation (see Berlinet and Thomas-Agnan 2004 for a collection of examples of spaces, norms, and kernels, and Andrieu 2013, Chap. 4 for more details). So, we propose to use the theory of thin-plate spline which leads to similar results for the form of the estimator (see [START_REF] Duchon | Splines minimizing rotation-invariant semi-norms in sobolev spaces[END_REF][START_REF] Meinguet | Multivariate interpolation at arbitrary points made simple[END_REF][START_REF] Wahba | Some new mathematical methods for variational objective analysis using splines and cross-validation[END_REF], but where the reproducing kernel is replaced by a semi-kernel with a simpler form that greatly simplifies the computation of the estimator.

The model ( 2) is a particular case of the general spline smoothing model defined in [START_REF] Wahba | Spline models for observational data[END_REF], since it can be rewritten as follows:

y j = L j F + ε j , j = 1, . . . , 2n, (3)
where -the observations y j are defined by: y j = y i with i = j for j = 1, . . . , n y j = v i with i = j -n for j = n + 1, . . . , 2n ,

the bounded linear functionals L j on W m [0, T ] are defined by:

L j F = F (t i ) with i = j for j = 1, . . . , n L j F = F ′ (t i ) with i = j -n for j = n + 1, . . . , 2n ,
the errors ε j are defined by:

ε j = ε x,i avec i = j pour j = 1, . . . , n ε j = ε v,i avec i = j -n pour j = n + 1, . . . , 2n .
A useful result state that the operator L defined by Lf = Berlinet and Thomas-Agnan 2004, Th. 133;[START_REF] Wahba | Some new mathematical methods for variational objective analysis using splines and cross-validation[END_REF]. So, we can deduce that linear functionals in the model (3) are bounded on W m [0, T ] if m > 1 since d = 1, k = 0 for j = 1, . . . , n and k = 1 for j = n + 1, . . . , 2n. Then, an estimator of F (t) is the minimizer of the following penalized least squares criterion in W m [0, T ]:

∂ k f ∂ α 1 x1...∂ α d x d for α 1 + . . . + α d = k (k, α 1 , . . . , α d ∈ N) is a continuous linear form if and only if 2m -2k -d > 0 (see
1 2n {σ -2 x n i=1 (y i -F (t i )) 2 + σ -2 v n i=1 (v i -F ′ (t i )) 2 } + λ T 0 (F (m) (t)) 2 dt. (4)
A solution to a more general minimization problem extended to dimension d ≥ 1 of which the minimization problem ( 4) is a special case with d = 1, is given in [START_REF] Wahba | Some new mathematical methods for variational objective analysis using splines and cross-validation[END_REF]. Before stating the result, we introduce the m polynomials functions φ 1 , . . . , φ m as a basis of the null space of the penalty functional

J m (f ) = T 0 (f (m) (t)) 2 dt, i.e. H 0 = {f : J m (f ) = 0}.
Thus, provided that the two hypothesis on the linear functionals (L 1 , . . . , L n linearly independent continuous linear functionals, and L k p ν=1 a ν φ ν = 0 implies that all the a ν are 0) are satisfied, which it's the case if the sampling points t 1 , . . . , t n are distincts, we can deduce that the minimizer of (4) can be written as:

F λ (t) = m ν=1 d ν φ ν (t) + n i=1 c i E m (t i , t) + n i=1 c ′ i ∂ ∂s E m (s, t)| s=ti , (5)
where

E m (s, t) = θ m,1 |s -t| 2m-1 , θ m,1 = Γ (1/2 -m) 2 2m π 1/2 (m -1)! . The coefficients c = (c 1 , . . . , c n , c ′ 1 , . . . , c ′ n ) T and d = (d 1 , . . . , d m )
T are solutions of the linear system

(K + nλW -1 )c + T d = y, (6)
T T c = 0, (7)
where

T = {L j φ ν } 2n m j = 1 ν = 1 , K = {L j(s) L k(t) E m (s, t)} 2n j,k=1 , W = diag(σ -2 x , . . . , σ -2 x , σ -2 v , . . . , σ -2 v ).
The solution F λ (t) is then a polynomial spline of order m with knots at the sampling time t 1 , . . . , t n . Note that the bivariate function E m (s, t), called semikernel, acts like a reproducing kernel in this approach, but its simple form is more appropriate for computational aspect and specifically for solving the linear system defined by Eqs. ( 6) and ( 7).

The error variances σ 2 x and σ 2 v are usually unknown in practice. In general, we use an estimator of the error variance corresponding to the criterion used for the selection of the smoothing parameter λ. Three scores are commonly used:

the UBR score ("Unbiased Risk") which is an extension of the Mallow's C p criterion ; -the GCV score ("Generalized Cross-Validation") which is a weighted version of the standard cross-validation ; -the GML score ("Generalized Maximum Likelihood") based on a Bayesian model, and that required a normality assumption on the errors.

The selection of the smoothing parameter results from the minimization of one of these criteria, and the error variances estimates depend on the smoothing parameter obtained, and therefore on the criterion chosen (see [START_REF] Gu | Smoothing spline ANOVA models[END_REF]. Thus, if in a first time, we consider only the position measurements, i.e. the data model (1), and if we denote A(λ) the hat matrix defined by

( F λ (t 1 ), . . . , F λ (t n )) T = A(λ)y, (8)
where F λ is the smoothing spline estimate of F for the smoothing parameter λ (that is a polynomial spline of order 2m), then the variance estimate of σ 2

x associated with the GCV criterion is

σ 2 gcv = y T (I -A(λ v )) 2 y tr(I -A(λ v )) , (9)
and the variance estimate associated with the GLM criterion is

σ 2 gml = y T (I -A(λ m ))y n -m . ( 10
)
Similarly, in a second time, we consider only the speed measurements and calculate the smoothing spline estimate of F ′ for the smoothing parameter λ, and then deduce an estimate of the variance error σ 2 v .

Smoothing under monotonicity constraint

The second constraint in the nonparametric model ( 1) is a monotonicity constraint since the function F representing the distance traveled as a function of time must be increasing. Various methods of smoothing under monotonicity constraint have been developed. The main approaches are based on kernel smoothers and splines. An overview of these methods can be found in Delecroix and Thomas-Agnan (2000). Among the main methods, we can cite the isotonic regression introduced by [START_REF] Brunk | Maximum likelihood estimates of monotone parameters[END_REF], the monotone splines (for example, Ramsay introduces the I-splines basis in [START_REF] Ramsay | Monotone regression splines in action[END_REF] for monotone regression splines) or the projection methods (e.g. Delecroix et al 1996 or [START_REF] Mammen | A general projection framework for constrained smoothing[END_REF].

In a previous study [START_REF] Sabatier Andrieu | Modélisation fonctionnelle et lissage sous contraintes d'un profil spatial de vitesse[END_REF], the method of homeomorphic splines developed by [START_REF] Bigot | Smoothing under diffeomorphic constraints with homeomorphic splines[END_REF] have been tested. However, if the monotonization step presented good results, we had difficulties in the implementation of the derivative. So, we propose to use a method developed by [START_REF] Ramsay | Estimating smooth monotone functions[END_REF] which has the advantage of being relatively simple to implement. The principle of this method is to transform the constrained smoothing problem to an unconstrained one. A monotone function has a positive first derivative. So the main idea is that any strictly monotonic function f satisfies the following differential equation:

D 2 f = wDf, (11)
where Df and D 2 f are respectively the first and second derivative of the function f , and w is an unconstrained function. Thus, any strictly monotonic function f can be written as following (as solution of the Eq. ( 11)):

f (t) = β 0 + β 1 t 0 exp[ u 0 w(v)dv]du, (12)
where β 0 and β 1 are arbitrary constants such that f (0) = β 0 and f ′ (0) = β 1 .

Then, the problem is to estimate the coefficients β 0 and β 1 and the unconstrained function w by minimizing the following criterion:

n i=1 (y i -β 0 -β 1 h(t i )) 2 + λ T 0 (w m (t)) 2 dt, (13)
where

h(t) = t 0 exp[ v 0 w(v)dv]du. (14)
The unconstrained function w is computed using an appropriate basis expansion (e.g. B-splines) and the coefficients β 0 and β 1 are estimated by numerical algorithms. However, due to a numerical optimization of the criterion (13), monotone smoothing spline involves considerably more computation than the usual smoothing spline process.

Thus, the monotonicity constraint is consider in a second smoothing step by applying the method of Ramsay described above, to the estimated values F λ (t i ) obtained at the first smoothing step for which we have used the speed measurements (Section 4.1). Therefore, this second smoothing step can be seen as a monotonization step and then is similar to the projection step in projection methods (see [START_REF] Mammen | A general projection framework for constrained smoothing[END_REF]. The new monotone estimator of F (t) is then denoted F mc (with "mc" for monotonicity constraint), and finally we deduce an estimator v S of the space-speed profile v S by the transformation

v S = F ′ mc • F -1 mc .

Simulation study

We propose to illustrate the performance of our smoothing procedure on three simulated examples. For each example, we study the monotone estimator F mc obtained after the two steps of the smoothing procedure (use of derivative information, and then monotonization), the derivative F ′ mc of this estimator, and the composite function

F ′ mc • F -1
mc . We investigate the regression model ( 2) with a fixed design made up of n = 50 points evenly distributed on [0, 1] for the first and second example, and n = 150 points evenly distributed on [0, 3] for the third example. The errors term ε x,i and ε v,i were simulated from centered gaussian distributions with σ x = 0.2 and σ v = 0.01. This choice of a smaller noise for the derivative of F than for F is motivated by the fact that speed measurements are generally more accurate than position measurements. The increasing regression functions chosen in these three examples are:

F 1 (t) = t 2 with t ∈ [0, 1], (15)
F 2 (t) = 1 2 (2t -1) 3 + 1 2 with t ∈ [0, 1], (16)
F 3 (t) =    (t -1) 3 + 1 si t ≤ 1 1 si 1 ≤ t ≤ 2 (t -2) 3 + 1 si t ≥ 2 with t ∈ [0, 3], (17)
and correspond respectively to a convex function, a function with a small plateau (inflection point), and a function with a large plateau (F 3 is an extension of F 2 with a plateau over [1,2]). The different functions, their derivatives, and the composite functions

F ′ i • F -1 i , i = 1, 2, 3, are displayed Fig.

(dashed blue lines).

The proposed estimator F mc of the the unknown regression function F for the three examples ( 15)-( 17) is the solid red line (left column in Fig. 5). As described in the previous section, this estimator was calculated in two steps. The first step corresponding to the construction of the estimator (5) that take into account the derivative information, was computed with the function ssr in the R (R Core Team 2013) package assist [START_REF] Wang | Assist: A suite of s functions implementing spline smoothing techniques[END_REF]). In the three examples, we have chosen m = 3 (i.e. a quintic spline) for the order of the Sobolev space, and the GML criterion was used for the selection of the smoothing parameter. The second step corresponding to the monotonization of the estimator (5) was computed with the function smooth.monotone in the R package fda. We have chosen m = 3 for the degree of the penalty in Eq. ( 13), and the smoothing parameter was chosen by trial and error.

The derivative F ′ mc of the estimator is represented in the middle column of the Fig. 5, and the composite function

F ′ mc • F -1
mc is in the right column (red lines). The computation of a point t 0 = F -1 (y 0 ) was made by the computation of the root (i.e. zero) of the function F (t)-y 0 (use of the R function uniroot). Due to computational problems of the inverse function at the edge of the interval, the composite function F ′ • F -1 and its estimator was computed over [0.1,0.9] for the first and second example ( 15)-( 16), and over [0.1,1.9] for the third example (17). We use 100 simulation runs to calculate the pointwise mean squared error (MSE) for the estimators F mc and F ′ mc evaluated on an equidistant grid of size 2n, and for the estimator

F ′ mc • F -1
mc evaluated on an equidistant grid with a step equal to 0.01. Curve of these MSE for the three examples are displayed Fig. 6. Table 1 Fig. 5 and Fig. 6 show the good results of the estimators F mc , F ′ mc and

F ′ mc • F -1
mc , even if we observe some boundary effects at the edge of the interval especially for the derivative estimator F ′ mc . As mentioned Sect. 3.2, the difficulty in the estimation of the composite function F ′ • F -1 is the estimation of the cusp when the function is equal to zero, and that corresponds to a plateau for the function F . The comparison between results obtained with the example (16) (regression function F 2 with a small plateau) and the example (17) (regression function F 3 with a large plateau) shows that the cusp is overestimate in the case of a small plateau (in the example ( 16), the estimator F ′ mc • F -1 mc does not cross the axis). This problem results from the mononization step, since the method of Ramsay which is used, provides a strictly increasing estimator F mc that is not appropriate for the estimation of plateaux even if in practice only small plateaux are not very good estimated. However, actually, we have not found a solution for this problem.

Application to the real data example

To illustrate the performance of the smoothing procedure in the estimation of space-speed profiles from noisy position and speed measurements, we applied the method on the real data set presented in Sect. 2. Thus, the smoothing procedure is applied to the data set of 78 individual space-speed profiles illustrated at Fig. 2 and reproduced Fig. 7.a, that are composed of noisy position and speed measurements, as follows:

1. A first smoothing step using derivative information, with for each path j, j = 1, . . . , 78, an estimation F λj (t) of each function F j (t) (representing the distance traveled as function of time) with the following parameters:

an estimation of the variance σ 2 x,j and σ 2 v,j for each path j, j = 1, . . . , 78; -m = 3 (quintic spline); -an automatic selection of each smoothing parameter λ j resulting from the minimization of the GML criterion. 2. A second smoothing step under monotonicity constraint, corresponding to a monotonization of each estimate F λj (t) obtained at the previous step with the following parameters:

m = 3 for the degree of the penalty ; -a selection of each smoothing parameter by trial and errors. Fig. 7.b illustrates the results of this smoothing procedure and shows the smooth individual space-speed profiles obtained with the transformation F ′ j,mc (t)• F -1 j,mc (t). Results are good since some peaks which appear in raw data, and that probably correspond to outliers, are reduced (e.g. blue and orange curves). Missing values are also be corrected by the smoothing procedure. The disadvantage of the step monotonisation mentioned in the simulation study, which causes an over-estimation of speed, appears mainly at the stop sign (short stop) but is less important at the traffic light (long stop). Finally, note that the main difference between Fig. 7.a and Fig. 7.b is that the smoothing procedure allows to reduce the study of these individual speed profiles to a functional framework since the estimated space-speed profiles belong to the space E SSP defined in Definition 1. Some advantages of studying a set of individual space-speed profiles in a functional framework are illustrated in the next sections.

Curve registration by landmarks alignment

The second step of our analysis is to summarize the information contained in a data set of individual space-speed profiles by the construction of the average profile. However, due to phase variation (i.e. horizontal variation) between the individual speed profiles, especially at stops (all vehicles do not stop at the same location), it is necessary to align them in order to obtain an aggregated profile which is representative of the set.

The curve registration or curve alignment problem appears in many areas such as biology, meteorology, pattern recognition... [START_REF] Ramsay | Curve registration[END_REF][START_REF] Bigot | Landmark-based registration of curves via the continuous wavelet transform[END_REF]). Indeed, frequently, observed curves exhibit two types of variability: amplitude variation which corresponds to vertical variation, and phase variation which corresponds to horizontal variation [START_REF] Ramsay | Functional Data Analysis, Second Edition[END_REF]. Then, to build a representative curve of a set of observed curves, it is necessary to correct the phase variation in order to obtain curves with similar features. The curve registration problem consists in finding, for each curve, a warping function and to deform all the curves in order to align them. If the literature about this problem is relatively large (e.g. [START_REF] Kneip | Statistical tools to analyze data representing a sample of curves[END_REF][START_REF] Wang | Alignment of curves by dynamic time warping[END_REF], it is not treated or treated in a simple way in velocity profiles studies.

We propose to use the method of landmarks alignment which consists to determine, for each curve, a deformation function so that specific points called landmarks of the registered curves are aligned. Specific points defined as landmarks are generally the positions of maxima, minima, inflection points, or zero crossings. Then, the landmarks registration of m signals f 1 , . . . , f m defined on the same interval [0, X] can be divided into the five following steps:

1. Definition of characteristic points to be used as landmarks (eg, minimum, maximum, zero crossing ...). 2. Extraction of landmarks x i,1 , . . . , x i,K from an observed sequence of each signal f i , i = 1 . . . , m. Note that since observed signals are noisy, the landmarks x i,1 , . . . , x i,K are usually extracted from a estimator f i of the signal f i . 3. Identify landmarks reference x 0,1 , . . . , x 0,K , i.e. the points at which the curves must match. 4. Determine deformation functions h 1 , . . . , h m so that corresponding landmarks are matched, i.e. for all i = 1, . . . , m, h i (x 0,j ) = x i,j , j = 1, . . . , K. 5. Deformation of the signals using transformations obtained in the previous step. The registered functions

f i (x) = f i [h -1 i (x)
], i = 1, . . . , m, are then aligned at each points x 0,1 , . . . , x 0,K .

The deformation functions h i (x), i = 1 . . . , m, called warping functions, must check the following properties:

-Initial conditions: h i (0) = 0, h i (X) = X. -Landmarks alignment: h i (x 0,j ) = x i,j .
-Strict monotonicity: x 1 < x 2 implies h i (x 1 ) < h i (x 2 ) (in order to respect the sequencing of points).

The method of landmarks alignment is applied to the set of speed profiles illustrated at Fig. 7.b. In order to compare similar speed profiles, we distinguish the two driving situations corresponding to the state of the traffic light (red or green light). Only the red light case will be studied in the following, that represents a sample of 36 individual profiles. We have chosen to define a. Unregistered space-speed profiles.

b. Registered space-speed profiles. landmarks as the positions of the two elements of the infrastructure that require a stop of the vehicle, namely the stop sign and the red light. Thus, the landmarks, corresponding to zero-crossing (or local minima) at the stop sign and the traffic light positions, are extracted from the estimated space-speed profiles obtained with the smoothing procedure, and are matched with the reference landmarks defined by the average position of vehicle stops at this two elements of the infrastructure. Then, monotone cubic spline interpolation have been determined as warping functions and have been computed with the R function splinefun and the option "monoH.FC". We also impose the condition that the warping functions are linear with a slope equal to one around the stops (we fix an interval of length 100 m around each stop) in order to not too distort the space-speed profiles in the neighborhood of each stop and to obtain "true" space-speed profiles as defined in Definition 1.

Fig. 8 compares the unregistered (Fig. 8.a) and the registered (Fig. 8.b) speed profiles in the red light case (36 curves). Fig. 8.a illustrates the fact that averaging unregistered profile results in an average profile (black curve) that is not representative of the set of the individual speed profiles. Indeed, this average profile doesn't equal to zero at the red light unlike all individual profiles. In contrast, Fig. 8.b shows that the average of the registered profiles tends to resemble much more closely most of the individual profiles, and then is a good aggregated speed profile of the sample.

6 Functional boxplot: a graphical tool to explore the variability of a functional data set Finally, the last step of our analysis is to explore the variability of a set of individual space-speed profiles. Indeed, if the construction of an aggregated speed profile, such as the average profile, leads to a good representation of the actual speeds on a road network section, such an aggregated profile does not reflect the variability between road users. The boxplot proposed by [START_REF] Tukey | Exploratory data analysis[END_REF] is a graphical method used to represent the distribution of univariate data, and can be used to represent speed variations between individuals at a given point. For example, Fig. 9.a represents pointwise boxplots calculated at a regular interval of 10 m in the red light case, with medians connected by a red line (V50 profile) and 85th percentiles connected by a blue line (V85 profile). However, this representation lost the continuous form of the individual profiles, and then the V50 and V85 profiles are not true space-speed profiles as defined in Definition 1 in contrast to the average speed profile obtained at Fig. 8.b.

So, we propose to use a graphical tool called functional boxplots, recently developed by Sun and Genton (2011), which extends the notion of boxplots to functional data. This tool is based on the notion of functional depth which generalizes order statistics or ranks to the functional setting. Indeed, the first step to construct a boxplot is the data ordering. But if the notion of order is obvious in the univariate setting, it is much more complicated in the functional setting. This problem has led to the emergence of the concept of functional depth, first introduced for multivariate data [START_REF] Zuo | General notions of statistical depth function[END_REF], that provides a measure of "centrality" and "outlyingness" for a function within a sample of curves and allows to order them from center-outward (López-Pintado and Romo 2009). The median curve is then the curve with the higher depth. Various examples of functional depth have been proposed in the literature such as the Fraiman and Muniz depth [START_REF] Fraiman | Trimmed means for functional data[END_REF], the random projection depth [START_REF] Cuevas | Robust estimation and classification for functional data via projection-based depth notions[END_REF] or the band depth (López-Pintado and Romo 2009). If Sun and Genton (2011) use the band depth and its modified version for the construction of its functional boxplots, a comparison of the results obtained with various functional depth led us to choose the h-mode depth introduced by [START_REF] Cuevas | A partial overview of the theory of statistics with functional data. J Stat Plan Inference Cuevas A[END_REF] and based on the concept of mode. The authors defined a functional mode as the curve most densely surrounded by the rest of curves of the dataset. Thus, the h-modal functional depth of a curve x i with respect the set of curves x 1 , . . . , x n is given by:

M D n (x i , h) = n k=1 K( x i -x k h ), (18)
where . is an appropriate norm, K is a kernel function, and h is a bandwidth. In practice, the L 2 norm and the truncated Gaussian kernel are used, and the bandwidth taken is the 15th percentile of the empirical distribution of { x i -x k , i, k = 1, . . . , n}. Functional boxplots create with the h-modal depth are illustrated in Fig. 9.b in the red light case. This functional boxplot is composed of the maximum enveloppe (blue curves), the median profile (black curve) which is the most central curve with the highest h-modal depth, the 25% central region (dark magenta region), the 50% central region (magenta region) and the 75% central region (pink region). The red dashed curves are the outlier candidates detected by the 1.5 times the 50% central region rule (see Sun and Genton 2011). This functional boxplot have been computed with the function fdepth of the R package rainbow and the function fbplot of the R package fda.

The advantage of this graphical tool is that it allows to represent the speed dispersion among individuals on a given road section. These speed corridors allow to distinguish road sections where the speed variability is large and those for which speeds are more homogeneous. Moreover, this tool leads to the extraction of the median profile which depends to the choice of a functional depth, and that can be used as a representative speed profile of the set of the individual speed profiles instead of the average profile.

Conclusion

In this paper, we have proposed a functional analysis of a set of space-speed profiles corresponding to speed as function of the distance traveled by the vehicle from an initial point. Thus, a definition of the functional space of these objects was proposed and the study of their mathematical properties has shown the remarkable property of non differentiability at points for which speed is zero and that corresponds to a cusp in the curve. Then, the first step of our analysis was the development of a smoothing procedure in order to be reduced to a functional framework. However, we have shown that the estimation of a space-speed profiles from noisy position and speed measurements was complex and can be reduced to a nonparametric regression problem taking into account two constraints: the use of the derivative information, and a monotonicity constraint. A two-step estimator (smooth, and then monotonize) based on the general theory of thin-plate spline was proposed and computed on simulation studies and on a real data set. If the proposed smoothing procedure presents good results, some limitations appear in the estimation of points for which the function is zero (speed tends to be overestimated at short stops). This point is an important challenge for the future, and the fusion of the two constraints in a single smoothing step will be subject to future research.

In a second time, a methodology has been proposed to summarize a set of individual space-speed profiles with an aggregated speed profile. The functional approach allows to use curve registration method in order to correct phase variation, and then to obtain a representative average speed profile with similar features of corresponding individual speed profiles. The method of landmarks alignment, which consists to align specific points of the curves, is applied on a data set where two driving situations corresponding to the state of the traffic light (red or green light) are distinguish. A comparison of the unregistered and the registered speed profiles at stops imposed by the infrastructure in the red light case, as well as the corresponding average profiles, illustrates the interest of the method. The development of unsupervised classification methods to distinguish traffic conditions (free vs congestion) or specific driving conditions (state of traffic lights) could also be subject to future works.

Finally, the variability of a set of individual space-speed profiles was explored by the use of functional boxplots, initially proposed by Sun and Genton (2011), which are an extension of the classical boxplots used in the univariate setting. This tool leads to the construction of speed corridors that reflect the variability between road users.

A Proofs of properties given in Section 3.2

Proof of Theorem 1: Let x 0 ∈ [0, x f ]. There are two distinct cases : 1 st case : x 0 is a point of continuity of F -1 . Then by composition of two continuous functions, we deduce that F ′ • F -1 is continuous at x 0 . 2 nd case : x 0 is a point of discontinuity of F -1 . We begin by demonstrating the following lemma:

Lemma 1 Let x 0 ∈ [0, x f ] a point of discontinuity of F -1 . Then the speed is zero at this point, i.e. v S (x 0 ) = F ′ • F -1 (x 0 ) = 0. This lemma can be proved easily. Indeed, if x 0 is a point of discontinuity of F -1 , then there is a close interval [t -, t + ] where F is constant and equal to x 0 , and by definition of F -1 , F -1 (x 0 ) = t -. This implies that F ′ + (t -) = 0 where F ′ + (t -) is the right derivative of F at t -, and as it was assumed that F was differentiable, we also have F ′ -(t -) = 0 where F ′ -(t -) is the left derivative of F at t -. Finally, F ′ (t -) = 0, and therefore F ′ • F -1 (x 0 ) = 0 which ends the proof of the lemma 1. Now, we study the one-sided limit of v S = F ′ • F -1 at x 0 . (v S ) -(x 0 ) = lim F ′ (t) = 0 since F ′ (t) = 0 on [t -, t + ]. So, we deduce that (v S ) + (x 0 ) = 0. Hence, using Lemma 1, we conclude that v S is continuous at x 0 .

⊓ ⊔ Proof of Theorem 2:

1 st case: Assume that F satisfies the assumptions (H 1 ).

Let x 0 such that t 0 = F -1 (x 0 ). Since F ′ (t 0 ) = 0, then v S (x 0 ) = 0, i.e. x 0 ∈ H 0 . Under the assumptions (H 1 ), we can apply the Taylor-Young's formula to F ′ : For all θ in a neighborhood of t 0 , F ′ (t 0 +θ) = F ′ (t 0 )+θF ′′ (t 0 )+ θ 2 2 F ′′′ (t 0 )+θ 2 ε(θ), where ε(θ) → 0 when θ → 0. But since F ′ (t 0 ) = 0, if we had F ′′ (t 0 ) = 0, then F ′ would change sign at t 0 , which contradicts the strict monotonicity of F . Therefore F ′′ 0 ) = 0. So, F ′ (t 0 + θ) ∼ θ→0 θ 2 2 F ′′′ (t 0 ) (since it is assumed that F ′′′ (t 0 ) = 0). Let h = F (t 0 + θ) -F (t 0 ). We apply the Taylor-Young's formula to F : h = F (t 0 + θ) -F (t 0 ) = θF ′ (t 0 ) + θ 2 2 F ′′ (t 0 ) + θ 3 6 F ′′′ (t 0 ) + θ 3 ε ′ (θ) where ε ′ (θ) → 0 when θ → 0. In order to study the differentiability of vs at x 0 , we define the following growth rates:

v S (x 0 +h)-v S (x 0 ) h = F ′ (t 0 +θ)-F ′ (t 0 ) F (t 0 +θ)-F (t 0 ) ∼ θ→0 θ 2 2 F ′′′ (t 0 ) θ 3 6 F ′′′ (t 0 )
= 3 θ . This growth rate has no limit when θ → 0, but this does not prove that it has also no limit when h → 0. We will prove this by contradiction. Assume that -ℓ| < ε. This means that the growth rate has a limit ℓ ∈ R when θ → 0, which is a contradiction. Hence, under the assumptions (H 1 ), v S is not differentiable at x 0 .

2 nd case: Assume that F satisfies the assumptions (H 2 ). As in the first case, we define x 0 such that t 0 = F -1 (x 0 ). The graph of G :

coincides with F on [0, t 0 ], -is deduced from the graph of F by the translation vector (t 0t 1 )

-→ i on [t 0 , T -(t 1t 0 )].

Thus, the graph of G is similar to the graph of F but removing the time period [t 0 , t 1 ] for which the function is constant. So, the same growth rate occurs at x 0 , and if v S is not differentiable at x 0 for one, it is not for the other. In other words, the results of the first case where F ′ = 0 at one point t 0 extend to the more general case where F ′ is zero on an interval [t 0 , t 1 ] (t 0 = t 1 ), subject to the assumptions (H 1 ) on G. ⊓ ⊔

Fig. 1

 1 Fig. 1 Map of the studied section.

Fig. 2

 2 Fig. 2 Sample of 78 space-speed profiles : raw data.

Fig. 3

 3 Fig. 3 Link between the three study areas : [distance × time, speed × time] and [speed × distance].

Fig. 4

 4 Fig. 4 Functional diagram illustrating the definition of space-speed profiles.

Fig. 5

 5 Fig. 5 Estimators of F , F ′ and F ′ • F -1 (left to right) for the three examples (15)-(17) (top to bottom). The unknown regression function F , its derivative F ′ and the composite function F ′ • F -1 are dashed blue lines. The noisy data are green points. The estimators Fmc, F ′ mc and F ′ mc • F -1 mc are solid red lines.

Fig. 6

 6 Fig. 6 Simulated mean squared error of the estimators Fmc, F ′ mc and F ′ mc • F -1 mc (left to right) computed over 100 simulations runs on the three examples (15)-(17) (top to bottom).

Fig. 7

 7 Fig. 7 Smoothing step on the 78 individual space-speed profiles.

Fig. 8

 8 Fig.8Registration of space-speed profiles in the red light case (36 curves). The black curve is the average profile.

Fig. 9

 9 Fig.9Pointwise boxplots ans functional boxplots in the red light case (36 curves).

FF

  ′ • F -1 (x). When x → x 0 by lower values, t → t -by lower values, and limt→t - t<t - F ′ (t) = 0 since F ′ (t) = 0 on [t -, t + ] and F ′ is continuous at t -. So, we deduce that (v S ) -(x 0 ) = 0. Similarly, (v S ) + (x 0 ) = lim ′ • F -1 (x). When x → x 0 by upper values, t → t -by upper values, and lim t→t - t>t -

  v S (x 0 +h)-v S (x 0 ) h → h→0 ℓ ∈ R. Then, by definition, ∀ε > 0, ∃α > 0 such that |h| < α ⇒ | v S (x 0 +h)-v S (x 0 ) h -ℓ| < ε. But since F is continuous at t 0 , ∃β > 0 such that |tt 0 | < β ⇒ |F (t) -F (t 0 )| < α, or similarly |θ| < β ⇒ | F (t 0 + θ) -F (t 0 ) h | < α. Hence, ∀ε > 0, ∃β > 0 such that |θ| < β ⇒ | v S (x 0 +h)-v S (x 0 ) h

  shows the results of the mean integrated squared error (MISE) for the estimators F mc , F ′ mc and F

							MSE of the estimator of F(t)											MSE of the estimator of F'(t)							MSE of the estimator of F'(F^(-1)(x))
		0.0013																0.06															0.006
		0.0012																0.05															0.005
		0.0011																0.04															0.004
	MSE	0.0010															MSE	0.03															MSE
		0.0009																0.02															0.003
		0.0008																0.01															0.002
		0.0007																0.00															0.001
		0.0			0.2			0.4			0.6			0.8			1.0	0.0			0.2			0.4			0.6			0.8			1.0	0.1	0.3		0.5	0.7		0.9
									t																t												x
							MSE of the estimator of F(t)											MSE of the estimator of F'(t)							MSE of the estimator of F'(F^(-1)(x))
		0.00085																0.025
																																	0.12
																		0.020
																																	0.10
		0.00084																0.015															0.08
	MSE																MSE																MSE
																		0.010															0.06
		0.00083																0.005															0.04
																		0.000															0.02
		0.0			0.2			0.4			0.6			0.8			1.0	0.0			0.2			0.4			0.6			0.8			1.0	0.2		0.4	0.6		0.8
									t																t												x
							MSE of the estimator of F(t)											MSE of the estimator of F'(t)
		6e-04																0.14															0.05
																		0.12															0.04
		5e-04																0.10
																		0.08															0.03
	MSE	4e-04															MSE
																		0.06															0.02
		3e-04																0.04															0.01
																		0.02
		2e-04																0.00															0.00
		0.0	0.2	0.4	0.6	0.8	1.0	1.2	1.4	1.6	1.8	2.0	2.2	2.4	2.6	2.8	3.0	0.0	0.2	0.4	0.6	0.8	1.0	1.2	1.4	1.6	1.8	2.0	2.2	2.4	2.6	2.8	3.0	0.1	0.4	0.7	1.0	1.3	1.6	1.9
									t																t												x

′ mc • F -1 mc in the three examples. MSE MSE of the estimator of F'(F^(-1)(x))

Table 1

 1 Mean integrated squared error (MISE) of the estimators Fmc, F ′ mc and F

	′ mc • F -1 mc ,

Acknowledgements The authors would like to gratefully thank Christine Thomas-Agnan from the University of Toulouse I for discussions about RKHS and the general smoothing problem, and for her help in the use of the derivative information in the smoothing procedure presented in this article. We thank Yuedong Wang from the University of California for his help in the implementation of the smoothing procedure with the R package assist, and Jeremie Bigot from the Toulouse Mathematics Institute for discussions about the present study. This work was funded by the French Institute of Science and Technology for Transport, Development and Networks (IFSTTAR).