The 14 GPCM Covariance Models¶

This notebook demonstrates the Phase 3 (v0.4 roadmap) addition: all 14 Celeux & Govaert (1995) Gaussian Parsimonious Clustering Models (GPCM) -- the same family mclust implements -- as covariance_type options on GaussianMixture, plus select_model's mclust-style default of sweeping all 14 codes.

See the GPCM Covariance Models user guide for what each of the 14 codes means.

In [1]:

import numpy as np
import torch
import matplotlib.pyplot as plt
import os

os.chdir('../')

from sklearn.datasets import load_iris
from tgmm import GaussianMixture, select_model
from tgmm.covariance.gpcm import GPCM_CODES
from tgmm.plotting import plot_gmm

random_state = 42

import numpy as np import torch import matplotlib.pyplot as plt import os os.chdir('../') from sklearn.datasets import load_iris from tgmm import GaussianMixture, select_model from tgmm.covariance.gpcm import GPCM_CODES from tgmm.plotting import plot_gmm random_state = 42

Part 1: the classic iris petal measurements (2-D)¶

Petal length and width, the two features that most cleanly separate the 3 iris species -- small enough to sweep all 14 models x several K quickly, and 2-D enough to plot the winning model's component ellipses directly.

In [2]:

iris = load_iris()
X_iris = torch.tensor(iris.data[:, [2, 3]], dtype=torch.float64)  # petal length, petal width

plt.figure(figsize=(5, 4))
plt.scatter(X_iris[:, 0], X_iris[:, 1], c=iris.target, cmap='viridis', s=15)
plt.xlabel('petal length (cm)')
plt.ylabel('petal width (cm)')
plt.title('Iris petal measurements (true species)')
plt.show()

iris = load_iris() X_iris = torch.tensor(iris.data[:, [2, 3]], dtype=torch.float64) # petal length, petal width plt.figure(figsize=(5, 4)) plt.scatter(X_iris[:, 0], X_iris[:, 1], c=iris.target, cmap='viridis', s=15) plt.xlabel('petal length (cm)') plt.ylabel('petal width (cm)') plt.title('Iris petal measurements (true species)') plt.show()

In [3]:

result_iris = select_model(
    X_iris, n_components=range(1, 6), criterion='bic',
    random_state=random_state, n_init=1, max_iter=200,
)

df_iris = result_iris.to_dataframe().sort_values('bic')
df_iris.head(10)

result_iris = select_model( X_iris, n_components=range(1, 6), criterion='bic', random_state=random_state, n_init=1, max_iter=200, ) df_iris = result_iris.to_dataframe().sort_values('bic') df_iris.head(10)

Out[3]:

	n_components	covariance_type	bic	aic	icl	log_likelihood	converged
38	3	VVE	347.054572	301.895043	357.870257	-0.906317	True
52	4	VVE	354.818447	294.605741	370.829914	-0.848686	True
35	3	VEE	354.977590	315.839331	369.030985	-0.966131	True
41	3	VVV	355.957404	304.776604	368.395953	-0.902589	True
24	2	VVE	361.178482	331.072129	361.250534	-1.036907	True
49	4	VEE	361.414820	310.234020	376.459484	-0.920780	True
39	3	VEV	364.349955	319.190426	376.750136	-0.963968	True
27	2	VVV	364.579644	331.462656	364.637312	-1.031542	True
63	5	VEE	367.825650	304.602308	393.042373	-0.875341	True
55	4	VVV	369.243180	299.998569	405.605136	-0.846662	True

BIC by model family and K¶

One line per covariance code, BIC against K -- lower is better.

In [4]:

plt.figure(figsize=(9, 5))
for code_name in GPCM_CODES:
    sub = df_iris[df_iris['covariance_type'] == code_name].sort_values('n_components')
    sub = sub.dropna(subset=['bic'])
    if len(sub) == 0:
        continue
    plt.plot(sub['n_components'], sub['bic'], marker='o', markersize=3, label=code_name, alpha=0.7)
