<html><head><meta name="color-scheme" content="light dark"></head><body><pre style="word-wrap: break-word; white-space: pre-wrap;">\documentclass[12pt]{article}
\usepackage{graphicx}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{amsfonts}
%\usepackage{eucal}

  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%         DENSITY ESTIMATION              %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}

\begin{center}
{\Huge \scshape Kernel Density Estimation Tutorial}\\
\vspace{15pt}

{\large\scshape Rohan Shiloh Shah}\\
\vspace{10pt}
{\large \scshape CS644A Pattern Recognition}

\end{center}

\vspace{30pt}

In  classification and regression, the primary goal is the estimation of a  prediction function. The likelihood or conditional density is one such function; for regression $p(y|x) = p( y, x)/ \int p(y, x) d y$ and similarily for classification $p(c|x) = p( c, x)/ \sum_c p(c, x) $ where $c$ is a class label from the set of labels $\mathfrak{C}$. These are supervised learning tasks since each training example is paired with a corresponding label or annotation; for regression $y \in \mathbb{R}^n$ and for classification $c \in \mathfrak{C}$.

Given an unannotated training data set, we seek to build a model, specifically an unconditional probability density function, that delineates the essential information contained in the observation space $\mathcal{X}$. This is \emph{un-supervised learning} since it is performed in the abscence of annotations (and hence without any cost or loss function) through direct interaction with `new experiences'.

One common approach is to assume the density has a fixed parametric form and then to estimate the parameters (using a maximum likelihood approach) associated with this form; for example using a mixture model we can decompose the unknown density as follows:
\begin{equation}
\hat{f}_{(\mu, \Sigma)}(x) = \sum_{i=1}^m \zeta_i \; P_{(\mu_i, \Sigma_i)}(x), \; \zeta_i \geq 0, \; \sum \zeta_i = 1
\label{parametricdensity}
\end{equation}
where the mixing coefficients $\zeta_i$ quantify the contribution of the $i^{th}$ model in the mixture to the generation of the estimate $\hat{f}_{(\mu, \Sigma)}$. However, in many instances \emph{parametric} mixture approximations do not converge in probability to the true density. The following sections present a regularized, non-parametric estimate that is a mixture of convolved kernel functions and is asymptotically both unbiased and consistent.

\section{Non-Parametric Density Estimation}
Provided with $n$ discrete observations of a random variable
$$ \mathcal{S} = \{ x_1,x_2, \cdots , x_n \} \subseteq \mathcal{X}$$ 
all of which are identically and independently distributed (iid) according to some \emph{unknown} probability distribution $F(x)$, we seek an estimate $\hat{f}(x)$ of the true \emph{probability density function} $f(x)$. The search for $\hat{f}(x)$ is usually performed in a restricted functional space $\mathcal{H}_f$; which is infact a Reproducing Kernel Hilbert Space (see Section \ref{RKHS}). The functional space is further restricted to only those functions that are non-negative and integrate to one. The probability density function (PDF) is simply the derivative of the cumulative distribution function (CDF):
\begin{equation}
f(x) = \frac{dx}{dF(x)} \; \; \Longleftrightarrow \; \; F(x)  = \int_{-\infty}^{x} f(\psi) d\psi
\label{cdf_integral}
\end{equation}
We can rewrite \ref{cdf_integral} as a linear mapping \cite{SVD_Vapnik}:
\begin{equation}
\int_{-\infty}^{_+\infty} \mathbb{I}_{\psi &lt; x} f(\psi) d\psi = F(x)  \; \; \Longleftrightarrow \; \; \mathcal{A} f(x) = F(x)
\label{cdf_linearoperator}
\end{equation}
where both integrals in \ref{cdf_integral} and \ref{cdf_linearoperator} are vector integrations and $\mathcal{A}$ is an injective mapping from $\mathcal{H}_f$ to the Hilbert Space where $F$ is defined; $\mathcal{H}_F$. Neither $f$ or $F$ are known (whereas the operator $\mathcal{A}$ and its inverse are well defined) so we begin by estimating $F$ using samples $\mathcal{S}$ generated by the random process and then proceed to deriving $\hat{f}$ from our estimate $\hat{F}$ using an approximation of the inverse of the linear transformation $\mathcal{A}$.

The \emph{empirical distribution} at $x$ can be estimated from the data by taking the ratio of the number of samples that are less than or equal to $x$ to the total number of samples:
\begin{equation}
\hat{F}(x) = \frac{1}{n}\sum_{i=1}^n \mathbb{I}_{( -\infty, \; x_i]} (x)\label{distribution_estimate1}
\end{equation}
and is an unbiased maximum likelihood estimate that is piece-wise constant. For the density to exist, the estimated distribution $\hat{F}$ must be differentiable and hence continuous and so to smooth out the estimate $\hat{F}$; a non-linear regression (Figure \ref{densityinterpolation}) is used to approximate the distribution in the regions where training samples are unavailable; specifically the regression is performed on the set of pairs 
$$\left\{ x_i,  \hat{F}(x_i) ; \; \; \; \; i=1,\cdots , n \right\} $$
and is parametrized by the vector $\rho$ which leads to our new estimate for the distribution: $\hat{F}_{\rho}(x)$. Support Vector Regression techniques may also be used to derive accurate regressions since regularization (see Section \ref{regression_regularization}) is then possible.

Now to estimate the density function at $x$ we need to estimate the derivative of $\hat{F}_{\rho}(x)$ which can roughly be done by taking the difference between two evaluations of the distribution function at fixed lengths $+h$ and $-h$ from $x$:

\begin{figure}
\includegraphics[scale=1]{image_density_interpolation.eps}
\caption{Estimating the Density and Distribution Functions: [left] A linear interpolation on $\hat{F}_{\rho}$ between the end points of the region $\mathcal{R}= [x-h, x+h]$ gives us an estimate for the slope $\hat{f}(x)$. [right] A non-linear regression on evaluations of \ref{distribution_estimate1} at all training samples, in this case on $(x_1, \hat{F}(x_1)), (x_2, \hat{F}(x_2)), \cdots , (x_6, \hat{F}(x_6))$, yields a smooth estimate for the distribution function.} 
\label{densityinterpolation}
\end{figure}

\begin{figure}
\includegraphics[scale=0.5]{image_histogram_original_curve_01.eps}
\includegraphics[scale=0.5]{image_histogram_random_samples_01.eps}
\includegraphics[scale=0.5]{image_histogram_02_01.eps}
\includegraphics[scale=0.5]{image_histogram_01_01.eps}
\includegraphics[scale=0.5]{image_histogram_001_01.eps}
\includegraphics[scale=0.5]{image_histogram_0001_01.eps}
\caption{Parametric Frequency Histogram Estimation: The true unknown density (top left) can be estimated by taking random samples (top right, 1000 random samples) and placing them in bins of fixed length to generate a histogram. Histograms with bin-size $h=0.2$, $h=0.1$, $h=0.01$ and $h=0.001$ are shown; the bin-size or bandwidth (as well as the actual placement of the bins) is an important parameter in estimating the density function; in this case only a bin-size of $h=0.01$ is able to capture the multi-modality of the true density. In the limit as the bandwidth goes to zero the histogram will converge to the true density provided the number of samples goes to infinity.}
\label{densitybinsize}
\end{figure}

\begin{eqnarray}
\label{density_estimate1}
\nonumber
\hat{f}(x) &amp;=&amp; \mathcal{A}^{-1} \hat{F}_{\rho}(x)\\
&amp;=&amp; \frac{1}{2h} \int_{x-h}^{x+h} d\hat{F}_{\rho}(\psi)\\ \nonumber
&amp;=&amp; \frac{1}{2h}\left( \hat{F}_{\rho}(x+h) - \hat{F}_{\rho}(x-h) \right)\\ \nonumber
&amp;=&amp; \frac{1}{2h}\left( \frac{1}{n}\sum_{i=1}^n \mathbb{I}_{x_i \leq x+h} - \frac{1}{n}\sum_{i=1}^n \mathbb{I}_{x_i \leq x-h} \right)\\ \nonumber
&amp;=&amp; \frac{1}{2h}\left( \frac{1}{n}\sum_{i=1}^n \mathbb{I}_{x_i \leq x+h} - \mathbb{I}_{x_i \leq x-h} \right)\\ \nonumber
&amp;=&amp; \frac{1}{2h}\left( \frac{1}{n}\sum_{i=1}^n \mathbb{I}_{x-h \; \leq \; x_i \; \leq \; x+h} \right)\\ \nonumber
&amp;=&amp; \frac{1}{V} \times \frac{k(x)}{n}
\end{eqnarray}
where $k(x) = \sum_{i=1}^n \mathbb{I}_{x-h \; \leq \; x_i \; \leq \; x+h}$ is the number of samples that fall in the region $\mathcal{R} = [x-h, x+h]$ and $V$ is the volume of $\mathcal{R}$ which in this case is simply $(2h)^{-1}$. 

The shape of the region $\mathcal{R}$ and hence its volume $V$ can be adjusted as more random samples become availible; it has been shown \cite{duda_hart_2001} that as the regions $\mathcal{R}$ get smaller ($\lim_{n \to \infty} V = 0$) and the samples in the region increases ($\lim_{n \to \infty} k = \infty$), the estimated density $\hat{f}(x)$ will converge to $f(x)$ provided that $\lim_{n \to \infty} k/n = \infty$ (that is to say the proportion of samples falling within the region to those outside it is very small); in the limit it will be a smooth density function.

\section{Kernel Density Estimation: Parzen Windows}

In the previous section we decomposed the CDF into regions or windows $\mathcal{R}$ and estimated the PDF for each window separately. There are several ways to choose the placement (i.e. the center) of the windows; for example they can be \emph{disjoint} and cover the entire domain as is the case with Frequency Histogram Estimation (Figure \ref{densitybinsize}); the benefit of using such a method is that the structure of the windows are independant of the number of random samples (although the estimation procedure is not). Alternatively, the approach we will consider henceforth will associate a window with each random sample $x_i$ so that we have a sequence of windows $\mathcal{R}_1, \mathcal{R}_2, \mathcal{R}_3, \cdots$ centered at $x_1, x_2, x_3, \cdots$. Since the number of windows (and their volume) is now dependent on the sample size we will re-write the result (\ref{density_estimate1}) as 
\begin{equation}
\label{density_estimate2}
\hat{f}_n(x) = \frac{1}{V_n} \times \frac{k_n(x)}{n}
\end{equation}
where $V_n$ is the volume of the region $\mathcal{R}_n$ and $k_n$ is the number of samples that fall inside it; as more samples are generated the regions can either be refined (as is the case with the Parzen method) or we can reverse this logic and gradually increase the size of the regions starting with an inconspicuously small $\mathcal{R}_1$ (which is the case with the nearest-neighbor estimation method).

Let us further generalize \cite{duda_hart_2001} the definition of the region $\mathcal{R}_n$ so that we can estimate multi-variate densities; assume the windows $\mathcal{R}_n$ are d-dimensional hypercubes with edges of length $h_n$; the volume is then simply $V_n = (h_n)^d$. 
\begin{figure}
\centering
\includegraphics[scale=1.25]{image_hypercube_window_function.eps}
\caption{Let the hypercube window $\mathcal{R}_n$ have dimension $d=2$; then $k_n(x)$ is simply a count of the random samples that fall in the square with sides of length $h_n$ and centered at $x$; in this case $k_n(x)=3$.}
\label{hypercubewindowfunction}
\end{figure}
Now that the window has changed we need to revise our definition of $k_n(x)$ which was previously defined using a simple indicator function $\mathbb{I}_{x-h \leq x_i \leq x+h}$; generalising this from counting random samples within an interval to counting within a hypercube can be done by first using a \emph{window function} $\omega$ of the form:
\[
\omega(s) = \left\{\begin{array}{cl}
	1, &amp; | s_j| \leq 0.5 \; \; \forall j = 1, \cdots , d\\
	0, &amp; otherwise 
	   \end{array}\right.
\]
which defines the boundary for a unit-hypercube centered at the origin and then defining $k_n$ as follows:
$$k_n(x) = \sum_{i=1}^n \omega \left( \frac{x-x_i}{h_n} \right)$$
So $\omega \left( \frac{x-x_i}{h_n} \right) = 1$ if and only if $\frac{x-x_i}{h_n} \leq (\frac{1}{2}, \cdots , \frac{1}{2})^T$ or $x-x_i \leq (\frac{1}{2} h_n, \cdots , \frac{1}{2} h_n )^T $, in other words if $x_i$ is less than half the length of an edge of the hypercube away from $x$ in all dimensions $i=1, \cdots , d$. So $\omega \left( \frac{x-x_i}{h_n} \right)$ defines a (d+1)-dimensional hypercube of volume $(h_n)^{d}$ (since the (d+1)th dimension has edges of length $1$), centered at $x$ which counts the number of random samples that fall within the d-dimensional hypercube window $\mathcal{R}$. Finally substituting this result into \ref{density_estimate2} gives:
\begin{eqnarray}
\label{density_estimate3}
\hat{f}_n(x) &amp;=&amp; \frac{1}{V_n} \times \frac{k_n}{n}\\ \nonumber
&amp;=&amp;  \frac{1}{n} \frac{1}{(h_n)^d} \sum_{i=1}^n \; \omega \left( \frac{x-x_i}{h_n} \right) \\ \nonumber
\end{eqnarray}
The resulting estimated density is `jagged' since the window (basis) functions in the above linear combination are hypercubes (Figure \ref{hypercubewindowfunction}) with abrupt edges; reducing the width $h_n$ as more samples are generated will smooth out the estimate. Infact if $\lim_{n \to \infty} h_n = 0$ then the estimated density function is proved \cite{Fukunaga} to be asymptotically unbiased; $\lim_{n \to \infty} E(\hat{f}_n) = f$.

