B4B33RPZ · Pattern Recognition and Machine Learning

Parzen windows

Lab 4: estimating a probability density from samples: histogram, k-NN and Parzen windows (kernel density estimation)

Klára JanouškováFEE CTU, Prague

Part I

Recap

What have we assumed so far?

	$p(x \mid k)$	priors $p(k)$	loss $W$
Bayesian decision	known	known	known
Neyman–Pearson, minimax, Wald	known	unknown	unknown
Today	unknown	we only have samples $x_1, \dots, x_n$

Probability density estimation: from samples

x_1, \dots, x_n

drawn from

p(x)

, build an estimate

\hat p(x)

. For a classifier we do it once per class: the samples of class

k

give

\hat p(x \mid k)

.

How can we estimate a density?

Parametric (next lab)

Assume a shape, e.g. a Gaussian; estimate its few parameters ( $\mu$ , $\sigma$ ).

MLE, MAP, Bayesian estimation

Non-parametric (today)

No fixed shape: the samples themselves define $\hat p$ .

histogram, k-NN, Parzen windows = kernel density estimation

Mixtures

A sum of a few simple densities: in between the two.

mixture of Gaussians, deep neural networks

Parametric: few numbers to learn, but if the true density has a different shape (two bumps, one Gaussian fitted), the error stays however many samples we get. Non-parametric: follows any shape given enough data, but needs more samples and keeps all of them to evaluate $\hat p(x)$ .

Part II

Histogram

We have samples. What is $p(x)$ ?

We have: $\{x_1, x_2, \dots, x_n\}$ , each $x_i \sim p(x)$ , independent

We want: $p(x)$

Where is $p(x)$ high? Where is it zero? How would you decide?

Count: many samples close together → $p(x)$ high; no samples → $p(x) \approx 0$ .

Simplest way to count: cut the axis into bins and count the samples in each bin.

Counts as bars = a histogram. But a count is not a density: how tall should a bar be?

How tall should a histogram bar be?

The same bins as before, now as bars. A bin $B$ holds $N_B$ of the $N$ samples and has width $w_B$ . Which height $\hat p$ ?

↓ what the data say: the fraction of samples that landed in the bin

\text{fraction in } B = \frac{\textcolor{blue}{N_B}}{N}

The bin has width $w_B$ .

↓ what our bar says: probability of the bin = area of the bar = height × width

\text{area of the bar} = \textcolor{purple}{\hat p(x)}\, \textcolor{orange}{w_B}

↓ we choose the height so that the two agree

\textcolor{purple}{\hat p(x)}\, \textcolor{orange}{w_B} = \frac{\textcolor{blue}{N_B}}{N}

↓ divide both sides by $\textcolor{orange}{w_B}$

\textcolor{purple}{\hat p(x)} = \frac{\textcolor{blue}{N_B}}{N\, \textcolor{orange}{w_B}}

In words: take each bin's relative frequency

\frac{\textcolor{blue}{N_B}}{N}

and divide it by the bin width

\textcolor{orange}{w_B}

. Bin

[-2, -1)

:

\frac{\textcolor{blue}{13}}{40} \div \textcolor{orange}{1} = \textcolor{purple}{0.325}

Histogram: pros and cons?

Discuss in pairs for a minute, then we collect.

✓ Pros

✓ simple, easy to visualise
✓ we don't need to keep the $x_i$ : the bin counts are enough
✓ fast to evaluate: find the bin, read the height

✗ Cons

✗ step-like, not smooth
✗ we must choose $w_B$ : too small → empty bins, too large → no detail
✗ depends on where the bins start
✗ sparse for dimension $> 2$ : $m$ bins per axis → $m^D$ cells

Try it: bin width and bin origin

bin width

w_B

bins start at show the true

p(x)

Part III

More dimensions: $K$ , $N$ and $V$

More dimensions?