plt.xlabel('n_components (K)')
plt.ylabel('BIC (lower is better)')
plt.title('BIC across all 14 GPCM models, iris petal data')
plt.legend(ncol=2, fontsize=8)
plt.show()

best = result_iris.best_
print(f'Winning model: {best.covariance_type_} (K={best.n_components}), BIC={best.bic(X_iris):.2f}')

plt.figure(figsize=(9, 5)) for code_name in GPCM_CODES: sub = df_iris[df_iris['covariance_type'] == code_name].sort_values('n_components') sub = sub.dropna(subset=['bic']) if len(sub) == 0: continue plt.plot(sub['n_components'], sub['bic'], marker='o', markersize=3, label=code_name, alpha=0.7) plt.xlabel('n_components (K)') plt.ylabel('BIC (lower is better)') plt.title('BIC across all 14 GPCM models, iris petal data') plt.legend(ncol=2, fontsize=8) plt.show() best = result_iris.best_ print(f'Winning model: {best.covariance_type_} (K={best.n_components}), BIC={best.bic(X_iris):.2f}')

Winning model: VVE (K=3), BIC=347.05

The winning model's component ellipses¶

In [5]:

plt.figure(figsize=(6, 5))
plot_gmm(X_iris, gmm=best, true_labels=torch.tensor(iris.target), color_by_cluster=True)
plt.title(f'{best.covariance_type_} (K={best.n_components}) fit to iris petal measurements')
plt.show()

plt.figure(figsize=(6, 5)) plot_gmm(X_iris, gmm=best, true_labels=torch.tensor(iris.target), color_by_cluster=True) plt.title(f'{best.covariance_type_} (K={best.n_components}) fit to iris petal measurements') plt.show()

Part 2: a 10-D latent-space-like dataset¶

Four clusters in 10 dimensions with deliberately varied covariance structure (different volumes, some axis-aligned, some correlated/elongated) -- the kind of thing a GPCM sweep is actually for: no single covariance type is obviously right, so the model-selection search has real work to do.

In [6]:

def make_latent_cluster(n, d, seed, volume, elongation=None, rotate=False):
    gen = torch.Generator().manual_seed(seed)
    base = torch.randn(n, d, generator=gen, dtype=torch.float64)
    scale = torch.ones(d, dtype=torch.float64)
    if elongation is not None:
        scale[: len(elongation)] = torch.tensor(elongation, dtype=torch.float64)
    X = base * scale.sqrt() * volume
    if rotate:
        Q, _ = torch.linalg.qr(torch.randn(d, d, generator=gen, dtype=torch.float64))
        X = X @ Q.T
    return X

torch.manual_seed(random_state)
d = 10
centers = torch.randn(4, d, dtype=torch.float64) * 6
clusters = [
    make_latent_cluster(150, d, seed=1, volume=1.0) + centers[0],
    make_latent_cluster(150, d, seed=2, volume=1.8) + centers[1],
    make_latent_cluster(150, d, seed=3, volume=1.0, elongation=[4, 4, 1, 1, 1, 1, 1, 1, 1, 1]) + centers[2],
    make_latent_cluster(150, d, seed=4, volume=1.0, elongation=[4, 4, 1, 1, 1, 1, 1, 1, 1, 1], rotate=True) + centers[3],
]
X_latent = torch.cat(clusters, dim=0)
print(X_latent.shape)

def make_latent_cluster(n, d, seed, volume, elongation=None, rotate=False): gen = torch.Generator().manual_seed(seed) base = torch.randn(n, d, generator=gen, dtype=torch.float64) scale = torch.ones(d, dtype=torch.float64) if elongation is not None: scale[: len(elongation)] = torch.tensor(elongation, dtype=torch.float64) X = base * scale.sqrt() * volume if rotate: Q, _ = torch.linalg.qr(torch.randn(d, d, generator=gen, dtype=torch.float64)) X = X @ Q.T return X torch.manual_seed(random_state) d = 10 centers = torch.randn(4, d, dtype=torch.float64) * 6 clusters = [ make_latent_cluster(150, d, seed=1, volume=1.0) + centers[0], make_latent_cluster(150, d, seed=2, volume=1.8) + centers[1], make_latent_cluster(150, d, seed=3, volume=1.0, elongation=[4, 4, 1, 1, 1, 1, 1, 1, 1, 1]) + centers[2], make_latent_cluster(150, d, seed=4, volume=1.0, elongation=[4, 4, 1, 1, 1, 1, 1, 1, 1, 1], rotate=True) + centers[3], ] X_latent = torch.cat(clusters, dim=0) print(X_latent.shape)

