PRML Kapitel 9 Mixed Gaussian Distribution Python-Implementierung

PRML In Kapitel 9 wird der EM-Algorithmus vorgestellt. Der EM-Algorithmus selbst ist eine Technik, die an verschiedenen Orten verwendet werden kann, und ich selbst habe den EM-Algorithmus mit einem Klassifikator kombiniert, um eine rauschresistentere Klassifizierung durchzuführen. Da ich mich jedoch nie mit der gemischten Gaußschen Verteilung gruppiert hatte, die das bekannteste Anwendungsbeispiel für die Anwendung des EM-Algorithmus ist, implementierte ich ** die wahrscheinlichste Schätzung der gemischten Gaußschen Verteilung mit dem EM-Algorithmus **.

Höchstwahrscheinlich Schätzung der gemischten Gaußschen Verteilung

Zum Beispiel eine normale Gaußsche Verteilung von Datenpunkten wie die blauen Punkte in der obigen Abbildung.

p({\bf x}) = \mathcal{N}({\bf x}|{\bf\mu},{\bf\Sigma})

Bei der Anpassung mit, Aus diesem Grund ist die Single-Peak-Gauß-Verteilung in diesem Fall kein gutes Modell. Konzentriert sich auf die Tatsache, dass die Datenpunkte in drei Gruppen unterteilt sind, in diesem Fall eine gemischte Gaußsche Verteilung unter Verwendung von drei Gaußschen Verteilungen.

p({\bf x}) = \sum_{k=1}^3 \pi_k\mathcal{N}({\bf x}|{\bf\mu}_k,{\bf\Sigma}_k)

Wird als geeignetes Modell angesehen. Es sei jedoch $ \ sum_k \ pi_k = 1 $. $ {\ Bf \ mu} \ _1 $ ist oben, $ {\ bf \ mu} \ _2 $ ist unten rechts und $ {\ bf \ mu} \ _3 $ ist der Durchschnitt der Datenpunktblöcke unten links. Ich werde. Daher ist es gut, jeden Block mit einer Gaußschen Verteilung zu versehen, aber für uns Menschen ist es offensichtlich, welcher Datenpunkt zu welchem Block gehört, aber die Maschine weiß das nicht.

Welcher Datenpunkt zu welchem der K Chunks gehört, mit den Koordinaten $ \ {{\ bf x} \ _ n \} \ _ {n = 1} ^ N $ von N Datenpunkten als beobachteten Variablen 1-of-k-Codierung der latenten Variablen $ \ {{\ bf z} \ _n \} \ _ {n = 1} ^ N $ und Parameter $ \ {{\ bf \ pi} \ _k , {\ bf \ mu} \ _ k, {\ bf \ Sigma} \ _ k \} \ _ {k = 1} ^ K $ werden gleichzeitig geschätzt. Wenn es solche latenten Variablen gibt, ist es üblich, den EM-Algorithmus zu verwenden.

EM-Algorithmus

Das Verfahren zur wahrscheinlichsten Schätzung der gemischten Gaußschen Verteilung durch den EM-Algorithmus (PRML-Gleichungen (9.23) bis (9.27)) ist in Abschnitt 9.2.2 von PRML zusammengefasst und wird hier weggelassen.

Code

Bibliothek

Verwenden Sie wie gewohnt Numpy und Matplotlib.

import matplotlib.pyplot as plt
import numpy as np

Gemischte Gaußsche Verteilung

class GaussianMixture(object):

    def __init__(self, n_component):
        #Anzahl der Gaußschen Verteilungen
        self.n_component = n_component

    #Höchstwahrscheinlich Schätzung unter Verwendung des EM-Algorithmus
    def fit(self, X, iter_max=10):
        #Datendimension
        self.ndim = np.size(X, 1)
        #Initialisierung des Mischungskoeffizienten
        self.weights = np.ones(self.n_component) / self.n_component
        #Durchschnittliche Initialisierung
        self.means = np.random.uniform(X.min(), X.max(), (self.ndim, self.n_component))
        #Initialisierung der Kovarianzmatrix
        self.covs = np.repeat(10 * np.eye(self.ndim), self.n_component).reshape(self.ndim, self.ndim, self.n_component)

        #Wiederholen Sie den Schritt E und M.
        for i in xrange(iter_max):
            params = np.hstack((self.weights.ravel(), self.means.ravel(), self.covs.ravel()))
            #E Schritt, berechnen Sie die Belastungsrate
            resps = self.expectation(X)
            #M Schritt, Parameter aktualisieren
            self.maximization(X, resps)
            #Überprüfen Sie, ob die Parameter konvergiert haben
            if np.allclose(params, np.hstack((self.weights.ravel(), self.means.ravel(), self.covs.ravel()))):
                break
        else:
            print("parameters may not have converged")

    #Gaußsche Funktion
    def gauss(self, X):
        precisions = np.linalg.inv(self.covs.T).T
        diffs = X[:, :, None] - self.means
        assert diffs.shape == (len(X), self.ndim, self.n_component)
        exponents = np.sum(np.einsum('nik,ijk->njk', diffs, precisions) * diffs, axis=1)
        assert exponents.shape == (len(X), self.n_component)
        return np.exp(-0.5 * exponents) / np.sqrt(np.linalg.det(self.covs.T).T * (2 * np.pi) ** self.ndim)

    #E Schritt
    def expectation(self, X):
        #PRML-Formel(9.23)
        resps = self.weights * self.gauss(X)
        resps /= resps.sum(axis=-1, keepdims=True)
        return resps

    #M Schritt
    def maximization(self, X, resps):
        #PRML-Formel(9.27)
        Nk = np.sum(resps, axis=0)

        #PRML-Formel(9.26)
        self.weights = Nk / len(X)

        #PRML-Formel(9.24)
        self.means = X.T.dot(resps) / Nk

        diffs = X[:, :, None] - self.means
        #PRML-Formel(9.25)
        self.covs = np.einsum('nik,njk->ijk', diffs, diffs * np.expand_dims(resps, 1)) / Nk

    #Wahrscheinlichkeitsverteilung p(x)Berechnung
    def predict_proba(self, X):
        #PRML-Formel(9.7)
        gauss = self.weights * self.gauss(X)
        return np.sum(gauss, axis=-1)

    #Clustering
    def classify(self, X):
        joint_prob = self.weights * self.gauss(X)
        return np.argmax(joint_prob, axis=1)

