Modeling Human Spatial Mobility Patterns with the L\'evy Flight Cluster Model

TL;DR

This paper introduces the Lévý Flight Cluster Model, a hierarchical Bayesian approach for analyzing and simulating human mobility patterns from irregular location data, with applications in understanding overlaps and anonymization.

Contribution

The paper presents a novel hierarchical Bayesian mixture model for human mobility analysis, capable of capturing movement characteristics and generating synthetic location data.

Findings

Accurately models human movement patterns

Effectively identifies activity overlaps

Generates realistic synthetic location data

Abstract

Despite the extensive collection of individual mobility data over the past decade, fueled by the widespread use of GPS-enabled personal devices, the existing statistical literature on estimating human spatial mobility patterns from temporally irregular location data remains limited. In this paper, we introduce the L\'{e}vy Flight Cluster Model (LFCM), a hierarchical Bayesian mixture model designed to analyze an individual's activity distribution. The LFCM can be utilized to determine probabilistic overlaps between individuals' activity patterns and serves as an anonymization tool to generate synthetic location data. We present our methodology using real-world human location data, demonstrating its ability to accurately capture the key characteristics of human movement.

Figures

Tables

Table 1. Table 1. Prior specification for the LFCM .

Parameter	Prior Distribution	Hyper-prior Distribution
$N_G$	Poiss($\lambda$)	Gamma($\frac{1}{2},\frac{1}{2}$)
c	Multin($\underset{¯}{\omega }$)	$\underset{¯}{\omega }\sim$ Dir($\underset{¯}{1}$)
$b$	Bern($\nu$)	$\nu \sim$ Beta$(2,2)$
${\Vert \underset{¯}{\Delta x}\Vert }_2|b=0$	$𝒩({\underset{¯}{\mu }}_g\Delta t,{𝛀}_g^{-1}\Delta t$)	${\underset{¯}{\mu }}_g,{𝛀}_g\sim NW(\underset{¯}{0},.01,\frac{1}{2},\frac{1}{2}𝐈)$
$z$	Multin($\underset{¯}{\xi }$)	$\underset{¯}{\xi }\sim \text{Dir}\left(\underset{¯}{1}\right)$
$\eta$	Bern($𝓅$)	$𝓅\sim$ Beta($2,2$)
${\Vert \underset{¯}{\Delta x}\Vert }_2|b=1$	Par(${\alpha }_g$) ($\epsilon$ known)	${\alpha }_g\sim$ Gamma$(\frac{1}{2},\frac{1}{2})$
${\tan }^{-1}\left(\frac{\Delta x_2}{\Delta x_1}\right)$	VM($\vartheta$, $\tau$)	$\vartheta \sim$ VMises(0, 1)

Table 2. Table 2. Average jump length, MSD, and ROG for observed data, simulated data using LFCM, and simulated data using the CPT estimator. Standard deviations are in parentheses.

	Jump Length	MSD	Radius of Gyration
Observed	1.556	2.889	1.250
	(0.052)	(0.041)	(0.0040)
Simulated (LFCM)	1.580	2.975	1.190
	(0.035)	(0.056)	(0.038)
Simulated (CPT), W=0.2	1.120	1.136	0.902
	(0.045)	(0.052)	(0.034)
Simulated (CPT), W=1.0	1.305	1.187	0.935
	(0.073)	(0.096)	(0.038)
Simulated (CPT), W=1.5	1.599	1.589	1.140
	(0.040)	(0.093)	(0.027)

Equations

$$\Delta r(t_i)$$

$$f_{\Delta r}(\Delta r(t))$$

$$\text{MSD}(t_n)$$

$$⟨r^2(t)⟩\propto t,$$

$$\sqrt{⟨r^2(t)⟩}\propto t^{1/2}.$$

$$r_g(t_n)$$

$$f_{r_g}(r_g(t))$$

$$f_{\Delta r}(\Delta r(t))={\int }_{r_0}^{\infty }{f}_{\Delta r}\left(\Delta r(t)|r_g(t)\right)f_{r_g}\left(r_g(t)\right)d{r}_g(t),$$

$$\underline{X}(t)$$

$$f_{\underline{X}(t)}(\underline{X}(t))$$

$$f(\lambda x)$$

$$B(t+h)-B(t)\perp \perp B(s),$$

$$\frac{\partial \rho (x,t)}{\partial t}$$

$$\underline{\mu }(t)\equiv \underline{a}+\frac{t-t_{i-1}}{t_i-t_{i-1}}(\underline{b}-\underline{a})$$

$$X(t_i)-X(t_{i-1})\perp \perp X(t_j)-X(t_{j-1}),$$

$$\underset{h\to 0}{\lim }P(|X(t+h)-X(t)|>\epsilon )=0.$$

$$1-F_r(r;\epsilon ,\alpha )$$

$$p(j\to k)=p_{d(j,k)\mid d(j,h)}\prod _{l\le d(j,k)}a(k_l),$$

$$J(t)$$

$$X(t)$$

$$J_1(t)+J_2(t)+\cdots J_n(t)\stackrel{d}{=}{a}_{n}J(t)+b_n,$$

$$X(t)$$

$$\equiv \mu t+{\sigma }^{\ast }B(t)+J(t).$$

$$\Delta X(t)$$

$$f_{\Delta x}(\Delta X(t)|\mu ,{\sigma }^{\ast },\vartheta ,\alpha ,\Delta t)$$

$$f_{\Delta x}(\Delta X(t)\mid N_G,\underline{\theta },\Delta t)$$

$$\displaystyle\phantom{\equiv\sum_{g=1}^{N_{G}}\omega_{g}\biggl{[}}+\nu f_{J^{*}}(\Delta r(t);\alpha_{g},\Delta t)\biggr{]},$$

$$f_{\underline{\Delta x}}\left(\underline{\Delta X}(t)\mid N_G,\underline{\theta },\Delta t\right)$$

$$\displaystyle\phantom{\equiv\sum\omega_{g}\ \biggl{[}}+\nu f_{J^{*}}(\underline{\Delta r}(t);\vartheta_{g},\alpha_{g},\Delta t)\bigr{]},$$

$$1-F_{\Delta r/\Delta t}(\Delta r/\Delta t)$$

Full text

KEYWORDS: Activity space; global positioning systems (GPS); human mobility; space-time geography; stochastic processes

CONTACT M. Wolff. Email: mlw32@uw.edu

August 30, 2025

1. 1 Introduction

2. 2 Mobile device GPS data

3. 3 Human mobility metrics

4. 4 Human mobility models

5. 5 The Lévy flight cluster model

6. 5.1 Memory Dependence

7. 5.2 The Lévy Flight Cluster Model and its Bayesian specification

8. 6 Simulation study

9. 7 Analysis of mobile device GPS data

10. 8 Discussion

11. A Collapsed MCMC sampler for the LFCM

12. B Derivations of posterior distributions

13. C Similarity metrics for activity areas

1. Introduction

The collection of human activity data from the early twentieth century was based on travel journals (Wolf et al., 2001; Schönfelder and Axhausen, 2003, 2004), surveys (Buliung et al., 2008), and social application queries (Luo et al., 2016). Advancements in technology now include high-resolution GPS data (Hato, 2006; Stopher et al., 2008; Bohte and Maat, 2009), with smartphones offering a wealth of information that enables deeper analysis of human activity patterns (e.g., Kwan, 2012; Matthews and Yang, 2013; Dobra et al., 2015; Williams et al., 2015). For example, agent-based modeling (ABM) frameworks simulate autonomous actions to analyze behavioral shifts within larger patterns. This can reveal how individuals face geographical mobility constraints such as access to transportation (Kenyon et al., 2002) and epidemic-related stay-at-home orders (Perry et al., 2021), linked to social exclusion and changes in employment opportunities (Goldhaber et al., 2015).

However, accurately estimating spatiotemporal trajectories from individuals’ GPS data presents its own set of challenges. Human mobility displays a combination of both stochastic and decidedly deterministic movement patterns; localized activities like being at home, in an office, or at a gym, as well as traveling, lead to distinctive behavior patterns, each suggesting different degrees of social overlap— probability of social interaction is high within a shared office space but practically nonexistent on a heavily trafficked highway. These behavioral patterns can also alter dramatically due to unobserved factors such as changes in financial status or the occurrence of a natural disaster. Since GPS data collected from mobile devices are primarily application dependent, the data consist of temporally irregular samples of individuals’ time-stamped geolocations, which are closely related to their current activities. This introduces further complexity in inferring where and how the unrecorded time is spent, thus presenting an added challenge to the process of activity space estimation.

Early models of individual activity often relied on assumptions of multivariate normality, resulting in broad elliptical regions of probabilistic activity (Sherman et al., 2005; Newsome et al., 1998). Other proposed models defined activity regions as the convex hull around a specific set of points (e.g. Buliung and Kanaroglou, 2006; Fan and Khattak, 2008), which often unrealistically described the true activity area and encompassed large areas of uncertainty. More recently, methods such as kernel density estimation (e.g., Chen and Dobra, 2020; Dong et al., 2020) have been employed, which assume uniform distributions of activity within arbitrarily defined grid cells—where the estimation of activity overlap can be significantly influenced by the chosen cell size, as we show in Section 6.

A separate body of research prominent in animal telemetry attempts to recover the latent position process that underlies the observed locations of individuals (see Hooten et al. (2017) for an extensive review). By modeling the latent position process, they addressed the temporal dependence present in the observed data. Typically, these models rely on GPS data and employ Gaussian interpolation to deduce a latent position based on observed points and elapsed time (Hooten et al., 2017). Over time, these models have seen various improvements, such as the inclusion of factors such as velocity, areas of attraction, and heterogeneous sampling times. For example, Scharf et al. (2017) introduce a hierarchical Bayesian multiple-imputation framework to fill gaps in sparse telemetry data, treating the true continuous path as a latent variable and demonstrating improved inference under measurement error. One notable model is the Brownian bridge movement model (BBMM) (Horne et al., 2007; Kranstauber et al., 2012, 2014), which determines individual latent positions using a sequence of Brownian bridges. However, BBMMs have been found to be inadequate in capturing the nuances of human mobility (Krumm, 2021); Brownian bridges do not allow “loitering” behavior or large movements for, for example, nonpedestrian travel, nor do they yield realistic inference when locations are detected not by GPS, but by a spatially static network of motion sensors that result in a highly irregular sampling rate that also varies between individuals (Ahmadian et al., 2025). Parallel to these advances, Scharf et al. (2015) use discrete regime-switching indicators to infer dynamic social networks from movement data, illustrating the power of latent indicators for regime changes.

