On the Inability of Markov Models to Capture Criticality in Human Mobility
Vaibhav Kulkarni Abhijit Mahalunkar Benoit Garbinato John D. Kelleher










We examine the non-Markovian nature of human mobility by exposing the inability of Markov models to capture criticality in human mobility. 
In particular, the assumed Markovian nature of mobility was used to establish a theoretical upper bound on the predictability of human mobility (expressed as a minimum error probability limit), based on temporally correlated entropy. 
Since its inception, this bound has been widely used and empirically validated using Markov chains.
We show that recurrent-neural architectures can achieve significantly higher predictability, surpassing this widely used upper bound.
In order to explain this anomaly, we shed light on several underlying assumptions in previous research works that has resulted in this bias.
By evaluating the mobility predictability on real-world datasets, we show that human mobility exhibits scale-invariant long-range correlations, bearing similarity to a power-law decay.
This is in contrast to the initial assumption that human mobility follows an exponential decay. 
This assumption of exponential decay coupled with Lempel-Ziv compression in computing Fano's inequality has led to an inaccurate estimation of the predictability upper bound. 
We show that this approach inflates the entropy, consequently lowering the upper bound on human mobility predictability.
We finally highlight that this approach tends to overlook long-range correlations in human mobility. 
This explains why recurrent-neural architectures that are designed to handle long-range structural correlations surpass the previously computed upper bound on mobility predictability. 










Exploring Sensitivity to NMSSM Signatures with Low Missing Transverse Energy at the LHC
    A. Titterton^1,2,3, U. Ellwanger^4, H.U. Flaecher^1, S. Moretti^2,3,
 C.H. Shepherd-Themistocleous^3
    December 30, 2023
========================================================================================================





§ INTRODUCTION



The rapid proliferation of mobile devices equipped with internet connectivity and positioning systems has resulted in the collection of massive amounts of human-mobility data.
Real-time user locations are typically collected using the Global Positioning System (GPS), Call Detail Record logs (CDR) and Wireless-LAN (WLAN); they can be used to study and model user mobility behaviours, beneficial to a variety of applications such as traffic management, urban planning and location-based advertisements. 
One of the applications of mobility modelling consists of formulating predictive models to forecast individual trajectories, for which various methods have been proposed, including Markov chains <cit.>, neural networks <cit.>, finite automata <cit.> and Naive Bayes <cit.>. 
Existing research have used several datasets differing with respect to their spatial and temporal granularity, resulting in vastly contrasting prediction accuracies ranging from over 90% to under 40% <cit.>.


In this context, the seminal paper of Song et al. <cit.> laid the foundations for computing a theoretical upper bound on the maximum predictability of human mobility.
This work establishes a benchmark for quantifying the performance of different algorithms and generalises its approach across various datasets.    
It defines the maximum predictability, noted π^max, as the temporally correlated entropy of information of an individual's trajectory.
π^max is estimated by first computing the entropy based on the Lempel-Ziv data compression <cit.> and then by solving the limiting case for Fano's inequality <cit.>, an information-theoretical result used to compute lower bounds on the minimum error probability in multiple-hypotheses testing problems.
The proposed theoretical upper bound (π^max = 93%) is empirically validated using a CDR dataset collected over 50,000 users by a telecommunications  operator for a duration of three months.
It should be noted however that CDRs are a rather rough approximation of  human mobility.
Song et al. <cit.> also show that π^max is independent of radius of gyration and movement periodicity, hence observe an insignificant level of variation across a heterogeneous population.
  

Several subsequent works either redefine this theoretical upper bound or perform empirical validations with different mobility datasets.
Lu et al. <cit.> re-estimate π^max to be 88% and use Markov chains to empirically verify this redefined upper bound. 
They analyse another CDR dataset consisting of 500,000 users, collected for a duration of five months and achieve an average predictability of 91% with an order-1 Markov chain.
They also show that higher-order Markov chain models do not improve prediction accuracy. 
Their interpretation behind surpassing their own estimated theoretical bound is that trajectories exceeding this bound are non-stationary, whereas the accuracy of stationary trajectories prevails within the bound. 
A trajectory is considered to be stationary when people tend to remain still during short time-spans.
This conclusion directly contradicts findings of Song et al. <cit.>, because non-stationary trajectories should by definition have a higher entropy.
Additionally, Cuttone et al. <cit.> show that the stationary nature of trajectories plays a significant role in the higher accuracies resulting from Markov models <cit.> as they often predict the user will remain in the previous location, i.e., self-transitions.
Lin et al. <cit.> also show that π^max is independent of the data sampling rate which was later questioned by Smith et al. <cit.> (π^max = 81 ± 4% ) and Cuttone et al. <cit.> (π^max = 65%).
Smith et al. <cit.> and Cuttone et al. <cit.> use mobility datasets <cit.> containing GPS trajectories and empirically show that predictability has a direct correlation with the temporal resolution and an inverse correlation with the spatial resolution.


