Generating Synthetic GPS Traces with Markov Models

Discretize space into a tessellation, estimate a smoothed first-order transition matrix from real traces, then sample new cell sequences and map each cell back to a representative coordinate — producing synthetic GPS traces whose aggregate structure matches the real data without copying any individual.

Core Calculation # permalink

Let the study area be partitioned into $m$ cells $\{s_1, \dots, s_m\}$. Each cell is a Markov state. From real traces you count $n_{ij}$, the number of transitions observed from cell $i$ to cell $j$. The maximum-likelihood transition probability is:

This raw estimate is unsafe for privacy: a count of $n_{ij} = 1$ means exactly one person made that move, and sampling it reproduces their behavior. Apply Laplace (additive) smoothing with a pseudocount $\lambda$:

where $n_{i\cdot} = \sum_k n_{ik}$ is the total outgoing count from cell $i$. Smoothing shifts probability mass onto never-observed transitions, so a lone rare move can no longer dominate its row. The long-run visit frequency is the stationary distribution $\boldsymbol{\pi}$, the left eigenvector satisfying:

Worked numeric example # permalink

Take a 4-cell tessellation $\{A, B, C, D\}$ with the observed transition-count matrix (rows = origin, columns = destination):

With $\lambda = 1$ and $m = 4$, each row denominator becomes $n_{i\cdot} + 4$. Row A gives:

The two unobserved moves out of A ($A \to C$, $A \to D$) now carry 8.3% each instead of zero — the smoothing that prevents memorization. The full smoothed matrix is:

Solving $\boldsymbol{\pi}\hat{P} = \boldsymbol{\pi}$ yields the stationary visit distribution:

So over a long synthetic walk, cells B, C, and D are visited slightly more than A — consistent with A being a low-connectivity origin. Any synthetic trace you sample will reproduce these aggregate visit frequencies without echoing a specific real path.

Python Implementation # permalink

The function builds an H3 tessellation, estimates the smoothed transition matrix from real cell sequences, samples synthetic sequences, and converts cells back to representative coordinates. Input coordinates are WGS84 (EPSG:4326); the representative-point step reprojects to a state plane CRS for metric fidelity.

from __future__ import annotations
import h3
import numpy as np
import geopandas as gpd
from shapely.geometry import Point

def synthesize_markov_traces(
    real_sequences: list[list[str]],
    n_traces: int,
    trace_length: int,
    smoothing: float = 1.0,
    target_crs: str = "EPSG:2229",   # e.g. California State Plane Zone 5 (feet)
    seed: int | None = None,
) -> gpd.GeoDataFrame:
    """Fit a smoothed first-order Markov model over H3 cells and sample traces.

    Args:
        real_sequences: real traces as ordered lists of H3 cell ids (states).
        n_traces: number of synthetic traces to generate.
        trace_length: number of steps per synthetic trace.
        smoothing: Laplace pseudocount lambda added to every transition. Larger
            values spread probability onto unobserved moves, so a single rare
            transition cannot be reproduced verbatim. This is the privacy knob.
        target_crs: projected CRS for the emitted representative coordinates.
        seed: RNG seed for reproducibility.

    Returns:
        GeoDataFrame with columns [trace_id, step, h3_cell, geometry] in target_crs.
    """
    rng = np.random.default_rng(seed)

    # 1. State space: every distinct H3 cell observed across all real traces.
    states = sorted({cell for seq in real_sequences for cell in seq})
    idx = {cell: i for i, cell in enumerate(states)}
    m = len(states)

    # 2. Count transitions, then Laplace-smooth so no single observed move
    #    (count == 1, i.e. one person) survives as a deterministic edge.
    counts = np.zeros((m, m), dtype=float)
    for seq in real_sequences:
        for a, b in zip(seq[:-1], seq[1:]):
            counts[idx[a], idx[b]] += 1.0
    smoothed = counts + smoothing
    P = smoothed / smoothed.sum(axis=1, keepdims=True)   # row-stochastic

    # 3. Initial-state distribution from smoothed first-step frequencies.
    start_counts = np.zeros(m)
    for seq in real_sequences:
        if seq:
            start_counts[idx[seq[0]]] += 1.0
    pi0 = (start_counts + smoothing) / (start_counts.sum() + smoothing * m)

    # 4. Sample synthetic cell sequences by walking the chain.
    rows = []
    for t in range(n_traces):
        cur = rng.choice(m, p=pi0)
        for step in range(trace_length):
            cell = states[cur]
            # Representative point = H3 cell centroid (lat, lng in WGS84).
            lat, lng = h3.cell_to_latlng(cell)
            rows.append({"trace_id": t, "step": step, "h3_cell": cell,
                         "geometry": Point(lng, lat)})
            cur = rng.choice(m, p=P[cur])

    gdf = gpd.GeoDataFrame(rows, geometry="geometry", crs="EPSG:4326")
    # Reproject cell centroids to a metric CRS for downstream distance work.
    return gdf.to_crs(target_crs)

