---
# Examples of 2-D classifications using QDA and KDE methods.
---

In this script, we show the differences in classification performances between the 
quadratic discriminant analysis (QDA) and the kernel density estimation (KDE) methods.

Both methods are applied to a series of 2-D datasets containing 2 classes with an increasing level 
spatial complexity. 

We estimate the probability distribution function (PDF) of an observation x using a mixture of PDF weighted with their 
<i>a priori</i> class probabilities $P(C_{i})$:

<blockquote>  $P(\bf{x}) = \sum_{i=0}^{N}P(\bf{x}|C_{i}) P(C_{i})$ </blockquote>

The two methods differ in the way they compute the $P(\bf{x}|C_{i})$ term.

The QDA method approximates the PDF of each class as a multi-normal gaussian distribution which parameters are 
adjusted to fit the observations (the datasets).

The KDE method approximates the PDF by convolving the datasets with spatial kernels. Hence, 
there is no <i>a priori</i> model to guide us; this is a non-parametric model.


## Discriminant methods

Once the $P(\bf{x}|C_{i})$ values are evaluated, discriminant methods are used to locate the influence zone 
of each class i over a spatial grid. This corresponds to regions where 

<blockquote> $P(\bf{x}|C_{i}) P(C_{i}) > P(\bf{x}|C_{j \ne i}) P(C_{j \ne i})$ </blockquote>

or 

<blockquote> $h_{i}(\bf{x}) > h_{j \ne i}(\bf{x})$ </blockquote>

with the discriminant function $h_{i}(\bf{x})$ defined as

<blockquote>  $h_{i}(\bf{x}) = \ln P(\bf{x}|C_{i}) P(C_{i})$ </blockquote>
 

## QDA method
The discriminant function for the 2-D gaussian PDF is:

<blockquote>  $h_{i}(\bf{x}) =\frac{1}{2}\ln |\Sigma_{i}| -\frac{1}{2}(\bf{x}-\bf{\mu_{i}})^\top \Sigma_{i}^{-1} (\bf{x}-\bf{\mu_{i}}) + \ln P(C_{i})$  </blockquote>

with 

<blockquote>  The observation:   $\bf{x}=[x_{1} x_{2}]^\top$ </blockquote> 

<blockquote>  The origin:   $\bf{\mu_{i}} = [\mu_{1} \mu_{2}]^\top$  </blockquote> 

<blockquote>  The covariance matrix: $\Sigma_{i} = \begin{pmatrix} \sigma_{x_{1}}^2 & \sigma_{x_{1,2}} \\ \sigma_{x_{1,2}} & \sigma_{x_{2}}^2 \end{pmatrix}$ </blockquote>

The parameters $\bf{\mu}$ and $\Sigma$ are estimated from the dataset.

## KDE method
The probability density $P(x|C_{i})$ can be modeled as the convolution of the data distribution with a kernel $K(u)$ designed to
smooth the estimated PDF

<blockquote>  $P(x|C_{i}) = \large{ \frac{1}{N_{i}h}\sum_{i=1}^{N_{i}} K(\frac{x-x_{i}}{h}) }$ </blockquote> 
                                 
where $N_{i}$ is the number of data points belonging to class i. The level of smoothing is controled by the bandwidth parameter <i>h</i>. 
<br><br>    


Gaussian and Epanechikov kernels are popular kernels. Their 2-D expressions are:
    
<blockquote>  Gaussian: $K(u)= \frac{1}{2\pi} \exp [-\frac{\|u\|^{2}}{2}]$ </blockquote>
 
<blockquote>  Epanechikov: $K(u)= \frac{15}{8\pi}(1-\|u\|^{2})$ for $\|u\| \in [0, 1]$ </blockquote>


The discriminant function becomes 

<blockquote> 
$
\begin{align}
P(x|C_{i})P(C_{i}) & = \frac{1}{N_{i}h}\sum_{i=1}^{N_{i}} K(\frac{x-x_{i}}{h}) P(C_{i}) \\
& =  \frac{1}{Nh}\sum_{i=1}^{N_{i}} K(\frac{x-x_{i}}{h})
\end{align}
$                
</blockquote>

                                 



```python
print(__doc__)

# Author: Pierre Gravel <pierre.gravel@iid.ulaval.ca>
# License: BSD

import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from matplotlib import colors
from sklearn.discriminant_analysis import QuadraticDiscriminantAnalysis

from sklearn import cluster, datasets
from sklearn.preprocessing import StandardScaler
from sklearn.neighbors import KernelDensity

import seaborn as sns
sns.set()

# Used for reproductibility of the results
np.random.seed(42)
```

    Automatically created module for IPython interactive environment


