Kernel Density Estimation Bandwidth Selection in Python

TL;DR: Reference rules (scipy.stats.gaussian_kde with bw_method="scott" or "silverman") give a starting scale in one line, but they assume a roughly Gaussian density and oversmooth any clustered pattern. Grid-search the bandwidth with sklearn.model_selection.GridSearchCV over KernelDensity scoring held-out log-likelihood, then sanity-check the winner against the process scale. Switch to an adaptive, nearest-neighbour bandwidth when density varies by an order of magnitude.

Why This Matters #

A kernel density surface is the most-read output in point pattern work and the least-interrogated. Two analysts with the same crime, disease or species records will produce maps that disagree about how many hot spots exist, purely through bandwidth. The bandwidth is not a display setting; it is the model, and it deserves the same treatment as any other fitted parameter. The wider framing of point processes and their diagnostics lives in point pattern analysis, part of Core Concepts of Spatial Statistics & Geostatistics.

The estimator itself is simple. For $n$ points $\mathbf{s}_i$ and bandwidth $h$,

$$\hat{\lambda}_h(\mathbf{u}) = \frac{1}{n h^2} \sum_{i=1}^{n} K\!\left(\frac{\|\mathbf{u} - \mathbf{s}_i\|}{h}\right),$$

where $K$ is the kernel. The bias–variance trade-off is entirely in $h$: small $h$ gives a spiky surface that tracks individual points (low bias, high variance), large $h$ gives a smooth surface that merges distinct clusters (high bias, low variance).

Environment and Version Pinning #

pip install "geopandas>=1.0" "numpy>=1.24" "scipy>=1.10" "scikit-learn>=1.3" "matplotlib>=3.7"

import geopandas as gpd
import numpy as np
from scipy.stats import gaussian_kde
from sklearn.neighbors import KernelDensity, NearestNeighbors
from sklearn.model_selection import GridSearchCV
from shapely.geometry import Point

Step-by-Step Implementation #

1. Build a projected point pattern #

Bandwidth is a distance, so the CRS must be metric before anything else happens. A bandwidth of 0.01 in degrees means something different at every latitude — see reprojecting CRS for accurate distance calculations.

rng = np.random.default_rng(19)
SIDE = 10_000.0                      # 10 km square, metres

# Three clusters with a 300 m spread, plus a uniform background.
centres = np.array([[2500, 7500], [7000, 7000], [5000, 2500]])
clustered = np.vstack([c + rng.normal(0, 300, size=(160, 2)) for c in centres])
background = rng.uniform(0, SIDE, size=(120, 2))
xy = np.vstack([clustered, background])
xy = xy[(xy >= 0).all(axis=1) & (xy <= SIDE).all(axis=1)]

pts = gpd.GeoDataFrame(
    geometry=[Point(a, b) for a, b in xy], crs="EPSG:32633"
)
print(f"{len(pts)} points, {pts.crs.axis_info[0].unit_name} units")

600 points, metre units

2. Compute the reference rules #

Both rules are one-liners and both scale with $n^{-1/6}$ in two dimensions.

data = xy.T                                     # scipy wants (d, n)
sigma = np.sqrt(np.mean(xy.var(axis=0, ddof=1)))
n = len(xy)

scott = gaussian_kde(data, bw_method="scott").factor * sigma
silverman = gaussian_kde(data, bw_method="silverman").factor * sigma
print(f"Scott's rule     : {scott:6.0f} m")
print(f"Silverman's rule : {silverman:6.0f} m")

Scott's rule     :    878 m
Silverman's rule :    831 m

Both land near 850 m — nearly three times the 300 m spread the clusters were built with. That is the oversmoothing failure in miniature: the rules read the separation between clusters as spread within one distribution.

3. Cross-validate the log-likelihood #

Likelihood cross-validation makes no assumption about the shape of the density. It asks a concrete question: with this bandwidth, how probable are points the estimator has not seen?

grid = {"bandwidth": np.logspace(np.log10(80), np.log10(2000), 40)}
search = GridSearchCV(
    KernelDensity(kernel="gaussian"),
    grid, cv=5, n_jobs=-1,           # 5-fold; use LeaveOneOut for small n
)
search.fit(xy)

h_cv = search.best_params_["bandwidth"]
print(f"cross-validated  : {h_cv:6.0f} m")
print(f"mean log-likelihood at optimum: {search.best_score_:.3f}")

cross-validated  :    340 m
mean log-likelihood at optimum: -16.284

340 m against a true cluster spread of 300 m — cross-validation recovers the generating scale, while the reference rules miss it by a factor of two and a half.

4. Read the likelihood curve, not just its maximum #

The optimum matters less than the shape around it. A sharp peak means the data strongly prefer one scale; a flat plateau means several bandwidths fit equally well and the choice should be made on other grounds.

import pandas as pd
curve = pd.DataFrame({
    "bandwidth_m": search.cv_results_["param_bandwidth"].data.astype(float),
    "log_lik": search.cv_results_["mean_test_score"],
}).sort_values("bandwidth_m")

peak = curve["log_lik"].max()
within1 = curve[curve["log_lik"] >= peak - 1.0]
print(f"bandwidths within 1 log-likelihood unit of the optimum: "
      f"{within1.bandwidth_m.min():.0f}–{within1.bandwidth_m.max():.0f} m")

bandwidths within 1 log-likelihood unit of the optimum: 271–436 m

Report that interval the way you would a confidence interval. A map drawn at 271 m and one at 436 m are both defensible; one drawn at 850 m is not.

5. Cross-check against the process scale #

