Computing travel times from filtered traffic states

This article experimentally assesses the influence of sensor data rates on travel time estimates computed from filtered traffic speed estimates. Using velocity data obtained from GPS smartphones and inductive loop detector data collected during the Mobile Century experiment near Berkeley, CA, and an evolution equation for average velocity along the roadway, an estimate of the traffic state is obtained via ensemble Kalman filtering. A large--scale batch of computations is run to produce estimates of traffic velocity with varying degrees of input data, and instantaneous and a posteriori dynamic travel times are compared to travel times recorded using license plate re-identification. We illustrate that dynamic travel time estimates can be computed with less than 10% error regardless of the data source, and that existing inductive loop detector data can significantly improve the accuracy of travel time estimates when GPS data is sparse.

Introduction.
1.1. Objective. Probe data will likely become ubiquitous in the not too distant future, due to the rapid expansion of consumer generated probe data from cellphones, personal navigation devices, and intelligent vehicles. As the use of probe data for traffic monitoring increases, so does the need to understand the benefits and trade-offs between GPS data and conventional data sources. Yet, a complete analysis of the trade-offs between probe data and fixed sensors is difficult, because the value of the data from any sensor (GPS equipped probe vehicles, inductive loop detectors, etc.), is dependent on the specifics of the sensing technology, the method used to process the data, and the specific traffic monitoring objective in question.
As navigation devices in vehicles become more common, either on smartphones or through integrated in-vehicle systems, the need to estimate travel times between any two arbitrary points on the network in real-time will become increasingly important. One approach to solving this problem is to construct explicit travel time estimators for all origin destination pairs on the network. Obviously, this approach quickly becomes intractable as the size of the network grows to scales relevant to commercial traffic monitoring companies. A computationally efficient alternative is to estimate the traffic speed throughout the network. Then, with the speeds on each link known, travel time estimates can be computed as requested by the navigation devices.
The latter approach to travel time estimation motivates our work, which is to study the following question. To what degree can GPS probe data act as a substitute for conventional traffic monitoring technologies such as inductive loop detectors, for the purpose of estimating travel times? Specifically, we focus our attention on estimating travel times by integrating various data volumes from inductive loop detectors and GPS equipped probe vehicles into a velocity flow model equivalent to the Cell Transmission Model [11,12], using an estimation technique known as ensemble Kalman filtering (EnKF). A preliminary version of this work was presented in conference form in [23].
This article emphasizes the experimental performance of a flow-model based estimator using data collected during a one day field experiment known as Mobile Century [17]. The data set collected during this experiment is unique because of the large number of GPS equipped probe vehicles representing 2-5% of the traffic flow, the dense coverage of working inductive loop detectors on the experiment site, and the availability of travel time data obtained from video license plate reidentification. Thus, although the results presented in this article are still limited in geographic scope and in time, to the best of the authors' knowledge, they are based on the most comprehensive publicly available GPS data set to date [30].
1.2. Related work. Several field experiments have been conducted to assess the applicability of cell phone-based measurements for traffic monitoring [1,2,3,8,22,25,28,29], including data generated from cell phone towers, which produces less accurate vehicle position and speed measurements compared to GPS. Bar-Gera [3] compared several months of network data from cellphones to inductive loop detector data on a 14 km freeway segment in Israel, and found them to be in good agreement. Liu et al. [22] evaluated a different network-based cell phone system in Minnesota, and compared travel times to license plate re-identification, and found the system generated results with varying accuracies. A summary of the major network-based cell phone experiments to date can be found in Liu et al. [22].
Several studies have also been conducted to assess the trade-offs between inductive loop detector data and data collected from GPS equipped probe vehicles. In Kwon et al. [20], it is shown that annual estimates of total delay, average duration of congestion, and average spatial extent of congestion can be made with less than 10% error by using either inductive loop detectors placed with half-mile spacing, or by using probe vehicle runs at a rate of about three vehicles an hour. Approximately four to six days of data is needed for reliable estimates from either data source.
The work of Herrera et al. [16] compares a nudging algorithm and a mixture Kalman filtering algorithm to examine how the addition of probe vehicle measurements sampled at a fixed time interval can decrease errors in estimating traffic velocity. On a 0.4 mile stretch of roadway, sampling 5% the traffic at 150 second intervals with inductive loops at both ends of the domain leads to a 16% improvement over the inductive loop detector data alone. The article also uses the Mobile Century experiment data to compare three scenarios of time-based sampling of probe vehicles, finding that probe data outperforms inductive loop detector data for estimating traffic velocity if a sufficient number of measurements can be obtained from probe vehicles. This work uses the same data set from Mobile Century, but we now consider nearly one thousand scenarios to compare probe data to inductive loop detector data.
The works [4,7,10] are also closely aligned with the present work. In Cristiani et al. [10], experimental mobile sensor data is used to calibrate a flux function and to determine an initial condition for a traffic flow model instantiated for a stretch of highway in Rome, Italy. The model is then simulated forward in time to make accurate travel time estimates. Blandin et al. [4] use forward simulation with boundary data obtained from sensors to compare estimates of the traffic state from a calibrated scalar traffic model with a calibrated 2 × 2 phase transition model [9]. Instead of open loop simulation with a calibrated model, the traffic prediction in [7] is updated by assimilating historical data in the forecast.
1.3. Methodology overview and organization of the article. In order to assess the trade-offs between velocity data collected from GPS smartphones and velocity data obtained from inductive loop detectors, it is necessary to define the process by which the data is transformed into an estimate of travel time. In this article, we rely on a velocity estimation algorithm developed as part of the Mobile Millennium project [30]. The algorithm combines velocity measurements from GPS smartphones or inductive loop detectors with a model of traffic evolution, using a technique known as ensemble Kalman filtering (EnKF) to produce an improved estimate of the velocity field, from which the travel time is computed. The resulting travel time computed from this process is then compared to the travel times recorded from the license plate re-identification video data.
With the data processing algorithm determined, we create a number of scenarios in which the volume of probe data and number of inductive loop detectors made available to the processing algorithm are adjusted. For example, this allows us to compare the accuracy of computing travel times when all of the probe data is made available, to travel times which are computed when only some of the probe data is available, to travel times when some probe data is available and some inductive loop detector data is available. In this way, we can quantify the trade-offs of various amounts of data from probes and inductive loop detector data in terms of increased or decreased accuracy of the computed travel times.
In order to describe and quantify what probe data is made available to the travel time processing algorithm, we detail two metrics of importance to probe data, namely the penetration rate and the sampling rate. The penetration rate is defined as the percentage of cars on the roadway reporting probe data compared to the overall traffic flow, including the vehicles which do not send data. In addition to increasing the number of measurements, as the penetration rate increases, the sample of vehicles which generate measurements are more likely to be representative of the total traffic flow. The sampling rate refers to the frequency at which data is collected from the probe vehicles, and can be used to increase or decrease the number of measurements made available for estimating travel times from the same vehicles. The probe data collection technique used in this work collects data from probe vehicles at fixed points in space using a technique known as Virtual Trip Lines (VTLs) [18] invented by Nokia. By decreasing spacing between the VTLs, the probe vehicles will send more measurements, with smaller spacing between measurements.
In order to modify the amount of data obtained from inductive loop detectors, the number of inductive loop detectors which are made available to the processing algorithm is adjusted. Because this article focuses on a real highway, it is not possible to modify the location of the inductive loop detectors. Instead, given a fixed number of inductive loop detectors to include for a given scenario, we select the specific loop detectors such that they achieve as uniform of a spacing along the highway as is possible.
The work presented in this article summarizes the findings of a large scale computational study in which an ensemble Kalman filtering estimation algorithm is used to produce traffic estimates, and to characterize the dependency of the solution in the amount of sensing data used for the computation. The contribution of the article is thus the method used to assess the potential gains provided by probe and inductive loop data.
The remainder of this article is organized as follows. The key features of the processing algorithm used for velocity estimation are given in Section 2, and the methods for computing travel times from the velocity field are described. The data collected from the Mobile Century experiment is detailed in Section 3.1, and in Section 3.2, the techniques for generating scenarios with various amounts of input data from probe vehicles and inductive loops are presented. In Section 4, the results of nearly one thousand scenarios using various amounts of inductive loop detector data and probe data for travel time estimation are analyzed. Finally, the discussion in Section 5 summarizes the results.
2. Algorithm for estimating travel times. The processing algorithm used in this work is based on a velocity estimation algorithm developed in the Mobile Millennium system. The algorithm takes velocity data from inductive loop detectors and probe vehicles as input, combines the data with a physical model of traffic evolution, and produces an improved estimate of the velocity along the corresponding stretch of roadway. Using this improved estimate of velocity, an estimated travel time is computed using an instantaneous method and a dynamic method, to compare against the travel times recorded from video data. A brief overview of this process is described in this section.
2.1. Mobile Millennium highway traffic estimation algorithm. The velocity estimation algorithm developed in the Mobile Millennium system is based on a discretization of a traffic flow model known as the Lighthill-Whitham-Richards (LWR) partial differential equation [21,24] which describes the evolution of traffic density on the highway. In its discrete form, this model is also known as the Cell Transmission Model [11,12]. In order to simplify the velocity estimation problem, this model is transformed into an equivalent velocity evolution equation [27] operating on a 30 second time step. The complete mathematical details of the employed traffic velocity evolution equation, and the fusion of velocity measurement data with the evolution equation using ensemble Kalman filtering (EnKF) are presented in [27]. We give an overview of the model and filtering algorithm next.
2.1.1. Velocity traffic dynamics. The seminal LWR equation, as proposed in [21,24], reads: ∂ρ(x, t) ∂t where ρ(x, t) and q(x, t) respectively denote the density and flow at location x and time t. Additionally, we denote v(x, t) the velocity field on the highway, and assume the velocity can be specified as a function V of the density only. This allows the flux q to be defined as a function of the density: The widely used triangular flux function assumes a constant velocity in free-flow and a hyperbolic velocity in congestion: where v max , ρ max , ρ c and w f are respectively the maximum velocity, maximum density, critical density at which the flow transitions from free-flow to congested, and the backwards propagating wave speed, respectively. Because the triangular velocity function is not strictly monotonic in free-flow, it cannot be inverted.
In order to use the triangular model in a velocity setting, we approximate it by the Smulders velocity function [26], with a linear expression in free-flow and a hyperbolic expression in congestion: For continuity of the flux at the critical density ρ c , the additional relation ρc ρmax = w f vmax must be satisfied. The Smulders velocity function (3) can be inverted to obtain the velocity as a function of density: where v c is the critical velocity: v c = V (ρ c ). The LWR PDE is discretized using a Godunov numerical scheme to obtain a discrete density evolution equation. After discretization, we can apply a variable change on the densities using the inverse velocity function V −1 . Let K and i max be two positive integers, we discretize time (indexed by k) and space (indexed by i) in K time steps of length ∆T = T K and i max + 1 space cells of length ∆x = L imax+1 . Then, according to the Godunov scheme, the velocity value v k+1,i of cell i at time step k + 1 can be computed as: where g (v 1 , v 2 ) represents the numerical flux between consecutive cells with respective velocities v 1 and v 2 . In the case of the Smulders model, we obtain: We note that the evolution of the velocity field at each discrete point on an edge is well defined by (5) and (6), except at the boundary points v k,0 and v k,imax . At these boundaries, the equations contain references to the ghost cells v k,−1 and v k,imax+1 , which are points which do not lie in the physical domain. The values of v k,−1 and v k,imax+1 are given by the prescribed boundary conditions to be imposed on the left and right side of the domain respectively. Extension of the model to networks follows a standard approach [15], and is detailed in [27].

