Tuning HMM Emission and Transition Probabilities

The Hidden Markov Model map matching cluster introduces sigma_z and beta as the two parameters governing emission and transition probability, and the OSRM-backed matcher implementation wires them into a runnable class — but neither page covers how to actually derive good values for a specific fleet, or why the right value shifts as GPS sampling rate changes. This page fills that gap: an estimator for sigma_z grounded in real position noise rather than a guessed constant, a transition model built explicitly around the route-length-vs-great-circle-distance gap, and a small grid-search harness for choosing beta against labelled data. For the broader question of when an HMM approach is worth this tuning effort at all versus a simpler geometric matcher, see choosing a map-matching algorithm.

Compatibility & Configuration Requirements

Production-Ready Implementation

This class estimates sigma_z from position data, exposes log-space emission and transition scoring, and includes a grid-search harness that scores candidate beta values against a labelled validation set without requiring live routing-engine calls during tuning.

import numpy as np
from dataclasses import dataclass


def _haversine_m(lat1: np.ndarray, lon1: np.ndarray, lat2: np.ndarray, lon2: np.ndarray) -> np.ndarray:
    """Vectorised great-circle distance in metres."""
    R = 6_371_000.0
    phi1, phi2 = np.radians(lat1), np.radians(lat2)
    dphi = np.radians(lat2 - lat1)
    dlambda = np.radians(lon2 - lon1)
    a = np.sin(dphi / 2) ** 2 + np.cos(phi1) * np.cos(phi2) * np.sin(dlambda / 2) ** 2
    return R * 2 * np.arctan2(np.sqrt(a), np.sqrt(1.0 - a))


class HMMProbabilityTuner:
    """
    Estimate and apply emission (sigma_z) and transition (beta) parameters
    for an HMM map matcher, entirely in log-space.

    Parameters
    ----------
    sigma_z : float | None
        GPS positional noise, metres. Pass None to require an explicit
        call to estimate_sigma_z_from_stationary before scoring.
    beta : float
        Transition decay rate, units m^-1. Larger beta penalises the gap
        between network routing distance and great-circle distance more
        harshly per metre of gap — see the tuning guidance below.
    max_candidates : int
        Candidate edges retained per GPS observation after sorting by
        perpendicular distance. Bounds the O(N^2) cost of building the
        transition matrix at each time step.
    min_sigma_z : float
        Floor applied to sigma_z to avoid a near-zero emission scale
        collapsing every candidate's log-probability toward -inf.
    """

    def __init__(
        self,
        sigma_z: float | None = None,
        beta: float = 0.02,
        max_candidates: int = 5,
        min_sigma_z: float = 3.0,
    ):
        self.sigma_z = sigma_z
        self.beta = beta
        self.max_candidates = max_candidates
        self.min_sigma_z = min_sigma_z

    @staticmethod
    def estimate_sigma_z_from_stationary(
        lat: np.ndarray, lon: np.ndarray
    ) -> float:
        """
        Estimate sigma_z as the RMS haversine deviation of GPS fixes from
        their centroid during a known-stationary period.

        Parameters
        ----------
        lat, lon : arrays of fixes recorded while the vehicle was parked
            or otherwise confirmed stationary (e.g. a depot dwell window
            identified by stop detection).

        Returns
        -------
        Estimated sigma_z in metres. Run this over several stationary
        windows across your fleet's device population and take the
        median, since GPS chipset quality varies by hardware model.
        """
        centroid_lat = np.mean(lat)
        centroid_lon = np.mean(lon)
        deviations_m = _haversine_m(
            lat, lon, np.full_like(lat, centroid_lat), np.full_like(lon, centroid_lon)
        )
        return float(np.sqrt(np.mean(deviations_m ** 2)))

    def log_emission(self, distances_m: np.ndarray) -> np.ndarray:
        """
        Log-space zero-mean Gaussian emission probability.

        log P(obs | edge) = -0.5 * (d / sigma_z)^2 - log(sigma_z * sqrt(2*pi))

        The normalising term matters when comparing scores across
        candidates scored with different sigma_z values (e.g. mixed
        device fleets); it can be dropped only if sigma_z is constant
        for every comparison in a single decode.
        """
        sigma = max(self.sigma_z or self.min_sigma_z, self.min_sigma_z)
        d = np.clip(np.asarray(distances_m, dtype=float), 0.0, None)
        return -0.5 * (d / sigma) ** 2 - np.log(sigma * np.sqrt(2 * np.pi))

    def log_transition(
        self, d_network_m: np.ndarray, d_gc_m: np.ndarray
    ) -> np.ndarray:
        """
        Log-space transition probability from the route-length vs
        great-circle distance gap, modelled as an exponential
        distribution with rate beta.

        log P(edge_j | edge_i) = log(beta) - beta * |d_network - d_gc|

        A larger beta concentrates probability mass near diff=0,
        strongly penalising any candidate whose real road-network path
        is much longer than the straight-line distance between fixes —
        exactly the shape of a genuine detour.
        """
        diff = np.abs(np.asarray(d_network_m, dtype=float) - np.asarray(d_gc_m, dtype=float))
        diff = np.clip(diff, 1e-6, None)
        return np.log(self.beta) - self.beta * diff

    def truncate_candidates(
        self, distances_m: np.ndarray, ids: np.ndarray
    ) -> tuple[np.ndarray, np.ndarray]:
        """Keep only the max_candidates closest candidates by distance."""
        order = np.argsort(distances_m)[: self.max_candidates]
        return distances_m[order], ids[order]

    def grid_search_beta(
        self,
        validation_samples: list[dict],
        beta_grid: np.ndarray,
    ) -> dict:
        """
        Score a range of beta values against labelled validation data
        and return the beta that best separates correct from incorrect
        candidate transitions.

        Parameters
        ----------
        validation_samples : list of dicts, each with keys
            'd_network_m', 'd_gc_m', 'is_correct_candidate' (bool),
            drawn from manually verified ground-truth traces.
        beta_grid : array of candidate beta values to evaluate, e.g.
            np.geomspace(0.002, 0.2, 20).

        Returns
        -------
        dict with 'scores' (beta -> margin) and 'best_beta'. The margin
        is the mean log-transition assigned to correct candidates minus
        the mean assigned to incorrect ones; a larger margin means beta
        discriminates better between right and wrong topology.
        """
        correct = [s for s in validation_samples if s["is_correct_candidate"]]
        incorrect = [s for s in validation_samples if not s["is_correct_candidate"]]

        correct_net = np.array([s["d_network_m"] for s in correct])
        correct_gc = np.array([s["d_gc_m"] for s in correct])
        incorrect_net = np.array([s["d_network_m"] for s in incorrect])
        incorrect_gc = np.array([s["d_gc_m"] for s in incorrect])

        scores: dict = {}
        for beta in beta_grid:
            self.beta = float(beta)

            correct_lp = (
                self.log_transition(correct_net, correct_gc)
                if len(correct)
                else np.array([-np.inf])
            )
            incorrect_lp = (
                self.log_transition(incorrect_net, incorrect_gc)
                if len(incorrect)
                else np.array([0.0])
            )
            margin = float(correct_lp.mean() - incorrect_lp.mean())
            scores[float(beta)] = margin

        best_beta = max(scores, key=scores.get)
        self.beta = best_beta
        return {"scores": scores, "best_beta": best_beta}