Characterizations of human mobility using non-Gaussian random-walk models have also been explored. For example, studies by Gonzalez et al. (2008) and Rhee et al. (2011) have highlighted that human mobility often shows heavy-tailed travel distances similar to a Lévy flight process. Song et al. (2010) employed a modified Lévy flight process to model mobility, thus encapsulating geographic return patterns that are not adequately represented by random movement. In a similar vein, Alessandretti et al. (2017) utilized logarithmic Gaussian mixtures in a discrete movement model to mirror the heavy-tailed nature of human mobility. Although these models represent characteristics of the latent process underlying observed individual mobility, they fail to provide uncertainty bounds for individuals’ regions of highest activity based on local and global movement patterns. This gap underlines the remaining need for models capable of capturing the complex dynamics of human mobility while providing a realistic notion of uncertainty.

To address these gaps, we introduce the Lévy flight cluster model (LFCM) — a generative Bayesian model for human mobility trajectories. Using the Lévy Itô decomposition, we model individuals’ movement as a Lévy flight exploration process that may exhibit a mixture of repetitive “activity clusters” represented by local Brownian motions. Our model can be estimated from temporally irregularly sampled location data and with an arbitrarily small or large number of data points. Unlike BBMMs, this model can represent non-movement (periods of time spent at the same location) as an individual’s “activity cluster”, non-pedestrian movement as a Lévy jump, and estimates propensity for an individual to return to each identified cluster. Samples from the joint posterior distribution of the model’s parameters can be used as a generative procedure for individuals’ activity distributions, which are highly comparable to the individuals’ original mobility patterns as measured by common mobility metrics. This provides a novel approach to human mobility to data anonymization. By simulating the spatial mobility patterns of individuals with the LFCM, researchers can fit accurate models of human mobility based on the resulting synthetic data while preserving the privacy of the individuals involved. Moreover, our generative model can be used an interpolative and extrapolative model, which allows one to use common metrics to compare distances between individuals even when individuals are not observed at the same times. This means that, as we show in Section 7, one may generate probabilistic sociomatrices across any number of individuals, which can be used for later statistical inference. Finally, our model can be efficiently estimated using analytical integration and Gibbs sampling methods and is completely parallelizable across individuals and trajectories. For these reasons, LFCM contributes significantly to understanding human mobility, providing not only a nuanced model that captures the complex dynamics of human mobility but also establishing a foundation for further research to estimate activity areas and their overlap.

The structure of the paper is as follows. In Section 2 we describe the characteristics of contemporary mobile device GPS data, specifically referring to GPS-generated geolocations collected at nonuniform time intervals from 293 anonymized devices within a metropolitan region for a duration of 12 weeks. In Section 3, we examine commonly used accuracy metrics related to the study of individual mobility. These include jump length, radius of gyration, mean-squared displacement, and measures of frequented locations. In Section 4 we review previous research on individual mobility models. In Section 5 we introduce the LFCM together with a Bayesian framework for inference and estimation of this model. In Sections 6 and 7 we evaluate the performance of the proposed model with simulated and real-world human location data. This assessment focuses on the model’s ability to generate individual-level probabilistic activity areas, representative mobility sample paths for each individual, and accurate representation through the accuracy metrics mentioned above. To demonstrate an immediate practical application of the proposed model, in Section 7 we illustrate the estimation of the probabilistic activity overlap between individuals. This investigation of the nuances of human mobility and its application in the definition of activity spaces serves to advance our understanding in the field and provides a solid foundation for future research on the relationships between human mobility and exogenous environmental factors.

2. Mobile device GPS data

Mobile devices and their applications have become ubiquitous in the lives of people living in metropolitan areas. These devices and apps gather a wealth of data including timestamps, geographic location, communication content, and application-specific details, often under the purview of the relevant privacy policies. Data provisioning companies collect device data from mobile application developers or their aggregators. The types of data they accumulate include mobile ad identifiers (like Apple iOS IDFAs or Google Android IDs), time-stamped geographic coordinates of devices (latitude and longitude), the precision of these coordinates, and the device’s direction and speed.

In this article, we use the location records provided by a data provisioning company for 293 mobile devices active in a large metropolitan area. These data span a 12-week period from November 2018 to January 2019. Figure 1 displays the distribution of records by day of the week. Apart from the dates surrounding the winter holiday season, the records are relatively uniformly distributed across the period.

Each record includes a unique device identifier, geographic coordinates, a timestamp, an estimated margin of error for the coordinates, and a reported speed of travel. This allows identification of daily movement patterns for individual devices over the 12-week period. Figure 2 presents the recorded trajectories of three different devices over a single week, starting November 12, 2018. It is evident from these trajectories that they contain both random and patterned components. They also highlight the intermittent temporal gaps in the recorded locations. Such gaps are common throughout this data set.

Numerous mobile device applications carry out location-sensitive data collection while in operation. This leads to a significant fluctuation in timestamps by hour of the day and day of the year. Figure 3 presents the average number of records by hour of the day for the same period (displayed at the bottom). Given this irregular rate of sampling, a temporally flexible method of estimating uncertainty around an individual’s observable activity patterns becomes essential. The unique attributes of mobile GPS data, including its high resolution and irregular sampling rate, present challenges that require robust and flexible analytical methods to fully capitalize on potential of the data to capture human mobility trajectories.

3. Human mobility metrics

Quantification of individuals’ daily living trajectories often uses four recognized metrics: jump length (Brockmann et al., 2006; Gonzalez et al., 2008), radius of gyration (Gonzalez et al., 2008), mean squared displacement (Alessandretti et al., 2017), and frequent locations (Song et al., 2010). These metrics, applied within defined time periods, can capture unique aspects of human mobility. Many studies find that human travel often follows a heavy-tailed distribution —- with many short trips and few long trips; jump length helps quantify this heterogeneity in travel behavior. The radius of gyration, as a measure of the typical distance an individual travels from their center of mass, is valuable for distinguishing between localized movement and more expansive travel patterns, indicating how confined or widespread an individual’s activity space is. The mean squared displacement serves a similar purpose but measures displacement from the process center, not a single origin point. The number of frequented locations, assessed through a frequency rank, captures the regularity and predictability in routines and can help identify meaningful activity centers. Each metric serves as a reference point for the validity of the movement model in the context of individual mobility. Details for each of these metrics are given below.

Jump Lengths

Let $\{\underline{X}(t)\in\mathbb{R}^{2}:t\geq 0\}$ denote a spatio-temporal process. The jump length of this process at time $t_{i}$ is the Euclidean distance between the locations $\underline{x}(t_{i-1})$ and $\underline{x}(t_{i})$, denoted as

[TABLE]

The jump length distribution $f_{\Delta r}(\Delta r(t))$ describes the likelihood that members of a population travel specific distances and is critical to accurately modeling population mobility. Brockmann et al. (2006) estimated the distribution of jump lengths for trajectories of U.S. bank notes according to the probability density function

[TABLE]

for the 95% confidence interval $\alpha\in(1.57,1.61)$. This “power law” distribution has been consistently observed in further research on human mobility. Gonzalez et al. (2008) modeled jump lengths between individual GPS location data using a truncated version of this power law,

[TABLE]

finding a 95% confidence interval for $\alpha\in(1.50,1.90)$. Using two separate data sets on the location of GPS of individuals and the mode of transportation, Zhao et al. (2015) found that the aggregate distribution of $\Delta r(t)$ followed a power law distribution with $\alpha=1.55$ and $\alpha=1.39$, respectively.

Figure 4 shows a histogram of jump length distributions for 293 individuals in the GPS mobility data described in Section 2 alongside the density function $f_{\Delta r}(\Delta r(t))\propto[\Delta r(t)]^{-1.59}$, where $1.59$ is the MLE for $\alpha$ found in Brockmann et al. (2006), showing that the observed jump lengths roughly match a power-law distribution.

However, this comparison assumes that jump lengths across individual/weeks are identically distributed. We also compute jump lengths between consecutive GPS locations within each (individual, week) group and analyze the upper tail of these lengths with a continuous power law model. For each group we select $x_{\min}$ by minimizing the Kolmogorov–Smirnov distance, then estimate the tail index $\alpha_{g}$ by maximum likelihood with asymptotic standard error $\mathrm{SE}(\hat{\alpha}_{g})\approx(\hat{\alpha}_{g}-1)/\sqrt{n_{g,\text{tail}}}$. To obtain a single overall effect while accounting for between–group heterogeneity, we pool the $\hat{\alpha}_{g}$ via random–effects meta–analysis (REML), weighting by $w_{g}=1/(\mathrm{SE}(\hat{\alpha}_{g})^{2}+\tau^{2})$, and find a 95% CI for the pooled $\hat{\alpha}=(1.71,1.74)$ with $\tau^{2}=0.075)$, suggesting jump lengths slightly above estimates in prior work.

Numerous papers have further corroborated the power law behavior of human mobility (e.g., Song et al., 2010; Rhee et al., 2011; Alessandretti et al., 2017). This behavior is evident at the individual level in our GPS mobility data as illustrated in Figure 5 that shows the empirical complementary cumulative distribution functions (ECCDF) on the log-log scale. Each of the three plots represents a different individual device and each curve represents the trajectory over a week of observation. The approximate linearity of the tail in these distributions suggests heavy-tailed, or approximate power-law behavior.

Mean Squared Displacement (MSD)

Characteristically, MSD quantifies the average squared distance an entity moves from its initial position up to a certain point in time; this can help infer the nature of the underlying movement process, e.g. whether it is diffusive or superdiffusive.

Definition 3.1** (Mean Squared Displacement).**

Let $\underline{x}(t_{0})$ denote the initial point observed in the spatio-temporal process $\{\underline{X}(t):t\geq 0\}$. The Mean-squared displacement of the process up to time $t_{n}$ is defined by

[TABLE]

the average squared distance from the starting location $\underline{x}(t_{0})$ up to time $t_{n}$.

