# Import some basic libraries
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set_context('paper')
# Fixing the random seed so that we all get the same data
np.random.seed(123456)
# Make synthetic dataset for clustering
num_clusters_true = 3
# The means of each cluster
mu_true = 5.0 * np.random.randn(num_clusters_true, 2)
# The variance of the observations around the cluster
# Covariance of each cluster
Sigmas = []
for i in range(num_clusters_true):
u = 0.5 * np.random.randn(2, 1)
v = 0.5 * np.random.randn(2, 1)
Sigma = np.dot(u, u.T) + np.dot(v, v.T)
Sigmas.append(Sigma)
# How many observations to generate per cluster
num_obs_cluster = [50, 50, 50]
#num_obs_cluster = [50] * num_clusters_true
# Generate the data
data = []
for i in range(num_clusters_true):
x_i = np.random.multivariate_normal(mu_true[i], Sigmas[i],
size=num_obs_cluster[i])
data.append(x_i)
data = np.vstack(data)
# Permute the data so that order info is lost
data = np.random.permutation(data)
Visualize the data:
fig, ax = plt.subplots(dpi=150)
ax.plot(data[:, 0], data[:, 1], '.')
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
Text(0, 0.5, '$x_2$')
We are going to use the sciki-learn implementation of Gaussian mixtures. Here is how you can set up a model and fit it:
# Import that Gaussian mixtures class
from sklearn.mixture import GaussianMixture
# Initialize and fit
model = GaussianMixture(n_components=3).fit(data)
Note that we picked 3 components for the mixture (which happens to be the true number of mixtures). Here are the three mean vectors identified by the model:
model.means_
array([[-7.46929254, -5.73926691],
[ 6.23874412, -0.93769159],
[ 2.40053388, -1.3058076 ]])
Compare these to the true mean vectors:
mu_true
array([[ 2.3455615 , -1.41431672],
[-7.54529252, -5.67816186],
[ 6.06056013, -0.86607325]])
The model identifies one covariance matrix per mixture component. These are stored in a 3D array (first dimension = mixture component). Here is how to access it:
model.covariances_
array([[[ 0.22431117, -0.33125352],
[-0.33125352, 0.53270818]],
[[ 0.46188069, -0.1988742 ],
[-0.1988742 , 0.13489411]],
[[ 0.20382394, 0.53776024],
[ 0.53776024, 1.94763635]]])
Remember that the Gaussian mixture model is essentially fitting the joint probability density function of the data, say $p(\mathbf{x})$. In 2D, you can plot the contours of this probability density function. As a matter of fact, we are going to plot the contours of the logarithm of the probability function, i.e., the contours of $\log p(\mathbf{x})$. Here is how to do it.
# Pick a bunch of x points
x = np.linspace(-9., 8.)
# Pick a bunch of y points
y = np.linspace(-9., 8.)
# Make a grid (X and Y are 50 x 50 matrices.
# The coordinates of the i, j grid point Are X[i, j], Y[i, j]
X, Y = np.meshgrid(x, y)
# Put the grid coordinates in a (50^2) x 2 array because this is what scikit likes
XX = np.array([X.ravel(), Y.ravel()]).T
# Get the log probability on the grid points
z = model.score_samples(XX)
# Respahe z
Z = z.reshape(X.shape)
fig, ax = plt.subplots(dpi=150)
c = ax.contour(X, Y, Z, levels=np.linspace(-400, 1.0, 30))
plt.colorbar(c)
plt.scatter(data[:, 0], data[:, 1])
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$');
You can also use the Gaussian mixture model to cluster the data. Clustering is performed as following. For each point $\mathbf{x}$ that you would like to cluster:
- calculate the probability $p_k$ that the point belongs to the $k$-th mixture component, i.e., calculate: $$ p_k = \frac{\pi_kN(\mathbf{x}|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)}{\sum_{k'=1}^K\pi_{k'}N(\mathbf{x}|\boldsymbol{\mu}_{k'},\boldsymbol{\Sigma}_{k'})}. $$
- assign the point to the cluster $k^*$ with the highest probability, i.e., $$ k^* = \arg\max_{k}p_k. $$ Scikit-learn already implements this functionality. Here is how to do it:
# Clustering with this model
labels = model.predict(data)
fig, ax = plt.subplots(dpi=150)
plt.scatter(data[:, 0], data[:, 1], c=labels)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$');
Another very useful thing that you can do is sample from the estimated probability density. This can be implemented through the following procedure:
- sample a mixture component $k$ in $1,\dots,K$, with probability $\pi_k$, i.e., sample from the Categorical random variable $\operatorname{Categorical}(\pi_1,\dots,\pi_K)$.
- then sample an $\mathbf{x}$ from the Gaussian $N(\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)$. Again, scikit-learn implements this for you. Here is how to access it:
# Sampling from the density
fig, ax = plt.subplots(dpi=150)
# We are going to sample 10,000 times
s_data, labels = model.sample(10000)
ax.scatter(s_data[:, 0], s_data[:, 1])
ax.set_xlabel('$x_1$');
ax.set_ylabel('$x_2$')
Text(0, 0.5, '$x_2$')
Alright! Everything works as expected when we use just the right number of components. Let's see what happens when we use the wrong number of components. We are going to start with fewer components than the right number. Let's try 2 components.
# Trying fewer components
model = GaussianMixture(n_components=2).fit(data)
x = np.linspace(-9., 8.)
y = np.linspace(-9., 8.)
X, Y = np.meshgrid(x, y)
XX = np.array([X.ravel(), Y.ravel()]).T
Z = model.score_samples(XX)
Z = Z.reshape(X.shape)
fig, ax = plt.subplots(dpi=150)
c = ax.contour(X, Y, Z, levels=np.linspace(-400, 0.0, 30))
plt.colorbar(c)
plt.scatter(data[:, 0], data[:, 1])
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$');
Observe that the model is grouping two of the components together. It puts them in the same Gaussian. That's the best it can do with two components. This is clearly underfitting.
Let's now try more components than the right number. We are going to do 5. Here we go:
# Trying fewer components
model = GaussianMixture(n_components=5).fit(data)
x = np.linspace(-9., 8.)
y = np.linspace(-9., 8.)
X, Y = np.meshgrid(x, y)
XX = np.array([X.ravel(), Y.ravel()]).T
Z = model.score_samples(XX)
Z = Z.reshape(X.shape)
fig, ax = plt.subplots(dpi=150)
c = ax.contour(X, Y, Z, levels=np.linspace(-400, 0.0, 30))
plt.colorbar(c)
plt.scatter(data[:, 0], data[:, 1])
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
Text(0, 0.5, '$x_2$')
The density doesn't look very bad even though we used more components. (The density will start looking bad if you use more components). How does the clustering look like when we use five componets. Let's see:
# Clustering with this model
labels = model.predict(data)
fig, ax = plt.subplots(dpi=150)
plt.scatter(data[:, 0], data[:, 1], c=labels)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$');
Okay. The model found five clusters (because we asked for them), but we can visually see in this example that it is overfitting.
Let's now apply the Bayesian information criterion (BIC) to select the number of cluster. First, we loop over all models with 1 to 10 components and we calculate the BIC for each one of them:
# Train with different number of components and estimate BIC
bics = []
models = []
for nc in range(1, 10):
m = GaussianMixture(n_components=nc).fit(data)
bics.append(m.bic(data))
models.append(m)
bics = np.array(bics)
Now let's visualize the BICS:
fig, ax = plt.subplots(dpi=150)
ax.bar(range(1, 10), bics)
ax.set_ylabel('BIC Score')
ax.set_xlabel('Number of components');
Remember that the model with the smallest BIC is the best. For this synthetic dataset the procedure selects 3 components which happens to be the right number.
Questions¶
Change the random seed a couple of times (to whatever integer you want) and generate a few more examples.
Change the number of true components in the synthetic dataset to 5 and repeat. Does BIC work?