1D: $\hat p = \textcolor{blue}{N_B} / (N\, \textcolor{orange}{w_B})$ . In 2D the bin is a square region $R$ around $x$ . Rename: count $\textcolor{blue}{N_B} \to \textcolor{blue}{K}$ , size $\textcolor{orange}{w_B} \to \textcolor{orange}{V}$ (area $s^2$ ; a volume in $D$ dims).

↓ what our estimate says: $\textcolor{purple}{\hat p}$ is constant over $R$ , so probability = height × volume

\text{probability of } R = \textcolor{purple}{\hat p(x)}\, \textcolor{orange}{V}

↓ what the data say: $\textcolor{blue}{K}$ of the $N$ samples fell into $R$

\text{fraction in } R = \frac{\textcolor{blue}{K}}{N}

↓ we choose $\textcolor{purple}{\hat p}$ so that the two agree, then divide by $\textcolor{orange}{V}$

\textcolor{purple}{\hat p(x)}\, \textcolor{orange}{V} = \frac{\textcolor{blue}{K}}{N} \;\Rightarrow\; \textcolor{purple}{\hat p(x)} = \frac{\textcolor{blue}{K}}{N\, \textcolor{orange}{V}}

Two ways to use it (

N

fixed):
fix $V$ , count $K$ : histogram, Parzen windows
fix $K$ , find $V$ : k-nearest neighbours

From now on $R$ is a ball around $x$ (a circle in 2D): one number, the radius $\textcolor{orange}{r}$ ; $\textcolor{orange}{V = \pi r^2}$ . k-NN grows $\textcolor{orange}{r}$ until $\textcolor{blue}{K}$ samples are inside.

Fix $K$ , find $V$ : k-nearest neighbours

$R$ is now a circle around the query. $N = 21$ samples, $K = 3$ . Question: is $\hat p$ higher at q1 or at q2?

q1: grow a circle until it holds 3 points. The 3rd nearest point sets the radius: $r = 2.42$

q2: the 3 nearest are much closer: $r = 1.12$

↓ in 2D the "volume" of the circle is its area

\textcolor{purple}{\hat p(q)} = \frac{\textcolor{blue}{K}}{N\, \textcolor{orange}{V}}, \qquad \textcolor{orange}{V = \pi r^2}

q1:

\textcolor{purple}{\hat p} = \frac{\textcolor{blue}{3}}{21 \cdot \textcolor{orange}{\pi \cdot 2.42^2}} = \frac{\textcolor{blue}{3}}{21 \cdot \textcolor{orange}{18.3}} = \textcolor{purple}{0.0078}

q2:

\textcolor{purple}{\hat p} = \frac{\textcolor{blue}{3}}{21 \cdot \textcolor{orange}{\pi \cdot 1.12^2}} = \frac{\textcolor{blue}{3}}{21 \cdot \textcolor{orange}{3.9}} = \textcolor{purple}{0.0364}

Same

K

, smaller circle

V

→ higher density

\hat p

at q2.

More: k-nearest neighbours

Fix $V$ , count $K$ : a window

Same 21 samples. Now every query gets a circle of the same radius $r = 2$ . Count the points inside.

q1: $K = 1$

q2: $K = 5$

↓ now $\textcolor{orange}{V}$ is the same for both queries, only $\textcolor{blue}{K}$ changes

\textcolor{purple}{\hat p(q)} = \frac{\textcolor{blue}{K}}{N\, \textcolor{orange}{V}}, \qquad \textcolor{orange}{V = \pi \cdot 2^2 = 12.6}

q1:

\textcolor{purple}{\hat p} = \dfrac{\textcolor{blue}{1}}{21 \cdot \textcolor{orange}{12.6}} = \textcolor{purple}{0.0038}

q2:

\textcolor{purple}{\hat p} = \dfrac{\textcolor{blue}{5}}{21 \cdot \textcolor{orange}{12.6}} = \textcolor{purple}{0.0189}

Both ways agree: q2 lies in the denser region. A window around

x

that counts points is the idea behind Parzen.

Quick check: fixed $K$ or fixed $V$ ?

