Unsupervised Learning and Clustering

Week 9 of 16 · Unsupervised · 7 days

Week 09 · Unsupervised

Unsupervised Learning and Clustering

Week 09 · Day 1 of 7

k-means From Scratch

The algorithm in four steps, and the three assumptions it makes silently

Eight weeks of supervised learning: you had the answers and searched for the rule. Now the answers are gone. Nobody labelled these customers, and the question is whether they fall into groups at all.

Unsupervised learning: Finding structure in data with no target column. There is no correct answer to check against, which changes everything about how you validate, and makes it far easier to convince yourself of something that is not there.

k-means, in full

Place k centres at random.
Assign every point to its nearest centre.
Move each centre to the mean of the points assigned to it.
Repeat 2 and 3 until nothing moves.

import numpy as np

rng = np.random.default_rng(0)
X = np.vstack([rng.normal([0, 0], 0.6, (60, 2)),
               rng.normal([4, 4], 0.6, (60, 2)),
               rng.normal([0, 5], 0.6, (60, 2))])

centres = X[rng.choice(len(X), 3, replace=False)]
for step in range(1, 11):
    d = ((X[:, None] - centres[None, :]) ** 2).sum(axis=2)
    labels = d.argmin(axis=1)
    new = np.array([X[labels == k].mean(axis=0) for k in range(3)])
    shift = np.abs(new - centres).max()
    centres = new
    if step <= 3 or shift < 1e-6:
        print('step %2d  largest centre movement %.6f' % (step, shift))
    if shift < 1e-6:
        break

print('\nfinal centres:')
print(np.round(centres[centres[:, 0].argsort()], 2))

step  1  largest centre movement 2.336408
step  2  largest centre movement 1.994041
step  3  largest centre movement 0.043246
step  4  largest centre movement 0.000000

final centres:
[[-0.    4.94]
 [-0.    0.1 ]
 [ 3.83  4.02]]

Converged in a handful of passes, and the centres landed on the three clusters the data was built from. Twelve lines, no library.

The library version, and the argument that matters

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')

NUM = ['tenure_months', 'support_calls', 'monthly_charges']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
from sklearn.pipeline import make_pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler

prep = make_pipeline(SimpleImputer(strategy='median'), StandardScaler())
Z = prep.fit_transform(df[NUM])
from sklearn.cluster import KMeans
import numpy as np

km = KMeans(n_clusters=4, n_init=10, random_state=42).fit(Z)
print('cluster sizes:', np.bincount(km.labels_))
print('inertia: %.1f' % km.inertia_)
print('iterations to converge:', km.n_iter_)

tiny = np.bincount(km.labels_).argmin()
print('\nthe smallest cluster, in full:')
print(df.loc[km.labels_ == tiny,
             ['customer_id', 'tenure_months', 'support_calls',
              'monthly_charges']].to_string(index=False))

cluster sizes: [1641  822  536    1]
inertia: 3012.8
iterations to converge: 8

the smallest cluster, in full:
customer_id  tenure_months  support_calls  monthly_charges
     C00008             48             97            76.42

One of your four segments has a single member

It is the customer with 97 support calls from week 2. k-means partitions everything, so a point far from all the others gets its own centre rather than being flagged as unusual, and you have spent a quarter of your segmentation on one person. This is the “every point belongs somewhere” assumption failing in the most visible way possible, and it is why day 4's DBSCAN, which can label a point as noise, matters.

n_init exists because k-means gets stuck

The result depends on where the centres started. A bad start converges to a bad local optimum, and the algorithm has no way to notice. n_init=10 runs the whole thing ten times from different starts and keeps the lowest inertia. Setting n_init=1 to save time is how you get a different segmentation every Monday.

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')

NUM = ['tenure_months', 'support_calls', 'monthly_charges']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
from sklearn.pipeline import make_pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler

prep = make_pipeline(SimpleImputer(strategy='median'), StandardScaler())
Z = prep.fit_transform(df[NUM])
from sklearn.cluster import KMeans
import numpy as np

print('%8s %14s %14s' % ('n_init', 'best inertia', 'worst inertia'))
for n_init in [1, 5, 25]:
    results = [KMeans(n_clusters=5, n_init=n_init, random_state=s).fit(Z).inertia_
               for s in range(8)]
    print('%8d %14.1f %14.1f' % (n_init, min(results), max(results)))

  n_init   best inertia  worst inertia
       1         2498.7         2628.4
       5         2498.7         2498.7
      25         2498.7         2498.7

Scaling decides the answer

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')

NUM = ['tenure_months', 'support_calls', 'monthly_charges']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
from sklearn.cluster import KMeans
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler
import numpy as np

PAIR = ['monthly_charges', 'total_charges']
raw = SimpleImputer(strategy='median').fit_transform(df[PAIR])
scaler = StandardScaler().fit(raw)

for name, data in [('unscaled', raw), ('scaled  ', scaler.transform(raw))]:
    centres = KMeans(n_clusters=3, n_init=10, random_state=0).fit(data).cluster_centers_
    if name.strip() == 'scaled':
        centres = scaler.inverse_transform(centres)
    order = centres[:, 1].argsort()
    print('%s centres (monthly, total):' % name)
    print(np.round(centres[order], 1))

