If kernel density estimation were a musical instrument, bandwidth would be its tuning—get it right, and the estimate sings with clarity and precision; get it wrong, and the result is either screeching noise (undersmoothing) or muffled silence (oversmoothing).

We established in the previous page that bandwidth profoundly impacts KDE quality: too small produces spiky, overfitted estimates; too large obliterates important features. The optimal bandwidth depends on the unknown true density $f$—creating a fundamental circular problem. How do we select a bandwidth without already knowing what we're trying to estimate?

This page is devoted entirely to answering this question. We'll explore multiple philosophies and methods: quick rule-of-thumb calculations, sophisticated plug-in estimators, principled cross-validation approaches, and adaptive techniques that vary bandwidth locally. Each method embodies different assumptions and tradeoffs, and understanding this landscape is essential for applying KDE effectively in practice.

Recall the fundamental result for optimal bandwidth. Minimizing the Mean Integrated Squared Error (MISE) yields:

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

where:

• $R(K) = \int K^2(u) , du$ — roughness of the kernel (known)
• $\mu_2(K) = \int u^2 K(u) , du$ — second moment of the kernel (known)
• $R(f'') = \int [f''(x)]^2 , dx$ — roughness of the true density's second derivative (unknown)

The problem is clear: the optimal bandwidth depends on $R(f'')$, which involves the unknown density $f$. This is not merely inconvenient—it's a fundamental epistemological challenge. Different approaches to bandwidth selection correspond to different strategies for breaking this circularity.

Evaluation Criteria for Bandwidth Selectors:

How do we judge a bandwidth selector? Key criteria include:

1. Accuracy: How close is selected $\hat{h}$ to optimal $h^*$?

   • Measured by relative error: $(\hat{h} - h^)/h^$
   • Or by ratio of achieved MISE to optimal MISE

2. Consistency: Does $\hat{h}/h^* \to 1$ as $n \to \infty$?

3. Convergence Rate: How fast does $\hat{h}$ converge to $h^*$?

   • Relative error $O(n^{-\alpha})$ for various $\alpha$

4. Robustness: Performance across different true densities?

5. Computational Cost: Time complexity as function of $n$?

6. Simplicity: Ease of implementation and understanding?

No method dominates on all criteria—hence the proliferation of approaches.

Reference rules assume the true density belongs to a specific family (typically Gaussian), compute the optimal bandwidth for that reference, and substitute sample estimates. They're simple, fast, and often sufficient for exploratory analysis.

Silverman's Rule (1986):

For a Gaussian density $f = N(\mu, \sigma^2)$: $$R(f'') = \frac{3}{8\sqrt{\pi}\sigma^5}$$

With Gaussian kernel ($R(K) = 1/(2\sqrt{\pi})$, $\mu_2(K) = 1$):

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

Robust Version (Default in most software): $$h = 0.9 \cdot \min\left(\hat{\sigma}, \frac{\text{IQR}}{1.34}\right) \cdot n^{-1/5}$$

The factor 0.9 adjusts for using the minimum of two scale estimates (slightly lower than 1.06 to prevent oversmoothing). The IQR/1.34 equals $\sigma$ for Gaussian but is robust to outliers.

Scott's Rule (1992):

$$h = 1.059 \cdot \hat{\sigma} \cdot n^{-1/5}$$

Essentially the same as Silverman's non-robust version. Sometimes used with the coefficient rounded to 1.06.

Comparison of Coefficients:

For Gaussian kernel and Gaussian reference:

• Exact optimal: $\left(\frac{4}{3n}\right)^{1/5}\sigma \approx 1.059 \cdot \sigma \cdot n^{-1/5}$
• Scott: $1.059 \cdot \hat{\sigma} \cdot n^{-1/5}$
• Silverman (robust): $0.9 \cdot \hat{\sigma} \cdot n^{-1/5}$ (with IQR adjustment)

For other kernels:

The coefficient changes based on $R(K)$ and $\mu_2(K)$. For Epanechnikov kernel: $$h = 2.34 \cdot \hat{\sigma} \cdot n^{-1/5}$$