As jump lengths are drawn from a heavy-tailed distribution, such as a power law or a truncated power law (see Equation (3)), the process initially exhibits large fluctuations with occasional long jumps. However, truncation that is caused, for example, by travel restrictions limits the maximum jump length, ensuring that the distribution has a finite variance. Thus the central limit theorem guarantees that the sum of many independent jumps converges to a Gaussian distribution. Consequently, while short-time dynamics may reflect heavy-tailed behavior, long-term aggregated motion will asymptotically resemble a Brownian motion.

For the Brownian motion, the mean squared displacement (MSD) grows linearly with time:

[TABLE]

where the operator $\langle\cdot\rangle$ is often used to represent expectation in the statistical mechanics literature, which implies that the typical displacement (i.e., the square root of the MSD) scales as

[TABLE]

Thus, despite the heavy-tailed nature of individual mobility events, the truncation forces the MSD to eventually exhibit the diffusive scaling characteristic of Brownian motion.

Radius of Gyration.

The heavy-tailed behavior of humans on certain scales suggests that they tend to move a characteristic distance away from their starting locations. This distance is often quantified by the radius of gyration, $r_{g}(t)$, defined as the root mean square distance of a set of points from a given axis. Traditional applications of this measure are rooted in physics (where it is related to the mass moment of inertia) and engineering (distribution of the cross-sectional area in a column).

Definition 3.2** (Radius of Gyration).**

The Radius of Gyration $r_{g}(t)$ is given by:

[TABLE]

where $\underline{x}(t_{i})$ are the $n$ geospatial coordinates of, e.g., an individual device trajectory, and $\overline{\underline{x}(t)}$ is mean of the trajectories $n$ points.

When applied to human mobility, the radius of gyration is often used to characterize the typical distance of an individual from the center of mass of their trajectory. Gonzalez et al. (2008) used mobile phone cell records to determine the radius of gyration for two sets of users; one set every two hours over one week and one set on a 6-month scale. The radius of gyration $r_{g}(t)$ was calculated for each user, and it was found that the distribution of $r_{g}(t)$ follows a truncated power law of the form

[TABLE]

where $r_{0}$ represents a minimum radius of the gyration cut-off point due to the spatial sensitivity of the data, and $\kappa_{r}$ represents the upper cut-off due to the finite size of the study area. The measured 95% confidence interval of the scaling exponent $\alpha_{r}\in(1.40,1.80)$—estimated from observations across both datasets—indicates a significant degree of heterogeneity in the travel habits of the observed population. Measuring the jump length distribution conditional on the radius of gyration $r_{g}$, $f_{\Delta r}(\Delta r(t)\ |r_{g}(t))$ revealed that users with small $r_{g}$ travel mostly over small distances, whereas those with large $r_{g}$ tend to display a combination of many small and a few larger jump sizes.

It is important to note that the distributions described above are not independent but instead are related by the equality

[TABLE]

where $f_{\Delta r}(\Delta r(t)|r_{g}(t))$ is the jump-length distribution introduced previously. If $f_{\Delta r}(\Delta r(t))$ has a power law scaling exponent $\alpha$, then we have $\alpha=\alpha_{c}+\alpha_{r}$ for conditional scaling exponent $\alpha_{c}$. This form of scaling corresponds to a type of random walk called a Lévy flight, and the results suggest that this may be the behavior of individuals up to their associated characteristic distance $r_{g}$.

Frequented Locations

A characteristic quality of human mobility is the frequent use of particular spaces or activity regions. Many people tend to travel to a specific region or along specific routes for work, return home on a daily basis, and frequent the same region to buy food, utilities, or clothing. In the literature on human mobility, it is common to quantify these frequent locations through the use of a classification system (Chen and Dobra, 2020), and it has been shown that the probability of an individual returning to a particular location is inversely proportional to its frequent rank (Gonzalez et al., 2008). The tendency to return to known areas has also been quantified in the literature through the number of new locations visited over time, where GPS mobility data have been found to grow proportionally to $t^{0.6}$ (Song et al., 2010).

4. Human mobility models

We give an overview of the statistical frameworks used in the development of human mobility models. These models are treated as stochastic processes, subject to variability due to local movement and more extensive habitual travel. They are frequently modeled as structured random walks, a method that borrows concepts derived from estimation of animal mobility (e.g. Liu et al., 2018; Meekan et al., 2017). These models also serve as the foundation for state-of-the-art mobility models. By examining these approaches, we can better understand the dynamics underlying individual-level mobility patterns.

Random Walks

The concept of a random walk, first put forth by Karl Pearson in 1905 (Pearson, 1905), is generally used to describe a sequence of random steps within a certain context, such as fluctuations in a gambler’s wealth, the movement of a molecule, variations in a stock price, or the foraging behavior of an animal. A well-known model of a random walk is one based on a regular lattice, in which movement occurs on a discrete grid of locations according to a predetermined probability distribution. This lattice-based model can be expanded into a discrete random walk by considering movements across more complex spaces. Continuous-time random walks take this concept a step further by incorporating movements made within a continuous space over random time periods. A basic illustration of a continuous-time random walk within a two-dimensional real coordinate space ($\mathbb{R}^{2}$) can be represented by the following stochastic process

[TABLE]

where $\underline{\Delta X}(t_{i})\sim f_{\underline{\Delta X}}(\underline{\Delta X}(t))$ are the first differences and $n(t)\sim f_{n|t}(n|t)$ is a random number of jumps in the interval $(t_{0},t)$. Then

[TABLE]

One property that plays an important role in the characterization of random-walk distributions is that of scale invariance, wherein any magnification shows self-similarity to the whole. A function $f(x)$ is considered scale invariant under all rescalings if for any scalar $\lambda$,

[TABLE]

for some exponent delta $\delta$. This property holds with $\delta=0$ for white noise (Gardiner, 2009; Van Kampen, 2007), and $\delta=-2$ for Brownian motion (Mandelbrot and Van Ness, 1968; Feller, 1971; Barabási and Stanley, 1995), for example.

The property of scale invariance has been empirically seen in human mobility data through the use of power-law distributions to model jump lengths (e.g., Gonzalez et al., 2008). However, random walks alone do not adequately describe numerous other aspects of human mobility such as a nonuniform distribution of jump angles and nonlinear growth of newly visited locations.

Brownian Motion

Brownian motion or the Wiener process is a stochastic process that is often used to model random movement in continuous time and space (Doob and Doob, 1953). Although the original observations of Brownian motion were rooted in physical phenomena, its mathematical representation has broad applications, such as financial modeling or physics.

Definition 4.1** (Brownian motion).**

A Brownian motion process $B(t)$ is characterized by the following properties:

• 1.

$B(0)=0$.

• 2.

Independent increments: for every $t>0$,

[TABLE]

for $u,s\geq 0$ and $s\leq t$.

• 3.

Gaussian increments: $B(t+h)-B(t)\sim\mathcal{N}(0,h)$.

• 4.

$B(t)$* is continuous in $t$.*

This process is also often described from a statistical mechanics perspective. For the probability density $\rho(x,t)$ of the particles at position $x$ and time $t$, Brownian motion satisfies the diffusion equation:

[TABLE]

where $D$ is the diffusion coefficient, a positive constant that quantifies the rate at which the probability mass spreads over time. If $n$ particles start at the origin at $t=0$, Eq. (10) for $\sigma^{2}=2Dt$ has the solution $\rho(x,t)=\frac{n}{\sqrt{2\pi\sigma^{2}}}\exp(-\frac{1}{2\sigma^{2}}x^{2})$, which is the density function of the normal distribution, suggesting that the mean square displacement of the Brownian motion is $\langle x^{2}\rangle=2Dt$. The form of the mean squared displacement in relation to time is often used to classify diffusion processes. For example, when the relationship between time and displacement is non-linear, $\langle x^{2}\rangle\propto t^{\gamma}$, it is called an anomalous diffusion process, said to exhibit super-diffusion if $\gamma>1$ and sub-diffusion when $\gamma<1$.

Brownian motion has a substantial history in describing movement processes in ecology. A seminal paper by Horne et al. (2007) introduced the Brownian bridge movement model (BBMM) to study individual animal movement and regional occupation times. The model infers gaps in the observed process $\{\underline{X}(t_{i}):i=1,\ldots,n\}$ with two-dimensional Brownian bridges. That is, conditional on $\underline{x}(t_{i-1})=\underline{a}$ and $\underline{x}(t_{i})=\underline{b}$, $\underline{X}(t)$ for $t\in(t_{i-1},t_{i})$ is distributed as $\mathcal{N}\left(\underline{\mu}(t),\sigma^{2}(t)\mathbf{I}\right)$ where

[TABLE]

The BBMM deals with serial correlation and unequal time intervals by directly incorporating this temporal independence.

Kranstauber et al. (2012) extended the BBMM to incorporate heterogeneous activity in animal movement, calling the model the dynamic Brownian bridge movement model (dBBMM). The authors augment the original BBMM by assuming local dispersion parameters on subsets of the observed process calculated by sliding windows. By introducing the additional parameters window size $w$ and break point $b$, they are able to identify sudden changes in the animal’s movement behavior. In separate work, Kranstauber et al. (2014) develop another extension called the bivariate Gaussian bridge (BGB) model to allow variation in the dispersion parameter in two dimensions, one parallel and the other orthogonal to the angle between the start and end points. Together, these models are widely used in the literature on animal movement (e.g., Kays et al., 2015; Allen and Singh, 2016; Rickbeil et al., 2019).

However, recent work by Krumm (2021) has suggested that BBMMs do not accurately capture human mobility patterns. Considering dispersion parameters estimated by triplets of points no more than 50 kilometers apart and 48 hours apart, the author finds a large variation in diffusion coefficients specified at the group, individual, and triplet level, suggesting ambiguity in the level at which the human mobility model should be considered, and finds that intermittent points rarely satisfy the assumption of normality after applying the model. The author concludes that the dBBMM and BGB models may address some of these concerns; however, the author leaves analysis of these models to future work.

Lévy processes and flights

Brownian motion processes are a subset of a larger class of stochastic processes with stationary independent increments called Lévy processes (Bertoin, 1996).

Definition 4.2** (Lévy process).**

A stochastic process $X(t)$ is said to be a Lévy process if it satisfies

• 1.

$X(0)=0$* almost surely.*

• 2.

Independent increments: for any $0\leq t_{1}\leq t_{2}\leq\dots t_{n}<\infty$,