Ganzer Code

gaussian_mixture.py

import matplotlib.pyplot as plt
import numpy as np


class GaussianMixture(object):

    def __init__(self, n_component):
        self.n_component = n_component

    def fit(self, X, iter_max=10):
        self.ndim = np.size(X, 1)
        self.weights = np.ones(self.n_component) / self.n_component
        self.means = np.random.uniform(X.min(), X.max(), (self.ndim, self.n_component))
        self.covs = np.repeat(10 * np.eye(self.ndim), self.n_component).reshape(self.ndim, self.ndim, self.n_component)
        for i in xrange(iter_max):
            params = np.hstack((self.weights.ravel(), self.means.ravel(), self.covs.ravel()))
            resps = self.expectation(X)
            self.maximization(X, resps)
            if np.allclose(params, np.hstack((self.weights.ravel(), self.means.ravel(), self.covs.ravel()))):
                break
        else:
            print("parameters may not have converged")

    def gauss(self, X):
        precisions = np.linalg.inv(self.covs.T).T
        diffs = X[:, :, None] - self.means
        assert diffs.shape == (len(X), self.ndim, self.n_component)
        exponents = np.sum(np.einsum('nik,ijk->njk', diffs, precisions) * diffs, axis=1)
        assert exponents.shape == (len(X), self.n_component)
        return np.exp(-0.5 * exponents) / np.sqrt(np.linalg.det(self.covs.T).T * (2 * np.pi) ** self.ndim)

    def expectation(self, X):
        resps = self.weights * self.gauss(X)
        resps /= resps.sum(axis=-1, keepdims=True)
        return resps

    def maximization(self, X, resps):
        Nk = np.sum(resps, axis=0)
        self.weights = Nk / len(X)
        self.means = X.T.dot(resps) / Nk
        diffs = X[:, :, None] - self.means
        self.covs = np.einsum('nik,njk->ijk', diffs, diffs * np.expand_dims(resps, 1)) / Nk

    def predict_proba(self, X):
        gauss = self.weights * self.gauss(X)
        return np.sum(gauss, axis=-1)

    def classify(self, X):
        joint_prob = self.weights * self.gauss(X)
        return np.argmax(joint_prob, axis=1)


def create_toy_data():
    x1 = np.random.normal(size=(100, 2))
    x1 += np.array([-5, -5])
    x2 = np.random.normal(size=(100, 2))
    x2 += np.array([5, -5])
    x3 = np.random.normal(size=(100, 2))
    x3 += np.array([0, 5])
    return np.vstack((x1, x2, x3))


def main():
    X = create_toy_data()

    model = GaussianMixture(3)
    model.fit(X, iter_max=100)
    labels = model.classify(X)

    x_test, y_test = np.meshgrid(np.linspace(-10, 10, 100), np.linspace(-10, 10, 100))
    X_test = np.array([x_test, y_test]).reshape(2, -1).transpose()
    probs = model.predict_proba(X_test)
    Probs = probs.reshape(100, 100)
    colors = ["red", "blue", "green"]
    plt.scatter(X[:, 0], X[:, 1], c=[colors[int(label)] for label in labels])
    plt.contour(x_test, y_test, Probs)
    plt.xlim(-10, 10)
    plt.ylim(-10, 10)
    plt.show()


if __name__ == '__main__':
    main()

Ergebnis

Die Parameter der gemischten Gaußschen Verteilung werden höchstwahrscheinlich unter Verwendung der Punkte als Trainingsdaten geschätzt, und die Wahrscheinlichkeitsverteilung wird durch Konturlinien dargestellt. Die Farbe der Punkte gibt auch an, zu welchem Cluster sie gehören. Dies ist das Ergebnis des Erfolgs, aber ** scheitert manchmal **. Wie in PRML beschrieben, ist die Maximierung der logarithmischen Wahrscheinlichkeitsfunktion dieses Mal ein schlechtes Einstellungsproblem und möglicherweise keine gute Lösung. Es gibt einige Heuristiken, die dies umgehen, aber diesmal schlägt dies fehl, weil wir keine Problemumgehung implementiert haben. Es sollte jedoch nicht sehr scheitern.

Am Ende

Unüberwachtes Clustering wurde durch Anpassen mit einer gemischten Gaußschen Verteilung durchgeführt. Hier wird die Anzahl der zu diesem Zeitpunkt verwendeten Gaußschen Verteilungen angegeben. Im nächsten Kapitel 10 werden wir eine Methode vorstellen, mit der die Anzahl der Elemente einer geeigneten gemischten Gaußschen Verteilung automatisch geschätzt werden kann. Daher werden wir beim nächsten Mal die dort verwendeten Variantenschächte implementieren.