---
# Example of classification using two gaussian distributions.
---

In this example, we generate a global probability distribution function (PDF) that is the sum of two gaussian PDF. 
Each one is weighted with its <i>a priori</i> class probability $P(C_{i})$:


<blockquote>  $P(\bf{x}) = P(\bf{x}|C_{1}) P(C_{1}) + P(\bf{x}|C_{2}) P(C_{2})$ </blockquote>



Our goal is to locate the influence zone of each class i. This corresponds to the zone where 

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


## The 2-D gausian PDF

In its matrix-form, the equation for the 2-D gausian PDF reads like this: 

<blockquote>  $P(\bf{x}) = \frac{1}{2\pi |\Sigma|^{0.5}} \exp{[-\frac{1}{2}(\bf{x}-\bf{\mu})^\top \Sigma^{-1} (\bf{x}-\bf{\mu})]}$  </blockquote>

where 

<blockquote> 
$
\begin{align}
\bf{x} &= [x_{1} x_{2}]^\top \\
\bf{\mu} &= [\mu_{1} \mu_{2}]^\top  \\
\Sigma &= \begin{pmatrix} \sigma_{x_{1}}^2 & \rho\sigma_{x_{1}}\sigma_{x_{2}} \\ \rho\sigma_{x_{1}}\sigma_{x_{2}} & \sigma_{x_{2}}^2 \end{pmatrix}
\end{align}
$                
</blockquote>

where $\rho$ is the correlation factor between the $x_{1}$ and $x_{2}$ data. 

## A useful trick to generate a covariance matrix $\Sigma$ with desired features 

Instead of guessing the values of $\sigma_{x_{1}}$, $\sigma_{x_{2}}$ and $\rho$ to create a 
covariance matrix $\Sigma$ for visualization purpose, we can design it with some desired characteristics. 
Those are the principal axis variances $\sigma_{1}^2$ and $\sigma_{2}^2$, and the rotation angle $\theta$.

First, we generate the covariance matrix of a correlated PDF with its principal axes oriented along the 
$x_{1}$ and $x_{2}$ axes.

<blockquote>  $\Sigma_{PA} = \begin{pmatrix} \sigma_{1}^2 & 0 \\ 0 & \sigma_{2}^2 \end{pmatrix}$ </blockquote>

Next, we generate the rotation matrix for the angle $\theta$:

<blockquote>  $R = \begin{pmatrix} \cos{\theta}  & -\sin{\theta} \\ \sin{\theta} & \cos{\theta} \end{pmatrix}$ </blockquote>

The covariance matrix we are looking for is 

<blockquote>  $\Sigma = R \Sigma_{PA} R^\top $ </blockquote>

Hence, we only have to specify the values of $\sigma_{1}$ and $\sigma_{2}$ along principal axes and the 
rotation angle $\theta$. The correlation coefficient $\rho$ depends on the value of $\theta$:

<blockquote>  $\rho>0$ when $\theta>0$ </blockquote>
<blockquote>  $\rho<0$ when $\theta<0$ </blockquote>

<br>
N.B. For ease of reading, we use below the variables x and y instead of x1 and x2. The final results are shown with x1 and x2.


```python
print(__doc__)

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

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import colors
from matplotlib import cm
from matplotlib.ticker import NullFormatter
from mpl_toolkits.mplot3d import Axes3D

from scipy import stats
from scipy.stats import multivariate_normal

from skimage.transform import resize

import seaborn as sns
sns.set(color_codes=True)

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

    Automatically created module for IPython interactive environment


Here are the characteristics of the 2-D PDF we want to generate:


```python
# Origins of the PDF
Mu = np.zeros((2,2))
Mu[0,:] = [6., 3.]
Mu[1,:] = [6., 7.]

# Standard deviations along the x1 and x2 directions (or x and y)
sigma = np.zeros((2,2))
sigma[0,:] = [2., 0.7]
sigma[1,:] = [1.2, 0.7]

# Rotation angles for the principal axis
theta = np.array([-30., 0.]) 

# A priori class probabilities
prob_C = np.array([0.5, 0.5]) 

```

Compute the covariance matrix $\Sigma$


```python
Sigma = np.zeros((2,2,2))
for i in range(2):
    angle = np.radians(theta[i])
    c, s = np.cos(angle), np.sin(angle)

    # Rotation matrix
    R = np.array(((c, -s), (s, c)))        

    # Covariance matrix for a PDF with its principal axes oriented along the x and y directions
    sig = np.array([[sigma[i,0]**2, 0.],[0., sigma[i,1]**2]])

    # Covariance matrix after rotation
    Sigma[:,:,i] = R.dot( sig.dot(R.T) )  

```

## Global PDF computation
Generate a spatial grid where the PDF will be evaluated locally.


```python
x_min, x_max = 0., 10.
y_min, y_max = 0., 10.
nx, ny = 300, 300
x = np.linspace(x_min, x_max, nx)
y = np.linspace(y_min, y_max, ny)
xx, yy = np.meshgrid(x,y) 
pos = np.dstack((xx, yy)) 
        