[TABLE]

when $i\neq j$.

• 3.

Stationary increments: For $t>s$, $X(t)-X(s)\stackrel{{\scriptstyle d}}{{=}}X(t-s)$.

• 4.

Continuity in probability: For any $\varepsilon>0$ and $t\geq 0$,

[TABLE]

Several key papers have found that several human mobility metrics show numerous statistical similarities with a particular type of Lévy process called a Lévy flight. The term “Lévy flight” was originally used by Benoît Mendelbrot (Cannon, 1984) as a description of a heavy-tailed stochastic process $r(t)$ with step sizes defined by the survival function

[TABLE]

commonly known as the Pareto distribution (Arnold, 2015). This distribution, which follows a power law according to its parameter $\alpha$, is widely used to model jump lengths in the human mobility literature.

Although traditional Lévy flights are memoryless, human mobility exhibits repetitive returns to specific areas. Recent modeling techniques of human mobility attempt to explicitly incorporate this. Song et al. (2010) develop an exploration and preferential return (EPR) model to match several intuitive aspects of human activity, modeling that the probability of visiting new locations decreases with time and a preference to return to a location with a probabilistic decision stage before moving: exploration with probability $P_{new}=\rho S^{-\gamma}$, or return with probability $P_{ret}=1-\rho S^{-\gamma}$. The conditional distance covered during an “exploration” is then modeled according to a heavy-tailed distribution $P(\Delta r)$, while during a “return” according to a distribution proportional to the number of previous visits. This alteration allows preservation of properties of Lévy flight while agreeing on the propensity to visit new locations and the number of returns to particular locations.

Alessandretti et al. (2017) offer an explanation for the apparent disconnect between the nonrandom intention of individuals and the scale-free nature of their movement, citing that mixtures of normal or log-normal distributions with different variances can generate a power law tail of movement often associated with scale-free distributions (Gheorghiu and Coppens, 2004). The authors model coordinate space as a nested structure of containers which humans traverse and movement as a two-step process in which the individual decides first which level of nested containers to move across and then which container to move in to. That is, the probability of transition from location $j$ to location $k$ is

[TABLE]

where $p_{d(j,k)|d(j,h)}$ represents the probability of traveling at level distance $d(j,k)$ given a current level distance $d(j,h)$ from a home location $h$, and $\prod_{l\leq d(j,k)}a(k_{l})$ is the probability of choosing a specific location $k$ at level distances $l\leq d(j,k)$. The size of each container, the maximum distance between any two points within a container, is described by a lognormal distribution, as is the time spent within each container. The authors find that this model is highly accurate in replicating certain characteristics of human mobility, such as the distribution of displacements, the radius of gyration, and visitation frequencies.

5. The Lévy flight cluster model

Motivated by the literature gap for a generative model for the estimation of individual activity, we develop the Lévy flight cluster model (LFCM), a novel statistical framework to estimate individual activity distributions from a collection of few temporally irregular GPS location samples. The joint posterior distribution of the LFCM parameters serve as a generative representation of individual activity distributions, which not only provides a novel approach to data anonymization, enabling researchers to simulate human mobility patterns while preserving the privacy of study participants, but also supports both interpolative and extrapolative analyses. As a result, comparisons of mobility between individuals are feasible even when data are collected at different times, allowing the creation of probabilistic sociomatrices for subsequent statistical inference. The remainder of this section develops the model progressively: we first lay the stochastic foundation via the Lévy-Itô decomposition and a jump-diffusion approximation; we then remove temporal dependence through a first-difference representation; next we introduce activity heterogeneity and memory via a finite mixture with returns; and we conclude with the Bayesian specification and priors for inference.

The LFCM leverages the Lévy-Itô decomposition theorem to model the latent stochastic process associated with an individual’s trajectory as a memory-aware mixture of Brownian motions and pure jump processes. Through the use of conjugate priors (Ryan et al., 2017) almost all nuisance parameters can be integrated out, allowing the development of efficient posterior sampling methods.

Lévy-Itô Decomposition

A fundamental result in the study of stochastic processes attributed to Paul Lévy (Lévy, 1938) and Kiyosi Itô (Itô, 1941) is the Lévy-Itô decomposition theorem, identifying a Lévy process $\{X(t):t\geq 0\}$ as the independent sum of a drifted Brownian motion process $\mu t+\sigma B(t)$, a compound Poisson process $J(t)$, and a square integrable martingale $\varsigma_{\varepsilon}(t)$ with almost surely countably many jumps with magnitude less than some specified $\varepsilon$. The compound Poisson process $J(t)$ is defined by

[TABLE]

where $N(t)$ is a Poisson process and $\mathfrak{J}_{i}$ are independent random variables with distribution function $F_{\mathfrak{J}}$. Hence, the Lévy process $X(t)$ can be written as

[TABLE]

If $\mathfrak{J}_{i}$ follows a symmetric distribution with power-law tails, decreasing as $|\mathfrak{J}|^{-\alpha-1}$ for $0<\alpha\leq 2$, the sum $J(t)$ will tend to a stable distribution (Gnedenko et al., 1954); that is, for every $n\in N$, there exist constants $a_{n}>0$ and $b_{n}\in\mathbb{R}$ such that

[TABLE]

where $J_{1}(t),J_{2}(t),\ldots,J_{n}(t)$ are i.i.d. copies of $J(t)$, and $\overset{d}{=}$ denotes equality in distribution.

Jump Approximation

For stable Lévy processes of index $\alpha\in(0,2]$, the martingale term $\varsigma(\varepsilon)$ is well approximated by Brownian motion in the sense that there exists $\sigma_{\varepsilon}$ such that the process $\sigma_{\varepsilon}^{-1}(\varsigma_{\varepsilon}(t)-\mathbb{E}[\varsigma_{\varepsilon}(t)])\stackrel{{\scriptstyle d}}{{\rightarrow}}B(t)$ as $\varepsilon\searrow 0$ (Asmussen and Rosiński, 2001). The Lévy process in Eq. (13) then follows the approximate equality

[TABLE]

Figure 6 illustrates the approximation using three instances of random walks in $\mathbb{R}^{2}$, each run for 7,000 steps. Brownian motion steps are shown in black and Lévy flight jumps are shown in red; the left process is a Brownian motion, the middle process is a Lévy flight, and the right process is a Lévy flight with jumps less than $\varepsilon$ replaced by a Brownian Motion step. While the Brownian motion and Lévy flight processes are distinct, with the latter exhibiting longer individual movements and the appearance of clustering, the right process in Figure 6 has similar characteristics to a Lévy flight.

Temporal Independence

The temporal dependence in Eq. (14) can be largely removed by considering the first difference representation

[TABLE]

since the components $B(\Delta t)\sim\mathcal{N}(\underline{0},\Delta t\mathbf{I})$ are independent. Here, $\Delta X(t)$ can represent a change in geospatial location, $\Delta t$ a change in time, and the jump diffusion process $J^{*}(\Delta t)$ is either zero (i.e., $N(\Delta t)=0$) or of magnitude strictly larger than some specified $\varepsilon$. Hence, Eq. (15) can be written as

[TABLE]

Letting $\nu=\nu_{\Delta t}\equiv P(N(\Delta t)>0)$, where we assume that the $N(\Delta t)$ are identically distributed, the process (16) has density

[TABLE]

for a pure jump process distribution $f_{J^{*}}(\Delta r;\alpha)$ with jump length parameter $\alpha$.

The above assumption highlights a key difference between Lévy flights and Lévy walks. In a Lévy flight, an individual (or object) makes instantaneous jumps between locations. These jumps can be of various lengths, but the probability of making a longer jump decreases as a power-law. In a Lévy walk, the individual (or object) still follows the power-law distribution for step lengths, but the motion within each step is not instantaneous. Instead, there is a defined speed, and longer steps will thus take more time to complete than shorter ones. Due to the sparsity of the data, and a lack of information on what the individual is doing, we make the assumption that the process is a Lévy flight, rather than a Lévy walk.

Activity Heterogeneity

To account for heterogeneity in an individual’s activities, the mixture formulation of $X(t)$ in Eq. (17) can be further extended to a finite mixture of $N_{G}$ processes:

[TABLE]

for mixture weights $\omega_{g}$ (e.g., the probability of a point belonging to activity $g$), where each $g^{th}$ mixture component is characterized by the tuple $(\mu_{g},\sigma_{g}^{*},\alpha_{g})$, and $\underline{\theta}$ is the vectorized collection of model parameters for parsimony. Eq. (18) allows the model to fit individuals with both a relatively low number of activities (e.g., living a sedentary lifestyle) and those with more active lifestyles. Eq. (18) is easily adapted to a multivariate framework described by

[TABLE]

where $N_{G}$ is the number of activity groups, $\underline{\mu}_{g}$ is a group-specific drift, $\mathbf{\Sigma}_{g}^{*}$ is the group-specific covariance matrix, and $f_{\Delta J^{*}}(\underline{\Delta r}(t);\vartheta_{g},\alpha_{g},\Delta t)$ is the distribution function of the pure jump process with jump angle parameter $\vartheta_{g}$ and jump length parameter $\alpha_{g}$. We assume the jump length and jump angle are independent, i.e. of the form $f_{J^{*}}(\underline{\Delta r}(t);\vartheta_{g},\alpha_{g},\Delta t)=\prod_{g=1}^{N_{G}}f_{\theta}(\theta;\vartheta_{g})f_{\Delta r}(\Delta r;\alpha_{g},\Delta t)$.

Pure Jump Size Distribution

Due to the complexity of aggregating a Poisson distributed number of arbitrary i.i.d. random variables, the compound Poisson process $J(t)$ is difficult to model in practice. Instead, we use a data-driven simplification of the density $f_{\Delta r/\Delta t}(\Delta r/\Delta t)$. A Lévy flight is characterized by the power-law tail of its step length, and the step lengths of human mobility as represented by GPS device location data observe a power-law tail distribution; Figure 7 depicts complementary cumulative distribution functions (CCDF) for three distinct individuals observed in the data, each of which presents a linear tail on the log scale. This follows because

[TABLE]

We model jump lengths as a Pareto distribution with common parameters $\alpha,\varepsilon$ between activity groups, that is,

[TABLE]