unscaled centres (monthly, total):
[[  53.1  665.4]
 [  69.5 2118.9]
 [  75.9 4450.9]]
scaled   centres (monthly, total):
[[  28.9  590.1]
 [  70.3 1101.6]
 [  75.2 3673. ]]

Total charges runs to 7,600 and monthly charges to 106, so unscaled the distance is almost entirely about the total. The three clusters differ hugely in total charges, 665, 2,119, 4,451, and barely at all in monthly charges, 53 to 76. Monthly charges has been ignored.

Scaled, the first cluster's monthly charge drops to 29, which is the no-internet group finally being allowed to separate itself. Neither answer is wrong arithmetic; only one of them is the question you meant to ask.

What k-means assumes

Assumption	Consequence when false
Clusters are roughly spherical	Elongated groups get cut in half
Clusters are similar in size	Large ones get split, small ones absorbed
Every point belongs somewhere	Outliers drag centres toward themselves
You know k	You do not, and the algorithm will not tell you

import numpy as np
from sklearn.cluster import KMeans
from sklearn.datasets import make_moons
from sklearn.metrics import adjusted_rand_score

X, truth = make_moons(n_samples=400, noise=0.06, random_state=0)
labels = KMeans(n_clusters=2, n_init=10, random_state=0).fit_predict(X)
print('two interleaved crescents, k-means agreement with the truth: %.3f'
      % adjusted_rand_score(truth, labels))
print('(1.0 is perfect, 0.0 is random)')

two interleaved crescents, k-means agreement with the truth: 0.274
(1.0 is perfect, 0.0 is random)

k-means fails completely here, because it can only draw straight boundaries between centres and a crescent is not a ball. Day 4's DBSCAN solves exactly this.

Day 1 takeaway

k-means alternates between assigning points to the nearest centre and moving centres to the mean. It depends on the starting position, so keep n_init high. It depends completely on scaling. And it assumes round, similarly sized clusters, when that is false it fails without warning you.

Week 09 · Day 2 of 7

Choosing k, and Whether There Are Clusters At All

Why inertia cannot answer it, and what silhouette and stability can

By Pawan Rama Mali911 words

k-means needs you to supply k, and the whole point of clustering is that you do not know it. There are three honest answers and one dishonest one.

The dishonest one: inertia

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')

NUM = ['tenure_months', 'support_calls', 'monthly_charges']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
from sklearn.pipeline import make_pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler

prep = make_pipeline(SimpleImputer(strategy='median'), StandardScaler())
Z = prep.fit_transform(df[NUM])
from sklearn.cluster import KMeans

print('%4s %14s' % ('k', 'inertia'))
for k in range(1, 11):
    km = KMeans(n_clusters=k, n_init=10, random_state=42).fit(Z)
    print('%4d %14.1f' % (k, km.inertia_))

   k        inertia
       9000.0
       6818.7
       4748.4
       3012.8
       2498.7
       2151.8
       1888.1
       1659.7
       1496.8
       1377.3

Inertia always falls, so it can never choose k

Adding a centre can only reduce the distance from points to their nearest centre. At k equal to the number of rows, inertia is zero and the clustering is useless. The “elbow method” asks you to spot where the curve bends, which on real data is usually a smooth arc with no bend at all. Look at the numbers above and try to name the elbow.

Silhouette: how well separated are the clusters

Silhouette score: For each point, compare its average distance to its own cluster against its average distance to the nearest other cluster. Runs from −1 to 1. Above 0.5 is strong structure; near 0 means the clusters overlap; negative means points are closer to another cluster than their own.

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')

NUM = ['tenure_months', 'support_calls', 'monthly_charges']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
from sklearn.pipeline import make_pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler

prep = make_pipeline(SimpleImputer(strategy='median'), StandardScaler())
Z = prep.fit_transform(df[NUM])
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score, davies_bouldin_score

print('%4s %14s %18s' % ('k', 'silhouette', 'davies-bouldin'))
for k in range(2, 9):
    labels = KMeans(n_clusters=k, n_init=10, random_state=42).fit_predict(Z)
    print('%4d %14.4f %18.4f'
          % (k, silhouette_score(Z, labels), davies_bouldin_score(Z, labels)))

   k     silhouette     davies-bouldin
       0.3201             1.2820
       0.3225             0.8410
       0.3706             0.7130
       0.2970             0.8996
       0.2970             0.9680
       0.2951             0.9170
       0.2965             0.9071

Higher silhouette is better; lower Davies-Bouldin is better. Read them together, and treat agreement between two different criteria as weak evidence rather than proof.

A low silhouette everywhere means there are no clusters

Which is a perfectly respectable finding, and the one most people refuse to accept. Our customers were generated from continuous distributions with no group structure at all, so any k-means result here is a partition of a cloud, not a discovery of groups. Partitions can still be useful for targeting, but calling them “segments we found” is a claim the data does not support.

Check against data that genuinely has clusters

import numpy as np
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs
from sklearn.metrics import silhouette_score

X, _ = make_blobs(n_samples=900, centers=4, cluster_std=0.8,
                  n_features=3, random_state=0)