```

Compute the global PDF


```python
# Weighted PDF
model = multivariate_normal(Mu[0,:], Sigma[:,:,0]) 
pdf0 = prob_C[0]*model.pdf(pos)

model = multivariate_normal(Mu[1,:], Sigma[:,:,1]) 
pdf1 = prob_C[1]*model.pdf(pos)

# Global PDF 
pdf = pdf0 + pdf1

# Define a mask to identify the influence zone of each class: Mask=False for class 0 and Mask=True for class 1
mask = pdf1 > pdf0

```

## Display the global PDF 


```python
z_offset = -0.1
view=[22., -30.]


fig = plt.figure(figsize = (15,10))
ax = fig.gca(projection='3d')


# ---------Generate overlying 3-D surface with colors indicating the influence zone of each class
rstride, cstride = 5, 5

s = ax.plot_surface(xx, yy, pdf, rstride=rstride, cstride=cstride, linewidth=.5, antialiased=True, color='gray', edgecolors='k')       
a1 = s.__dict__['_original_facecolor']
b1 = s.__dict__['_facecolors']
c1 = s.__dict__['_facecolors3d']

s = ax.plot_surface(xx, yy, pdf, rstride=rstride, cstride=cstride, linewidth=.5, antialiased=True, color='w', edgecolors='k')
a2 = s.__dict__['_original_facecolor']
b2 = s.__dict__['_facecolors']
c2 = s.__dict__['_facecolors3d']

Lx = int(nx/rstride)
Ly = int(ny/cstride)

mask2 = resize(mask, (Lx,Ly), order=0)
indx = np.argwhere(mask2)
idx = indx[:,0]*Lx + indx[:,1]

a = a1
b = b1
c = c1
for i in idx:
    a[i,:] = a2[i,:]
    b[i,:] = b2[i,:]
    c[i,:] = c2[i,:]
    
s.__dict__['_original_facecolor'] = a
s.__dict__['_facecolors'] = b
s.__dict__['_facecolors3d'] = c


# --------- Generate underlying contours with borders between each class
cset = ax.contourf(xx, yy, pdf, zdir='z', offset=z_offset, cmap='viridis')
ax.contour(xx, yy, mask, [0.5], offset=z_offset, linewidths=2., colors='white') 

# N.B. For convenience, we use classes 1 and 2 instead of classes 0 and 1 in the figure
ax.text(2, 1, z_offset,'$\mathcal{R}_{1}$', zdir='y', color='w',fontsize=18)        
ax.text(9, 8, z_offset,'$\mathcal{R}_{2}$', zdir='y', color='w',fontsize=18)    
ax.text(6, 2.5, 0.05, '$p(\mathbf{X}|C_{1})p(C_{1})$', zdir='y', horizontalalignment='right', verticalalignment='bottom', 
        color='k',fontsize=18)    
ax.text(6, 6., 0.08, '$p(\mathbf{X}|C_{2})p(C_{2})$', zdir='y', horizontalalignment='right', verticalalignment='bottom', 
        color='gray',fontsize=18)    

ax.set_zlim(z_offset, np.max(pdf))
ax.view_init(view[0], view[1])

ax.set_ylabel('$x_{2}$', fontsize=18)
ax.xaxis.set_rotate_label(False)  
ax.set_xlabel('$x_{1}$', rotation=10, fontsize=18)
ax.set_zlabel('$Z$', rotation=0, fontsize=18)
ax.set_title('$Z = PDF(x_{1},x_{2})$', fontsize=20)
    
fig.tight_layout()

plt.savefig('Influence_zones_in_gaussian_PDF_classification.png')
plt.savefig('Influence_zones_in_gaussian_PDF_classification.pdf')

plt.show()

```


```python

```

75925852beae19ced419ba1288cb01ea2df63c20
ipynb
Jupyter Notebook
generate_explanatory_diagram_for_influence_zones_in_two_class_segmentation.ipynb
AstroPierre/Scripts-for-figures-courses-GIF-4101-GIF-7005
a38ad6f960cc6b8155fad00e4c4562f5e459f248
BSD-2-Clause
generate_explanatory_diagram_for_influence_zones_in_two_class_segmentation.ipynb
AstroPierre/Scripts-for-figures-courses-GIF-4101-GIF-7005
a38ad6f960cc6b8155fad00e4c4562f5e459f248
BSD-2-Clause
generate_explanatory_diagram_for_influence_zones_in_two_class_segmentation.ipynb
AstroPierre/Scripts-for-figures-courses-GIF-4101-GIF-7005
a38ad6f960cc6b8155fad00e4c4562f5e459f248
BSD-2-Clause
Qwen/Qwen-72B
1. YES
2. YES