Emitting the cell centroid rather than a real observed coordinate is deliberate: every synthetic point snaps to the H3 lattice, so the output cannot leak a precise real location. If even the centroid is too precise for the smoothing you chose, add calibrated jitter using Laplace and Gaussian noise for coordinate data, and account for it against your privacy budget (epsilon) for spatial queries.

Verification # permalink

Confirm the synthetic set does not reproduce any individual. For every real trace, find its nearest synthetic trace; an implausibly small distance flags memorization.

import numpy as np
from scipy.spatial import cKDTree

def membership_inference_check(
    real_xy: np.ndarray,        # (N, 2) real points, projected metric CRS
    synth_xy: np.ndarray,       # (M, 2) synthetic points, same CRS
    leak_threshold_m: float = 1.0,
) -> dict:
    """Flag real points with a near-exact synthetic match (memorization)."""
    tree = cKDTree(synth_xy)
    dist, _ = tree.query(real_xy, k=1)          # nearest synthetic per real point
    leaked = float(np.mean(dist < leak_threshold_m))
    return {
        "min_distance_m": float(dist.min()),
        "median_distance_m": float(np.median(dist)),
        "leak_fraction": leaked,                 # target: ~0.0
    }

# A leak_fraction near zero and a min distance comfortably above the tessellation
# resolution indicate the Markov model generalized rather than memorized.

Pair this with a utility check: the synthetic stationary visit frequencies should match the real per-cell visit counts (correlation close to 1), and synthetic trip-length and radius-of-gyration distributions should track the real ones.

Edge Cases & Adjustments # permalink

• Sparse transitions. Rows whose entire outgoing count comes from one or two users memorize those users even after light smoothing. Raise $\lambda$, coarsen the tessellation so each cell aggregates more movement, or prune states with fewer than $k$ contributing users before fitting.
• Higher-order dependence. A first-order chain forgets routines longer than one step (home → work → gym). Second-order models capture more realism but memorize faster because each state pair is supported by fewer people; prefer adding a diary/EPR layer over raising the order.
• Temporal dynamics. A single transition matrix ignores time of day and weekday/weekend structure. Fit separate matrices per time bin, or drive location choice with a Markov diary, so peak-hour flows are realistic without encoding one person’s schedule.
• CRS mismatch. H3 centroids are WGS84 (EPSG:4326) degrees; any distance-based verification or jitter must run in a projected metric CRS. Reproject before measuring, or the leak-threshold check silently compares degrees to metres.

FAQ # permalink

How many real traces do I need to fit a Markov model?

Enough that every retained transition is supported by many different users — a practical floor is that no kept transition count comes from a single individual. Coarser tessellations reach adequate support faster because each cell aggregates more movement; smoothing and pruning handle the residual sparse states.

Why apply Laplace smoothing to the transition matrix?

A transition observed once encodes exactly one person’s move. Adding a pseudocount $\lambda$ to every cell spreads probability onto unobserved transitions, so a lone rare move no longer stands out and cannot be reproduced verbatim. It is the discrete-space analogue of adding noise to a numeric statistic.

Is a first-order Markov model enough for realistic traces?

First-order chains capture cell-to-cell structure but forget longer routines and time-of-day effects. Higher-order models are more realistic yet memorize individuals faster. For temporal realism, add a Markov diary or EPR layer rather than raising the chain order.

Related

• Synthetic Mobility Data Generation — the full generative pipeline this technique sits inside
• Trajectory Anonymization Techniques — perturb real traces instead of synthesizing new ones
• Privacy Budget Allocation for Spatial Queries — account for any noise you add on top of smoothing
• Laplace & Gaussian Noise for Coordinate Data — jitter representative coordinates when centroids are too precise
• Estimating Uniqueness of Mobility Traces — quantify how memorable the real traces are before fitting

← Back to Synthetic Mobility Data Generation