[Python] Ich habe versucht, marginalisiertes Gibbs-Sampling zu implementieren

Dieser Artikel ist eine Fortsetzung meines Artikels [Python] Implementieren von Clustering mit gemischten Gaußschen Modellen. Professionelle Reihe zum maschinellen Lernen "Nichtparametrischer Punktpunktprozess und statistische Mathematik des maschinellen Lernens" wird im Buch gelesen und geschrieben Ich habe das Modell in Python implementiert. Dieses Mal werden wir das zuletzt implementierte Modell verbessern und Clustering mit ** marginalisierter Gibbs-Stichprobe ** durchführen.

In diesem Artikel verwendete Symbole

Hier sind einige der in diesem Artikel verwendeten mathematischen Symbole, die anderswo selten zu finden 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.

Überprüfung der Gibbs-Probenahme

Bevor wir über Marginalisierung sprechen, lassen Sie uns die vorherige Gibbs-Stichprobe etwas genauer betrachten.

Beim letzten Mal haben wir angenommen, dass der Datensatz $ x_ {1: N} $ durch ein gemischtes Gauß-Modell von $ K $ -Clustern dargestellt wird und die durchschnittliche Gauß-Verteilung für jeden Cluster $ \ mu_ {1: K beträgt. } $ Und der Cluster $ z_ {1: N} $ jeder Daten sollte geschätzt werden. Die gleichzeitige Verteilung von $ z_ {1: N} $ und $ \ mu_ {1: K} $ ergibt die folgenden bedingten Wahrscheinlichkeiten, die analytisch schwer zu handhaben sind.

python

P(\mu_{1: K}, z_{1:N} | x_{1: N}, \Sigma_{1: K}, \pi)

In der obigen Gleichung ist $ \ Sigma_ {1: K} $ die Kovarianzmatrix der Gaußschen Verteilung, und $ \ pi: = (\ pi_1, \ cdots, \ pi_K) ^ T $ ist die kategoriale Verteilung, der jedes $ z_i $ folgt. Es ist ein Parameter. [^ 2] Anstatt $ \ mu_ {1: K}, z_ {1: N} $ gleichzeitig zu schätzen, wird daher jede Variable $ \ mu_1, \ cdots, \ mu_K, z_1, \ cdots, z_N $ nacheinander nacheinander geschätzt. Das Ganze wird geschätzt, indem dies abgetastet und wiederholt wird. Diese Technik wird als ** Gibbs-Abtastung ** bezeichnet. Daher wird es in den folgenden zwei Schritten geschätzt.

• Für jedes $ k = 1, \ cdots, K $, $ P (\ mu_k | \ mu_ {1: K} ^ {\ Backslash k}, z_ {1: N}, x_ {1: N}, \ Sigma_ {1: K}, \ pi, \ mu_ {\ rm pri}, \ Sigma_ {\ rm pri}) Beispiel $ \ mu_k $ basierend auf $
• Für jedes $ i = 1, \ cdots, N $, $ P (z_i | z_ {1: N} ^ {\ Backslash i}, \ mu_ {1: K}, x_ {1: N}, \ Sigma_ {1: K}, \ pi) Beispiel $ z_i $ basierend auf $

$ \ Mu_ {\ rm pri}, \ Sigma_ {\ rm pri} $ sind Parameter der vorherigen Verteilung von $ \ mu_k $. [^ 1]

Was ist periphere Gibbs-Abtastung?

In der obigen Gibbs-Stichprobe ist der Teil, der den Datencluster schätzt, der zweite Schritt, und der erste Schritt ist der zusätzliche Prozess, der zur Berechnung des Clusters erforderlich ist. Hier wird unter Verwendung einer Technik namens ** Marginalisierung ** plötzlich der zweite Schritt "$ z_ {1: N} $" ohne den ersten Schritt "Abtastung von $ \ mu_ {1: N} $" ausgeführt. Erwägen Sie "Sampling". Diese Technik wird als ** marginalisierte Gibbs-Abtastung ** bezeichnet. Die periphere Gibbs-Abtastung wiederholt einen Schritt:

• Für jedes $ i = 1, \ cdots, N $, $ P (z_i | z_ {1: N} ^ {\ Backslash i}, x_ {1: N}, \ Sigma_ {1: K}, \ pi , \ mu_ {\ rm pri}, \ Sigma_ {\ rm pri}) Beispiel $ z_i $ basierend auf $