\subsection{Kernel Basis Functions}
\begin{figure}
\centering
\includegraphics[scale=0.75]{image_gaussian_window_function.eps}
\caption{Product Kernel  Window Functions: instead of counting the number of random samples within a hypercube centered at $x$, we can associate a single-variate kernel function with each dimension and weight the count for each random sample by the product of its kernelized distances from $x$ in each dimension. More generally a multi-variate kernel function may be used.}
\label{gaussianwindowfunction}
\end{figure}
Instead of simply counting the number of random samples that fall within a fixed volume surronding $x$, we can weight the count \cite{duda_hart_2001} for each random sample by its \emph{kernelised distance} from $x$. This can be achieved by replacing the unit hypercube window function $\omega(s)$ with a smooth, symmetric \emph{kernel density function} $K(s)$ satisfying $V_n = \int_{-\infty}^{+\infty} K(s) \; ds = 1$ and $K(s) \geq 0$ and then rewriting \ref{density_estimate3} as:
\begin{equation}
\label{density_estimate4}
\hat{f}_n(x) = \frac{1}{n}  \sum_{i=1}^n \; K_{h_n} \left( x - x_i \right)
\end{equation}
where the bandwidth $h_n$ is shifted into the definition of the kernel as the standard-deviation so that $K_{h_n}(s) = K(s/h_n)$ and the term involving the volume disappears since $V_n=1$. The gaussian kernel is most often used;
\begin{equation}
\label{multivariategaussian}
K_\Sigma(x-x_i) = \frac{1}{(2 \pi)^{n/2} | \Sigma |^{1/2}} \exp \left( - \frac{1}{2} (x - x_i)^T \Sigma^{-1} (x - x_i) \right)
\end{equation}
where $\Sigma$ is the covariance or \emph{bandwidth matrix}. The key difference between the parametric density estimate \ref{parametricdensity} and non-parametric kernel density estimation \ref{density_estimate4} is that in the former the models that define the mixture have means or centers that are estimated from the data, while the latter makes use of kernel functions that are centered at the various samples in the training data.