Traffic state estimation. Given the velocity field
T on all space cells at time step k, the velocity field at time k + 1 can be obtained by applying: where f represents the update algorithm as described in (5), (6), (7) and (8). The term w k ∼ (0, Q) is the white, zero mean model noise with covariance Q. At every time step k, the measurements y k are related to the traffic state through the observation equation: where H k ∈ {0, 1} z k ×(imax+1) encodes the location of the z k data sources that sent measurements during that time interval, and ν k represents the measurement errors, which are assumed to follow a zero mean distribution with a measurement error covariance matrix R k . Note that H k depends on the time due to the fact that the location where measurements are recorded changes over time [27]. The recursive state estimation problem for linear time invariant systems is commonly solved using the Kalman filtering [19] algorithm. New measurements can be integrated optimally at every time step, using only the estimate at the previous time step. An a priori estimate at time k given data through time k − 1, denoted v k|k−1 , is computed by evolving the a posteriori estimate v k−1|k−1 through a linear evolution equation. Then, measurements at time k are collected and used to construct a corrected estimate v k|k by minimizing the trace of the posterior error covariance matrix P k|k .
Due to the nonlinearity of (9), a Kalman filter cannot be used, and an ensemble Kalman filter (EnKF) approach is chosen instead, as described in [13] and [14]. The EnKF circumvents the need for linearizing the evolution equation to fit the Kalman filtering framework, and instead approximates the evolution of the error covariance matrix P k|k through an ensemble (sample) approximation. The algorithm is detailed in Algorithm 1.
A few remarks on the performance of the velocity estimation algorithm described above are in order. First, it is noted that the algorithm was designed as part of the Mobile Millennium system, where it is not possible to track probe vehicles for privacy reasons, and thus no continuous GPS records from probes are assumed to Algorithm 1 Ensemble Kalman filter 1: Initialization: Draw J ensemble realizations v j 0|0 (with j ∈ {1, · · · , J}) from a process with a meanv 0|0 and covariance P 0|0 . 2: Model prediction: Update each of the J ensemble members according to the CTM-v forward simulation algorithm f : Then update the ensemble meanv k|k−1 and covariance P k|k−1 according to: 3: Measurement update: Obtain measurements, compute the Kalman gain G, and update the estimate of the state: 4: Return to 2.
be available for the estimation algorithm. In practice, it is expected that the performance of the estimation algorithm could be improved when tracking of individual probe vehicles is allowed. Second, it should be noted that the flow model requires some historical flow information to calibrate the model. In this study, historical inductive loop detector data from PeMS was used to estimate constant mainline inflow and outflows, and thus all results presented use inductive loop data in this way.
Next, the methods for computing the instantaneous and dynamic travel times from an estimated velocity field are described.

2.2.
Methods for computing travel times. The position of a vehicle p(t) moving at the average speed of traffic evolves according to: where v(x, t) = V (ρ (x, t)) . Thus, the dynamic travel time τ dyn (t 0 ) of a vehicle starting a trip at time t 0 at position p 0 is computed by solving (16) for the vehicle position, then solving: where L is the distance of the trip. Note that in general, the dynamic travel time may be infinite if the velocity drops to zero and remains zero indefinitely, although this is not an issue in practice.
A limitation of (16) is that it relies on knowledge of the velocity field up to τ dyn into the future. This means that the velocity field requires a short term estimate into the future, which further introduces computational overhead. Moreover, when the future forecast is an open loop model prediction without data to correct the estimate, the error in the velocity estimate increases. It is also worth noting that the dynamic travel time for individual vehicles is measurable, for example via license plate re-identification.
As an alternative, a common simplifying assumption is to assume the velocity field does not evolve during the travel time computation. This simplification gives rise to the instantaneous travel time, τ inst . To compute the instantaneous travel time, the position of the vehicle starting at p 0 at time t 0 is assumed to evolve according to: Then, the instantaneous travel time can be computed by solving (17) for the vehicle position, then solving: The instantaneous travel time has several advantages, most notably that it does not require the overhead of computing the velocity field forward in time, which can be intensive for large network predictions with an accurate discretization. This comes at the cost of decreased accuracy when the traffic state is changing quickly. Moreover, the issue of an infinite travel time is a more serious issue for the instantaneous travel time, as it can be generated if at any point the velocity drops to zero. To keep the estimates bounded, a minimum velocity greater than zero can be used during the integration of the vehicle's position.
Since the velocity field is the a posteriori filtered estimate, the instantaneous and dynamic travel times are approximated by integrating through the discrete velocity field v k|k,i . Specifically, the velocity field in (16) or (17) is approximated by: where i = p(t)/∆x with · denoting the floor operator. When computing the dynamic travel time, the time index k is updated according to k = t/∆T , and fixed as k = t 0 /∆T for the instantaneous travel time.
Because the estimated velocity field is piecewise constant in each cell over each time step, the solution of (16) or (17) can be computed exactly. Assuming the current position of the vehicle is in cell i, the vehicle updates its position each ∆T by computing the time τ required to reach the boundary between cells i and i + 1 when traveling at the speed v k|k,i . If τ < ∆T , the vehicle reaches the boundary within the time step. The vehicle travels at the speed v k|k,i until the boundary is reached, and travels at the speed v k|k,i+1 afterward. If the boundary is not reached, the vehicle travels at v k|k,i for the full time step. The position update is summarized as After ∆T , the position of the vehicle is updated, k is incremented, the new cell index i is computed, and the process is repeated.
Note that higher resolution numerical reconstructions are available for the dynamic travel time reconstructions [5] if discretization errors become significant. 3. Experimental data.
3.1. The Mobile Century experiment. In this section, the key features of the data collected during the Mobile Century experiment are presented. For a complete description of the experiment, the interested reader is referred to [17]. The Mobile Century field experiment was a one-day test in the San Francisco Bay Area which collected GPS data from cell phones in probe vehicles, inductive loop detector data, and travel time data from license plate re-identification video data. The experiment took place on February 8th, 2008, and involved 100 probe vehicles equipped with Nokia N95 cell phones which repeatedly drove a stretch of the I-880 freeway near Union City, CA, generating 2,200 vehicle trajectories.
The experiment site is also covered with 17 working inductive loop detector stations which feed measurements into the PeMS system [6]. The inductive loop detectors record the sensor occupancy and vehicle counts every 30 seconds, which is processed by a Mobile Millennium filtering algorithm in order to obtain the 30 second average velocity at the sensor. At 5 minute intervals, the PeMS system produces an estimate of the 5 minute average velocity at the sensor, which is shown in Figure 1b for the northbound traffic. The locations of the inductive loop detector stations are shown in Figure 1a.
Finally, as part of the experiment, high definition video cameras were temporarily installed to record license plates of northbound traffic. The travel times recorded from the re-identified vehicles traveling northbound is shown in Figure 2 most re-identified drivers experienced heavy evening congestion with travel times increasing to 15-20 minutes.

