Applying Huff Model Probabilities to Store Catchments

This page solves one task: computing the Huff-model probability that a customer at each demographic origin patronizes each candidate store, then rolling those probabilities up into expected captured demand that a ranking model can consume.

The Huff model is the workhorse gravity formulation of retail catchment analysis. It replaces the fiction of hard trade-area boundaries with a continuous probability field: a customer is more likely to shop at a store that is larger, closer, or both, and that pull decays with distance. The output — a probability per origin-store pair — is exactly the capture input that ranking and shortlisting candidate sites folds into its Huff-adjusted rank score.

The model #

The probability that a customer at origin $i$ patronizes store $j$ is that store’s attractiveness-over-distance pull, divided by the summed pull of every store the customer could choose:

$$P_{ij} = \frac{A_j \, d_{ij}^{-\beta}}{\sum_{k} A_k \, d_{ik}^{-\beta}}$$

Here $A_j$ is the attractiveness of store $j$ (commonly floor area, but any positive utility proxy works), $d_{ij}$ is the distance from origin $i$ to store $j$, and $\beta > 0$ is the distance-decay exponent that controls how sharply pull falls with distance. The denominator sums over the full competitive set $k$, which is what makes each $P_{ij}$ a proper probability: for any origin $i$, the probabilities over all stores sum to exactly 1.

Expected captured demand for store $j$ is then the probability-weighted sum of demand across origins:

$$D_j = \sum_{i} P_{ij} \, C_i$$

where $C_i$ is the buying power at origin $i$ — population, households, or dollar demand. The decay exponent $\beta$ is the parameter that most changes the result; a grocery trip has a steeper decay (customers will not drive far for milk) than a furniture showroom. Calibrating it belongs with the broader work of calibrating distance-decay functions for trade areas; this page takes $\beta$ as a given input.

Two modeling choices shape the result as much as $\beta$ does. The first is the attractiveness proxy $A_j$: floor area is the conventional default because larger stores carry more assortment, but any positive utility measure works — anchor-tenant count, parking capacity, or a blended index. Because attractiveness enters linearly in both numerator and denominator, only the ratio of one store’s attractiveness to its competitors’ matters, so consistent units across the store layer are what count, not absolute magnitude. The second is the competitive set: the denominator must include every store a customer at origin $i$ would realistically consider. Truncating that set to a fixed radius is a common shortcut, but it biases probabilities upward near the radius edge, so prefer a distance cutoff generous enough that omitted stores contribute negligible pull.

A useful way to read $D_j$ is as an expected market share: dividing a store’s captured demand by total demand in the study area gives the share of buying power the gravity field assigns to it. That share, computed for a candidate against the incumbent competitive field, is precisely the quantity that separates a well-located but contested site from one that actually wins its trade area.

Prerequisites #

• Python packages: geopandas, pandas, numpy, shapely, and pyproj. Install with pip install geopandas pandas numpy shapely pyproj.
• A store layer as a GeoDataFrame of points with an attractiveness column (floor area in m² is the usual proxy). Include both candidate and incumbent competitor stores — omitting competitors inflates every capture probability.
• A demographic origin layer of block-group centroids with a demand column (population or households), joined from ACS via how to join ACS 5-year estimates to custom trade area polygons. Use centroids, not polygons, so distance is point-to-point.
• A calibrated decay exponent $\beta$. Typical retail values run 1.5–2.2; start at 2.0 if you have no local calibration.

Configuration and execution parameters #

The distance floor matters: raw $d_{ij}^{-\beta}$ diverges as $d \to 0$, so an origin centroid coincident with a store produces an infinite pull. Clamping $d_{ij}$ to a small positive minimum keeps the ratio finite and well-behaved.

Annotated implementation #

The function computes a full origin-by-store probability matrix, then aggregates expected demand per store. All distance math happens in EPSG:5070, and the CRS is asserted before any geometric operation.

import geopandas as gpd
import numpy as np
import pandas as pd
from pyproj import CRS

EQUAL_AREA_CRS = CRS.from_epsg(5070)   # NAD83 / Conus Albers, meters