print('data built from four real clusters:')
for k in range(2, 8):
    labels = KMeans(n_clusters=k, n_init=10, random_state=0).fit_predict(X)
    print('  k=%d  silhouette %.4f' % (k, silhouette_score(X, labels)))

data built from four real clusters:
  k=2  silhouette 0.5886
  k=3  silhouette 0.6695
  k=4  silhouette 0.6581
  k=5  silhouette 0.5398
  k=6  silhouette 0.4460
  k=7  silhouette 0.3641

A clear peak at the true k, and values far above anything the customer data produced. That contrast is the calibration you need: this is what real structure looks like.

Stability: would you get the same answer next month?

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')

NUM = ['tenure_months', 'support_calls', 'monthly_charges']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
from sklearn.pipeline import make_pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler

prep = make_pipeline(SimpleImputer(strategy='median'), StandardScaler())
Z = prep.fit_transform(df[NUM])
from sklearn.cluster import KMeans
from sklearn.metrics import adjusted_rand_score
import numpy as np

rng = np.random.default_rng(0)
print('%4s %22s' % ('k', 'agreement across samples'))
for k in [2, 3, 4, 6]:
    scores = []
    for _ in range(6):
        idx_a = rng.choice(len(Z), int(0.8 * len(Z)), replace=False)
        idx_b = rng.choice(len(Z), int(0.8 * len(Z)), replace=False)
        shared = np.intersect1d(idx_a, idx_b)
        la = KMeans(k, n_init=10, random_state=0).fit(Z[idx_a])
        lb = KMeans(k, n_init=10, random_state=0).fit(Z[idx_b])
        scores.append(adjusted_rand_score(la.predict(Z[shared]),
                                          lb.predict(Z[shared])))
    print('%4d %22.4f' % (k, np.mean(scores)))

   k agreement across samples
               0.8236
               0.8786
               0.7375
               0.7785

Stability is the most practical criterion

Cluster two overlapping samples and see whether the same customers end up together. A segmentation that reshuffles when you add a month of data cannot be the basis of a marketing programme, whatever its silhouette says. This test costs nothing and almost nobody runs it.

Day 2 takeaway

Inertia always falls with k and cannot choose it. Use silhouette and Davies-Bouldin together, compare against what genuine clusters score, and test whether the clustering survives resampling. If every score is weak, the honest conclusion is that there are no clusters, and saying so is a real result.

Week 09 · Day 3 of 7

Hierarchical Clustering

Every k at once, and how linkage decides the shape of what you find

By Pawan Rama Mali687 words

k-means makes you choose k in advance. Hierarchical clustering builds the whole nested family of solutions and lets you cut it wherever you like.

How it merges

Agglomerative clustering: Start with every point as its own cluster. Repeatedly merge the two closest clusters until one remains. The record of merges is a tree, a dendrogram, and cutting it at any height gives you a clustering.

import numpy as np
from scipy.cluster.hierarchy import linkage, dendrogram, fcluster

rng = np.random.default_rng(0)
X = np.vstack([rng.normal([0, 0], 0.4, (12, 2)),
               rng.normal([3, 3], 0.4, (12, 2))])

Zl = linkage(X, method='ward')
print('the last five merges (cluster a, cluster b, distance, size):')
print(np.round(Zl[-5:], 3))

for k in [2, 3, 4]:
    labels = fcluster(Zl, k, criterion='maxclust')
    print('cut into %d: sizes %s' % (k, np.bincount(labels)[1:]))

the last five merges (cluster a, cluster b, distance, size):
[[20.    41.     1.004  9.   ]
 [38.    39.     1.054  9.   ]
 [40.    43.     1.557 12.   ]
 [36.    42.     1.687 12.   ]
 [44.    45.    15.417 24.   ]]
cut into 2: sizes [12 12]
cut into 3: sizes [12  3  9]
cut into 4: sizes [3 9 3 9]

Look at the distance column: the final merge happens at a far larger distance than the ones before it. That jump is the signature of two genuinely separate groups being forced together, and it is how you read a dendrogram.

Linkage decides the shape of what you find

Linkage	Distance between clusters	Tends to produce
`ward`	Increase in within-cluster variance	Compact, similar-sized clusters
`complete`	Furthest pair	Compact, sensitive to outliers
`average`	Mean over all pairs	A compromise
`single`	Closest pair	Long chains, can follow shapes, often produces one giant blob

import numpy as np
from sklearn.cluster import AgglomerativeClustering
from sklearn.datasets import make_moons
from sklearn.metrics import adjusted_rand_score

X, truth = make_moons(n_samples=400, noise=0.06, random_state=0)
for method in ['ward', 'complete', 'average', 'single']:
    labels = AgglomerativeClustering(n_clusters=2, linkage=method).fit_predict(X)
    print('%-9s agreement with the truth %.3f'
          % (method, adjusted_rand_score(truth, labels)))

ward      agreement with the truth 0.239
complete  agreement with the truth 0.285
average   agreement with the truth 0.329
single    agreement with the truth 1.000

Single linkage solves the crescents

