Gaussian Mixture Models und Expectation-Maximization

Probabilistisches Clustering mit scikit-learn

Ein Gaussian Mixture Model (GMM) beschreibt eine Wahrscheinlichkeitsverteilung als gewichtete Mischung mehrerer Gauß-Verteilungen. GMMs werden häufig für unüberwachtes Lernen, Clustering und Dichteschätzung eingesetzt.

Anders als ein rein hartes Clustering kann ein GMM für jeden Datenpunkt angeben, mit welcher Wahrscheinlichkeit er zu den einzelnen Komponenten gehört. Diese Wahrscheinlichkeiten werden häufig Responsibilities genannt.

In scikit-learn wird ein klassisches GMM durch sklearn.mixture.GaussianMixture bereitgestellt.

Dokumentation: GaussianMixture

GMM und k-means: wichtige Unterschiede

GMM wird gelegentlich mit k-means verwechselt. Beide Verfahren können Clusterstrukturen entdecken, verfolgen aber unterschiedliche Modelle.

k-means

  • ordnet jeden Punkt genau einem Cluster zu (hard assignment),
  • minimiert quadrierte euklidische Abstände zu Zentren,
  • führt dadurch zu Voronoi-Zellen um die Clusterzentren.

Gaussian Mixture Model

  • modelliert eine vollständige Wahrscheinlichkeitsdichte,
  • liefert für jeden Punkt Wahrscheinlichkeiten für alle Komponenten (soft assignment),
  • kann mit covariance_type="full" unterschiedlich orientierte ellipsoidale Komponenten beschreiben,
  • berücksichtigt neben Lage und Kovarianz auch das Mischungsgewicht jeder Komponente.

Ein GMM ist daher kein „kNN mit Ellipsen“. k-nearest neighbors (kNN) ist ein anderes, überwiegend überwachtes Verfahren und gehört konzeptionell nicht in diesen Vergleich.

Mathematisches Modell

Für (K) Komponenten lautet die Mischungsverteilung

[ p(\mathbf{x}) = \sum_{k=1}^{K} \pi_k,\mathcal{N}(\mathbf{x}\mid \boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k), ]

wobei (\pi_k \ge 0), (\sum_k \pi_k = 1), (\boldsymbol{\mu}_k) den Mittelwertvektor und (\boldsymbol{\Sigma}_k) die Kovarianzmatrix bezeichnet.

Die Zugehörigkeitswahrscheinlichkeit (Responsibility) eines Datenpunkts (\mathbf{x}_i) zur Komponente (k) ist

[ \gamma_{ik} = \frac{\pi_k,\mathcal{N}(\mathbf{x}_i\mid\boldsymbol{\mu}_k,\boldsymbol{\Sigma}k)} {\sum{j=1}^{K}\pi_j,\mathcal{N}(\mathbf{x}_i\mid\boldsymbol{\mu}_j,\boldsymbol{\Sigma}_j)}. ]

Für jeden Datenpunkt gilt (\sum_k \gamma_{ik}=1).

Expectation-Maximization (EM)

Die Parameter eines GMM werden typischerweise iterativ mit dem Expectation-Maximization-Algorithmus geschätzt.

E-Schritt

Mit den aktuellen Parametern werden die Responsibilities (\gamma_{ik}) berechnet.

M-Schritt

Mit (N_k = \sum_i \gamma_{ik}) werden die Parameter neu geschätzt:

[ \pi_k = \frac{N_k}{N}, \qquad \boldsymbol{\mu}k = \frac{1}{N_k}\sum_i \gamma{ik}\mathbf{x}_i, ]

[ \boldsymbol{\Sigma}k = \frac{1}{N_k}\sum_i \gamma{ik} (\mathbf{x}_i-\boldsymbol{\mu}_k) (\mathbf{x}_i-\boldsymbol{\mu}_k)^T. ]

Danach beginnt erneut der E-Schritt. EM erhöht die Log-Likelihood bei jedem regulären Iterationsschritt beziehungsweise lässt sie unverändert, kann aber in einem lokalen Optimum enden. Deshalb sind Initialisierung und gegebenenfalls mehrere Initialisierungen (n_init) wichtig.

EM in einer Dimension selbst implementiert