torch.Size([600, 10])

In [7]:

result_latent = select_model(
    X_latent, n_components=range(2, 6), criterion='bic',
    random_state=random_state, n_init=1, max_iter=100,
)

df_latent = result_latent.to_dataframe().sort_values('bic')
df_latent.head(10)

result_latent = select_model( X_latent, n_components=range(2, 6), criterion='bic', random_state=random_state, n_init=1, max_iter=100, ) df_latent = result_latent.to_dataframe().sort_values('bic') df_latent.head(10)

Out[7]:

	n_components	covariance_type	bic	aic	icl	log_likelihood	converged
33	4	VVI	21907.941776	21542.996615	21907.941776	-17.814164	True
31	4	VEI	21916.855605	21670.627544	21916.855605	-17.965523	True
45	5	VEI	21983.885007	21684.893791	22065.915103	-17.957411	True
29	4	VII	21986.613174	21779.957480	21986.613174	-18.071631	True
38	4	VVE	21990.115924	21427.308929	21990.115924	-17.642757	True
47	5	VVI	22013.273801	21555.993117	22080.064980	-17.789994	True
43	5	VII	22056.274831	21796.855982	22140.860137	-18.065713	True
52	5	VVE	22083.877563	21428.735045	22120.824224	-17.608946	True
35	4	VEE	22113.303741	21669.213846	22113.303741	-17.889345	True
49	5	VEE	22180.738906	21683.885855	22254.017093	-17.881572	True

BIC vs ICL for the winning K¶

Both criteria, compared across models at the BIC-winning K.

In [8]:

best_latent = result_latent.best_
winning_k = best_latent.n_components
print(f'Winning model: {best_latent.covariance_type_} (K={winning_k}), '
      f'BIC={best_latent.bic(X_latent):.2f}, ICL={best_latent.icl(X_latent):.2f}')

at_k = df_latent[df_latent['n_components'] == winning_k].dropna(subset=['bic', 'icl']).sort_values('bic')
at_k[['covariance_type', 'bic', 'icl', 'converged']].head(8)

best_latent = result_latent.best_ winning_k = best_latent.n_components print(f'Winning model: {best_latent.covariance_type_} (K={winning_k}), ' f'BIC={best_latent.bic(X_latent):.2f}, ICL={best_latent.icl(X_latent):.2f}') at_k = df_latent[df_latent['n_components'] == winning_k].dropna(subset=['bic', 'icl']).sort_values('bic') at_k[['covariance_type', 'bic', 'icl', 'converged']].head(8)

Winning model: VVI (K=4), BIC=21907.94, ICL=21907.94

Out[8]:

	covariance_type	bic	icl	converged
33	VVI	21907.941776	21907.941776	True
31	VEI	21916.855605	21916.855605	True
29	VII	21986.613174	21986.613174	True
38	VVE	21990.115924	21990.115924	True
35	VEE	22113.303741	22113.303741	True
30	EEI	22442.458830	22442.458830	True
32	EVI	22446.458565	22446.458565	True
28	EII	22474.234238	22474.234238	True

Takeaways¶

• select_model's default sweep (all 14 GPCM codes) finds a sensible winner without having to guess a covariance structure up front -- VVE-ish results on the iris data (shared orientation, elongated per-species shapes) are a well-known finding for this dataset.
• On genuinely mixed-structure data (axis-aligned vs. correlated clusters, Part 2), no single covariance type dominates at every K -- exactly the situation the 14-model family is for.
• volume_/shape_/orientation_ on the winning model expose why it won, not just that it did -- see the user guide for how to read them.