The use of kernel basis functions has several advantages, the most significant of which is that the resulting estimate $\hat{f}(x)$ is also a smooth density function. It has been shown \cite{Fukunaga} that provided that $\lim_{n \to \infty} h_n = 0$ and $lim_{n \to \infty} n h_n = \infty$  the estimated kernel density estimate pointwise converges in probability to the true density - this is asymptotic consistency; uniform convergence in probability is also proved under the additional condition $\lim_{n \to \infty} n(h(n))^2 = \infty$.

\begin{figure}
\centering
\includegraphics[scale=0.65]{image_kernel_density_comparison.eps}
\caption{Comparing the Gaussian and Epanetchnikov Kernels \cite{kde-matlab}: a bandwidth of $0.0215$ is used - the entropy for the Gaussian and Epanetchnikov Kernels are $0.0439$ and $0.0430$ respectively. Notice how even though the original or true density is defined only on the interval $[0,1]$ so that random samples are also only generated on this interval, the resulting estimated density extends outside this interval; this can be good if there are regions of missing values so that an implicit non-linear interpolation estimates the density in these regions; it can be bad when the estimation extends into regions for which the density is meant to be undefined.}
\label{kerneldensitycomparison}
\end{figure}

\subsection{Regularization by Convolution}
The estimate \ref{density_estimate4} can be re-written as a convolution of the kernel with the true density function; 
\begin{eqnarray}
f(s) \star K(s) &amp;\equiv&amp; \int_{-\infty}^{+\infty} K(s - x) f(x) dx \\ \nonumber
&amp;=&amp; E_{f(x)} \left( K(s - x) \right)  \\ \nonumber
&amp;\approxeq&amp; \frac{1}{n} \sum_{i=1}^n K(s - x_i) \\ \nonumber
&amp;=&amp; \hat{f} (s)\\ \nonumber
\end{eqnarray}
and so in a sense, the kernel density estimate is approximately a deconvolution from the true density; in \cite{duda_hart_1973} we see that in the limit as the number of random samples approaches infinity, $\hat{f}(s)$ converges to $f(s) \star K(s)$. We can also write the density estimate as a smoothing (we refer to smooth and smoothing in two contexts; \emph{smooth functions} have continuous derivatives and \emph{smoothing operations} remove localized features) convolution of the \emph{impulsive density function} (which assigns a probability mass of $\frac{1}{n}$ to each of the $n$ random samples; $\hat{f}^i(x) = 1/n$) with the kernel function:
\begin{eqnarray}
 \hat{f}^i(s) \star K(s) &amp;=&amp; \sum_{i = 1}^{n} \hat{f}^i(x_i) \; K( s - x_i)\\ \nonumber
 &amp;=&amp; \sum_{i = 1}^{n} \frac{1}{n} \; K( s - x_i)\\ \nonumber
 &amp;=&amp; \hat{f}(s)\\ \nonumber