We have $N \ge K$ samples. With which estimate can some $x$ get $\hat p(x) = 0$ ?

$\hat p = K/(NV)$ . A fixed window in an empty area holds $K = 0$ , so $\hat p = 0$ . A k-NN circle grows until it holds $K$ samples, so $\hat p > 0$ everywhere (with fewer than $K$ samples k-NN is not defined at all).

In the k-NN estimate, a query far from all points gets…

$K$ is fixed, so $V$ grows with the distance: $\hat p \sim 1/V$ is small but never 0.

k-NN density in 1D: what does $K$ do?

$K = 1$ : a spike wherever two samples are close: overfits.

$K = 5$ : follows both bumps, still noisy.

$K = 20$ : smooth, but flat and with heavy tails.

Notice the tails: far from the data the k-NN estimate stays much too high. Why, and what it means: next slide.

What happens far from the data?

A query right of all samples, $x = 6$ . Its $\textcolor{blue}{K = 5}$ nearest samples are the rightmost ones.

↓ in 1D the "ball" is an interval of half-width $\textcolor{orange}{r}$ ; it must reach back to the data, so $\textcolor{orange}{r}$ grows with the distance from the data

\textcolor{orange}{V} = 2\textcolor{orange}{r}, \quad \textcolor{orange}{r} = x - 2.21 \;\approx\; |x| \text{ for large } |x|

here $\textcolor{orange}{r} = 6 - 2.21 = 3.79$ , $\textcolor{orange}{V} = 7.58$ (2.21: the 5th sample from the right)

↓ into $\textcolor{blue}{K}/(N\textcolor{orange}{V})$ : far out, the estimate falls only like $1/|x|$

\textcolor{purple}{\hat p(x)} \approx \frac{\textcolor{blue}{K}}{2N|x|}

At $x = 6$ : $\textcolor{purple}{\hat p} = 0.016$ , true $p = 0.00005$ : about 308× too high. A Gaussian tail falls like $e^{-x^2/2}$ .

↓ total area of the right tail, from $a$ to $b$

\int_a^b \frac{\textcolor{blue}{K}}{2N x}\,dx = \frac{\textcolor{blue}{K}}{2N} \ln\frac{b}{a} \xrightarrow{\;b \to \infty\;} \infty

TL;DR: k-NN can never say "nothing here": it always stretches until it finds

\textcolor{blue}{K}

samples, so the tails never die out and the total area is infinite. Parzen: each bump has area

\frac{1}{n}

, total exactly 1.

Part IV

Parzen windows

A window around $x$

The fixed circle from before, in 1D: a window of width $h$ centred at $x$ . Five samples, $h = 1.4$ .

At $x = 0.6$ : how many samples fall into the window?

\textcolor{blue}{K = 2}

↓ $\textcolor{blue}{K}/(N\textcolor{orange}{V})$ with $\textcolor{orange}{V = h}$

\textcolor{purple}{\hat p(x)} = \frac{\textcolor{blue}{K}}{n\,\textcolor{orange}{h}}

\textcolor{purple}{\hat p(0.6)} = \dfrac{\textcolor{blue}{2}}{5 \cdot \textcolor{orange}{1.4}} = \textcolor{purple}{0.29}

Seen from one sample $x_i$ : it adds $\frac{1}{nh}$ to every $x$ within $h/2$ of it: a box of width $h$ , height $\frac{1}{nh}$ , so area $\frac{1}{n}$ .

Every sample adds its box: $n \cdot \frac{1}{n}$ = total area 1, and $\hat p \ge 0$ : a proper probability density.

↓ the estimate at $x$ = the boxes that cover $x$ ; $[\cdot]$ is 1 if true, else 0

\hat p(x) = \frac{1}{n} \sum_{i=1}^{n} \frac{1}{h}\, [\, |x - x_i| < h/2 \,]

Smooth bumps instead of boxes