Because it only asks whether two clusters have any close pair, it can follow a curved shape point by point along its length. Ward and complete linkage insist on compactness and cut the crescents in half, exactly as k-means did. The same chaining behaviour makes single linkage fragile on noisy data, where one bridging point merges two real clusters.

On the customers

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')

NUM = ['tenure_months', 'support_calls', 'monthly_charges']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
from sklearn.pipeline import make_pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler

prep = make_pipeline(SimpleImputer(strategy='median'), StandardScaler())
Z = prep.fit_transform(df[NUM])
from sklearn.cluster import AgglomerativeClustering, KMeans
from sklearn.metrics import silhouette_score, adjusted_rand_score
import numpy as np

sample = Z[np.random.default_rng(0).choice(len(Z), 1200, replace=False)]

for k in [3, 4, 5]:
    ward = AgglomerativeClustering(n_clusters=k, linkage='ward').fit_predict(sample)
    km = KMeans(k, n_init=10, random_state=0).fit_predict(sample)
    print('k=%d  ward silhouette %.4f   kmeans %.4f   they agree %.3f'
          % (k, silhouette_score(sample, ward), silhouette_score(sample, km),
             adjusted_rand_score(ward, km)))

k=3  ward silhouette 0.3314   kmeans 0.3173   they agree 0.452
k=4  ward silhouette 0.2958   kmeans 0.3727   they agree 0.400
k=5  ward silhouette 0.2636   kmeans 0.3017   they agree 0.371

It does not scale

Agglomerative clustering computes distances between every pair, which is quadratic in memory and worse in time. Three thousand rows is fine; a million is impossible. That is why the snippet above takes a sample, and why k-means remains the default for large data despite its assumptions.

Day 3 takeaway

Hierarchical clustering builds every k at once and you cut the tree where the merge distances jump. Linkage is the real choice: ward and complete find compact blobs, single linkage follows shapes but chains through noise. It is quadratic, so it is a tool for thousands of rows, not millions.

Week 09 · Day 4 of 7

DBSCAN and Density

Clusters as dense regions, arbitrary shapes, and outliers named as noise

By Pawan Rama Mali765 words

k-means partitions everything, including the points that belong to no group. DBSCAN starts from a different idea: a cluster is a dense region, and everything else is noise.

DBSCAN: Density-Based Spatial Clustering of Applications with Noise. A point is a core point if at least min_samples points lie within eps of it. Core points that are close enough form a cluster; points near a cluster join it; everything else is labelled noise, as −1.

It solves the crescents

import numpy as np
from sklearn.cluster import DBSCAN, KMeans
from sklearn.datasets import make_moons
from sklearn.metrics import adjusted_rand_score

X, truth = make_moons(n_samples=400, noise=0.06, random_state=0)

km = KMeans(2, n_init=10, random_state=0).fit_predict(X)
db = DBSCAN(eps=0.25, min_samples=5).fit_predict(X)
print('k-means agreement %.3f' % adjusted_rand_score(truth, km))
print('DBSCAN  agreement %.3f' % adjusted_rand_score(truth, db))
print('clusters found: %d, noise points: %d'
      % (len(set(db)) - (1 if -1 in db else 0), (db == -1).sum()))

k-means agreement 0.274
DBSCAN  agreement 1.000
clusters found: 2, noise points: 0

No k required, arbitrary shapes handled, and outliers named rather than forced into a group. The cost is that eps now has to be chosen, and it is less forgiving than k.

eps is the whole difficulty

import numpy as np
from sklearn.cluster import DBSCAN
from sklearn.datasets import make_moons
from sklearn.metrics import adjusted_rand_score

X, truth = make_moons(n_samples=400, noise=0.06, random_state=0)
print('%8s %10s %10s %12s' % ('eps', 'clusters', 'noise', 'agreement'))
for eps in [0.05, 0.12, 0.25, 0.4, 1.0]:
    lab = DBSCAN(eps=eps, min_samples=5).fit_predict(X)
    n_clusters = len(set(lab)) - (1 if -1 in lab else 0)
    print('%8s %10d %10d %12.3f'
          % (eps, n_clusters, (lab == -1).sum(),
             adjusted_rand_score(truth, lab)))

     eps   clusters      noise    agreement
05         28        218        0.022
12          2          2        0.990
25          2          0        1.000
4          1          0        0.000
0          1          0        0.000

Too small and everything is noise; too large and everything is one cluster

At 0.05 almost every point fails the density test. At 1.0 the two crescents merge. The usable window here runs from about 0.12 to 0.4, and on real data it can be much narrower. This sensitivity is DBSCAN's real weakness.

Choosing eps from the data

import numpy as np
from sklearn.neighbors import NearestNeighbors
from sklearn.datasets import make_moons

X, _ = make_moons(n_samples=400, noise=0.06, random_state=0)

# Distance to the 5th nearest neighbour, sorted. The 'knee' in that curve
# is roughly where points stop being in dense regions.
d, _ = NearestNeighbors(n_neighbors=5).fit(X).kneighbors(X)
fifth = np.sort(d[:, 4])
print('percentile of 5th-neighbour distance:')
for q in [50, 75, 90, 95, 99]:
    print('  %3d%%  %.4f' % (q, np.percentile(fifth, q)))