The CDR datasets used in these works are known to have inherent gaps due to the short bursts of calls masking the user's true entropy.
Human mobility varies under time translations, therefore the entropy not only depends on the duration of past observations but also on number of visited locations, these factors tend to be hidden in such datasets <cit.>. 
Additional inconsistencies become evident due to the fact that the authors in <cit.> group the user locations into one hour bins when constructing the historical trajectory of a user.
Further inspection suggests that these models can foresee future locations at π^max, only when an individual is present in one of the top n bins <cit.>.
The first two works <cit.> thus consider the last location of each day, consequently predicting only the user's home place.
Under such a scenario, Ikanovic et al. <cit.> and Cuttone et al. <cit.> show that the predictability of the true next location is significantly lower (π^max=71.1 ± 4.7%) than the predictability of the location in the subsequent bin.
They further show that an individual's mobility entropy is directly proportional to the number of visited locations. 
The authors also point out that the generating function behind the stochastic mobility behaviour is often unknown.
Therefore the bounds cannot be estimated theoretically and require empirical derivation. 
Cuttone et al. <cit.> achieve an even lower bound on π^max of 65% on the same datasets with the same methods as Ikanovic et al. <cit.>.


Zhao et al. <cit.> argue that lower entropy and higher predictability in mobility behaviour originates from its non-Markovian character.
To highlight the non-Markovian character of human mobility, they analyse the rank distribution of the visited locations and the associated dwelling times.
Using a CDR dataset, they show the presence of the scaling law (Zipf's) behaviour in one-point statistics <cit.>. 
However, human mobility is better described by two-point statistics as it involves implicit behavioural dynamics while traveling from a given location to another with a specific intent <cit.>.


Based on the research literature discussed above, we observe a discrepancy regarding the maximum predictability bound, π^max and disagreements on the impact of entropy, the number of uniquely visited locations and the spatiotemporal resolution of the trajectory on π^max.   
In order to gain a deeper understanding about this discrepancy, we construct next-place prediction models by using seven different approaches.
We compute the empirical maximum accuracy and compare it with the theoretically derived π^max, considering three large-scale real-world datasets containing GPS trajectories. 
We find that recurrent-neural architectures <cit.> significantly surpass π^max on datasets compiled for long timespans. 



The current approaches <cit.> estimate π^max by first computing the true entropy of user mobility, denoted by S^real using Lempel-ziv compression <cit.>, followed by the computation of the minimum error probability by leveraging Fano's inequality <cit.>.
Fano's inequality is based on the assumption that the system is governed by a Markovian process and computes the conditional entropy of a random variable X, relative to the correlated variable Y, to the probability of incorrectly estimating X from Y, thus yielding the minimum error probability, noted p_e.
In practice, p_e is computed by segmenting the entire trajectory into sub-strings, where the length of the shortest substring beginning at index i does not appear previously <cit.>.
The total predictability π^max is thus the complement of the average of error probabilities on individual substrings. 



We show that when a trajectory is split into substrings, the entropy associated with the individual substring increases, thus increasing S^real.
This occurs as Fano's inequality, rooted in information theory <cit.>, is intended for a data source with known probability distribution <cit.>.
Human mobility prediction however, is based on the discretization of the trajectories, where the probability distribution is not known a priori.
Furthermore, the estimation of entropy by using Lempel-Ziv coding <cit.> was originally constructed to provide a complexity measure for finite sequences, i.e., input sequence displaying exponential decay in long-range correlations (memoryless structure). 
Such sequences, when further split, changes the true distribution of the data and increase the associated entropy; and the derived π^max thus acts as a limit on the Markov model.
We argue that this is due to ignoring the presence of long-range structural correlations that are present in human trajectories.



 
The current π^max computation is thus based on the widely used assumption that human mobility is Markovian (memoryless), i.e., the movements are independently distributed, as testified by the numerous research works relying on Markov models to characterize mobility <cit.>.
In opposition to this widely used assumption, we show the presence of non-Markovian character in human mobility dynamics by empirically showing that the drop in the mutual information <cit.> follows a power law function of the distance between any two points in a trajectory.
We also show the state-dependent nature of mobility by illustrating the presence of scaling laws in the distribution of location dwell times and the rank associated with location visits. 
Using real-datasets, we show that these underlying mechanisms that govern human mobility are the same across disparate mobility behaviours.




§ BACKGROUND





 §.§ Complexity and Criticality


Self-organized complexity is a phenomenon which combines self-organization and criticality to describe complexity <cit.>.
It is a property of dynamic systems to regulate their microscopic behavior to be spatial/temporal scale independent <cit.>.
This resembles their behavior at a critical point of phase transition.  
However, unlike conventional phase transitions, self-organized complexity governed systems do not depend on any external tuning of control parameters, i.e., the evolution of dynamic complex systems is self organized into the critical behavior <cit.>.



The critical point of a statistical system is defined as the point in space, parameterized by intensive quantities (for example, temperature and pressure), at which there exist no boundaries between the phase transitions.
Simply put, it is the end of the co-existence curves of two phases <cit.>.
On the exterior of this critical point, the different phases are cramped due to a correlation in the phase (or other quantifiable properties) of adjacent elements. 
Such a correlation prevails at come distinctive length scale (the correlation length).
This can be also visualized as the density of the boundary between the phases. 
When a dynamic system arrives at the critical point, this correlation length becomes infinite, resulting in fuzzing the boundaries. 
Such an infinite correlation length does not contain any information regarding a finite scale length.
As one cannot use an infinite length as a unit of measurement, the theory describing the system behavior at the critical point is scale-invariant.
Therefore, it is clear that critical behavior implies scale invariance, and therefore when critical behavior is involved, the effects at distances much longer than microscopic lengths are imperative to study.




In general, a dynamic system self organizes into a complex state but with a fairly general structure. 
The complexity occurs due to an absence of a single well defined characteristic event size. 
However, despite the complexity, the system exhibits some statistical properties governed by power laws <cit.>.
For instance, the number of events, E as a function of its size s (wherein a major event is less likely to occur as compared to a small one) can be described as in Equation <ref>.



    E(s) = A× s^-α


Where A is some constant and α describes the statistical features of a self-organized critical (SOC) state.
In general SOC has been observed in slowly-driven and non-equilibrium systems which posses extended degrees of freedom and high level of non-linearity <cit.>.
In this paper, we show that human mobility patterns follow the above characteristics.





 §.§ Criticality in Practice


We give some examples to connect such a physical phenomenon to instances in the practical world. 
There has not been a plethora of work to characterize practical systems as SOC.
Amongst the limited work, we look at how Carreras et al. <cit.> analyze a 15-year long time series of electric power transmission system blackouts.
The authors present a possible quantitative explanation of the complex dynamics in a power systems which contains properties which leads to dynamic equilibrium with some properties of SOC.
Another work by Bartolozzi et al. <cit.> consider the stock markets which are complex self-interacting systems, characterized by intermittent behaviors. 
Here periods of high activity alternate with periods of relative calm. 
On the basis on empirical evaluations the authors conclude that stock markets depict a near-SOC state. 


In physical systems, criticality occurs when there exists correlation length-scale or time-scale that diverges to infinity. 
Therefore, the relationships in the sequence are not local in time, but are arbitrarily non-local.




 §.§ Exponential Decay and Power-law




In SOC dynamic systems, the probability of occurrence of major disruptive events decreases as a power function of the event size <cit.>. 
This is in contrast to many conventional systems in which this probability decays exponentially with event size.
The difficulty arises as the correlation length is infinite.
However, it is well known and established that in a self-organized critical system (SOC), the long range interactions occur due to power-law decay of correlation.


The power law is detected and characterized, by combining maximum-likelihood fitting methods with goodness-of-fit tests based on Kolmogorov-Smirnov statistic and likelihood ratios <cit.>. 

At this point, we summarize the difference between exponential decay and power law.
Although visually both might look similar as both are positive and go asymptotically to 0, the difference is pointed out in Equation <ref> and Equation <ref>.


    Power   Law: y = x^k



    Exponential   Decay: y = k^x


Where k is some constant.
Thus, it is evident that exponential decay goes to zero much faster than power law.
Furthermore, exponential probability distribution has an inherent property of memoryless-ness, whereas power law follows a polynomial relationship exhibiting the property of scale invariance.  




 §.§ Entropy & Mutual Information





Formally, the entropy of a discrete random variable X with probability mass function (pmf) p_X(x) is:


    H(X) = -∑_x p(x)log p(x) = -E[log(p(x))]


In a nutshell, the entropy measures the expected uncertainty in X.
Therefore, H(X) is equal to the amount of information learnt on an average from one instance of the random variable X and p(x), the probability distribution.
Furthermore, the entropy does not depend on the value that the random variable takes, but only on the probability distribution p(x) <cit.>.

In case two random variables, X,Y jointly distributed according to the pmf p(x,y), the joint entropy is given by:


    H(X,Y) = -∑_x,yp(x,y)log p(x,y)


Furthermore, the conditional entropy of X given Y can be given by:


    H(X|Y) = -∑_x,ylog p(x|y) = -E[log(p(x|y))]


Simply put, the conditional entropy is a measure of how much uncertainty remains about a random variable X given Y. 

The Mutual Information I, between two discrete random variables X,Y jointly distributed according to pmf p(x,y) is given by:


    I(X;Y)    = ∑_x,yp(x,y)logp(x,y)/p(x).p(y)
        = H(X) - H(X|Y) 
        = H(Y) - H(Y|X) 
        = H(X) + H(Y) - H(X,Y)



Mutual information can be analogously defined for continuous variables and conditional mutual information.
Simply, mutual information can be explained at follows.
Consider two people Alice and Bob, living in Switzerland. 
Alice goes out in a t-shirt if the temperature is more than 20 degrees celsius.
Bob always wears a t-shirt irrespective of the temperature.
Now, notice that Alice's actions give information about the temperature in Switzerland.
However, Bob's actions give no information whatsoever. 
This is due to Alice's actions are random and correlated with the temperature, whereas Bob's actions are deterministic. 
Mutual information I mathematically quantifies this notion.




 §.§ Markov Processes



Memoryless Markov Models. A Markov process is defined by a matrix M containing conditional probabilities M_ab such that M_ab = P(X_t+1 = a | X_t = b).
The Markov matrices (or stochastic matrices) are bounded by M_ab≥ 0 and ∑_aM_ab = 1.
Therefore the complete dynamic of the model can be specified as provided in equation below:


    p_t+1 = M.p_t


where p_t is a vector with components P(X_t = a) that specifies the probability distribution at time t.
If λ_i denotes the eigenvalues of M, sorted in decreasing magnitude such that: |λ_1| ≥ |λ_2| ≥ |λ_3|.
The Markov matrices have |λ_i| ≤ 1, with corresponding eigenvector resulting in a stationary probability distribution.

Furthermore, two additional conditions are imposed on Markov matrices: i) M is irreducible, i.e. every state is accessible from every other state, and ii) M is aperiodic, to avoid processes that will never converge (for example a → b → a → b). Thus it has been shown that M being irreducible and aperiodic implies that the individual eigenvalues are less than one and the stationary probability distribution is unique <cit.>.