def huff_probabilities(
    origins: gpd.GeoDataFrame,
    stores: gpd.GeoDataFrame,
    beta: float = 2.0,
    attractiveness_col: str = "floor_area_m2",
    demand_col: str = "population",
    min_distance_m: float = 100.0,
) -> tuple[np.ndarray, pd.DataFrame]:
    """
    Compute Huff patronage probabilities for every origin-store pair.

    Returns:
        P     : (n_origins, n_stores) probability matrix; each row sums to 1.
        demand: per-store expected captured demand, indexed like `stores`.
    """
    # Distances are only meaningful in a projected, equal-area CRS.
    assert origins.crs is not None, "origins missing CRS"
    assert stores.crs is not None, "stores missing CRS"
    o = origins.to_crs(EQUAL_AREA_CRS)
    s = stores.to_crs(EQUAL_AREA_CRS)

    if (s[attractiveness_col] <= 0).any():
        raise ValueError("attractiveness must be strictly positive")

    # Coordinate arrays for a vectorized pairwise distance matrix.
    ox = o.geometry.x.to_numpy()[:, None]   # (n_origins, 1)
    oy = o.geometry.y.to_numpy()[:, None]
    sx = s.geometry.x.to_numpy()[None, :]   # (1, n_stores)
    sy = s.geometry.y.to_numpy()[None, :]

    dist = np.hypot(ox - sx, oy - sy)       # meters, shape (n_origins, n_stores)
    dist = np.maximum(dist, min_distance_m)  # clamp so d^-beta stays finite

    attract = s[attractiveness_col].to_numpy()[None, :]
    pull = attract * np.power(dist, -beta)   # A_j * d_ij^-beta

    # Normalize each origin row so probabilities sum to 1 across all stores.
    row_totals = pull.sum(axis=1, keepdims=True)
    P = pull / row_totals

    # Expected captured demand: P_ij weighted by origin demand, summed over i.
    demand_i = o[demand_col].to_numpy()[:, None]
    captured = (P * demand_i).sum(axis=0)    # shape (n_stores,)
    demand = pd.DataFrame({
        "store_id": stores.index,
        "expected_demand": captured,
    }).set_index("store_id")

    return P, demand

A candidate site’s mean capture probability — the scalar the ranking stage consumes — is the demand-weighted average of its column in P:

def mean_capture_for_store(P, origins, store_col_idx, demand_col="population"):
    """Demand-weighted mean patronage probability for one store column."""
    w = origins[demand_col].to_numpy()
    return float(np.average(P[:, store_col_idx], weights=w))

Failure modes and debugging #

Omitting competitor stores is the most consequential error: because the denominator sums over the full choice set, dropping incumbents makes every candidate look like it captures the whole market. The competitor footprint that populates the store layer is the subject of competitor mapping and cannibalization analysis, and the sales a new site pulls from existing ones follows in measuring sales transfer between nearby stores.

Verification #

Confirm the model is well-formed before trusting its output:

1. Row-stochastic check: every origin row of P sums to 1 within floating-point tolerance.
2. Bounded probabilities: all entries lie in [0, 1].
3. Demand conservation: summed expected demand across stores equals summed origin demand.

assert np.allclose(P.sum(axis=1), 1.0), "rows must sum to 1"
assert (P >= 0).all() and (P <= 1).all(), "probabilities out of range"
total_in = origins["population"].sum()
total_out = demand["expected_demand"].sum()
assert np.isclose(total_in, total_out), "demand not conserved"

The conservation check is the strongest single test: because the probabilities partition each origin’s demand across the full store set, every unit of buying power must land somewhere. If the totals disagree, an origin has zero total pull (all distances or attractiveness invalid) and its row silently dropped out of the normalization.

Related #

• Ranking & Shortlisting Candidate Sites — consumes these capture probabilities in the Huff-adjusted rank score.
• Measuring Sales Transfer Between Nearby Stores — quantifies demand a new store pulls from incumbents.
• How to Join ACS 5-Year Estimates to Custom Trade Area Polygons — builds the demographic origin layer.

← Back to Ranking & Shortlisting Candidate Sites