[Python] Implementierung von Clustering mit einem gemischten Gaußschen Modell

Professionelle Serie zum maschinellen Lernen "Nichtparametrischer Punktpunktprozess und statistische Mathematik des maschinellen Lernens" wird im Buch gelesen und geschrieben Ich habe das Modell in Python implementiert. Nichtparametrische Felder sind eine Erweiterung des gemischten Gaußschen Modells. Ich möchte wirklich nichtparametrische Felder implementieren, aber in Vorbereitung darauf werde ich ein gemischtes Gauß-Modell implementieren.

• Es gab einen Artikel über Qiita, in dem das gleiche Modell und die gleiche Methode sowie die Implementierung erläutert wurden. Daher möchte ich ihn vorstellen. Beide waren sehr hilfreich.

• "Einführung in das maschinelle Lernen durch Bayes'sche Inferenz" Ungefähre Inferenz eines Poisson-Mischmodells, das nur mit Python numpy implementiert wurde
• Gemischte Gaußsche Verteilung und logsumexp

• Die in diesem Artikel beschriebenen Formeln enthalten diejenigen, die Sie anhand der Bücher berechnet haben.

In diesem Artikel verwendete Symbole

Bevor wir zur Hauptgeschichte kommen, werfen wir einen Blick auf einige der in diesem Artikel verwendeten mathematischen Symbole, die anderswo selten zu sehen sind.

• $ x_ {1: N} $: $ x_1, \ cdots, x_N $
• $ x_ {1: N} ^ {\ Backslash i} $: $ x_ {1: N} $ minus $ x_i $

Darüber hinaus ist der Vektor in diesem Artikel ein Spaltenvektor und wird nicht fett dargestellt.

Was ist ein gemischtes Gaußsches Modell?

Ein Gaußsches Mischungsmodell (GMM) ist ein Modell eines Datensatzes, in dem Daten, die aus mehreren Gaußschen Verteilungen generiert wurden, gemischt werden. Nehmen wir als Beispiel Ayame Dataset. Dies ist ein Datensatz der Eigenschaften und Arten von Iris, und drei Arten von Iris "setosa", "versicolor" und "virginica" werden gemischt. Viele von Ihnen sind damit vertraut, da es häufig für Übungen zur Klassifizierung des maschinellen Lernens verwendet wird. Es gibt vier erklärende Variablen in diesem Datensatz, aber der Einfachheit halber haben wir zwei herausgenommen, "petal_length" und "petal_width", und ein Streudiagramm der "farbcodierten Version" und "farbcodierten Version für jeden Iris-Typ" erstellt. Ich werde.

import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
df = sns.load_dataset('iris')
ax1.scatter(df['petal_length'], df['petal_width'], color='gray')
ax1.set_title('without label')
ax1.set_xlabel('petal_length')
ax1.set_ylabel('petal_width')
for sp in df['species'].unique():
    x = df[df['species'] == sp]['petal_length']
    y = df[df['species'] == sp]['petal_width']
    ax2.scatter(x, y, label=sp)
ax2.legend()
ax2.set_title('with label')
ax2.set_xlabel('petal_length')
ax2.set_ylabel('petal_width')
plt.show()

In der Abbildung rechts kann interpretiert werden, dass alle Daten gemäß der multivariaten Gaußschen Verteilung [^ 1] generiert werden, die für jeden Iris-Typ einen anderen Mittelwert hat. Dies ist das gemischte Gaußsche Modell. Ein grober Durchschnitt ist in der folgenden Tabelle dargestellt.

Unter Verwendung des gemischten Gaußschen Modells kann der Durchschnitt mehrerer Gaußscher Verteilungen aus dem links gezeigten unbeschrifteten Datensatz geschätzt werden, und die Daten und jede Gaußsche Verteilung können zum Clustering verknüpft werden.

Mathematische Erklärung

Das gemischte Gaußsche Modell wird als probabilistisches Modell beschrieben. Der Datencluster $ x_i $ sei $ z_i $ und die multivariate Gaußsche Verteilung, die jedem Cluster $ 1, \ cdots, k $ entspricht, sei $ N (\ mu_k, \ Sigma_k) $. Mit anderen Worten, wenn $ x_i $ zum Clasper $ k $ gehört, kann seine Generierungswahrscheinlichkeit wie folgt geschrieben werden: $ D $ ist die Dimension von $ x_i $ und $ N $ ist die Anzahl der Daten. ..

