# A Probabilistic Approach for Discovering Daily Human Mobility Patterns with
Mobile Data

Weizhu Qian
CIAD, Univ. Bourgogne Franche-Comté,
UTBM, F-90010
Belfort, France
Email<EMAIL_ADDRESS>
Fabrice Lauri
CIAD, Univ. Bourgogne Franche-Comté,
UTBM, F-90010
Belfort, France
Email<EMAIL_ADDRESS>
Franck Gechter
CIAD, Univ. Bourgogne Franche-Comté,
UTBM, F-90010
Belfort, France
Email<EMAIL_ADDRESS>
###### Abstract

Discovering human mobility patterns with geo-location data collected from
smartphone users has been a hot research topic in recent years. In this paper,
we attempt to discover daily mobile patterns based on GPS data. We view this
problem from a probabilistic perspective in order to explore more information
from the original GPS data compared to other conventional methods. A non-
parameter Bayesian modeling method, Infinite Gaussian Mixture Model, is used
to estimate the probability density for the daily mobility. Then, we use
Kullback-Leibler divergence as the metrics to measure the similarity of
different probability distributions. And combining Infinite Gaussian Mixture
Model and Kullback-Leibler divergence, we derived an automatic clustering
algorithm to discover mobility patterns for each individual user without
setting the number of clusters in advance. In the experiments, the
effectiveness of our method is validated on the real user data collected from
different users. The results show that the IGMM-based algorithm outperforms
the GMM-based algorithm. We also test our methods on the dataset with
different lengths to discover the minimum data length for discovering mobility
patterns.

###### Index Terms:

Probabilistic model; Infinite Gaussian Mixture Model; Kullback-Leibler
divergence; Human mobility.

##  I Introduction

Smartphone devices are equipped with multiple sensors that can record user
behavior on the handsets. With the help of a large-scale smartphone usage
data, researchers are able to study human behavior in the real world. Since
location information is one of the crucial aspects of human behaving,
investigating human mobility from mining mobile data has become recently a
popular research topic.

Previous research in this filed mainly only focus on discovering the
significant places or predicting the transition among the significant places
[1], [2], [3]. However, these research neglect the data sampled at the places
where one stay for a relatively short time, for instance, in the middle of
transitions. As opposed to this point of view, we believe that these data is
important for revealing human mobility patterns as well.

To characterize human mobility, we should realize that there are multiple
patterns lying in the mobility data even for the same individuals. For
example, normally, on workdays, one goes to work or school at daytime while on
weekends, he/she may prefer to stay at home. As a result, the mobility on
weekends is different from the one on workdays. Therefore, for each user, it
is reasonable to depict the trajectory for each day and then discover the
common mobility patterns shared by the trajectories of different days.

In our work, the human mobility is recorded by the GPS module embedded on the
smartphone devices. It should be emphasized that the GPS data (longitudes and
latitudes) is not evenly distributed spatially because one may stay longer at
a significant place (i.e, home or workplace/school) than at a less significant
place (i.e, restaurants or the roads). Thus, an appropriate description for
human mobility is to treat the location of an individual as a set of data
points randomly distributed in the space with respect to different
probabilities. Moreover, in practice, the data collecting procedure may not be
continuous all the time because the GPS module is turned off or does not
function sometimes. As a consequence, it arises the issue of data sparsity. In
general, the human mobility GPS data has the following properties:

  * •

The data has latent structures.

  * •

The data is not evenly distributed in space.

  * •

The data is sparse and noisy.

These unique data characteristics prevent researchers adopting some
conventional methods. Therefore, in our work, we adopt a probabilistic
approach to describe the daily human mobility. As compared to the conventional
methods, we believe our approach can explore more information from the
original GPS data and decrease the impact of data sparsity. The approach
presented in this paper is aimed at investigating this hypothesis and is
structured into three steps as shown in Fig. 1.

Figure 1: Overview of our method.

The first stage of the method is to estimate the probability density for each
day’s trajectories. For such task, Gaussian Mixture Model [4] is a possible
solution. However, the standard Gaussian Mixture Model needs to set the number
of components in advance, which is not practical to implement because
trajectory data can be statistically heterogeneous and a fixed component
number for every daily trajectories is not appropriate. To handle this
problem, we adopted the Infinite Gaussian Mixture Model [5], in which the
Dirichlet process prior is used to modify the mixed weights of components.

To measure the difference between different mobility probability densities,
the Kullback-Leibler divergence [6] estimator is used. The Kullback-Leibler
divergence is an asymmetric metric, which means the distance from distribution
p to distribution q is not the same as the distance from distribution q to
distribution p unless they are identical distributions. We exploit the
inequality property of KL divergence to reveal the subordinate relationship of
one trajectory to another.

Finally, we devise a clustering algorithm using the Infinite Gaussian Mixture
Models with Kullback-Leibler divergence to discover the mobility patterns
existing in human mobility data. More importantly, as compared to traditional
methods, our clustering algorithm is automatic because it does not require a
preset of the pattern number.

The work presented in this paper has then 3 main contributions:

  * •

For estimating probability density of daily mobility, we illustrate that the
Infinite Gaussian Mixture Models outperform the Gaussian Mixture Models.

  * •

We prove that Kullback-Leibler divergence is an appropriate metrics to measure
the closeness of mobility probability densities.

  * •

We develop a clustering algorithm based on Infinite Gaussian Mixture Model and
Kullback-Leibler divergence to find the human mobility patterns.

The reminder of the paper is organized as follows. Section II surveys the
related work. Section III addresses the problem we are tackling in this paper.
In Section IV, the proposed method is depicted. In Section V presents the
conducted experiment and its results to evaluate our method with real user
data. Finally, we conclude our paper and discuss about the future work in
Section VI.

##  II Related Work

In literature, previous research such as [7], [8], [9], [10], [11], [12] which
were studying human mobility with mobile data, mainly focused on tasks such as
extracting significant places, predicting next visiting places, predicting
visit duration or clustering trajectories.