In specifying the Pareto distribution as a conjugate prior, we must provide an estimate of $\varepsilon$. Since the estimation of the parameter $\varepsilon$ involves a domain search, traditional MLE methods for joint estimation do not apply. Naive grid searches of the entire algorithm may be prohibitively expensive. Although one could still specify a prior distribution on $\varepsilon$, this removes the benefit of integrating the closed-form expression for the Pareto posterior, speeding up computational efficiency.

Empirical Estimation of $\varepsilon$

We consider three separate methods of empirical estimation for $\varepsilon$—using the Kolmogorov-Smirnov (KS) statistic, the Anderson-Darling (AD) statistic, and the Kuiper’s (K) statistic. The KS statistic is the maximum absolute difference between the empirical cumulative distribution function (ECDF) of a random variable $\varepsilon$, $\widehat{F}_{n}(\varepsilon)$, and a corresponding fitted parameterized distribution $\widehat{F}_{\varepsilon}(\varepsilon;\theta)$:

[TABLE]

As an alternative to the KS statistic, the AD statistic is defined as

[TABLE]

The AD statistic tends to be more sensitive to the tails of the distribution (Engmann and Cousineau, 2011). The KS statistic takes the pointwise maximum absolute difference of the ECDF and a fitted parameterized CDF, where this difference is unlikely to be largest at the tails. However, very small values of $\widehat{F}_{\varepsilon}(x_{i})$ have a large weight in the AD statistic.

Kuiper’s statistic (K) considers both positive and negative differences between the distributions. This makes the statistic equally sensitive to the tails as at the median and also makes it invariant under cyclic transformations of the independent variable. The statistic is defined by

[TABLE]

Other formulations have been used, such as the Bayesian information criterion (BIC), but many of these do not perform well in low sample regimes (Clauset et al., 2009).

Next we propose a new methodology for estimating $\varepsilon$ for small samples from power law generative distributions— a neighborhood averaging procedure that improves estimates in low sample sizes without degrading estimates for larger samples. Let $\mathcal{E}$ be a range of potential values of $\varepsilon$. For a ball $B_{r}({\varepsilon})$ of radius $r$ around ${\varepsilon}$, let the $k$-smoothed neighborhood of $\varepsilon$ be defined as

[TABLE]

with corresponding estimates $\widehat{\varepsilon}^{\mathcal{N}}=\tfrac{1}{|\mathcal{E}_{\mathcal{N}}|}\Sigma_{\varepsilon\in\mathcal{E}_{\mathcal{N}}}\ \varepsilon$ and $\widehat{\alpha}_{\text{MLE}}^{\mathcal{N}}=\tfrac{1}{|\mathcal{E}_{\mathcal{N}}|}\Sigma_{\varepsilon\in\mathcal{E}_{\mathcal{N}}}\ \widehat{\alpha}_{\text{MLE}}(\varepsilon)$.

To show the implication of estimation bias when computing $\widehat{\varepsilon}^{\mathcal{N}}$, Figure 8 illustrates the mean estimate of $\alpha$ for 5,000 samples drawn from the Pareto distribution, with $\alpha=2.5$, $\varepsilon=100$, and n = 2,500, plotted as a function of $\varepsilon$ for power law (purple) and truncated power law (green) distributions.

To show differences in the methods for computing $\widehat{\varepsilon}^{\mathcal{N}}$, we performed the following simulation study. For each $N\in\{20,30,50,100,500\}$, we draw $N$ samples from the distribution

[TABLE]

with $\rho(\cdot)$ denoting the density function for a Pareto distribution, $\alpha=2.5$ and $x_{\min}=1$. Then we estimate $\varepsilon$ and $\alpha_{\text{MLE}}$ using the KS statistic, the AD statistic, and Kuiper’s statistic.

The average values of $\varepsilon$ and $\alpha_{\text{MLE}}$ by sample size and metric are shown in Figure 9. We find that the Kuiper statistic outperforms the KS statistic in estimation of power law parameters for all sample sizes, with more significant differences with smaller samples. Using neighborhood averaging with the Kuiper statistic does at least as well, producing better parameter estimates at the smallest sample sizes.

Pure Jump Angle Distribution

A standard Lévy flight assumes a uniform angular distribution of movement, but in practice the angular distribution of movement must also be specified. Figure 10 shows radial histograms of jump angles for three example individuals over the 12-week period, evidencing non-uniform traversal directions within individuals.

Due to its favorable conjugacy properties (Guttorp and Lockhart, 1988), a convenient choice to model angular movement is the von Mises distribution with density

[TABLE]

where $\mu$ is the mean direction, $\kappa$ is a concentration parameter, and

[TABLE]

is the modified Bessel function of order zero (Abramowitz and Stegun, 1965). When $\kappa=0$ the von Mises distribution (21) becomes a uniform distribution in the unit circle and often serves as a reasonable and convenient prior. Guttorp and Lockhart (1988) showed that the conjugate prior of (21) is proportional to

[TABLE]

Given a set of observed angles $\underline{\vartheta}=(\vartheta_{1},\ldots,\vartheta_{n})$, the posterior distribution of $(\mu,\kappa)$ is

[TABLE]

where

[TABLE]

In what follows, we model the jump angle distribution $f_{\theta}(\theta;\vartheta)$ as a von Mises distribution with common parameter $\vartheta$.

5.1. Memory Dependence

Identifying clusters of repeated activity is an important aspect of understanding the long-term activity distribution of individuals. While the mixture model in Eq. (19) captures heterogeneous activities and distinguishes those activities from travel, it does not identify the regions in which the activities occur, nor it incorporates locations in the future likelihood of travel. In this section, we propose a probabilistic method for conditional identification of activity regions and returns to those regions.

Given a location $\mu_{0}$, the average Brownian movement over the time period $[0,T]$ is

[TABLE]

We assume that the distribution of $\underline{x}(t_{i})$ being a return to the activity region $z$ is proportional to $\mathcal{N}(\underline{\widetilde{\mu}}_{z},\mathbf{\Sigma}_{z})$, where $z$ represents a contiguous Brownian motion in the activity group $g$, and

[TABLE]

Figure 11 provides an illustrative example. Suppose that at some point in the past an individual was active in a region $z$ with center $\underline{\widetilde{\mu}}_{z}$ and covariance $\mathbf{\Sigma}_{z}$, denoted by the red sample path in the figure. At time $t_{i}$ an individual moves to $\underline{x}(t_{i})$ according to the jump distribution $f_{J^{*}}(\Delta r;\vartheta,\alpha)$. The conditional likelihood that the jump $\underline{x}(t_{i})$ is then classified as a return to activity $z$, i.e. $f_{z(t_{i})}(z(t_{i})|\underline{\Delta x}(t_{i}),\underline{\theta})$, is proportional to $\mathcal{N}(\widetilde{\underline{\mu}}_{z},\mathbf{\Sigma}_{z})$. Specifically,

[TABLE]

where $\mathcal{Z}_{i}$ is the set of all activities identified at time $t_{i}$.

This serves as a computationally convenient conditional classification procedure for capturing nuances of individual mobility. Returns to activity regions must be defined explicitly since, although individuals do exhibit nonrandom routines, GPS spatial data measurements are on a continuous scale and unlikely to be repeated at multiple time points. When viewed in continuous space, the number of visited locations can resemble a Lévy flight, but as a more intuitive notion of activity regions, the number of new locations visited decreases over time. This conditional identification procedure incorporates historical movement variability in activity areas during the estimation rather than relying on the assumption of uniform likelihoods within grid cells used in previous work (Chen and Dobra, 2020; Alessandretti et al., 2017). Once identified, these distributions, characterized by the parameters $(\widetilde{\underline{\mu}}_{z},\mathbf{\Sigma}_{z})$, can be used as known locations to simulate realistic probability-based return behavior observed in true mobility (e.g., using the 95% CI of the distribution as a location of activity, see Figure 17).

5.2. The Lévy Flight Cluster Model and its Bayesian specification

We are now ready to introduce our main modeling framework for human mobility patterns.

Definition 5.1** (Lévy Flight Cluster Model).**

For a stochastic process $\{\underline{X}(t):t\geq 0\}$, the Lévy Flight Cluster Model (LFCM) is given by the mixture

[TABLE]

where $\nu$ is the probability of a jump occurring, $N_{G}$ is the number of activity groups, $\omega_{g}$ is the probability that $\underline{\Delta x}(t)$ is part of activity $g$, $\mathfrak{p}$ is the probability of a return occurring, and $\xi_{z}$ is the conditional probability of returning to activity region $z$. Note that $N_{Z}$ is implicitly defined by the number of identified contiguous Brownian motion sample paths that belong to one of the $N_{G}$ activities.

We define the $n\times 1$ latent vectors $\underline{b},\underline{c},\underline{\eta}$, and $\underline{z}$: the assignment of observation $i$ to activity cluster $g$ is denoted by the latent vector $\underline{c}\in\{1,\ldots,N_{G}\}^{n}$; assignment to a jump step is denoted by the latent vector $\underline{b}\in\{0,1\}^{n}$, the assignment of a return step is denoted by the latent vector $\underline{\eta}\in\{0,1\}^{n}$, and the assignment to an activity area is denoted by $\underline{z}\in\{(g,l)\ :\ g=1,\ldots N_{G};\ l=1,\ldots,L_{g}\}$. The classification $\underline{z}$ represents a return to the activity space associated with the $l^{\text{th}}$ contiguous Brownian motion in activity $g$. Let $\underline{\theta}\equiv\{\underline{\mu}_{g},\mathbf{\Sigma}_{g}^{*},\alpha,\vartheta,\omega_{g}\}_{g=1}^{G}$ and $\underline{\mathfrak{h}}\equiv\{\underline{b},\underline{c},\underline{\eta},\underline{z}\}$. The mixture likelihood in Eq. (23) becomes

[TABLE]

where

[TABLE]

A Bayesian framework for inference in the LFCM is defined by the prior specification given in Table 1 and follows in part the priors for finite latent mixtures parameters of Ryan et al. (2017). These priors are chosen to be conjugate and non-informative to lead to efficient MCMC samplers from the corresponding joint posterior distributions. Nuisance parameters can be integrated out as needed. Ryan et al. (2017) find the prior distribution $N_{G}\sim\text{Gamma}(\tfrac{1}{2},\tfrac{1}{2})$ works well in practice for absorbing existing components and ejecting new components in a finite mixture. The minimal jump distance $\varepsilon$ for classification as a jump governs the scale of the model, and is assumed to be known to obtain conjugacy. The determination of appropriate values for $\varepsilon$ is discussed in Sections 6 and 7.