Hidden Markov Models (HMM). In an HMM, in addition to the observed sequences X_1, X_2,...,X_n, there are hidden latent variables Y_1, Y_2...,Y_n that form a Markov chain as well.
These internal hidden dynamics are never observed, but at each time-step, an output is produced such that Y_i → X_i form the observed sequence.



Markov Model with Memory. Markov processes are able to model and represent higher-order dependencies between successive observations of a state variable.
One way to increase the process memory is by using high-order Markov models.
In such a model, the conditional probability distribution of future states in the process not only depends on the current state but also on the past states.
A Markov chain with memory m is a process satisfying,


    P(X_t+1 = s_t+1|X_t = s_t,...,X_1=s_1) = P(X_t+1 = s_t+1|X_t=s_t,...,X_t-m+1 = s_t-m+1)


for all t ≥ m =. By defining


    Y_t = (X_t, X_t-1,...,X_t-m+1)


and by taking the ordered m-tuples of X values as its product space so that the chain Y_t with suitable starting values satisfies the Markov property.



§ RESULTS


In this section, we first present the values of S^real and π^max computed using the approach mentioned in the works of Song et al. <cit.> and Lu et al. <cit.>.
We then discuss the accuracy results estimated using seven algorithms and compare them with respect to the theoretical upper bound.
To investigate the presence of memory in human mobility, we conduct several experiments and use the results to illustrate the reason for surpassing the upper bound. 
Finally, we investigate the existing approach and discuss their failure to compute the true entropy (S^real) of mobility trajectories.



