80  Kernel Density Estimation

80.1 General Definition

Kernel Density Estimation (Rosenblatt 1956; Parzen 1962) is a non-parametric method for estimating the probability density function of a random variable. If \((x_1, x_2, …, x_n)\) is an independently and identically distributed sample from any density function \(f\) then the Kernel Density Estimator is

\[ \hat{f}_h(x) = \frac{1}{n} \sum_{i=1}^{n} K_h \left( x - x_i \right) = \frac{1}{n h} \sum_{i=1}^{n} K \left( \frac{x - x_i}{h} \right) \]

where \(K_h()\) is the Kernel function with bandwidth \(h > 0\). The function \(K_h(z) \geq 0\) and has a unit integral \(\left( \int_\mathbb{R} K_h(z)\text{d}z = 1 \right)\) and a mean of 0. In other words, each sample point \(x_i\) is replaced by a Kernel function and the estimated Kernel Density is equal to the sum of all the Kernel functions.

80.1.1 Horizontal axis of the Kernel Density Plot

The horizontal axis represents the values of the variable under investigation.

80.1.2 Vertical axis of the Kernel Density Plot

The vertical axis represents the estimated density.

80.2 Evaluation Distance

80.2.1 Definition

\[ z_{ij} = \frac{1}{h} \left( v_j - x_i \right) \]

This is the distance of observation \(x_i\) from the point \(v_j\) at which the density will be evaluated.

80.3 Local Density

80.3.1 Definition

\[ \hat{f}\left( v_j \right) = \frac{1}{nh} \sum_{i=1}^{n} K \left[ z_{ij} \right] \]

for \(j=1, 2, …, m\) where

  • \(n\) is the number of observations in the sample
  • \(m\) is the number of equally spaced evaluation points
  • \(h\) is the bandwidth of the sliding window
  • \(z_{ij}\) is the evaluation distance
  • \(K\) is the weight function or Kernel

80.4 Rectangular Kernel

80.4.1 Definition

\[ \begin{aligned}\left| z_{ij} \right| \leq 1 &\Rightarrow K_{h,R} \left[ z_{ij} \right] = 0.5 \\\left| z_{ij} \right| > 1 &\Rightarrow K_{h,R} \left[ z_{ij} \right] = 0\end{aligned} \]

where

\[ z_{ij} = \frac{1}{h} \left( v_j - x_i \right) \]

and \(h\) is the bandwidth of the sliding window.

80.5 Triangular Kernel

80.5.1 Definition

\[ \begin{aligned}\left| z_{ij} \right| \leq 1 &\Rightarrow K_{h,T} \left[ z_{ij} \right] = 1 - \left| z_{ij} \right| \\\left| z_{ij} \right| > 1 &\Rightarrow K_{h,T} \left[ z_{ij} \right] = 0\end{aligned} \]

where

\[ z_{ij} = \frac{1}{h} \left( v_j - x_i \right) \]

and \(h\) is the bandwidth of the sliding window.

80.6 Gaussian Kernel

80.6.1 Definition

\[ K_{h,G} \left[ z_{ij} \right] = \frac{1}{\sqrt{2 \pi}} e^{-\frac{z_{ij}^2}{2} } \]

where

\[ z_{ij} = \frac{1}{h} \left( v_j - x_i \right) \]

and \(h\) is the bandwidth of the sliding window.

80.7 Epanechnikov Kernel (Epanechnikov 1969)

80.7.1 Definition

\[ \begin{aligned}\left| z_{ij} \right| \leq \sqrt{5} &\Rightarrow K_{h,E} \left[ z_{ij} \right] = \frac{3}{4 \sqrt{5}} \left( 1 - \frac{z_{ij}^2}{5} \right) \\\left| z_{ij} \right| > \sqrt{5} &\Rightarrow K_{h,E} \left[ z_{ij} \right] = 0\end{aligned} \]

where

\[ z_{ij} = \frac{1}{h} \left( v_j - x_i \right) \]

and \(h\) is the bandwidth of the sliding window.

80.8 R Module

80.8.1 Public website

The Kernel Density Plot is available on the public website:

80.8.2 RFC

The Kernel Density Plot is also available in RFC under the “Descriptive / Kernel Density Estimation”.

To compute the Kernel Density Plot on your local machine, the following script can be used in the R console:

x <- runif(400)
par1 = 0
par2 = 'no'
par3 = 512
ylab = 'Density'
xlab = 'Value'
if (par1 == 0) bw <- 'nrd0'
plot(density(x, bw=bw, kernel='gaussian', na.rm=TRUE), main='Gaussian Kernel', xlab=xlab, ylab=ylab)
grid()

plot(density(x, bw=bw, kernel='epanechnikov', na.rm=TRUE), main='Epanechnikov Kernel', xlab=xlab, ylab=ylab)
grid()

plot(density(x, bw=bw, kernel='rectangular', na.rm=TRUE), main='Rectangular Kernel', xlab=xlab, ylab=ylab)
grid()

plot(density(x, bw=bw, kernel='triangular', na.rm=TRUE), main='Triangular Kernel', xlab=xlab, ylab=ylab)
grid()

plot(density(x, bw=bw, kernel='biweight', na.rm=TRUE), main='Biweight Kernel', xlab=xlab, ylab=ylab)
grid()

plot(density(x, bw=bw, kernel='cosine', na.rm=TRUE), main='Cosine Kernel', xlab=xlab, ylab=ylab)
grid()

plot(density(x, bw=bw, kernel='optcosine', na.rm=TRUE), main='Optcosine Kernel', xlab=xlab, ylab=ylab)
grid()

To compute the Kernel Density Plot, the R code uses the density function with a parameter that defines the bandwidth and another parameter for the kernel.

80.9 Purpose

The Kernel Density Plot is used to visualize the empirical density of the variable under investigation. Sometimes the Kernel Densities are used to transform the data before computing some statistic which is used in subsequent analysis. The transformation which is induced can be interpreted as having a smoothing effect.

80.10 Bandwidth Selection

The bandwidth \(h\) controls the smoothness of the estimated density:

  • small \(h\): more detail, but noisier estimates (higher variance)
  • large \(h\): smoother estimates, but potentially oversmoothed (higher bias)

In practice, software often uses a default rule such as Silverman’s rule of thumb (Silverman 1986) or similar plug-in rules. These defaults are usually reasonable as a starting point, but it is good practice to inspect how conclusions change under slightly smaller and larger bandwidths.

80.11 Pros & Cons

80.11.1 Pros

The Kernel Density Plot has the following advantages:

  • Unlike the Histogram, it does not require to specify the number of bins.
  • It is easily interpreted and provides a lot of (detailed) information about the empirical density function.

80.11.2 Cons

The Kernel Density Plot has the following disadvantages:

  • It requires a bandwidth parameter to be specified. In most cases, however, the software will use a fairly appropriate bandwidth.
  • Not all software packages allow the computation of the Kernel Density Plot.

80.12 Example

Interactive Shiny app (click to load).
Open in new tab

The analysis shows the Gaussian Kernel Density Plot for the monthly marriages time series. It can be concluded that the time series exhibits a “camel shape” which indicates a bimodal density (there are two local maxima). The left hump of the camel shape is higher than the right one -- this implies that the left hump is displayed in the Table as a “global” maximum.

80.13 Task

Try to explain why the marriages time series is bimodal. Also, compare this result with the Kernel Density Plot for the monthly divorces time series (same country, same period). Why do both time series seem to have different distributions?