Execution & Tuning Guidelines

import numpy as np

tuner = HMMProbabilityTuner(max_candidates=5)

# Step 1 — estimate sigma_z from a known-stationary window per device
tuner.sigma_z = tuner.estimate_sigma_z_from_stationary(
    lat=stationary_df["lat"].to_numpy(),
    lon=stationary_df["lon"].to_numpy(),
)

# Step 2 — grid search beta against a labelled validation set
result = tuner.grid_search_beta(
    validation_samples=labelled_transitions,   # list[dict]
    beta_grid=np.geomspace(0.002, 0.2, 20),
)
print(f"best beta: {result['best_beta']:.4f}")

# Step 3 — score a live candidate pair in log-space
log_p = tuner.log_emission(np.array([8.5])) + tuner.log_transition(
    np.array([420.0]), np.array([310.0])
)

Sampling-rate adjustment. beta should not be a single fixed constant across a fleet with mixed sampling rates. As the interval between fixes grows, the physical distance between consecutive true positions grows too, and with it the normal amount by which network distance can exceed great-circle distance even for a correctly matched path. A practical heuristic: scale beta inversely with the expected inter-fix distance, beta ≈ c / (v_avg * Δt) with c around 3–8, so 1 Hz urban delivery traces use a stricter beta than 0.05 Hz long-haul traces with 20-second gaps between fixes.

Common Pitfalls

Up: Hidden Markov Model Map Matching in Python | Trajectory Analysis & Map Matching Techniques

Related

• Hidden Markov Model Map Matching in Python — the parent probabilistic framework these parameters plug into
• Building an HMM-based Map Matcher with OSRM and Python — the runnable matcher class that consumes tuned sigma_z and beta values
• Choosing a Map-Matching Algorithm: HMM vs Geometric vs Viterbi — when the tuning cost in this guide is worth paying versus a simpler matcher
• Kalman Filtering for GPS Noise Reduction — upstream denoising that changes the sigma_z you should estimate downstream
• Speed Profiling from Raw GPS Coordinates — per-segment speed estimates feed the sampling-rate-aware beta heuristic above