Selbst wenn die Varianz von $ \ mu_ {1: K} $ groß ist [^ 6], wird $ z_ {1 bei der marginalisierten Gibbs-Abtastung nicht durch das instabile Abtastergebnis von $ \ mu_ {1: K} $ beeinflusst. : N} $ kann abgetastet werden. Sie fragen sich vielleicht, wie Sie ohne $ \ mu_ {1: K} $ schätzen können, aber die Informationen für $ \ mu_ {1: K} $ lauten $ z_ {1: N} ^ {\ backslash i}, x_ {1: N}, \ Sigma_ {1: K}, \ pi, \ mu_ {\ rm pri}, \ Sigma_ {\ rm pri} $ haben sie, also benutze sie Bild zum Schätzen von $ z_i $ (ohne Verwendung von $ \ mu_ {1: K} $). Lassen Sie uns nun von der Gibbs-Sampling-Geschichte weggehen und die Marginalisierung diskutieren.

Peripherisierung der gleichzeitigen Verteilung

Angenommen, $ x \ in X, y \ in Y $ ist eine kontinuierliche stochastische Variable und die gleichzeitige Verteilung $ P (x, y) $ ist bekannt. Zu diesem Zeitpunkt wird das Eliminieren der stochastischen Variablen durch Integration wie folgt als ** Marginalisierung ** bezeichnet.

python

P(x) = \int_{Y}P(x, y)dy

Die Transformation der rechten Seite unter Verwendung des Bayes-Theorems und das Schreiben mit bedingter Wahrscheinlichkeit ergibt: Verwenden Sie diese Option, wenn Sie die bedingte Wahrscheinlichkeit und nicht die gleichzeitige Verteilung kennen.

python

P(x) = \int_{Y}P(x|y)P(y)dy

Übrigens, wenn die Wahrscheinlichkeitsvariable $ y $ diskret ist, lautet die Marginalisierungsformulierung wie folgt. Dies ist möglicherweise leichter zu verstehen.

python

P(x) = \sum_{y \in Y}P(x, y) = \sum_{y \in Y}P(x|y)P(y)

Anwendung auf Gibbs-Probenahme

Kehren wir zur Geschichte von Gibbs Sampling zurück. Peripherie, um die Wahrscheinlichkeitsverteilung von $ z_i $ ohne $ \ mu_ {1: K} $ zu finden. $ D $ ist die Dimension von $ x_i $ und $ \ mu_k $. Da es viele bedingte Variablen mit bedingter Wahrscheinlichkeit gibt, werden sie weggelassen.

python

P(z_i=k | \cdots)
= \int_{\mathbb{R}^D} P(z_i=k | \mu_{1:K}, \cdots)P(\mu_{1:K} | \cdots)d\mu_{1:K}

Wie beim letzten Mal wird dies unter der Annahme von $ \ Sigma_k = \ sigma_k ^ 2 I_D, \ Sigma_ {\ rm pri}: = \ sigma_ {\ rm pri} ^ 2I_D $ wie folgt berechnet. [^ 3]

python

\begin{align}
P(z_i=k | \cdots)
&\propto \pi_k N\left(x_i | a_k\left(\frac{n_k^{\backslash i}}{\sigma^2}\overline{x_k}^{\backslash i} + \frac{1}{\sigma_{\rm pri}^2} \mu_{\rm pri}  \right),
 \left( \sigma^2 + a_k \right)I_D\right) \\
a_k &:= \left( \frac{n_k^{\backslash i}}{\sigma^2} + \frac{1}{\sigma_{\rm pri}^2} \right)^{-1}
\end{align}

$ n_k ^ {\ backslash i} $ ist die Anzahl der Daten, die zum Cluster $ k $ im Datensatz gehören, die durch Subtrahieren von $ x_i $ von $ x_ {1: N} $, $ \ override {x_k} ^ {\ erhalten werden Backslash i} $ ist ihr Durchschnitt.

Implementierung

Basierend auf der vorherigen Implementierung werden wir marginalisierte Gibbs-Stichproben implementieren. Es gibt zwei wesentliche Änderungen gegenüber dem letzten Mal.

• Löschen Sie die Klassenvariable mu (unnötig, da sie nicht abgetastet wird)
• Änderung der Methode "log_deformed_gaussian"

Die erste ist die Kernidee der peripheren Gibbs-Abtastung und die wichtigste. Lassen Sie mich den zweiten erklären. Das letzte Mal wurde die Berechnung durchgeführt, um das Verhältnis der Wahrscheinlichkeitsdichtefunktionswerte der Normalverteilung zu erhalten, wobei die Faktoren ignoriert wurden, die nicht vom Cluster "k" abhängen. Ich werde die vorherige Formel erneut veröffentlichen. [^ 4]

python

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

Hier wie die Formel- \frac{1}{2}\log|\Sigma_k|Ich konnte den Faktor von ignorieren\Sigma_kIst ein ClusterkDies liegt daran, dass es sich um einen festen Wert handelte, von dem nichts abhing. Diesmal ändert sich jedoch $ \ Sigma_k = (\ sigma ^ 2 + a_k) I_D $ und $ a_k $ mit $ k $, sodass dieser Faktor nicht ignoriert werden kann. Daher lautet die Formel für diese Zeit wie folgt.

python