3.2.
Algorithms for data selection. To assess the trade-offs between different amounts of probe data and inductive loop detector data for the purpose of estimating travel times, we algorithmically select different subsets of loop and GPS probe data from the Mobile Century experiment [17], and use these subsets as inputs to the estimation process described in the previous section. We now describe the scenarios which modify the type and amount of the data which is made available for estimation, and the selection criteria which are used to generate the scenarios.

Selection of inductive loop detector data.
In order to modify the number of inductive loop detector stations which are made available for computing travel times, a simple selection criterion is developed. Specifically, given a fixed number of stations to include, the loops are selected in order to minimize the variance of the distance between consecutive sensors, resulting in approximately uniformly spaced sensors.
We consider a stretch of highway of length L, starting at x = 0 and ending at x = L, with n inductive loop detector stations located at x 1 , x 2 , · · · , x n , as shown in Figure 3. Let S i denote the spacing between sensor i and i + 1. In order to treat the boundaries without explicit knowledge of sensors outside the domain x ∈ [0, L], it is assumed that only half of the first inter-station spacing S 0 and of the last interstation spacing S n lie in the domain of interest. The weighted average spacing between the n sensors is given by: where the first and last spacings have a weight 1 2 , since only half of these spacings actually lie within the [0, L] domain. Note that the average spacing is independent of the specific locations of the sensors x i , and consequently cannot be used as a selection criterion.
Instead, we use a selection criterion which explicitly takes the uniformity of the inter-station distances S i into account. This is achieved by minimizing the variance σ 2 of the inter-station spacings S i , given by: Again, the first and last spacings have a weight 1 2 , since only half of these spacings actually lie within the [0, L] domain.
In practice, rather than minimizing the variance σ 2 , it is convenient to minimize an equivalent loop detector placement criterion denotedS: which shares the same minimizer as σ 2 . The best set of m inductive loop detector stations is then given by: where |U | represents the number of elements in the set U . The resulting selections for the inductive loop detector stations are shown in Table 1.
In the case when the chosen inductive loop detector stations are uniformly spaced within the section of interest, the criterionS is equal to the average spacingS. BecauseS serves as a lower bound forS, the difference betweenS andS indicates the degree of non uniformity of the sensor spacings caused by the fixed set from which the sensors are selected. Table 1 shows the difference between the inductive loop detector placement criterionS(U * (m)) and its lower bound, the average inductive loop detector spacingS(m), is small, indicating that the sensor spacing is relatively uniform.

Penetration rate for probe data.
A few remarks about probe penetration rates are important to make, before criteria to modify the penetration rate are discussed. In general, the penetration rate is difficult to determine for probe vehicles specifically because it depends on the number of equipped probe vehicles, the total traffic flow, and the evolution of the traffic flow in space and time. Typically, only the total number of equipped probe vehicles is known to probe data providers. Similarly, the total traffic flow can only be estimated from counts recorded by inductive loop detectors at predefined locations. Finally, because the evolution of the traffic flow Table 1. Inductive loop detector selection results. Given a number m, the selection algorithm returns the set U * (m) of m inductive loop detector stations which minimizes the inductive loop detector placement indexS(U * (m)). The labels in U * (m) correspond to the labels of the inductive loop detectors in Figure 1a. is not under the control of the probe vehicles, it is nearly impossible to a priori specify a penetration rate which is both uniform in space, and uniform in time.
Because of the inherent difficulty in specifying the penetration rate a priori, we instead elect to directly modify the number of equipped probe vehicles as a proxy for modifying the penetration rate. The number of equipped probe vehicles in this study varies from 0% to 100% of the 2,200 Mobile Century vehicle trajectories, increasing by increments of 10%. Over the eight hour experiment, this corresponds to an average rate of probe vehicles between 27.5 veh/hr and 275 veh/hr. At a probe rate of 275 veh/hr, the 20 minute average penetration rate at the center of the experiment site ranges between 1.5% and 3% (Figure 4b). When a subset of the vehicle trajectories is required, the subset is determined by selecting the trajectories at random before the simulation.

3.2.3.
Space-based sampling. In order to modify the number of measurements used from each probe vehicle trajectory under spatial sampling, the number of locations where measurements are collected are modified. The locations where measurements are obtained are encoded through the placement of virtual trip lines (VTLs), which can be viewed as virtual geographic markers which trigger vehicles to send measurements when the vehicle trajectory intersects the VTL. A complete description of the VTL sampling strategy is described in detail in Hoh et al. [18].
Because the VTLs are virtual, it is possible to place them anywhere on the experiment site. The determination of optimal VTL placement is complex, so instead we elect to place the VTLs uniformly across the experiment site. The number of VTLs, denoted by n VTL , tested in our scenarios varies from 9 VTLs to 99 VTLs, increasing by increments of 10 VTLs. This corresponds to an average spacing between 0.79 to 8.68 VTL/mi. Note the number of VTLs used on the experiment site is significantly higher than the number of inductive loop detector stations. This is possible because unlike inductive loop detector stations, the marginal cost of virtual trip lines is small. The number of probe vehicle measurements used in each simulation is shown in Figure 4a.  3.3. Summary of scenarios considered. Below, we give a summary of the various combinations of input data used for computing travel times in this work.
• Number of inductive loop detectors. Nine sets of inductive loop detector data, ranging from scenarios with zero to 16 stations, increasing by increments of two. • Number of probe data measurements. The amount of probe data is modified in two ways.
-Penetration rate. Eleven penetration rates are considered, ranging from scenarios with no probe data, to scenarios when 100% of the 2,200 probe vehicle trajectories are used, increasing by increments of 10%. This corresponds to an average rate of probe vehicles between 27.5 veh/hr and 275 veh/hr. -Number of measurements per vehicle. Ten sets of VTL locations are considered, ranging from scenarios with 9 (0.79 VTL/mi) to 99 VTLs (8.68 VTL/mi), increasing by increments of 10 VTLs. In total, 917 scenarios are created by instantiating scenarios with all combinations of the 9 sets of inductive loop detector data sets, 11 probe penetration rates, and the 10 VTL sets. In the remainder of this section, we describe the specific algorithms which select the data for each scenario. The scenarios tested are summarized in Table 2.  4. Results and discussion. In this section, we present the results of the 917 runs with varying amounts of probe and inductive loop detector data. We also vary the type of travel time computed (instantaneous or dynamic). First, the quantification of error is described, then the computational results are presented.

