Kernel density estimation

Section of Hierarchical clustering and kernel density estimation, part of Statistical learning

Type	Section
Taught in	Hierarchical clustering and kernel density estimation
Part	Statistical learning

Place probability mass 1/n at each observation. That empirical distribution remembers the sample exactly, and it has no ordinary density with respect to volume. A kernel estimate replaces each point mass by a small density bump. Its width controls how far the influence of an observation extends. With Gaussian bumps this is a mixture with one equally weighted component per observation. The bandwidth is chosen directly, in place of a component count.

Definition 1 (Kernel density estimator) For a sample \{\mathbf{x}_1,\ldots,\mathbf{x}_n\} in p dimensions, take K to be a symmetric probability density and H a symmetric positive definite bandwidth matrix. The kernel density estimate is

\widehat f(\mathbf{x};H)=\frac1n\sum_{i=1}^{n}K_H(\mathbf{x}-\mathbf{x}_i), \qquad K_H(\mathbf{x})=|H|^{-1/2}K(H^{-1/2}\mathbf{x}),

with H^{-1/2} obtained from the spectral decomposition. These assumptions give a nonnegative estimate that integrates to one. More general kernel constructions exist. The formulas below use this density-kernel setting.

Article: Kernel density estimator →

The factor |H|^{-1/2} is the Jacobian of \mathbf{u}=H^{-1/2}\mathbf{x}. It is there so that K_H integrates to one whenever K does. In one term,

\int K_H(\mathbf{x}-\mathbf{x}_i)\,d\mathbf{x} =|H|^{-1/2}\int K\bigl(H^{-1/2}(\mathbf{x}-\mathbf{x}_i)\bigr)\,d\mathbf{x}.

Put \mathbf{u}=H^{-1/2}(\mathbf{x}-\mathbf{x}_i), so d\mathbf{x}=|H|^{1/2}\,d\mathbf{u}, because the Jacobian determinant of \mathbf{x}=H^{1/2}\mathbf{u}+\mathbf{x}_i is |H^{1/2}|=|H|^{1/2}. The integral collapses to \int K(\mathbf{u})\,d\mathbf{u}=1. Each of the n terms integrates to one, and the factor 1/n gives total mass one. That is the change of variables asked for in the exercises.

One kernel is the standard multivariate normal density,

K(\mathbf{x})=\frac{1}{(2\pi)^{p/2}}e^{-\frac12\mathbf{x}^\top\mathbf{x}}, \qquad K_H(\mathbf{x})=\frac{1}{(2\pi)^{p/2}\,|H|^{1/2}}e^{-\frac12\mathbf{x}^\top H^{-1}\mathbf{x}},

so K_H is the \operatorname{N}_p(\mathbf{0},H) density. Each sample point is the centre of one such normal. The mixture reading is then exact: n components, weights 1/n, means \mathbf{x}_i, and one common covariance H. The components, the weights, and the means are fixed by the sample. The choice that remains is H.

The Epanechnikov kernel has bounded support, so a sample point contributes nothing beyond a finite radius. The normal kernel is the usual multivariate default, and its support is the whole space. A full matrix H can follow an elliptical cloud. A diagonal H, or a product of separate univariate smoothers, has axis-aligned contours. On a tilted cloud those contours place part of each kernel’s mass off the data. One H everywhere is a further restriction: a cloud whose shape changes from one region to another is smoothed with the same orientation in every region.

Choosing the bandwidth

Main article: Choosing the bandwidth

The matrix H may be constrained to be diagonal, or left unconstrained. Either way the target is approximation of the density f that generated the sample. There is no unique choice, and the multivariate problem is not a coordinate-wise copy of the univariate one. Likelihood cross-validation, as on Model selection, fits the estimate on one part of the sample and scores the log-likelihood on the left-out part. The bandwidth with the highest held-out score is preferred. With a short sample the split is repeated over folds.

A standard distance between f and \widehat f, different from Kullback–Leibler divergence, is the mean integrated squared error,

\mathrm{MISE}(H) =\mathbb{E}\int_{\mathbb{R}^p}\bigl(\widehat f(\mathbf{x};H)-f(\mathbf{x})\bigr)^2\,d\mathbf{x},

regarded as a function of H. The optimal bandwidth minimises MISE over symmetric positive definite p\times p matrices. The expectation is over the sample that builds \widehat f. MISE can be estimated in several ways, including cross-validation. Most methods optimise H numerically. A full matrix has p(p+1)/2 free entries, which is already a difficult optimisation at moderate p. Large n adds a computational cost. Likelihood cross-validation, or a single validation set, is one practical selector. In high dimension one usually reduces the coordinates first, and that reduction discards information.

Reducing dimension before density estimation can make the smoothing problem manageable, and it changes the target to the projected distribution. Discarded directions cannot be reconstructed from the density of the retained scores alone. There is no universal dimension, such as six, at which estimation becomes impossible: difficulty depends on sample size, smoothness, and effective dimension. The software’s supported dimension is a separate constraint.

Without a smoothness restriction, many irregular densities are compatible with the same finite sample, and the sample does not choose among them. Smoothness is the assumption that lets the estimate take values between the observed points.

Bias and variance in one bandwidth

Main article: Bias and variance in one bandwidth

In one dimension, write H=h^2. A very small bandwidth gives a spiky estimate with one bump per observation: low bias, high variance. A very large bandwidth gives an over-smoothed estimate that hides structure such as two modes: high bias, low variance. MISE balances the two. On iris petal length, three Gaussian estimates with bandwidths h=0.23, h=0.58, and h=1.46 show the three regimes, and a two-dimensional estimate can use both petal measurements at once. The likelihood on the training data favours h\to 0. Each component then sits on its observation and the density at that point diverges, the same singularity as a mixture component collapsing on a data point on the mixture page. The bandwidth is chosen by a held-out score or by an estimate of MISE, rather than by the training likelihood.

For a Gaussian kernel the estimate is a mixture of bumps centred at the data. Small bandwidths retain fine detail and sampling noise. Large bandwidths erase both. Choose the scale using a criterion matched to the goal, then inspect how nearby choices change the picture. A visible mode is a feature of this smoothed estimate. It is evidence of a separate population when something beyond the smoother says so.

Kernel density estimation

Choosing the bandwidth

Bias and variance in one bandwidth

Used in

See also