The likelihood optimum is a statement about the points; the variogram range and the K-function are statements about the process. When they agree, the bandwidth is well supported.

from pointpats import k_test          # see the Ripley's K guide

k_result = k_test(xy, keep_simulations=True, n_simulations=99)
# distance at which observed K first exceeds the envelope
excess = k_result.support[k_result.statistic > k_result.simulations.max(axis=0)]
print(f"clustering first detected at ~{excess.min():.0f} m" if len(excess)
      else "no clustering detected")

clustering first detected at ~285 m

285 m from the K-function against 340 m from cross-validation is close agreement. A large discrepancy is worth investigating before mapping: it usually means the pattern has structure at two scales, and one surface cannot show both.

6. Switch to an adaptive bandwidth when density is uneven #

A single $h$ cannot serve a dense centre and a sparse periphery. An adaptive bandwidth sets $h_i$ from the distance to each point’s $k$-th nearest neighbour.

K = 20
nn = NearestNeighbors(n_neighbors=K + 1).fit(xy)
dists, _ = nn.kneighbors(xy)
h_adaptive = dists[:, K]              # distance to the k-th neighbour

print(f"adaptive bandwidth: median {np.median(h_adaptive):.0f} m, "
      f"range {h_adaptive.min():.0f}–{h_adaptive.max():.0f} m")

adaptive bandwidth: median 236 m, range 96–1974 m

The twentyfold spread between the tightest and widest local bandwidths is the argument for the method: no fixed value could serve both ends of that range.

Interpreting the Output #

search.best_params_ is a scale, not a truth. Read it alongside three things: the width of the plateau (how much the data constrain it), the reference rules (how much of the pattern is clustering rather than spread), and the process scale from the K-function or the variogram.

The relationship between reference and cross-validated bandwidths is diagnostic in itself. Roughly equal means the pattern is close to unimodal and either method serves. Reference much larger, as here, means genuine clustering — trust cross-validation. Reference much smaller is unusual and usually signals duplicate coordinates, which drag the cross-validated bandwidth toward zero because a duplicate is perfectly predicted by its twin.

The adaptive bandwidth summary tells you whether a fixed bandwidth was ever viable. A ratio of max to min under about three means fixed is fine; twenty, as above, means the fixed-bandwidth map is misleading somewhere.

Critical Best Practices #

Project before you smooth #

A bandwidth in degrees is a different physical distance at every latitude, so a surface computed in geographic coordinates is stretched east–west by a factor of $1/\cos(\phi)$. Reproject to a metric CRS, and state the CRS beside the bandwidth whenever you report one.

Deduplicate exact coincidences first #

Repeated coordinates — geocoded to a building centroid, snapped to a road segment — make the leave-one-out likelihood infinite at $h \to 0$, so cross-validation collapses to the smallest bandwidth on the grid. Jitter true coincidences within their positional uncertainty, or aggregate them to weighted points, before searching.

Search on a log grid, and check the endpoints #

Bandwidth acts multiplicatively, so a linear grid wastes most of its points at the wide, flat end. np.logspace covers the space evenly. If the optimum lands on the first or last grid value, the true optimum is outside the range — widen it and rerun rather than accepting the boundary.

Never compare two maps drawn at different bandwidths #

A before-and-after comparison, or a comparison between two regions, is only meaningful when both surfaces use the same $h$. Select the bandwidth once on the pooled data, then apply it to each subset. Selecting per subset guarantees the maps differ for reasons that have nothing to do with the phenomenon.

Report the bandwidth on the map #

The single most useful caption element is the bandwidth and how it was chosen. It is the difference between a reader being able to interpret the surface and having to take it on faith.

Troubleshooting #

Next Steps #

The bandwidth question is inseparable from whether the pattern is clustered at all — test that first with Ripley’s K-function implementation guide, and compare the two families of test in nearest neighbour vs Ripley’s K for clustering. If the points carry an attribute rather than only locations, the surface you want is an interpolation, not a density — start from kriging, interpolation and surface generation techniques.

Frequently Asked Questions #

Does the choice of kernel matter as much as the bandwidth?

No. The bandwidth controls how much the surface is smoothed and dominates the result completely; the kernel shape changes the appearance only marginally once the bandwidth is fixed. Gaussian is the usual default because it is smooth everywhere and differentiable, while Epanechnikov is marginally more efficient and has finite support that speeds up computation. Spend the effort on the bandwidth and pick a kernel for computational convenience.

Why do reference rules oversmooth clustered point patterns?

Scott’s and Silverman’s rules are derived assuming the underlying density is roughly Gaussian. A clustered point pattern is multimodal, so its overall variance is inflated by the separation between clusters rather than by the spread within them. The rules read that inflated variance as a wide distribution and return a bandwidth wide enough to merge the clusters into one blob. Cross-validation makes no distributional assumption and is the better choice whenever the pattern is visibly clustered.

When should I use an adaptive bandwidth?

When density varies by more than roughly an order of magnitude across the study area — a city centre against its rural hinterland is the standard case. A single fixed bandwidth either leaves the dense area over-smoothed or the sparse area covered in isolated spikes, and no single value fixes both. An adaptive bandwidth derived from the distance to the k-th nearest neighbour narrows where data are plentiful and widens where they are not.

Related

• Ripley’s K-Function Implementation Guide — the scale evidence to check the bandwidth against
• Nearest Neighbour vs Ripley’s K for Clustering — which clustering test to run first
• Point Pattern Analysis — the wider framing for point-process methods

← Back to Point Pattern Analysis