Error metric.
Because validation data is available for dynamic travel times via license plate re-identification, an error metric is used to compare the velocity estimation algorithm output that has been converted to travel times with the travel time measured from video recordings. By using the travel time error as a performance metric, estimation algorithm results can be compared with the results obtained when using different types and quantities of the input data.
Since the license plate re-identification data provides a distribution of individual vehicle travel times (see Figure 2), we define the true travel time as a one minute moving average of the recorded travel times. Figure 2 also shows the division of the experiment into four time periods that represent the different phases of the traffic during the experiment. These periods are (i ) the morning accident, were travel times are decreasing as an incident clears, (ii ) a free flow period during the middle of the day when travel times are low, (iii ) a congestion building period before the evening rush hours, and (iv ) full congestion during the evening rush hours. Because of the different traffic conditions present in these time intervals, in addition to computing the error across the full day, the error is also computed for each time interval.
The travel time error is computed as follows. Let K be the number of estimates given in a period for which the error is to be computed, with each estimate indexed by k. Let τ v (t k ) be the mean travel time from the video data for vehicles entering at time t k , and let τ inst (t k ) and τ dyn (t k ) be the estimated mean travel time computed with the instantaneous and dynamic methods at time t k , respectively. The mean absolute percent error (MAPE) for the travel time computed with the instantaneous method is: while the MAPE for the travel time computed with the dynamic method is computed similarly.
4.2.1. Implementation. The estimates were computed using the existing Mobile Millennium [27] highway model summarized in Section 2, and implemented in the Java programming language. The model was run 917 times with various data inputs. Each run consisted in the computation of the mean speed field evolution and computation of both instantaneous and dynamic travel time every 30 seconds, from 10 am to 6 pm. The runs took 320 CPU-hours, and were distributed on 8 servers equipped with 2.2 GHz dual core AMD Opteron CPUs and 8 GB RAM, which reduced the computation time to 39 hours. Despite taking nearly two weeks of CPU time, we note that our computational implementation is actually quite efficient, which enables computational experiments of this scale. Specifically, each simulation requires propagation of 100 ensembles through the discretized network flow model, and correcting the model predictions with the EnKF algorithm repeatedly over a period of eight hours, before the instantaneous and dynamic travel time estimates are computed. In total, each eight-hour numerical experiment takes 20 min to complete, which means it runs about 24 times faster than real-time.

4.2.2.
Using only inductive loop detector data in the model. The first analysis of the traffic estimates is based on the results obtained when using inductive loop detector data as the only input to the model. These results give us a baseline for the comparison between probe and loop detector data. A total number of 17 runs were conducted based only on the inductive loop detector data, by varying the number of sensors according to the selection algorithm in Section 3.2. Both instantaneous and dynamic travel times were computed. The labels of the inductive loop detector stations used in the estimation are presented in Table 1 as a function of the number of stations selected (see also Figure 1a).
The results of these runs are shown in Figure 5. The subfigures show the estimation error broken down by time of the day, as defined in Figure 2. During the morning accident (Figure 5a), the dynamic travel times converge to estimates with 7% error, while the instantaneous estimates remain above 20% error. The instantaneous and dynamic estimates have between 6% and 7% error during the free flow and congestion building periods (Figure 5b and 5c), and 13% error during the full congestion period (Figure 5d), with the instantaneous and dynamic estimates performing similarly.
The number of inductive loop detector stations used tends to have a positive impact on the quality of the estimate when less than eight inductive loop detector stations are used. Note that the curve is not monotonically decreasing. This is because when only a few sensors are deployed, the error becomes highly dependent on the placement of the sensors. It is expected that an optimal sensor placement  (Table 1). However, using data from more than eight inductive loop detector stations does not improve the quality of the estimates. If fewer than three inductive loop detectors are used, the estimation error is unacceptably high, at some points reaching as high as 100% error.

4.2.3.
Using only VTL data in the model. The second part of the analysis consists of the travel time estimates obtained when using VTL data only. The changing parameters of the input data are the number of VTLs deployed on the experiment site and the rate of the probe vehicles used to produce speed measurements at the locations of VTLs. The estimation errors of the travel times obtained with the dynamic method are shown in Figure 6.
In each of the time periods, estimates of the travel time can be achieved with less than 6% MAPE, with sufficient probe vehicles and virtual trip lines. However, when more than 137.5 veh/hr are used with more than 2.54 VTL/mi, only small improvements in the accuracy of the estimates can be achieved. When compared with inductive loop detectors, using 137.5 veh/hr and 2.54 VTL/mi performs as well as the estimates using more than eight inductive loop detector stations during the morning accident, free flow, and congestion building periods, but has less than half the error of inductive loops during the full congestion period.

4.2.4.
Fusing VTL and loop detector data. The dynamic travel time estimation error using both VTL and loop detector data simultaneously is assessed in Figure 7, where the change in the dynamic travel time MAPE due to the addition of data from six inductive loop detectors is computed. The results shown are a representative subset of all the runs performed when mixing the two data types.
At low probe data rates during the morning accident, free flow, and congestion building periods, adding inductive loop detector data increases the accuracy of the dynamic travel time estimates. For example, during the morning accident (Figure 7a), with a probe rate of 27.5 veh/hr and a VTL spacing of 0.79 VTL/mi, adding inductive loop detector data reduced the error from 29% to 8%. During the full congestion period, the dynamic travel time estimate accuracy decreased when inductive loop detector data was added at low probe rates (27.5 veh/hr). This is likely due to the fact that the estimates based on virtual trip line data only were unusually accurate, even performing better than simulations with more probe vehicles. At higher penetration rates (above 137.5 veh/hr) adding data from the six inductive loops has negligible effect, increasing or decreasing the accuracy only slightly. The exception is during the free flow period, when the MAPE increased (between 0.05 and 0.08) even at high probe rates, when 0.79 VTL/mi were used. The errors in the free flow period are magnified due to the small base travel time, which is under 10 minutes, and it is in fact not constant during the period (see Figure 2). Moreover, it is clear from Figure 1 that there is an area of heavy congestion around postmile 26 even during the free flow period, which is difficult to capture correctly with sparse sampling. can be determined at any time on any route using the speed estimates, and used as a proxy for dynamic travel times. As was shown for the inductive loop detector data in Figure 5b and 5d, instantaneous and dynamic travel time estimates are very similar when traffic conditions change sufficiently slowly. The same holds when estimating travel times from probe data. By looking at the instantaneous travel time errors in Figure 8a, an interesting result can be seen. The results suggest that adding more probe data results in an increased travel time estimation error. However, this result is expected, and can be explained by focusing on the scenarios (in Figure 8a) in which the penetration rate of the probe vehicles is low and no loop detectors are used. Here, the instantaneous travel time estimate performs well, and may seem like a valid estimate of the true travel time during the incident. However, this gives a misleading indication of the quality of these travel time estimates. The good performance of the instantaneous estimate is caused by the fact that the current state of the traffic (speed field) is very poorly captured in the underlying scenario and the speed of the traffic is heavily overestimated. This causes the instantaneous travel time estimate to perform as a good predictor of the future traffic conditions, namely, as a predictor of the clearing incident. When the number of probe measurements increases, the speed field estimate is captured more accurately, and the increased error in the travel time estimate is caused by the instantaneous approximation. 5. Summary. In this study, trade-offs between velocity data collected from GPS smartphones in probe vehicles, and velocity data obtained from inductive loop detectors, for the purpose of computing travel times on a stretch of highway were studied.
This work used experimental probe data obtained from the Mobile Century field experiment and loop detector data from PeMS. The measurements were combined with a mathematical traffic model in a highway traffic estimation algorithm using a data assimilation technique called ensemble Kalman filtering, developed as a part of the Mobile Millennium project. The results of the algorithm were compared against the true travel times experienced by the drivers, obtained through license plate re-identification. A number of scenarios were created in which the volume of the probe data and number of inductive loop detector stations available for the estimation algorithm could be adjusted.
There are several practical results which were obtained from the extensive experimental studies including: • Achieving 10% error for dynamic travel times. In this study, it was found that the dynamic travel time estimates can be achieved with less than 10% error when using a flow model with data assimilation, by using either inductive loop detector data, probe data, or a mixture of both inductive loop detector data and probe data. Moreover, the estimates from virtual trip line-based probe data can achieve a higher degree of accuracy when all available probe data is used compared to estimates from inductive loop detectors when all inductive loops on the experiment site are used, although in general the performance is similar. • Minimum loop detector spacing for travel time estimation. In this study, using data from more than eight inductive loop detector stations (average spacing 0.83 miles) did not give additional benefit in the travel time estimation. The error remains constant between 6-13% depending on the time of day, regardless of the added loop detector stations. • Diminishing travel time accuracy improvement. When sampling probe vehicles at a rate of 137.5 veh/hr with more than 2.54 VTL/mi, increasing the number of probe measurements by adding more probe vehicles or additional trip lines causes only a small improvement on the travel time accuracy. • A mixture of probe and loop detector data in travel time estimation. It was found that when complementing loop detector data with probe vehicle data, better estimates for travel times are obtained, especially at low penetration rates. For example, if using loop detectors spaced more than 2.11 miles apart, probe data can give over 50% increase in the travel time accuracy. The above findings may assist practitioners to identify the value of flow model based traffic monitoring systems for computing travel times.


extfieldsofstudy: Engineering Medicine

text:
A novel channel selection method for multiple motion classification using high-density electromyography

