Density Estimation Fundamentals

Kernel Density Estimation (KDE) is the premier nonparametric method for estimating probability density functions. First proposed by Murray Rosenblatt in 1956 and independently by Emanuel Parzen in 1962, KDE elegantly addresses the limitations of histograms while maintaining the flexibility to capture arbitrary distributional shapes.

The core insight is profound in its simplicity: rather than counting points that fall into predefined bins, we place a smooth "bump" at each data point and sum these bumps to create a continuous density estimate. This approach eliminates the histogram's discontinuities, removes dependence on arbitrary bin origins, and achieves faster convergence rates.

KDE has become ubiquitous in modern data science, underpinning applications from visualization tools (density plots, violin plots) to sophisticated machine learning methods (mean shift clustering, anomaly detection, generative modeling). Understanding KDE deeply is essential for any practitioner working with probability distributions.

Given observations $\mathcal{D} = {x_1, x_2, \ldots, x_n}$ drawn i.i.d. from unknown density $f(x)$, the kernel density estimator is defined as:

$$\hat{f}h(x) = \frac{1}{nh} \sum{i=1}^{n} K\left(\frac{x - x_i}{h}\right)$$

where:

• $K(\cdot)$ is the kernel function — a symmetric, non-negative function integrating to 1
• $h > 0$ is the bandwidth (also called smoothing parameter or window width)
• The argument $(x - x_i)/h$ measures how far $x$ is from $x_i$ in units of bandwidth

Intuition: At each data point $x_i$, we place a kernel with height $1/(nh)$ and spread $h$. As we move along the $x$-axis, we sum contributions from all kernels. Near clusters of data, many kernels contribute, yielding high density. In sparse regions, few kernels contribute, yielding low density.

Equivalent Form:

Define the scaled kernel $K_h(u) = \frac{1}{h}K\left(\frac{u}{h}\right)$. Then:

$$\hat{f}h(x) = \frac{1}{n} \sum{i=1}^{n} K_h(x - x_i)$$

This shows KDE as an average of kernels centered at each data point.

Why This Works:

Consider the expected value of $\hat{f}_h(x)$ at a point $x$: $$\mathbb{E}[\hat{f}_h(x)] = \mathbb{E}\left[K_h(x - X)\right] = \int K_h(x - t) f(t) , dt = (K_h * f)(x)$$

This is the convolution of the kernel with the true density. For small $h$, $K_h$ approximates a Dirac delta, so $(K_h * f)(x) \approx f(x)$. The kernel acts as a smoothing filter on the density.

Verification of Valid Density:

$$\int \hat{f}h(x) , dx = \frac{1}{n} \sum{i=1}^{n} \int K_h(x - x_i) , dx = \frac{1}{n} \sum_{i=1}^{n} 1 = 1$$

Since each kernel integrates to 1, so does their average. Also, $\hat{f}_h(x) \geq 0$ everywhere (assuming $K \geq 0$), so KDE produces a valid probability density.

A kernel function $K(u)$ must satisfy the following conditions:

1. Non-negativity: $K(u) \geq 0$ for all $u$ (ensures $\hat{f}_h \geq 0$)
2. Normalization: $\int_{-\infty}^{\infty} K(u) , du = 1$ (ensures $\int \hat{f}_h = 1$)
3. Symmetry: $K(-u) = K(u)$ for all $u$ (ensures unbiased location)

Additionally, most kernels satisfy: 4. Zero mean: $\int u \cdot K(u) , du = 0$ (from symmetry) 5. Finite variance: $\int u^2 \cdot K(u) , du = \mu_2(K) < \infty$

The quantity $\mu_2(K) = \int u^2 K(u) , du$ is the second moment of the kernel, which appears in bias calculations.

The Gaussian Kernel:

The Gaussian kernel is by far the most widely used:

$$K(u) = \frac{1}{\sqrt{2\pi}} e^{-u^2/2}$$

Why Gaussian is Popular:

1. Infinite differentiability: $\hat{f}_h$ inherits all derivatives of $K$, making KDE with Gaussian kernel infinitely smooth.

2. Infinite support: Every data point contributes (albeit negligibly) to the estimate everywhere. No arbitrary cutoffs.

3. Mathematical convenience: Convolution of Gaussians is Gaussian. Many theoretical results simplify.

4. Connection to diffusion: KDE with Gaussian kernel can be viewed as heat diffusion from point sources.

5. Ubiquitous availability: Every statistical software implements Gaussian KDE.

Disadvantage: Unbounded support means computing $\hat{f}_h(x)$ technically requires summing $n$ terms, though in practice far-away points contribute negligibly.

The Epanechnikov Kernel:

$$K(u) = \frac{3}{4}(1 - u^2) \cdot \mathbf{1}(|u| \leq 1)$$