When interest lies only in activity regions, a collapsed MCMC sampler can be derived by integrating almost all parameters of the LFCM – see Appendix A.

6. Simulation study

We evaluate the efficacy of the LFCM model to represent the distribution of activity of an individual by simulating the trajectory of a simple routine in which, on a daily basis, an individual goes from home, to work, to a public space, and back home. With varying amounts of observed data, we show estimates of activity distributions for this example individual using the LFCM and a competing estimator called the conservative proportional time (CPT) estimator that is constructed conditional on a spatial grid (Dong et al., 2020). We show that the LFCM consistently captures the three key locations visited by the example individual, while the estimates produced by the CPT estimator exhibit more variability because they depend on the size of the cells in the spatial grid used. In addition, we show that the Pareto hyperparameter associated with minimal jump size does effect the behavior of the LFCM—in particular, it may incorrectly identify a jump process as an activity space when set to large values.

Data Generation

We simulate the spatiotemporal trajectory of an individual during their activities of daily living. They first spend their time at a “home” location from 0:00 to 8:40 represented by the coordinates $(0,0)$. They then travel for 20 minutes to a “work” location at the coordinates $(1,1)$, and stay there from 9:00 to 17:00. They then they travel for 15 minutes to a public space at coordinates $(1,0)$, and stay there from 17:15 to 23:45. Finally, they take 15 minutes to return home. When at home and at work, the individual moves around the area according to a zero-drift Brownian motion $\sigma B(t)$ with standard deviation $\sigma=0.1$. When in public space, the individual moves similarly with $\sigma=0.2$. The travel between these three key locations occurs on the shortest line between them.

We generate a new GPS location for the example individual every 10 minutes when at home and at work, every 5 minutes when in public space, and every 60 seconds when traveling between any two key locations. This reflects a real-world location data collection process in which individuals are more likely to use GPS tracking applications during high-movement activities (e.g., fitness tracking) or travel compared to time periods with low or no spatial movement. GPS locations during travel are simulated from a uniform distribution between the source and destination to simulate the possibility of traffic. Figure 12 displays the generated data and illustrates the simulated routine.

We generate three subsets of the simulated data shown in Figure 12: the entire set of simulated data, a random subset consisting of 50% of the simulated data, and a random subset consisting of 25% of the simulated data. We used these three subsets to estimate the activity distributions with LFCM together with the CPT estimator of Dong et al. (2020).

The Conservative Proportional Time Estimator (CPT)

The CPT estimator is a discrete estimator of an individual’s activity distribution. We assume that a spatial window $\mathcal{A}$ in which the individual is active is divided into cells $\mathcal{C}$ of equal size. Denote by $\mathcal{X}=\{\underline{x}(t_{j}):j=1,\ldots,n\}$ the GPS locations of this individual. The probability that the individual is located in the grid cell $a\in\mathcal{C}$ is proportional to the time spent in this grid cell (Dong et al., 2020):

[TABLE]

This estimator is conservative because the time spent traveling between two grid cells is excluded from the calculation. Chen and Dobra (2020) define the density ranking $\widehat{\alpha}(\underline{x})$ as

[TABLE]

where $\widehat{p}(\underline{x})$ is defined as $\Sigma_{a\in\mathcal{C}}\pi_{a}\mathbf{1}\{\underline{x}\in a\}$. The density ranking function $\widehat{\alpha}(\underline{x})$ is the fraction of observations in $\mathcal{X}$ whose estimated density is lower than the estimated density of the given point $\underline{x}$. The density ranking of a point $\underline{x}$ is interpreted as $\underline{x}$ being in the top $100\times[1-\widehat{\alpha}(\underline{x})]\%$ of locations that an individual is observed visiting. This can be used to generate level sets of the form

[TABLE]

corresponding to the top $\gamma\times 100\%$ activity regions.

Results

We used the MCMC Algorithm 2 described in Appendix A to estimate activity types, jump points, and activity regions with the LFCM. As a point of comparison, the level set $\widehat{A}_{0.9}$ is estimated using the CPT in Eq. (25). Figure 13 displays the results of LFCM (left) where Brownian motion is identified by square points, jumps are identified by triangular points, and shaded areas are identified activity regions; and the CPT method of Dong et al. (2020) (right). The top row corresponds to all data, the middle row corresponds to a random selection of 50% of the data, and the bottom row corresponds to a random selection of 25% of the data.

While the LFCM and the CPT recover the activity regions associated with home $(0,0)$, work $(1,1)$ and public space $(1,0)$, LFCM achieves a Jensen-Shannon divergence of $1.06\times 10^{-5}$ relative to the generative distribution, while CPT has a Jensen-Shannon divergence of 0.013. Figure 13 shows the sensitivity of CPT (Dong et al., 2020) to the implicit grid size. Indeed, changes in grid size can have substantial effects on the activity regions identified. Figure 14 illustrates three distinct grid sizes that substantively change the level set $\widehat{A}_{0.9}$. From left to right, the grid sizes used are 0.15, 0.20, and 0.25. A grid size of 0.15 identifies additional jump points as activity regions while shrinking the area of true activity, while a grid size of 0.25 removes one of the three areas of activity altogether.

LFCM is subject to less sensitivity, as the majority of the parameters in the LFCM can be analytically integrated out. However, LFCM has sensitivity to the a priori minimal jump length parameter $\varepsilon$. The simulation above uses a minimal jump length value $\varepsilon$ of $0.1$ based upon visual inspection of the data. If this value is set low, a number of points in the activity space may be labeled as jumps, while if the value is high, origin-destination travel may be labeled as Brownian motion activity. Figure 15 illustrates this behavior. The left panel shows LFCM results on the simulated data when $\varepsilon=0.01$, the middle panel shows LFCM results with $\varepsilon=0.10$, and the right panel shows LFCM results with $\varepsilon=0.20$. In the left panel, a substantial number of the Brownian motion points at $(1,0)$ are instead classified as jump points, while in the right panel, a portion of travel is classified as an activity space. For this reason, a careful calibration of the minimal jump length is needed by fitting the LFCM on a grid of possible values for $\varepsilon$.

Persistence Curves

The activity regions estimated by the LFCM naturally provide a probability density from the mixture of Brownian motions which generate them. We use this estimated density, along with the discrete density estimated with the CPT estimator to generate the persistence curve in Figure 16, which illustrates the level sets $\widehat{A}_{\gamma}$ for $\gamma\in[0,1]$.

The persistence of a connected component is defined by its birth and death times, which are determined by the corresponding level of density ranking. A decrease in the number of connected components suggests that either a connected component disappeared or two components merged. The persistence curve also provides evidence of mode stability; modes that persist over a larger range of $\gamma$ values are more robust, suggesting that they represent a true feature of the underlying data distribution. While the CPT identifies 4 components, the LFCM identifies three individual connected components similarly to the original data generation mechanism. Moreover, each of the modes identified by the CPT has a lower persistence than those identified by the LFCM, suggesting the latter provides a more stable estimate of the true activity spaces.

7. Analysis of mobile device GPS data

In this section we provide an application of the LFCM to real GPS data of the individual devices described in Section 2. First, to qualitatively illustrate the estimation output, we use the LFCM to show the types of activity and the activity activity spaces for a single device. Second, to demonstrate the LFCM’s efficacy of anonymized path generation, for the same device we display the generated sample path using the identified activity regions and draws of the posterior parameters returned by the model. Third, we provide a quantitative comparison across all weeks and individual devices between the generated sample paths and the true mobility data on the mobility metrics reviewed in Section 3: jump length, mean squared displacement (MSD), radius of gyration and frequented locations. Fourth, we show that the proportion of new locations visited over time of the simulated trajectories are aligned with prior findings in the literature. Fifth, we qualitatively and quantitatively assess the stability of activity spaces identified by the LFCM by estimating the model over progressively longer time intervals. Finally, conditional on the known trajectories, we use the LFCM estimate to simulate and interpolate a continuous trajectory for each individual, allowing for direct comparison across non-uniformly recorded trajectories, showing that the LFCM can be used to estimate activity overlap.

Activity region estimation

To obtain draws of posterior parameters for an example individual device’s activity space, we obtain posterior samples of the parameters of the LFCM using 10,000 iterations of Algorithm 2 described in Appendix A based on their 1,116 GPS locations recorded in four weeks. We generate activity regions using uncertainty bounds for the individual’s unobserved time given the posterior Brownian motion parameters $(\underline{\mu}_{g},\mathbf{\Sigma}_{g}^{*})$ of an activity. Figure 17 displays the entire movement of the example individual device over the four-week period. The colored points represent contiguous Brownian paths within the same activity group $g$, and the shaded regions represent the activity regions.

Each of the gray points in Figure 17 represents a location record identified as a non-Brownian jump in the observed process, i.e. as having a low probability of being a member of any of the activity groups after accounting for the time between records, usually occurring when a point is spread spatially far from neighboring data. The number of jumps is controlled by the hyperparameters $(\alpha,\varepsilon)$ of the Pareto distribution. The parameter $\varepsilon$ controls the minimum distance from the prior location that can be identified as a jump, and the parameter $\alpha$ controls the probability concentration that a far point is classified as a jump relative to a Brownian motion.

The shaded regions in Figure 17 are activity regions using LFCM representing the potential range of the device’s movement between records while in the area after accounting for drift and dispersion. Specifically, the region represents a 95% confidence region assuming the Brownian bridge movement model in Eq. (11) with dispersion $\mathbf{\Sigma}_{g}^{*}$ and time $\Delta t$ between each point identified as in the same activity regions. Activity regions may appear large for two reasons: the magnitude of the dispersion parameter within the activity group $g$ is high, or the time spent in the area is high. The constructed activity regions need not be elliptical, but activities with low dispersion over a long period of time tend to create elliptical regions of uncertainty.

Model record generation

We show that simulations from the LFCM can be used to generate representative trajectories with characteristics that match the mobility of the individual device. Based on posterior samples from the LFCM, we estimate the jump length distribution, the activity type distribution, the mean and variance of a Brownian motion describing the activities’ movement. Areas of frequent return are a function of these quantities.

