Resampling, Likelihood and Bayes
Kernel density estimation
KDE turns a pile of data points into a smooth curve by stacking a small bump on each point — with one smoothing knob that decides what you see.
- 8 min read
- 3 reading levels
- Published
Read these first
On this page 4
One lesson, three depths. Pick the one that fits you today — you can switch any time.
Beginner — No maths. Plain English.
Kernel density estimation draws a smooth curve of where your data is dense, by piling a little mound of sand on every data point.
Drop one fistful of sand wherever each data point sits on a line. Where points crowd together, the mounds overlap into a hill. Where points are rare, the sand lies thin.
Step back and look at the sand's outline. That outline is the kernel density estimate: a smooth landscape of "how common are values around here?".
Why it exists
The histogram does this job crudely, and it cheats in two ways. Its bars jump when you nudge the bin edges — the same data can look single-humped or double-humped depending on where the bins fall. And it claims density changes in sudden steps, which nature rarely does.
KDE fixes both. Every point gets its own mound centred exactly on it — no bins, no edges, nothing to nudge. Statisticians Rosenblatt and Parzen formalised it in the late 1950s.
One honest choice remains, and it is the entire craft: how wide to make each mound. That width is called the bandwidth.
narrow mounds: ▲▲ ▲▲▲ ▲ ▲▲ spiky — every quirk shows,
noise included
wide mounds: ▁▂▄▅▄▂▁ one smooth blur — real
structure smoothed away
right width: ▂▅▂ ▁▄▂ two groups, honestly drawnToo narrow, and you draw the noise. Too wide, and you erase the truth. There are automatic rules for choosing, and they are decent starting points — not verdicts.
A real example you have seen
Population-density heat maps on Google Maps glow by this logic — each person a mound, cities become hills. Crime "hot-spot" maps in the news, wildlife-territory maps from sighting records, and the smooth violin-shaped distribution plots in data dashboards are all KDE at work.
Remember this
- KDE = a mound on every point; overlaps add up into a smooth curve.
- No bins, no bin-edge lottery — the histogram's two cheats removed.
- The bandwidth decides everything: too narrow draws noise, too wide erases structure.
What to learn next
- Q-Q plots and normality checks — comparing a shape against a reference family instead of drawing it freehand.
- The bootstrap — put uncertainty bands on your density curve.
- Matplotlib — draw the curves this lesson printed as text.
Developer — Code and libraries.
Setup
pip install numpy scipyOutputs verified with numpy 1.26 and scipy 1.14, CPU.
Two kinds of users, one bandwidth lesson
App session lengths in minutes: quick checkers and long bingers, almost nothing between. Can KDE see the gap?
import numpy as np
from scipy import stats
data = np.array([2, 3, 1, 4, 2, 3, 2, 45, 50, 38, 55, 42, 3, 2, 48])
default = stats.gaussian_kde(data) # bandwidth by Scott's rule
narrow = stats.gaussian_kde(data, bw_method=0.25)
print("minutes default narrow")
for x in [0, 3, 15, 25, 47, 60]:
d, n = default(x)[0], narrow(x)[0]
print(f"{x:7d} {d:.4f} {n:.4f} {'#' * int(n * 300)}")minutes default narrow
0 0.0179 0.0383 ###########
3 0.0183 0.0418 ############
15 0.0125 0.0037 #
25 0.0078 0.0004
47 0.0112 0.0194 #####
60 0.0071 0.0048 #The walkthrough
Read the narrow column's bars first. Two clear hills — around 3 minutes and around 47 — with a genuine valley at 15–25. That matches the data: checkers and bingers, nothing between.
Now read the default column. Density at 25 minutes: 0.0078 — nearly half its peak value, for a region containing zero data points. The automatic bandwidth (Scott's rule) assumes one gentle hump; handed two distant clusters, it pours sand into the valley. Automatic rules are starting points, not truths — when the picture matters, try a few bandwidths and watch which features persist.
bw_method=0.25 scales the automatic bandwidth down to a quarter. You can also pass a callable or a fixed value; the scaling form is the everyday one.
kde(x) returns density, not probability. Values are "probability per minute" — only areas under the curve are probabilities. Integrate for those:
print(f"P(session < 10 min) = {narrow.integrate_box_1d(-np.inf, 10):.2f}")P(session < 10 min) = 0.54
That 0.54 hides a confession: 9 of 15 sessions (60%) are under 10 minutes, but part of their sand spilled below zero, where sessions cannot exist. The next section's second mistake covers the cure.
KDE also generates. narrow.resample(5, seed=1) draws new synthetic sessions from the estimated shape — KDE is a legitimate (if humble) generative model, and the conceptual ancestor of fancier ones.
Common mistakes
Trusting the default bandwidth on multi-humped data. You watched it happen above. Scott's and Silverman's rules are derived assuming near-normal data; clusters break the assumption. Vary the bandwidth; report features that survive.
Reading density spillover as data. Gaussian mounds have infinite tails, so the estimate is positive at 0 minutes and even at -5 — where sessions cannot exist. For strictly positive quantities, estimate the density of log(data) and transform back, or use a boundary-corrected method.
Using KDE on ten points and presenting the curve with a straight face. The smooth line looks authoritative regardless of sample size. Below a few dozen points, show the points too, or bootstrap the curve to reveal its wobble.
Feeding KDE high-dimensional data. In more than a few dimensions, all points sit far apart and the sand never overlaps — the curse of dimensionality bites hard. Beyond 2–3 dimensions, use mixture models or normalising flows.
Try it yourself
Plot the full curves with matplotlib: evaluate both KDEs on np.linspace(-10, 70, 300) and overlay a histogram (density=True). Then find, by trial, the largest bandwidth scaling that still shows the valley at 20 minutes.
What to learn next
- Q-Q plots and normality checks — comparing a shape against a reference family instead of drawing it freehand.
- The bootstrap — put uncertainty bands on your density curve.
- Matplotlib — draw the curves this lesson printed as text.
Researcher — Mathematics and papers.
The estimator
Given i.i.d. samples $x_1, \dots, x_n$ from unknown density $f$:
$$ \hat{f}h(x) = \frac{1}{nh} \sum{i=1}^{n} K!\left(\frac{x - x_i}{h}\right) $$
Where:
- $K$ — the kernel: a symmetric density (Gaussian, Epanechnikov, …).
- $h$ — the bandwidth, the smoothing scale.
- $\hat f_h$ — a valid density: non-negative, integrating to 1.
Kernel choice is second-order: the Epanechnikov kernel $K(u) = \tfrac{3}{4}(1 - u^2)_+$ minimises asymptotic error, but the Gaussian loses only about 5% efficiency. Bandwidth choice is first-order — it is the estimator.
Bias–variance and the optimal rate
Taylor expansion gives the pointwise asymptotics:
$$ \operatorname{Bias} \approx \frac{h^2}{2} \sigma_K^2 f''(x), \qquad \operatorname{Var} \approx \frac{f(x) R(K)}{n h} $$
with $\sigma_K^2 = \int u^2 K(u)\,du$ and $R(K) = \int K^2(u)\,du$. Bias grows with $h$ (smoothing flattens curvature — precisely the valley-filling seen in the developer tab); variance shrinks with $h$. Minimising the integrated MSE yields
$$ h_{opt} = \left(\frac{R(K)}{\sigma_K^4\, R(f'')\, n}\right)^{1/5}, \qquad \text{IMSE} = O!\left(n^{-4/5}\right) $$
— slower than the parametric $n^{-1}$: the price of assuming nothing about shape. $h_{opt}$ depends on the unknown $R(f'') = \int f''(x)^2 dx$, hence the plug-in industry:
- Silverman's rule: $h = 0.9 \min(\hat\sigma, IQR/1.34)\, n^{-1/5}$ — exact for Gaussian truth, oversmooths multimodal data.
- Scott's rule: $h = \hat\sigma\, n^{-1/5}$ (scipy's default via
bw_method='scott', applied as a covariance scaling). - Sheather–Jones (1991): solve a fixed-point plug-in for $R(f'')$ — the accepted default among specialists.
- Least-squares / likelihood cross-validation: minimise estimated risk directly; LSCV is unbiased for IMSE but high-variance.
Multivariate form and its limits
$$ \hat f_H(x) = \frac{1}{n \left|H\right|^{1/2}} \sum_i K!\left(H^{-1/2}(x - x_i)\right) $$
with bandwidth matrix $H$. The optimal IMSE degrades to $O!\left(n^{-4/(4+d)}\right)$ in dimension $d$ — at $d = 10$, millions of samples buy what hundreds buy in 1-D. This rate is minimax over smoothness classes: not an algorithmic failure but an information-theoretic wall (Stone, 1982).
Boundary and shape corrections
Near a support boundary the kernel spills mass outside, biasing $\hat f$ downward by up to half. Remedies: reflection, boundary kernels, or transformation KDE ($\log$ for positive data). Adaptive/variable-bandwidth KDE ($h_i \propto f(x_i)^{-1/2}$, Abramson 1982) sharpens peaks while calming tails.
Connections
KDE with a Gaussian kernel is the predictive density of a Dirichlet-process-like mixture with one component per point; the mean-shift clustering algorithm ascends the KDE gradient; Parzen-window classifiers apply Bayes' rule to per-class KDEs; modern normalising flows and diffusion models occupy the same niche — density estimation — with learned, high-dimensional smoothers. Foundational papers: Rosenblatt (1956), Remarks on some nonparametric estimates of a density function; Parzen (1962), On estimation of a probability density function and mode.
What to learn next
- Q-Q plots and normality checks — comparing a shape against a reference family instead of drawing it freehand.
- The bootstrap — put uncertainty bands on your density curve.
- Matplotlib — draw the curves this lesson printed as text.