# Session 10 - Unsupervised Learning 

## Contents

- [Principal Components Analysis](#Principal-Components-Analysis)
- [Clustering Methods](#Clustering-Methods)


```python
# Import
import pandas as pd
import numpy as np
import seaborn as sns
import time

from sklearn.preprocessing import scale
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans
from scipy.cluster import hierarchy
from scipy.cluster.hierarchy import linkage, dendrogram, cut_tree
```


```python
# Import matplotlib for graphs
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
from IPython.display import clear_output

# Set global parameters
%matplotlib inline
plt.style.use('seaborn-white')
plt.rcParams['lines.linewidth'] = 3
plt.rcParams['figure.figsize'] = (10,6)
plt.rcParams['figure.titlesize'] = 20
plt.rcParams['axes.titlesize'] = 18
plt.rcParams['axes.labelsize'] = 14
plt.rcParams['legend.fontsize'] = 14
```

The difference between *supervised learning* and *unsupervised learning* is that in the first case we have a variable $y$ which we want to predict, given a set of variables $\{ X_1, . . . , X_p \}$.

In *unsupervised learning* are not interested in prediction, because we do not have an associated response variable $y$. Rather, the goal is to discover interesting things about the measurements on $\{ X_1, . . . , X_p \}$. Question that we can answer are:

- Is there an informative way to visualize the data? 
- Can we discover subgroups among the variables or among the observations?

We are going to answe these questions with two main methods, respectively:

- Principal Components Analysis
- Clustering

## Principal Component Analysis

Suppose that we wish to visualize $n$ observations with measurements on a set of $p$ features, $\{X_1, . . . , X_p\}$, as part of an exploratory data analysis.

We could do this by examining two-dimensional scatterplots of the data, each of which contains the n observations’ measurements on two of the features. However, there are $p(p−1)/2$ such scatterplots; for example,
with $p = 10$ there are $45$ plots! 

PCA provides a tool to do just this. It finds a low-dimensional represen- tation of a data set that contains as much as possible of the variation. 

The **first principal component** of a set of features $\{X_1, . . . , X_p\}$ is the normalized linear combination of the features

$$
Z_1 = \phi_{11} X_1 + \phi_{21} X_2 + ... + \phi_{p1} X_p
$$

that has the largest variance. By normalized, we mean that $\sum_{i=1}^p \phi^2_{i1} = 1$.

In other words, the first principal component loading vector solves the optimization problem

$$
\underset{\phi_{11}, \ldots, \phi_{p 1}}{\operatorname{max}} \ \left\{\frac{1}{n} \sum_{i=1}^{n}\left(\sum_{j=1}^{p} \phi_{j 1} x_{i j}\right)^{2}\right\} \quad \text { subject to } \quad \sum_{j=1}^{p} \phi_{j 1}^{2}=1
$$

The objective that we are maximizing is just the sample variance of the $n$ values of $z_{i1}$.

After the first principal component $Z_1$ of the features has been determined, we can find the second principal component $Z_2$. The **second principal component** is the linear combination of $\{X_1, . . . , X_p\}$ that has maximal variance out of all linear combinations that are *uncorrelated* with $Z_1$.

We illustrate the use of PCA on the `USArrests` data set. For each of the 50 states in the United States, the data set contains the number of arrests per $100,000$ residents for each of three crimes: `Assault`, `Murder`, and `Rape.` We also record the percent of the population in each state living in urban areas, `UrbanPop`.


```python
# Load crime data
df = pd.read_csv('data/USArrests.csv', index_col=0)
df.head()
```




<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>Murder</th>
      <th>Assault</th>
      <th>UrbanPop</th>
      <th>Rape</th>
    </tr>
    <tr>
      <th>State</th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>Alabama</th>
      <td>13.2</td>
      <td>236</td>
      <td>58</td>
      <td>21.2</td>
    </tr>
    <tr>
      <th>Alaska</th>
      <td>10.0</td>
      <td>263</td>
      <td>48</td>
      <td>44.5</td>
    </tr>
    <tr>
      <th>Arizona</th>
      <td>8.1</td>
      <td>294</td>
      <td>80</td>
      <td>31.0</td>
    </tr>
    <tr>
      <th>Arkansas</th>
      <td>8.8</td>
      <td>190</td>
      <td>50</td>
      <td>19.5</td>
    </tr>
    <tr>
      <th>California</th>
      <td>9.0</td>
      <td>276</td>
      <td>91</td>
      <td>40.6</td>
    </tr>
  </tbody>
</table>
</div>




```python
# Scale data
X_scaled = pd.DataFrame(scale(df), index=df.index, columns=df.columns).values
```


```python
# Fit PCA with 2 components
pca2 = PCA(n_components=2).fit(X_scaled)
```


```python
# Get weights
weights = pca2.components_.T
df_weights = pd.DataFrame(weights, index=df.columns, columns=['PC1', 'PC2'])
df_weights
```




<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>PC1</th>
      <th>PC2</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>Murder</th>
      <td>0.535899</td>
      <td>0.418181</td>
    </tr>
    <tr>
      <th>Assault</th>
      <td>0.583184</td>
      <td>0.187986</td>
    </tr>
    <tr>
      <th>UrbanPop</th>
      <td>0.278191</td>
      <td>-0.872806</td>
    </tr>
    <tr>
      <th>Rape</th>
      <td>0.543432</td>
      <td>-0.167319</td>
    </tr>
  </tbody>
</table>
</div>




```python
# Transform X to get the principal components
X_dim2 = pca2.transform(X_scaled)
df_dim2 = pd.DataFrame(X_dim2, columns=['PC1', 'PC2'], index=df.index)
df_dim2.head()
```




<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>PC1</th>
      <th>PC2</th>
    </tr>
    <tr>
      <th>State</th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>Alabama</th>
      <td>0.985566</td>
      <td>1.133392</td>
    </tr>
    <tr>
      <th>Alaska</th>
      <td>1.950138</td>
      <td>1.073213</td>
    </tr>
    <tr>
      <th>Arizona</th>
      <td>1.763164</td>
      <td>-0.745957</td>
    </tr>
    <tr>
      <th>Arkansas</th>
      <td>-0.141420</td>
      <td>1.119797</td>
    </tr>
    <tr>
      <th>California</th>
      <td>2.523980</td>
      <td>-1.542934</td>
    </tr>
  </tbody>
</table>
</div>




```python
fig, ax1 = plt.subplots(figsize=(10,10))
ax1.set_title('Figure 10.1');

# Plot Principal Components 1 and 2
for i in df_dim2.index:
    ax1.annotate(i, (df_dim2.PC1.loc[i], -df_dim2.PC2.loc[i]), ha='center', fontsize=14)
    
# Plot reference lines
m = np.max(np.abs(df_dim2.values))*1.2
ax1.hlines(0,-m,m, linestyles='dotted', colors='grey')
ax1.vlines(0,-m,m, linestyles='dotted', colors='grey')
ax1.set_xlabel('First Principal Component')
ax1.set_ylabel('Second Principal Component')
ax1.set_xlim(-m,m); ax1.set_ylim(-m,m)

# Plot Principal Component loading vectors, using a second y-axis.
ax1b = ax1.twinx().twiny() 
ax1b.set_ylim(-1,1); ax1b.set_xlim(-1,1)
for i in df_weights[['PC1', 'PC2']].index:
    ax1b.annotate(i, (df_weights.PC1.loc[i]*1.05, -df_weights.PC2.loc[i]*1.05), color='orange', fontsize=16)
    ax1b.arrow(0,0,df_weights.PC1[i], -df_weights.PC2[i], color='orange', lw=2)
```

### PCA and spectral analysis

In case you haven't noticed, calculating principal components, is equivalent to calculating the eigenvectors of the design matrix $X'X$, i.e. the variance-covariance matrix of $X$. Indeed what we performed above is a decomposition of the variance of $X$ into orthogonal components.

The constrained maximization problem above can be re-written in matrix notation as

$$
\max \ \phi' X'X \phi \quad \text{ s. t. } \quad \phi'\phi = 1
$$

Which has the following dual representation

$$
\mathcal L (\phi, \lambda) = \phi' X'X \phi - \lambda (\phi'\phi - 1)
$$

Which gives first order conditions

$$
\begin{align}
& \frac{\partial \mathcal L}{\partial \lambda} = \phi'\phi - 1 \\
& \frac{\partial \mathcal L}{\partial \phi} = 2 X'X \phi - 2 \lambda \phi
\end{align}
$$

Setting the derivatives to zero at the optimum, we get

$$
\begin{align}
& \phi'\phi = 1 \\
& X'X \phi = \lambda \phi
\end{align}
$$

Thus, $\phi$ is an **eigenvector** of the covariance matrix $X'X$, and the maximizing vector will be the one associated with the largest **eigenvalue** $\lambda$. 

We can now double-check it using `numpy` linear algebra package.


```python
eigenval, eigenvec = np.linalg.eig(X_scaled.T @ X_scaled)
data = np.concatenate((eigenvec,eigenval.reshape(1,-1)))
idx = list(df.columns) + ['Eigenvalue']
df_eigen = pd.DataFrame(data, index=idx, columns=['PC1', 'PC2','PC3','PC4'])

df_eigen
```




<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>PC1</th>
      <th>PC2</th>
      <th>PC3</th>
      <th>PC4</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>Murder</th>
      <td>0.535899</td>
      <td>0.418181</td>
      <td>0.649228</td>
      <td>-0.341233</td>
    </tr>
    <tr>
      <th>Assault</th>
      <td>0.583184</td>
      <td>0.187986</td>
      <td>-0.743407</td>
      <td>-0.268148</td>
    </tr>
    <tr>
      <th>UrbanPop</th>
      <td>0.278191</td>
      <td>-0.872806</td>
      <td>0.133878</td>
      <td>-0.378016</td>
    </tr>
    <tr>
      <th>Rape</th>
      <td>0.543432</td>
      <td>-0.167319</td>
      <td>0.089024</td>
      <td>0.817778</td>
    </tr>
    <tr>
      <th>Eigenvalue</th>
      <td>124.012079</td>
      <td>49.488258</td>
      <td>8.671504</td>
      <td>17.828159</td>
    </tr>
  </tbody>
</table>
</div>



The spectral decomposition of the variance of $X$ generates a set of orthogonal vectors (eigenvectors) with different magnitudes (eigenvalues). The eigenvalues tell us the amount of variance of the data in that direction.

If we combine the eigenvectors together, we form a projection matrix $P$ that we can use to transform the original variables: $\tilde X = P X$


```python
X_transformed = X_scaled @ eigenvec
df_transformed = pd.DataFrame(X_transformed, index=df.index, columns=['PC1', 'PC2','PC3','PC4'])

df_transformed.head()
```




<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>PC1</th>
      <th>PC2</th>
      <th>PC3</th>
      <th>PC4</th>
    </tr>
    <tr>
      <th>State</th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>Alabama</th>
      <td>0.985566</td>
      <td>1.133392</td>
      <td>0.156267</td>
      <td>-0.444269</td>
    </tr>
    <tr>
      <th>Alaska</th>
      <td>1.950138</td>
      <td>1.073213</td>
      <td>-0.438583</td>
      <td>2.040003</td>
    </tr>
    <tr>
      <th>Arizona</th>
      <td>1.763164</td>
      <td>-0.745957</td>
      <td>-0.834653</td>
      <td>0.054781</td>
    </tr>
    <tr>
      <th>Arkansas</th>
      <td>-0.141420</td>
      <td>1.119797</td>
      <td>-0.182811</td>
      <td>0.114574</td>
    </tr>
    <tr>
      <th>California</th>
      <td>2.523980</td>
      <td>-1.542934</td>
      <td>-0.341996</td>
      <td>0.598557</td>
    </tr>
  </tbody>
</table>
</div>



This is exactly the dataset that we obtained before.

### More on PCA

#### Scaling the Variables

The results obtained when we perform PCA will also depend on whether the variables have been individually scaled. In fact, the variance of a variable depends on its magnitude.


```python
# Variables variance
df.var(axis=0)
```




    Murder        18.970465
    Assault     6945.165714
    UrbanPop     209.518776
    Rape          87.729159
    dtype: float64



Consequently, if we perform PCA on the unscaled variables, then the first principal component loading vector will have a very large loading for `Assault`, since that variable has by far the highest variance.


```python
# Fit PCA with unscaled varaibles
X = df.values
pca2_u = PCA(n_components=2).fit(X)
```


```python
# Get weights
weights_u = pca2_u.components_.T
df_weights_u = pd.DataFrame(weights_u, index=df.columns, columns=['PC1', 'PC2'])
df_weights_u
```




<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>PC1</th>
      <th>PC2</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>Murder</th>
      <td>0.041704</td>
      <td>0.044822</td>
    </tr>
    <tr>
      <th>Assault</th>
      <td>0.995221</td>
      <td>0.058760</td>
    </tr>
    <tr>
      <th>UrbanPop</th>
      <td>0.046336</td>
      <td>-0.976857</td>
    </tr>
    <tr>
      <th>Rape</th>
      <td>0.075156</td>
      <td>-0.200718</td>
    </tr>
  </tbody>
</table>
</div>




```python
# Transform X to get the principal components
X_dim2_u = pca2_u.transform(X)
df_dim2_u = pd.DataFrame(X_dim2_u, columns=['PC1', 'PC2'], index=df.index)
df_dim2_u.head()
```




<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>PC1</th>
      <th>PC2</th>
    </tr>
    <tr>
      <th>State</th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>Alabama</th>
      <td>64.802164</td>
      <td>11.448007</td>
    </tr>
    <tr>
      <th>Alaska</th>
      <td>92.827450</td>
      <td>17.982943</td>
    </tr>
    <tr>
      <th>Arizona</th>
      <td>124.068216</td>
      <td>-8.830403</td>
    </tr>
    <tr>
      <th>Arkansas</th>
      <td>18.340035</td>
      <td>16.703911</td>
    </tr>
    <tr>
      <th>California</th>
      <td>107.422953</td>
      <td>-22.520070</td>
    </tr>
  </tbody>
</table>
</div>




```python
fig, (ax1,ax2) = plt.subplots(1,2,figsize=(18,9))

# Scaled PCA
for i in df_dim2.index:
    ax1.annotate(i, (df_dim2.PC1.loc[i], -df_dim2.PC2.loc[i]), ha='center', fontsize=14)
ax1b = ax1.twinx().twiny() 
ax1b.set_ylim(-1,1); ax1b.set_xlim(-1,1)
for i in df_weights[['PC1', 'PC2']].index:
    ax1b.annotate(i, (df_weights.PC1.loc[i]*1.05, -df_weights.PC2.loc[i]*1.05), color='orange', fontsize=16)
    ax1b.arrow(0,0,df_weights.PC1[i], -df_weights.PC2[i], color='orange', lw=2)
ax1.set_title('Scaled')
    
# Unscaled PCA
for i in df_dim2_u.index:
    ax2.annotate(i, (df_dim2_u.PC1.loc[i], -df_dim2_u.PC2.loc[i]), ha='center', fontsize=14)
ax2b = ax2.twinx().twiny() 
ax2b.set_ylim(-1,1); ax2b.set_xlim(-1,1)
for i in df_weights_u[['PC1', 'PC2']].index:
    ax2b.annotate(i, (df_weights_u.PC1.loc[i]*1.05, -df_weights_u.PC2.loc[i]*1.05), color='orange', fontsize=16)
    ax2b.arrow(0,0,df_weights_u.PC1[i], -df_weights_u.PC2[i], color='orange', lw=2)
ax2.set_title('Unscaled')

# Plot reference lines
for ax,df in zip((ax1,ax2), (df_dim2,df_dim2_u)):
    m = np.max(np.abs(df.values))*1.2
    ax.hlines(0,-m,m, linestyles='dotted', colors='grey')
    ax.vlines(0,-m,m, linestyles='dotted', colors='grey')
    ax.set_xlabel('First Principal Component')
    ax.set_ylabel('Second Principal Component')
    ax.set_xlim(-m,m); ax.set_ylim(-m,m)
```

As predicted, the first principal component loading vector places almost all of its weight on `Assault`, while the second principal component loading vector places almost all of its weight on `UrpanPop`. Comparing this to the left-hand plot, we see that scaling does indeed have a substantial effect on the results obtained. However, this result is simply a consequence of the scales on which the variables were measured. 

#### The Proportion of Variance Explained

We can now ask a natural question: how much of the information in a given data set is lost by projecting the observations onto the first few principal components? That is, how much of the variance in the data is not contained in the first few principal components? More generally, we are interested in knowing the **proportion of variance explained (PVE)** by each principal component. 


```python
# Four components
pca4 = PCA(n_components=4).fit(X_scaled)
```


```python
# Variance of the four principal components
pca4.explained_variance_
```




    array([2.53085875, 1.00996444, 0.36383998, 0.17696948])




```python
# As a percentage of the total variance
pca4.explained_variance_ratio_
```




    array([0.62006039, 0.24744129, 0.0891408 , 0.04335752])



In the `Arrest` dataset, the first principal component explains $62.0\%$ of the variance in the data, and the next principal component explains $24.7\%$ of the variance. Together, the first two principal components explain almost $87\%$ of the variance in the data, and the last two principal components explain only $13\%$ of the variance.


```python
fig, (ax1, ax2) = plt.subplots(1,2, figsize=(12,5))
fig.suptitle('Figure 10.2');

# Relative 
ax1.plot([1,2,3,4], pca4.explained_variance_ratio_)
ax1.set_ylabel('Prop. Variance Explained')
ax1.set_xlabel('Principal Component');

# Cumulative
ax2.plot([1,2,3,4], np.cumsum(pca4.explained_variance_ratio_))
ax2.set_ylabel('Cumulative Variance Explained');
ax2.set_xlabel('Principal Component');
```

#### Deciding How Many Principal Components to Use

In general, a $n \times p$ data matrix $X$ has $\min\{n − 1, p\}$ distinct principal components. However, we usually are not interested in all of them; rather, we would like to use just the first few principal components in order to visualize or interpret the data. 

We typically decide on the number of principal components required to visualize the data by examining a *scree plot*.

However, there is no well-accepted objective way to decide how many principal com- ponents are enough.

## Clustering

Clustering refers to a very broad set of techniques for finding subgroups, or clusters, in a data set. When we cluster the observations of a data set, we seek to partition them into distinct groups so that the observations within each group are quite similar to each other, while observations in different groups are quite different from each other.

In this section we focus on perhaps the two best-known clustering approaches: 

1. **K-means clustering**: we seek to partition the observations into a pre-specified clustering number of clusters
2. **Hierarchical clustering**: we do not know in advance how many clusters we want; in fact, we end up with a tree-like visual representation of the observations, called a dendrogram, that allows us to view at once the clusterings obtained for each possible number of clusters, from 1 to n.

### K-Means Clustering

The idea behind K-means clustering is that a good clustering is one for which the within-cluster variation is as small as possible. Hence we want to solve the problem

$$
\underset{C_{1}, \ldots, C_{K}}{\operatorname{minimize}}\left\{\sum_{k=1}^{K} W\left(C_{k}\right)\right\}
$$

where $C_k$ is a cluster and $ W(C_k)$ is a measure of the amount by which the observations within a cluster differ from each other.

There are many possible ways to define this concept, but by far the most common choice involves **squared Euclidean distance**. That is, we define

$$
W\left(C_{k}\right)=\frac{1}{\left|C_{k}\right|} \sum_{i, i^{\prime} \in C_{k}} \sum_{j=1}^{p}\left(x_{i j}-x_{i^{\prime} j}\right)^2
$$

where $|C_k|$ denotes the number of observations in the $k^{th}$ cluster. 

### Algorithm

1. Randomly assign a number, from $1$ to $K$, to each of the observations. These serve as initial cluster assignments for the observations.

2. Iterate until the cluster assignments stop changing:

    a) For each of the $K$ clusters, compute the cluster centroid. The kth cluster centroid is the vector of the $p$ feature means for the observations in the $k^{th}$ cluster.
    
    b) Assign each observation to the cluster whose centroid is closest (where closest is defined using Euclidean distance).


```python
np.random.seed(123)

# Simulate data
X = np.random.randn(50,2)
X[0:25, 0] = X[0:25, 0] + 3
X[0:25, 1] = X[0:25, 1] - 4
```


```python
fig, ax = plt.subplots(figsize=(6, 5))
fig.suptitle("Baseline")

# Plot
ax.scatter(X[:,0], X[:,1], s=50, alpha=0.5, c='k') 
ax.set_xlabel('X0'); ax.set_ylabel('X1');
```


```python
np.random.seed(1)

# Init clusters
K = 2
clusters0 = np.random.randint(K,size=(np.size(X,0)))
```


```python
fig, ax = plt.subplots(figsize=(6, 5))
fig.suptitle("Random assignment")

# Plot
ax.scatter(X[clusters0==0,0], X[clusters0==0,1], s=50, alpha=0.5) 
ax.scatter(X[clusters0==1,0], X[clusters0==1,1], s=50, alpha=0.5)
ax.set_xlabel('X0'); ax.set_ylabel('X1');
```


```python
# Compute new centroids
def compute_new_centroids(X, clusters):
    K = len(np.unique(clusters))
    centroids = np.zeros((K,np.size(X,1)))
    for k in range(K):
        if sum(clusters==k)>0:
            centroids[k,:] = np.mean(X[clusters==k,:], axis=0)
        else:
            centroids[k,:] = np.mean(X, axis=0)
    return centroids
```


```python
# Print
centroids0 = compute_new_centroids(X, clusters0)
print(centroids0)
```

    [[ 1.35725989 -2.15281035]
     [ 1.84654757 -1.99437838]]



```python
# Plot assignment
def plot_assignment(X, centroids, clusters, d, i):
    clear_output(wait=True)
    fig, ax = plt.subplots(figsize=(6, 5))
    fig.suptitle("Iteration %.0f: inertia=%.1f" % (i,d))

    # Plot
    ax.clear()
    colors = plt.rcParams['axes.prop_cycle'].by_key()['color'];
    K = np.size(centroids,0)
    for k in range(K):
        ax.scatter(X[clusters==k,0], X[clusters==k,1], s=50, c=colors[k], alpha=0.5) 
        ax.scatter(centroids[k,0], centroids[k,1], marker = '*', s=300, color=colors[k])
        ax.set_xlabel('X0'); ax.set_ylabel('X1');
    
    # Show
    plt.show();
```


```python
# Plot
plot_assignment(X, centroids0, clusters0, 0, 0)
```


```python
# Assign X to clusters
def assign_to_cluster(X, centroids):
    K = np.size(centroids,0)
    dist = np.zeros((np.size(X,0),K))
    for k in range(K):
        dist[:,k] = np.mean((X - centroids[k,:])**2, axis=1)
    clusters = np.argmin(dist, axis=1)
    
    # Compute inertia
    inertia = 0
    for k in range(K):
        if sum(clusters==k)>0:
            inertia += np.sum((X[clusters==k,:] - centroids[k,:])**2)
    return clusters, inertia
```


```python
# Get cluster assignment
[clusters1,d] = assign_to_cluster(X, centroids0)
```


```python
# Plot
plot_assignment(X, centroids0, clusters1, d, 1)
```

We now have all the components to proceed iteratively.


```python
def kmeans_manual(X, K):

    # Init
    i = 0
    d0 = 1e4
    d1 = 1e5
    clusters = np.random.randint(K,size=(np.size(X,0)))

    # Iterate until convergence
    while np.abs(d0-d1) > 1e-10:
        d1 = d0
        centroids = compute_new_centroids(X, clusters)
        [clusters, d0] = assign_to_cluster(X, centroids)
        plot_assignment(X, centroids, clusters, d0, i)
        i+=1
```


```python
# Test
kmeans_manual(X, K)
```

Here the observations can be easily plotted because they are two-dimensional.
If there were more than two variables then we could instead perform PCA
and plot the first two principal components score vectors.

In this example, we knew that there really were two clusters because
we generated the data. However, for real data, in general we do not know
the true number of clusters. We could instead have performed K-means
clustering on this example with `K  =  3`. If we do this, K-means clustering will split up the two "real" clusters, since it has no information about them:


```python
# K=3
kmeans_manual(X, 3)
```

The automated function in `sklearn` to persorm $K$-means clustering is `KMeans`. 


```python
# SKlearn algorithm
km1 = KMeans(n_clusters=3, n_init=1, random_state=1)
km1.fit(X)
```




    KMeans(n_clusters=3, n_init=1, random_state=1)




```python
# Plot
plot_assignment(X, km1.cluster_centers_, km1.labels_, km1.inertia_, km1.n_iter_)
```

As we can see, the results are different in the two algorithms? Why? $K$-means is susceptible to the initial values. One way to solve this problem is to run the algorithm multiple times and report only the best results

To run the `Kmeans()` function in python with multiple initial cluster assignments, we use the `n_init` argument (default: 10). If a value of `n_init` greater than one is used, then K-means clustering will be performed using multiple random assignments, and the `Kmeans()` function will report only the best results.


```python
# 30 runs
km_30run = KMeans(n_clusters=3, n_init=30, random_state=1).fit(X)
plot_assignment(X, km_30run.cluster_centers_, km_30run.labels_, km_30run.inertia_, km_30run.n_iter_)
```

It is generally recommended to always run K-means clustering with a large value of `n_init`, such as 20 or 50 to avoid getting stuck in an undesirable local optimum.

When performing K-means clustering, in addition to using multiple initial cluster assignments, it is also important to set a random seed using the `random_state` parameter. This way, the initial cluster assignments can be replicated, and the K-means output will be fully reproducible.

### Hierarchical Clustering

One potential disadvantage of K-means clustering is that it requires us to pre-specify the number of clusters $K$. Hierarchical clustering is an alternative approach which does not require that we commit to a particular choice of $K$


#### The Dendogram

Hierarchical clustering has an added advantage over K-means clustering in that it results in an attractive tree-based representation of the observations, called a **dendrogram**.


```python
plt.figure(figsize=(25, 10))
plt.title('Hierarchical Clustering Dendrogram')

# calculate full dendrogram
plt.xlabel('sample index')
plt.ylabel('distance')
dendrogram(
    linkage(X, "complete"),
    leaf_rotation=90.,  # rotates the x axis labels
    leaf_font_size=8.,  # font size for the x axis labels
)
plt.show()
```

Each leaf of the *dendrogram* represents one observation. 

As we move up the tree, some leaves begin to fuse into branches. These correspond to observations that are similar to each other. As we move higher up the tree, branches themselves fuse, either with leaves or other branches. The earlier (lower in the tree) fusions occur, the more similar the groups of observations are to each other.

We can use de *dendogram* to understand how similar two observations are: we can look for the point in the tree where branches containing those two obse rvations are first fused. The height of this fusion, as measured on the vertical axis, indicates how different the two observations are. Thus, observations that fuse at the very bottom of the tree are quite similar to each other, whereas observations that fuse close to the top of the tree will tend to be quite different.

The term **hierarchical** refers to the fact that clusters obtained by cutting the dendrogram at a given height are necessarily nested within the clusters obtained by cutting the dendrogram at any greater height.

#### The Hierarchical Clustering Algorithm

1. Begin with $n$ observations and a measure (such as Euclidean distance) of all the $n(n − 1)/2$ pairwise dissimilarities. Treat each 2 observation as its own cluster.

2. For $i=n,n−1,...,2$

    a) Examine all pairwise inter-cluster dissimilarities among the $i$ clusters and identify the **pair of clusters that are least dissimilar** (that is, most similar). Fuse these two clusters. The dissimilarity between these two clusters indicates the height in the dendrogram at which the fusion should be placed.
    
    b) Compute the new pairwise inter-cluster dissimilarities among the $i−1$ remaining clusters.

