# Exercise 4 - Clustering Algorithms

## a) Clustering using K-means

### Program 1 - Implementing K-means Clustering from scratch and using scikit learn

#### AIM:
To implement K-means clustering algorithm from scratch and also implement it with scikit learn.

#### ALGORITHM:
##### K-means:
1. Initialize randomly selected points as centroids.
2. Find the clusters for each point.
3. Compute new centroids by computing the mean of each cluster.
4. Repeat 1-3 until convergence.

##### Distrotion(metric):
Distrotion is calculated as the average of the squared distances from the cluster centers of the respective clusters. which can be mathematically defined as:
$$
\newcommand{\vect}{\mathbf}
\text{Distrotion} = \sum_{\vect{c}_k,X_k \in C} \sum_{\vect{x}_i \in X_k} (\vect{x}_i-\vect{c}_k)^2
$$
where,   
$C$ - Set of clusters  
$\vect{c}_k$ - centroid of the $k$th cluster  
$X_i$ - Points in the $k$th cluster

##### Sillhoutte Score(metric):
Silhoette Score is defined as the mean of the sillhouette coefficents of the follows:
Let $A$ be the cluster where the $i$th points is assigned and $C$ be any cluster.
$$
\text{SilhouetteScore} = \frac{1}{m}\sum_{i=1}^{m} s(i)\\
\text{silhouetteCoeffient. } s(i) = \frac{ b(i) - a(i)}{\max(a(i),b(i))}\\
b(i) = \min_{C\neq A} d(i,C)
$$
where,  
$a(i)$ - average distance from $i$th point to all points in $A$  
$d(i,C)$ - average distaance from $i$th point to all points in cluster $C$

#### Part 1 - Defining class for KMeans


```python
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from scipy.spatial.distance import cdist
```


```python
class KMeans:
    def __init__(self, n_clusters):
        self.n_clusters = n_clusters
        
    def fit(self, train_x, random_state = 0):
        m = train_x.shape[0]
        np.random.seed(random_state)
        self.centroids = train_x[np.random.choice(m,self.n_clusters,replace=False)]
        while True:
            old_centroids = self.centroids.copy()
            self.labels = self.predict(train_x)
            for cluster_id in range(self.n_clusters):
                cluster = train_x[self.labels == cluster_id]
                if len(cluster):
                    self.centroids[cluster_id] = cluster.mean(axis=0)
            if np.all(old_centroids != self.centroids) : break
                
    def predict(self, points):
        return np.argmin(cdist(points,self.centroids),axis=-1)
```

#### Part 2 - Defining evaluation metrics


```python
def distortion(X,labels,centroids):
    _, label_codes = np.unique(labels,return_inverse=True)
    return np.linalg.norm((X-centroids[label_codes]),axis=-1).mean()

def silhouette_score(X, labels):
    m = X.shape[0]
    clusters = np.unique(labels)
    n_clusters = clusters.shape[0]
    if n_clusters == 1: return -1
    cluster_codes = np.arange(n_clusters)
    dist_mask = labels.reshape(-1,1) == \
        ((labels + cluster_codes.reshape(-1,1))%n_clusters)[:,np.newaxis,:]
    dist_mask[0][np.eye(m).astype(bool)] = 0
    dist_mat = cdist(X,X)
    mean_dist = (dist_mask*dist_mat).sum(axis=-1)/dist_mask.sum(axis=-1)
    A = mean_dist[0]
    B = np.min(mean_dist[1:],axis=0)
    return ((B-A)/np.max(np.dstack([A,B]),axis=-1)).mean()
```

#### Part 3 - Loading the iris dataset


```python
names = ['s_len', 's_wid', 'p_len', 'p_wid', 'species']
iris_df = pd.read_csv('./datasets/iris.csv', header=None, names=names)
display(iris_df.head())
print(
    f"The data set contains {iris_df.shape[0]} records "
    f"and {iris_df.shape[1]} features.",
    iris_df['species'].value_counts(), sep="\n"
)

```