python

\begin{align}
P(x_i | z_i=k, x_{1:N}^{\backslash i}, z_{1:N}^{\backslash i},  \mu_{1:K}, \Sigma_{1:K}) 
&= N(x | \mu_k, \Sigma_k)
\end{align}

Andererseits wird $ z_i $ auch probabilistisch ausgedrückt. Angenommen, $ z_i $ wird aus einer kategorialen Verteilung mit $ \ pi: = (\ pi_1, \ cdots, \ pi_K) ^ T $ als Parametern generiert. $ K $ ist die Anzahl der Cluster.

python

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

Aus dem Obigen kann die Wahrscheinlichkeit, dass der Cluster $ z_i $ der Daten $ x_i $ $ k $ ist, wie folgt geschrieben werden. [^ 2]

python

P(z_i = k | x_{1:N}^{\backslash i}, z_{1:N}^{\backslash i},  \mu_{1:K}, \Sigma_{1:K}, \pi)
= \frac{\pi_k N(x_i| \mu_k, \Sigma_k)}{\sum_{l=1}^K \pi_l N(x_i| \mu_l, \Sigma_l)} \qquad\cdots (1)

Das Clustering erfolgt durch Schätzen jedes $ z_i $ aus dem Datensatz unter Verwendung dieser Formel.

So schätzen Sie den Cluster

Um den Cluster $ z_ {1: N} $ aller Daten aus der Formel $ (1) $ zu schätzen, entsprechen die Parameter der multivariaten Gaußschen Verteilung jedem Cluster $ \ mu_ {1: K}, \ Sigma_ { Sie müssen 1: K} $ und den kategorialen Verteilungsparameter $ \ pi $ kennen. Diese werden nicht explizit angegeben und sollten aus dem Datensatz geschätzt oder als geeignete Werte angenommen werden. Dieses Mal schätzen wir den Mittelwert $ \ mu_ {1: K} $ der multivariaten Gaußschen Verteilung aus dem Datensatz und nehmen die verbleibenden Parameter auf die folgenden Werte an: $ I_D $ ist die Einheitsmatrix von $ D \ times D $ und $ \ sigma ^ 2 $ ist der Hyperparameter. [^ 6]

python

\Sigma_1 = \cdots = \Sigma_K = \sigma^2 I_D \\
\pi = (1 / K, \cdots, 1 / K)^T

Um dann $ \ mu_ {1: K} $ und $ z_ {1: N} $ zu schätzen, werden wir die Methode anwenden, diese nacheinander probabilistisch abzutasten. Diese Methode wird als ** Gibbs-Abtastung ** bezeichnet. [^ 9] [^ 10]

Da die Abtastung von $ z_i $ gemäß der Formel $ (1) $ erfolgen kann, berücksichtigen Sie die Wahrscheinlichkeitsverteilung von $ \ mu_k $ unten. Ich möchte die Wahrscheinlichkeitsverteilung schätzen, daher scheint es, dass die Bayes'sche Schätzung verwendet werden kann. Der Durchschnitt $ \ mu_ {\ rm pri} $ der konjugierten vorherigen Verteilung von $ \ mu_k $ und die Kovarianzmatrix $ \ Sigma_ {\ rm pri} $ sind wie folgt definiert. [^ 7] Diese sind allen $ k = 1, \ cdots, K $ gemeinsam. $ \ Sigma_ {\ rm pri} ^ 2 $ ist ein Hyperparameter.

python

\begin{align}
\mu_{\rm pri} &= (0, \cdots, 0)^T \\
\Sigma_{\rm pri} &= \sigma_{\rm pri}^2I_D
\end{align}

Zu diesem Zeitpunkt ist die hintere Verteilung von $ \ mu_k $ wie folgt. $ n_k $ ist die Anzahl der Daten, die zum Cluster $ k $ in $ x_ {1: N} $ gehören, und $ \ bar {x} _k $ ist der Durchschnitt von ihnen.

python