With boxes a sample counts fully or not at all, so $\hat p$ jumps. Let its weight fade with the distance instead.

	boxes (before)	bumps (now)
one sample $x_i$ adds	$\dfrac{1}{n\textcolor{orange}{h}}\,[\,\|x - x_i\| < \tfrac{\textcolor{orange}{h}}{2}\,]$	$\dfrac{1}{n\textcolor{orange}{h}}\,\varphi\!\left(\dfrac{x - x_i}{\textcolor{orange}{h}}\right)$
$\textcolor{orange}{h}$ is the	width of the box	std. deviation
area of one	$\textcolor{purple}{1/n}$	$\textcolor{purple}{1/n}$
all $n$ samples	staircase, area $\textcolor{purple}{1}$	smooth, area $\textcolor{purple}{1}$

Every sample gets its bump; $\varphi(u) = \frac{1}{\sqrt{2\pi}} e^{-u^2/2}$ is the standard normal density, and $\frac{1}{\textcolor{orange}{h}}\varphi(\frac{u}{\textcolor{orange}{h}}) = K_h(u)$ , a normal with std. $\textcolor{orange}{h}$ .

\hat p_h(x) = \frac{1}{n} \sum_{i=1}^{n} K_h(x - x_i)

Parzen windows: the general recipe

Given samples $x_1, \dots, x_n$ , a kernel $K$ and a bandwidth $\textcolor{orange}{h} > 0$ :

\textcolor{purple}{\hat p_h(x)} = \frac{1}{n} \sum_{i=1}^{n} \frac{1}{\textcolor{orange}{h}}\, K\!\left(\frac{x - x_i}{\textcolor{orange}{h}}\right)

$K$ sets the shape of one bump, $\textcolor{orange}{h}$ its width; the $\frac{1}{\textcolor{orange}{h}}$ keeps its area unchanged when we stretch it. In $D$ dimensions: $\frac{1}{\textcolor{orange}{h}^D}$ .

Box: $K(u) = [\,|u| < \frac12\,]$ . Gaussian: $K(u) = \varphi(u)$ , as in the assignment. More shapes in the Optional section.

What must

K

look like?

$K(u) \ge 0$ : no negative density
$\int K(u)\,du = 1$ : total area 1
symmetric, $K(u) = K(-u)$

The first two are needed for a density; the third is a sensible choice.

More: Kernel density estimation

Kernel symmetry

↓ the Gaussian is symmetric: $K_h(u) = K_h(-u)$

K_h(x - x_i) = K_h(x_i - x)

So $\hat p_h(x)$ is also one bump on $x$ , averaged over the samples: a soft count of the samples near $x$ , like the fixed-radius circle.

Try it: the bandwidth $h$

bandwidth

h

show the 40 bumps show the true

p(x)

Part V

Choosing $h$

How big should $h$ be?

The same trade-off as in every model: too flexible fits the noise in the samples, too rigid misses the structure.

Can the training data choose $h$ ?

Idea: pick the $h$ under which our samples are most probable (their log-likelihood; why the log: next slide)

L_{\text{trn}}(h) = \sum_{j=1}^{n} \log \hat p_h(x_j)

As $h \to 0$ , $L_{\text{trn}}(h)$ …

Every training sample sits under its own bump, whose height $\frac{1}{nh}\varphi(0)$ grows without limit.

↓ $x_j$ is itself a sample: its own bump sits on it, with height $\frac{1}{n}\frac{1}{h}\varphi(0)$

\hat p_h(x_j) \ \ge\ \frac{1}{n\,h}\,\varphi(0) \ \xrightarrow{\ h \to 0\ }\ \infty

The "best" $h$ would be 0: one spike per sample, a perfect score on these samples, useless for new ones.

Score

\hat p_h

on samples it has not seen.

More: likelihood function

How do we score $\hat p$ on new samples?

A good estimate puts high density where new samples actually land. Two estimates, three held-out samples:

	$x'_1$	$x'_2$	$x'_3$	product	sum of logs
