🔔 Gaussian Mixture Models: Probabilistic Clustering

📚 What You'll Learn

By the end of this lesson, you will be able to:

Explain how a Gaussian Mixture Model represents data as a weighted blend of Gaussian components
Describe how the Expectation-Maximization (EM) algorithm fits a GMM
Contrast GMM's soft, probabilistic assignments with K-means' hard cluster assignments
Choose an appropriate covariance type and select the number of components using AIC/BIC
Use a fitted GMM for density estimation and probabilistic cluster membership

⏱️ Estimated Time: 45–60 minutes

🎯 Project: Fit a Gaussian Mixture Model to a dataset, choose the number of components with BIC, and compare its soft clustering against K-means.

Introduction

Gaussian Mixture Models (GMM) are probabilistic models that assume data points are generated from a mixture of several Gaussian distributions with unknown parameters. Unlike K-means which performs hard clustering, GMM provides soft clustering where each point has a probability of belonging to each cluster. This makes GMM more flexible for handling clusters with different shapes, sizes, and densities.

GMM vs K-means Comparison

graph TD A[Clustering Task] --> B{Method} B --> C[K-means] C --> C1[Hard Assignment] C --> C2[Spherical Clusters] C --> C3[Equal Variance] C --> C4[Distance Based] B --> D[GMM] D --> D1[Soft Assignment] D --> D2[Elliptical Clusters] D --> D3[Different Variances] D --> D4[Probability Based] C1 --> E[Point belongs to ONE cluster] D1 --> F[Point has probability for EACH cluster] style C fill:#ffebee style D fill:#e8f5e9 style E fill:#fff3e0 style F fill:#e3f2fd

How GMM Works

flowchart LR A[Input Data] --> B[Initialize Parameters] B --> B1[Means μ] B --> B2[Covariances Σ] B --> B3[Weights π] B --> C[E-Step] C --> C1[Calculate Responsibilities] C1 --> C2[P(cluster|data point)] C --> D[M-Step] D --> D1[Update μ] D --> D2[Update Σ] D --> D3[Update π] D --> E{Converged?} E -->|No| C E -->|Yes| F[Final Model] F --> G[Soft Clusters] F --> H[Density Estimation] style A fill:#e3f2fd style F fill:#c8e6c9

Setting Up GMM

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.mixture import GaussianMixture
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import make_blobs, make_moons
from sklearn.model_selection import train_test_split
from scipy.stats import multivariate_normal
import warnings
warnings.filterwarnings('ignore')

# Set style for better visualizations
plt.style.use('seaborn-v0_8-darkgrid')
sns.set_palette("husl")

# Set random seed
np.random.seed(42)

print("Gaussian Mixture Models Tutorial")
print("=" * 50)

Basic GMM Implementation