print('\na reasonable eps sits near the 90th to 95th percentile')

percentile of 5th-neighbour distance:
   50%  0.0642
   75%  0.0837
   90%  0.1051
   95%  0.1179
   99%  0.1489

a reasonable eps sits near the 90th to 95th percentile

On the customers, where density is uniform

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')

NUM = ['tenure_months', 'support_calls', 'monthly_charges']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
from sklearn.pipeline import make_pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler

prep = make_pipeline(SimpleImputer(strategy='median'), StandardScaler())
Z = prep.fit_transform(df[NUM])
from sklearn.cluster import DBSCAN
import numpy as np

print('%8s %10s %10s %12s' % ('eps', 'clusters', 'noise', 'largest'))
for eps in [0.3, 0.5, 0.8, 1.2]:
    lab = DBSCAN(eps=eps, min_samples=10).fit_predict(Z)
    n_clusters = len(set(lab)) - (1 if -1 in lab else 0)
    sizes = np.bincount(lab[lab >= 0]) if n_clusters else np.array([0])
    print('%8s %10d %10d %12d'
          % (eps, n_clusters, (lab == -1).sum(), sizes.max()))

     eps   clusters      noise      largest
3         11        411          905
5          1         43         2957
8          1          4         2996
2          1          2         2998

Either almost everything is noise, or almost everything is one enormous cluster. There is no setting that produces several balanced groups, because there are no dense regions separated by sparse ones. DBSCAN is telling you the same thing the silhouette scores did in day 2, and it is telling you more clearly.

	k-means	Hierarchical	DBSCAN
Need k in advance	Yes	No	No
Cluster shapes	Spherical	Depends on linkage	Any
Handles outliers	No	Poorly	Yes, explicitly
Scales to millions	Yes	No	Moderately
Main difficulty	Choosing k	Choosing linkage	Choosing eps

Day 4 takeaway

DBSCAN defines clusters as dense regions, so it needs no k, handles any shape, and labels outliers as noise instead of absorbing them. In exchange eps is sensitive and can be chosen from the k-nearest-neighbour distance curve. When no eps gives sensible groups, that is evidence about your data.

Week 09 · Day 5 of 7

Gaussian Mixture Models

Soft assignment, cluster shape, and a selection criterion with a real minimum

By Pawan Rama Mali786 words

k-means gives every point one label. A Gaussian mixture gives every point a probability of belonging to each cluster, which is both more honest and more useful.

Gaussian mixture model: Assumes the data was generated by several Gaussian distributions mixed together, and fits their means, covariances and weights. Each point gets a membership probability per component rather than a single assignment, soft rather than hard clustering.

Soft assignment

import numpy as np
from sklearn.mixture import GaussianMixture

rng = np.random.default_rng(0)
X = np.vstack([rng.normal([0, 0], 1.0, (300, 2)),
               rng.normal([3, 3], 1.0, (300, 2))])

gm = GaussianMixture(n_components=2, random_state=0).fit(X)
probs = gm.predict_proba(X)

confident = (probs.max(axis=1) > 0.95).mean()
print('points assigned with >95%% confidence: %.1f%%' % (100 * confident))
print('points that are genuinely ambiguous:  %.1f%%'
      % (100 * (probs.max(axis=1) < 0.7).mean()))

border = np.argsort(np.abs(probs[:, 0] - 0.5))[:3]
print('\nmost ambiguous points:')
for i in border:
    print('  at %s  probabilities %s' % (X[i].round(2), probs[i].round(3)))

points assigned with >95% confidence: 92.3%
points that are genuinely ambiguous:  1.8%

most ambiguous points:
  at [0.48 2.69]  probabilities [0.51 0.49]
  at [1.79 1.27]  probabilities [0.527 0.473]
  at [0.45 2.66]  probabilities [0.466 0.534]

Those borderline customers are real. k-means would have assigned each of them firmly to one side and told you nothing about the doubt.

Covariance shape is the reason to prefer it

import numpy as np
from sklearn.mixture import GaussianMixture
from sklearn.cluster import KMeans
from sklearn.metrics import adjusted_rand_score

rng = np.random.default_rng(0)
# Two elongated, differently oriented clouds -- not spheres.
a = rng.normal(0, 1, (400, 2)) @ np.array([[3.0, 0.0], [0.0, 0.3]])
b = rng.normal(0, 1, (400, 2)) @ np.array([[0.3, 0.0], [0.0, 3.0]]) + [4, 0]
X = np.vstack([a, b])
truth = np.r_[np.zeros(400), np.ones(400)]

km = KMeans(2, n_init=10, random_state=0).fit_predict(X)
gm = GaussianMixture(2, covariance_type='full', random_state=0).fit_predict(X)
print('k-means agreement %.3f' % adjusted_rand_score(truth, km))
print('GMM     agreement %.3f' % adjusted_rand_score(truth, gm))

k-means agreement 0.393
GMM     agreement 0.851

k-means is a special case of a GMM

Fix every covariance to the same spherical shape and give every point a hard assignment, and the mixture model reduces exactly to k-means. covariance_type='full' lets each component have its own shape and orientation, which is what makes it fit the two elongated clouds above.