$\hat p_A$	0.30	0.25	0.20	0.015	−4.2
$\hat p_B$	0.50	0.01	0.40	0.002	−6.2

$B$ is higher at two samples, but nearly misses $x'_2$ : the product punishes that. $A$ wins.

↓ samples are independent, so the probability of seeing all of them is a product: the likelihood

\prod_{x \in X_{\text{tst}}} \hat p(x)

↓ take the log: products of many small numbers underflow (100 samples at 0.05: $10^{-130}$ ); the log turns the product into a sum and keeps the same order, so the same $h$ wins

\log \prod \hat p(x) = \sum_{x \in X_{\text{tst}}} \log \hat p(x)

Higher log-likelihood on held-out samples = better estimate. A sample where

\hat p \approx 0

gives

\log \hat p \to -\infty

: a strong penalty for spikes and gaps.

More: likelihood and log-likelihood

Choosing $h$ with 10-fold cross-validation

In the assignment

crossval(num_data, 10): index sets itrn[i] ( $\textcolor{blue}{X^i_{\text{trn}}}$ ), itst[i] ( $\textcolor{orange}{X^i_{\text{tst}}}$ )
compute_Lh: $L(h)$ for one $h$ ; fminbound: $h^*$ (it minimises, so pass $-L$ )

Split the $n$ samples into 10 equal parts. Fold $i$ holds out part $i$ :

$\textcolor{orange}{X^i_{\text{tst}}}$ : the samples in part $i$ ; $\;\textcolor{blue}{X^i_{\text{trn}}}$ : the other 9 parts

↓ Parzen built from $\textcolor{blue}{X^i_{\text{trn}}}$ only (notebook: $\hat p_h(x \mid X^i_{\text{trn}})$ , " $\mid$ " = "built from")

\hat p^{(i)}_h(x) = \frac{1}{|\textcolor{blue}{X^i_{\text{trn}}}|} \sum_{x_j \in \textcolor{blue}{X^i_{\text{trn}}}} K_h(x - x_j)

↓ score it on the held-out samples it never saw

L^i(h) = \sum_{x \in \textcolor{orange}{X^i_{\text{tst}}}} \log \hat p^{(i)}_h(x)

↓ average over the 10 folds: every sample is held out exactly once

L(h) = \frac{1}{10} \sum_{i=1}^{10} L^i(h)

↓ choose the $h$ that predicts unseen samples best

h^* = \operatorname*{argmax}_h L(h)

More: Cross-validation · scikit-learn: cross-validation · scipy fminbound

Cross-validation, one fold at a time

fold

i

bandwidth

h

Try it: training vs. cross-validated log-likelihood

bandwidth

h

Part VI

Assignment

From densities to a classifier

Images of letters A and C, one measurement $x$ per image. We don't know $p(x \mid A)$ or $p(x \mid C)$ .

↓ one Parzen estimate per class, from that class's training samples, each with its own $h$ from cross-validation

\hat p(x \mid A) = \hat p_{h_A}(x \mid x_A), \quad \hat p(x \mid C) = \hat p_{h_C}(x \mid x_C)

↓ priors: the fraction of training samples in each class

p_K(A) = \frac{N_A}{N_A + N_C}, \quad p_K(C) = \frac{N_C}{N_A + N_C}

↓ 0/1 loss: the most probable class (lab 2), with $\hat p$ in place of $p$

q(x) = \operatorname*{argmax}_{k \in \{A, C\}} p_K(k)\, \hat p(x \mid k)

Then: classify the test images and compute the test error. Different $h_A$ and $h_C$ are fine: each class has its own spread.

The assignment in three parts

1 · Parzen estimate

parzen(x, trn_x, h); $x$ scalar or vector
unit test: the hand exercise
plots for $h \in \{100, 500, 1000, 2000\}$ vs. a normalised histogram

2 · Best

h

compute_Lh(itrn, itst, x, h): 10-fold $L(h)$
h_range = linspace(100, 1000, 19)
np.random.seed(42) before crossval
fminbound for $h_A$ , $h_C$