#### The linkage function

We have a concept of the dissimilarity between pairs of observations, but how do we define the dissimilarity between two clusters if one or both of the clusters contains multiple observations?

The concept of dissimilarity between a pair of observations needs to be extended to a pair of groups of observations. This extension is achieved by developing the notion of **linkage**, which defines the dissimilarity between two groups of observations.

The four most common types of linkage are:

1. **Complete**: Maximal intercluster dissimilarity. Compute all pairwise dissimilarities between the observations in cluster A and the observations in cluster B, and record the largest of these dissimilarities.
2. **Single**: Minimal intercluster dissimilarity. Compute all pairwise dissimilarities between the observations in cluster A and the observations in cluster B, and record the smallest of these dissimilarities. 
3. **Average**: Mean intercluster dissimilarity. Compute all pairwise dissimilarities between the observations in cluster A and the observations in cluster B, and record the average of these dissimilarities.
4. **Centroid**: Dissimilarity between the centroid for cluster A (a mean vector of length p) and the centroid for cluster B. Centroid linkage can result in undesirable inversions.

Average, complete, and single linkage are most popular among statisticians. Average and complete linkage are generally preferred over single linkage, as they tend to yield more balanced dendrograms. Centroid linkage is often used in genomics, but suffers from a major drawback in that an inversion can occur, whereby two clusters are fused at a height below either of the individual clusters in the dendrogram. This can lead to difficulties in visualization as well as in interpretation of the dendrogram.