(The larger coefficient reflects the kernel's shape, not that it requires wider bandwidth.)

Oversmoothing Rule:

Sheather and Jones (1991) proposed using an oversmoothed reference to prevent undersmoothing: $$h_{\text{os}} = \left(\frac{243 R(K)}{35 [\mu_2(K)]^2 n}\right)^{1/5} \hat{\sigma}$$

This is designed to be larger than optimal, ensuring that features aren't spuriously created. One can then reduce from this starting point.

Multivariate Reference Rules:

For $d$-dimensional data with bandwidth matrix $H \propto \widehat{\Sigma}$:

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

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

These become increasingly unreliable as $d$ grows, because the Gaussian reference becomes less representative of complex multivariate structure.

Plug-in methods estimate $R(f'') = \int [f''(x)]^2 dx$ directly from the data, then substitute ("plug in") this estimate into the optimal bandwidth formula. The key insight is that $R(f'')$ can itself be estimated by kernel methods.

The Idea:

To estimate $R(f'') = \int [f''(x)]^2 dx$, we can:

1. Estimate $f''(x)$ using a kernel derivative estimator
2. Integrate the squared estimate

But the kernel estimator of $f''$ has its own bandwidth $g$, which itself depends on $R(f^{(4)})$. This creates a hierarchy of estimation problems.

Sheather-Jones Plug-in (1991):

The most sophisticated and widely-used plug-in method. It solves the bandwidth selection problem by:

1. Estimate $R(f^{(6)})$ using a normal reference
2. Use this to set bandwidth $g_6$ for estimating $f^{(4)}$
3. Estimate $R(f^{(4)})$ using bandwidth $g_6$
4. Use this to set bandwidth $g_4$ for estimating $f''$
5. Estimate $R(f'')$ using bandwidth $g_4$
6. Plug $\widehat{R(f'')}$ into the $h^*$ formula

Mathematical Details:

To estimate $\psi_r = \int f^{(r)}(x)^2 dx$ (where $f^{(r)}$ is the $r$-th derivative), use:

$$\hat{\psi}r = \frac{1}{n^2 g^{r+1}} \sum{i=1}^n \sum_{j=1}^n K^{(r)}\left(\frac{x_i - x_j}{g}\right)$$

where $K^{(r)}$ is the $r$-th derivative of the kernel.

The bandwidth $g$ for estimating $\psi_r$ depends on $\psi_{r+2}$. This creates the hierarchy:

$$g_4 = C_4(K) \cdot \psi_6^{-1/(5+6)} \cdot n^{-1/(5+6)}$$ $$g_2 = C_2(K) \cdot \psi_4^{-1/(5+4)} \cdot n^{-1/(5+4)}$$ $$h = C_0(K) \cdot \psi_2^{-1/(5+2)} \cdot n^{-1/(5+2)}$$

where $C_r(K)$ are constants depending on the kernel.

The chain is stopped by using a Gaussian reference for the highest derivative.

Properties of Plug-in Methods:

Advantages:

• Faster convergence rate than cross-validation: $\hat{h}/h^* - 1 = O_p(n^{-5/14})$ vs $O_p(n^{-1/10})$ for CV
• Computationally efficient: $O(n \log n)$ with FFT
• Adapts to non-Gaussian densities better than reference rules

Disadvantages:

• Still relies on normal reference at the top of the hierarchy
• Can fail for very complex densities (many modes)
• More complex to implement than reference rules

Direct Plug-in (DPI) Variant:

Simpler than Sheather-Jones, stops the hierarchy earlier:

1. Estimate $\psi_4$ using normal reference for $\psi_6$
2. Estimate $\psi_2$ using the estimated $\psi_4$
3. Compute $h$ from estimated $\psi_2$

Slightly less accurate but easier to implement.

Cross-validation (CV) methods select bandwidth by optimizing an estimate of the density estimator's error that doesn't require knowing the true density. Unlike plug-in methods, CV makes no assumptions about the shape of $f$—it is fully nonparametric.

Least-Squares Cross-Validation (LSCV):

Also called Unbiased Cross-Validation (UCV). We want to minimize the Integrated Squared Error: $$\text{ISE}(h) = \int [\hat{f}_h(x) - f(x)]^2 dx = \int \hat{f}_h^2 - 2\int \hat{f}_h f + \int f^2$$

The last term $\int f^2$ doesn't depend on $h$, so we minimize: $$J(h) = \int \hat{f}_h^2 - 2\int \hat{f}_h f$$

The first integral can be computed directly. The second term $\int \hat{f}_h f = \mathbb{E}_f[\hat{f}_h(X)]$ can be estimated by leave-one-out:

$$\widehat{\int \hat{f}h f} = \frac{1}{n} \sum{i=1}^n \hat{f}_{h,-i}(x_i)$$

where $\hat{f}_{h,-i}(x_i)$ is the KDE at $x_i$ using all points except $x_i$.

LSCV Score:

$$\text{LSCV}(h) = \int \hat{f}h^2(x) dx - \frac{2}{n} \sum{i=1}^n \hat{f}_{h,-i}(x_i)$$

For Gaussian kernel, both terms have closed-form expressions:

$$\int \hat{f}h^2 dx = \frac{1}{n^2 h \sqrt{4\pi}} \sum{i=1}^n \sum_{j=1}^n \exp\left(-\frac{(x_i - x_j)^2}{4h^2}\right)$$

$$\hat{f}{h,-i}(x_i) = \frac{1}{(n-1)h\sqrt{2\pi}} \sum{j eq i} \exp\left(-\frac{(x_i - x_j)^2}{2h^2}\right)$$

Minimize LSCV$(h)$ over a grid of $h$ values.

Properties of LSCV:

• Consistent: $\hat{h}_{\text{LSCV}} / h^* \xrightarrow{P} 1$
• Slow convergence: relative error $O_p(n^{-1/10})$
• High variability: the LSCV curve is often wiggly, making optimization difficult
• Tends to undersmooth for small samples

Biased Cross-Validation (BCV):

A variant that reduces the variability of LSCV at the cost of some bias:

$$\text{BCV}(h) = \int \hat{f}h^2 + \frac{R(K)}{nh} - \frac{2}{n} \sum{i=1}^n \hat{f}_{h,-i}(x_i) - \frac{2K(0)}{nh}$$

The additional terms stabilize the objective function, reducing variance but introducing small bias.

Likelihood Cross-Validation (LCV):

Instead of minimizing ISE, maximize the leave-one-out log-likelihood:

$$\text{LCV}(h) = \sum_{i=1}^n \log \hat{f}_{h,-i}(x_i)$$

This measures how well the density predicts held-out points—maximizing average log-density of truly sampled points.

Equivalent form (Pseudo-likelihood): $$\text{LCV}(h) = n \cdot \frac{1}{n}\sum_{i=1}^n \log \hat{f}_{h,-i}(x_i)$$

Related to minimizing KL divergence $D_{KL}(f | \hat{f})$ rather than ISE.

How do the various bandwidth selection methods compare in practice? Extensive simulation studies (e.g., Park & Marron 1990, Jones et al. 1996) have evaluated methods across different density types.

Key Findings from Simulation Studies:

1. For Unimodal Densities: Reference rules and plug-in methods perform similarly well. Simpler is often better.

2. For Multimodal Densities: Cross-validation and plug-in methods outperform reference rules, which systematically oversmooth.

3. For Small Samples (n < 100): All methods have high variability. Reference rules are often as good as anything.

4. For Large Samples (n > 1000): Plug-in methods show their theoretical advantage in convergence rate.

5. Sheather-Jones Dominance: Across many scenarios, Sheather-Jones offers the best balance of speed and accuracy.

6. CV Variability: Cross-validation can be erratic, especially for small samples or smooth densities. Multiple minima in the CV curve cause problems.

Practical Recommendations:

1. Start with Silverman's robust rule for quick exploration. It's rarely catastrophic.

2. For publication/final analysis: Use Sheather-Jones (if available) or LSCV.

3. When you suspect multimodality: Cross-validation or plug-in; avoid reference rules.

4. When outliers are present: Use robust scale estimate (IQR-based) in reference rules.

5. When computational cost matters: Reference rules or plug-in (both are fast).

6. When you have no prior knowledge: Cross-validation makes the fewest assumptions.

7. Always validate visually: Plot multiple bandwidths spanning 0.5× to 2× the selected value. Features that appear across the range are robust.

Fixed-bandwidth KDE uses the same smoothing everywhere. But optimal smoothing may vary: high-density regions (modes) benefit from smaller bandwidth (capture detail), while low-density regions (tails) need larger bandwidth (reduce variance). Adaptive (variable) bandwidth methods address this.

Two Philosophies:

1. Sample-point adaptive: Each data point $x_i$ has its own bandwidth $h_i$. The kernel placed at $x_i$ has width $h_i$.

2. Design-point adaptive: The bandwidth at evaluation point $x$ depends on local density, $h(x) = h \cdot \lambda(x)$.

Both aim to achieve the same effect: narrower kernels where data is dense, wider kernels where data is sparse.

Abramson's Square-Root Rule (1982):

A widely-used adaptive approach. The bandwidth at point $x_i$ is: $$h_i = h \cdot \left(\frac{\bar{f}}{\tilde{f}(x_i)}\right)^{1/2}$$

where:

• $h$ is a global bandwidth (e.g., from Silverman's rule)
• $\tilde{f}(x_i)$ is a pilot density estimate at $x_i$ (fixed-bandwidth KDE)
• $\bar{f} = \left(\prod_{i=1}^n \tilde{f}(x_i)\right)^{1/n}$ is the geometric mean of pilot densities

Intuition: If $\tilde{f}(x_i)$ is high (dense region), then $h_i$ is small. If $\tilde{f}(x_i)$ is low (sparse region), then $h_i$ is large.

The adaptive KDE is: $$\hat{f}{\text{adapt}}(x) = \frac{1}{n} \sum{i=1}^n \frac{1}{h_i} K\left(\frac{x - x_i}{h_i}\right)$$

Normalization: The geometric mean $\bar{f}$ ensures that $\prod h_i$ remains approximately constant, preserving the overall smoothing level.

Properties of Adaptive KDE:

Advantages:

• Better captures peaks and troughs simultaneously
• Reduces bias near modes
• Reduces variance in tails
• Achieves faster convergence for some density classes

Disadvantages:

• More complex to implement
• Requires pilot estimate (introduces additional tuning)
• May not integrate to exactly 1
• Harder to analyze theoretically

Balloon Estimator:

A design-point adaptive approach: $$\hat{f}_{\text{balloon}}(x) = \frac{k}{n \cdot V_d \cdot r_k(x)^d}$$

where $r_k(x)$ is the distance to the $k$-th nearest neighbor of $x$, and $V_d$ is the volume of a unit $d$-sphere. Larger $r_k$ (sparse region) gives smaller density. This is k-NN density estimation.

Sample-Smoothing Estimator:

A sample-point adaptive approach where $h_i \propto r_k(x_i)$ (distance to $k$-th nearest neighbor of $x_i$).

Armed with theoretical understanding, let's consolidate practical implementation guidance.

Decision Flowchart:

Start
  │
  ├── Is this exploratory analysis?
  │   Yes → Use Silverman's Rule
  │   No ↓
  │
  ├── Do you suspect multiple modes?
  │   Yes → Use Cross-Validation or Plug-in
  │   No ↓
  │
  ├── Is sample size > 500?
  │   Yes → Sheather-Jones Plug-in
  │   No → Silverman's Rule (robust)
  │
  └── Always: Visualize with multiple bandwidths!

Common Mistakes to Avoid:

1. Using default without checking: Software defaults are sensible but not always appropriate for your data.

2. Ignoring visual inspection: Numbers can't tell you if modes are real or artifacts.

3. Over-relying on CV: For small samples, CV can be very noisy.

4. Reporting single bandwidth: Always report which method was used and sensitivity to choice.

5. Ignoring domain knowledge: If you know the density should be unimodal (or bimodal, etc.), use that information.

Bandwidth selection is the critical practical challenge in kernel density estimation. We've explored the major approaches in depth:

What's Next:

With univariate density estimation fully explored, we turn to the challenge of multivariate density estimation. The next page addresses the curse of dimensionality head-on, exploring techniques for estimating joint distributions and understanding how sample requirements scale with dimension.