3 · Classifier

classify_bayes_parzen(x_test, xA, xC, pA, pC, h_bestA, h_bestC)
test error, show_classification

Upload to BRUTE (zip, no folders): parzen.py, parzen.ipynb, parzen_estimates.png, optimal_h_classA.png, optimal_h_classC.png, parzen_classif.png. Bonus: 2D measurements (Optional section).

Optional

Extra material

Other kernels

Any symmetric bump with area 1 works: $\hat p_h(x) = \frac{1}{n} \sum_i K_h(x - x_i)$ with $K_h(u) = \frac{1}{h} K\!\left(\frac{u}{h}\right)$ .

The box gives the window count from before (steps); the Gaussian is smooth everywhere.
"Epanechnikov is optimal": with the best $h$ for each kernel, it gives the smallest error (the AMISE from the rules-of-thumb slide) of all non-negative kernels. The others are close: for the same error the Gaussian needs only about 5% more samples, the box about 7%.
So in practice $h$ matters far more than the shape of the kernel.

More: Kernel (statistics): common kernels and their efficiency

Where do rules of thumb for $h$ come from?

The error of $\hat p_h$ has two parts that pull $h$ in opposite directions.

↓ bias: a wide kernel flattens peaks and fills valleys, by about $\frac{1}{2} h^2 p''(x)$

\int \text{bias}^2 \approx \frac{h^4}{4} \int p''(x)^2\,dx

↓ variance: a narrow window holds few samples, about $n h\, p(x)$ of them

\int \text{var} \approx \frac{1}{n h} \int K(u)^2\,du

↓ add them, set the derivative in $h$ to 0

h^{*5} = \frac{\int K^2}{n \int p''^2}

More: Kernel density estimation: bandwidth selection, AMISE

Silverman's rule: pretend $p$ is Gaussian

$h^*$ needs $\int p''^2$ , which depends on the unknown $p$ .

↓ plug in a Gaussian with the sample's standard deviation $\sigma$ , and a Gaussian kernel

\int p''^2 = \frac{3}{8\sqrt{\pi}\,\sigma^5}, \qquad \int K^2 = \frac{1}{2\sqrt{\pi}}

↓ into $h^{*5} = \int K^2 / (n \int p''^2)$

h = \left(\tfrac{4}{3}\right)^{1/5} \sigma\, n^{-1/5} \approx 1.06\, \sigma\, n^{-1/5}

$n^{-1/5}$ : 10× more data → $h$ only 0.63× smaller.

running example, $n = 40$	$h$
Silverman, $1.06\,\sigma n^{-1/5}$ ( $\sigma = 1.62$ )	0.82
Scott (scipy default), $\sigma n^{-1/5}$	0.77
robust: $0.9 \min(\sigma, \text{IQR}/1.34)\, n^{-1/5}$	0.70
cross-validation	0.57
best, using the true $p$	0.43

The rules oversmooth here: two bumps curve more than one Gaussian with the same $\sigma$ . Fine as a quick default, or as a starting point for CV.

More: Silverman's rule · scipy gaussian_kde (Scott, Silverman) · interquartile range

Adaptive bandwidth: one width per sample

One $h$ is a compromise: too wide where samples are crowded, too narrow where they are sparse. Idea: every sample $x_i$ gets its own width $\textcolor{orange}{h_i}$ .

$h$ : one global width, as before (Silverman or cross-validation).

Pilot $\textcolor{blue}{\tilde p}$ : ordinary Parzen with that $h$ . $\textcolor{blue}{\tilde p(x_i)}$ says how crowded it is around $x_i$ .

$\textcolor{purple}{g}$ : a typical pilot value (geometric mean of all $\tilde p(x_i)$ ). Ratio $\textcolor{blue}{\tilde p(x_i)}/\textcolor{purple}{g}$ : $>1$ crowded, $<1$ sparse.