<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>s_len</th>
      <th>s_wid</th>
      <th>p_len</th>
      <th>p_wid</th>
      <th>species</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>5.1</td>
      <td>3.5</td>
      <td>1.4</td>
      <td>0.2</td>
      <td>setosa</td>
    </tr>
    <tr>
      <th>1</th>
      <td>4.9</td>
      <td>3.0</td>
      <td>1.4</td>
      <td>0.2</td>
      <td>setosa</td>
    </tr>
    <tr>
      <th>2</th>
      <td>4.7</td>
      <td>3.2</td>
      <td>1.3</td>
      <td>0.2</td>
      <td>setosa</td>
    </tr>
    <tr>
      <th>3</th>
      <td>4.6</td>
      <td>3.1</td>
      <td>1.5</td>
      <td>0.2</td>
      <td>setosa</td>
    </tr>
    <tr>
      <th>4</th>
      <td>5.0</td>
      <td>3.6</td>
      <td>1.4</td>
      <td>0.2</td>
      <td>setosa</td>
    </tr>
  </tbody>
</table>
</div>


    The data set contains 150 records and 5 features.
    setosa        50
    versicolor    50
    virginica     50
    Name: species, dtype: int64


#### Part 4 - Implementing K-means Clustering


```python
train_x = iris_df.drop("species", axis=1).values
k = 3
kmeans = KMeans(n_clusters=k)
kmeans.fit(train_x)
all_df = iris_df.copy(deep = True)
centroids = pd.DataFrame(kmeans.centroids, columns=names[:-1])
centroids["cluster"] = "centroid"
all_df['cluster'] = kmeans.labels.astype("str")
all_df = pd.concat([all_df, centroids])
```

#### Part 5 - Comparing the Clustering results with original classes


```python
plt.figure(figsize=(15,6))
plt.suptitle("classes vs clusters")
plt.subplot(121)
plt.title("Original Classes")
sns.scatterplot(data=all_df,x="s_len",y="s_wid",hue="species",s=80)
plt.legend(loc="upper right")

plt.subplot(122)
plt.title("Predicted Clusters and Centroids")
sns.scatterplot(
    data=all_df,x="s_len",y="s_wid",
    hue="cluster",style = "cluster",
    markers="osPD",s=80
)
plt.legend(loc="upper right")

plt.show()
```

#### Part 6 - Testing with multiple values of k(Hyperparameter tuning)


```python
sil_coefs = []
distortions = []
K = np.arange(2,6)
for k in K:
    kmeans = KMeans(n_clusters=k)
    kmeans.fit(train_x)
    avg_sil_coef = silhouette_score(train_x,kmeans.labels)
    dist = distortion(train_x,kmeans.labels,kmeans.centroids)
    print(f"For k={k:<4} Avg.Sil.Coef: {avg_sil_coef:<10.5f} Distortion: {dist:.5f}")
    distortions.append(dist)
    sil_coefs.append(avg_sil_coef)
```

    For k=2    Avg.Sil.Coef: 0.41657    Distortion: 1.30604
    For k=3    Avg.Sil.Coef: 0.49127    Distortion: 0.67576
    For k=4    Avg.Sil.Coef: 0.45307    Distortion: 0.59402
    For k=5    Avg.Sil.Coef: 0.32360    Distortion: 0.56721


#### Part 7 - Using Elbow method to Find the optimal value of k


```python
plt.plot(K, distortions, 's-', markersize=8)
plt.xlabel('k')
plt.xticks(K)
plt.ylabel('Distortion')
plt.title('Elbow Method for Finding Optimal k')
plt.show()
```

#### Part 8 -  Implementing K-means clustering with scikit learn


```python
from sklearn.cluster import KMeans as sklKMeans
from sklearn.metrics import silhouette_score as sk_silhouette_score
```


```python
k = 3
kmeans = sklKMeans(n_clusters=k, init='random').fit(train_x)
ss = sk_silhouette_score(train_x, kmeans.labels_)
dist = distortion(train_x,kmeans.labels_, kmeans.cluster_centers_)
print(f"For k={k}\tAvg. Sil. Coef: {ss:.5f}\tDistortion: {dist:.5f}")
```

    For k=3	Avg. Sil. Coef: 0.55259	Distortion: 0.64884


## b) Clustering using Gaussian Mixture Model(GMM)

### Program 1 - Implementing GMM from scratch using Expectation and Maximization algorithm 

#### AIM:

To implement GMM from scratch in python using Expectation and Maximization algorithm.