Zunächst implementieren wir den Kern des Verfahrens für drei eindimensionale Gauß-Komponenten. Das Beispiel ist bewusst kompakt gehalten: Es soll die E- und M-Schritte sichtbar machen, nicht GaussianMixture nachbauen.

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)

x = np.concatenate([
    rng.normal(-4.0, 0.8, 180),
    rng.normal(0.5, 1.1, 260),
    rng.normal(5.0, 0.6, 160),
])
rng.shuffle(x)

def normal_pdf(x, mean, variance):
    return np.exp(-0.5 * (x - mean) ** 2 / variance) / np.sqrt(
        2 * np.pi * variance
    )

K = 3
weights = np.full(K, 1 / K)
means = np.array([-5.5, -0.5, 4.0], dtype=float)
variances = np.full(K, 2.0)

log_likelihood = []

for _ in range(40):
    weighted_densities = np.column_stack([
        weights[k] * normal_pdf(x, means[k], variances[k])
        for k in range(K)
    ])

    total_density = weighted_densities.sum(axis=1, keepdims=True)
    responsibilities = weighted_densities / total_density
    log_likelihood.append(np.log(total_density[:, 0]).sum())

    N_k = responsibilities.sum(axis=0)
    weights = N_k / len(x)
    means = (responsibilities * x[:, None]).sum(axis=0) / N_k
    variances = (
        responsibilities * (x[:, None] - means) ** 2
    ).sum(axis=0) / N_k

order = np.argsort(means)

print("Gewichte:", np.round(weights[order], 3))
print("Mittelwerte:", np.round(means[order], 3))
print("Standardabweichungen:", np.round(np.sqrt(variances[order]), 3))
print("Zeilensummen der Responsibilities:",
      np.round(responsibilities[:5].sum(axis=1), 6))

Gewichte: [0.295 0.44  0.265]
Mittelwerte: [-4.076  0.472  4.988]
Standardabweichungen: [0.655 1.164 0.607]
Zeilensummen der Responsibilities: [1. 1. 1. 1. 1.]

Die Zeilensummen der Responsibility-Matrix sind 1. Jeder Datenpunkt verteilt seine Zugehörigkeit also probabilistisch über alle Komponenten.

Die geschätzten Mittelwerte und Standardabweichungen sollten nahe bei den Parametern liegen, mit denen die Daten erzeugt wurden. Die Reihenfolge der Komponenten ist grundsätzlich bedeutungslos: „Komponente 0“ besitzt keine fest vorgegebene inhaltliche Bedeutung.

grid = np.linspace(x.min() - 2, x.max() + 2, 600)

plt.figure(figsize=(9, 5))
plt.hist(x, bins=45, density=True, alpha=0.35)

mixture_density = np.zeros_like(grid)

for k in order:
    component = weights[k] * normal_pdf(grid, means[k], variances[k])
    mixture_density += component
    plt.plot(grid, component, linewidth=2)

plt.plot(grid, mixture_density, linewidth=3)
plt.xlabel("x")
plt.ylabel("Dichte")
plt.title("Geschätzte Komponenten und Gesamtmischung")
plt.show()
No description has been provided for this image 
plt.figure(figsize=(8, 4))
plt.plot(log_likelihood)
plt.xlabel("EM-Iteration")
plt.ylabel("Log-Likelihood")
plt.title("Entwicklung der Log-Likelihood")
plt.show()

diffs = np.diff(log_likelihood)
print("Kleinste Änderung:", diffs.min())
No description has been provided for this image 
Kleinste Änderung: -4.547473508864641e-13

Numerisch können durch Rundung winzige Abweichungen auftreten. Im idealen EM-Ablauf sinkt die Log-Likelihood jedoch nicht. Das Verfahren stoppt, wenn sich die Zielfunktion beziehungsweise die Parameter nur noch vernachlässigbar ändern.

Zweidimensionales GMM mit scikit-learn

Nun verwenden wir GaussianMixture auf einem zweidimensionalen Datensatz. Die Daten werden linear transformiert, damit die Cluster nicht nur kreisförmig, sondern anisotrop und unterschiedlich orientiert erscheinen.

from sklearn.datasets import make_blobs
from sklearn.mixture import GaussianMixture

X, _ = make_blobs(
    n_samples=700,
    centers=3,
    cluster_std=(0.8, 1.0, 0.7),
    random_state=12,
)