\textcolor{orange}{h_i} = h \left(\frac{\textcolor{blue}{\tilde p(x_i)}}{\textcolor{purple}{g}}\right)^{-1/2}

4× more crowded than typical → half the width; 4× sparser → twice the width (Abramson's square-root law)

$\hat p(x) = \frac{1}{n}\sum_i K_{\textcolor{orange}{h_i}}(x - x_i)$ : the same recipe; every bump still has area $\frac1n$ , so $\hat p$ is still a density.

Adaptive bandwidth on a skewed sample

60 samples from a skewed density: a sharp peak near 1 and a long thin tail to the right.

Fixed $h$ (from CV): a compromise. The peak comes out too low, and lonely tail samples still make small bumps (around 4.3).

Adaptive $h_i$ : narrow bumps at the peak, so it is closer to the truth; wide bumps in the tail, so it is smooth. Grey bars: the widths $h_i$ of the far samples.

Used in: t-SNE (each point's Gaussian width $\sigma_i$ is set by the perplexity, i.e. by its neighbours) · particle physics, Cranmer 2001, in ROOT as RooKeysPdf · Abramson 1982

The curse of dimensionality, numerically

Standard Gaussian in $D$ dimensions. Window: a cube of side 1 around the peak. How many samples land in it?

$D$	probability inside	$n$ for $K \approx 20$ inside
1	0.38	52
2	0.15	136
3	0.056	356
5	0.0082	2 400
10	0.000068	295 000
20	4.6·10⁻⁹	4.4·10⁹

↓ the coordinates are independent; each one is inside $[-\frac12, \frac12]$ with probability 0.383

P(\text{inside}) = 0.383^D, \qquad n \approx \frac{20}{0.383^D}

A 28×28 image has $D = 784$ . Parzen (and histograms) work in 1–3 dimensions; beyond that we need models that assume structure.

More: Curse of dimensionality

Try it: the same window in 1, 2 and 3 dimensions

samples

n

window side

Bonus: two measurements

$x$ = (left half) − (right half) intensity
$y$ = (top half) − (bottom half) intensity

↓ a round 2D Gaussian bump of width $h$ on every sample, covariance $h^2 I$

K_h(\mathbf{u}) = \frac{1}{2\pi h^2}\, e^{-\|\mathbf{u}\|^2 / (2h^2)}

↓ the same average as in 1D

\hat p_h(\mathbf{x}) = \frac{1}{n} \sum_{i=1}^{n} K_h(\mathbf{x} - \mathbf{x}_i)

Watch the scales: one $h$ for both axes only makes sense if $x$ and $y$ have similar ranges. Rescale them first (e.g. divide each by its standard deviation), or use a separate $h$ per axis.

In $D$ dimensions $V \sim h^D$ : to keep $K$ samples in the window we need far more data. Parzen works well only in low dimensions.

Bonus task: show the estimated

\hat p(\mathbf{x} \mid k)

for both letters, e.g. as contour plots.

More: scipy multivariate_normal

Optional

Where Parzen windows show up

You have seen kernel density estimates before

A violin plot is a Parzen estimate drawn on both sides of an axis: the running example, $h = 0.57$ .

seaborn: sns.kdeplot(x), sns.violinplot(…); bw_adjust scales the default (Scott's) $h$
pandas: df.plot.kde()
maps of hotspots (accidents, crime, GPS traces): a 2D Parzen estimate on a map

Next time a plot looks suspiciously smooth or bumpy: ask which

h

was used.

More: seaborn kdeplot · violin plot

Background subtraction in video

A static camera. Which pixels belong to a moving object?

Per pixel: keep its last $N$ values, fit a 1D Parzen estimate. A new value with low $\hat p$ is foreground.

Why non-parametric: a pixel behind a swaying tree sees leaf or sky: two bumps. One Gaussian would call both "unusual".

Elgammal, Harwood, Davis: Non-parametric model for background subtraction, ECCV 2000

Mean shift: climb $\hat p$ to a mode

x \leftarrow \frac{\sum_i K_h(x - x_i)\, x_i}{\sum_i K_h(x - x_i)}

Move to the kernel-weighted mean of the samples. Every step goes uphill on $\hat p_h$ .

start

x

bandwidth

h

Comaniciu & Meer 2002: clustering, segmentation, tracking · Wikipedia

Parzen windows inside transformers

Parzen regression: to predict $y$ at a new $x$ , average the training outputs $y_i$ , weighting each by how close its $x_i$ is to $x$ .

\hat y(x) = \sum_i w_i\, y_i, \qquad w_i = \frac{K_h(x - x_i)}{\sum_j K_h(x - x_j)}

The weights are $\ge 0$ and sum to 1: close samples count more.

Attention in a transformer does the same: each token's output is a weighted average of the other tokens' values, with larger weights for tokens that are more "similar" to it.

The difference: "similar" is not a fixed kernel with a bandwidth $h$ , but learned during training.

Explained step by step, from this exact formula to attention: Dive into Deep Learning, 11.2: Attention pooling by similarity · the kernel view of attention: Tsai et al. 2019, Transformer dissection

Cautionary tale: Parzen for generative models

2014: the original GAN paper could not compute the likelihood of its samples. So it fitted a Parzen estimate to generated images, chose $h$ by validation, and reported the log-likelihood of the test images.

2016: Theis, van den Oord, Bethge showed that in image dimensions this number is far from the true log-likelihood even with very many samples, and can rank models in the wrong order.

Why: images have hundreds of dimensions. The curse of dimensionality from the Optional section: no sample lands near a test image, so

\hat p

measures distance to the nearest sample, not probability.

The metric was dropped. Lesson: a density estimate is only as good as the dimension it lives in.

Parzen windows

Recap

What have we assumed so far?

How can we estimate a density?

Histogram

We have samples. What is p(x)p(x)p(x)?

How tall should a histogram bar be?

Histogram: pros and cons?

Try it: bin width and bin origin

More dimensions: KKK, NNN and VVV

More dimensions?

Fix KKK, find VVV: k-nearest neighbours

Fix VVV, count KKK: a window

Quick check: fixed KKK or fixed VVV?

k-NN density in 1D: what does KKK do?

What happens far from the data?

Parzen windows

A window around xxx

Smooth bumps instead of boxes

Parzen windows: the general recipe

Kernel symmetry

Try it: the bandwidth hhh

Choosing hhh

How big should hhh be?

Can the training data choose hhh?

How do we score p^\hat pp^​ on new samples?

Choosing hhh with 10-fold cross-validation

Cross-validation, one fold at a time

Try it: training vs. cross-validated log-likelihood

Assignment

From densities to a classifier

The assignment in three parts

Extra material

Other kernels

Where do rules of thumb for hhh come from?

Silverman's rule: pretend ppp is Gaussian

Adaptive bandwidth: one width per sample

Adaptive bandwidth on a skewed sample

The curse of dimensionality, numerically

Try it: the same window in 1, 2 and 3 dimensions

Bonus: two measurements

Where Parzen windows show up

You have seen kernel density estimates before

Background subtraction in video

Mean shift: climb p^\hat pp^​ to a mode

Parzen windows inside transformers

Cautionary tale: Parzen for generative models

Further reading

We have samples. What is $p(x)$ ?

More dimensions: $K$ , $N$ and $V$

Fix $K$ , find $V$ : k-nearest neighbours

Fix $V$ , count $K$ : a window

Quick check: fixed $K$ or fixed $V$ ?

k-NN density in 1D: what does $K$ do?

A window around $x$

Try it: the bandwidth $h$

Choosing $h$

How big should $h$ be?

Can the training data choose $h$ ?

How do we score $\hat p$ on new samples?

Choosing $h$ with 10-fold cross-validation

Where do rules of thumb for $h$ come from?

Silverman's rule: pretend $p$ is Gaussian

Mean shift: climb $\hat p$ to a mode