[Python] Clustering mit einem unendlich gemischten Gaußschen Modell

Eine der Clustering-Methoden ist die Schätzung des gemischten Gaußschen Modells. Dies ist eine Methode zum Ausdrücken und Clustering eines Datensatzes mit mehreren Gaußschen Verteilungen. Normalerweise nehmen diese Gaußschen Verteilungen eine endliche Zahl an, aber wenn es eine unendliche Zahl gäbe, welche Art von Clusterbildung wäre das?

In diesem Artikel wird auf Maschinelles Lernen Professional Series "Nichtparametrischer Bayes-Punkt-Prozess und statistische Mathematik des maschinellen Lernens" verwiesen. Einführung in die Gliederung und Implementierung von "unendlich gemischtem Gaußschen Modell" und "nicht parametrisch". Ich empfehle dieses Buch, weil es sehr leicht zu verstehen ist.

Ich möchte die detaillierte Theorie weglassen und mich auf den Umriss und die Geschichte konzentrieren. Außerdem sollte die Implementierung auf den Quellcode beschränkt sein.

• Vorheriger Artikel [Python] Ich habe marginalisiertes Gibbs-Sampling implementiert Ich bin mir jedoch bewusst, dass der Fluss vermittelt werden kann, auch wenn ich nur diesen Artikel lese.

Einführung eines unendlichen gemischten Gaußschen Modells

In einem normalen gemischten Gaußschen Modell wird der Datensatz durch eine endliche Anzahl von Gaußschen Verteilungen dargestellt, und jedes Daten wird als mit der darin enthaltenen Gaußschen Verteilung assoziiert angesehen. Die Anzahl der Clasper (= Anzahl der Gaußschen Verteilungen) muss im Voraus bestimmt werden, aber die entsprechende Anzahl von Clustern variiert je nach Datensatz und ist nicht leicht zu erfassen.

Wenn die Anzahl der Cluster während des Clusters automatisch ermittelt wird, müssen Sie sich nicht um die oben genannten Probleme kümmern. Das ** unendlich gemischte Gaußsche Modell ** macht dies möglich.

Das unendlich gemischte Gaußsche Modell basiert auf der Idee, dass es unendlich viele Gaußsche Verteilungen gibt, und die Daten im Datensatz sind einer endlichen Anzahl von ihnen zugeordnet. Die folgende Abbildung zeigt das Ergebnis der Clusterbildung, das später gezeigt wird. (Da es problematisch ist, werde ich die Abbildung erneut verwenden. Bitte haben Sie Verständnis.)

Dieser Datensatz wird durch eine Gaußsche Verteilung von drei dargestellt. Während ein normales gemischtes Gauß-Modell glaubt, dass es genau drei Gauß-Verteilungen auf dieser Welt gibt, glaubt ein unendliches gemischtes Gauß-Modell, dass es unendlich viele potenzielle Gauß-Verteilungen gibt, die nicht an einen Datensatz gebunden sind. Auf diese Weise können Sie neben der Wahrscheinlichkeit, dass Daten mit den drei vorhandenen Clustern verbunden werden, auch die Wahrscheinlichkeit berechnen, mit der die Daten mit einem unbekannten Cluster verbunden werden, bei dem Sie nicht wissen, wo der Durchschnitt liegt.

Mit anderen Worten, beim Clustering mit dem unendlich gemischten Gaußschen Modell werden Daten einem neuen Cluster zugewiesen, während vorhandene Cluster verschwinden, ohne mit Daten verbunden zu sein, die Anzahl der Cluster erhöhen oder verringern und schließlich zu einer geeigneten Anzahl konvergieren. Sie können.

Berechnung der Wahrscheinlichkeitsverteilung

Ich werde die Berechnungsformel von hier aus veröffentlichen. In einer normalen gemischten Gaußschen Verteilung wird die *** kategoriale Verteilung *** auf die Wahrscheinlichkeitsverteilung von $ z_i $ angewendet (wenn $ x_ {1: N} $ [^ 4] nicht beobachtet wird).