A widespread topic is to predict human mobility with the smartphone usage
contextual information, e.g., temporal information, application usage, call
logs, WiFi status, Cell ID, etc. In [2] and [13] for instance, the researchers
applied various machine learning techniques to accomplish prediction tasks
such as next-time slot location prediction and next-place prediction. In
particular, they exploited how different combinations of contextual features
related to smartphone usage can affect prediction accuracy. Meanwhile, they
also compared the predicting performance of individual models and generic
models.

Another frequently-used method for such tasks is to use probabilistic models.
Through calculating the conditional probabilities between contextual features,
[14] developed the contextual conditional models for the next-place prediction
and visit duration prediction. In [15] and [16], the researchers presented the
probabilistic prediction frameworks based on kernel density estimation. [15]
utilized conditional kernels density estimation to predict the mobility events
while [16] devised different kernels for different context information types.
And in [3], the authors developed a location Hierarchical Dirichlet Process
(HDP) based approach to model heterogeneous location habits under data
sparsity.

In addition, generative models can also be applied to predict human mobility.
This type of models can be Naïve Bayes [17], Markov Model [18], Hidden Markov
Model (HMM) [19] or Dynamic Bayesian Network (DBN) [20], [21], etc. These
generative models attempt to predict the future states of human behavior by
computing the state transition probabilities. However, when the number of
states expands, the calculation grows exponentially.

Among the other possible approaches, [22] proposed a Hypertext Induced Topic
Search-based inference model for mining interesting locations and travel
sequences using a large GPS dataset in certain region. In [1], the authors
employed the random forests classifiers to label different places without any
geo-location information. [23] made use of nonlinear time series analysis of
the arrival time and residence time for location prediction.

In particular, for clustering user trajectories, there exists several
different methods. However, these conventional algorithms are not applicable
to our objectives. For example, some researchers used K-means [24], [25] in
their work, whereas K-means can not handle the trajectories with complex
shapes or noisy data because it is based on Euclidean distance. Besides, it
also need the pre-knowledge of cluster number, which is not acquirable in many
cases.

Though DBSCAN [26], [18], a density based clustering techniques, can deal with
data with arbitrary shapes and does not require the number of cluster in
advance. However, it still needs to set the minimum points number and
neighbourhood radius to recognize core areas and it treats the non-core data
points as noise. From our study, we argue that the trajectory parts with less
data density are also essential to demonstrate the human mobility
trajectories. And the grid searching algorithm [14] focus on detecting the
stay points within a set of square regions and fails to reveal the mobility at
a larger scale.

In this paper we decided to focus on a probabilistic point of view. As
compared to aforementioned previous works, we aim at describing the daily
trajectories using their probability densities. Moreover, to discover the
common mobility patterns shared among these trajectories, we devise an
automatic clustering algorithm. As opposed to traditional clustering
algorithms, our method is able to exploit more information from the sparse and
noisy original GPS data and free from pre-defining clusters number.

##  III Problem Formulation

As expressed in introduction, our purpose is to discover the mobility patterns
for each individual from their GPS location data.

As shown in Fig. 2, the mobility for one individual consists of many different
trajectories (the data is from the MDC dataset, the detailed data description
will be in following experiments). We believe that one’s daily mobility is
rather regular and there are common mobility patterns shared among different
daily trajectories. From a common sense, one may follow the regular daily
itineraries, for instance, home-work place/school-home. Yet, on different days
the daily itineraries may not be the same, for instance, on the way to home,
one may take a detour to do shopping in a supermarket sometimes. Hence, our
objective is to discover all the potential daily mobility from the data with
location information.

Figure 2: GPS data for a randomly selected user. Different colors represent
different days.

We extract each day’s trajectory from the whole dataset as shown in Fig. 3.
Fig. 3 reveals that a daily trajectories recorded by GPS data is not
distributed evenly in space, and is even not continuous in some areas. It may
be caused by the data collecting procedure: some data collecting time range is
actually relatively short (less than 24 hours, in fact, only few hours
sometimes), which leads to the data sparsity problem.

Figure 3: One randomly selected daily trajectories for a given user.

In order to overcome this problem and exploit as much information as possible
from the GPS data, we argue that a reasonable way to describe the daily
trajectories is to estimate the probability density of the location data. And
the relationship among the trajectories can be represented by their
probability densities. As a result, we can discover all the mobility patterns
for each user.

The tasks in this paper will be as follows:

  * •

Task 1: Estimate the probability density for mobility for each day. We will
compare results of GMM and IGMM.

  * •

Task 2: Measure the closeness between different trajectories. We will use the
KD divergence as metrics.

  * •

Task 3: Discover the similar mobility patterns among all the recorded daily
trajectories. This can be regarded as a clustering problem.

  * •

Task 4: Compare the IGMM algorithm with the GMM based algorithms.

  * •

Task 5: Identify the minimum data length for discovering all mobility
patterns.

##  IV Proposed Method

###  IV-A Estimate Daily Trajectories Probability Density

We assume that the GPS location data points are distributed randomly
spatially. Besides, the distribution of each day consists of unknown number of
heterogeneous sub-distributions. Therefore, it is reasonable to adopt mixed
Gaussian models for estimating probability density of daily mobility.

####  IV-A1 Gaussian Mixture Model

A Gaussian Mixture Model (GMM) is composed of a fixed number $K$ of sub
components. The probability distribution of a GMM can be described as follows:

$P(x)=\sum_{k=1}^{K}{\pi_{k}P(x|\theta_{k})}$ (1)

where, $x$ is the observable variable, $\pi_{k}$ is the assignment probability
for each model, with $\sum_{k=1}^{K}\pi_{k}=1,(0<\pi_{k}<1)$, and $\theta_{k}$
is the internal parameters of the base distribution.

Let $z_{n}$ be the latent variables for indicating categories.

$z_{n}\sim Categorical(z_{n}|\pi)\quad\sum_{k=1}^{K}z_{nk}=1$ (2)

where, $z_{n}=[z_{n1},z_{n2},...,z_{nk},...,z_{nK}]$, in which only one
element $z_{nk}=1$. It means $x_{n}$ is correspondent to $\theta_{k}$.

If the base distribution is a Gaussian, then:

$P(x|\theta_{k})=N(x|\mu_{k},\Lambda^{-1}_{k})$ (3)

where, $\mu_{k}$ is the mean vector and $\Lambda_{k}$ is the precision matrix.