class GMMDemo:
    """
    Demonstrate Gaussian Mixture Models
    """
    
    def __init__(self):
        self.model = None
        self.data = None
        
    def generate_mixed_data(self, n_samples=1000):
        """
        Generate data from multiple Gaussian distributions
        """
        # Create data with different cluster properties
        np.random.seed(42)
        
        # Cluster 1: Tight, circular
        cluster1 = np.random.multivariate_normal(
            mean=[0, 0],
            cov=[[0.5, 0], [0, 0.5]],
            size=n_samples // 3
        )
        
        # Cluster 2: Elongated, diagonal
        cluster2 = np.random.multivariate_normal(
            mean=[3, 3],
            cov=[[1, 0.7], [0.7, 1]],
            size=n_samples // 3
        )
        
        # Cluster 3: Wide, horizontal
        cluster3 = np.random.multivariate_normal(
            mean=[-2, 3],
            cov=[[2, 0], [0, 0.3]],
            size=n_samples // 3
        )
        
        # Combine clusters
        X = np.vstack([cluster1, cluster2, cluster3])
        
        # Create labels for visualization
        y_true = np.hstack([
            np.zeros(len(cluster1)),
            np.ones(len(cluster2)),
            np.ones(len(cluster3)) * 2
        ])
        
        # Shuffle data
        indices = np.random.permutation(len(X))
        X = X[indices]
        y_true = y_true[indices]
        
        return X, y_true
    
    def fit_gmm(self, X, n_components=3, covariance_type='full'):
        """
        Fit Gaussian Mixture Model
        
        covariance_type options:
        - 'full': each component has its own covariance matrix
        - 'tied': all components share the same covariance matrix
        - 'diag': diagonal covariance (axes aligned)
        - 'spherical': single variance per component (like K-means)
        """
        self.model = GaussianMixture(
            n_components=n_components,
            covariance_type=covariance_type,
            max_iter=100,
            random_state=42
        )
        
        self.model.fit(X)
        
        # Get cluster assignments
        labels = self.model.predict(X)
        
        # Get soft assignments (probabilities)
        proba = self.model.predict_proba(X)
        
        return labels, proba
    
    def visualize_results(self, X, labels, proba, y_true=None):
        """
        Visualize GMM clustering results
        """
        fig, axes = plt.subplots(2, 2, figsize=(12, 10))
        
        # Plot 1: True clusters (if available)
        if y_true is not None:
            axes[0, 0].scatter(X[:, 0], X[:, 1], c=y_true, cmap='viridis',
                             edgecolor='black', linewidth=0.5, s=30)
            axes[0, 0].set_title('True Clusters')
        else:
            axes[0, 0].scatter(X[:, 0], X[:, 1], c='gray',
                             edgecolor='black', linewidth=0.5, s=30)
            axes[0, 0].set_title('Input Data')
        
        # Plot 2: GMM hard clustering
        axes[0, 1].scatter(X[:, 0], X[:, 1], c=labels, cmap='viridis',
                         edgecolor='black', linewidth=0.5, s=30)
        axes[0, 1].set_title('GMM Hard Clustering')
        
        # Plot 3: Probability of dominant cluster
        max_proba = np.max(proba, axis=1)
        scatter = axes[1, 0].scatter(X[:, 0], X[:, 1], c=max_proba, 
                                    cmap='YlOrRd', vmin=0, vmax=1,
                                    edgecolor='black', linewidth=0.5, s=30)
        axes[1, 0].set_title('Cluster Confidence\n(Max Probability)')
        plt.colorbar(scatter, ax=axes[1, 0])
        
        # Plot 4: Gaussian ellipses
        axes[1, 1].scatter(X[:, 0], X[:, 1], c=labels, cmap='viridis',
                         alpha=0.3, s=20)
        
        # Draw ellipses for each component
        self.draw_ellipses(axes[1, 1], self.model)
        axes[1, 1].set_title('Gaussian Components')
        
        # Set labels
        for ax in axes.flat:
            ax.set_xlabel('Feature 1')
            ax.set_ylabel('Feature 2')
            ax.grid(True, alpha=0.3)
        
        plt.suptitle('Gaussian Mixture Model Results', fontsize=14)
        plt.tight_layout()
        plt.show()
    
    def draw_ellipses(self, ax, gmm):
        """
        Draw ellipses representing Gaussian components
        """
        from matplotlib.patches import Ellipse
        
        colors = plt.cm.viridis(np.linspace(0, 0.8, gmm.n_components))
        
        for i in range(gmm.n_components):
            # Get mean and covariance
            mean = gmm.means_[i]
            
            if gmm.covariance_type == 'full':
                cov = gmm.covariances_[i]
            elif gmm.covariance_type == 'tied':
                cov = gmm.covariances_
            elif gmm.covariance_type == 'diag':
                cov = np.diag(gmm.covariances_[i])
            elif gmm.covariance_type == 'spherical':
                cov = np.eye(2) * gmm.covariances_[i]
            
            # Compute eigenvalues and eigenvectors
            eigenvalues, eigenvectors = np.linalg.eig(cov)
            angle = np.degrees(np.arctan2(eigenvectors[1, 0], eigenvectors[0, 0]))
            
            # Draw ellipse at 2 standard deviations
            width = 2 * np.sqrt(eigenvalues[0]) * 2  # 2 std
            height = 2 * np.sqrt(eigenvalues[1]) * 2
            
            ellipse = Ellipse(mean, width, height, angle=angle,
                             facecolor=colors[i], alpha=0.3,
                             edgecolor=colors[i], linewidth=2)
            ax.add_patch(ellipse)
            
            # Mark center
            ax.scatter(mean[0], mean[1], c='red', s=100, marker='x')