python

P(z_i = k | \pi) = \pi_k

Wobei $ \ pi: = (\ pi_1, \ cdots, \ pi_K) $ ein Punkt allein in der Dimension $ (K-1) $ ist. [^ 1] Als Methode zur Bestimmung des Wertes von $ \ pi $ gibt es eine Methode zur Annahme einer Gleichverteilung oder eine Methode zur Bayes'schen Schätzung unter Annahme einer vorherigen Verteilung [^ 2].

Andererseits wendet das unendlich gemischte Gaußsche Modell einen stochastischen Prozess an, der als ** chinesischer Kochprozess ** bezeichnet wird. Dies ist ein Mechanismus zum Berechnen der Wahrscheinlichkeitsverteilung des Clusters $ z_i $ der $ i $ -ten Daten aus dem Cluster $ z_ {1: N} ^ {\ backslash i} $ der Daten außer dem $ i $ th.

python

P(z_i = k | z_{1:N}^{\backslash i}, \alpha) =
\left\{
\begin{array}{ll}
\frac{n_k}{N - 1 + \alpha} & (k = 1, \cdots, K) \\
\frac{\alpha}{N - 1 + \alpha} & (k = K + 1)
\end{array}
\right.

$ k = 1, \ cdots, K $ entspricht dem vorhandenen Cluster und $ k = K + 1 $ entspricht dem neuen Cluster. Hier ist $ n_k $ die Anzahl von $ j $, so dass $ i \ neq j, z_j = k $ und $ \ alpha $ ein Hyperparameter ist, der einen positiven Wert annimmt. Wie Sie der Formel entnehmen können, ist es umso wahrscheinlicher, dass ein neuer Cluster zugewiesen wird, je größer $ \ alpha $ ist. Von nun an werde ich diese Wahrscheinlichkeitsverteilung als $ {\ rm CRP} (k | z_ {1: N} ^ {\ backslash i}, \ alpha) $ schreiben.

Betrachten Sie als nächstes die Wahrscheinlichkeitsverteilung, wenn Sie $ x_ {1: N} $ beobachten. In einem normalen gemischten Gaußschen Modell kann es wie folgt berechnet werden.

python

\begin{align}
P(z_i = k | x_{1:N}, \mu_{1:K}, \Sigma_{1:K}, \pi) 
&= \frac{P(z_i = k | \pi) N(x_i | \mu_k, \Sigma_k)}{\sum_{l=1}^{K}P(z_i = l | \pi) N(x_i | \mu_l, \Sigma_l)} \\
&= \frac{\pi_k N(x_i | \mu_k, \Sigma_k)}{\sum_{l=1}^{K}\pi_l N(x_i | \mu_l, \Sigma_l)}
\end{align}

$ N (x | \ mu, \ Sigma) $ ist die Wahrscheinlichkeitsdichtefunktion der multivariaten Gaußschen Verteilung. In dem unendlich gemischten Gaußschen Modell kann die Wahrscheinlichkeitsverteilung erhalten werden, indem $ P (z_i = k | \ pi) $ in der obigen Gleichung in die Gleichung für den chinesischen Kochprozess geändert wird.

python

\begin{align}
P(z_i = k | x_{1:N}, \mu_{1:K}, \Sigma_{1:K}, z_{1:N}^{\backslash i}, \alpha) 
&= \frac{{\rm CRP}(k | z_{1:N}^{\backslash i}, \alpha) N(x_i | \mu_k, \Sigma_k)}{\sum_{l=1}^{K+1}{\rm CRP}(l | z_{1:N}^{\backslash i}, \alpha) N(x_i | \mu_l, \Sigma_l)}
\end{align}

Hier gibt es eine Einschränkung. Für vorhandene Cluster $ k = 1, \ cdots, K $ ist es möglich, den Durchschnitt $ \ mu_k $ und die verteilte gemeinsam verteilte Matrix $ \ Sigma_k $ aus den zu jedem Cluster gehörenden Daten zu schätzen. Für den neuen Cluster $ k = K + 1 $ ist es jedoch nicht möglich, den Mittelwert und die verteilte, gemeinsam verteilte Matrix direkt zu schätzen, da keine Daten zum Cluster gehören. Daher werden wir dieses Mal die folgende Methode anwenden.

• Nehmen Sie eine vorherige Verteilung für $ \ mu_k $ an, marginalisieren Sie die gleichzeitige Verteilung $ P (x_i, \ mu_k | \ cdots) $ und ermitteln Sie die Wahrscheinlichkeitsverteilung, ohne $ \ mu_k $ explizit zu schätzen
• Lassen Sie $ \ Sigma_k $ die gemeinsame Matrix $ \ sigma ^ 2 I $ für alle Cluster sein. ($ \ Sigma ^ 2 $ ist ein Hyperparameter)

Für den Mittelwert berechnen die Bayes'sche Schätzung und Marginalisierung die Wahrscheinlichkeitsverteilung nur aus der vorherigen Verteilung. In Bezug auf die Verteilungs-Co-Verteilungsmatrix ist es für das Clustering ausreichend, sie wie oben zu fixieren, sodass wir sie diesmal nicht schätzen werden. [^ 5] Natürlich ist es auch möglich, die Varianz abzuschätzen und zu gruppieren.

Da es unsinnig ist, Fälle zu trennen, wird diese Methode auch angewendet, wenn $ k = 1, \ cdots, K + 1 $. Informationen zur Marginalisierung finden Sie in Vorheriger Artikel.

Clustering mit nicht parametrischen Feldern

Um den Cluster $ z_ {1: N} $ von $ N $ Daten zu schätzen, werden wir wahrscheinlich eins nach dem anderen von $ z_1 $ bis $ z_N $ abtasten. Diese Methode wird als ** Gibbs-Abtastung ** bezeichnet, und die Methode, die auf der unendlich dimensionalen Wahrscheinlichkeitsverteilung wie dieser Zeit basiert, wird insbesondere als ** nichtparametrische Buchten ** bezeichnet.

Um diese Methode auszuführen, müssen wir einige Annahmen über das anfängliche Clustering $ z_ {1: N} $ treffen. Dieses Mal wiederholen wir die Stichprobe von $ z_1 $ bis $ z_N $ mehrmals mit demselben Cluster $ k = 1 $ wie das anfängliche Clustering. Anfangs gibt es nur einen Cluster, aber es ist ein Bild, bei dem Daten, die weit vom Durchschnitt entfernt sind, mit einem neuen Cluster kombiniert werden und der Klassenring fortschreitet, während die Anzahl der Cluster erhöht wird.

Implementierung

Implementiertes Clustering mit einem unendlich gemischten Gaußschen Modell. Es basiert auf der Implementierung von Last time, wurde jedoch insgesamt aufgefrischt und der Unterschied kann groß sein.

from collections import Counter
from functools import partial
from typing import Sequence, Tuple
import numpy as np
from scipy.special import logsumexp


def _log_gaussian(x: np.ndarray, mu: np.ndarray, var: float) -> np.ndarray:
    #2 im vorherigen Artikel*np.Ich habe pi weggelassen, aber der Klarheit halber werde ich es diesmal berücksichtigen.
    d = x.shape[0]
    norm_squared = ((x - mu) * (x - mu)).sum(axis=1)
    return - d * np.log(2 * np.pi * var) / 2 - norm_squared / (2 * var)


#Array-Klasse zum Speichern des Clusters
class InfiniteClusterArray(object):
    def __init__(self, array: Sequence[int]):
        self._array = np.array(array, dtype=int)
        self._counter = Counter(array)
        self._n_clusters = max(array) + 1

    def update(self, i: int, k: int):
        assert k <= self._n_clusters

        pre_value = self._array[i]
        if pre_value == k:
            return

        if self._counter[pre_value] > 0:
            self._counter[pre_value] -= 1
        self._array[i] = k
        self._counter[k] += 1

        if k == self._n_clusters:
            self._n_clusters += 1

        if not self._counter[pre_value]:
            self._array -= 1 * (self._array > pre_value)
            self._counter = Counter(self._array)
            self._n_clusters -= 1

    @property
    def array(self) -> np.ndarray:
        return self._array.copy()

    @property
    def n_clusters(self) -> int:
        return self._n_clusters

    def count(self, k: int, excluded=None) -> int:
        assert (excluded is None) or isinstance(excluded, int)
        if excluded is None:
            return self._counter[k]

        counted = self._counter[k] - bool(self._array[excluded] == k)
        return counted

    def __getitem__(self, i: int) -> int:
        return self._array[i]

            

#Array-Klasse zum Speichern des Datasets
class DataArray(object):
    def __init__(self, array: np.ndarray):
        assert array.ndim == 2
        self._array = array

    @property
    def n(self) -> int:
        return self._array.shape[0]

    @property
    def d(self) -> int:
        return self._array.shape[1]

    def mean(self, excluded=None) -> np.ndarray:
        assert (excluded is None) or isinstance(excluded, Sequence)

        if excluded is None:
            return self._array.mean(axis=0)
        else:
            return self._array[excluded].mean(axis=0)

    def __getitem__(self, i):
        return self._array.__getitem__(i)

    def __iter__(self):
        for i in range(self.n):
            yield self._array[i]

    def __str__(self) -> str:
        return f'DataArray({self._array})'


#Unendlich gemischte Gaußsche Gibbs-Stichprobenklasse
class IGMMGibbsSampling(object):
    def __init__(self, var: float, var_pri: float, alpha: float, seed=None):
        #Der vorherige Verteilungsdurchschnitt wird auf den Durchschnitt für den gesamten Datensatz festgelegt.
        assert (var > 0) and (var_pri > 0)
        assert alpha > 0

        self.var = var
        self.var_pri = var_pri
        self.alpha = alpha
        self._random = np.random.RandomState(seed)

    def predict_clusters(self, X: np.ndarray, n_iter: int, init_z=None) -> np.ndarray:
        #Methode zur Schätzung des Clusters
        assert X.ndim == 2
        X = DataArray(X)
        mu_pri = X.mean()

        if init_z is None:
            init_z = np.zeros(X.n, dtype=int)
        z = InfiniteClusterArray(init_z)

        for _ in range(n_iter):
            compute_probs = partial(self._compute_pmf_zi, X, z, mu_pri)
            for i in range(X.n):
                probs = compute_probs(i)
                z_i = self._random.multinomial(1, probs)
                z_i = np.where(z_i)[0][0]  # One-hot to index
                z.update(i, z_i)

        return z.array

    def _compute_pmf_zi(self, X: DataArray, z: InfiniteClusterArray, mu_pri: np.ndarray, i: int) -> np.ndarray:
        mu, var = self._compute_marginal_distribution(X, z, mu_pri, i)
        pi = np.array([
            z.count(k, excluded=i) / (X.n - 1 + self.alpha)
            if k < z.n_clusters
            else self.alpha / (X.n - 1 + self.alpha)
            for k in range(z.n_clusters + 1)
        ])

        log_pdf_xi = _log_gaussian(X[i], mu, var)
        with np.errstate(divide='ignore'):
            pmf_zi = np.exp(np.log(pi) + log_pdf_xi - logsumexp(np.log(pi) + log_pdf_xi)[np.newaxis])

        return pmf_zi

    def _compute_marginal_distribution(self, X: DataArray, z: InfiniteClusterArray,
                                       mu_pri: np.ndarray, i: int) -> Tuple[np.ndarray, float]:

        n_vec = np.array([
            z.count(k, excluded=i)
            for k in range(z.n_clusters + 1)
        ], dtype=np.int)

        tmp_vec = 1 / (n_vec / self.var + 1 / self.var_pri)

        mu = np.tile(mu_pri / self.var_pri, (z.n_clusters + 1, 1))  # Shape (z.n_clusters + 1, X.d)
        for k in range(z.n_clusters):
            excluded_list = [
                j for j in range(X.n)
                if (j != i) and (z[j] == k)
            ]
            if excluded_list:
                mean = X.mean(excluded=excluded_list)
                mu[k] += n_vec[k] * mean / self.var
        mu *= tmp_vec[:, np.newaxis]
        var = tmp_vec + self.var

        return mu, var

Clustering-Praxis

Versuchen Sie als Test das Clustering mit der implementierten selbst erstellten Klasse. Diesmal gilt dies für den folgenden Datensatz.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
%matplotlib inline

np.random.seed(0)
cmap = plt.rcParams['axes.prop_cycle'].by_key()['color']

pi = [0.6, 0.2, 0.2]
mu = np.array([[1, 2], [-2, 0], [3, -2]])
var = [0.8, 0.4, 0.2]
size = 100
ndim = 2

def make_data():
    z = np.random.multinomial(1, pvals=pi)
    z = np.where(z==1)[0][0]
    cov = var[z] * np.eye(2)
    
    return np.random.multivariate_normal(mu[z], cov)

X = np.array([make_data() for _ in range(100)])

fig = plt.figure()
ax = plt.axes()
ax.scatter(X[:, 0], X[:, 1], color='gray')
ax.set_aspect('equal')
plt.show()

Dies ist die nicht gruppierte Version des Datensatzes, den ich Ihnen zuvor gezeigt habe. Die Anzahl der Daten beträgt 100 und wird durch zufällige Auswahl aus 3 Gaußschen Verteilungen erstellt. Daher ist es angebracht, es in drei Cluster zu unterteilen.

Lassen Sie uns jetzt Cluster.

python

model = IGMMGibbsSampling(var=0.5, var_pri=4, alpha=0.1, seed=0)  #Übergeben Sie Hyperparameter an den Konstruktor
clusters = model.predict_clusters(X, n_iter=10)  #Clustering durchführen

#Visualisierung
fig = plt.figure()
ax = plt.axes()
for k in range(max(clusters) + 1):
    data = np.array([
        X[i, :] for i in range(X.shape[0])
        if clusters[i] == k
    ])
    ax.scatter(data[:, 0], data[:, 1], label=f'cluster {k + 1}', color=cmap[k])

ax.set_aspect('equal')
plt.legend()
plt.show()

Die übliche Verteilung var im Cluster ist ein kleiner Datensatz von $ 0,5 $, die Verteilung der mittleren vorherigen Verteilung var_pri ist ein großer $ 4 $, damit sie wie eine gleichmäßige Verteilung aussieht, und der Parameter alpha des chinesischen Kochprozesses beträgt $ 0,1. Es ist auf $ gesetzt. Übrigens wird der Mittelwert der vorherigen Verteilung des Mittelwerts als Mittelwert des gesamten Datensatzes implementiert. Der Clustering-Methode "Predict_Clusters" werden der Datensatz und die Anzahl der Iterationen "n_iter" übergeben. Hier wird die Anzahl der Iterationen auf 10 gesetzt.

Das Ergebnis ist wie folgt.

Obwohl ich die Anzahl der Cluster nicht angegeben und mit demselben Clusterstatus begonnen habe, wurden sie in drei Cluster unterteilt.

Da dies eine große Sache ist, untersuchen wir den Übergang des Clusters für jede Anzahl von Iterationen. Die Clustering-Ergebnisse bei "n_iter = 2, 4, 6, 8, 10" werden angezeigt.

• n_iter = 2

• n_iter = 4

• n_iter = 6

• n_iter = 8

• n_iter = 10

Es scheint, dass eine große Anzahl von Clustern aus dem Anfangszustand erstellt wurde, und dann verringerte sich die Anzahl von Clustern allmählich auf drei. Übrigens, selbst wenn ich "n_iter" von hier aus erhöht habe, hat die Anzahl der Cluster im Grunde 3 Zustände beibehalten.

• Gelegentlich wurde ein Cluster mit nur einem Daten erstellt, und die Anzahl der Cluster wurde zu vier.

Hyperparameteranpassung

Wie oben erwähnt, ist es für neue Cluster umso einfacher, neue Cluster zu bilden, je größer der Hyperparameter $ \ alpha $ im chinesischen Kochprozess ist. Ursprünglich gab es ein Motiv, die Schätzung der Anzahl der Cluster wegzulassen, und wenn Sie Schwierigkeiten haben, dieses $ \ alpha $ anzupassen, wird die Bedeutung der Einführung verringert. Als ich versuchte, $ \ alpha $ zu ändern, hatte ich den Eindruck, dass es im Grunde kein Problem gibt, wenn es kleiner gemacht wird. Die folgende Abbildung zeigt das Ergebnis von $ \ alpha = 0,0001 $ und 10 Iterationen.

Es ist fest in 3 Cluster unterteilt. Selbst wenn $ \ alpha $ ziemlich klein gemacht wird, scheint es, dass es richtig geclustert werden kann, ohne einen großen Einfluss auf die Anzahl der Iterationen zu haben.

Übrigens, als $ \ alpha = 1 $ gesetzt wurde, wurden zum Zeitpunkt von 10 Iterationen 5 Cluster gebildet, und die Anzahl verringerte sich danach nicht auf 3. Es scheint, dass die Anzahl der Cluster nicht gut bestimmt werden kann, wenn $ \ alpha $ zu groß ist. Es hat die gleichen Eigenschaften wie die Schrittweite des beim maschinellen Lernen verwendeten Optimierungsalgorithmus.

abschließend

Bisher habe ich versucht, Clustering mit einem unendlich gemischten Gaußschen Modell durchzuführen. Es war interessant, weil es wie beabsichtigt in die Anzahl der Cluster unterteilt wurde.

[^ 1]: Eine Menge von $ (\ pi_1, \ cdots, \ pi_K) $, die $ \ pi_1, \ cdots, \ pi_K \ ge 0 $ und $ \ sum_ {k = 1} ^ K \ pi_k = 1 $ erfüllt Wird als $ (K-1) $ -Dimensionseinheit bezeichnet. Ursprünglich in der Geometrie verwendet, repräsentiert die Dimensionseinheit $ 0 $ einen Punkt, die Dimensionseinheit $ 1 $ ein Liniensegment und die Dimensionseinheit $ 2 $ ein Dreieck. Wenn Sie das Bayes'sche Schätzbuch lesen, wird dieser Satz im Kontext der Wahrscheinlichkeit verwendet. [^ 2]: Es gibt auch eine Methode, um $ \ pi $ am meisten aus der Beobachtung von $ z_ {1: N} ^ {\ backslash i} $ zu schätzen, aber diese Methode ist immer $ , wenn der Cluster $ k $ verschwindet. Beachten Sie, dass pi_k = 0 $ ist und der Cluster nicht wiederbelebt wird. [^ 4]: $ N $ Variablen $ x_1, \ cdots, x_N $ werden als $ x_ {1: N} $ abgekürzt. Außerdem schreibe ich $ x_ {1: N} ^ {\ backslash i} $ für $ N-1 $ -Variablen mit Ausnahme von $ i $ th. [^ 5]: Auch wenn $ \ Sigma_k $ ein fester Wert ist, der Clustern gemeinsam ist, denke ich, dass Clustering mit der gleichen Genauigkeit wie k-means durchgeführt werden kann. Dies liegt daran, dass die k-means-Methode der wahrscheinlichsten Schätzung von jedem $ \ mu_k $ und $ z_i $ im gemischten Gaußschen Modell mit einer festen Varianz wie oben beschrieben entspricht.