\begin{align}
P(z_i = k | x_{1:N}^{\backslash i}, z_{1:N}^{\backslash i},  \mu_{1:K}, \Sigma_{1:K}, \pi) 
&\propto \pi_k \exp\left\{ \frac{D}{2}\log \sigma_k^2 -\frac{1}{2 \sigma^2}\ \| x_i - \mu_k \|_2^2\right\}
\end{align}

Die Implementierung der Methode "log_deformed_gaussian", die den Inhalt dieser "exp" berechnet, ist wie folgt.

python

def log_deformed_gaussian(x, mu, var):
    D = x.shape[0]
    norm_squared = ((x - mu) * (x - mu)).sum(axis=1)
    return -D * np.log(var) / 2 - norm_squared / (2 * var)

Basierend auf dem oben Gesagten habe ich versucht, es zu implementieren. Es ist lang, also ist es gefaltet.

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 -np.log(var) / 2 - 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]


#Klasse für Clustering mit marginalisierter Gibbs-Stichprobe
class GMMClustering(object):
    def __init__(self, X, K, var=1, var_pri=1, seed=None):
        self.X = X.copy()
        self.K = K
        self.N = X.shape[0]
        self.D = X.shape[1]
        self.z = None

        #Parametereinstellung der Wahrscheinlichkeitsverteilung
        # self.mu wird nicht verwendet, ist also nicht definiert
        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, n_iter):
        init_z = self._random.randint(0, self.K, self.N)
        self.z = ClusterArray(init_z)

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

    def _sample_zi(self, i):
        n_vec = np.array([
            self.z.count(k) - bool(k == self.z.array[i])
            for k in range(self.K)
        ])
        x_bar_vec = np.array([self._mean_in_removed(k, i) for k in range(self.K)])

        tmp = 1/(n_vec/self.var + 1/self.var_pri)
        mu = np.tile(tmp, (self.D, 1)).T * (np.tile(n_vec, (self.D, 1)).T * x_bar_vec / self.var + self.mu_pri / self.var_pri)
        var = tmp + self.var

        log_probs_xi = log_deformed_gaussian(self.X[i], mu, 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 _mean_in_removed(self, k, i):
        # self.Berechnen Sie den Durchschnitt der zu Cluster k gehörenden Daten aus dem Datensatz, wobei das i-te von X subtrahiert wird
        mean = np.array([
            x for j, x in enumerate(self.X)
            if (j != i) and (self.z[j] == k)
        ]).mean(axis=0)
        return mean

Ergebnis

Wie beim letzten Mal habe ich versucht, Clustering für Ayame Dataset durchzuführen. Die Hyperparameter $ \ sigma ^ 2 = 0.1, \ sigma_ {\ rm pri} ^ 2 = 1 $ und die Anzahl der Gibbs-Abtastiterationen beträgt $ 10 $. Dies sind die gleichen Einstellungen wie beim letzten Mal. 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(X, K=3, var=0.05, seed=1)
gmc.fit(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()

Als ich das Ergebnis mit "Pair Plot" von "Seaborn" visualisierte, bekam ich Folgendes.

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

Selbst mit marginalisierter Gibbs-Stichprobe, die nicht $ \ mu_ {1: K} $ abtastet, können Sie gut gruppieren.

abschließend

Implementierte periphere Gibbs-Stichproben aus Machine Learning Professional-Reihe "Nichtparametrischer Bayes-Punkt-Prozess und statistische Mathematik des maschinellen Lernens" Es wurde bestätigt, dass dasselbe Clustering wie Letztes Mal durchgeführt werden kann. Das nächste Mal planen wir, den Hauptinhalt dieses Buches, die nicht parametrischen Felder, zu implementieren.

[^ 1]: Da die Wahrscheinlichkeitsverteilung von $ \ mu_k $ Bayesed ist, muss die vorherige Verteilung festgelegt werden. [^ 2]: Dies sind nach wie vor feste Werte, keine stochastischen Variablen. [^ 3]: Das Buch beschreibt diese Ableitungsmethode, aber ich habe die Berechnung selbst überprüft. Es war ziemlich schwierig. [^ 4]: Letztes Mal habe ich die Wahrscheinlichkeitsverteilung selbst beschrieben, aber diesmal habe ich die Proportionalformel beschrieben. Wenn Sie sich daran gewöhnen, wird dieser leichter zu verstehen sein. [^ 6]: Die Verteilung von $ \ mu_k $ ist groß, wenn die Anzahl der zum Cluster $ k $ gehörenden Daten klein ist. In der marginalisierten Gibbs-Stichprobe können Sie als extremes Beispiel die Wahrscheinlichkeit ermitteln, dass $ z_i = k $ auftritt, selbst wenn die zum Cluster $ k $ gehörenden Daten $ 0 $ sind. Diese Eigenschaft ist wichtig für nicht parametrische Felder.