Calibrating Distance-Decay Functions for Trade Areas

This page solves one exact task: fitting a distance-decay function to observed customer or loyalty-card data so the reachable-population term in a suitability score reflects how patronage actually falls off with distance, rather than assuming every reachable customer counts equally.

A raw drive-time catchment treats a household 2 minutes away and one 14 minutes away as identical members of the reachable population. Real patronage does not work that way: the probability a customer shops at a site decays with travel cost, and the shape of that decay is empirically estimable from data you already hold. Calibrating it converts a flat catchment count into a distance-weighted expected-customer count, which is the honest input to building weighted site suitability scores.

Concept: exponential and power decay #

Two functional forms dominate retail gravity modeling. The exponential form assumes patronage falls at a constant proportional rate per unit distance:

$$f(d) = e^{-\lambda d}$$

where $d$ is travel distance (or drive time) and $\lambda > 0$ is the decay rate — larger $\lambda$ means a tighter, more local trade area. The power form assumes decay is steep near the store and flattens with distance:

$$f(d) = d^{-\beta}$$

with exponent $\beta > 0$. Exponential decay generally fits convenience and grocery formats where customers refuse long trips; power decay often fits destination or comparison retail where shoppers tolerate distance for selection. The calibration task is to estimate $\lambda$ (or $\beta$) from observed data and pick the form that fits best.

The estimated decay function then multiplies each demographic increment by its distance weight, so the decay-adjusted reachable population is $\sum_k P_k \, f(d_k)$ over block groups $k$ at distance $d_k$ — the term that feeds the scoring model instead of a raw sum. This is the same $f_d$ that appears in the accessibility score discussed across the Isochrone Generation & Network Analysis guides.

A third form, the exponential-power (or gamma) hybrid $f(d) = d^{-\beta} e^{-\lambda d}$, is worth fitting when neither pure form captures the data: it decays steeply near the store like a power law but suppresses the far tail like an exponential. It costs a second free parameter and needs more observations to fit stably, so reach for it only when the residuals of both simpler forms show systematic curvature. For most grocery and convenience formats the plain exponential wins on parsimony; for destination retail the power form usually edges it. Whichever wins, the decision should be made on held-out residuals, not on which curve looks prettier over the fitted range.

Prerequisites #

Before running this task you need:

• Python packages: numpy and scipy for the curve fit, pandas for the observation table, and geopandas / pyproj to compute distances in a projected CRS. Install with pip install numpy scipy pandas geopandas pyproj.
• Observed patronage-by-distance data. A loyalty-card or point-of-sale extract giving, for an existing store, the share of customers (or trip counts normalized to a share) originating in each distance band. The method needs the distance of each origin and its observed patronage.
• Store and customer-origin coordinates in a known CRS. Distances must be computed in a projected, distance-preserving CRS — this task uses EPSG:5070 (NAD83 / Conus Albers) — never in raw WGS84 degrees.

Configuration and execution parameters #

The distance metric is the load-bearing choice: if the suitability model weights population by drive time, calibrate decay against drive time, not straight-line km. Mixing metrics produces a decay curve that is internally consistent but wrong for the pipeline that consumes it.

Sample size per band matters as much as the number of bands. A distance band built from a handful of loyalty customers carries a noisy patronage share that the optimizer will nonetheless try to fit exactly. When band sizes vary widely, pass the per-band standard error to curve_fit via its sigma argument so the fit weights well-populated bands more heavily and the covariance it returns reflects real uncertainty. Bands with fewer than a minimum customer count are better dropped than trusted, because a single atypical shopper can swing a thin band’s share by tens of percent and drag the whole curve with it.

Annotated implementation #

The function below computes network-consistent distances in EPSG:5070, fits both an exponential and a power decay to the observed patronage shares with scipy.optimize.curve_fit, and returns the better-fitting model with its estimated parameter and goodness-of-fit statistics.

import numpy as np
import pandas as pd
import geopandas as gpd
from pyproj import CRS
from scipy.optimize import curve_fit

DIST_CRS = CRS.from_epsg(5070)   # NAD83 / Conus Albers: distances in meters


def exp_decay(d, lam):
    """Exponential distance decay f(d) = exp(-lambda * d)."""
    return np.exp(-lam * d)


def power_decay(d, beta):
    """Power distance decay f(d) = d^(-beta), guarded near zero."""
    return np.power(np.maximum(d, 1e-6), -beta)


def observed_distances(store: gpd.GeoDataFrame,
                       origins: gpd.GeoDataFrame) -> np.ndarray:
    """Distance (km) from a store to each customer origin, in a projected CRS."""
    if store.crs is None or origins.crs is None:
        raise ValueError("store/origins missing CRS; refusing to measure distance")
    store_p = store.to_crs(DIST_CRS)
    origins_p = origins.to_crs(DIST_CRS)
    assert CRS.from_user_input(store_p.crs) == DIST_CRS   # never measure in degrees
    pt = store_p.geometry.iloc[0]
    return origins_p.geometry.distance(pt).to_numpy() / 1000.0   # meters -> km