\begin{align}
\mu_{k, {\rm pos}} &= \Sigma_{k, {\rm pos}}(n_k \Sigma_k^{-1} \overline{x}_k + \Sigma_{\rm pri}^{-1}\mu_{\rm pri})  \\
\Sigma_{k, {\rm pos}} &= (n_k \Sigma_{k}^{-1} + \Sigma_{\rm pri}^{-1})^{-1} \\
\mu_k &= N(\mu | \mu_{k, {\rm pos}}, \Sigma_{k, {\rm pos}}) \qquad\cdots (2)
\end{align}

Diesmal sind $ \ Sigma_k $ und $ \ Sigma_ {\ rm pri} ^ {-1} $ konstante Vielfache von $ I_D $, also $ \ mu_ {k, {\ rm pos}} $ und $ \ Sigma_ { k, {\ rm pos}} $ kann wie folgt transformiert werden.

python

\begin{align}
\mu_{k, {\rm pos}} &= \Sigma_{k, {\rm pos}}\left(\frac{n_k}{\sigma^2} \overline{x}_k + \frac{1}{\sigma_{\rm pri}^2}\mu_{\rm pri}\right)  \\
\Sigma_{k, {\rm pos}} &= \left(\frac{n_k}{\sigma^2} + \frac{1}{\sigma_{\rm pri}^{2}}\right)^{-1}
\end{align}

Aus dem Obigen ist ersichtlich, dass die Abtastung von $ \ mu_k $ gemäß Gleichung (2) durchgeführt werden sollte.

Implementierung

Hier sind einige wichtige Punkte bei der Implementierung.

Implementierung eines Arrays, das die Anzahl der Daten im Cluster zählen kann

Die Anzahl der Daten, über die jeder Cluster verfügt, muss aus der Formel $ (2) $ berechnet werden. Implementieren Sie daher die Array-Klasse "ClusterArray" mit einer "count" -Methode, die die Nummer zurückgibt. [^ 8] Zusätzlich zur Methode "count" werden die Methoden, die wahrscheinlich verwendet werden, von "numpy.ndarray" delegiert.

python

import numpy as np
from collections import Counter

class ClusterArray(object):
    def __init__(self, array):
        #Array ist eine eindimensionale Liste, Array
        self._array = np.array(array, dtype=np.int)
        self._counter = Counter(array)

    @property
    def array(self):
        return self._array.copy()

    def count(self, k):
        return self._counter[k]

    def __setitem__(self, i, k):
        #Selbst wenn ausgeführt._Zähler wird ebenfalls aktualisiert
        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

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

Entfernen Sie zusätzliche Berechnungen

Sie können die Wahrscheinlichkeitsdichtefunktion direkt in der Formel $ (1) $ berechnen, aber Sie können den Rechenaufwand reduzieren, indem Sie Faktoren löschen, die nicht mit $ k $ zusammenhängen.

python

\begin{align}
P(z_i = k | x_{1:N}^{\backslash i}, z_{1:N}^{\backslash i},  \mu_{1:K}, \Sigma_{1:K}, \pi)
&= \frac{\pi_k \exp\{- \frac{1}{2}\log|\Sigma_k| - \frac{1}{2} (x_i - \mu_k)^T\Sigma_k(x_i - \mu_k)\}}{\sum_{l=1}^K \pi_l \exp\{- \frac{1}{2}\log|\Sigma_l|- \frac{1}{2} (x_i - \mu_l)^T\Sigma_l(x_i - \mu_l)\}}  \\
&= \frac{\pi_k \exp\{- \frac{1}{2 \sigma^2}\ \| x_i - \mu_k \|_2^2\}}{\sum_{l=1}^{K}\pi_l \exp\{- \frac{1}{2 \sigma^2}\ \| x_i - \mu_l \|_2^2\}}
\end{align}

Implementieren Sie die Methode log_deformed_gaussian, um den Inhalt von $ \ exp $ in dieser Formel zu berechnen und für die Berechnung zu verwenden.

python

def log_deformed_gaussian(x, mu, var):
    norm_squared = ((x - mu) * (x - mu)).sum(axis=1)
    return -norm_squared / (2 * var)

Über logsumexp

Für Berechnungen wie $ \ log (\ sum_ {i = 1} \ exp f_i (x)) $ ist logsumexp eine Technik, um Über- und Unterlauf zu verhindern. Dies ist jedoch nicht sehr effizient, da $ \ log $ und $ \ exp $ häufig verwendet werden. Daher werden wir dieses Mal kein logsumexp verwenden. [^ 5]

