# RUN THIS CELL FIRST
# Imports
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
# sklearn
from sklearn.datasets import load_iris
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
# Better plots
plt.rcParams.update(
{
"lines.markersize": 10, # Big points
"font.size": 15, # Larger font
"xtick.major.size": 5.0, # Bigger xticks
"ytick.major.size": 5.0, # Bigger yticks
}
)
4.1 Plots¶
For plotting in Python there are two standard options. The first is to use matplotlib, and the second is to used seaborn which is based on matplotlib. Seaborn typically creates better-looking plots than matplotlib with less effort, but matplotlib is a lot more customisable.
I have used matplotlib for the vast majority of the plots in the lectures and labs, but you are more than welcome to use seaborn for the labs and tests if you prefer it.
It is also possible to plot through pandas itself if you are working with a dataframe (some examples are here) using matplotlib as a backend. This is quite neat but less flexible than using matplotlib directly.
I'll now provide code to recreate different plots from the lecture. The fig, ax = plt.subplots() that we use to create the figure and its axes looks a bit esoteric, but you can find a good explanation of why we're doing it here.
4.1.1 Bar plots¶
These are useful for visualising the distribution of a categorical variable.
### Bar plot of lecture ratings
fig, ax = plt.subplots()
# Synthesise data
opinions = ["1 Awful", "2 Bad", "3 OK", "4 Good", "5 Great"]
tallies = [12, 5, 2, 10, 20]
# Plot
ax.bar(opinions, tallies)
plt.show()
4.1.2 Histograms¶
These are useful for visualising the distribution of a continuous variable, although you may need to tweak the bin size.
### Histogram of exam results
# Package for creating skewed data
from scipy.stats import skewnorm
# The number of bins.
num_bins = 10
fig, ax = plt.subplots()
# Setting a seed so the plot looks the same each time
np.random.seed(4)
# Skewed data
x = skewnorm.rvs(1, loc=60, scale=10, size=200)
ax.hist(x, bins=num_bins)
ax.set_xlabel("Marks")
ax.set_ylabel("Frequencies")
ax.set_title(f"Exam marks (%) - {num_bins} bins", fontsize=20)
plt.show()
4.1.3 Scatter plots¶
I would argue these are the most important plot of all. They let you compare values for different features in either 2D or 3D.
### Petal length vs. petal width for irises (2D)
data = load_iris()
df = pd.DataFrame(data.data, columns=data.feature_names)
# Add a column with the species as an integer 0,1,2
# We can use different colours for different species
df.insert(4, "species", data.target)
fig, ax = plt.subplots(figsize=[8,8])
ax.scatter(
df["petal length (cm)"][df["species"] == 0],
df["petal width (cm)"][df["species"] == 0],
color="r",
edgecolor='k',
)
ax.scatter(
df["petal length (cm)"][df["species"] == 1],
df["petal width (cm)"][df["species"] == 1],
color="g",
edgecolor='k',
)
ax.scatter(
df["petal length (cm)"][df["species"] == 2],
df["petal width (cm)"][df["species"] == 2],
color="b",
edgecolor='k',
)
ax.legend(data.target_names)
# Tidy up ticks (optional)
ax.set_yticks([0, 1, 2, 3])
# Label
ax.set_title("Petal length vs. petal width of irises")
ax.set_xlabel("Petal length (cm)")
ax.set_ylabel("Petal width (cm)")
plt.show()
### Petal length vs. petal width vs. sepal length for irises (3D)
# We have to make fig, ax in a slightly different way for 3D
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(projection="3d")
ax.scatter(
df["petal length (cm)"][df["species"] == 0],
df["petal width (cm)"][df["species"] == 0],
df["sepal length (cm)"][df["species"] == 0],
color="r",
edgecolor='k',
)
ax.scatter(
df["petal length (cm)"][df["species"] == 1],
df["petal width (cm)"][df["species"] == 1],
df["sepal length (cm)"][df["species"] == 1],
color="g",
edgecolor='k',
)
ax.scatter(
df["petal length (cm)"][df["species"] == 2],
df["petal width (cm)"][df["species"] == 2],
df["sepal length (cm)"][df["species"] == 2],
color="b",
edgecolor='k',
)
ax.legend(data.target_names)
# Label
ax.set_title("Sepal Length vs. Petal length vs. petal width of irises")
ax.set_xlabel("Petal length (cm)")
ax.set_ylabel("Petal width (cm)")
ax.set_zlabel("Sepal length (cm)")
plt.show()
4.1.4 Box plots¶
I do use seaborn for boxplots, as they look a lot better.
### Box plot of petal length for the iris dataset
# import seaborn
import seaborn as sns
fig = plt.figure(figsize=(6, 4))
sns.boxplot(x=df['petal length (cm)'])
plt.show()
### Box plot of petal length for the iris dataset by species
fig = plt.figure(figsize=(7, 8))
sns.boxplot(y=df["petal length (cm)"], x=df["species"])
plt.show()
4.1.5 KDE plots¶
Kernel density estimation (KDE) plots are easiest to do in seaborn.
### KDE plot of petal length for the iris dataset
fig = plt.figure(figsize=(7, 8))
sns.kdeplot(x=df['petal length (cm)'])
plt.show()
4.1.6 Violin plots¶
Also seaborn!
### Violin plot of petal length for the iris dataset by species
fig = plt.figure(figsize=(7, 8))
sns.violinplot(y=df["petal length (cm)"], x=df["species"])
plt.show()
4.1.7 Heat maps¶
Heat maps are a nice way to visualise the magnitude of quantities in a matrix. Here we are going to use seaborn to see how different attributes correlate for a penguin dataset.
### Penguin correlation heatmap
# Read pandas dataset from seaborn package
df = sns.load_dataset("penguins", cache=False)
# Compute correlation matrix and plot using seaborn
# Notice have used .corr() from pandas
matrix = df.corr(numeric_only=True).round(2)
sns.heatmap(matrix, annot=True, vmax=1, vmin=-1, center=0, cmap="vlag").set(
title="Correlations for penguin attributes"
)
plt.show()
4.1.8 Pair plots¶
These are great for a quick overview of a dataset you are considering.
### Penguin pair plot
# The hue argument below does the class separation in the pair plot.
sns.pairplot(df, hue='species')
4.2 Dimensionality reduction¶
Unfortunately, a lot of datasets are high dimensional, and can't really be visualised using the plots we've looked at. One way to resolve this is to use dimensionality reduction techniques to bring a high dimensional dataset down to e.g. 2D where we can just view it on a scatter plot. In the lectures we looked at two dimensionality reduction techniques: PCA and t-SNE.
4.2.1 Principal Component Analysis (PCA)¶
PCA takes a standardised dataset matrix $\mathbf{X}\in \mathbb{R}^{N\times D} $ and returns an orthonormal matrix $\mathbf{W}_{PC}\in \mathbb{R}^{D\times D}$. Note that orthonormal just means that $\mathbf{W}_{PC}^\top \mathbf{W}_{PC} = \mathbf{I}$.
The columns of $\mathbf{W}_{PC}$, which we will denote as $\mathbf{w}_1, \mathbf{w}_2, \dots,\mathbf{w}_D$ are the eigenvectors of the empirical covariance matrix $\mathbf{\Sigma }= \frac{1}{N}\mathbf{X^\top}\mathbf{X}$ sorted from highest to lowest eigenvalue $\lambda_1 , \lambda_2, \dots , \lambda_D$. These eigenvectors are known as the principal components of the data. These are unit vectors (i.e. of length 1).
We are interested in PCA for the purpose of dimensionality reduction. We often have a dataset that is too high dimensional to visualise. We would like a simple linear transform $\mathbf{Z=XW}$ that turns our high dimensional data points $\{\mathbf{x}^{(n)}\}_{n=1}^{N}$(the rows of $\mathbf{X}$) into low dimensional data points $\{\mathbf{z}^{(n)}\}_{n=1}^{N}$ (the rows of $\mathbf{Z}$). Here, $\mathbf{x}\in\mathbb{R}^D$ and $\mathbf{z}\in \mathbb{R}^d$ where $d<D$.
In order to quantify how good this transform is, we create a setup where its transpose can be used to (attempt to) reverse the transformation. This means we can compute $\mathbf{\hat{X}}=\mathbf{ZW^\top}=\mathbf{X}\mathbf{W}\mathbf{W}^\top$ where the rows of $\mathbf{\hat{X}}$ are reconstructions of the original data points $\{\hat{\mathbf{x}}^{(n)}\}_{n=1}^{N}$.
We can now assess how good a transformation $\mathbf{W}$ is by computing the average (squared) distance between points and their reconstructions.
\begin{equation} \frac{1}{N}\sum_{n} \lVert{\mathbf{x}^{(n)}-\mathbf{\hat{x}}^{(n)}}\rVert^2=\frac{1}{N}\sum_{n} \lVert{\mathbf{x}^{(n)}-\mathbf{W}\mathbf{W}^\top\mathbf{x}^{(n)}}\rVert^2 \end{equation}We can write this as an optimisation problem. We want to solve
\begin{equation} \underset{\mathbf{W}}{\mathrm{minimise}}\,\, \, \frac{1}{N}\sum_{n} \lVert{\mathbf{x}^{(n)}-\mathbf{W}\mathbf{W}^\top\mathbf{x}^{(n)}}\rVert^2 \end{equation}As it happens, we actually want to solve the constrained optimisation problem:
\begin{equation}\underset{\mathbf{W}}{\mathrm{minimise}}\,\, \, \frac{1}{N}\sum_{n} \lVert{\mathbf{x}^{(n)}-\mathbf{W}\mathbf{W}^\top\mathbf{x}^{(n)}}\rVert^2 \, \, \, \text{subject to} \, \, \, \mathbf{W}^\top \mathbf{W} = \mathbf{I}\end{equation}The constraint makes the columns of $\mathbf{W}$ orthonormal.
What has this got to do with PCA? Well, the solution to this optimisation problem when projecting down to $d$ dimensions is the matrix whose $d$ columns are the first $d$ principal components of the data.
\begin{equation}\mathbf{W} = \begin{bmatrix}\mathbf{w}_1&\mathbf{w}_2&\dots &\mathbf{w}_d\end{bmatrix}\end{equation}Let's visit the ellipse example in the lecture to see PCA in action.
# Generate points on an unrotated ellipse
a = 10
b = 1
theta = np.linspace(0,2*np.pi,100)
y = a*np.sin(theta)
z = b*np.cos(theta)
x = np.zeros(100)
# Stack into dataset matrix
X = np.vstack((x,y,z)).T
# Rotate
t = np.pi / 4
R = np.array([[np.cos(t), np.sin(t),0],[-np.sin(t), np.cos(t),0],[0,0,1]])
X = X @ R
# Standardise in a single line. This is an alternative to what we did in the last cell
scaler = StandardScaler()
X = scaler.fit_transform(X)
# Finally, plot
fig = plt.figure(figsize=(12,6))
ax = fig.add_subplot(projection='3d')
ax.grid()
ax.set_xlim([-2,2])
ax.set_ylim([-2,2])
ax.set_zlim([-2,2])
ax.set_xticks([])
ax.set_yticks([])
ax.set_zticks([])
ax.scatter(X[:,0],X[:,1],X[:,2],edgecolor='k')
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
# Annoyingly, this x_3 doesn't seem to show below. My apologies.
ax.set_zlabel('$x_3$')
fig.tight_layout()
Here, we have a dataset of points that maps out an ellipse in 3D. We know that ellipses are only 2D so we would hope that if we projected this dataset down to 2D we would get something that looks an awful lot like an ellipse.
Let's use the PCA class from sklearn to compute $\mathbf{W}$ and then $\mathbf{Z}=\mathbf{XW}$. When creating an object from this class, it takes in an argument n_components. This refers to the number of principal components we want to keep. In this example, we want to go down to 2D so we need to keep 2 principal components.
# Import PCA package
# Create a pca object
pca = PCA(n_components=2)
# Calling pca.fit(X) computes W behind the scenes
pca.fit(X)
# Calling pca.transform(X) then computes X@W
Z = pca.transform(X)
# Let's now plot our transformed dataset to see if it makes sense
# Plot Z
fig, ax = plt.subplots(figsize=[6,6])
ax.set_aspect("equal", "box")
ax.grid()
ax.set_xlim([-2.2, 2.2])
ax.set_ylim([-2.2, 2.2])
ax.scatter(Z[:, 0], Z[:, 1],color='g',edgecolor='k')
ax.set_xlabel("$z_1$", size=20)
ax.set_ylabel("$z_2$", size=20)
plt.show()
Great! It's an ellipse! In addition to giving you a transformation that minimises reconstruction error, the principal components (i.e. the features in the transformed space) correspond to (orthogonal) directions of maximum variance. This is why $z_1$ above is the major axis of the ellipse, and $z_2$ is the minor axis because the variance along the first axis is greatest, followed by the variance along the second.
We can quantify how much of the variance has been captured in these dimensions by considering the ratio of the eigenvalues corresponding to the principal components we've used relative to the sum of all eigenvalues of $\boldsymbol\Sigma$. Sklearn computes these automatically:
# The amount of variance captured in the 1st and 2nd principal component
# It captures all of the variance because the third dimension is redundant in describing a 2D shape
pca.explained_variance_ratio_
4.2.2 t-distributed Stochastic neighbour embedding (t-SNE)¶
We saw in the lecture that there were some limitations to using PCA in that it was linear. t-SNE is a very popular non-linear technique that produces a low-dimensional mapping that preserves local structure. Instead of learning some transformation matrix, we directly learn the values of each data point in this mapping.
This means that you can't use t-SNE to map a new point. You would have to run the algorithm again incorporating the new point.
To do this, we first devise an expression for $p_{i,j}$: the probability that two point $\mathbf{x}^{(i)}$ and $\mathbf{x}^{(j)}$ are paired in the original space. In t-SNE we use something Gaussian-like:
$$p_{i,j}= \frac{p_{i|j}+p_{j|i}}{2N}$$where
$$p_{j|i}=\frac{\exp (-\frac{1}{2\sigma_i^{2}}\lVert \mathbf{x}^{(i)}-\mathbf{x}^{(j)}\lVert ^2)}{\sum_{k\neq i}\exp (-\frac{1}{2\sigma_i^{2}}\lVert \mathbf{x}^{(i)}-\mathbf{x}^{(k)}\lVert ^2)}$$and each bandwidth $\sigma_i$ is chosen to match a user-selected perplexity, which controls the effective neighbourhood size. We set $p_{i,i}=0$. This gives high pair probabilities to nearby points.
The next step is to devise an expression for $q_{i,j}$: the probability that the low dimensional embeddings of these points $\mathbf{z}^{(i)}$ and $\mathbf{z}^{(j)}$ are paired. In t-SNE we use a Cauchy distribution which has long tails:
$$q_{i,j}=\frac{(1+ \lVert \mathbf{z}^{(i)}-\mathbf{z}^{(j)}\lVert ^2)^{-1}}{\sum_{k\neq l}(1+ \lVert \mathbf{z}^{(k)}-\mathbf{z}^{(l)}\lVert ^2)^{-1}}$$We also set $q_{i,i}=0$.
If you think these expressions are a bit random and unprincipled, then you are absolutely right. Often empirical effectiveness trumps all else in ML.
Given some pair, we would like the probabilities of pairing to be preserved in both the original space and the embedded space. Given what we know about linear regression, it is tempting to write a loss function that is the sum (or mean) of the distances between probabilities. However, we can do better. Given that we are working with probabilities, we can use the KL divergence:
$$ D_{KL}(P\lVert Q) = \sum_i \sum _j p_{i,j}\log \frac{p_{i,j}}{q_{i,j }}$$which is a measure of how one probability distribution diverges from another. Specifically, we can find the exact embeddings $\mathbf{z}^{(1)},\mathbf{z}^{(2)},\dots,\mathbf{z}^{(N)}$ that minimise this quantity by solving $\underset{\mathbf{z}^{(1)},\mathbf{z}^{(2)},\dots,\mathbf{z}^{(N)}}{\mathrm{minimise}}\ D_{KL}(P\lVert Q)$.
This can be done using good old Gradient descent or something similar. You will compare t-SNE to PCA in the lab.