# A few useful functions

Let us first define several functions that will be used to generate each classification example.

## Data normalization
Normalize dataset for easier parameter selection for the classification methods.


```python
def normalize_dataset(X, nbins):
    X = StandardScaler().fit_transform(X)
    x_min = 1.1 * X[:,0].min()
    x_max = 1.1 * X[:,0].max()
    y_min = 1.1 * X[:,1].min()
    y_max = 1.1 * X[:,1].max()

    # Generate a spatial grid where the PDF will be evaluated locally.
    xx, yy = np.meshgrid(np.linspace(x_min, x_max, nbins), np.linspace(y_min, y_max, nbins))
    Xgrid = np.c_[xx.ravel(), yy.ravel()]
    
    return (X, Xgrid, xx, yy, x_min, x_max, y_min, y_max)
```

## KDE classification method
The KDE method is used to estimate the local probability density function associated to each class. 
Both PDF and their <i>a priori</i> probabilities are then used to find 
the local winning class at each point of the grid. 

Because the gaussian kernel has an infinite spatial extent, the winning class is defined everywhere. 
This is not the case with the Epanechikov kernel which drops to zero at distances larger than the
bandwidth from any data point. This problem can be alleviated by using larger bandwidths but at the price
of losing spatial resolution as we will see in the examples below.   


```python
def KDE_classifications(h_gaussian, h_epanechikov, X, y, Xgrid, nbins, prob_C):

    pdf_g = np.zeros((nbins,nbins,2))
    pdf_e = np.zeros((nbins,nbins,2))
    
    kde_g = KernelDensity(bandwidth=h_gaussian, kernel='gaussian', leaf_size=100) 
    kde_e = KernelDensity(bandwidth=h_epanechikov, kernel='epanechnikov', leaf_size=100) 
   
    # Normalisation for each data cluster
    for i in range(2):
        kde_g.fit(X[y == i])  
        Z = np.exp(kde_g.score_samples(Xgrid))
        pdf_g[:,:,i] = Z.reshape([nbins, nbins]) * prob_C[i]

        kde_e.fit(X[y == i])  
        Z = np.exp(kde_e.score_samples(Xgrid))
        pdf_e[:,:,i] =  Z.reshape([nbins, nbins]) * prob_C[i]
        
    # Find winning class
    C_g = np.argmax(pdf_g, axis=2)
    C_e = np.argmax(pdf_e, axis=2)
    
    # Find the zones on the grid where none of the Epanechikov PDF is defined.  
    eps = np.finfo(float).eps
    C_undef = np.logical_and((pdf_e[:,:,0]<eps), (pdf_e[:,:,1]<eps))
    
    # Assign a class 2 to those pixels
    C_e[C_undef] = 2  
    
    return (C_g, C_e)
        
```

## Classification results displaying function 


```python
def display_classification_results(C_g, C_e, C_qda, X, y, Xgrid, figName):

    cm_dim = ListedColormap(['#00BB00', '#0000BB', '#000000'])
    cm_bright = ListedColormap(['#00FF00', '#0000FF'])
    
    fig, axs = plt.subplots(1,3, figsize=(15,5), sharex=True, sharey=True) 
    for k in range(3):
        ax = axs.ravel()[k]

        if (k==0):
            # QDA classification results
            ax.pcolormesh(xx, yy, C_qda, cmap=cm_dim, norm=colors.Normalize(0., 2.), zorder=0) 
            ax.scatter(X[:, 0], X[:, 1], c=y, cmap=cm_bright, edgecolors='w')
            ax.set_title('QDA', fontsize=18, color='k')
        elif (k==1):
            # Gaussian KDE classification results            
            ax.pcolormesh(xx, yy, C_g, cmap=cm_dim, norm=colors.Normalize(0., 2.), zorder=0) 
            ax.scatter(X[:, 0], X[:, 1], c=y, cmap=cm_bright, edgecolors='w')
            ax.set_title('KDE (Gaussian Kernel)', fontsize=18, color='k')
        else:
            # Epanechikov KDE classification results    
            # A black color is assigned to zones where the PDF are undefined
            ax.pcolormesh(xx, yy, C_e, cmap=cm_dim, norm=colors.Normalize(0., 2.), zorder=0) 
            ax.scatter(X[:, 0], X[:, 1], c=y, cmap=cm_bright, edgecolors='w')
            ax.set_title('KDE (Epanechikov Kernel)', fontsize=18, color='k')

        ax.set_xlim([x_min,x_max])
        ax.set_ylim([y_min,y_max])
        ax.set_xlabel('$x_{1}$', fontsize=14)
        ax.set_ylabel('$x_{2}$', fontsize=14, rotation=0) 
        ax.xaxis.set_label_coords(0.5, -0.05) 
        ax.yaxis.set_label_coords(-0.05, 0.5) 
        ax.set_aspect('equal')
        ax.set_xticks([])
        ax.set_yticks([])
        
        ax.grid(None)  
        ax.grid(None)  
    
    fig.tight_layout()
        
    plt.savefig(figName + '.png')
    plt.savefig(figName + '.pdf')
```