covariance_type	Each component gets	Parameters
`full`	Its own shape and orientation	Most
`tied`	A shared shape	Fewer
`diag`	Axis-aligned, own scale per axis	Fewer still
`spherical`	One radius	Fewest, nearly k-means

Choosing the number of components properly

Unlike k-means, a mixture model has a likelihood, so you can use information criteria, which penalise complexity and therefore do not improve forever.

import numpy as np
from sklearn.mixture import GaussianMixture
from sklearn.datasets import make_blobs

X, _ = make_blobs(n_samples=900, centers=4, cluster_std=0.9,
                  n_features=3, random_state=0)
print('%4s %14s %14s' % ('k', 'BIC', 'AIC'))
for k in range(1, 9):
    gm = GaussianMixture(k, covariance_type='full', random_state=0).fit(X)
    print('%4d %14.1f %14.1f' % (k, gm.bic(X), gm.aic(X)))

   k            BIC            AIC
      12623.6        12580.3
      10613.2        10521.9
       9782.0         9642.7
       9638.6         9451.3
       9695.3         9460.0
       9761.8         9478.5
       9798.7         9467.3
       9870.6         9491.2

BIC has a minimum, unlike inertia

This is the real advantage over the elbow method. BIC adds a penalty proportional to the number of parameters, so it falls while extra components genuinely explain the data and rises once they are just fitting noise. The minimum is an actual answer rather than a judgement call about where a curve bends.

On the customers

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')

NUM = ['tenure_months', 'support_calls', 'monthly_charges']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
from sklearn.pipeline import make_pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler

prep = make_pipeline(SimpleImputer(strategy='median'), StandardScaler())
Z = prep.fit_transform(df[NUM])
from sklearn.mixture import GaussianMixture

print('%4s %14s' % ('k', 'BIC'))
best = None
for k in range(1, 9):
    gm = GaussianMixture(k, covariance_type='full', random_state=0).fit(Z)
    bic = gm.bic(Z)
    if best is None or bic < best[1]:
        best = (k, bic)
    print('%4d %14.1f' % (k, bic))
print('\nlowest BIC at k = %d' % best[0])

   k            BIC
      25540.7
      21961.8
      20710.0
      20479.7
      19982.6
      19816.1
      19807.6
      19895.4

lowest BIC at k = 7

Day 5 takeaway

A Gaussian mixture gives membership probabilities rather than hard labels, and lets each component have its own shape, so it fits elongated clusters k-means cannot. Because it has a likelihood, BIC can choose the number of components with a genuine minimum, which inertia never has.

Week 09 · Day 6 of 7

Clustering Mixed Data

Why one-hot encoding makes your categories outvote everything else

By Pawan Rama Mali906 words

Every algorithm this week needs distances, and our data is half categorical. What is the distance between a fibre optic customer and a DSL one? The answer you choose changes everything.

The one-hot trap

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')

NUM = ['tenure_months', 'support_calls', 'monthly_charges']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
from sklearn.preprocessing import OneHotEncoder, StandardScaler
from sklearn.impute import SimpleImputer
from sklearn.pipeline import make_pipeline
import numpy as np

cat = make_pipeline(SimpleImputer(strategy='most_frequent'),
                    OneHotEncoder(sparse_output=False)).fit_transform(df[CAT])
num = make_pipeline(SimpleImputer(strategy='median'),
                    StandardScaler()).fit_transform(df[NUM])

print('numeric columns     %d' % num.shape[1])
print('one-hot columns     %d' % cat.shape[1])
print('\ntypical distance contributed by each block:')
rng = np.random.default_rng(0)
i, j = rng.integers(0, len(num), (2, 400))
print('  numeric  %.3f' % np.linalg.norm(num[i] - num[j], axis=1).mean())
print('  one-hot  %.3f' % np.linalg.norm(cat[i] - cat[j], axis=1).mean())

numeric columns     3
one-hot columns     12

typical distance contributed by each block:
  numeric  1.954
  one-hot  2.111

Four categorical variables became twelve columns

Three numeric columns contribute 1.95 of typical distance between two customers; twelve one-hot columns contribute 2.11. So slightly more than half the notion of “similar customer” is now categorical, and worse, it is unevenly spread: contract type gets three columns and payment method four, so payment method counts more than tenure, monthly charges or support calls individually.

Nobody decided that. It fell out of how many levels each variable happened to have.

Three honest options

Weight the blocks so numeric and categorical contribute comparably.
Use a distance built for mixed data, such as Gower's, which handles each column type on its own terms.
Cluster within categories: segment fibre customers separately from DSL ones. Often the most interpretable.

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')

NUM = ['tenure_months', 'support_calls', 'monthly_charges']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
from sklearn.preprocessing import OneHotEncoder, StandardScaler
from sklearn.impute import SimpleImputer
from sklearn.pipeline import make_pipeline
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
import numpy as np

cat = make_pipeline(SimpleImputer(strategy='most_frequent'),
                    OneHotEncoder(sparse_output=False)).fit_transform(df[CAT])
num = make_pipeline(SimpleImputer(strategy='median'),
                    StandardScaler()).fit_transform(df[NUM])