Using one week of data from the individual device described above, we draw samples from the joint posterior distribution of the LFCM. We use these samples to determine MAP parameter estimates, identify activity regions and generate sample paths. Figure 18 displays the trajectory of the device in 500 points simulated using the MAP estimates, where Brownian motion movement is shown in red, Lévy jumps are shown in gray, and the identified activity regions are shaded in blue. This simulation can be compared with Figure 17, which is the true mobility of the device in the same week—the simulated points appear to be very similar to those in Figure 17, which are the true mobility data.

This simulation process also results in reasonable mobility metrics mentioned in Section 3. Table 2 displays the mean and standard deviation of jump length, MSD, and radius of gyration across 12 weeks of movement and 293 individual devices for the observed data, data simulated from LFCM, and data simulated from the conservative proportional time estimator in Eq. (25). Whereas the LFCM closely reproduces jump length, MSD, and radius of gyration, simulations from the CPT estimator tend to exhibit under-bias on mean jump length and MSD.

The noticably better performance of the LFCM compared to the CPT estimator can be explained by a couple of factors. First, the likelihood of movement is constant within grid cell, and some cells are much more likely than others. This causes a large amount of the trajectory to be within-cell movement. Second, region size is uncorrelated with time spent in a region for the CPT estimator. In contrast, for the LFCM, with the exception of certain key locations, the likelihood of a region is correlated with the size of the region; nearby areas of return tend to merge if there is sufficient Brownian motion activity in those areas.

Frequented Locations

One of the primary difficulties in the estimation of individual mobility using Lévy flight models is faithful replication of an individual’s frequented locations. This is often visually assessed in the literature as the number of new locations an individual visits over the time period of data collection. Similarly, we plot the proportion of new locations visited over a week-long period of data collection–temporally normalized to the unit interval $[0,1]$–estimated by the LFCM in Figure 19; using MAP estimates estimated by the LFCM for six individual devices, we generate new trajectories conditioning on the observed times seen in the data. The figure shows a reference line of $t^{0.6}$, corresponding to the observed rates of exploration in the literature (Song et al., 2010), while the average across our simulations is $t^{0.52}$.

Activity Space Stability

While individuals might continually visit new locations, the number of places they spend most of their time stays relatively stagnant. We illustrate this in Figure 20, which displays an individual device’s activity spaces estimated by LFCM after 1 week (left), 4 weeks (middle), and 12 weeks (right) of observation. While the individual’s activity distribution becomes more concentrated, and the areas in which the individual spends the majority of their time becomes more obvious over time, the most frequented activity spaces appear to be stable after a single week; the activity regions around $(-1.5,1)$, $(0,0)$, and the multiple regions between $(-0.5,-2)$ and $(0.5,-1)$ are present across all estimates.

Figure 20 also shows that along with persistent areas, new areas are identified over time. However, adults living in urban areas tend to have two or three most relevant places: work, home, and a primary social activity (Oldenburg and Brissett, 1982; Finlay et al., 2019). Although work may require significant travel, the few areas of highest activity are likely to be the most relevant to understanding patterns of habitual activity. To assess the extent LFCM can capture this behavior, we empirically evaluate convergence of the “Top 1” and “Top 3” activity areas as a function of the observed time.

Specifically, let $\widehat{p}_{i,w}$ denote the estimated activity distribution for device $i$, estimated on points $\underline{x}_{1},\ldots,\underline{x}_{w}$. First we define the $\gamma$ upper level set of distribution $\widehat{p}_{i,w}$ as

[TABLE]

For each device we use the LFCM to estimate $\widehat{\mathcal{A}}_{i,w}^{0.8}$ in Eq. (27) from $w\in\{2,4,\ldots,12\}$ weeks of data points. For each device and week, we select the “Top 1” and “Top 3” disconnected subsets $A\in\widehat{\mathcal{A}}_{i,w}^{0.8}$, according to highest probability mass. Figure 21 displays activity areas estimated with increasingly longer temporal data for an individual device. Both the “Top 1” and “Top 3” activity area estimates anecdotally corroborate the patterns shown in Figure 20; while some spatial error appears in the “Top 3” activity areas at week 4, the activity areas are generally consistent with even two weeks of data.

Finally, to construct a quantitative measure of convergence for the LFCM activity space estimates, we use the LFCM to estimate $\widehat{\mathcal{A}}_{i,12}^{0.8}$ on individual devices’ full 12 weeks of data, and their corresponding “Top 1” and “Top 3” activity areas, as the “true” activity space. Then, for each device, we estimate the similarity between identified activity areas at week 12 and at week $w\in\{2,4,6,8,10\}$. Several similarity metrics described in Appendix C can be applied to assess convergence in the estimated activity areas: the KL divergence, the Hausdorff distance, the Jaccard distance, and a modification of the Jaccard distance which we call the Overlap. The choice of metric depends on the representation and modeling of these areas. Generally, the Overlap and Jaccard distances are likely to be the least sensitive, especially to minor alterations in large activity spaces.

The empirical mean and 95% confidence band of the Overlap distance across devices for each week $w$ are shown in Figure 22. We found that, on average, “Top 1” and “Top 3” activity areas estimated on 2 weeks of data have approximately 50% overlap with the corresponding estimates on 12 weeks of data. Using 4 weeks of data yields an average 10% increase in overlap for both the “Top 1” and “Top 3” activity areas — showing approximately 60% overlap with 12-week activity areas. Finally, plots for both “Top 1” and “Top 3” show confidence interval lower bounds of approximately 20% for 2-week estimates, suggesting substantial overlap with top activity areas estimated on the full 12 week data even in the worst case.

Probabilistic Contact Networks

In social and behavioral sciences, understanding the probabilistic contact networks can aid in studying social interaction patterns, shared behaviors, and commonalities in lifestyle choices; in public health, these networks may imply similar exposures to health risks (e.g., pollution, noise), which can be important when studying disease spread or common health issues in a community; in urban planning and development, understanding these networks can inform the planning of public amenities, recreational areas, and even residential or commercial zones based on the population’s behaviors and patterns of movement. In what follows, we leverage Bayesian parameter estimates of the LFCM to calculate these probabilistic contact networks.

First, across all individual devices $i$ and weeks $w$, we find the minimum and maximum times across individuals for that week.

[TABLE]

and the final observation time as

[TABLE]

Then, the minimal initial observation time across all individuals in week $w$ is given by

[TABLE]

and the maximal end observation time is

[TABLE]

For each individual, interpolation and extrapolation are done by generating conditional sample paths from posterior scans of the LFCM—sampling randomly from the sequence of scans

[TABLE]

from the target posterior distribution $p(\theta\mid D)$, after convergence is achieved. Figure 23 illustrates the resulting trajectory for one individual and week.

As a baseline comparison to this method, we also generate a simple interpolation and extrapolation strategy, which linearly interpolates between observed records, and extended to $t_{w}^{\min}$ and $T_{w}^{\max}$ by assuming that the individual remains still in their first and last observed record, respectively.

For each scan, we then calculate the $L_{2}$ distance between individual trajectories across all devices for a given week, resulting in an activity overlap matrix. We generate 500 such matrices using 500 randomly sampled posterior scans of the individuals, which corresponds to a potential trajectory that the individual might follow over the observed period.

The left panel of Figure 24 shows the mean pairwise distance matrices for 293 devices estimated on 500 generated sample paths from individual/week LFCM models. The right panel shows the same matrix after linear interpolation.

There is relative agreement between both interpolation/extrapolation strategies, which reveal significant structure in the distance matrices, holding valuable implications for social and behavioral science. The clear vertical and horizontal lines suggest that some individuals remain consistently far from others - potentially indicating social isolation or distinct spatial domains - while small dark patches point to local clusters, which may reflect subgroups exhibiting confined or periodic movement. Note that the LFCM-based distance matrix has relatively fewer outliers, which may be due to the poor extrapolation behavior of the linear method. This spatial patterning may offer critical insights to understand community structure, mobility, and the design of targeted interventions and policies.

8. Discussion

We have shown that the Lévy flight cluster model (LFCM) robustly captures the intricate dynamics of human mobility, as evidenced by both simulation studies and real-world GPS analyses. Our proposed model successfully delineates probabilistic activity spaces by distinguishing between localized activity clusters characterized by stationary repetitive behavior and sporadic long-range movements consistent with Lévy flight dynamics. These dual components were reliably estimated despite irregular sampling, and the posterior distributions aligned well with established mobility metrics such as jump length, radius of gyration, and mean square displacement. In addition, the ability of the model to generate probabilistic sociomatrices via interpolation and extrapolation further underscores its efficacy in representing temporal and spatial overlap among individuals. This capacity is pivotal for examining phenomena such as social isolation, geographic mobility constraints, and the emergence of localized subgroups, thus offering nuanced insights into patterns of social exclusion and differential access to resources.

Overall, the LFCM provides a substantial methodological advancement over traditional approaches by accurately modeling human mobility through a generative Bayesian framework that integrates both local Brownian motion and long-range Lévy jumps. The success of the model in reproducing individual-level activity distributions, while simultaneously offering uncertainty bounds, highlights its potential as a tool for both rigorous scientific inquiry and practical applications in urban planning, epidemiology, and data anonymization. By allowing researchers to simulate realistic activity spaces and derive probabilistic overlap measures even under irregular sampling conditions, the LFCM lays a comprehensive foundation for future investigations into the behavioral and social dimensions of human mobility.

Acknowledgment

AHW and GSC thank the following university departments for their in-kind support through their affiliate faculty positions: Department of Statistical Sciences and Operations Research at the Virginia Commonwealth University (AHW, GSC), Department of Statistics at the University of Washington (GSC), and the Department of Statistics and Actuarial Science at the University of Waterloo (GSC).

Appendix A Collapsed MCMC sampler for the LFCM

The LFCM can be “collapsed” through integrating out nuisance parameters, significantly lessening computational complexity. Below we derive the collapsed MCMC sampling algorithm for the LFCM. Let $n_{g}=\sum_{i=1}^{n}\ \mathbf{1}\{c_{i}=g\}$, $n_{b}=\sum_{i=1}^{n}b_{i}$, $n_{\eta}=\sum_{i=1}^{n}\eta_{i}$, $\mathfrak{Z}=\{(g,l)\ \mid\ g=1,\ldots,N_{G};\ \ l=1,\ldots,L_{g}\},$ and $n_{z}=\sum_{i=1}^{n}\mathbf{1}\{z_{i}=(g,l)\}$. Let $\mathbf{\Delta X}$ represent $n$ observed spatial movements and let $\mathfrak{H}$ represent the collection of $n\times 1$ latent vectors $\underline{\mathfrak{h}}$. The joint distribution of the data and model parameters can be written as