Background Selecting an appropriate number of surface electromyography (EMG) channels with desired classification performance and determining the optimal placement of EMG electrodes would be necessary and important in practical myoelectric control. In previous studies, several methods such as sequential forward selection (SFS) and Fisher-Markov selector (FMS) have been used to select the appropriate number of EMG channels for a control system. These exiting methods are dependent on either EMG features and/or classification algorithms, which means that when using different channel features or classification algorithm, the selected channels would be changed. In this study, a new method named multi-class common spatial pattern (MCCSP) was proposed for EMG selection in EMG pattern-recognition-based movement classification. Since MCCSP is independent on specific EMG features and classification algorithms, it would be more convenient for channel selection in developing an EMG control system than the exiting methods. Methods The performance of the proposed MCCSP method in selecting some optimal EMG channels (designated as a subset) was assessed with high-density EMG recordings from twelve mildly-impaired traumatic brain injury (TBI) patients. With the MCCSP method, a subset of EMG channels was selected and then used for motion classification with pattern recognition technique. In order to justify the performance of the MCCSP method against different electrode configurations, features and classification algorithms, two electrode configurations (unipolar and bipolar) as well as two EMG feature sets and two types of pattern recognition classifiers were considered in the study, respectively. And the performance of the proposed MCCSP method was compared with that of two exiting channel selection methods (SFS and FMS) in EMG control system. Results The results showed that in comparison with the previously used SFS and FMS methods, the newly proposed MCCSP method had better motion classification performance. Moreover, a fixed combination of the selected EMG channels was obtained when using MCCSP. Conclusions The proposed MCCSP method would be a practicable means in channel selection and would facilitate the design of practical myoelectric control systems in the active rehabilitation of mildly-impaired TBI patients and in other rehabilitation applications such as the multifunctional myoelectric prostheses for limb amputees.