print('%10s %14s %28s' % ('cat weight', 'silhouette', 'largest cluster is'))
for w in [0.0, 0.25, 0.5, 1.0]:
    M = np.hstack([num, cat * w])
    lab = KMeans(4, n_init=10, random_state=0).fit_predict(M)
    big = np.bincount(lab).argmax()
    top = df.loc[lab == big, 'internet_service'].value_counts(normalize=True)
    print('%10.2f %14.4f %20s %6.0f%%'
          % (w, silhouette_score(M, lab), top.index[0], 100 * top.iloc[0]))

cat weight     silhouette           largest cluster is
00         0.3706          Fibre optic     66%
25         0.3403          Fibre optic     67%
50         0.2856          Fibre optic     68%
00         0.2083                  DSL     60%

At weight 0, numeric columns only. The silhouette is 0.37. At full weight, which is what plain one-hot encoding gives you by default, it falls to 0.21 and the dominant service in the largest cluster flips from fibre optic to DSL. Letting the categories in has made the clusters measurably worse and changed what they are about, and the default is the worst setting on the table.

Do not read this as a search for the best weight

The weight that maximises silhouette here is zero, which just says the numeric columns cluster more cleanly on their own. If the categorical variables genuinely matter to your question you cannot discard them because they score badly. You need a distance that treats them fairly, which is the next section. Optimising silhouette over the weight would be choosing what to measure by what scores well.

Gower distance, implemented

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')

NUM = ['tenure_months', 'support_calls', 'monthly_charges']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
import numpy as np

def gower(A, B, numeric_cols, cat_cols, ranges):
    """Mean per-column dissimilarity: scaled absolute difference for
    numeric columns, 0 or 1 for categorical ones."""
    total = np.zeros(len(A))
    for c in numeric_cols:
        total += np.abs(A[c].to_numpy() - B[c].to_numpy()) / ranges[c]
    for c in cat_cols:
        total += (A[c].to_numpy() != B[c].to_numpy()).astype(float)
    return total / (len(numeric_cols) + len(cat_cols))

d = df.dropna(subset=NUM + CAT).reset_index(drop=True)
ranges = {c: d[c].max() - d[c].min() for c in NUM}

a = d.iloc[[0, 0, 0]].reset_index(drop=True)
b = d.iloc[[1, 2, 3]].reset_index(drop=True)
print('customer 0 against three others:')
print(np.round(gower(a, b, NUM, CAT, ranges), 4))
print('\n0 is identical, 1 is maximally different.')
print('Every column contributes exactly once, whatever its type.')

customer 0 against three others:
[0.4508 0.3574 0.7095]

0 is identical, 1 is maximally different.
Every column contributes exactly once, whatever its type.

One column, one vote

That is the whole idea, and it is what one-hot encoding breaks. Gower scales each numeric column by its range so it lands in 0 to 1, scores each categorical column as match or mismatch, and averages. A three-level variable and a seventeen-level variable then carry the same weight, which is almost always what you meant.

Day 6 takeaway

One-hot encoding before clustering gives categorical variables one vote per level, so they dominate the distance and the clusters become a restatement of your categories, while the silhouette score improves, which makes it worse. Weight the blocks, use Gower distance, or cluster within categories.

Week 09 · Day 7 of 7

A Customer Segmentation Project

Features, k, stability, business description, and an honest validation

By Pawan Rama Mali1077 words

A segmentation, done properly, including the part where you check whether it means anything.

1. Decide what the segmentation is for

This decides the features. Segmenting for a retention campaign means behaviour and value; segmenting for network planning means service and usage. There is no general-purpose segmentation, and starting with “cluster the customers” is how you end up with groups nobody can act on.

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')

NUM = ['tenure_months', 'support_calls', 'monthly_charges']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
import numpy as np

# For a retention programme: how long, how much, how unhappy.
FEATURES = ['tenure_months', 'monthly_charges', 'support_calls']
prep = make_pipeline(SimpleImputer(strategy='median'), StandardScaler())
Z = prep.fit_transform(df[FEATURES])
print('clustering on %s' % FEATURES)
print('rows %d, columns %d' % Z.shape)

clustering on ['tenure_months', 'monthly_charges', 'support_calls']
rows 3000, columns 3

2. Choose k with more than one criterion

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')

NUM = ['tenure_months', 'support_calls', 'monthly_charges']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
from sklearn.cluster import KMeans
from sklearn.mixture import GaussianMixture
from sklearn.metrics import silhouette_score, davies_bouldin_score

FEATURES = ['tenure_months', 'monthly_charges', 'support_calls']
Z = make_pipeline(SimpleImputer(strategy='median'),
                  StandardScaler()).fit_transform(df[FEATURES])

print('%4s %12s %14s %14s' % ('k', 'silhouette', 'davies-b', 'BIC'))
for k in range(2, 8):
    lab = KMeans(k, n_init=10, random_state=42).fit_predict(Z)
    bic = GaussianMixture(k, covariance_type='full', random_state=0).fit(Z).bic(Z)
    print('%4d %12.4f %14.4f %14.1f'
          % (k, silhouette_score(Z, lab), davies_bouldin_score(Z, lab), bic))

   k   silhouette       davies-b            BIC
     0.3201         1.2820        21961.8
     0.3225         0.8410        20710.0
     0.3706         0.7130        20479.7
     0.2970         0.8996        19982.6
     0.2970         0.9680        19816.1
     0.2951         0.9170        19807.6

