In [9]:
plt.scatter(X[:, 0], X[:, 1])
plt.show()
import os
import pandas as pd
import numpy as np
from typing import NamedTuple, Tuple
from scipy.special import logsumexp
# Graphic part
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib.patches import Circle, Arc
from IPython.display import HTML
import geopandas as gpd
K-Means and Gaussian Mixtures are two different approaches, which can be applied to clustering problems. K-Means searchs for the best distribution of K clusters, in order to minimize the loss function among all points in dataset. Gaussian Mixtures work on different conception, it searchs for best probability distribution, along K gaussian distributions, each one with its own mean and variance.
We will analyze their outcomes on different cases. We start with a 2 dimensional random dataset, to compare both techniques, with different cluster numbers.
Next, we apply both techniques to real world countries data, in order to separate into clusters and analyze them. For this we will apply two features, in a 2-D dimension dataset, but we will try to analyze different clusters conceptions.
We start by declaring the functions and classes used.
We build classes for K-Means and Gaussian Mixtures and calculate all their mathematics.
For our mixture models, we define their mean, variance and weight components.
We also define a function to create the animations and the other plotting functions.
class GaussianMixture(NamedTuple):
"""Tuple holding a gaussian mixture"""
mu: np.ndarray # (K, d) array - each row corresponds to a gaussian component mean
var: np.ndarray # (K, ) array - each row corresponds to the variance of a component
p: np.ndarray # (K, ) array = each row corresponds to the weight of a component
class Kmeans:
"""Mixture model based on kmeans"""
def __init__(self, X: np.ndarray, K: int, seed: int = 0) -> Tuple[GaussianMixture, np.ndarray]:
"""Initializes the mixture model with random points as initial means and uniform assingments
Args:
X: (n, d) array holding the data
K: number of components
seed: random seed
Returns:
mixture: the initialized gaussian mixture
post: (n, K) array holding the soft counts for all components for all examples
"""
np.random.seed(seed)
n, d = X.shape
p = np.ones(K) / K
# select K random points as initial means
mu = X[np.random.choice(n, K, replace=False)]
var = np.zeros(K)
# Compute variance
for j in range(K):
var[j] = ((X - mu[j])**2).mean()
mixture = GaussianMixture(mu, var, p)
post = np.ones((n, K)) / K
self.X = X
self.K = K
self.n = n
self.d = d
self.cost = 0
self.mixture = mixture
self.post = post
def estep(self) -> np.ndarray:
"""E-step: Assigns each datapoint to the gaussian component with the closest mean
Args:
self.X: (n, d) array holding the data
self.mixture: the current gaussian mixture
Returns:
np.ndarray: (n, K) array holding the soft counts for all components for all examples
"""
post = np.zeros((self.n, self.K))
for i in range(self.n):
tiled_vector = np.tile(self.X[i, :], (self.K, 1))
sse = ((tiled_vector - self.mixture.mu)**2).sum(axis=1)
j = np.argmin(sse)
post[i, j] = 1
self.post = post
def mstep(self) -> Tuple[GaussianMixture, float]:
"""M-step: Updates the gaussian mixture. Each cluster yields a component mean and variance.
Args:
self.X: (n, d) array holding the data
self.post: (n, K) array holding the soft counts for all components for all examples
Returns:
GaussianMixture: the new gaussian mixture
float: the distortion cost for the current assignment
"""
n_hat = self.post.sum(axis=0)
p = n_hat / self.n
cost = 0
mu = np.zeros((self.K, self.d))
var = np.zeros(self.K)
for j in range(self.K):
mu[j, :] = self.post[:, j] @ self.X / n_hat[j]
sse = ((mu[j] - self.X)**2).sum(axis=1) @ self.post[:, j]
cost += sse
var[j] = sse / (self.d * n_hat[j])
self.mixture = GaussianMixture(mu, var, p)
self.cost = cost
class Gmixtures:
"""Mixture model using EM"""
def __init__(self, X: np.ndarray, K: int, seed: int = 0) -> Tuple[GaussianMixture, np.ndarray]:
"""Initializes the mixture model with random points as initial means and uniform assingments
Args:
X: (n, d) array holding the data
K: number of components
seed: random seed
Returns:
mixture: the initialized gaussian mixture
post: (n, K) array holding the soft counts for all components for all examples
"""
np.random.seed(seed)
n, d = X.shape
p = np.ones(K) / K
# select K random points as initial means
mu = X[np.random.choice(n, K, replace=False)]
var = np.zeros(K)
# Compute variance
for j in range(K):
var[j] = ((X - mu[j])**2).mean()
mixture = GaussianMixture(mu, var, p)
post = np.ones((n, K)) / K
self.X = X
self.K = K
self.n = n
self.d = d
self.cost = 0
self.mixture = mixture
self.post = post
def estep(self) -> Tuple[np.ndarray, float]:
"""E-step: Softly assigns each datapoint to a gaussian component
Args:
self.X: (n, d) array holding the data
self.mixture: the current gaussian mixture
Returns:
np.ndarray: (n, K) array holding the soft counts for all components for all examples
float: log-likelihood of the assignment
"""
post = np.zeros((self.n, self.K))
def log_N(x, mixture):
d = x.shape[0]
log_P = -1/(2*mixture.var)*np.linalg.norm(x-mixture.mu, axis=1)**2 - \
d/2 * np.log(2*np.pi*mixture.var)
return log_P
ll = 0
for i in range(self.n):
log_P = log_N(X[i], self.mixture)
f = np.log(self.mixture.p + 1e-16) + log_P
post[i] = np.exp(f - logsumexp(f))
ll += logsumexp(f).sum()
self.post = post
self.cost = ll
def mstep(self, min_variance: float = .25) -> GaussianMixture:
"""M-step: Updates the gaussian mixture by maximizing the log-likelihood of the weighted dataset
Args:
self.X: (n, d) array holding the data
self.post: (n, K) array holding the soft counts for all components for all examples
self.mixture: the current gaussian mixture
min_variance: the minimum variance for each gaussian
Returns:
GaussianMixture: the new gaussian mixture
"""
n_hat = self.post.sum(axis=0)
p = n_hat / self.n
mu = self.mixture.mu
cu = self.X != 0
den = self.post.T @ cu
num = self.post.T @ (cu * self.X)
for k in range(self.K):
for dim in range(self.d):
if den[k, dim] >= 1:
mu[k, dim] = num[k, dim] / den[k, dim]
var = np.zeros(self.K)
p = np.zeros(self.K)
for k in range(self.K):
num = 0
den = 0
for u in range(self.n):
cu = self.X[u, :] != 0
x_cu = self.X[u, cu]
mu_cu = mu[k, cu]
num += self.post[u, k] * np.linalg.norm(x_cu - mu_cu) ** 2
den += self.post[u, k] * sum(cu)
p[k] += self.post[u, k]
var[k] = max(min_variance, num / den)
p[k] /= self.n
self.mixture = GaussianMixture(mu, var, p)
def create_animation(model, frames, interval, axis, title_pos):
# First we set up the figure, the axis, and the plot element for animation
fig, ax = plt.subplots(figsize=(8,8))
ax.set_xlim((axis[0], axis[1]))
ax.set_ylim((axis[2], axis[3]))
line, = ax.plot([], [])
plt.close()
def plot_init():
# initialization function: plot the background of each frame
line.set_data([], [])
return (line)
def plot_circle(model, line):
percent = model.post / model.post.sum(axis=1).reshape(-1, 1)
color = {0: "r",
1: "b",
2: "k",
3: "y",
4: "m",
5: "c"}
colors = list(map(lambda x: color[x], np.argmax(percent, axis=1)))
ax.scatter(X[:,0], X[:,1], c=colors)
for j in range(model.K):
mu = model.mixture.mu[j]
sigma = np.sqrt(model.mixture.var[j])
circle = Circle(mu, sigma, color=color[j], fill=False)
ax.add_patch(circle)
def animate(i):
# animation function. This is called sequentially
ax.clear()
line, = ax.plot([], [])
ax.set_xlim((axis[0], axis[1]))
ax.set_ylim((axis[2], axis[3]))
global model
model.estep()
model.mstep()
plot_circle(model, line)
ax.text(title_pos[0], title_pos[1], f'Step: {i+1}', fontsize=20)
return (line,)
# call the animator
ani = FuncAnimation(fig, animate, init_func=plot_init, frames=frames, interval=interval, repeat=True)
return ani
def plot_models(model_1, model_2=None, df=False, world=False, title=None):
""" Usages of this function:
- model_1, model_2: plot both models outcome
- model_1, title: plot only model_1 outcome
- model_1, df, world, title: plot model_1 in world map
"""
color = {0: "r",
1: "b",
2: "k",
3: "y",
4: "m",
5: "c"}
titles = ['K-Means', 'Gaussian Mixtures']
if world == False:
if model_1 and model_2:
models = [model_1, model_2]
fig, ax = plt.subplots(1, 2, figsize = (16, 7))
for i, model in enumerate(models):
percent = model.post / model.post.sum(axis=1).reshape(-1, 1)
colors = list(map(lambda x: color[x], np.argmax(percent, axis=1)))
ax[i].scatter(model.X[:,0], model.X[:,1], c=colors)
ax[i].set_xlabel('gdp (log10)')
ax[i].set_ylabel('pop (log10)')
ax[i].set_title(titles[i])
for j in range(model.K):
mu = model.mixture.mu[j]
sigma = np.sqrt(model.mixture.var[j])
circle = Circle(mu, sigma, color=color[j], fill=False)
ax[i].add_patch(circle)
elif model_1 and not model_2:
fig, ax = plt.subplots(figsize = (8, 8))
percent = model_1.post / model_1.post.sum(axis=1).reshape(-1, 1)
colors = list(map(lambda x: color[x], np.argmax(percent, axis=1)))
ax.scatter(model_1.X[:,0], model_1.X[:,1], c=colors)
ax.set_xlabel('gdp (log10)')
ax.set_ylabel('pop (log10)')
ax.set_title(title)
for j in range(model_1.K):
mu = model_1.mixture.mu[j]
sigma = np.sqrt(model_1.mixture.var[j])
circle = Circle(mu, sigma, color=color[j], fill=False)
ax.add_patch(circle)
else:
fig, ax = plt.subplots(figsize = (12, 6))
countries = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres"))
countries.plot(color="lightgrey", ax=ax)
percent = model_1.post / model_1.post.sum(axis=1).reshape(-1, 1)
colors = list(map(lambda x: color[x], np.argmax(percent, axis=1)))
ax.scatter(df['long'], df['lat'], c=colors)
ax.set_title(title)
def plot_ground_truth(df, how, model=None, world=False):
if how == 'hemisphere':
color = {'north': "b",
'south': "r"}
elif how == 'continent':
color = {'america': 'y',
'europe': 'm',
'africa': 'k',
'asia': 'r',
'oceania': 'b'}
colors = df.apply(lambda row: color.get(row[how]), axis=1).values
if world == False:
fig, ax = plt.subplots(figsize = (8, 8))
ax.scatter(df['gdp'], df['population'], c=colors)
ax.set_xlabel('gdp (log10)')
ax.set_ylabel('pop (log10)')
else:
fig, ax = plt.subplots(figsize = (12, 6))
countries = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres"))
countries.plot(color="lightgrey", ax=ax)
ax.scatter(df['long'], df['lat'], c=colors)
ax.set_title('Ground truth')
if model:
if how == 'hemisphere':
color = ["r", "b"]
elif how == 'continent':
color = ["r", "b", "k", "y", "m"]
for j in range(model.K):
mu = model.mixture.mu[j]
sigma = np.sqrt(model.mixture.var[j])
circle = Circle(mu, sigma, color=color[j], fill=False)
ax.add_patch(circle)
X = np.loadtxt("toy_data.txt")
plt.scatter(X[:, 0], X[:, 1])
plt.show()
K = 2
model = Kmeans(X, K, 2)
ani = create_animation(model, frames=6, interval=1000, axis = [-8, 12, -8, 8], title_pos=[1, 7])
HTML(ani.to_html5_video())
K = 2
model = Gmixtures(X, K, 2)
ani = create_animation(model, frames=14, interval=500, axis = [-8, 12, -8, 8], title_pos=[1, 7])
HTML(ani.to_html5_video())
K = 4
model = Kmeans(X, K, 0)
ani = create_animation(model, frames=20, interval=200, axis = [-8, 12, -8, 8], title_pos=[1, 7])
HTML(ani.to_html5_video())
K = 4
model = Gmixtures(X, K, 0)
ani = create_animation(model, frames=80, interval=100, axis = [-8, 12, -8, 8], title_pos=[1, 7])
HTML(ani.to_html5_video())
With K = 4, we see a difference between the two techniques. While K-Means tries to distribute evenly the points, among all clusters, Gaussian Mixture foccuses on explaining better the points and their density, and we can have overlapping distributions, each with their own mean and variance. Sometimes Gaussian Mixtures can explain better real problems.
We now proceed to analyze real world data. We will use 2 features, GDP (gross domestic product) and population, of 87 countries, from 2020. These data were obtained from Kaggle.
df_gdp = pd.read_csv('GDP.csv')
df_gdp = df_gdp[['Country Name', '2020']]
df_gdp.set_index('Country Name', inplace=True)
df_pop = pd.read_csv('Countries/Pupulation density by countries.csv')
df_pop = df_pop[['Country', 'Population']]
df_pop.set_index('Country', inplace=True)
df_gdp = df_gdp.rename(index={'Bahamas, The': 'Bahamas',
'Congo, Dem. Rep.': 'Democratic Republic of the Congo',
'Congo, Rep.': 'Republic of the Congo',
'Hong Kong SAR, China': 'Hong Kong',
'St. Kitts and Nevis': 'Saint Kitts and Nevis'})
df_pop = df_pop.rename(index={'Bermuda (UK)': 'Bermuda',
'Puerto Rico (US)': 'Puerto Rico',
'Eswatini (Swaziland)': 'Eswatini',
'Uruguay[note 5]': 'Uruguay',
'Saint Vincent and the Grenadines': 'St. Vincent and the Grenadines'})
df_info = pd.read_csv('others_info.csv')
df_info.set_index('Country Name', inplace=True)
df = df_gdp.copy()
df = df.merge(df_pop, how='left', left_index=True, right_index=True)
df = df.dropna()
df['Population'] = (df['Population'].str.replace(",", "")).astype(int)
df = np.log10(df)
df.columns = ['gdp', 'population']
df = df.merge(df_info, how='left', left_index=True, right_index=True)
df
| gdp | population | hemisphere | division | continent | lat | long | |
|---|---|---|---|---|---|---|---|
| Country Name | |||||||
| Australia | 12.123852 | 7.409086 | south | first | oceania | -25.274398 | 133.775136 |
| Austria | 11.636488 | 6.949517 | north | first | europe | 47.516231 | 14.550072 |
| Burundi | 9.453591 | 7.049822 | south | third | africa | -3.373056 | 29.918886 |
| Belgium | 11.717671 | 7.061620 | north | first | europe | 50.503887 | 4.469936 |
| Benin | 10.194557 | 7.069411 | north | third | africa | 9.307690 | 2.315834 |
| ... | ... | ... | ... | ... | ... | ... | ... |
| United States | 13.322219 | 8.517817 | north | first | america | 37.090240 | -95.712891 |
| St. Vincent and the Grenadines | 8.907129 | 5.043441 | north | third | america | 12.984305 | -61.287228 |
| South Africa | 11.525045 | 7.769193 | south | first | africa | -30.559482 | 22.937506 |
| Zambia | 10.257934 | 7.214982 | south | third | africa | -13.133897 | 27.849332 |
| Zimbabwe | 10.256505 | 7.180688 | south | third | africa | -19.015438 | 29.154857 |
87 rows × 7 columns
fig, ax = plt.subplots(figsize=(8,8))
plt.scatter(df['gdp'], df['population'])
plt.xlabel('gdp (log10)')
plt.ylabel('pop (log10)')
plt.show()
K = 2
X = np.array(df.iloc[:,0:2])
model = Kmeans(X, K, 9)
ani = create_animation(model, frames=7, interval=800, axis = [8, 14, 4, 10], title_pos=[11, 9.5])
HTML(ani.to_html5_video())
model_1 = model
K = 2
X = np.array(df.iloc[:,0:2])
model = Gmixtures(X, K, 9)
ani = create_animation(model, frames=14, interval=500, axis = [8, 14, 4, 10], title_pos=[11, 9.5])
HTML(ani.to_html5_video())
model_2 = model
As stated before, K-Means and Gaussian Mixtures show similar clustering when K equals 2.
By separating the data into two clusters, one interpretation we can make is dividing the globe into north and south hemispheres. Below we show both models representation, as well as real hemisphere distributions.
plot_models(model_1, model_2)
plot_ground_truth(df, how='hemisphere', model=model_1, world=False)
By looking at the pictures above, we can make some assumptions.
Both techniques group the dataset in two clusters, by similarities: one at the top right and one at the bottom left.
By looking at the real hemispheres division, both are spread over the values. However the blue cluster spreads more over extreme sides, while the red cluster is smaller and concentrate more in the center.
For a better visualization, we can plot these points in a world map, by their geographic location.
plot_ground_truth(df, how='hemisphere', world=True)
We can also plot the models outcome in the world map. By analyzing the clusterings outcome, we make interesting insights. Although we did not give information about geographic location, we see that clearly there are similarities between countries in the south hemisphere, as well between countries in the north hemisphere.
Overall, we know that countries in the north are more developed and richer, like USA, European and Asian countries. On the other side, south countries are poorer, like South America and African countries. And we can see this represented by the clusters outcome.
By looking at it we can also see some details. Although Australia is located in the south hemisphere, it is more developed and have similar features of the north hemisphere group. Also there are some mixed countries in South/Central America and Africa. This mean that, althouth, they are located in the south hemisphere, they may show features closer to north hemisphere group.
plot_models(model_1, df=df, world=True, title='K-Means')
K = 5
X = np.array(df.iloc[:,0:2])
model = Kmeans(X, K, 8)
ani = create_animation(model, frames=10, interval=500, axis = [8, 14, 4, 10], title_pos=[10.5, 9.5])
HTML(ani.to_html5_video())
model_3 = model
K = 5
X = np.array(df.iloc[:,0:2])
model = Gmixtures(X, K, 8)
ani = create_animation(model, frames=40, interval=200, axis = [8, 14, 4, 10], title_pos=[10.5, 9.5])
HTML(ani.to_html5_video())
model_4 = model
Now that we have 5 clusters, we can analyze if the model outcomes correspond to real world continents. We start by plotting the models outcome.
plot_models(model_3, model_4)
We can see that both techniques have different distributions over the data. In the Gaussian Mixtures the top group has only 4 points, more isolated from the others.
But how do they correspond to real continents? Let's check in the world map.
plot_ground_truth(df, how='continent', world=True)
plot_models(model_3, df=df, world=True, title='K-Means')
plot_models(model_4, df=df, world=True, title='Gaussian Mixtures')
These results are very interesting. First we shall say that the Gaussian Mixtures showed a better outcome, closer to the real continents. Since it represents all data with gaussian distributions, this method cluster better similar and closer points. The outcome of K-Means technique showed more mixed countries, all over America, Africa and Europe.
Also we can see that many African countries share similarities between themselves, same as European countries. Gaussian Mixtures classified both USA and Australia belonging to the European group, which is understandable, due to both being developed.
It is very important to emphasize here that we are investigating clustering techniques and their information. We are not applying proper classification algorithms, neither training or predicting on different datasets. We are using our ground truth variables only to correlate to achieved outcomes. In real problems we may not have the ground truth clusters and it may be very hard to find real groups or their meaning.
Concluding, it was very interesting to discover all these real similarities, in this study case, based only on two features. For future works we could analyze more features, and their representation by applying clustering on PCA components.