[TABLE]

If one is only interested in the latent vectors and hyperparameters associated with the Brownian motion, the corresponding marginal joint distribution is:

[TABLE]

The hyperparameters associated with the Brownian motion can also be integrated out of the marginal joint distribution – Appendix B:

[TABLE]

where \mathbf{K}_{\accentset{\smash{\raisebox{-0.75346pt}{\scriptscriptstyle\sim}}}{n}_{g}}\equiv(\kappa\mathbf{I}+n\overline{\mathbf{B}}^{\top}\overline{\mathbf{B}}), \nu_{\accentset{\smash{\raisebox{-0.75346pt}{\scriptscriptstyle\sim}}}{n}_{g}}\equiv\nu+\tilde{n}_{g}, and

[TABLE]

Here $\tilde{n}_{g}=\#\{\underline{\Delta x}_{i}:(c_{i},b_{i})=(g,0)\text{ or }(b_{i},\eta_{i},z_{i})=(1,0,(g,l))\}$ is the number of spatial movements that are classified as a Brownian motion movement or a return, and the matrices $\mathbf{A},\mathbf{B},\mathbf{D}$ are given by $\mathbf{B}_{i}=\mathbf{A}_{i}^{1/2}\mathbf{D}_{i}$,

[TABLE]

The marginal joint distribution becomes

[TABLE]

with $\pi(N_{G})$ the marginal distribution on the number of groups $N_{G}$,

[TABLE]

and

[TABLE]

Hence, the joint distribution $f(\Delta\mathbf{X},\mathfrak{H},N_{G},\Delta\underline{t})$ is discrete in the latent vectors and the number of groups—$(\underline{c},\underline{b},\underline{z},\underline{\eta},N_{G})$—and stochastic searches over the discrete space of finite mixtures can be used to sample from the posteriors without reversible jumps (Ryan et al., 2017; Richardson and Green, 1997). This discrete search technique has been successful in a number of finite mixture applications (e.g., Ryan et al., 2017; Wyse and Friel, 2012; Nobile and Fearnside, 2007). The exact probability of the proposed movement can be calculated as

[TABLE]

and often reduces to computationally simpler forms through cancelation of the non-varying terms from Eq. (30), depending on the update being considered.

Conditional distribution of $\eta$

The conditional distribution of $\eta_{i}$ is constructed as a Beta-Bernoulli($\mathpzc{p,b_{1},b_{2}}$). Hence,

[TABLE]

Conditional distribution of $N_{G}$

Conditional sampling of the number of activity groups $N_{G}$ is performed using the absorption / ejection method of Ryan et al. (2017). Let $p^{\text{ej.}}=p^{\text{ab.}}=0.5$. If $N_{G}>1$ and $N_{G}<M_{G}$ for some maximum number of groups $M_{G}$, then with probability $p^{\text{ej.}}$ a “eject” move is made, and with probability $p^{\text{ab.}}$ an “absorb” move is made. If $N_{G}=1$, the eject move is made with probability $1$, and if $N_{G}=M_{G}$, the absorb move is made with probability $1$.

The eject move chooses one of the $N_{G}$ clusters at random and attempts to reallocate a random number of the $L_{g}$ continuous paths in $g$. The number of paths with probability $\varrho\sim\text{Beta}(\mathpzc{a},\mathpzc{a})$ is selected to move from $g$ to a new cluster, $N_{G}+1$. The proposal, denoted below by the superscript $`*^{\prime}$, is accepted with probability:

[TABLE]

where $\pi(N_{G})$ represents the marginal distribution on the number of groups$N_{G}$.

The absorb move is performed by randomly selecting two components $g_{1},g_{2}$ from the $N_{G}+1$ groups and combining them into a single component. This move is accepted with probability:

[TABLE]

MCMC Sampling

Sampling of the marginal posterior distribution

[TABLE]

in Eq. (30) is conducted using a Markov chain Monte Carlo (MCMC) sampling algorithm. A single sweep of parameter updates (see Table 1 in the main text) are described in Algorithm 2 below. We start with two types of movement data. First, we have the spatial difference data, which tells us how far an individual has moved. Second, we have the data for the time difference, indicating the time it took for these movements to take place. We are interested in grouping movements into what we call “activity groups” or categories of movement such as walking, running, or jumping. We also want to understand specific characteristics of these movements, for example, whether they correspond to a trip into a new spatial region or whether they represent a return trip to a spatial region that has been previously visited.

We start with an initial assignment of movements into groups and their classifications (e.g., exploratory vs. return movements) that is generated by sampling from a hypergeometric distribution. Next we sequentially cycle through the following steps:

Update the group classifications. For each group, we look at every movement and decide if its current group assignment is good or if it needs to be updated. This decision is based on the spatial and temporal data and the characteristics of the movements.

• 2.

Adjust the number of groups. We consider if we need more or fewer groups to accurately represent the variety of movements.

• 3.

Refine the movement characteristics. For each movement, we revisit our initial guesses about its characteristics (like whether it’s an exploration or a return) and update these guesses to better match the data.

• 4.

Classify the activity areas. We decide if each movement is towards a new area (exploration) or a known area (return). This helps in understanding how individuals or objects interact with their environment over time.

With each repetition of these steps, the movement groupings and their classification are refined. This algorithm is a systematic way to sift through data on movements, categorizing them into different types of activities, and understanding their nature.

Appendix B Derivations of posterior distributions

For positive definite symmetric matrices $\mathbf{D},\ \mathbf{A}$, define the prior normal and Wishart distributions

[TABLE]

The corresponding densities are proportional to

[TABLE]

The exponent of Eq. (B) can be written as

[TABLE]

Letting $\mathbf{B}_{i}:=\mathbf{A}_{i}^{1/2}\mathbf{D}_{i}$, the quadratic expression above is expanded to

[TABLE]

where the overscore of a matrix, e.g. $\overline{\mathbf{A}}$, represents the mean operator. The term in Eq. (B) can be combined with the exponent of Eq. (B) to become

[TABLE]

Since $\overline{\mathbf{B}}\mathbf{\Lambda}=\mathbf{\Lambda}\overline{\mathbf{B}}$ when $\mathbf{B}_{i}=\mathbf{A}_{i}^{1/2}\mathbf{D}_{i}$, this becomes

[TABLE]

Adding and subtracting the term

[TABLE]

allows for the factorization

[TABLE]

Adding and subtracting the term $N\overline{\mathbf{A}^{1/2}\underline{x}}^{\top}\mathbf{\Lambda}\overline{\mathbf{A}^{1/2}\underline{x}}$ allows factorization

[TABLE]

Combining the other two terms:

[TABLE]

This is a scalar and is equal to its own trace. Moreover, since the trace has a cyclic commutative property, we have

[TABLE]

and

[TABLE]

These can be combined with Eq. (B) so that the exponent becomes

[TABLE]

The posterior parameters are then

[TABLE]

When $\mathbf{A}=\mathbf{D}=\mathbf{I}$, the posterior distribution becomes the standard posterior parameters of a Normal-Wishart distribution. The normalizing constant is

[TABLE]

For $\{\underline{\Delta x}_{i}\ :\ c_{i}=g\}$,

[TABLE]

For $\{\underline{\Delta x}_{i}\ :\ b_{i}=1,\eta_{i}=0,z_{i}=(g,l)\}$,

[TABLE]

Appendix C Similarity metrics for activity areas

Several similarity metrics can be applied in the context of evaluating the convergence of estimated activity areas. In this section we discuss the Kullback-Leibler (KL) divergence, the Hausdorff distance, the Jaccard distance, and finally, a modification of the latter called the Overlap.

The Kullback-Leibler (KL) divergence provides a measure of the difference between two probability distributions. Let $P$ and $Q$ be two probability distributions defined on the same measurable space. The measure quantifies how much information is lost when the distribution $P$ is approximated by $Q$. For discrete and continuous variables, respectively, it is defined as

[TABLE]

Its formula is essentially an expectation of the logarithmic difference between the probabilities $P$ and $Q$, where $P$ and $Q$ are the two distributions being compared. This metric is appropriate when “activity areas” are conceptualized as probability distributions. However, its sensitivity is notably pronounced when one distribution assigns a significant probability to an event that the other largely disregards, leading to potential large divergences even for small differences. As seen in Figure 20, the number of activity areas naturally increases over time, leading to difficulty in interpreting this metric.

The Hausdorff distance gauges the maximum distance between two point sets, encapsulating the idea of shape dissimilarity. It captures the greatest distance from a point in one set to the closest point in the other set. This metric is inherently sensitive to outliers, as even well-aligned datasets can yield a high Hausdorff distance if a single point pair exhibits a large distance. Due to the sparsity of the data, activity areas in early weeks are likely to have far greater spread than those in later weeks.

The Jaccard distance, derived from the Jaccard index, evaluates the dissimilarity between two discrete sets. Defined as

[TABLE]

this metric quantifies the proportion of non-overlapping elements in two sets. It is particularly suitable when “activity areas” are represented as discrete locations. Given its foundational basis on ratios, the Jaccard distance offers a balanced sensitivity to changes, especially in large datasets where individual differences might have proportionally smaller impacts.

Finally, we propose the Overlap distance $1-\frac{|A\cap B|}{|B|}$ which computes the proportion of elements in set $B$ that are absent from set $A$, serving as an asymmetric dissimilarity measure. Contrasting this with the symmetric Jaccard distance, which considers mutual overlaps, the asymmetric measure emphasizes the extent to which $A$ covers $B$. Specifically, it captures the dissimilarity value when observing $A$ relative to $B$, focusing on the coverage of $B$ by $A$. In applications such as evaluating “activity areas”, this measure provides insight into the degree of coverage of the activity areas estimated with 12 weeks of data by the activity areas estimated with a fewer number of weeks, offering a practical perspective of the coverage of “true” activity areas, with less sensitivity than the Jaccard distance.