# Create and demonstrate GMM
gmm_demo = GMMDemo()

# Generate data
print("Generating mixed Gaussian data...")
X, y_true = gmm_demo.generate_mixed_data(n_samples=600)

# Fit GMM
print("Fitting Gaussian Mixture Model...")
labels, proba = gmm_demo.fit_gmm(X, n_components=3)

# Visualize results
print("Visualizing results...")
gmm_demo.visualize_results(X, labels, proba, y_true)

# Print model parameters
print("\nModel Parameters:")
print(f"Means:\n{gmm_demo.model.means_}")
print(f"\nWeights: {gmm_demo.model.weights_}")
print(f"\nConverged: {gmm_demo.model.converged_}")

Covariance Types Comparison

def compare_covariance_types(X):
    """
    Compare different covariance types in GMM
    """
    covariance_types = ['spherical', 'diag', 'tied', 'full']
    n_components = 3
    
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    axes = axes.flatten()
    
    for idx, cov_type in enumerate(covariance_types):
        # Fit GMM
        gmm = GaussianMixture(
            n_components=n_components,
            covariance_type=cov_type,
            random_state=42
        )
        gmm.fit(X)
        
        # Predict
        labels = gmm.predict(X)
        
        # Plot
        axes[idx].scatter(X[:, 0], X[:, 1], c=labels, cmap='viridis',
                         alpha=0.5, s=30)
        
        # Draw ellipses
        gmm_demo.draw_ellipses(axes[idx], gmm)
        
        # Calculate BIC for model selection
        bic = gmm.bic(X)
        aic = gmm.aic(X)
        
        axes[idx].set_title(f'{cov_type.capitalize()}\nBIC: {bic:.0f}, AIC: {aic:.0f}')
        axes[idx].set_xlabel('Feature 1')
        axes[idx].set_ylabel('Feature 2')
        axes[idx].grid(True, alpha=0.3)
    
    plt.suptitle('GMM Covariance Types Comparison', fontsize=14)
    plt.tight_layout()
    plt.show()
    
    # Print comparison
    print("\nCovariance Types Explanation:")
    print("-" * 50)
    print("Spherical: Each component has single variance (like K-means)")
    print("Diagonal: Each component has diagonal covariance (axis-aligned)")
    print("Tied: All components share the same covariance matrix")
    print("Full: Each component has its own full covariance matrix")
    print("\nLower BIC/AIC indicates better model fit")

# Compare covariance types
print("\n" + "=" * 50)
print("COMPARING COVARIANCE TYPES")
print("=" * 50)
compare_covariance_types(X)

Model Selection: Finding Optimal Components

flowchart TB A[Dataset] --> B[Try Different K] B --> C[K=2] B --> D[K=3] B --> E[K=4] B --> F[K=5] C --> G[Calculate BIC/AIC] D --> G E --> G F --> G G --> H{Select Best K} H --> I[Lowest BIC] H --> J[Elbow Method] H --> K[Domain Knowledge] I --> L[Final Model] style A fill:#e3f2fd style L fill:#c8e6c9

