Kernel smoothing

Kernel smoothing is a non parametric estimation method of the probability density function of a distribution.

Introduction

Let X be a random variable with probability density function p. Given a sample of independent observations x_1, ..., x_n of X and any point x \in \Rset, the kernel smoothing estimator provides an approximation \widehat{p}(x) of p(x).

In dimension 1, the kernel smoothed probability density function \widehat{p} has the following expression, where K is the univariate kernel, n the sample size and (X_1, \cdots, X_n) \in \Rset^n the univariate random sample with \forall i, \, \, X_i \in \Rset ([wand1994] eq. 2.2 page 12):

(1)\widehat{p}(x) = \frac{1}{nh}\sum_{i=1}^{n} K\left(\frac{x-X_i}{h}\right)

The kernel K is a function satisfying:

\int K(x)\, dx=1.

The parameter h is called the bandwidth.

Multivariate kernels

In dimension d>1, the kernel may be defined as a product kernel K_d, as follows where \vect{x} = (x_1, \cdots, x_d)\in \Rset^d ([wand1994] eq. 2.2 page 91):

K_d(\vect{x}) = \prod_{j=1}^{d} K(x_j).

The kernel smoothed probability density function in dimension d is:

\widehat{p}(\vect{x}) = \frac{1}{n \prod_{j=1}^{d}h_j} \sum_{i=1}^{n} K_d\left(\frac{x_1 - (X_{i})_1 }{h_1}, \dots, \frac{x_d - (X_{i})_d}{h_d}\right)

where (\vect{X}_1, \cdots, \vect{X}_n) is the d-variate random sample which components are denoted \vect{X}_i = (X_{i1}, \dots, X_{id}).

In the multivariate case, the bandwidth is the vector \vect{h} = (h_1, \cdots, h_d).

Asymptotic error and asymptotically optimal bandwidth

The quality of the approximation may be controlled by the AMISE (Asymptotic Mean Integrated Square error) criteria defined as follows. By definition, the mean squared error (MSE) at point \vect{x} is (see [scott2015] eq. page 40, [wand1994] pages 14-21):

\operatorname{MSE}(\widehat{p}, \vect{x}) = \mathbb{E}_\vect{X}\left[\left(\widehat{p}(\vect{x}) - p(\vect{x})\right)^2 \right].

It can be proved that the mean squared error is the sum of the variance and the squared bias:

\operatorname{MSE}(\widehat{p}, \vect{x}) = \operatorname{Var}_\vect{X}\left[\widehat{p}(\vect{x})\right] + \left(\operatorname{Bias}\left(\widehat{p}(\vect{x})\right)\right)^2

where the bias is:

\operatorname{Bias}\left(\widehat{p}(\vect{x})\right) = \mathbb{E}_\vect{X}\left[\widehat{p}(\vect{x})\right] - p(\vect{x}).

The MSE depends on the point where the density is evaluated. The mean integrated squared error (MISE) is (see [scott2015] eq. page 41):

\operatorname{MISE}(\widehat{p}) = \mathbb{E}_\vect{X}\left[||\widehat{p} - p||^2_{L_2}\right] = \int \, \operatorname{MSE}(\widehat{p}, \vect{x}) \, d\vect{x} \\

Finally, the asymptotic mean integrated squared error (AMISE), denoted \operatorname{AMISE}(\widehat{p}) is defined as the two first terms in the Taylor series expansion of \operatorname{MISE}(\widehat{p}) when n tends to infinity.

The quality of the estimation essentially depends on the value of the bandwidth h. In dimension 1, the bandwidth that minimizes the AMISE criteria is (see [wand1994] eq 2.13 page 22):

(2)h_{\operatorname{AMISE}}(K) = \left( \frac{R(K)}{\mu_2(K)^2 R\left(p^{(2)}\right)}\right)^{\frac{1}{5}}n^{-\frac{1}{5}}

where the rugosity of the kernel is:

R(K) = \int K(\vect{x})^2\, d\vect{x}

and the second raw moment of the kernel is:

\mu_2(K) = \int \vect{x}^2K(\vect{x})\, d\vect{x} = \sigma_K^2.

In the equation (2), the expression R\left(p^{(2)}\right) is the rugosity of the second derivative of the density probability function p that we wish to approximate. Since, by hypothesis, the true density p is unknown, its second derivative is also unknown. Hence the equation (2) cannot be used directly to compute the bandwidth.

We have ([wand1994] page 67):

R\left(p^{(r)}\right) = (-1)^r\Psi_{2r}

where:

\Psi_r(p) = \int p^{(r)}p(x)\, dx = \mathbb{E}_\vect{X}\left[p^{(r)}\right].

Therefore:

(3)h_{\operatorname{AMISE}}(K) = \left( \frac{R(K)}{\mu_2(K)^2\Psi_4(p)}\right)^{\frac{1}{5}}n^{-\frac{1}{5}}

Several methods exist to evaluate the optimal bandwidth h_{\operatorname{AMISE}}(K) based on different approximations of \Psi_4:

  • Silverman’s rule in dimension 1,

  • the plug-in bandwidth selection,

  • Scott’s rule in dimension d.

Efficiency of a kernel

Let K be a kernel. We may be interested if a particular kernel may be able to reduce the estimation error. The efficiency of a kernel is (see [scott2015] page 151):

\operatorname{eff}(k) = \frac{\sigma_k R(k)}{\sigma_{k_E} R(k_E)}

where k_E is Epanechnikov’s kernel. The AMISE error is proportional to the efficiency (see [scott2015] eq. 6.25 page 151):

\operatorname{AMISE} \propto \operatorname{eff}(k)