This parabola (inverted) is theoretically optimal in the sense of minimizing MISE among all kernels when bandwidth is optimally chosen. The proof, due to Epanechnikov (1969), shows that the optimal kernel shape for minimizing $\int \text{Var}[\hat{f}_h(x)] dx$ subject to a fixed $\int \text{Bias}^2[\hat{f}_h(x)] dx$ is a parabola.

In Practice:

Despite its optimality, the efficiency gain over Gaussian is small—approximately 5% in MISE. Since kernel choice matters much less than bandwidth selection, the convenience of Gaussian often outweighs the marginal efficiency of Epanechnikov.

Bounded Support Advantage: For compact-support kernels (Epanechnikov, Uniform, etc.), computing $\hat{f}_h(x)$ only requires summing over points within distance $h$ of $x$, enabling $O(n)$ or better algorithms with spatial data structures.

The statistical properties of KDE reveal precisely how bandwidth controls the bias-variance tradeoff.

Expected Value (Bias Derivation):

Recall $\mathbb{E}[\hat{f}_h(x)] = \int K_h(x-t) f(t) , dt$. Using Taylor expansion of $f(t)$ around $x$:

$$f(t) = f(x) + (t-x)f'(x) + \frac{(t-x)^2}{2}f''(x) + O((t-x)^3)$$