def find_optimal_components(X, max_components=10):
    """
    Find optimal number of components using BIC and AIC
    """
    n_components_range = range(1, max_components + 1)
    bic_scores = []
    aic_scores = []
    
    for n in n_components_range:
        gmm = GaussianMixture(
            n_components=n,
            covariance_type='full',
            random_state=42
        )
        gmm.fit(X)
        bic_scores.append(gmm.bic(X))
        aic_scores.append(gmm.aic(X))
    
    # Plot scores
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
    
    # BIC scores
    ax1.plot(n_components_range, bic_scores, 'bo-', linewidth=2, markersize=8)
    ax1.set_xlabel('Number of Components')
    ax1.set_ylabel('BIC Score')
    ax1.set_title('BIC Score vs Number of Components')
    ax1.grid(True, alpha=0.3)
    
    # Mark minimum
    min_bic_idx = np.argmin(bic_scores)
    ax1.axvline(x=n_components_range[min_bic_idx], color='red', 
                linestyle='--', label=f'Best: {n_components_range[min_bic_idx]}')
    ax1.legend()
    
    # AIC scores
    ax2.plot(n_components_range, aic_scores, 'go-', linewidth=2, markersize=8)
    ax2.set_xlabel('Number of Components')
    ax2.set_ylabel('AIC Score')
    ax2.set_title('AIC Score vs Number of Components')
    ax2.grid(True, alpha=0.3)
    
    # Mark minimum
    min_aic_idx = np.argmin(aic_scores)
    ax2.axvline(x=n_components_range[min_aic_idx], color='red',
                linestyle='--', label=f'Best: {n_components_range[min_aic_idx]}')
    ax2.legend()
    
    plt.suptitle('Model Selection: Optimal Number of Components', fontsize=14)
    plt.tight_layout()
    plt.show()
    
    print(f"Optimal components by BIC: {n_components_range[min_bic_idx]}")
    print(f"Optimal components by AIC: {n_components_range[min_aic_idx]}")
    
    return n_components_range[min_bic_idx]

# Find optimal number of components
print("\n" + "=" * 50)
print("MODEL SELECTION")
print("=" * 50)
optimal_k = find_optimal_components(X, max_components=8)

GMM for Density Estimation

class DensityEstimation:
    """
    Use GMM for probability density estimation
    """
    
    def __init__(self):
        self.gmm = None
        
    def fit_density_model(self, X, n_components):
        """
        Fit GMM for density estimation
        """
        self.gmm = GaussianMixture(
            n_components=n_components,
            covariance_type='full',
            random_state=42
        )
        self.gmm.fit(X)
        
    def plot_density_contours(self, X):
        """
        Plot density contours from GMM
        """
        # Create mesh
        h = 0.1
        x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
        y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
        xx, yy = np.meshgrid(np.arange(x_min, x_max, h),
                            np.arange(y_min, y_max, h))
        
        # Calculate log probability for each point in mesh
        mesh_points = np.c_[xx.ravel(), yy.ravel()]
        log_prob = self.gmm.score_samples(mesh_points)
        log_prob = log_prob.reshape(xx.shape)
        
        # Plot
        fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
        
        # Density contours
        contour = ax1.contourf(xx, yy, np.exp(log_prob), levels=20, cmap='YlOrRd')
        ax1.scatter(X[:, 0], X[:, 1], c='black', s=10, alpha=0.5)
        ax1.set_title('Probability Density Estimation')
        ax1.set_xlabel('Feature 1')
        ax1.set_ylabel('Feature 2')
        plt.colorbar(contour, ax=ax1, label='Density')
        
        # Log probability
        contour2 = ax2.contour(xx, yy, log_prob, levels=15, colors='black', alpha=0.4)
        ax2.scatter(X[:, 0], X[:, 1], c=self.gmm.predict(X), 
                   cmap='viridis', s=30, alpha=0.7)
        ax2.set_title('Log Probability Contours')
        ax2.set_xlabel('Feature 1')
        ax2.set_ylabel('Feature 2')
        ax2.clabel(contour2, inline=True, fontsize=8)
        
        plt.suptitle('GMM for Density Estimation', fontsize=14)
        plt.tight_layout()
        plt.show()
    
    def detect_outliers(self, X, threshold_percentile=5):
        """
        Detect outliers using GMM density estimation
        """
        # Get log probabilities
        log_prob = self.gmm.score_samples(X)
        
        # Set threshold
        threshold = np.percentile(log_prob, threshold_percentile)
        
        # Identify outliers
        outliers = log_prob < threshold
        
        # Visualize
        plt.figure(figsize=(10, 6))
        plt.scatter(X[~outliers, 0], X[~outliers, 1], 
                   c='blue', label='Normal', alpha=0.6, s=30)
        plt.scatter(X[outliers, 0], X[outliers, 1],
                   c='red', label='Outliers', s=50, marker='x')
        plt.xlabel('Feature 1')
        plt.ylabel('Feature 2')
        plt.title(f'Outlier Detection using GMM\n({np.sum(outliers)} outliers detected)')
        plt.legend()
        plt.grid(True, alpha=0.3)
        plt.show()
        
        return outliers