```python
# Init
linkages = [hierarchy.complete(X), hierarchy.average(X), hierarchy.single(X)]
titles = ['Complete Linkage', 'Average Linkage', 'Single Linkage']
```


```python
fig, (ax1,ax2,ax3) = plt.subplots(1,3, figsize=(18,10))

# Plot
for linkage, t, ax in zip(linkages, titles, (ax1,ax2,ax3)):
    dendrogram(linkage, ax=ax)
    ax.set_title(t)
```

For this data, complete and average linkage generally separates the observations into their correct groups.

c9d85aa3b6c0835835356452dbe0efaafcd32ac2
ipynb
Jupyter Notebook
10_unsupervised.ipynb
albarran/Machine-Learning-for-Economic-Analysis-2020
eaca29efb8347e51178bdb7fcf90f934d7fee144
MIT
2021-05-14T17:23:16.000Z
2021-05-14T17:23:16.000Z
10_unsupervised.ipynb
albarran/Machine-Learning-for-Economic-Analysis-2020
eaca29efb8347e51178bdb7fcf90f934d7fee144
MIT
10_unsupervised.ipynb
albarran/Machine-Learning-for-Economic-Analysis-2020
eaca29efb8347e51178bdb7fcf90f934d7fee144
MIT
Qwen/Qwen-72B
1. YES
2. YES