Substituting and using symmetry of $K$ (which zeros the $f'$ term):

$$\mathbb{E}[\hat{f}_h(x)] = f(x) \int K(u) , du + \frac{h^2 f''(x)}{2} \int u^2 K(u) , du + O(h^4)$$

$$= f(x) + \frac{h^2 f''(x)}{2} \mu_2(K) + O(h^4)$$

where $\mu_2(K) = \int u^2 K(u) , du$ is the kernel's second moment.

Bias: $$\text{Bias}[\hat{f}_h(x)] = \mathbb{E}[\hat{f}_h(x)] - f(x) \approx \frac{h^2}{2} f''(x) \mu_2(K)$$

Bias is proportional to:

• $h^2$: smaller bandwidth → smaller bias
• $f''(x)$: high curvature regions → more bias
• $\mu_2(K)$: wider kernels → more bias

Variance Derivation:

Since $\hat{f}h(x) = \frac{1}{n}\sum{i=1}^n K_h(x-x_i)$ is an average of i.i.d. random variables:

$$\text{Var}[\hat{f}_h(x)] = \frac{1}{n} \text{Var}[K_h(x-X)] = \frac{1}{n}\left(\mathbb{E}[K_h^2(x-X)] - \mathbb{E}[K_h(x-X)]^2\right)$$

For small $h$, the second term is approximately $f(x)^2$, and:

$$\mathbb{E}[K_h^2(x-X)] = \int K_h^2(x-t) f(t) , dt = \frac{1}{h^2} \int K^2\left(\frac{x-t}{h}\right) f(t) , dt$$

$$= \frac{f(x)}{h} \int K^2(u) , du + O(1) = \frac{f(x) R(K)}{h}$$

where $R(K) = \int K^2(u) , du$ is the roughness of the kernel.

Variance: $$\text{Var}[\hat{f}_h(x)] \approx \frac{f(x) R(K)}{nh}$$

Variance is proportional to:

• $1/(nh)$: more data or wider bandwidth → less variance
• $f(x)$: higher density → more variance
• $R(K)$: rougher kernels → more variance

Mean Squared Error (MSE):

$$\text{MSE}[\hat{f}_h(x)] = \text{Bias}^2 + \text{Var} \approx \frac{h^4}{4} [f''(x)]^2 [\mu_2(K)]^2 + \frac{f(x) R(K)}{nh}$$

Mean Integrated Squared Error (MISE):

$$\text{MISE}(\hat{f}_h) = \int \text{MSE}[\hat{f}_h(x)] , dx$$

$$\approx \frac{h^4}{4} [\mu_2(K)]^2 R(f'') + \frac{R(K)}{nh}$$

where $R(f'') = \int [f''(x)]^2 , dx$ is the roughness of the true density's second derivative.

Deriving the Optimal Bandwidth:

To find the MISE-optimal $h$, differentiate MISE with respect to $h$ and set to zero:

$$\frac{d}{dh} \text{MISE} = h^3 [\mu_2(K)]^2 R(f'') - \frac{R(K)}{nh^2} = 0$$

Solving:

$$h^* = \left(\frac{R(K)}{n [\mu_2(K)]^2 R(f'')}\right)^{1/5}$$

Order Notation: $$h^* = O(n^{-1/5})$$

The optimal bandwidth shrinks as the fifth root of sample size.

Optimal MISE:

$$\text{MISE}^* = \frac{5}{4} \left(\frac{[\mu_2(K)]^2 R(K)^4}{n^4} R(f'')\right)^{1/5} = O(n^{-4/5})$$

Thus KDE converges at rate $n^{-4/5}$, faster than histograms ($n^{-2/3}$) but slower than parametric ($n^{-1}$).

Interpretation of Optimal Bandwidth:

The formula for $h^*$ reveals key insights:

1. Sample Size Effect: $h^* \propto n^{-1/5}$ — more data allows finer resolution.

   • 32× more data → bandwidth halved
   • Slow reduction: even 100× data only reduces $h$ by ~60%

2. Density Roughness Effect: $h^* \propto R(f'')^{-1/5}$ — rough densities need smaller bandwidths.

   • High curvature → more rapid changes → need fine resolution
   • Flat density → large bandwidth acceptable

3. Kernel Effect: Through $R(K)$ and $\mu_2(K)$.

   • Epanechnikov minimizes the relevant constant
   • But differences between kernels are minor

The Chicken-and-Egg Problem:

$R(f'') = \int [f''(x)]^2 dx$ involves the unknown density $f$. We cannot compute the optimal bandwidth exactly without already knowing $f$. This motivates:

1. Plug-in rules: estimate $R(f'')$ from a preliminary density estimate
2. Reference rules: assume $f$ is Gaussian and use the known formula
3. Cross-validation: optimize a criterion that doesn't require $f$

Silverman's Rule of Thumb (1986):

Assuming $f$ is Gaussian $N(\mu, \sigma^2)$, we have $R(f'') = \frac{3}{8\sqrt{\pi}\sigma^5}$.

For the Gaussian kernel with $R(K) = (2\sqrt{\pi})^{-1}$ and $\mu_2(K) = 1$:

$$h^* = \left(\frac{4\hat{\sigma}^5}{3n}\right)^{1/5} \approx 1.06 \cdot \hat{\sigma} \cdot n^{-1/5}$$

To handle non-normal data, Silverman suggested using a robust scale estimate:

$$h = 0.9 \cdot \min\left(\hat{\sigma}, \frac{\text{IQR}}{1.34}\right) \cdot n^{-1/5}$$

The factor 0.9 adjusts for the different roughness assumptions. IQR/1.34 equals $\sigma$ for Gaussian but is more robust.

Properties of Silverman's Rule:

• Simple to compute
• Works well for unimodal, roughly symmetric distributions
• Oversmooths multimodal distributions (assumes single mode)

When the density has bounded support (e.g., positive data, proportions in [0,1]), standard KDE suffers from boundary bias—systematic errors near the edges of the support.

The Problem:

Consider data on $[0, \infty)$. At $x = 0$, the kernel extends into negative values where no data can exist. The integral:

$$\int_{-\infty}^{\infty} K_h(0 - t) f(t) , dt$$

is evaluated over $[0, \infty)$, but the kernel expects contribution from $(-\infty, 0)$ as well. Since $f(t) = 0$ for $t < 0$, the estimate is biased downward near 0.

Quantifying Boundary Bias:

At a boundary point $x = 0$, the relative bias is approximately $O(1)$ rather than $O(h^2)$ in the interior. This is a dramatic increase—the estimate can be biased by 50% or more!

Correction Methods:

1. Reflection Method:

Reflect data across the boundary and include both original and reflected points: $$\hat{f}h(x) = \frac{1}{nh}\sum{i=1}^n \left[K\left(\frac{x-x_i}{h}\right) + K\left(\frac{x+x_i}{h}\right)\right]$$

This implicitly assumes $f'(0) = 0$ (slope zero at boundary).

2. Boundary Kernels:

Use modified kernels near boundaries that integrate to 1 over the restricted domain: $$K_b(u; x, h) = \frac{K(u)}{\int_{L(x)}^{R(x)} K(v) dv}$$

where $L(x)$ and $R(x)$ are the effective kernel limits given the boundary.

3. Pseudodata Methods:

Generate pseudodata beyond the boundary using local linear extrapolation.

4. Local Polynomial Methods:

Fit a local polynomial (rather than constant) at boundary regions, which automatically corrects for boundary effects.

5. Transformation:

Map data to an unbounded domain (e.g., log transform for positive data), apply KDE, and back-transform.

Beta Kernel for [0,1] Data:

For data on [0, 1], a specialized solution uses Beta distribution kernels:

$$K_{\text{Beta}}(t; x, h) = \frac{t^{x/h}(1-t)^{(1-x)/h}}{B(x/h + 1, (1-x)/h + 1)}$$

The kernel is centered at $x$ but naturally confined to [0, 1]. As $h \to 0$, it concentrates around $x$. This eliminates boundary bias entirely by construction.

Gamma Kernel for Positive Data:

Similarly, Gamma kernels can be used for data on $[0, \infty)$: $$K_{\text{Gamma}}(t; x, h) = \text{Gamma distribution centered near } x$$

These specialized kernels achieve optimal rates at boundaries but require different bandwidth selection procedures.

Extending KDE to multiple dimensions is conceptually straightforward but introduces significant computational and statistical challenges.

Product Kernel:

The simplest multivariate extension uses a product of univariate kernels: $$K_H(\mathbf{u}) = \prod_{j=1}^d K(u_j / h_j) / \prod_{j=1}^d h_j$$

with different bandwidths $h_j$ for each dimension. The estimate is: $$\hat{f}H(\mathbf{x}) = \frac{1}{n} \sum{i=1}^n K_H(\mathbf{x} - \mathbf{x}_i)$$

Full Bandwidth Matrix:

More generally, use a bandwidth matrix $H$ (positive definite $d \times d$ matrix): $$K_H(\mathbf{u}) = |H|^{-1/2} K(H^{-1/2}\mathbf{u})$$

where $K(\cdot)$ is a multivariate kernel (often multivariate Gaussian).

For multivariate Gaussian kernel: $$K_H(\mathbf{u}) = (2\pi)^{-d/2} |H|^{-1/2} \exp\left(-\frac{1}{2}\mathbf{u}^T H^{-1} \mathbf{u}\right)$$

The bandwidth matrix $H$ controls smoothing in all directions and correlations.

The Curse of Dimensionality in KDE:

Recall the optimal MISE convergence rate: $$\text{MISE}^* = O(n^{-4/(d+4)})$$

Beyond $d \approx 5$, KDE becomes impractical without enormous datasets. This is the "curse of dimensionality" in action—high-dimensional spaces are effectively empty, and data become too sparse for nonparametric estimation.

Bandwidth Matrix Selection:

In $d$ dimensions, the full bandwidth matrix $H$ has $d(d+1)/2$ free parameters. This is usually too many to estimate reliably. Common simplifications:

1. Scalar bandwidth: $H = h^2 I$ — same bandwidth in all directions. Simplest but doesn't adapt to different scales or correlations.

2. Diagonal bandwidth: $H = \text{diag}(h_1^2, \ldots, h_d^2)$ — different bandwidths per dimension but no correlation. Often standardize data first.

3. Full matrix: Allow arbitrary positive definite $H$. Most flexible but requires more data to estimate. Often set $H \propto \widehat{\Sigma}$, the sample covariance matrix.

Scott's Rule for multivariate: $$H = n^{-2/(d+4)} \widehat{\Sigma}$$

Silverman's Rule for multivariate: $$H = \left(\frac{4}{(d+2)n}\right)^{2/(d+4)} \widehat{\Sigma}$$

Both scale the sample covariance by a factor depending on $n$ and $d$.

Naive KDE evaluation requires $O(n)$ operations per query point—for $m$ query points, total cost is $O(nm)$. With modern datasets having millions of points, this becomes prohibitive.

Computational Complexity:

For $n$ data points and $m$ evaluation points:

• Naive: $O(nm)$ kernel evaluations
• With Gaussian kernel: Each evaluation involves exponentiating and summing
• For visualization on grid: $m$ can be large (e.g., 1000×1000 = 10⁶)

Fast Algorithms:

FFT-Based KDE:

The key insight is that KDE at grid points is a discrete convolution: $$\hat{f}(x_k) = \frac{1}{n} \sum_{i=1}^n K_h(x_k - x_i)$$

If we bin the data onto the same grid:

1. Create histogram counts $c_j$ on grid ${g_1, \ldots, g_m}$
2. KDE at grid points is approximately $\hat{f}(g_k) \approx \sum_j c_j K_h(g_k - g_j) / n$
3. This is a convolution: $\hat{f} \approx (c * K_h) / n$
4. Apply FFT: $\mathcal{F}(\hat{f}) = \mathcal{F}(c) \cdot \mathcal{F}(K_h) / n$
5. Inverse FFT to get $\hat{f}$

Total complexity: $O(m \log m)$ for the FFT operations, plus $O(n)$ for binning.

Error Analysis: Binning introduces bias proportional to $(\Delta g)^2$ where $\Delta g$ is grid spacing. With grid spacing much smaller than bandwidth, error is negligible.

Kernel density estimation represents the gold standard for nonparametric density estimation. We've covered its theoretical foundations and practical aspects in depth:

What's Next:

Bandwidth selection is arguably the most critical practical concern in KDE—more important than kernel choice. The next page focuses entirely on bandwidth selection methods, including Silverman's rule, plug-in estimators, and cross-validation approaches. We'll develop both theoretical understanding and practical guidance for choosing bandwidths that balance smoothness against fidelity to data.