Für logsumexp ist der folgende Artikel sehr hilfreich. Gemischte Gaußsche Verteilung und logsumexp

Gesamtumsetzung

Ich habe versucht, es zu implementieren, wobei ich die oben genannten Punkte berücksichtigt habe. Unter Berücksichtigung von "Scikit-Learn" kann das Clustering mit der "Fit" -Methode ausgeführt werden. Der Quellcode ist lang, also falte ich ihn.

import numpy as np
from collections import Counter


def log_deformed_gaussian(x, mu, var):
    norm_squared = ((x - mu) * (x - mu)).sum(axis=1)
    return -norm_squared / (2 * var)


class ClusterArray(object):
    def __init__(self, array):
        #Array ist eine eindimensionale Liste, Array
        self._array = np.array(array, dtype=np.int)
        self._counter = Counter(array)

    @property
    def array(self):
        return self._array.copy()

    def count(self, k):
        return self._counter[k]

    def __setitem__(self, i, k):
        #Selbst wenn ausgeführt._Zähler wird ebenfalls aktualisiert
        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

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


class GaussianMixtureClustering(object):
    def __init__(self, K, D, var=1, var_pri=1, seed=None):
        self.K = K  #Anzahl der Cluster
        self.D = D  #Erklärende variable Abmessungen(Zum Zeitpunkt der Erstellung des Konstruktors festgelegt, um die Implementierung zu vereinfachen)
        self.z = None

        #Parametereinstellung der Wahrscheinlichkeitsverteilung
        self.mu = np.zeros((self.K, self.D))
        self.var = var  #Fest, allen Clustern gemeinsam
        self.pi = np.full(self.K, 1 / self.K)  #Fest, allen Clustern gemeinsam

        #Festlegen der vorherigen Verteilung
        self.mu_pri = np.zeros(self.D)
        self.var_pri = var_pri

        self._random = np.random.RandomState(seed)

    def fit(self, X, n_iter):
        init_z = self._random.randint(0, self.K, X.shape[0])
        self.z = ClusterArray(init_z)

        for _ in range(n_iter):
            for k in range(self.K):
                self.mu[k] = self._sample_mu_k(X, k)
            for i, x_i in enumerate(X):
                self.z[i] = self._sample_zi(x_i)

    def _sample_zi(self, x_i):
        log_probs_xi = log_deformed_gaussian(x_i, self.mu, self.var)

        probs_zi = np.exp(log_probs_xi) * self.pi
        probs_zi = probs_zi / probs_zi.sum()

        z_i = self._random.multinomial(1, probs_zi)
        z_i = np.where(z_i)[0][0]
        return z_i

    def _sample_mu_k(self, X, k):
        xk_bar = np.array([x for i, x in enumerate(X) if self.z[i] == k]).mean(axis=0)
        var_pos = 1 / (self.z.count(k) / self.var + 1 / self.var_pri)
        mu_pos = var_pos * (xk_bar * self.z.count(k) / self.var + self.mu_pri / self.var_pri)

        mu_k = self._random.multivariate_normal(mu_pos, var_pos * np.eye(self.D))
        return mu_k

Versuchen Sie es mit Clustering

Lassen Sie uns den öffnenden Ayame-Datensatz mithilfe des implementierten gemischten Gaußschen Modells gruppieren. Die Hyperparameter sind $ \ sigma ^ 2 = 0,1, \ sigma_ {\ rm pri} ^ 2 = 1 $, und die Anzahl der Gibbs-Abtastiterationen beträgt 10. Wenn wir die Beschriftungen des tatsächlichen Datensatzes mit den Clustering-Ergebnissen vergleichen, erhalten wir:

import seaborn as sns
import matplotlib.pyplot as plt
%matplotlib inline


#Datensatz laden
df = sns.load_dataset('iris')
X = df[['sepal_length', 'sepal_width', 'petal_length', 'petal_width']].values

#Clustering mit gemischtem Gaußschen Modell
gmc = GMMClustering(K=3, D=4, var=0.1, seed=1)
gmc.fit(X, n_iter=10)
df['GMM_cluster'] = gmc.z.array