# Comparison of classification results obtained with QDA and KDE methods 
## Initialisations


```python
# Number of 2-D points per dataset.
npts = 200

# A priori probability of the clusters
prob_C = np.array([0.5, 0.5]) 

# Number of data points allocated to each cluster 
nc = (npts*prob_C).astype(int)

# Number of spatial bins used in each spatial direction
nbins = 500

# Quadratic discriminant analysis
qda = QuadraticDiscriminantAnalysis(store_covariance=True)
```

## Two isotropic gaussians distributions with identical variance.
This is an example where the QDA produces a linear boundary. This is shown in the first panel of the figure below.
There is no <i>a priori</i> reason for the KDE method to produce linear boundaries. 

Small values of KDE bandwidths were selected to generate noisy boundaries in the next two panels. 

The third panel shows a limitation of using the Epanechikov kernel; the PDF associated to a cluster 
does not extend beyond a bandwidth distance from its outermost data points. The black zones are such regions 
where no PDF can be computed and thus no winning class can be found. 


```python
X, y = datasets.make_blobs(n_samples=nc)

# Normalize dataset
(X, Xgrid, xx, yy, x_min, x_max, y_min, y_max) = normalize_dataset(X, nbins)

```


```python
# QDA classification
C_qda = qda.fit(X, y).predict(Xgrid) 
C_qda = C_qda.reshape(xx.shape)

# KDE classifications
h_gaussian = .01
h_epanechikov = .3

(C_g, C_e) = KDE_classifications(h_gaussian, h_epanechikov, X, y, Xgrid, nbins, prob_C)

# Compare classification results
figName = '2D_KDE_classif_gaussian_1'
display_classification_results(C_g, C_e, C_qda, X, y, Xgrid, figName)
```

What happens when the bandwidths are much larger? In this case, the PDF associated to each cluster is smoothed out 
and extends beyond the panels' field of view. This creates linear boundaries as can be seen in the figure below. 
The Epanechikov-resulting PDF are now defined everywhere and the local winning class can be evaluated everywhere; 
there are no black zones left. 


```python
# KDE classifications
h_gaussian = 5
h_epanechikov = 10

(C_g, C_e) = KDE_classifications(h_gaussian, h_epanechikov, X, y, Xgrid, nbins, prob_C)

# Compare classification results
figName = '2D_KDE_classif_gaussian_2'
display_classification_results(C_g, C_e, C_qda, X, y, Xgrid, figName)
```

## Noisy concentric circles 

This example shows more complex boundaries. The influence zones of the QDA and the gaussian-KDE are similar.

The Epanechikov-KDE results are very different for the same reason as before; its bandwidth was
selected for that reason. Notice how the distance between the black boundaries and the nearest points are 
equal to the bandwidth, as expected. 


```python
X,y = datasets.make_circles(n_samples=npts, factor=.5, noise=.05)

# Normalize dataset
(X, Xgrid, xx, yy, x_min, x_max, y_min, y_max) = normalize_dataset(X, nbins)

```


```python
C_qda = qda.fit(X, y).predict(Xgrid) 
C_qda = C_qda.reshape(xx.shape)

h_gaussian = .25
h_epanechikov = .5

(C_g, C_e) = KDE_classifications(h_gaussian, h_epanechikov, X, y, Xgrid, nbins, prob_C)

# Compare classification results
figName = '2D_KDE_classif_circles_1'
display_classification_results(C_g, C_e, C_qda, X, y, Xgrid, figName)
```

The following figure shows the results when the Epanechikov-KDE bandwiths is adjusted to clear the black zones. 
All 3 methods now produce similar-looking influence zones.