\end{eqnarray}
Essentially the convolution is a regularization of the estimate through the addition of `smoothing' noise in the regions where the impulsive density is undefined.

It is interesting to note that the Parzen Method has been shown \cite{parzen_regularisation} to be equivalent to a Regularized Least Squares Method (or Tikhonov Regularization) where the regularizing functional is taken to be the norm of the resulting estimated density. One benefit \cite{regularisationbyconvolution} of `regularizing by convolving' is that there are no explicit regularization parameters that need to be estimated and then re-trained.
\subsection{Bandwidth Matrix Selection}
Changes in the bandwidth $h_n$ of the Kernel function, can severely effect the resulting estimate; choosing a large bandwidth will produce biased estimates that hide localised features whereas smaller bandwidths will increase the estimates variability by introducing sharp modulations. The simplest choice of bandwidth is as some function of $n$, for example $k_n = 1/n$, so that as the number of samples increases the kernel gets smaller. A more complex method minimizes the integrated mean squared error between $\hat{f}_n$ and $f_n$ with respect to $k_n$. In Figure \ref{bandwidthselection}, three bandwidth selection methods are compared. 

\subsection{Computation and Sparsity}
To form an estimation of the kernel density using Parzen Windows, \ref{density_estimate4} (or \ref{density_estimate3}) must be evaluated for \emph{all} $n$ random samples. The sum in \ref{density_estimate4} (or \ref{density_estimate3}) is another $n$ evaluations of the window function, which in itself is a function of $d$ the number of dimensions. In subsequent sections on support vector density estimation, an equally robust estimate involving considerably less computation will be considered. 

\begin{figure}
\begin{center}
\includegraphics[scale=0.3]{image_gaussian_density_estimation.eps}
\end{center}
\caption{Kernel Bandwidth Selection \cite{kde-matlab}: Three bandwitdh selection methods are compared: the entropy and bandwidth of each are: ROT (0.1037, 0.0761), LCV (0.0439, 0.0215), HALL(0.0491, 0.0268). Notice the positive correlation between the entropy and the bandwidth; intuitively as the bandwidth increases the kernel density estimate gets smoother and closer to a uniform distribution which has the maximal entropy.}
\label{bandwidthselection}
\end{figure}

\bibliographystyle{annotate}
\bibliography{bibliography}

\end{document}


</pre></body></html>