#Visualisierung der Ergebnisse
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
for sp in df['species'].unique():
    x = df[df['species'] == sp]['petal_length']
    y = df[df['species'] == sp]['petal_width']
    ax1.scatter(x, y, label=sp)
ax1.legend()
ax1.set_title('species')
ax1.set_xlabel('petal_length')
ax1.set_ylabel('petal_width')
for k in range(gmc.K):
    x = df[df['GMM_cluster'] == k]['petal_length']
    y = df[df['GMM_cluster'] == k]['petal_width']
    ax2.scatter(x, y, label=k)
ax2.legend()
ax2.set_title('GMM cluster')
ax2.set_xlabel('petal_length')
ax2.set_ylabel('petal_width')
plt.show()

Die Grenze zwischen "versicolor" und "virginica" ist verdächtig, aber die Clusterbildung gemäß dem Etikett ist fast vollständig. Ich habe auch versucht, das Clustering-Ergebnis mit "Pair Plot" von "Seaborn" zu visualisieren.

sns.pairplot(
    df.drop(columns=['species']),
    vars=['sepal_length', 'sepal_width', 'petal_length', 'petal_width'],
    hue='GMM_cluster'
)

Es ist schön gruppiert. Wenn Sie sich die diagonalen Diagramme von "Pairplot" ansehen, können Sie sehen, dass der Datensatz durch ein gemischtes Gaußsches Modell dargestellt wird.

abschließend

Professionelle Serie zum maschinellen Lernen "Nichtparametrischer Punktpunktprozess und statistische Mathematik des maschinellen Lernens" implementiert ein gemischtes Gaußsches Modell und ist einfach Ich habe versucht, mit verschiedenen Datensätzen zu gruppieren. Dieses Mal haben wir uns mit gemischten Gaußschen Modellen befasst, aber Sie können sich auch Modelle wie "Mixed Bernoulli model" und "Mixed Poisson model" vorstellen, die für das Clustering verwendet werden können. Das nächste Mal werde ich über marginalisierte Gibbs-Sampling schreiben. Dies ist eine notwendige Technik zum Durchführen nicht parametrischer Felder.

(Hinzugefügt am 21. Januar 2020) Ich konnte den nächsten Artikel schreiben. [Python] Ich habe marginalisiertes Gibbs-Sampling implementiert

[^ 1]: Die Verteilung ist unterschiedlich, der Einfachheit halber wird jedoch nur der Durchschnitt berücksichtigt. [^ 2]: Kann mit dem Bayes-Theorem abgeleitet werden. [^ 5]: Ich habe entschieden, dass es in Ordnung ist, logsumexp nicht zu verwenden, da ich Über- und Unterlauf vermeiden kann, indem ich $ \ sigma ^ 2 $ setze, das weder zu groß noch zu klein für den Datensatz ist. [^ 6]: $ \ sigma ^ 2 $ wird verteilt. Alle Kovarianzen sind $ 0 $. Es mag wie eine grobe Annahme erscheinen, aber wenn Sie $ \ sigma ^ 2 $ richtig entscheiden, können Sie immer noch genug Cluster bilden, und wenn Sie es vereinfachen, können Sie den Rechenaufwand reduzieren. Das Schätzen dieser Daten aus dem Datensatz sollte natürlich eine genauere Clusterbildung ermöglichen. [^ 7]: Die konjugierte vorherige Verteilung von $ \ mu_k $ ist eine multivariate Gaußsche Verteilung. [^ 8]: Ich denke, Sie können die Counter-Instanz draußen lassen, ohne eine solche Klasse zu erstellen. Ich mag es, diese Dinge in Klassen zusammenzufassen. [^ 9]: Gibbs-Sampling ist eine Technik zur Approximation simultaner Verteilungen, die analytisch schwer zu berechnen sind. Dieses Mal approximieren wir $ P (\ mu_ {1: K}, z_ {1: N} | \ cdots) $. [^ 10]: Andere Methoden können verwendet werden, um die gemischte Gaußsche Verteilung abzuschätzen, aber wir erkennen, dass Gibbs-Abtastung erforderlich ist, um sich auf nicht parametrische Buchten auszudehnen. Um ehrlich zu sein, verstehe ich es nicht richtig. Bitte lassen Sie mich wissen, ob Sie es verstehen können.