transformation = np.array([
    [0.7, -0.6],
    [0.35, 1.1],
])
X = X @ transformation

gmm_model = GaussianMixture(
    n_components=3,
    covariance_type="full",
    n_init=5,
    reg_covar=1e-6,
    random_state=12,
)

gmm_model.fit(X)

print("Konvergiert:", gmm_model.converged_)
print("Iterationen:", gmm_model.n_iter_)
print("Gewichte:", np.round(gmm_model.weights_, 3))
print("Mittelwerte:")
print(np.round(gmm_model.means_, 3))

Konvergiert: True
Iterationen: 3
Gewichte: [0.336 0.331 0.333]
Mittelwerte:
[[-3.226  9.379]
 [-2.991  3.574]
 [-3.923 15.016]]

GaussianMixture verwendet EM zur Maximum-Likelihood-Schätzung. n_init=5 startet die Parameterschätzung mehrfach und behält die beste Lösung. Das ist hilfreich, weil EM von der Initialisierung abhängt.

reg_covar addiert einen kleinen nichtnegativen Wert auf die Diagonale der Kovarianzmatrizen. Dadurch werden die Kovarianzen numerisch stabiler und positiv definit gehalten.

probabilities = gmm_model.predict_proba(X[:5])
hard_labels = gmm_model.predict(X[:5])

print("Responsibilities der ersten fünf Punkte:")
print(np.round(probabilities, 3))
print("\nZeilensummen:", np.round(probabilities.sum(axis=1), 6))
print("Harte Zuordnung:", hard_labels)

Responsibilities der ersten fünf Punkte:
[[0. 0. 1.]
 [0. 0. 1.]
 [0. 1. 0.]
 [1. 0. 0.]
 [0. 1. 0.]]

Zeilensummen: [1. 1. 1. 1. 1.]
Harte Zuordnung: [2 2 1 0 1]

predict_proba liefert die Responsibilities. predict reduziert diese Information auf die jeweils wahrscheinlichste Komponente.

Die numerischen Clusterlabels sind beliebig. Bei unüberwachtem Lernen darf man nicht erwarten, dass die Nummern automatisch mit eventuell vorhandenen Klassenlabels übereinstimmen.

from matplotlib.patches import Ellipse

def add_covariance_ellipse(ax, mean, covariance, n_std=2.0):
    eigenvalues, eigenvectors = np.linalg.eigh(covariance)
    order = np.argsort(eigenvalues)[::-1]
    eigenvalues = eigenvalues[order]
    eigenvectors = eigenvectors[:, order]

    angle = np.degrees(
        np.arctan2(eigenvectors[1, 0], eigenvectors[0, 0])
    )
    width, height = 2 * n_std * np.sqrt(eigenvalues)

    ellipse = Ellipse(
        xy=mean,
        width=width,
        height=height,
        angle=angle,
        fill=False,
        linewidth=2,
    )
    ax.add_patch(ellipse)

labels = gmm_model.predict(X)

fig, ax = plt.subplots(figsize=(8, 6))
ax.scatter(X[:, 0], X[:, 1], c=labels, s=20, alpha=0.65)

for mean, covariance in zip(gmm_model.means_, gmm_model.covariances_):
    add_covariance_ellipse(ax, mean, covariance)
    ax.scatter(mean[0], mean[1], marker="x", s=100)

ax.set_title("Gaussian Mixture Model mit vollen Kovarianzmatrizen")
ax.set_xlabel("$x_1$")
ax.set_ylabel("$x_2$")
plt.show()
No description has been provided for this image

Wie viele Komponenten?

Die Anzahl der Komponenten ist nicht immer bekannt. Die reine Trainings-Log-Likelihood ist dafür ungeeignet, weil ein komplexeres Modell mit mehr Komponenten die Daten meist mindestens ebenso gut anpassen kann.

Für GMMs können beispielsweise AIC und BIC verwendet werden. Beide bestrafen zusätzliche Modellkomplexität. Bei BIC gilt: kleiner ist besser.

Das folgende Beispiel vergleicht Modelle mit einer bis sechs Komponenten.

component_counts = range(1, 7)
bic_values = []
models = []