Mobility datasets. We conduct the experiments by using three mobility datasets whose specifications, along with the estimated values of S^real and π^max, are shown in Table <ref>. 
The PrivaMov dataset <cit.> was collected through GPS, WiFi and GSM in the city of Lyon (France) and includes university students, staff and their family  members.
The Nokia mobile dataset <cit.> (NMDC) was collected in the Lake Geneva region of Switzerland and consists young individuals' trajectories, collected through GPS, WLAN, GSM and Bluetooth.
The GeoLife dataset <cit.> was collected in Beijing (China) and contains trajectories recorded through GPS loggers and GPS-phones.
The computation of S^real and π^max at the aggregate level for the dataset is based on our observation of independence of predictability on travel distance (radius of gyration r_g) in human mobility, which is consistent with previous studies <cit.>.






Prediction Algorithms. We estimate the empirical predictability using seven different approaches: (1) Markov chains <cit.> (order 1-5), (2) Hidden Markov model <cit.>(HMM), (3) Vanilla Recurrent Neural Network <cit.> (Vanilla-RNN), (4) Recurrent Neural Network with Long Short-Term Memory <cit.> (RNN-LSTM), (5) Dilated Recurrent Neural Network <cit.> (Dilated-RNN), (6) Recurrent Highway Network <cit.> (RHN), and (7) Pointer Sentinel Mixture Model <cit.> (PSMM).
We find that higher order Markov chains (typically >3) do not contribute to increase prediction accuracy, as also observed by Lu et al. <cit.>.
The prediction accuracy for Markov chain models and the recurrent-neural architectures for all datasets is shown in Figure <ref> and Figure <ref>, respectively. 
 