Therefore, an observable sample $x_{n}$ is drawn from GMM according to:

$\ x_{n}\sim\prod_{k=1}^{K}{{N(x_{n}|\nu_{k},\Lambda_{k})}^{z_{nk}}}$ (4)

As it is illustrated above, one crucial issue of GMM is to pre-define the
number of components $K$. This is tricky because the probability distribution
for each day’s mobility is not identical. Thus, to define a fixed $K$ for all
mobility GMM models is not suitable in our case.

####  IV-A2 Infinite Gaussian Mixture Model

Alternatively, we resort to the Infinite Gaussian Mixture Model (IGMM) [5]. As
compared to finite Gaussian Mixture Model, by using a Dirichlet process (DP)
prior, IGMM does not need to specify the number of components in advance. Fig.
4 presents the graphical structure of the Infinite Gaussian Mixture Model.

Figure 4: The plate representation of Infinite Gaussian Mixture Model.

In Fig. 4, the nodes represents the random variables and especially, the
shaded node is observable and the unshaded nodes are unobservable. The edges
represent the conditional dependencies between variables. And the variables
are within the plates means that they are drawn repeatedly.

According to Fig. 4, the Dirichlet process can be depicted as:

$G\sim DP(\alpha,G_{0})$ (5)

where, $G$ is a random measure, which consists of infinite base measure
$G_{0}$ and $\lambda$ is the hyper-parameter of $G_{0}$. In our case, it is a
series of Gaussian distributions. And $\alpha\sim Gamma(1,1)$ is the
concentration parameter. $N$ is the total samples number. $\theta_{k}$ is the
parameters of base distribution. $X_{k}$ is the observable data for
$\theta_{k}$. $Z_{k}$ is the latent variables that indicates the category of
$X_{k}$.

Alternatively, $G$ can be explicitly depicted as follow:

$G(\theta)=\sum_{k=1}^{\infty}\pi_{k}\delta_{\theta_{k}}$ (6)

where, $\theta_{k}\sim G_{0}(\lambda)$, and $\delta$ is Dirac function.
$\pi_{k}$ determines the proportion weights of the clusters and the
$\delta_{\theta_{k}}$ is the prior of the $\theta_{k}$ to determine the
location of clusters in space.

We choose the Stick-breaking process (SBP) [27] to implement the Dirichlet
process as the prior for $\pi_{k}$. The the Stick-breaking process can be
described as follow:

$\pi_{k}=\nu_{k}\prod_{j=1}^{k-1}{(1-\nu_{j})}\quad k\geq 2$ (7)

where, $\nu_{k}\sim Beta(1,\alpha)$.

Since $P(x|\theta)$ is Gaussian, $\theta=\\{\mu,\Lambda\\}$. Further, let
$G_{0}$ be a Gaussian-Wishart distribution, then, $\mu_{k},\Lambda_{k}\sim
G_{0}(\mu,\Lambda)$. Therefore, similarly, draw an observable sample $x_{n}$
from IGMM:

$\
x_{n}\sim\prod_{k=1}^{\infty}{{N(x_{n}|\nu_{k},\Lambda^{-1}_{k})}^{z_{nk}}}$
(8)

Then, Variational Inference is used to solve the IGMM models. As compared to
Gibbs sampling or to a Markov chain Monte Carlo (MCMC) method which consumes a
large mount of calculating time, Variational Inference is relatively fast
[28]. The results will be demonstrated in the later experiments.

###  IV-B Measure Daily Trajectories Similarities

The Kullback-Leibler (KL) divergence is a metric to evaluate the closeness
between two distributions. For continuous variables, the KL divergence
$D_{KL}(p||q)$ the expectation of the logarithmic difference between the $p$
and $q$ with respect to probability $p$ and vice versa. From (9) and (10), it
can be seen that the KL divergence is non-negative and asymmetric. In many
occasions, the inequality of the KL divergence is notorious. However, in our
methodology, on the contrary, we take advantage of the characteristics of
inequality to reveal the similarities among different trajectories instead of
the Jensen-Shannon divergence which is a symmetric metrics.

$D_{KL}(p||q)=\int_{-\infty}^{\infty}p(x)log(\frac{p(x)}{q(x)})dx$ (9)

$D_{KL}(q||p)=\int_{-\infty}^{\infty}q(y)log(\frac{q(y)}{p(y)})dy$ (10)

There is no closed form to implement the KL divergence by the definition of
(9) and (10) for Gaussian Mixture Models. Instead, we resort to the Monte
Carlo simulation method proposed in [29]. Then, the KL divergence can be
caculated by:

$D_{KL_{MC}}(p||q)=\frac{1}{n}\sum_{i=1}^{n}log(\frac{p(x_{i})}{q(x_{i})})$
(11)

$D_{KL_{MC}}(q||p)=\frac{1}{n}\sum_{i=1}^{n}log(\frac{q(y_{i})}{p(y_{i})})$
(12)

This method is to draw a large amount of i.i.d samples $x_{i}$ from
distribution $p$ to calculate $D_{KL_{MC}}(p||q)$ according to (11) and
$D_{KL_{MC}}(p||q)\to D_{KL}(p||q)$ as $n\to\infty$. And it is the same for
implementing (10) by using (12). The results will be demonstrated in the later
experiments. Furthermore, if we define a representative trajectory for a
mobility pattern then we can distinguish whether a new trajectory belong to
this cluster by comparing it to the representative trajectory. To do so, we
need to set a threshold with a lower bound and an upper bound for the KL
divergence, then it can be used as the metrics to cluster mobility patterns.

###  IV-C Discover Mobility Patterns

As mentioned before, our task is to find the trajectories which are mutually
similar. For this reason, we treat the different mobility patterns as
different clusters in which the daily trajectories are their sub-members. Even
so, the trajectories within the same clusters still can not be treated as
identically distributed as other conventional clustering methods because of
different trajectory lengths. Hence, we need to devise a algorithm that is
able to cluster the trajectories based on the distribution similarity and the
aforementioned KL divergence can be applicable as closeness metrics. Note that
due to the large data scale and the number of the potential clusters, a high
accuracy solution is intractable sometimes. Therefore, instead of pursuing a
very accurate result, our purpose is to reach a relative accurate result in a
reasonable amount of calculating time.