for n_components in component_counts:
    model = GaussianMixture(
        n_components=n_components,
        covariance_type="full",
        n_init=5,
        reg_covar=1e-6,
        random_state=12,
    ).fit(X)
    models.append(model)
    bic_values.append(model.bic(X))

best_index = int(np.argmin(bic_values))
best_components = list(component_counts)[best_index]

print("Nach BIC bevorzugte Komponentenanzahl:", best_components)

plt.figure(figsize=(8, 4))
plt.plot(list(component_counts), bic_values, marker="o")
plt.xlabel("Anzahl Komponenten")
plt.ylabel("BIC")
plt.title("GMM-Modellwahl mit BIC")
plt.show()

Nach BIC bevorzugte Komponentenanzahl: 3
No description has been provided for this image

BIC ist ein Auswahlkriterium, kein Beweis dafür, dass eine bestimmte „wahre“ Clusterzahl existiert. Ob die statistischen Komponenten fachlich sinnvolle Gruppen darstellen, muss immer im Kontext der Daten beurteilt werden.

Scikit-learn zeigt in seinem offiziellen Beispiel ebenfalls die Auswahl von Komponentenanzahl und Kovarianztyp mit Informationskriterien:

Gaussian Mixture Model Selection

Kovarianztypen

GaussianMixture unterstützt mehrere Annahmen für die Kovarianzmatrizen:

  • "full": jede Komponente besitzt eine eigene vollständige Kovarianzmatrix,
  • "tied": alle Komponenten teilen sich eine vollständige Kovarianzmatrix,
  • "diag": jede Komponente besitzt nur diagonale Kovarianzen,
  • "spherical": pro Komponente wird nur eine Varianz verwendet.

"full" ist am flexibelsten, benötigt aber auch die meisten Parameter. Bei wenig Daten oder hoher Dimensionalität können einfachere Kovarianzmodelle stabiler sein.

Numerische Singularitäten und Regularisierung

Ein klassisches Problem von Gauß-Mischungen entsteht, wenn eine Komponente auf sehr wenige Punkte oder praktisch auf einen einzelnen Punkt kollabiert. Dann kann ihre Kovarianzmatrix nahezu singulär werden und die Likelihood ohne geeignete Regularisierung degenerieren.

Scikit-learn begegnet dem unter anderem mit reg_covar. Der Parameter fügt einen kleinen Wert zur Diagonale der Kovarianzmatrix hinzu. Der Standardwert ist in der aktuellen API 1e-6.

Mehrere Initialisierungen (n_init) und eine sachgerechte Wahl der Komponentenanzahl sind ebenfalls wichtig. Eine extrem große Anzahl von Komponenten ist nicht automatisch ein besseres Modell.

BayesianGaussianMixture

Wenn die Zahl der relevanten Komponenten nicht fest vorgegeben werden soll, bietet scikit-learn außerdem BayesianGaussianMixture. Dieses Modell verwendet variationale Bayes-Schätzung und kann bei ausreichend vielen initial erlaubten Komponenten unwichtige Komponenten stark heruntergewichten.

Das ist konzeptionell ein anderes Schätzverfahren als das klassische Maximum-Likelihood-GMM mit EM und sollte deshalb nicht einfach als austauschbarer Parameter betrachtet werden.

Dokumentation: BayesianGaussianMixture

Zusammenfassung

Gaussian Mixture Models sind besonders nützlich, wenn

  • Cluster probabilistisch statt nur hart zugeordnet werden sollen,
  • Komponenten unterschiedliche Streuungen und Korrelationen besitzen können,
  • eine Dichteschätzung gewünscht ist.

Die zentralen Punkte sind:

  1. GMM modelliert eine Mischung von Wahrscheinlichkeitsverteilungen.
  2. EM wechselt zwischen Responsibility-Berechnung und Parameterschätzung.
  3. Initialisierung ist wichtig, weil EM lokale Optima finden kann.
  4. predict_proba liefert die weichen Clusterzugehörigkeiten.
  5. AIC/BIC können bei der Modellwahl helfen.
  6. reg_covar verbessert die numerische Stabilität der Kovarianzmatrizen.

Damit ist GMM deutlich mehr als nur ein alternatives Verfahren zum Zeichnen von Clustergrenzen: Es ist ein probabilistisches generatives Modell der Datenverteilung.