Location Intelligence Retail Site Selection Automation

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−λdf(d) = e^{-\lambda d}

where dd is travel distance (or drive time) and λ>0\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−βf(d) = d^{-\beta}

with exponent β>0\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 ∑kPk f(dk)\sum_k P_k \, f(d_k) over block groups kk at distance dkd_k — the term that feeds the scoring model instead of a raw sum. This is the same fdf_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−βe−λdf(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

Parameter Value for this task Notes
Decay form exponential / power Fit both; select by lower residual sum of squares.
Distance metric drive time (min) or network km Must match the metric used downstream in scoring.
Distance CRS EPSG:5070 Equal-area / distance-preserving for CONUS; assert before any distance math.
p0 (initial guess) [1.0] for λ\lambda curve_fit needs a starting point; 1.0 is a safe scale for km/min.
bounds (0, inf) Decay rate must be strictly positive.
Min observations ≥ 8 distance bands Fewer points give an unstable, over-confident fit.

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.

python
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:

python
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.

Fitting the decay to observed patronage rather than assuming it Loyalty-card trips binned by drive time give an observed patronage share that falls from about 0.62 at two minutes to 0.04 at twenty-eight minutes. A fitted exponential with lambda 0.11 tracks the bins closely to twenty minutes and then sits slightly above them, which is the signature of a tail better described by a power law. The decay parameter is measured, not chosen 0.7 0.5 0.3 0.1 0 8 16 24 drive time to the store (minutes) observed patronage share fitted curve: exp(−0.11 t), R² = 0.97 bars are observed trips per drive-time bin Fit on the log of the share so the regression is linear; check the tail separately, where an exponential systematically over-predicts the long trips that a power law captures.

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:

python
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.

The same catchment, three decay rates, three different demand estimates For a catchment holding 18,000 people within five minutes, 41,000 between five and fifteen, and 76,000 between fifteen and thirty, a slow decay of 0.05 credits 84,900 effective people, a fitted 0.11 credits 52,600, and a fast 0.20 credits 31,800. The ranking of two candidate sites can invert between the first and the last. One catchment, three decay rates, a 2.7× spread in credited demand band people λ = 0.05 λ = 0.11 (fitted) λ = 0.20 0–5 min 18,000 16,000 13,900 11,000 5–15 min 41,000 25,800 13,700 5,600 15–30 min 76,000 43,100 25,000 15,200 effective demand 135,000 84,900 52,600 31,800 Why this decides rankings A site whose population sits mostly in the far band loses two thirds of its credited demand at λ = 0.20 and very little at λ = 0.05 — so an unfitted λ silently picks winners by geography rather than by evidence. Band weights use the midpoint travel time: 2.5, 10 and 22.5 minutes.

Failure modes and debugging

Symptom Cause Fix
curve_fit raises RuntimeError Optimizer did not converge in maxfev iterations Raise maxfev, supply a better p0, or rescale distance to a saner unit.
Negative or absurd lambda Missing bounds, or patronage not monotonically decreasing Enforce bounds=(0, inf); check the data is a decay, not noise.
R² near zero Wrong functional form, or distance measured in degrees Fit both forms; confirm distances came from a projected CRS.
Fit dominated by one point An outlier band with few underlying customers Weight the fit by sample size (sigma= in curve_fit) or drop thin bands.
Decay too flat / too steep vs. reality Straight-line distance used where drive time was needed Recompute distances against the drive-time metric used downstream.

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.

Choosing the decay family on held-out trips, not on the fit Exponential, power-law and Gaussian decays are compared on in-sample fit and on trips held out from the calibration. The exponential fits the near bands best, the power law wins on the held-out long trips, and the Gaussian looks strong in sample but degrades most out of sample. Pick the family on held-out trips — in-sample fit always flatters family in-sample R² held-out error where it goes wrong exponential 0.97 0.041 over-predicts the long tail power law 0.95 0.028 unstable below 2 minutes Gaussian 0.98 0.067 too flat near the store The Gaussian has the best in-sample R² and the worst held-out error — the classic shape of a family with enough freedom to trace the calibration bins rather than the behaviour underneath them. A defensible calibration therefore reports both columns and the split rule that produced them, so the choice can be re-examined when a new season of trip data arrives. Held-out error is mean absolute deviation in patronage share on trips excluded from the fit.

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.
python
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.

Frequently Asked Questions

How much trip data is enough to calibrate a decay?

Enough that every drive-time band has trips in it, which usually matters more than the total. A hundred thousand trips concentrated within ten minutes of the store constrains the near bands beautifully and says nothing about the tail, so the fitted parameter is decided by whichever handful of long trips happened to land in the sample. Check the per-band counts before trusting a fit, and widen the bands at the far end rather than fitting a curve through bins holding single-digit observations.

Can one decay parameter serve the whole chain?

Rarely. Decay is a property of the trip purpose and the format, not of the company, so a dense urban store and a highway-adjacent one in the same chain genuinely behave differently. Fitting per format and per market type — and falling back to a pooled parameter only where the local sample is too thin — keeps the model honest, and the fallback should be recorded on the site so a reviewer can see which sites were scored on borrowed evidence.

What if no trip data exists for a new market?

Borrow from the most similar market you have measured, and mark the borrowing explicitly. A parameter fitted in a comparable metro is far better than an invented one, and it comes with a provenance note that keeps its uncertainty visible. Then plan the measurement: once the site opens, the loyalty data from its own catchment is the calibration set that replaces the borrowed value, and the recalibration should be a scheduled job rather than an act of remembering.

Should the decay be fitted on distance or on travel time?

On travel time, because that is what shoppers actually experience and what the routing layer already computes. Distance and time diverge exactly where the decision matters — a five-kilometre trip across a river crossing and a five-kilometre trip on an arterial are not the same trip — and a decay fitted on distance quietly absorbs the network into a parameter that then cannot transfer to another market. Fit on time, and let the routing engine own the geography.

← Back to Building Weighted Site Suitability Scores