```python
h_gaussian = .25
h_epanechikov = 1.2

(C_g, C_e) = KDE_classifications(h_gaussian, h_epanechikov, X, y, Xgrid, nbins, prob_C)

# Compare classification results
figName = '2D_KDE_classif_circles_2'
display_classification_results(C_g, C_e, C_qda, X, y, Xgrid, figName)
```

# Noisy crescents

This example shows more complex boundaries. This is a case where the QDA cannot well delineate the influence zones since 
the boundaries cannot be linear nor quadratic. The gaussian-KDE results are excellent and so are the Epanechikov-KDE results except 
for the black zones.



```python
X,y = datasets.make_moons(n_samples=npts, noise=.15)

# Normalize dataset
(X, Xgrid, xx, yy, x_min, x_max, y_min, y_max) = normalize_dataset(X, nbins)

```


```python
C_qda = qda.fit(X, y).predict(Xgrid) 
C_qda = C_qda.reshape(xx.shape)

h_gaussian = .25
h_epanechikov = .5


(C_g, C_e) = KDE_classifications(h_gaussian, h_epanechikov, X, y, Xgrid, nbins, prob_C)

# Compare classification results
figName = '2D_KDE_classif_crescents_1'
display_classification_results(C_g, C_e, C_qda, X, y, Xgrid, figName)
```

Can we improve again the Epanechikov-KDE results by using a larger bandwidth as we did before? 

The next figure shows that it is not the case. The Epanechikov-kernel bandwidth was increased until no more black
zones were found in the third panel. Unfortunately, this led to many data poits being misclassified. 

With this complex dataset, the best classification was obtained using the KDE method with a gaussian kernel. 


```python
h_gaussian = .25
h_epanechikov = 2.

(C_g, C_e) = KDE_classifications(h_gaussian, h_epanechikov, X, y, Xgrid, nbins, prob_C)

# Compare classification results
figName = '2D_KDE_classif_crescents_2'
display_classification_results(C_g, C_e, C_qda, X, y, Xgrid, figName)
```

# Unstructured random distribution

This example shows a 2-class distribution with no natural boundary. Once again, this is a case where the 
QDA method performs poorly since no linear or quadratic boundaries can be found. 

The gaussian-KDE results display overfitting to the noisy data. This 
shows that kernel methods can nevertheless be used to approximate complex distributions since most green and blues points 
are well separated. The spatial resolution is high.

The Epanechikov-KDE results, on the other hand, have a lower spatial resolution. They were obtained using the 
minimum bandwith value to prevent black zones. The results are less impressive than the previous ones 
eventhough they showed overfitting. Hence, the limited extent of the Epanechikov kernel has a strong effect 
on the final spatial resolution. 


```python
X0,y0 = np.random.rand(nc[0], 2), np.zeros((nc[0],), dtype=int)
X1,y1 = np.random.rand(nc[1], 2), np.ones((nc[1],), dtype=int)
X = np.r_[X0, X1]
y = np.r_[y0, y1]

# Normalize dataset
(X, Xgrid, xx, yy, x_min, x_max, y_min, y_max) = normalize_dataset(X, nbins)
```


```python
C_qda = qda.fit(X, y).predict(Xgrid) 
C_qda = C_qda.reshape(xx.shape)

h_gaussian = .07
h_epanechikov = .7

(C_g, C_e) = KDE_classifications(h_gaussian, h_epanechikov, X, y, Xgrid, nbins, prob_C)

# Compare classification results
figName = '2D_KDE_classif_random'
display_classification_results(C_g, C_e, C_qda, X, y, Xgrid, figName)
```


```python

```

6d18d1b2a415f6fa2cbf6a0d3ad60035e059f215
ipynb
Jupyter Notebook
2D_KDE_classification.ipynb
AstroPierre/Scripts-for-figures-courses-GIF-4101-GIF-7005
a38ad6f960cc6b8155fad00e4c4562f5e459f248
BSD-2-Clause
2D_KDE_classification.ipynb
AstroPierre/Scripts-for-figures-courses-GIF-4101-GIF-7005
a38ad6f960cc6b8155fad00e4c4562f5e459f248
BSD-2-Clause
2D_KDE_classification.ipynb
AstroPierre/Scripts-for-figures-courses-GIF-4101-GIF-7005
a38ad6f960cc6b8155fad00e4c4562f5e459f248
BSD-2-Clause
Qwen/Qwen-72B
1. YES
2. YES