The algorithm we devise is shown in Algorithm 1 and its variables are
described in Table I.

The first step of the clustering algorithm is to calculate the probability
densities using the Infinite Gaussian Mixture Models. At this step, we create
a list, in which the members are the probability densities of each. Then, the
first cluster is created with one trajectory as its first member and it also
will be compared with other trajectories.

Afterwards, we select another daily trajectory in the list and calculate the
KL divergences, both $D_{KL(p||q)}$ and $D_{KL(q||p)}$. And the new trajectory
is added to the current cluster if the minimum and maximum of the KL-
divergences are smaller than the lower bound and upper bound of the thresholds
at the same time, respectively. And if the $D_{KL(p||q)}$ is smaller than
$D_{KL(q||p)}$, the new trajectory become the benchmark for the current
cluster. An alternative way to do this is to compute the probability density
of the current cluster using all the data of the discovered trajectories,
however, the calculation will be massive.

This step will be repeated until all the trajectories belonging to the current
cluster are discovered at the end of this iteration. Then, all the members of
the current cluster are removed from iteration because, we assume that each
trajectories can only be a member of one mobility pattern. At the start of new
iteration, a new cluster is created, repeat the above steps until the list is
empty. Finally, all the mobility patterns are discovered.

TABLE I: Variables Description Variable | Domain | Description
---|---|---
$d$ | $\\{1,2,\dots,D\\}$ | Number of data collecting day
$X$ | $\\{X_{1},X_{2},\dots,X_{d},\dots,X_{D}\\}$ | Total GPS data (longitudes, latitudes
$P$ | $\\{P_{1},P_{2},\dots,P_{d},\dots,P_{D}\\}$ | Probability density for $X$
$M$ | $\\{M_{1},M_{2},\dots,M_{k},\dots M_{K}\\}$ | Total mobility patterns
$M_{k}$ | $\\{X_{k1},X_{k2},\dots,X_{kn}\\}$ | Discovered mobility pattern
$Th$ | $\\{lowerbound,upperbound\\}$ | Threshold for distinguishing patterns
$D_{KL}$ | $\\{D_{KL(p||q)},D_{KL(q||p)}\\}$ | KL divergences
Algorithm 1 Mobility Pattern Discovering Algorithm

1:$X$

2:$M$

3:

4:$P$ $\leftarrow$ IGMM($X$) $\triangleright$ probability density estimation

5:

6:Initialize:$M=\\{M_{k}\\}$ $\triangleright$ create the mobility patterns set

7:while $P\neq\emptyset$ do

8: $X_{s}=X_{1}$ $\triangleright$ set the baseline mobility for $M_{k}$

9: $M_{k}=\\{X_{s}\\}$ $\triangleright$ create current pattern $M_{k}$

10: for $d=2,\dots,D$ do

11: $D_{KL}$ $\leftarrow$ $(P_{s},P_{d})$$\triangleright$ measure similarity

12: if ($\min(D_{KL})<Th[0]$) & ($\max(D_{KL})<Th[1]$) then$\triangleright$
two patterns are similar

13: add $P_{d}$ to $M_{k}$$\triangleright$ add new member

14: if $D_{KL}[0]>D_{KL}[1]$ then

15: $P_{s}$ $\leftarrow$ $P_{d}$ $\triangleright$ change the baseline mobility

16: end if

17: end if

18: end for

19: remove $P_{d}\in M_{k}$ from $P$$\triangleright$ current pattern is
finished

20:

21: create $M_{k+1}$ $\triangleright$ find new mobility pattern

22: add $M_{k+1}$ to $M$

23:end while

24:

25:return $M$

As it can be seen that our algorithm is designed to discover the latent
mobility patterns automatically without the pre-knowledge of the numbers of
existing patterns.

##  V Experiments and Results

###  V-A Dataset Description

We use the Mobile Data Challenge (MDC) dataset [30], [31] to validate our
method. This dataset records comprehensive smartphone usages with fine
granularity of time. The participants of the MDC dataset are up to nearly
$200$ and the data collection campaign lasts more than $18$ months. This
abundant information thus can be used to investigate individual mobility
patterns.

To collect the individual location information, as compared to other methods,
for instance, through stand-alone GPS devices, using GPS-equipped smartphones
is a more practical way to have a larger group of participants without
affecting their daily life.

In our study, we attempt to find the trajectories that belong to the same
mobility patterns, thus we focus the spatial information of the GPS records,
namely, the latitudes and longitudes and the time-stamps of the data are not
considered. Meanwhile, since we consider not only the significant places but
all location records, we use the unlabeled data without any semantic
information.

###  V-B Experimental Setup

In the conducted experiment, we randomly select $20$ users with sufficient
data. And each user’s is segmented by the time range of one day. Fig. 5
demonstrates the number of data collecting days for each user. It can be seen
that the data collecting days for most users are more than $200$. And with
such amount of data, we believe that it is possible to discover individual’s
mobility patterns from it.

Figure 5: Number of data collecting days for each user

However, as it is illustrated in Fig. 6, the data length of each day varies
from less than $4$ hours to $24$ hours. And most of them is less than $8$
hours. Hence, we also should be aware that some data can be missing because
the GPS modules were turned off or were not functioning. Consequently, it is
one of the reasons that cause the data sparsity problem. In the following
part, we will prove that our method can mitigate the impact of data sparsity.

Table II summarizes the temporal information about the GPS data for conducting
the experiments.

TABLE II: Data collecting time  | Average | Total
---|---|---
Collecting days for all users | 300.25 days | 6005.0 days
Collecting hours per day for all users | 6.93 hours | 41595.0 hours
Collecting hours per day for each user | 6.67 hours | 2084.65 hours
Figure 6: Empirical cumulative distribution of hours per data collecting days.

To test the performance of our method, we will conduct three experiments from
different perspectives:

  * •

We compare the IGMM model with the GMM model on estimating the daily
trajectories probability density.

  * •

We use that the KL divergence to measure the closeness of different
trajectories.

  * •

We test our method on each selected user data so as to find the daily mobility
patterns for each individual.

  * •

We compare the results of the IGMM models to a series of fixed-number
components GMM models.

  * •

We run the algorithm on the varying-length datasets, in the aim to find the
minimum data length for discovering most mobility patterns of one individual.

###  V-C Experimental Results

####  V-C1 Probability Density Estimation

Figure 7: Distribution estimation by GMM (negative log-likelihood) Figure 8:
Distribution estimation by IGMM (negative log-likelihood)

Fig. 7 and Fig. 8 show the density estimation results obtained by GMM and
IGMM, respectively. It can be seen that, compared to the GMM model, the result
of the IGMM model is more smooth. It suggests that IGMM is not affected by the
number of components and it infers more information from the original data and
it is less influenced by data sparsity. That is to say, on the same dataset,
the computational results of IGMM have higher fidelity. Hence, in our
approach, we chose IGMM to estimate probability density of daily mobility.

####  V-C2 Measuring Daily Trajectories Similarities

As shown in Fig. 9, we select $5$ daily trajectories from the data of one
random user to present the KD divergences between different trajectories. The
baseline trajectory is the Trajectory 1 and the rest of trajectories are
chosen to make comparisons.

(a) Trajectory 1 and Trajectory 2

(b) Trajectory 1 and Trajectory 3

(c) Trajectory 1 and Trajectory 4

(d) Trajectory 1 and Trajectory 5

Figure 9: Comparison Between Different Trajectories. TABLE III: KL-Divergences for Different Trajectories. p | q | $D_{KL(p||q)}$ | $D_{KL(p||q)}$
---|---|---|---
Trajectory 1 | Trajectory 2 | 7.21 | 2.82
Trajectory 1 | Trajectory 3 | 1.28 | 1.83
Trajectory 1 | Trajectory 4 | 19.07 | 1269.47
Trajectory 1 | Trajectory 5 | 3.08 | 996.17

The combinations are shown in Figure 9 and the results are illustrated in
Table III:

Trajectory 2 is nearly a subset of Trajectory 1 and thus $D_{KL}(p||q)$ is
larger than $D_{KL(p||q)}$. And their values are both small, thus Trajectory 2
and Trajectory 1 can be regarded to belong to the same mobility pattern.
Trajectory 3 is very similar to Trajectory 1 and $D_{KL}(p||q)$ almost equals
to $D_{KL}(q||p)$. Hence, they also are the members of the same mobility
pattern. Trajectory 4 share a small part with Trajectory 1 whereas generally
they are very different. $D_{KL}(p||q)$ and $D_{KL}(q||p)$ are both very
large. Therefore, it is reasonable to recognize Trajectory 4 and Trajectory 1
as different patterns. And Trajectory 5 is totally different from Trajectory
1. And $D_{KL}(p||q)$ is small but $D_{KL}(p||q)$ are very large. So they
naturally are not in the same pattern. According to the trajectories in the
Fig. 9 and the results in Table III, it shows that the KL divergence is able
to illustrate the difference among trajectories and can be the metrics for
clustering.

####  V-C3 Discovering Daily Mobility Patterns

We run our algorithm on the data of the $20$ users to discover their daily
mobility patterns.

  * •

Discovered Patterns: The partial results for different randomly selected user
data are demonstrated in Fig. 10. It shows that, after clustered by our
proposed algorithm, the data is split into different mobility patterns. Each
cluster is composed of trajectories close to each other even if they are not
distributed with the same density in the space. That proves our methodology is
able to find the different mobility patterns even under the condition of noisy
data and discontinuous trajectories.

Figure 10: Discovered mobility patterns from three random selected users.
Different colors represents different days.

Fig. 11 shows that our methodology is not only able to identify the different
patterns in the daily trajectories data but is also able to find the most
representative trajectories for each mobility pattern.

Figure 11: Representative trajectories for each discovered mobility patterns.

  * •

Number of Patterns and Trajectories: Fig. 12 shows the number of discovered
mobility pattern for all the user in our experiments. We can see that the
number of mobility patterns varies from $5$ to more than $30$ and most of them
are about $10$ to $15$. It also can be found that the lengths of data
collecting days are not proportional to the number of discovered mobility
patterns, which indicates that the results rely more on the individual
behavior rather than the data length.

Figure 12: Number of discovered mobility patterns for each user.

  * •

Number of members for each patterns: Fig. 13 depicts the number of members for
each discovered mobility patterns for all users. We can see that most mobility
patterns consist of less than $50$ trajectories. And nearly $40\%$ of the
patterns have only one trajectory, whereas few patterns have more than $100$
trajectories.

Figure 13: Empirical cumulative distribution of discovered patterns’ members.

One needs to notice that the number of discovered patterns depends on the
Kullback-Leibler divergence thresholds we set in the clustering algorithm.
When the thresholds are small, it means that the condition to be in the same
mobility pattern is more strict and naturally the discovered mobility patterns
are more and the member of each patterns are less, and vice versa.

####  V-C4 Comparison to GMM

To compare the Infinite Gaussian Mixture Models, we use a group of Gaussian
Mixture Models with different numbers of components to estimate the daily
mobility probability densities in our proposed clustering algorithm. The
metrics we adopt to evaluate the results is the mean log-likelihood. The
results show in Table IV that changing the fixed number of component Gaussian
Mixture Models can not enhance the clustering performance. On the contrary,
the Infinite Gaussian Mixture Models can improve the clustering performance.

TABLE IV: Overall Mean Log-likelihood For Different Models Model | Mean log-likelihood
---|---
GMM-1 | -26078.15
GMM-2 | -38514.32
GMM-3 | -52431.62
GMM-4 | -63794.70
GMM-5 | -73508.10
IGMM | -24871.78

####  V-C5 Varying Data Length

To investigate how the data length, namely, the number of days of the data,
affects the results, we utilize different data lengths which varies from $50$
days to $350$ days. The results are shown in Fig. 14. It can be seen that,
from $50$-day data length to $200$-day data length, the average discovered
mobility pattern numbers increase as the data length grows. While, when the
data length is larger than $200$ days, the patterns numbers change marginally.
Thus, according to the results, we can say that, generally, a $200$-day GPS
dataset is large enough to discover most of the mobility patterns of an
individual.

Figure 14: Average number of discovered patterns for different data collecting
day length. Error bars represent the standard deviation.

##  VI Conclusion and Perspective

In this work, we presented a probabilistic approach to discover human daily
mobility patterns based on GPS data collected by smartphones.

In our approach, the human daily mobility is considered as sets of probability
distributions. The proposed approach is divided into three parts. The first
step is to estimate the probability densities. We argue that the Infinite
Gaussian Mixture Model is more appropriate than the standard Gaussian Mixture
Model to this issue, this argument being besides validated by the experimental
results. Further, in order to find the similar trajectories, one needs to
measure the closeness between the trajectories. For this task, we chose the
Kullback-Leibler divergence as distance metrics. According to the
computational results from the selected trajectories, we validated, on test
sets, that KL divergence is able to measure the similarities among the
trajectories. Finally, we devised a novel automatic clustering algorithm
combining the advantages of both IGMM and the KL divergence so as to discover
human daily mobility patterns without having the knowledge of the cluster
number in advance.

For validation, we select $20$ random individual data from the MDC dataset to
conduct the different experiments. The results obtained show that our proposed
approach can discern different mobility patterns and select the most
representative trajectories for each mobility patterns from the GPS data. In
addition, we also compared the IGMM based algorithm with a group of GMM based
algorithms with various fix-number components, the results reveal that the
IGMM model performs better. Finally, testing varying-length dataset on our
methods leads to results which suggest that a $200$-day GPS is generally
sufficient enough to discover most of the individual daily mobility patterns.

We are aware of that human mobility is also a time-related behavior. Thus, as
future work, we plan to take into account the temporal information, for
example, hour of day and day of week. Based on that, we will try to build a
spatial-temporal probabilistic model to predict human mobility. In addition,
for further study, we may exploit other smartphone usage information (i.g.,
application usage) in the dataset to obtain more knowledge about human
behavior.

## Acknowledgment

The research in this paper used the MDC Database made by Idiap Research
Institute, Switzerland and owned by Nokia. The authors would like to thank the
MDC team for providing the access to the database.

## References

  * [1] T. M. T. Do and D. Gatica-Perez, “The places of our lives: Visiting patterns and automatic labeling from longitudinal smartphone data,” _IEEE Transactions on Mobile Computing_ , vol. 13, no. 3, pp. 638–648, 2014.
  * [2] P. Baumann, C. Koehler, A. K. Dey, and S. Santini, “Selecting individual and population models for predicting human mobility,” _IEEE Transactions on Mobile Computing_ , vol. 17, no. 10, pp. 2408–2422, 2018.
  * [3] J. McInerney, J. Zheng, A. Rogers, and N. R. Jennings, “Modelling heterogeneous location habits in human populations for location prediction under data sparsity,” in _Proceedings of the 2013 ACM international joint conference on Pervasive and ubiquitous computing_. ACM, 2013, pp. 469–478.
  * [4] D. Reynolds, “Gaussian mixture models,” _Encyclopedia of biometrics_ , pp. 827–832, 2015.
  * [5] C. E. Rasmussen, “The infinite gaussian mixture model,” in _Advances in neural information processing systems_ , 2000, pp. 554–560.
  * [6] S. Kullback, _Information theory and statistics_. Courier Corporation, 1997.
  * [7] X. Lu, E. Wetter, N. Bharti, A. J. Tatem, and L. Bengtsson, “Approaching the limit of predictability in human mobility,” _Scientific reports_ , vol. 3, p. 2923, 2013.
  * [8] J. Ye, S. Dobson, and S. McKeever, “Situation identification techniques in pervasive computing: A review,” _Pervasive and mobile computing_ , vol. 8, no. 1, pp. 36–66, 2012.
  * [9] M. Lin and W.-J. Hsu, “Mining gps data for mobility patterns: A survey,” _Pervasive and mobile computing_ , vol. 12, pp. 1–16, 2014.
  * [10] P. Pirozmand, G. Wu, B. Jedari, and F. Xia, “Human mobility in opportunistic networks: Characteristics, models and prediction methods,” _Journal of Network and Computer Applications_ , vol. 42, pp. 45–58, 2014.
  * [11] Y. Zheng, “Trajectory data mining: an overview,” _ACM Transactions on Intelligent Systems and Technology (TIST)_ , vol. 6, no. 3, p. 29, 2015.
  * [12] H. Cao, N. Mamoulis, and D. W. Cheung, “Discovery of periodic patterns in spatiotemporal sequences,” _IEEE Transactions on Knowledge and Data Engineering_ , vol. 19, no. 4, pp. 453–467, 2007.
  * [13] T. M. T. Do and D. Gatica-Perez, “Where and what: Using smartphones to predict next locations and applications in daily life,” _Pervasive and Mobile Computing_ , vol. 12, pp. 79–91, 2014.
  * [14] T.-M.-T. Do and D. Gatica-Perez, “Contextual conditional models for smartphone-based human mobility prediction,” in _Proceedings of the 2012 ACM conference on ubiquitous computing_. ACM, 2012, pp. 163–172.
  * [15] T. M. T. Do, O. Dousse, M. Miettinen, and D. Gatica-Perez, “A probabilistic kernel method for human mobility prediction with smartphones,” _Pervasive and Mobile Computing_ , vol. 20, pp. 13–28, 2015.
  * [16] A. Peddemors, H. Eertink, and I. Niemegeers, “Predicting mobility events on personal devices,” _Pervasive and Mobile Computing_ , vol. 6, no. 4, pp. 401–423, 2010.
  * [17] M. Muhlenbrock, O. Brdiczka, D. Snowdon, and J.-L. Meunier, “Learning to detect user activity and availability from a variety of sensor data,” in _Second IEEE Annual Conference on Pervasive Computing and Communications, 2004. Proceedings of the_. IEEE, 2004, pp. 13–22.
  * [18] C. Yu, Y. Liu, D. Yao, L. T. Yang, H. Jin, H. Chen, and Q. Ding, “Modeling user activity patterns for next-place prediction,” _IEEE Systems Journal_ , vol. 11, no. 2, pp. 1060–1071, 2017.
  * [19] S.-B. Cho, “Exploiting machine learning techniques for location recognition and prediction with smartphone logs,” _Neurocomputing_ , vol. 176, pp. 98–106, 2016.
  * [20] V. Etter, M. Kafsi, E. Kazemi, M. Grossglauser, and P. Thiran, “Where to go from here? mobility prediction from instantaneous information,” _Pervasive and Mobile Computing_ , vol. 9, no. 6, pp. 784–797, 2013.
  * [21] D. J. Patterson, L. Liao, D. Fox, and H. Kautz, “Inferring high-level behavior from low-level sensors,” in _International Conference on Ubiquitous Computing_. Springer, 2003, pp. 73–89.
  * [22] Y. Zheng, L. Zhang, X. Xie, and W.-Y. Ma, “Mining interesting locations and travel sequences from gps trajectories,” in _Proceedings of the 18th international conference on World wide web_. ACM, 2009, pp. 791–800.
  * [23] S. Scellato, M. Musolesi, C. Mascolo, V. Latora, and A. T. Campbell, “Nextplace: a spatio-temporal prediction framework for pervasive systems,” in _International Conference on Pervasive Computing_. Springer, 2011, pp. 152–169.
  * [24] S. Jiang, J. Ferreira, and M. C. González, “Clustering daily patterns of human activities in the city,” _Data Mining and Knowledge Discovery_ , vol. 25, no. 3, pp. 478–510, 2012.
  * [25] D. Ashbrook and T. Starner, “Using gps to learn significant locations and predict movement across multiple users,” _Personal and Ubiquitous computing_ , vol. 7, no. 5, pp. 275–286, 2003.
  * [26] J. Tang, F. Liu, Y. Wang, and H. Wang, “Uncovering urban human mobility from large scale taxi gps data,” _Physica A: Statistical Mechanics and its Applications_ , vol. 438, pp. 140–153, 2015.
  * [27] J. Sethuraman, “A constructive definition of dirichlet priors,” _Statistica sinica_ , pp. 639–650, 1994.
  * [28] D. M. Blei, M. I. Jordan _et al._ , “Variational inference for dirichlet process mixtures,” _Bayesian analysis_ , vol. 1, no. 1, pp. 121–143, 2006\.
  * [29] J. R. Hershey and P. A. Olsen, “Approximating the kullback leibler divergence between gaussian mixture models,” in _2007 IEEE International Conference on Acoustics, Speech and Signal Processing-ICASSP’07_ , vol. 4. IEEE, 2007, pp. IV–317.
  * [30] N. Kiukkonen, J. Blom, O. Dousse, D. Gatica-Perez, and J. Laurila, “Towards rich mobile phone datasets: Lausanne data collection campaign,” _Proc. ICPS, Berlin_ , vol. 68, 2010.
  * [31] J. K. Laurila, D. Gatica-Perez, I. Aad, O. Bornet, T.-M.-T. Do, O. Dousse, J. Eberle, M. Miettinen _et al._ , “The mobile data challenge: Big data for mobile computing research,” Tech. Rep., 2012.


contents=  ©2023 IEEE. Personal use of this material is permitted. Permission
from IEEE must be obtained for all other uses, in any current or future media,
including reprinting/republishing this material for advertising or promotional
purposes, creating new collective works for resale or redistribution to
servers or lists, or reuse of any copyrighted component of this work in other
works. , scale=1, angle=0, position=0.38515mm, color=black, anchor=below

# Relating EEG recordings to speech using envelope tracking and the speech-FFR

###### Abstract

During speech perception, a listener’s electroencephalogram (EEG) reflects
acoustic-level processing as well as higher-level cognitive factors such as
speech comprehension and attention. However, decoding speech from EEG
recordings is challenging due to the low signal-to-noise ratios of EEG
signals. We report on an approach developed for the ICASSP 2023 ‘Auditory EEG
Decoding’ Signal Processing Grand Challenge. A simple ensembling method is
shown to considerably improve upon the baseline decoder performance. Even
higher classification rates are achieved by jointly decoding the speech-evoked
frequency-following response and responses to the temporal envelope of speech,
as well as by fine-tuning the decoders to individual subjects. Our results
could have applications in the diagnosis of hearing disorders or in
cognitively steered hearing aids.

Index Terms—  EEG decoding, deep learning, speech

##  1 Introduction

When people listen to speech, their brainwaves synchronise with acoustic
features such as the speech envelope. The degree of this neural tracking
reflects cognitive factors such as attention to speech, speech comprehension,
and intelligibility [1, 2]. Accurately decoding speech from
electroencephalography (EEG) is, however, a challenging task owing to the low
signal-to-noise ratios of EEG recordings and the limited availability of EEG
data recorded during speech perception.

Neural speech tracking in a particular subject is best assessed when EEG
responses to speech from that particular subject are available and can be used
to train a subject-specific decoder. Decoding EEG responses to speech from
unseen subjects is a harder task, since EEG signals vary greatly between
individuals. The ICASSP 2023 ‘Auditory EEG Decoding’ Signal Processing Grand
Challenge involves both types of decoding. Its task is to develop population
match-mismatch decoders: given a temporal segment of EEG data and two
candidate speech segments, the decoder should predict which of the speech
segments corresponds to the EEG signal. This should be done both for EEG
recordings from subjects that have been included in the training set, and for
others whose EEG data has not been seen in the training stage.

Here we describe the development of our decoders which placed first in this
match-mismatch task. Besides the speech envelope, we relate a second acoustic
feature to the EEG recordings. This feature is the temporal fine structure of
the voiced parts of speech (which consist of a fundamental frequency and many
higher harmonics). The electroencephalogram displays a strong response at the
fundamental frequency of speech, termed the speech-frequency-following
response or speech-FFR [3]. The speech-FFR is driven both by the fundamental
frequency itself as well as by higher harmonics, and can be decoded from EEG
responses to continuous speech [3, 4, 5, 6].

##  2 Materials and methods

2.1. Dataset. A large training dataset was provided by the ICASSP competition
organisers, which consisted of EEG recordings from 71 subjects who listened to
speech material [7]. Decoders were evaluated against a heldout dataset
comprising EEG from 70 subjects included in the training dataset, and 15 new
unseen subjects.

We used two pre-processed versions of the dataset that contained the two
speech features of interest together with the corresponding EEG signals. The
first dataset consisted of speech envelopes and EEG recordings sampled at
$64\text{\,}\mathrm{Hz}$. The second dataset contained the envelope
modulations of the higher harmonics of the fundamental frequency of the voiced
parts of speech, together with EEG sampled at $512\text{\,}\mathrm{Hz}$ [6].
Pre-processing and methodological details will be described in a forthcoming
publication [8].

2.2. Baseline decoder. The baseline decoder is a deep neural network (DNN)
proposed by Accou et al. [2]. This DNN brings the speech envelopes and EEG
recordings into a space where the matched envelope and EEG segments are
maximally similar. The output layer consists of a single Sigmoid neuron. The
matched and mismatched envelopes are presented as an ordered pair, and the
Sigmoid neuron predicts the probability that the first envelope is matched to
the EEG. The decoder is trained with the binary crossentropy loss function and
Adam optimizer without regularisation.

2.3. Baseline + speech-FFR decoder. The baseline decoder relates the EEG
recordings to the temporal envelope of speech only. We were also interested in
relating the EEG signals to the temporal fine structure of speech, since this
can improve classification performance [9]. We retained the architecture of
the baseline decoder, but swapped the speech envelopes for the high-frequency
envelope-modulations feature [6]. The Sigmoid outputs of this speech-FFR
decoder and the baseline decoder were combined via linear discriminant
analysis (LDA) to produce the final predicted label.

2.4. Decoder training and fine-tuning. Training examples were presented to the
decoders as temporal segments of $3\text{\,}\mathrm{s}$ in duration (the same
duration was used for evaluation). Sources of randomness in the training
procedure include the decoder initialisation, and the order in which training
examples were presented. The effects of these were marginalised by averaging
the Sigmoid outputs of several trained instances of the decoders. For the
population decoders, the hop length between the onsets of the training
examples was $1\text{\,}\mathrm{s}$. When fine-tuning the population decoder
to individual subjects, this was reduced to $0.125\text{\,}\mathrm{s}$, and
regularisation (batch normalisation and spatial dropout) was applied to the
input layers.

##  3 Results

3.1. Averaging of decoder outputs. We trained 100 instances of the baseline
decoder. By averaging the Sigmoid outputs of the instances, the classification
accuracy was improved (Figure 1). Therefore, we formed two averaged population
decoders: these used 50 instances of the baseline decoder, and 30 instances of
the speech-FFR decoder, respectively. The averaged decoders were used in
Section 3.2.

Fig. 1: Classification accuracies against number of averaged decoders (log
scale). The bootstrapped mean and range of the accuracies are denoted by the
dots and errorbars respectively.

3.2. Combining averaged baseline and speech-FFR decoders. The baseline decoder
generally achieves higher accuracies than the speech-FFR decoder (Figure 2,
left). The LDA classifier was fitted on an unseen portion of the training
dataset, for which the correlation between the Sigmoid outputs of the decoders
was moderate ($R=0.229$, Figure 2, right). For the heldout dataset the
correlation was similar ($R=0.206$), confirming that there was no severe
overfitting or distributional shift. This composite decoder achieved an
accuracy of 81.18% on the heldout dataset for unseen subjects.

Fig. 2: Comparison between the averaged baseline and speech-FFR decoders.
(Left) The accuracies of the two decoders are shown for the 71 subjects.
(Right) The Sigmoid outputs of the averaged decoders are shown. A red (blue)
point indicates an EEG segment that matches the first (second) speech segment
(recall that these are grouped as an ordered pair). The LDA decision boundary
is shown.

3.3. Decoder fine-tuning. Figure 3 shows the effect of fine-tuning the
decoders to individual subjects. For both the baseline decoders and the
speech-FFR decoders, fine-tuning significantly improved the classification
accuracy ($p<0.001$, signed-rank tests). The fine-tuned baseline decoders
achieved an accuracy of 82.71% on the heldout dataset for seen subjects.

Fig. 3: Effect of fine-tuning the decoders to individual subjects. Each
datapoint represents the classification accuracy for a single subject.

##  4 Conclusions

The performance of the baseline decoder could be improved by averaging the
outputs of several trained decoder instances. Combining the averaged baseline
and speech-FFR decoders enhanced the decoding accuracies further. The best
results were obtained by using versions of the baseline decoder which were
fine-tuned to individual subjects where possible, and by using the composite
decoder for unseen subjects. This approach won the match-mismatch task of the
ICASSP 2023 Signal Processing Grand Challenge ‘Auditory EEG Decoding’. Future
work will establish which aspects of the fine-tuning procedure led to such
high classification accuracies.

## References

  * [1] O. Etard and T. Reichenbach,  “Neural speech tracking in the theta and in the delta frequency band differentially encode clarity and comprehension of speech in noise,”  J. Neurosci., vol. 39, no. 29, pp. 5750–5759, 2019.
  * [2] B. Accou et al.,  “Predicting speech intelligibility from EEG in a non-linear classification paradigm,”  J. Neural Eng., vol. 18, no. 6, pp. 066008, 2021.
  * [3] N. Kraus et al.,  “The frequency-following response: a window into human communication,”  in The Frequency-Following Response, pp. 1–15. Springer, 2017.
  * [4] A. E. Forte et al.,  “The human auditory brainstem response to running speech reveals a subcortical mechanism for selective attention,”  eLife, vol. 6, pp. e27203, 2017.
  * [5] O. Etard et al.,  “Decoding of selective attention to continuous speech from the human auditory brainstem response,”  NeuroImage, vol. 200, pp. 1–11, 2019.
  * [6] P. Kulasingham et al.,  “High gamma cortical processing of continuous speech in younger and older listeners,”  NeuroImage, vol. 222, pp. 117291, 2020.
  * [7] L. Bollens et al.,  “A Large Auditory EEG decoding dataset,”, doi: 10.48804/K3VSND 2023\.
  * [8] M. Thornton et al.,  “Relating EEG with speech using envelope tracking and the speech-FFR,”  in preparation.
  * [9] C. Puffay et al.,  “Relating the fundamental frequency of speech with EEG using a dilated convolutional network,”  in Proc. Interspeech 2022, pp. 4038–4042.