def calibrate_decay(distance_km: np.ndarray,
                    patronage: np.ndarray) -> dict:
    """
    Fit exponential and power decay to observed (distance, patronage-share) data.

    `patronage` should be a normalized share in (0, 1], highest near the store.
    Returns the better-fitting model, its parameter, and R^2.
    """
    d = np.asarray(distance_km, dtype=float)
    y = np.asarray(patronage, dtype=float)
    if d.size < 8:
        raise ValueError(f"only {d.size} observations; need >= 8 for a stable fit")
    if np.isnan(d).any() or np.isnan(y).any():
        raise ValueError("NaN in distance or patronage input")

    results = {}
    for name, fn, p0 in (("exponential", exp_decay, [1.0]),
                         ("power", power_decay, [1.0])):
        params, cov = curve_fit(fn, d, y, p0=p0,
                                bounds=(0, np.inf), maxfev=10000)
        resid = y - fn(d, *params)
        ss_res = float(np.sum(resid ** 2))
        ss_tot = float(np.sum((y - y.mean()) ** 2))
        r2 = 1.0 - ss_res / ss_tot if ss_tot > 0 else float("nan")
        stderr = float(np.sqrt(np.diag(cov))[0])   # 1-sigma on the parameter
        results[name] = {"param": float(params[0]), "stderr": stderr,
                         "ss_res": ss_res, "r2": r2, "fn": fn}

    best = min(results, key=lambda k: results[k]["ss_res"])
    return {"best_form": best, "models": results}

Driving the calibration on a loyalty extract is then direct:

obs = pd.DataFrame({
    "distance_km": [0.6, 1.2, 2.1, 3.4, 4.8, 6.5, 8.1, 10.3, 13.0],
    "share":       [0.98, 0.82, 0.63, 0.44, 0.31, 0.19, 0.12, 0.07, 0.03],
})
fit = calibrate_decay(obs["distance_km"].to_numpy(), obs["share"].to_numpy())

best = fit["best_form"]
model = fit["models"][best]
print(f"best form: {best}")
print(f"param (lambda/beta) = {model['param']:.4f} +/- {model['stderr']:.4f}")
print(f"R^2 = {model['r2']:.3f}")

With curve_fit returning both the parameter and its covariance, the reported standard error is a genuine uncertainty band on $\lambda$, not a point estimate presented as fact — essential when the decay rate propagates into a capital decision.

Applying the fitted decay #

Once $\lambda$ is estimated, weight each demographic increment by its decay factor to produce the decay-adjusted reachable population that enters the scoring matrix:

def decay_weighted_population(blockgroups: gpd.GeoDataFrame,
                              store: gpd.GeoDataFrame,
                              lam: float) -> float:
    """Sum block-group population weighted by exp(-lambda * distance)."""
    d_km = observed_distances(store, blockgroups.centroid.to_frame("geometry")
                              .set_crs(blockgroups.crs))
    weights = exp_decay(d_km, lam)
    return float((blockgroups["population"].to_numpy() * weights).sum())

This decay-weighted count replaces the flat reach_pop sum from a raw catchment, so a captive nearby household contributes nearly its full weight while a distant one contributes a small fraction — the behavior the suitability model needs.

Failure modes and debugging #

The most consequential silent error is measuring distance in WGS84 degrees. A degree of longitude is not a fixed distance, so a decay fit on degree “distances” produces a $\lambda$ that is wrong and location-dependent. The observed_distances function guards against this by asserting the projected CRS before any distance is taken.

Verification #

Confirm the calibration before letting $\lambda$ drive a score:

1. Positivity and plausibility: assert model["param"] > 0; a fitted decay rate must be strictly positive, and for grocery formats $\lambda$ typically implies patronage halving every few km.
2. Goodness of fit: require R² above a threshold (0.85 is a reasonable bar for clean loyalty data); a low R² means the chosen form is wrong or the data is too noisy to calibrate.
3. Residual structure: plot or test residuals for systematic curvature. Residuals that arc consistently above then below the fit signal the wrong functional form — switch exponential to power or vice versa.
4. Monotonicity of weights: assert (np.diff(exp_decay(sorted_d, lam)) <= 0).all() — decay weights must never increase with distance.

sorted_d = np.sort(obs["distance_km"].to_numpy())
w = exp_decay(sorted_d, model["param"])
assert (np.diff(w) <= 1e-12).all(), "decay weights must be non-increasing"
assert model["r2"] > 0.85, f"weak fit: R^2={model['r2']:.3f}"

A calibration that passes these checks yields a decay curve you can defend to an investment committee: an empirically fitted, uncertainty-quantified weighting that ties the reachable-population term to observed customer behavior rather than an assumed radius.

Related #

• Building Weighted Site Suitability Scores — the parent stage that consumes the decay-weighted population term.
• Normalizing Mixed-Scale Site Attributes for Scoring — put the decay-weighted count on the same scale as other criteria.
• Weighting Demographic Variables for Target Audiences — build the demographic term that decay then discounts by distance.

← Back to Building Weighted Site Suitability Scores