#### Formula:

##### Log-likelihood
$$
l(\vect{\theta}) = \sum_i\log\left( \sum_k p(\vect{x}_i,Z=k|\theta)\right) = \sum_i\log \left(\sum_k \pi_k \mathcal{N}(\vect{x}_i|\vect{\mu}_k.\vect{\Sigma}_k))\right)
$$

##### Auxilary Function
$$
\newcommand{\vect}[1]{\boldsymbol{\mathbf{#1}}}
$$
$$
Q(\vect{\theta},\vect{\theta}^{(t-1)})= \sum_i\sum_k r_{ik} \log \pi_k +\sum_i\sum_k r_{ik} \log p(\vect{x}_i |\vect{\theta}_k )
$$
where $r_{ik} \triangleq p(z_i = k|x_i , \vect{\theta}^{(t−1)} )$ is the responsibility that cluster $k$ takes for data point $i$.

##### Responibility
$$
r_{ik}= \frac{\pi_k p(\vect{x}_i|\vect{\theta}_k^{(t-1)})}
{\sum_{k'}\pi_{k'} p(\vect{x}_i|\vect{\theta}_{k'}^{(t-1)})}
=\frac{\pi_k \mathcal{N}(\vect{x}_i|\vect{\mu}_k,\vect{\Sigma}_k)}
{\sum_{k'}\pi_{k'} \mathcal{N}(\vect{x}_i|\vect{\mu}_{k'},\vect{\Sigma}_{k'}))}
$$

##### Parameter updates
$$
\begin{align}
\pi_k &= \frac{1}{N}\sum_i r_{ik} = \frac{r_k}{N}\\
\vect{\mu}_k &= \frac{\sum_i r_{ik} \vect{x}_i}{r_k}\\
\vect{\Sigma}_k & = \frac{\sum_i r_{ik}(\vect{x}_i-\vect{\mu}_k)(\vect{x}_i-\vect{\mu}_k)^T}{r_k}
\end{align}
$$
where   
$r_k \triangleq \sum_i r_{ik}$ - the weighted number of points.  
$N$ - number of Points

#### ALGORITHM:
##### GMM Using EM:
1. Initialize the parmeters of the gaussian mixture distribution $\theta^0 = (\vect{\theta_k},\vect{\pi})$
2. Until convergence, do the following:
    1. E-step - Compute the terms inside the auxilary function $Q(\theta,\theta^{(t-1)})$, i.e, the responsibilities $r_{ik}$
    2. M-step - Compute the updated parameters $(\pi_k,\vect{\mu}_k,\vect{\Sigma}_k)$ for each cluster $k$ of the Gaussian Mixture distribution

#### Part 1 - Defining the class for GMM


```python
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
from sklearn.cluster import KMeans as SklKMeans
```


```python
class GMM:
    def __init__(self, 
        n_components: int, n_iters: int, tol: float,
        random_state: int=0, init_params="random"
    ):
        self.n_components = n_components
        self.n_iters = n_iters
        self.tol = tol
        self.init_params = init_params
        self.random_state = random_state


    def fit(self, X,plot=False,plot_params={}):
        m,n= X.shape
        self.X = X
        k = self.n_components
        
        self.resp = np.zeros((m, k))
        self.weights = np.full(k, 1 / k)
        
        
        if self.init_params == "kmeans":
            kmeans = SklKMeans(k,random_state=self.random_state)
            kmeans.fit(X)
            self.means = kmeans.cluster_centers_
        if self.init_params == "random" :
            np.random.seed(self.random_state)
            self.means = X[np.random.choice(m, k , replace=False)]
            

        self.covs = np.full((k, n, n), np.cov(X, rowvar=False))
    
        self.converged = False
        self.log_likelihood_trace = []
        
        if plot:
            fig,ax = plt.subplots(1,4,figsize=(20,5))
            fig.suptitle("GMM Using Estimation and Maximization")
            self.draw(ax[0],"Initial Clusters", **plot_params)
        
        for i in range(self.n_iters):
            self._do_estep(X)
            self._do_mstep(X)
            self.log_likelihood_trace.append(self.log_likelihood)
            if i == 0:
                if plot:
                    self.draw(ax[1],"Clusters after 1 iteration", **plot_params)
            else:
                ol,nl = self.log_likelihood_trace[-2:]
                if nl-ol <= self.tol:  #likelihood difference less than tol
                    self.converged = True
                    break
        if plot:
            self.draw(ax[2],f"clusters after {i+1} iterations", **plot_params)
            ax[3].plot(self.log_likelihood_trace, "-o")
            ax[3].set_title("Log Likelihood")
            plt.show()

    def _do_estep(self, X):
        for k in range(self.n_components):
            self.resp[:, k] = self.weights[k] * multivariate_normal(self.means[k], self.covs[k]).pdf(X)
        p_x = np.sum(self.resp, axis=1,keepdims=1)
        self.resp = self.resp / p_x
        self.log_likelihood = np.sum(np.log(p_x))
        
    def _do_mstep(self, X):
        resp_weights = self.resp.sum(axis=0)
        self.weights = resp_weights / X.shape[0]
        self.means = self.resp.T @ X/ resp_weights.reshape(-1, 1)
        diff = X[:,np.newaxis,:] - self.means
        self.covs = np.einsum('ik,ikj,ikl->kjl',self.resp,diff,diff)/resp_weights.reshape(-1,1,1)
        

    def draw(self,ax,  title="", **plot_params):
        ax.set_title(title)
        ax.scatter(self.X[:, 0], self.X[:, 1],**plot_params)
        
        delta = 0.05
        k = self.n_components
        x = np.arange(*ax.get_xlim(), delta)
        y = np.arange(*ax.get_ylim(), delta)
        x, y = np.meshgrid(x, y)
        col = [f"C{(i+1)%10}" for i in range(k)]
        for i in range(k):
            mean, cov = self.means[i], self.covs[i]
            z = multivariate_normal(mean, cov).pdf(np.dstack([x, y]))
            ax.scatter(mean[0], mean[1], color=col[i])
            ax.contour(x, y, z, levels=[.01], colors=col[i])

```

#### Part 2 - Generating Random data from a Gaussian Mixture


```python

def gen_data(k=3, dim=2, points_per_cluster=200, lim=[-10, 10],seed = 1):
    np.random.seed(seed)
    x = np.ndarray((k,points_per_cluster,dim))
    mean = np.random.rand(k, dim)*(lim[1]-lim[0]) + lim[0]
    for i in range(k):
        cov = np.random.rand(dim, dim+10)
        cov = cov @ cov.T
        x[i] = np.random.multivariate_normal(mean[i], cov, points_per_cluster)
    x = x.reshape(-1,dim)

    if(dim == 2):
        plt.figure()
        plt.title("generated random gaussian distributions")
        plt.scatter(x[:, 0], x[:, 1], s=1, alpha=0.4)
        plt.show()
        
    return x

```


```python
X = gen_data(k=3, dim=2, points_per_cluster=1000, seed=3)
```

#### Part 3 - Implementing GMM on generated random data


```python
gmm = GMM(n_components =3,n_iters=30,tol=10e-4, random_state=5, init_params="random")
gmm.fit(X, plot= True,plot_params={"s":1,"alpha":.2})
```

#### Part 4 - Implementing GMM by using kmeans to initialize means


```python
k = 3
gmm = GMM(n_components = k,n_iters = 30,tol=10e-4,random_state=1, init_params="kmeans")
gmm.fit(X,plot=True,plot_params={"s":1,"alpha":.2})

```

efe9f36feda9d8ef4b62192626b5124e137d918e
ipynb
Jupyter Notebook
notebooks/Exercise 4.ipynb
msuriya2001/machine-learning-lab-exercises
1f638ba437a918aafec95bf26087c39059c9b81d
MIT
2022-02-03T03:23:10.000Z
2022-02-03T03:23:10.000Z
notebooks/Exercise 4.ipynb
msuriya2001/machine-learning-lab-exercises
1f638ba437a918aafec95bf26087c39059c9b81d
MIT
notebooks/Exercise 4.ipynb
msuriya2001/machine-learning-lab-exercises
1f638ba437a918aafec95bf26087c39059c9b81d
MIT
2021-08-24T17:02:11.000Z
2022-03-09T13:19:09.000Z
Qwen/Qwen-72B
1. YES
2. YES