Introduction
Surface electromyography (EMG) is an electrical manifestation of muscle contractions, which has been used as a convenient and proprioceptive control signals for the operations of limb prostheses [1], assistive exoskeletons [2], and rehabilitation robots [3]. In recent years, one of the most popular applications of EMG signals would be to use EMG-based pattern recognition (EMG-PR) algorithms to control multifunctional powered prostheses [1,[4][5][6][7]. More recently, the EMG-PR also was adopted in the neural rehabilitation of the motor impaired patients who suffer from neurological injuries for development of an active rehabilitation system. The previous studies have proved that the active rehabilitation is promising in enhancing the therapeutic effect [8] and accelerating the brain plasticity [9]. In light of these findings, the EMG-PR algorithm was proposed to identify the motor intentions of stroke survivors and patients with incomplete cervical spinal cord injury [10,11]. The results of these previous studies have suggested that EMG-PR algorithms would be feasible and useful in active rehabilitation of patients suffering from neurological injuries.
High density EMG recordings have been increasingly applied in EMG-PR algorithms [12][13][14] with an attempt to capture more temporal and spatial information about the muscle activities and electrophysiology. With more channels of temporal information that can reflect the activities of some small forearm muscles, a better motion classification performance could be achieved when using high density EMG recordings [13,14]. And the topographical maps of EMG amplitude can be used to examine the exact electrode locations where the strong myoelectric activity is experienced during a motion task [12][13][14]. However, the data processing of high-density EMG recordings from a large number of EMG electrodes is computationally expensive, making it impractical in real-time myoelectric control. Additionally, the high dimensional EMG recordings may cause the classifier to over-fit the training data due to the irrelevant or redundant information. Therefore, selecting an appropriate number of EMG channels with desired classification performance and determining the optimal placement locations of electrodes would be necessary in the practical myoelectric control for a potential user.
Currently, there are two commonly used ways to determine the number and the locations of EMG electrodes in the applications of EMG-based motion classification. One direct way is based on the physiologically known anatomical knowledge of skeletal muscles [15][16][17], in which the electrode sites would be chosen by viewing the intensity, repeatability and consistency of multi-channel EMG signal recordings. For the persons with the intact functions of skeletal muscles, using the muscle-physiology-based way we may quickly and easily determine the appropriate sites of EMG electrodes since their almost consistent muscle anatomy could ensure producing very similar muscle contractions when they actuate a limb activity. However, for the persons with some muscular issues such as the patients with a post-stroke, traumatic brain injury or limb amputation, when doing a limb activity, their muscle contractions would be different from healthy people and may be various between patients since the extents, the positions and the causes of their muscular damages might be inconsistent. Thus it would be hardly to use this kind of clinical method for determining the appropriate electrode sites in patients.
Another way to determine a subset of appropriate EMG channels is based on a certain optimizing criterion. A commonly used method is the sequential forward selection (SFS)-based channel selection algorithm [14,16], in which the most informative EMG channels were determined by an optimization procedure with an optimizing objective of high motion classification accuracy. Since SFS method requires to repeatedly perform the searching procedure for the optimal channels until obtaining a given number of channels, this classification-accuracy-based method would take a long computation time to get a relatively optimal subset of EMG channels. Moreover, SFS method must rely on a specific set of EMG features and a specific EMG pattern-recognition algorithm to do the selection of the optimal channels for an EMG control application. When the EMG features and/or pattern-recognition algorithm would be changed for some reasons such as trying to use different features or algorithms for higher classification performance, the selected EMG channels might be not optimal yet. Thus the optimization procedure would be executed again to search the optimal channels accordingly. Some feature-dependent algorithms such as the quasi-optimal channel selection method based on partial Kullback-Leibler information [18] and coefficient-based channel selection method using principal component analysis [19] also have been proposed and used for channel selection. In addition, another two feature selection methods, Fisher-Markov selector (FMS) based on Fisher criterion [20] and the minimal-redundancy-maximal-relevance based on mutual information [21], also fall into the group of feature-dependent channel selection methods. It is obvious that these feature-dependent channel selection methods would rely on the EMG features. Similarly, when changing the set of EMG features such as increasing/decreasing a feature to/from it, the selected channels might not work well again, thus another channel selection procedure would be required to re-determine the appropriate channels accordingly.
Alternatively, the direct channel selection method (i.e., variable selection) which directly works on the raw EMG data also has been proposed to determine a set of appropriate channels for EMG recordings. The Monte Carlo method reported by Nagata et al. for hand motion classification is one of such kinds of the direct channel selection approaches [22]. Unlike the classification-accuracy-based and the feature-dependent channel selection methods, this raw-data-based method is independent of EMG features and classification algorithms. Thus a unique set of optimal EMG channels might be determined from the multi-channel EMG recordings. As a result, when the features and/or classification algorithm are changed, the selected EMG channels would be retained. In this study, we proposed a novel direct channel selection approach named multi-class common spatial pattern (MCCSP) for channel selection in EMG patternrecognition-based movement classification. Since MCCSP would be independent on EMG features and classification algorithms, it would be more convenient for channel selection in developing an EMG control system than the exiting methods. The performance of the proposed MCCSP method in selecting an optimal set of EMG channels was assessed with high-density EMG recordings from twelve mildly-impaired traumatic brain injured (TBI) patients. Generally, the TBI patients would suffer from the physical and behavioral disorders, so appropriately physical therapy should be employed for the rehabilitation of their physical functions. The EMG-control active rehabilitation systems would be a promising way for TBI patients' neural function recovery [23]. In addition, the performance of the proposed MCCSP method was compared with that of two exiting channel selection methods, SFS and FMS.

Data collection
Twelve mildly-impaired male patients with TBI participated in this study. They were chosen based on the upper limb motion impairment level assessed by a physical therapist. According to the definition of the stages in the Brunnstrom Assessment Scale, they were in stage IV-V and got the scores of 49 to 61 with the Fugl-Meyer Assessment of Sensorimotor Recovery after stroke, in which a zero score denotes no any function and a score of 66 designates as normal function. All of the subjects did not have any experience of attending this kind of research study before. In the experiment, they were asked to use their unilateral arm with severer motor impairment to perform 21 forearm and hand movements ( Figure 1) plus one "no movement". Each movement was maintained for 6 s with a moderate force and repeated 6 times. A rest time of 8 s was set between two successive movements in each trial. All subjects could choose to finish all or part of the 22 movements based on their own motor ability. The Research Ethics Board of the Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences, approved the experimental protocol of this study, and each subject gave written informed consent and provided permission for publication of photographs for a scientific and educational purpose.
The high-density EMG acquisition system (Refa-128, TMS International BV, Netherlands) was used to record the EMG signals during the experiment. The 56 monopolar electrodes (5 mm in diameter) were placed on the forearm and hand of subjects, as shown in Figure 2. The 48 of 56 electrodes were placed on the forearm in an 8 × 6 grid from 1 cm proximal to the elbow crease to 1/3 distal to the wrist joint with an electrode inter-distance of around 2 cm and other eight electrodes were placed on the hand muscles with two electrodes on the first dorsal interosseous, three on the thenar group muscles, and three on the hypothenar group muscles. A reference electrode was fixed on a nylon bracelet that was worn on subject's wrist. The sampling rate of EMG signals was set as 1024 Hz.

EMG data preprocessing
The EMG signals recorded with 56 electrodes were digitally filtered with a five-order Butterworth high-pass filter at 30 Hz, and then the 50 Hz power line interference was reduced with a notch filter from EMG recordings. Two electrode configurations, monopolar mode and bipolar mode, were considered in the study with an attempt to see if the two modes of EMG signals would provide different performance in classifying the movements. The raw recordings were 56 monopolar EMG signals (Figure 3(a)). The bipolar EMG recordings were formed with the differential values between two adjacent monopolar channels along the orientation of muscle fibers from the 56 monopolar EMG signals, resulting in 45-channel bipolar EMG signals, as shown in Figure 3(b).

MCCSP algorithm
MCCSP-based channel selection algorithm considers the EMG data from all the channels simultaneously. To clearly explain its principle and implementation, the principle of two-class common spatial pattern (CSP) algorithm, the CSP-based channel selection algorithm, and the MCCSP-based channel selection algorithm were successively introduced as follows.

CSP algorithm
CSP is a supervised two-class filter algorithm by calculating the linear spatial filters that maximize the variance of one class and meanwhile minimize the variance of another class [24,25]. In this way, the two classes can be maximally separated by their variances. Suppose we want to classify two motion tasks, class A and class B, and use X A and X B to denote their corresponding signal matrix with a dimension n × d, where d is the number of channels and n is the number of samples per channel. The object becomes to find a spatial filter matrix W which can maximize the variance of class A and minimize the variance of class B, which would correspond to an optimization problem that can be formulated as: By solving a generalized eigenvalue problem, the linear spatial filter matrix W can be obtained by simultaneously diagonalizing the covariance matrix Σ A and Σ B : With the constraint condition D A + D B = I, the eigenvectors with the largest eigenvalues for D A have the smallest eigenvalues for D B and vice versa. Applying the filter matrix W to the raw signals X would give d output signals Y = W * X that were also called components. The variance of each component is indicated by its corresponding eigenvalue on the principal diagonal of D A , for class B of D B . In the binary classification problem of EEG signals, the first component that has the largest variance for class A and the smallest variance for class B and the last component that has the smallest variance for class A and the largest variance for class B were commonly used [25].  4 5 6 7 8   9 10 11 12 13 14 15 16   17 18 19 20 21 22 23 24   25 26 27 28 29 30 4 5 6 7 8   9 10 11 12 13 14 15 16   17 18 19 20 21 22 23 24   25 26 27 28 29 30

CSP-based channel selection
With W −1 (inverse matrix of W) that is called spatial patterns, the original signal matrix X can be reconstructed by X = W −1 * Y. The first and last columns of W −1 would be the most important spatial patterns that explain the largest variance of one motion and the smallest variance of the other. If considering the coefficients of spatial patterns as projection weights, the two channels corresponding to the maximal coefficients of spatial pattern vectors may be the channels that are most correlated with the motion specific sources. This method has been used by Wang et al. to select four EEG channels from 118 channels for classifying three motor imaginary tasks and achieved a high classification accuracy [24].

MCCSP-based channel selection
Given a problem with multiple classes of movements in the study, we firstly extended the two-class CSP into multi-class form by using one-versus-rest (OvR) scheme, as done in [26]. For each motion class, a corresponding filter matrix W C −1 was calculated by maximizing the variance of the class and minimizing the sum of the variances of all other classes: Accordingly, multiple spatial pattern matrices corresponding to all the included motion classes would be obtained, where c denotes the class label. Then we used these spatial pattern matrices to select the optimal channels for subsequent pattern recognition analysis. In the study, we adopted similar coefficients-based channel selection method as did in the two-class CSP. For each spatial pattern matrix corresponding to a specific motion class, the two channels corresponding to the maximal coefficients of the first and the last spatial pattern vectors were selected, and then all the selected two-channel from each spatial pattern matrix W C −1 were combined together, which were then used as the optimal channels in the subsequent pattern recognition analysis. Note that for the two channels selected from a specific spatial pattern matrix, e.g., W A −1 , one or both of which may be overlapped with that obtained from another spatial pattern matrix W B −1 . So the selected channels for all movements were ranked in the order of their occurrence frequency from high to low. The number of selected EMG channels might vary for different subjects.

SFS algorithm
SFS is an iterative searching procedure, in which one optimal channel that produces the highest classification accuracy was firstly selected among all the channels, and then another channel that can achieve the maximum classification accuracy with the combination with the selected channels was added [14,16]. Suppose we require to choose an optimal subset with n channels from all the 56 channels, the classification computation procedure would be repeatedly implemented by (56 + 56 + 1 − n) * n/2 times.

FMS algorithm
FMS is essentially a feature selection algorithm, which was proposed to select the globally optimal subset of features from high-dimensional feature space in the spirit of Fisher's class separation criterion [20]. By using specific kernel functions and the Markov random field optimization techniques, the FMS is capable of selecting features and has efficient computational complexity. In this study, we firstly applied FMS to the EMG feature matrix to get the indices of all the features that were ranked in the order of their importance coefficients, and then the feature indexes were mapped to their corresponding channel indexes. In this way, the optimal EMG channels could be selected.

Feature extraction and motion classification
With a selected subset of EMG signals, a shifting analysis window with a time length of 150 ms and an increment of 100 ms (50 ms overlapping) was used to segment the EMG signals into a series of analysis windows. Then two commonly used time-domain feature sets were extracted from each analysis window, respectively. They were (1) four time-domain features (TD) [6], mean absolute value (MAV), number of zeros crossings (ZC), number of slope sign changes (SSC), and waveform length (WL), and (2) the six order autoregressive (AR) model coefficients [27] plus the root mean square (RMS) amplitude of EMG signals (TDAR). For each analysis window, the features extracted from all the selected EMG channels was concatenated to form a feature matrix, which was then fed into a classifier for motion classification.
Two classifier algorithms, the linear discriminant analysis (LDA) and the k-nearest neighbor (KNN), were used in this study. The LDA classifier was used for its merit of algorithm simplicity and low computational cost [5], and the KNN classifier was used for its tolerance to arbitrary data distribution [28]. In order to find the best featureclassifier combination for each selected EMG subset, each of the two classifiers was separately combined with each of the two feature sets, yielding four different featureclassifier combinations.

Performance evaluation
The classification accuracy, defined as the percentage of the number of correct classification decisions over the total number of classification decisions, was used to assess the motion classification performance [7]. Five-fold cross validation was performed in motion classification when using the EMG subset selected via MCCSP and FMS for the purpose of evaluating their performance in EMG-PR. The EMG feature matrix extracted from the selected EMG recordings was randomly divided into five subsets with equal length. Four of the five subsets were used as a training set and the remaining one subset was considered as a testing set in each classification process to compute the classification accuracy. The classification accuracies were averaged over the five validation results from the folds. For the EMG subset selected by SFS method, the first half of the feature matrix was used as the training set and its second half was used as the testing set. No cross validation was conducted herein considering the newly added EMG channel was determined according to the maximum classification accuracy during last round repetition.

Statistical analysis
To assess the statistical difference among the three channel selection methods, oneway ANOVA was performed in terms of motion classification accuracy with SPSS Statistical Modeling Software (SPSS 17.0 IBM Corp., Chicago, IL). In addition, we also analyzed the effect of feature set, classifier, and combination of feature set and classifier on the classification performance by performing two-way ANOVA. The interaction effect of feature set and classifier on classification performance was firstly assessed. If the effect was insignificant, the main effects of features set and classifiers were examined. The level of statistical significance was set to p < 0.05.

Results
Number of selected EMG channels using MCCSP Figure 4 illustrates the numbers of selected EMG channels by MCCSP for all the 12 subjects when using EMG signals from monopolar and bipolar channel configurations, respectively. It can clearly observe from Figure 4 that the number of selected EMG channels varied for different subjects, ranging from 18 to 30 for the two channel configurations. And we can see that the numbers of selected channels by using monololar and bipolar EMG signals, respectively, were similar (in 10 subjects) or same (in two subjects). In order to make the following results comparable among all the subjects, the maximum number of selected EMG channels was set as 18. Note that for a subject the number of selected channels were corresponding to the movement classes that could be completed by him/her. Since the motor impairment of all the subjects was different, the number of completed movement classes (nClass designated by the yellow squares in Figure 4) would vary with a range of nClass = 19-22.

Optimal feature-classifier combination
To find the optimal feature-classifier combination for each selected EMG subset, the average classification accuracy across all subjects were calculated and shown in Figure 5 when using TD-LDA, TDAR-LDA, TD-KNN, and TDAR-KNN, respectively. The x-axis denoted the number of involved EMG channels that were selected by using MCCSP (Figure 5 EMG subset selected by MCCSP, TDAR-KNN and TD-KNN had almost the same classification performance. In addition, when using the EMG subset determined by SFS, the difference among all the four feature-classifier combinations became smaller. We further investigated whether the feature set, classifier, or the combination of feature and classifier had a significant impact on the classification performance when using the EMG subset determined by MCCSP (Figure 6(a) and (b)), SFS (Figure 6(c) and (d)), and FMS (Figure 6(e) and (f )), respectively. The results demonstrate that when using the monopolar EMG channels determined by MCCSP (Figure 6(a)), the monopolar EMG channels determined by SFS (Figure 6(c)), and the bipolar EMG channels determined by SFS (Figure 6(d)), the classification performance was relatively stable with respect to different feature sets and classifiers. However, when using the bipolar EMG channels selected by MCCSP (Figure 6(b)), the monopolar EMG channels selected by FMS (Figure 6(e)), and the bipolar EMG channels selected by FMS (Figure 6(f )), the classification performance became sensitive to the choice of feature set and/or classifier.  Figure 7 shows the comparison of average classification accuracy across all subjects when using the EMG subsets (1-18 channels) selected by MCCSP, SFS, and FMS and their corresponding optimal feature-classifier combinations. The results indicate that applying TD-KNN to the EMG subset determined by MCCSP brought the highest classification accuracy and the best convergence when 4 to 13 optimal monopolar electrodes (Figure 7(a)) and 2 to 18 optimal bipolar electrodes (Figure 7(b)) were utilized. But the difference was not significant (p-value > 0.05). Figure 8 compares the two electrode configurations in terms of motion classification accuracy, which was obtained by using the EMG subset selected via MCCSP and the TD-KNN feature-classifier combination. The results show that the bipolar configuration was consistently better than the monopolar configuration with the increase of the number of included EMG channels. In addition, 18 EMG channels might be sufficient  to get a proper motion classification performance. Because the average classification accuracy was 93.03% when using 18 optimal monopolar EMG channels, only 1.47% lower than that when using all 56 monopolar EMG channels (94.50%), and 18 optimal bipolar EMG channels brought an average classification accuracy of 95.58%, only 2.59% lower than that when using all 45 bipolar EMG channels (98.17%).

Comparison of distribution of selected EMG channels
We also analyzed the distribution of selected EMG channels when using different feature-classifier combinations for each subject. Figure 9 demonstrated a representative example from a subject (BI01), where 20 monopolar EMG channels (Figure 9(a)) and 10 bipolar EMG channels (Figure 9(b)) were selected, respectively. It was found that for both electrode configurations, the distribution of EMG channels selected by SFS varied with the choice of feature-classifier combination. When using FMS, the distribution of selected EMG channels varied slightly following the change of feature set. When using MCCSP, however, a fixed channel distribution was obtained independent of feature and classifier. The channel distributions for other subjects almost presented the similar patterns as the subject.

Discussion
Generally speaking, using more EMG electrodes could capture more electrophysiological information that may improve the performance of movement classification in EMG-PR-based control system such as multifunctional myoelectric prostheses and active rehabilitation robots. However, this would simultaneously increase the complexity and cost of the EMG controlled systems. Thus it becomes necessary and important to find an appropriate number of electrodes and their locations for the high performance of a myoelectric control system before it is clinically viable. Theoretically, to get the optimal channels from high density channels, all the combinations of desired number of channels should be involved in the pattern recognition analysis. But it often be impractical to search a global optimal solution when a large number of channels are involved. For example, searching 6 from 56 channels would result in 32,468,436 different combinations. So some suboptimal methods have been proposed and used in the previous studies [14,16,18,19]. Two commonly used methods for channel selection are SFS and FMS, which are dependent on either classification accuracy or EMG features. While these methods might be applicable for channel selection, they have some limitations, as mentioned above in Introduction section. In this study, a novel direct channel selection algorithm was proposed and its performance in selecting the appropriate channels has been evaluated by a high-density EMG recordings from twelve mildly-impaired TBI patients. Since the newly proposed method (MCCSP) for channel selection is independent on EMG features and classification algorithms, it should be more convenient for channel selection in comparison to the SFS and FMS, especially when changing the EMG features or the classification algorithms is needed.
By examining two different types of EMG features and two pattern-recognition algorithms, we found that the TD features combining with a KNN classifier would be a better configuration for the selected EMG subset via MCCSP, and TDAR features with a LDA classifier would be a better one for the determined channels by using SFS (Figure 5). In addition, the electrode configuration is also a factor that would affect the motion classification performance. When the number of selected electrodes was increasing, the motion classification accuracy in bipolar electrode configuration reached a plateau (Figures 5 and 7) more quickly than that in monopolar electrode configuration. Furthermore, using same number of EMG channels, the bipolar electrode configuration always achieved higher classification accuracy than the monopolar electrode configuration (Figure 8). These findings are similar as the conclusions obtained by Zhou and Huang et al. [14,16]. Generally, a monopolar electrode can provide more electrophysiological information than a bipolar electrode. However, different EMG applications may need different information involved in EMG signals.
For the EMG-based motion classification, the useful information should be these components of EMG signals that are characteristic of discrepancy for different movements, which could make a classifier to easily identify a movement. This may be one reason why the bipolar electrodes would outperform the monopolar electrodes in movement classification. Another possible reason underlying these findings would be that compared to a monopolar electrode, a bipolar electrode acts as a differential operator that can remove some common components such as direct-current part and common-mode noise in EMG signals. These common components in different channels might not provide any discriminate information for movement classification, even making the classification performance worse. In future study, we are interested in exploring the reasons. As a commonly used channel selection method in EMG-PR based classification, the previous studies have shown that SFS would be an effective method in selecting the appropriate EMG channels [7,13,14,16]. However, in comparison to the proposed MCCSP, SFS was consistently worse just with one exception of using only one EMG channel (Figure 7). With SFS method, the first chosen EMG channel would be the global optimum. However, when more channels were considered, it was impossible to get the globally optimal EMG subset by using the SFS method. Compared with FMS, the proposed MCCSP performed slightly better (Figure 7) in classifying different arm/hand movements. In addition, the computational cost should be also an interesting performance metric in channel selection. In this study, we found that the average cost time taken by searching 18 optimal monopolar EMG channels was about 0.98 s when using the proposed MCCSP, significantly lower than that when using SFS (with TD feature and LDA classifier) (about 2096 s), FMS (with TD feature) (10.56 s), or FMS (with TDAR feature) (18.67 s), respectively. Additionally, a fixed combination of the selected EMG channels could be achieved when using MCCSP (Figure 9). These outcomes of this study may suggest that as an alternative, the MCCSP might be an effective and practicable choice for channel selection in the design of a practical myoelectric control system. Note that the channels selected by the three different methods in the study were inconsistent, but the similar classification performance was achieved by these selected channels. This would be because the three methods depend on different principles in channel selection. As mentioned above, as a direct channel selection algorithm the proposed MCCSP uses the raw EMG data to select the appropriate channels, the FMS and the SFS rely on EMG features and classification accuracy for searching the optimal channels, respectively. Thus the different channel combinations would be obtained by using the three methods. However, this does not mean there is no need to carefully select the channels. In the EMG-PR movement classification, it is the patterns of EMG signals from the selected channels that are considered as a whole to train and test a classifier for movement classification. Although different channel combinations were determined by the three methods, their EMG patterns all could provide identifiable information to properly classify different movements. So a similar classification performance could be achieved by the three methods. This may indicate that the channel combination for similar classification performance would be not only one.
It should be noted that all the results of this study were conducted by means of postprocessing and evaluated with offline classification accuracy which might not reflect the real-time performance directly [29]. In real-time applications, the EMG signals would be easily affected by some factors like muscle contraction force [30], skin impedance [31], electrodes shift [30,32], arm position variation [33,34], and muscle fatigue [30,35]. The stability of the proposed channel selection algorithm was not considered in current study. Note that Huang et al. investigated the temporal stability of the SFS-selected electrodes by validating the EMG data recorded from follow-up experiment in terms of classification performance [16]. The repeatability of EMG signals recorded from the selected electrode sites and their inter-parameter agreement [15] would be also necessary to examine the stability of channel selection algorithm. With these important issues, we will examine the stability of the channels selected by the proposed MCCSP in our future work.


extfieldsofstudy:
Computer Science
Mathematics

text:
Stochastic Gradient Estimate Variance in Contrastive Divergence and Persistent Contrastive Divergence

Contrastive Divergence (CD) and Persistent Contrastive Divergence (PCD) are popular methods for training the weights of Restricted Boltzmann Machines. However, both methods use an approximate method for sampling from the model distribution. As a side effect, these approximations yield significantly different biases and variances for stochastic gradient estimates of individual data points. It is well known that CD yields a biased gradient estimate. In this paper we however show empirically that CD has a lower stochastic gradient estimate variance than exact sampling, while the mean of subsequent PCD estimates has a higher variance than exact sampling. The results give one explanation to the finding that CD can be used with smaller minibatches or higher learning rates than PCD.

One known feature of Contrastive Divergence (CD) learning is that it yields a biased estimate of the gradient (Bengio & Delalleau, 2009;Carreira-Perpinan & Hinton, 2005). On the other hand, it is known to be fast for reaching good results (Carreira-Perpinan & Hinton, 2005;Tieleman, 2008). In addition to the computationally light sampling procedure in CD, it is claimed to benefit from a low variance of the gradient estimates (Hinton, 2002;Carreira-Perpinan & Hinton, 2005). However, the current authors are not aware of any rigorous research on whether this claim holds true, and what the magnitude of the effect is 2 .
On the other hand, Persistent Contrastive Divergence (PCD) has empirically been shown to require a lower learning rate and longer training than CD 3 (Tieleman, 2008). The authors propose that the low learning rate is required since the model weights are updated while the Markov chain runs, which means that in order to sample from a distribution close to the stationary distribution the weight cannot change too rapidly. However, for similar reasons that CD updates are assumed to have low variance, subsequent PCD updates are likely to be correlated leading to a possibly undesirable "momentum" in the updates. This behavior would effectively increase the variance of the mean of subsequent updates, requiring either larger minibatches or smaller learning rates.
In this paper we explore the variances of CD, PCD and exact stochastic gradient estimates. By doing so, we hope to shed light on the observed fast speed of CD learning, and on the required low learning rate for PCD learning compared to CD learning. Thereby we hope to contribute to the understanding of the difference between CD and PCD beyond the already well documented bias of CD.

Contrastive Divergence and Persistent Contrastive Divergence
A restricted Boltzmann machine (RBM) is a Boltzmann machine where each visible neuron x i is connected to all hidden neurons h j and each hidden neuron to all visible neurons, but there are no edges between the same type of neurons. An RBM defines an energy of each state (x, h) by −E(x, h | θ) = b x + c h + x Wh, and assigns the following probability to the state via the Boltzmann distribution: W} is a set of parameters and Z(θ) normalizes the probabilities to sum up to one.
The log likelihood of one training data point is hence Sampling the positive phase of the gradient of the log likelihood ∂φ + ∂W is easy, but sampling the negative phase ∂φ − ∂W is intractable. A popular method to solve sampling of the negative phase is Contrastive Divergence (CD). In CD, the negative particle is sampled only approximately by running a Markov Chain a limited number of steps (often only one step) from the positive particle (Hinton, 2002). Another method, called Persistent Contrastive Divergence (PCD) solves the sampling with a related method, only that the negative particle is not sampled from the positive particle, but rather from the negative particle from the last data point (Tieleman, 2008).

Experiments
In order to examine the variance of CD and PCD gradient estimates, we use an empirical approach. We train an RBM and evaluate the variance of gradient estimates from different sampling strategies at different stages of the training process. The sampling strategies are CD-k with k ranging from 1 to 10, PCD, and CD-1000 that is assumed to correspond to an almost unbiased stochastic gradient. In addition, we test CD-k with independent samples (I-CD), where the negative particle is sampled from a random training example. The variance of I-CD separates the effect of the negative particle being close to the data distribution in general, and the effect of the negative particle being close to the positive particle in question.
We use three different data sets. The first is a reduced size MNIST (LeCun et al., 1998) set with 14 × 14 pixel images of the first 1 000 training set data points of each digit, totaling 10 000 data points. The second data set are the center 14 × 14 pixels of the first 10 000 CIFAR (Krizhevsky & Hinton, 2009) images converted into gray scale. The third are the Caltech 101 Silhouettes , with 8 641 16 × 16 pixel black and white images. We binarize the grayscale images by sampling the visible units with activation probabilities equal to the pixel intensity.
We set the number of hidden neurons equal to the number of visible neurons. The biases are initialized to zero, while the weights are initially sampled from a zero-mean normal distribution with standard deviation 1/ √ n v + n h where n v and n h are the number of visible and hidden neurons, respectively. We train the model with CD-1, and evaluate the variance of the gradient estimates after 10, and 500 epochs. We use Adaptive learning rate (Cho et al., 2011) with an initial learning rate of 0.01. We do not use weight decay.
In all of the gradient estimates, the final sampling step for the probabilities of the hidden unit activations is omitted. The gradient estimate is therefore based on sampled binary visible unit activations, but continuous hidden unit activation probabilities conditional on these visible unit activations. This process is called Rao-Blackwellisation , and is often used in practice. The variance is calculated on individual gradient estimates based on only one positive and negative particle each. In practice, the gradient is usually estimated by averaging over a mini-batch of N independent samples, which diminishes the variance Nfold. We ignore the bias gradient estimates.
When analyzing subsequent PCD gradient estimates, the negative particles of the first estimate are sampled 1 000 steps from a random training example. Subsequent k estimates are then averaged, where the positive particle is randomly sampled from the data for each step while the negative particle is sampled from the previous negative particle. No learning occurs between the subsequent estimates. We can therefore disentangle the effects of weight updates from the effect of correlation between subsequent estimates.
We iterate all results for 10 different random initializations of the weights, and evaluate the variance by sampling gradient estimates of individual training examples 10 times for each training example in the data set. The variance is calculated for each weight matrix element separately, and the variances of the individual weights are then averaged.

Results
As we can see from Figure 1, the variance of Contrastive Divergence is indeed smaller than for exact sampling of the negative particle. We also see that the variance of CD estimates quickly increases with the number of CD steps. However, this effect is significant only in later stages of training. This phenomenon is expected, as the model is expected not to mix as well in later stages of training as when the weights are close to the small initial random weights.
If we sample the negative particle from a different training example than the positive particle (I-CD), in Figure 3 we see that the variance is similar or even larger compared to the variance with exact sampling. Although it is trivial that the variance of the I-CD estimates is higher than for CD, the interesting result is that I-CD loses all of the variance advantage against exact sampling. The result supports the hypothesis that the low variance of CD precisely stems from the fact that the negative particle is sampled from the positive particle, and not from that the negative particle is sampled only a limited number of steps from a random training example.  For subsequent PCD updates, we see in Figure 2 that the variance indeed is considerably higher than for independent sampling. Again, as expected this effect is stronger the later during training the evaluation is done.
When looking at the magnitude of the variance difference, we see that for CD-1, the mean of 10 subsequent updates have a multiple times smaller variance than PCD. In effect, this means that ignoring any other effects and the effect of weight updates, PCD would need considerably smaller learning rates or larger minibatches to reach the same variance per minibatch. This magnitude is substantial, and might explain the empirical finding that PCD performs best with smaller learning rates than CD.

Conclusions
Contrastive Divergence or Persistent Contrastive Divergence are often used for training the weights of Restricted Boltzmann machines. Contrastive Divergence is claimed to benefit from low variance of the gradient estimates when using stochastic gradients. Persistent Contrastive Divergence could on the other hand suffer from high correlation between subsequent gradient estimates due to poor mixing of the Markov chain estimating the model distribution.
In this paper, we have empirically confirmed both of these findings. In experiments on three data sets, we find that the variance of CD-1 gradient estimates are considerably lower than when independently sampling with many steps from the model distribution. Conversely, the variance of the mean of subsequent gradient estimates using PCD is significantly higher than with independent sampling. This effect is mainly observable towards the end of training. In effect, this indicates that from a variance perspective, PCD would require considerably lower learning rates or larger minibatches than CD. As CD is known to be a biased estimator, it therefore seems that the choice between CD and PCD is a trade-off between bias and variance.
Although the results in this paper are practically significant, the approach in this paper is purely empirical. Further theoretical analysis of the variance of PCD and CD gradient estimates would therefore be warranted to confirm these findings. In addition, we intend to repeat the experiments with a larger data set, and train the models and replace the CD-1000 baseline with better approximations of exact sampling using e.g. Enhanced gradient (Cho et al., 2011) and parallel tempering (Cho et al., 2010).