3. Check it is stable

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')

NUM = ['tenure_months', 'support_calls', 'monthly_charges']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
from sklearn.cluster import KMeans
from sklearn.metrics import adjusted_rand_score
import numpy as np

FEATURES = ['tenure_months', 'monthly_charges', 'support_calls']
Z = make_pipeline(SimpleImputer(strategy='median'),
                  StandardScaler()).fit_transform(df[FEATURES])
rng = np.random.default_rng(0)

for k in [3, 4, 5]:
    agree = []
    for _ in range(8):
        ia = rng.choice(len(Z), int(0.8 * len(Z)), replace=False)
        ib = rng.choice(len(Z), int(0.8 * len(Z)), replace=False)
        shared = np.intersect1d(ia, ib)
        a = KMeans(k, n_init=10, random_state=0).fit(Z[ia])
        b = KMeans(k, n_init=10, random_state=0).fit(Z[ib])
        agree.append(adjusted_rand_score(a.predict(Z[shared]),
                                         b.predict(Z[shared])))
    print('k=%d  stability %.4f  +/- %.4f' % (k, np.mean(agree), np.std(agree)))

k=3  stability 0.7134  +/- 0.4151
k=4  stability 0.8612  +/- 0.2177
k=5  stability 0.9153  +/- 0.0507

4. Describe the segments in business language

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')

NUM = ['tenure_months', 'support_calls', 'monthly_charges']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
from sklearn.cluster import KMeans
import pandas as pd

FEATURES = ['tenure_months', 'monthly_charges', 'support_calls']
Z = make_pipeline(SimpleImputer(strategy='median'),
                  StandardScaler()).fit_transform(df[FEATURES])

seg = df.copy()
seg['segment'] = KMeans(4, n_init=10, random_state=42).fit_predict(Z)

profile = seg.groupby('segment').agg(
    customers=('customer_id', 'size'),
    tenure=('tenure_months', 'median'),
    monthly=('monthly_charges', 'median'),
    calls=('support_calls', 'mean'),
    churn_rate=('churned', 'mean'),
)
profile['monthly_value'] = (profile['customers'] * profile['monthly']).round(0)
print(profile.round(2).to_string())

         customers  tenure  monthly  calls  churn_rate  monthly_value
segment                                                              
0             1641    15.0    74.22   1.41        0.40       121787.0
1              822    16.0    27.98   0.82        0.16        23004.0
2              536    52.0    68.40   1.33        0.04        36660.0
3                1    48.0    76.42  97.00        0.00           76.0

Churn rate was not a clustering feature

It was deliberately left out, so the fact that it differs across segments is real information rather than something you built in. If it had been a feature, finding that the segments differ in churn would be circular. Always hold back the variable you want to explain and check it afterwards.

5. Test whether the segments predict anything

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')

NUM = ['tenure_months', 'support_calls', 'monthly_charges']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.pipeline import make_pipeline
from sklearn.cluster import KMeans
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score, StratifiedKFold
import numpy as np

FEATURES = ['tenure_months', 'monthly_charges', 'support_calls']
prep = make_pipeline(SimpleImputer(strategy='median'), StandardScaler())
Z = prep.fit_transform(df[FEATURES])
seg = KMeans(4, n_init=10, random_state=42).fit_predict(Z)
y = df['churned']
cv = StratifiedKFold(5, shuffle=True, random_state=0)

onehot = OneHotEncoder(sparse_output=False).fit_transform(seg.reshape(-1, 1))
print('segment label alone      AUC %.4f'
      % cross_val_score(LogisticRegression(max_iter=1000), onehot, y,
                        cv=cv, scoring='roc_auc').mean())
print('raw features alone       AUC %.4f'
      % cross_val_score(LogisticRegression(max_iter=1000), Z, y,
                        cv=cv, scoring='roc_auc').mean())
print('raw features + segment   AUC %.4f'
      % cross_val_score(LogisticRegression(max_iter=1000),
                        np.hstack([Z, onehot]), y, cv=cv,
                        scoring='roc_auc').mean())

segment label alone      AUC 0.6977
raw features alone       AUC 0.7721
raw features + segment   AUC 0.7735

The segment label carries real information, but less than the features it was derived from, and adding it to those features gains little. That is the usual finding, and it is worth stating: a segmentation is a communication tool, a way to give four groups names that a marketing team can act on. It is rarely a way to improve a model.

Your assignment

Rebuild the segmentation with contract one-hot encoded and added to the features. The silhouette will improve. Then print the contract breakdown of each segment, and decide whether you have discovered anything that df.groupby('contract') would not have told you.

Day 7 takeaway

Choose features from what the segmentation is for. Pick k with several criteria and check stability under resampling. Hold back the outcome you care about so it can validate the segments rather than defining them. And be clear about what a segmentation is: a way to describe customers to people, not usually a way to predict better.