# Demonstrate density estimation
print("\n" + "=" * 50)
print("DENSITY ESTIMATION WITH GMM")
print("=" * 50)

density_estimator = DensityEstimation()
density_estimator.fit_density_model(X, n_components=3)

print("Plotting density contours...")
density_estimator.plot_density_contours(X)

print("\nDetecting outliers...")
outliers = density_estimator.detect_outliers(X, threshold_percentile=5)
print(f"Found {np.sum(outliers)} outliers ({np.sum(outliers)/len(X)*100:.1f}% of data)")

GMM vs K-means Practical Comparison

from sklearn.cluster import KMeans

def compare_gmm_kmeans():
    """
    Compare GMM and K-means on different datasets
    """
    # Generate different types of data
    datasets = {}
    
    # Dataset 1: Well-separated spherical clusters
    X1, y1 = make_blobs(n_samples=300, centers=3, n_features=2,
                        cluster_std=0.5, random_state=42)
    datasets['Well-separated'] = (X1, y1)
    
    # Dataset 2: Overlapping clusters
    X2, y2 = make_blobs(n_samples=300, centers=3, n_features=2,
                        cluster_std=1.5, random_state=42)
    datasets['Overlapping'] = (X2, y2)
    
    # Dataset 3: Different variances
    X3 = np.vstack([
        np.random.multivariate_normal([0, 0], [[0.1, 0], [0, 0.1]], 100),
        np.random.multivariate_normal([2, 2], [[1, 0], [0, 1]], 100),
        np.random.multivariate_normal([-2, 2], [[0.5, 0.3], [0.3, 0.5]], 100)
    ])
    y3 = np.hstack([np.zeros(100), np.ones(100), np.ones(100)*2])
    datasets['Different Variances'] = (X3, y3)
    
    # Compare on each dataset
    fig, axes = plt.subplots(3, 3, figsize=(15, 12))
    
    for row, (name, (X_data, y_data)) in enumerate(datasets.items()):
        # Standardize data
        scaler = StandardScaler()
        X_scaled = scaler.fit_transform(X_data)
        
        # True clusters
        axes[row, 0].scatter(X_scaled[:, 0], X_scaled[:, 1], 
                           c=y_data, cmap='viridis', s=30, alpha=0.7)
        axes[row, 0].set_title(f'{name}\nTrue Clusters')
        
        # K-means
        kmeans = KMeans(n_clusters=3, random_state=42)
        kmeans_labels = kmeans.fit_predict(X_scaled)
        axes[row, 1].scatter(X_scaled[:, 0], X_scaled[:, 1],
                           c=kmeans_labels, cmap='viridis', s=30, alpha=0.7)
        axes[row, 1].set_title('K-means')
        
        # GMM
        gmm = GaussianMixture(n_components=3, covariance_type='full', random_state=42)
        gmm_labels = gmm.fit_predict(X_scaled)
        axes[row, 2].scatter(X_scaled[:, 0], X_scaled[:, 1],
                           c=gmm_labels, cmap='viridis', s=30, alpha=0.7)
        axes[row, 2].set_title('GMM')
    
    # Set labels
    for ax in axes.flat:
        ax.set_xlabel('Feature 1')
        ax.set_ylabel('Feature 2')
        ax.grid(True, alpha=0.3)
    
    plt.suptitle('GMM vs K-means Comparison', fontsize=14)
    plt.tight_layout()
    plt.show()