We observe that the accuracy of Markov models lie in the vicinity of the theoretical π^max. 
It is also clearly evident that recurrent-neural architectures significantly outperform Markov models with respect to their average accuracies.
Recurrent-neural architectures are a class of artificial neural networks, which use their hidden memory representation to process input sequences. 
The variants of these architectures differ in their capacity to manipulate this memory and propagate gradients along the network.
For instance, RHN's are built to account for short and long-range correlations present in a sequence, which explains their superior performance as compared to the other architectures. 
Whereas, PSMM's weigh long-range dependencies much higher than short-distance correlations in the sequence.
In Table <ref>, we show the maximum predictability achieved by using the best performing models from each algorithm, and in Figure <ref> we compare their performance with the theoretical upper bound.







The prediction accuracies of recurrent-neural architectures also surpass the theoretical upper bound for the respective dataset.  
This anomaly in computing π^max is puzzling, even more so considering the diversity of the datasets with respect to their collective time spans, visited number of locations, demographics and spatiotemporal granularity.
This lets us question the assumption that human mobility follows a Markov process. 
Here, we conduct extensive analysis to empirically prove its non-Markovian nature.




Location Rank Distribution. In order to gain insight into the datasets, we first analyse the rank distribution of the locations, according to the visit frequency at individual and aggregated levels.
An individual visits different locations depending on a perceived priority attached to the location <cit.>; this results in a heterogenous location frequency distribution <cit.>.
To study the location-rank distribution, we follow the approach stated in Zhao et al. <cit.> in order to rank locations according to their collective magnitude at the aggregate level. 
Figure <ref> shows the existence of power-law scaling  (Zipf's law) in the rank distribution of visited locations in human mobility.
We also observe a convergence and robustness at the individual level, which clearly indicates non-uniform mobility behaviour and its effect on entropy, hinting at the non-Markovian nature of human mobility <cit.>.






Dwell Time Distribution. To further confirm this non-Markovian nature, we check the distribution of the dwell-times associated with the individual locations.  
The current mobility models are based on an assumption that human movements are randomly distributed in space and time, hence are approximated by a Poisson process <cit.>.
However, Barabasi <cit.> shows that human activities are non-Poissonian, by showing that inter-event timings depict long-tailed distribution.
We observe a similar behaviour when considering  human mobility in all the datasets, when examining the dwell-times associated with each location; most locations are visited at high periodicity, while few locations encounter long waiting times.
The current models assume that inter-event time follows exponential distribution <cit.>, rather, we observe an emergence of power-law as shown in Figure 5, 6, 7.
The spikes in the plot correspond to delays and display the visit regularity, which indicates a long-tailed process.
The delay-time distribution depicts the priority list model in human mobility, bearing similarity to other activities as remarked by Barabasi <cit.>.
When an individual is presented with multiple events under the context of mobility, the next location is determined on a perceived priority, thus resulting in power-law dynamics in inter-location waiting times <cit.>.
This shows that the dwell-times associated with human mobility are not memoryless, hence cannot be considered as Markovian. 
In the above analysis, we also observe a convergence between individual mobility patterns and aggregated datasets, which concurs with the observations of Yan et al. <cit.>. 

 










Criticality and Mutual Information. Criticality is the property of dynamic systems to regulate their microscopic behavior to be spatiotemporally scale-independent <cit.>.
As a result, critical behaviour implies scale invariance; and when criticality is involved, the effects at distances much longer than microscopic lengths are crucial to study.
In practice, criticality is observed when there exists a correlation time-scale that diverges to infinity, thus the relationships in the sequence are arbitrarily non-local.
 

Mutual information I quantifies the measure of information provided by a symbol/location coordinate (Y) about another symbol/location coordinate (X).
In case of mobility, mutual information between two location instances is the realisation of a discrete stochastic process, with separation τ in time <cit.>.
Lin et at. <cit.> express I, as a function of the number of symbols (locations) between any two symbols and state that it would decay with a power-law for any context-free grammar and hence must be non-Markovian.
In order to perform this validation on human mobility, we first estimate mutual information as a function of distance for the GeoLife dataset <cit.>.
This choice is based on uniformly sampled location points in the dataset.
To validate the emergence of power-law at distinct sampling rates, we undersample and oversampled (spatial semivariance interpolation <cit.>) the dataset by a factor of two and four and show the trend in Figure <ref>.




As suggested by Lin et al. <cit.>, we observe a power-law decay at all the sampling rates.
Contrary to our assumption that I would increase and decrease by the factor of under/over sampling, we observe a decrease in I for all the contexts in which the true distribution of the data is altered.
We also observe that the reduction is proportional to the Kullback-Leibler divergence <cit.> between their respective distributions.
The reduction in I stems from the fact that a change in the distribution results in the alteration of the true correlation between the location pairs.
The true distribution will therefore show maximum I, compared to the cases when either artificial pairs are introduced (oversampling) or true pairs are removed (undersampling) from the dataset.
To verify our hypothesis, we also calculate the joint entropy for all the cases and observe an increase in H(X,Y) for the altered distributions as shown in Figure <ref>. 
We see that the increased entropy is due to an increase in the ratio between unique pairs in the dataset over the total number of pairs.
The introduction of spurious pairs scrambles the true distribution as it leads to introduction of data points in the true sequence, thereby changing the random variables sampled at distance D, hence reducing I.
This occurrence was confirmed after computing the area under the curve (ROC), which was maximum for the true data distribution in the first quartile as compared to the rest as shown in Figure <ref>. 
This explains our observation of higher joint entropy for the oversampled and the undersampled case. 












To analyse the long-range correlations present in each of the datasets, we compute their respective mutual information decay.
This information will serve as basis for the difference in accuracy for each dataset and the performance difference between the prediction algorithms.
Again, we observe a power-law decay across all the datasets and their respective joint entropy, as shown in Figure <ref> and Figure <ref>.
We further explore the Markov transition matrices for these datasets and observe that they are reducible and periodic, resulting in the decay of I to a constant.
It has been shown that such a characteristic of the transition matrix cannot result in an exponential decay by Lin et al. <cit.>.  
They show that an irreducible and aperiodic Markov process, with non-degenerate eigenvalues, cannot produce critical behaviour because I decays exponentially. 
This phenomenon is seen in a number of cases, including hidden and semi-Markov models <cit.>.
In the literature, such behaviour is superficially dealt with by increasing the state space to include symbols from the past, which does not address the main issue <cit.> with Markov models; lack of memory.
This analysis shows that GeoLife dataset consists of considerably higher number of long-range correlations, compared to the PrivaMov dataset and the NMDC dataset.
This should be self-evident from their respective data collection durations. 
However, the lower dependencies in the NMDC dataset, compared to PrivaMov, is due to the smaller area of the data collection region, which generally results in lower entropy of movement <cit.>.




Here, we reason about the accuracy variation within and between the datasets and about the performance differences between the prediction algorithms. 
We observe that the NMDC dataset provides higher accuracy as compared to the other datasets, and witness a lower variation within the accuracies of different algorithms. 
This stems from the presence of very short dependencies in the individual trajectories present in the dataset, as seen in Figure <ref>.
The lower correlations also result in roughly equivalent prediction accuracies within the predictive models. 
The lower accuracies of recurrent-neural architectures, compared to Markov chain at some time-steps are due to the models tendency to actively seek for long-range dependencies.
However, if the dataset does not contain such dependencies, the model underperforms, unless it is weighted to account for such an existence. 
This underperformance is evident from the behaviour of dilated-RNN's, where an increase in dilations (to account for longer dependencies) results in dropping accuracy.
Such a phenomenon has also been observed in language modelling tasks, which suggests that this is not a domain specific occurrence <cit.>.
The performance drop in the recurrent-neural architectures at different time steps is due to capturing the long-distance dependencies to different degrees, resulting in either under/over fitting.
An additional reason for higher accuracy in NMDC dataset is due to a lower number of unique locations and smaller variations in the dwell-times, as compared to the PrivaMov and GeoLife datasets, as shown in Figure <ref> and Figure <ref>.
These aspects directly correlate with the entropy and thus affect predictability <cit.>.
We also observe that PSMMs, perform better on GeoLife dataset, compared to other two, due to its ability to search for dependencies at longer distances.




Entropy and Predictability Estimation. The current method <cit.> uses Lempel-Ziv data compression scheme <cit.> to compute the mobility entropy. 
This approach segments the complete trajectory into substrings, where a substring is defined as the shortest length element subsequence yet to be encountered.
As observed by Lesne et al. <cit.>, a vast majority of substrings are of length one or two, which are the dominant contributors to the entropy.
The estimated entropy is thus the outcome of finite-size fluctuations; and the total count of the substrings and of the elements in a substring does not represent the true probability distribution. 
Furthermore, in this process the structural correlation between the individual substrings is ignored, based on the argument that the probability of joint occurrences is weak <cit.>.
This argument stems from the reasoning that the parsed substrings are independently and identically distributed (iid) according to Gaussian distribution, that does not apply to mobility trajectories.
Finally, the correlated features can be compressed only by memorising all the cases of intervening random variables between the correlated instances. <cit.>.
It has thus been proved that Lempel-Ziv approach fails to capture redundancies in the data sources with long-range correlations <cit.>. 



More importantly, Storer et al. <cit.> shows that standard data compression approaches, such as the Lempel-Ziv approach cannot truly capture long-range dependencies, as the information carriers of a sequence lie in its structural origin. 
However, these approaches limit the entropy estimation process at the sub-string level.
Given that entropy is the complete quantitative measure of the dependency relations (including many point correlations), the computation of higher-order entropy is non-trivial. 
Therefore, it is flawed to assume that the π^max derived from such an approximate estimation of S^real should act as an upper bound of predictability on trajectories compiled for long time-spans.  





§ DISCUSSION



The previous research <cit.> estimated S^real and pi^max by using CDR datasets spanning a period of three to five months.
Such datasets do not truly capture features such as the total number of unique locations visited by an individual, due to its low granularity (typically 4-5km <cit.>).
This results in a dataset with a masked entropy and mobility patterns ignoring long-range correlations. 
An important point to note is that for very short distances, power-law decay and exponential decay may not be trivial to differentiate <cit.>.
This was due in part to the unavailability of high granularity GPS datasets and the fact that that previous works <cit.> were only studied for short distances of human mobility. 
Therefore, the assumptions underlying the computation of S^real and π^max would have been fairly easy to overlook.




The aforementioned inadequacies would reinforce the empirical validation of π^max using Markov chains, however, as mentioned above this would result in an error-prone estimation of the predictability. 
As seen in other domains of sequential-data modelling such as natural language processing, it has been shown that Markov chains are fundamentally unsuitable for modelling such processes <cit.>.
Our empirical observations, backed by theoretical foundations, indicate that human mobility will be poorly approximated by Markov chains.
This is particularly afflictive for trajectories that satisfy criteria such as long time-span of collection and large radius of gyration of movement.    




Our choice to rely on mutual information was based on its triviality and its domain independence.
As shown by Lin et al. <cit.>, the mutual information decay offers some insights into why recurrent-neural architectures exceed probabilistic models in terms of capturing criticality. 
As for a Markov process, the observations at t_n depends only on events at previous time step t_n-1 or on previous n time-steps for an n-order Markov chain.
Under such a context, the maximum possible predictive information is given by the entropy of the distribution of states at one time step, which is in turn bounded by the logarithm of the number of accessible states.
Unlike Markov chains, the recurrent-neural architectures, such as RHN's, approach this bound while maintaining the memory long enough, that the predictive information is reduced by the entropy of transition probabilities.




We expect to provide a more sophisticated description of the underlying phenomenon as more of the trajectory is observed. 
Consequently, increasing the number of parameters in the model. 
That is, when we examine trajectories on the scale of individual coordinates, we learn about the rules of combining these points into points of interest and the transition paths between them.
At the next level, if we consider several of these points of interest and the paths, we learn the rules for combining these points into semantic patterns.
Similarly, when we look at semantic patterns, we learn about the visitation periodicities and circadian rhythms associated with the mobility behaviours.
Therefore, longer traces have increasing number of long range structural correlations that are non-trivial to be captured by the currently available entropy measure.
Moreover, the current approximation implies that the substrings have the same compressibility factor <cit.>, hence the results derived from this approach would coincide with the average.   
Thus, the current computation will result in higher estimates of entropy, consequently resulting in a lower predictability bound. 




Even though Markov models tend to underperform in modelling human mobility, their use in human mobility prediction is not entirely without precedent.   
In fact, considering their low computational complexity, it might be advantageous to opt for a Markov model when a dataset contains short-distance dependencies and low number of unique locations.
However, in datasets exhibiting criticality, long-range correlations appear in the vicinity of the critical point, which necessitate recurrent-neural architectures to accurately model human mobility. 
In this paper, we have shown that human mobility exhibits scale-invariant long-range correlations which can be quantitatively measured by a power-law decay of mutual information.
We highlight that the exponent characterising the power-law decay of the correlations is well defined for infinite sequences. 
For mobility trajectories, however, the accuracy of the analysis is restricted by the length of the substrings and their entropy, which results in an incorrect estimation of maximum predictability.
This explains why the empirical results surpass the theoretical upper bound in several previous research works and in our own experiments.







§ METHODS



Mobility prediction. We define mobility prediction as forecasting the transitions between places, after eliminating all self-transitions <cit.>.
A preliminary step in achieving this consists of transforming the raw GPS locations into a sequences of points of interest <cit.>.
A point of interest is defined as any location where an individual visits with an intentional purpose with a perceived priority for e.g., home/work place, gym, train station etc. 
Among plethora of existing works dedicated to the problem of extracting these points, we rely on our approach that is independent of a priori assumptions regarding the data and individual mobility behaviors <cit.>.
We then convert the raw GPS trajectory of a user u, T_u = ⟨(lat_1,lon_1,t_1), (lat_2,lon_2,t_2)...(lat_n,lon_n,t_n) ⟩, where lat_i, lon_i are the latitude and longitude coordinates respectively and t_i is the timestamp such that t_i+1 > t_i into a sequence of temporally ordered points of interest, s(t) = ⟨ (poi_1,t_1), (poi_2 t_2)...(poi_n,t_n)⟩, where poi_i is the point of interest at index i.
The mobility prediction task is thus formulated as: given a sequence s(t) up to a timestamp n, predict the next point of interest at timestamp n+1. 
The prediction accuracy is then estimated by following the approach stated by Lu et al. <cit.>.

  

Predictive models. We use the standard implementations of the predictive algorithms as described in their respective papers. 
Markov chains <cit.> and hidden Markov models <cit.> are implemented using the standard python libraries (hmmlearn).
We use hyper-parameters stated in these papers.
Vanilla-RNN <cit.>, RNN-LSTM <cit.> and dilated-RNN <cit.> are based on predicting the next character (language modeling) in the text, whereas RHN <cit.> and PSMM <cit.> model the prediction task as multivariate classification. 
For dilated-RNN <cit.> we use the dilations of 1, 2, 4, 8, 16, 32 and 64 and provided the results for dilation 32 after which we observe a drop in the accuracy.




Computing S^real and π^max. These values are computed as stated in Song et al. <cit.> as per Equation <ref> which is based on Lempel-Ziv data compression <cit.>.


    S^real = (1/n∑_i=1^nλ_i )^-1ln(n)



Where, n is the length of the trajectory (total number of locations) and λ is defined as the length of the shortest substring at an index i not appearing previously from index 1 to i-1.
Note that we use the same base (2) in entropy estimation as for the logarithm in Fano's inequality.
Furthermore, their length is set to zero upon reaching index i, when no more unique substrings can be computed using the above method. 
π^max is then estimated by solving the limiting case of Fano's inequality <cit.>.



Mutual information. Equation <ref> gives mutual information, I, between two discrete random variables X,Y jointly distributed according to the probability mass function p(x,y).
Mutual information I thus quantifies the number of bits of information provided by a symbol (Y) about another symbol (X).



    I(X,Y) = ∑_x,yp(x,y)log(p(x,y)/p(x).p(y))
 

Mutual Information can also be estimated by computing the entropy of the marginal distribution of discrete random variables X and Y, and the joint entropy of discrete random variables X and Y.


    I(X,Y) = H(X) + H(Y) - H(X,Y) = D_KL(p(XY)||p(X)p(Y))


where H(X) is the entropy of a random variable X and H(X, Y) is the joint entropy of X and Y. 
D_KL is the Kullback-Liebler Divergence <cit.>. 
Thus, Mutual Information is same as the Kullback-Leibler Divergence between distributions of X and Y.
In order to compensate for insufficient samplings, we use the following adjustment proposed by Grassberger et al. <cit.> to compute H(X), H(Y), H(X,Y).

    H(X) = log N - 1/N ∑_i=1^k Niψ((Ni)



In order to compute the mutual information in a mobility trajectory, we first estimate the distribution of a random variable from index 0 followed by the distribution of the second random variable in the data at index D. 
The random variables X and Y are sampled from the individual trajectory sequence. 
X denotes sequence starting at index 0 and Y is sequences starting at offset D from 0. 
D is then varied to compute long-distance dependencies at every separation by creating displacements between the random variables. 
Once the contextual dependence limit is reached, this approach starts sampling noise, which sets the termination criterion.
Finally, the average similarity between the two symbols is quantified.




§ DATA AVAILABILITY


The GeoLife dataset that supports the findings of this study is public and is made available by Microsoft Asia <cit.>. 
The PrivaMov dataset is collected by Universite de Lyon and can be obtained by submitting an online form <cit.>. 
The Nokia Mobile Data Challenge (NMDC) can be obtained in a similar fashion <cit.>.
We provide links in the citations to avail these datasets.
Furthermore, our source codes are made public <cit.> and further clarifications will be provided upon request.



abbrv