The next table presents several kernels available in the library and their associated variance, rugosity and efficiency. We see that the best kernel is Epanechnikov’s kernel. We also see that there is not much difference between the different kernels. This is one of the reasons why the normal kernel is often used.

Kernel

\operatorname{Var}(k)

R(k)

\sigma_k R(k)

\operatorname{eff}(k)

Epanechnikov

0.2000

0.6000

0.2683

100.00 %

Biweight

0.1429

0.7143

0.2700

99.39 %

Quartic

0.1429

0.7143

0.2700

99.39 %

Triweight

0.1111

0.8159

0.2720

98.67 %

Triangular

0.1667

0.6667

0.2722

98.59 %

Normal

1.0000

0.2821

0.2821

95.12 %

Uniform

0.3333

0.5000

0.2887

92.95 %

Logistic

3.2899

0.1667

0.3023

88.76 %

Table 1. Efficiency of several order 2 kernels.

Silverman’s rule (dimension 1)

In this section, we consider a random variable i.e. d = 1. If the density p is normal with standard deviation \sigma, then the term \Psi_4 can be exactly evaluated. By definition, the Silverman rule for the bandwidth is the optimal bandwidth of the AMISE criteria when the true density p is normal (see [silverman1986] page 45):

(4)h^{Silver}(K) := h^{p = normal}_{\operatorname{AMISE}}(K) = \left( \frac{8\sqrt{\pi} R(K)}{3\mu_2(K)^2}\right)^{\frac{1}{5}} \sigma n^{-\frac{1}{5}}.

The Silverman rule is based on the hypothesis that the true density p is close to the normal density, even if the density p is not necessarily normal.

The equation (4) is accurate when the density is not far from a normal one. In the special case where we use the normal kernel, the Silverman rule is (see [silverman1986] eq 3.28 page 45):

(5)h^{Silver}(\mathcal{N}) = 1.06 \sigma n^{-\frac{1}{5}}.

Choice for the standard deviation

We estimate the standard deviation \sigma by its sample counterpart \hat{\sigma}, evaluated from the sample (x_1, \dots, x_n):

(6)h^{Silver}(K) = \left( \frac{8\sqrt{\pi} R(K)}{3\mu_2(K)^2}\right)^{\frac{1}{5}} \hat{\sigma} n^{-\frac{1}{5}}

The estimator \hat{\sigma} of the true standard deviation can be estimated using the sample standard deviation based on the sample (x_1, \dots, x_n). This is:

\hat{\sigma} = \sqrt{\frac{1}{n - 1} \sum_{i = 1}^n (x_i - \bar{x})^2 }

where \bar{x} is the sample mean:

\bar{x} = \frac{1}{n} \sum_{i = 1}^n x_i.

Another method is to use the standardized interquartile range ([wand1994] page 60):

\widehat{\sigma}_{\operatorname{IQR}} = \frac{\widehat{q}(3/4) - \widehat{q}(1/4)}{\Phi^{-1}(3/4) - \Phi^{-1}(1/4)}

where \Phi^{-1} is the quantile function of the standard normal distribution and \widehat{q}(3/4) and \widehat{q}(1/4) are the sample quartiles at levels 75% and 25% respectively. The previous estimator is robust against outliers that might be in the sample.

Plug-in bandwidth selection method (dimension 1)

The plug-in bandwidth selection method improves the estimation of the rugosity of the second derivative of the density. Instead of making the Gaussian assumption, the method uses a kernel smoothing method in order to make an approximation of higher derivatives of the density. This method is due to [sheather1991] who used ideas from [park1990]. The algorithm is presented in [wand1994], page 74 under the “Solve the equation rule” name. The implementation uses ideas from [raykar2006], but the fast selection is not implemented.

The equation (3) requires the evaluation of the quantity \Psi_4. We use the estimator \hat{\Psi}_r of \Psi_r, using a kernel density estimator of the r-th derivative of the density. The estimator is (see [wand1994] page 64):

(7)\hat{\Psi}_r(K) = \frac{1}{n}\sum_{i=1}^{n} \widehat{p}^{(r)}(X_i)

where \hat{\Psi}_r(K) is the estimator based on the kernel K. Deriving equation (1) leads to:

(8)\widehat{p}^{(r)}(x) = \frac{1}{nh^{r+1}}\sum_{i=1}^{n} K^{(r)}\left(\frac{x-X_i}{h}\right)

and then the estimator \hat{\Psi}_r(h) is defined as:

(9)\hat{\Psi}_r(h) = \frac{1}{n^2h^{r+1}}\sum_{i=1}^{n}\sum_{j=1}^{n} K^{(r)}\left(\frac{X_i-X_j}{h}\right)

We note that \hat{\Psi}_r(h) depends of the parameter h which can be taken in order to minimize the Asymptotic Mean Square Error (AMSE) criteria evaluated between \Psi_r and \hat{\Psi}_r(h). The optimal parameter h is:

(10)h^{(r)}_{\operatorname{AMSE}} = \left(\frac{-2K^{(r)}(0)}{\mu_2(K) \Psi_{r+2}}\right)^{\frac{1}{r+3}}n^{-\frac{1}{r+3}}

The previous equation states that the bandwidth h^{(r)} required to compute \widehat{p}^{(r)} depends on \Psi_{r+2}. But to compute \Psi_{r+2}, we need h^{(r + 2)}, etc. The goal of the method is to break this infinite set of equations at some point by providing a pilot bandwidth. The \ell-stage plug-in bandwidth method uses \ell different intermediate bandwidths before evaluating the final one.

In this document, we present the two stage solve-the-equation plug-in method.

  • The equation (3) defines h_{\operatorname{AMISE}}(K) as a function of \Psi_4. Let t be the function defined by the equation:

    (11)