# Run comparison
print("\n" + "=" * 50)
print("GMM VS K-MEANS COMPARISON")
print("=" * 50)
compare_gmm_kmeans()

Best Practices

🎯 GMM Guidelines

Feature Scaling: Standardize features before fitting GMM
Initialization: Use k-means initialization for stability
Covariance Type: Start with 'full', try simpler types for speed
Model Selection: Use BIC/AIC for choosing number of components
Convergence: Check convergence and increase iterations if needed
Singularities: Add regularization if covariance becomes singular
Interpretability: Examine component parameters for insights
Outlier Detection: Use low probability threshold for anomalies

Practice Exercises

Exercise 1: Customer Segmentation with GMM

Build a customer segmentation system using GMM:

Load customer purchase data
Select optimal number of segments using BIC
Compare different covariance types
Interpret segment characteristics
Calculate customer probability for each segment

Exercise 2: Image Segmentation

Use GMM for image color segmentation:

Load an image and extract RGB values
Fit GMM to color space
Segment image into regions
Compare with K-means segmentation

Summary

✅ You've Learned

GMM fundamentals and probabilistic clustering
Difference between hard and soft clustering
EM algorithm for parameter estimation
Different covariance types and their effects
Model selection using BIC and AIC
Using GMM for density estimation and outlier detection
When to use GMM vs K-means

📓 Learning Journal

Keep a learning journal — digital or physical. After this lesson, take a few minutes to write down:

Key concepts you learned
Techniques that clicked for you
Questions or confusion points to revisit
Ideas you want to try
Your progress and feelings about learning this

✍️ This lesson's prompt: GMM gives each point a probability of belonging to each cluster rather than a single hard label. When has a "soft" answer ("70% this, 30% that") described reality better than forcing a single category — and why might that matter for clustering?

📝 Lesson Summary

🎓 Key Takeaways

A Gaussian Mixture Model represents data as a weighted sum of Gaussian components, each with its own mean and covariance.
It is fit with the Expectation-Maximization (EM) algorithm, which alternates soft assignments and parameter updates.
Unlike K-means' hard assignments, GMM gives each point a probability of membership in every cluster, and its covariances allow elliptical clusters.
Choose the number of components with information criteria (AIC/BIC) rather than guessing.

🎉 What You've Accomplished

You can now fit a Gaussian Mixture Model, interpret its soft probabilistic assignments, select the right number of components and covariance type, and use the fitted model for both clustering and density estimation.

❓ Common Questions at This Stage

How is GMM different from K-means?

K-means makes hard assignments and effectively assumes spherical, equal-size clusters. GMM makes soft, probabilistic assignments and models each cluster's shape via a covariance matrix, so it handles elliptical, overlapping clusters.

How many components should I use?

Fit GMMs across a range of component counts and compare AIC or BIC; the value that minimizes the criterion balances fit against complexity. Domain knowledge and validation on downstream tasks also help.

What does the covariance type control?

It sets how flexible each cluster's shape can be — from spherical (circular) through diag and tied to full (arbitrary ellipses). More flexibility fits complex data but needs more samples and risks overfitting.

🔭 Looking Ahead

GMM rounds out your clustering toolkit alongside K-means and DBSCAN. With hard, density-based, and probabilistic clustering in hand, you're ready to apply unsupervised learning to real datasets and richer applications.

✅ Before the Next Lesson

Fit GMMs with 1–10 components and plot AIC/BIC to pick the best number of components.
Compare GMM's soft cluster probabilities against K-means' hard labels on the same overlapping data.
Write your Learning Journal entry for this lesson

🌟 Encouragement for the Journey

Probabilistic thinking is a superpower in machine learning, and GMM is a beautiful example of it. You've added a flexible, principled clustering method to your toolkit — keep applying it and trust your growing intuition.