PRML Kapitel 10 Variante Mixed Gaussian Distribution Python-Implementierung

Dieses Mal habe ich die variable gemischte Gaußsche Verteilung implementiert, die am Ende des vorherigen Artikels erwähnt wurde. Letztes Mal haben wir den EM-Algorithmus verwendet, um die Parameter der gemischten Gaußschen Verteilung am meisten zu schätzen, aber dieses Mal werden wir die variante Bayes-Methode verwenden, um die Parameter der gemischten Gaußschen Verteilung zu schätzen. In der vorherigen Implementierung wurde die Anzahl der gemischten Elemente in der gemischten Gaußschen Verteilung festgelegt, aber unter Verwendung der varianten Bayes-Methode wird die Anzahl der gemischten Elemente automatisch bestimmt. Es gibt einige andere Vorteile, wenn man es wie einen Bayesianer behandelt.

Variantengemischte Gaußsche Verteilung

Bei der wahrscheinlichsten Schätzung einer normalen gemischten Gaußschen Verteilung wurden der Mischungskoeffizient $ \ pi_k $, der Mittelwert $ \ mu_k $ und die Kovarianz $ \ Sigma_k $ (oder die Genauigkeit $ \ Lambda_k $) höchstwahrscheinlich geschätzt, diesmal jedoch alle Schätzen Sie die Wahrscheinlichkeitsverteilung als stochastische Variable.

Die folgende Abbildung (aus PRML Abbildung 10.5 extrahiert) ist ein grafisches Modell des diesmal verwendeten Modells. graphical_model.png Die gleichzeitige Verteilung all dieser stochastischen Variablen

p({\bf X}, {\bf Z}, {\bf \pi}, {\bf \mu}, {\bf \Lambda}) = p({\bf X}|{\bf Z},{\bf\mu},{\bf\Lambda})p({\bf Z}|{\bf\pi})p({\bf\pi})p({\bf\mu}|{\bf\Lambda})p({\bf\Lambda})

Es wird sein.

Der Punkt dieser Bayes-Variante ist die posteriore Verteilung $ p ({\ bf Z}, {\ bf \ pi}, {\ bf \ mu}, {nach Beobachtung von ** $ {\ bf X} $. \ bf \ Lambda} | {\ bf X}) $ ist eine differenzielle Annäherung an die Wahrscheinlichkeitsverteilung, die in die latente Variable $ {\ bf Z} $ und die folgenden Parameter zerlegt ist **.

\begin{align}
p({\bf Z}, {\bf \pi}, {\bf \mu}, {\bf \Lambda}|{\bf X}) &\approx q({\bf Z}, {\bf \pi}, {\bf \mu}, {\bf \Lambda})\\
&= q({\bf Z})q({\bf\pi},{\bf\mu},{\bf\Lambda})
\end{align}

Da es schwierig ist, die ursprüngliche posteriore Verteilung genau zu berechnen, verwenden wir eine variante Näherung.

Der allgemeine Ablauf dieser Variante der Bayes-Methode ist

  1. Initialisieren Sie die Wahrscheinlichkeitsverteilungen $ q ({\ bf Z}) $ und $ q ({\ bf \ pi}, {\ bf \ mu}, {\ bf \ Lambda}) $
  2. Aktualisieren Sie $ q ({\ bf Z}) $ (entspricht dem E-Schritt des EM-Algorithmus)
  3. Aktualisieren Sie $ q ({\ bf \ pi}, {\ bf \ mu}, {\ bf \ Lambda}) $ (entspricht dem M-Schritt des EM-Algorithmus)
  4. Wiederholen Sie die Schritte 2 und 3 bis zur Konvergenz Es sieht aus wie das. Wenn Sie die Formeln in den Schritten 2 und 3 aufschreiben, wird es eine beträchtliche Menge sein, also werde ich es weglassen.

Code

Bibliothek

Zusätzlich zu den üblichen Matplotlib und Numpy verwende ich auch Scipy. Da wir eine Digammafunktion und eine Gammafunktion benötigen, werden wir auch für den Algorithmus eine andere Bibliothek als Numpy verwenden.

import matplotlib.cm as cm
import matplotlib.pyplot as plt
import numpy as np
from scipy.special import digamma, gamma

Variante gemischte Gauß

class VariationalGaussianMixture(object):

    def __init__(self, n_component=10, alpha0=1.):

        #Anzahl gemischter Elemente in einer gemischten Gaußschen Verteilung
        self.n_component = n_component

        #Parameter der vorherigen Verteilung des Mischungskoeffizienten pi
        self.alpha0 = alpha0

    def init_params(self, X):
        self.sample_size, self.ndim = X.shape

        #Stellen Sie vorherige Verteilungsparameter ein
        self.alpha0 = np.ones(self.n_component) * self.alpha0
        self.m0 = np.zeros(self.ndim)
        self.W0 = np.eye(self.ndim)
        self.nu0 = self.ndim
        self.beta0 = 1.

        self.component_size = self.sample_size / self.n_component + np.zeros(self.n_component)

        #Initialisieren Sie die Parameter der Wahrscheinlichkeitsverteilung
        self.alpha = self.alpha0 + self.component_size
        self.beta = self.beta0 + self.component_size
        indices = np.random.choice(self.sample_size, self.n_component, replace=False)
        self.m = X[indices].T
        self.W = np.tile(self.W0, (self.n_component, 1, 1)).T
        self.nu = self.nu0 + self.component_size

    #Gibt die Parameter der Wahrscheinlichkeitsverteilung zurück
    def get_params(self):
        return self.alpha, self.beta, self.m, self.W, self.nu

    #Variante Bayes-Methode
    def fit(self, X, iter_max=100):
        self.init_params(X)

        #Aktualisieren, bis die Wahrscheinlichkeitsverteilung konvergiert
        for i in xrange(iter_max):
            params = np.hstack([array.flatten() for array in self.get_params()])

            #Wahrscheinlichkeitsverteilung q(z)Aktualisieren
            r = self.e_like_step(X)

            #Wahrscheinlichkeitsverteilung q(pi, mu, Lambda)Aktualisieren
            self.m_like_step(X, r)

            #Wenn es konvergiert, endet es
            if np.allclose(params, np.hstack([array.ravel() for array in self.get_params()])):
                break
        else:
            print("parameters may not have converged")

    #Wahrscheinlichkeitsverteilung q(z)Aktualisieren
    def e_like_step(self, X):
        d = X[:, :, None] - self.m
        gauss = np.exp(
            -0.5 * self.ndim / self.beta
            - 0.5 * self.nu * np.sum(
                np.einsum('ijk,njk->nik', self.W, d) * d,
                axis=1)
        )
        pi = np.exp(digamma(self.alpha) - digamma(self.alpha.sum()))
        Lambda = np.exp(digamma(self.nu - np.arange(self.ndim)[:, None]).sum(axis=0) + self.ndim * np.log(2) + np.linalg.slogdet(self.W.T)[1])

        #PRML-Formel(10.67)Berechnung der Belastungsrate
        r = pi * np.sqrt(Lambda) * gauss

        #PRML-Formel(10.49)Normalisieren Sie die Belastungsrate
        r /= np.sum(r, axis=-1, keepdims=True)

        #Falls 0% auftritt
        r[np.isnan(r)] = 1. / self.n_component
        return r

    #Wahrscheinlichkeitsverteilung q(pi, mu, Lambda)Aktualisieren
    def m_like_step(self, X, r):

        #PRML-Formel(10.51)
        self.component_size = r.sum(axis=0)

        #PRML-Formel(10.52)
        Xm = X.T.dot(r) / self.component_size

        d = X[:, :, None] - Xm
        #PRML-Formel(10.53)
        S = np.einsum('nik,njk->ijk', d, r[:, None, :] * d) / self.component_size

        #PRML-Formel(10.58) q(pi)Aktualisieren Sie die Parameter für
        self.alpha = self.alpha0 + self.component_size

        #PRML-Formel(10.60)
        self.beta = self.beta0 + self.component_size

        #PRML-Formel(10.61)
        self.m = (self.beta0 * self.m0[:, None] + self.component_size * Xm) / self.beta

        d = Xm - self.m0[:, None]
        #PRML-Formel(10.62)
        self.W = np.linalg.inv(
            np.linalg.inv(self.W0)
            + (self.component_size * S).T
            + (self.beta0 * self.component_size * np.einsum('ik,jk->ijk', d, d) / (self.beta0 + self.component_size)).T).T

        #PRML-Formel(10.63)
        self.nu = self.nu0 + self.component_size

    # p(Testdaten|Trainingsdaten)Pi,mu,Ungefähr mit Lambdas Schätzung
    def predict_proba(self, X):

        #PRML-Formel(B.80)Erwarteter Wert der Präzisionsmatrix Lambda
        covs = self.nu * self.W
        precisions = np.linalg.inv(covs.T).T

        d = X[:, :, None] - self.m
        exponents = np.sum(np.einsum('nik,ijk->njk', d, precisions) * d, axis=1)
        #PRML-Formel(10.38)
        gausses = np.exp(-0.5 * exponents) / np.sqrt(np.linalg.det(covs.T).T * (2 * np.pi) ** self.ndim)

        #PRML-Formel(10.69)Erwarteter Wert des Mischungskoeffizienten pi
        gausses *= (self.alpha0 + self.component_size) / (self.n_component * self.alpha0 + self.sample_size)

        #Peripherisierung der latenten Variablen z
        return np.sum(gausses, axis=-1)

    #Clustering
    def classify(self, X):
        #Wahrscheinlichkeitsverteilung q(Z)Argmax
        return np.argmax(self.e_like_step(X), 1)

    #Berechnen Sie die Verteilung der Schüler
    def student_t(self, X):
        nu = self.nu + 1 - self.ndim

        #PRML-Formel(10.82)
        L = nu * self.beta * self.W / (1 + self.beta)

        d = X[:, :, None] - self.m
        #PRML-Formel(B.72)
        maha_sq = np.sum(np.einsum('nik,ijk->njk', d, L) * d, axis=1)

        #PRML-Formel(B.68)
        return (
            gamma(0.5 * (nu + self.ndim))
            * np.sqrt(np.linalg.det(L.T))
            * (1 + maha_sq / nu) ** (-0.5 * (nu + self.ndim))
            / (gamma(0.5 * nu) * (nu * np.pi) ** (0.5 * self.ndim)))

    #Voraussichtliche Verteilung p(Testdaten|Trainingsdaten)Berechnung
    def predict_dist(self, X):
        #PRML-Formel(10.81)
        return (self.alpha * self.student_t(X)).sum(axis=-1) / self.alpha.sum()

Ganzer Code

import matplotlib.cm as cm
import matplotlib.pyplot as plt
import numpy as np
from scipy.special import digamma, gamma


class VariationalGaussianMixture(object):

    def __init__(self, n_component=10, alpha0=1.):
        self.n_component = n_component
        self.alpha0 = alpha0

    def init_params(self, X):
        self.sample_size, self.ndim = X.shape
        self.alpha0 = np.ones(self.n_component) * self.alpha0
        self.m0 = np.zeros(self.ndim)
        self.W0 = np.eye(self.ndim)
        self.nu0 = self.ndim
        self.beta0 = 1.

        self.component_size = self.sample_size / self.n_component + np.zeros(self.n_component)
        self.alpha = self.alpha0 + self.component_size
        self.beta = self.beta0 + self.component_size
        indices = np.random.choice(self.sample_size, self.n_component, replace=False)
        self.m = X[indices].T
        self.W = np.tile(self.W0, (self.n_component, 1, 1)).T
        self.nu = self.nu0 + self.component_size

    def get_params(self):
        return self.alpha, self.beta, self.m, self.W, self.nu

    def fit(self, X, iter_max=100):
        self.init_params(X)
        for i in xrange(iter_max):
            params = np.hstack([array.flatten() for array in self.get_params()])
            r = self.e_like_step(X)
            self.m_like_step(X, r)
            if np.allclose(params, np.hstack([array.ravel() for array in self.get_params()])):
                break
        else:
            print("parameters may not have converged")

    def e_like_step(self, X):
        d = X[:, :, None] - self.m
        gauss = np.exp(
            -0.5 * self.ndim / self.beta
            - 0.5 * self.nu * np.sum(
                np.einsum('ijk,njk->nik', self.W, d) * d,
                axis=1)
        )
        pi = np.exp(digamma(self.alpha) - digamma(self.alpha.sum()))
        Lambda = np.exp(digamma(self.nu - np.arange(self.ndim)[:, None]).sum(axis=0) + self.ndim * np.log(2) + np.linalg.slogdet(self.W.T)[1])
        r = pi * np.sqrt(Lambda) * gauss
        r /= np.sum(r, axis=-1, keepdims=True)
        r[np.isnan(r)] = 1. / self.n_component
        return r

    def m_like_step(self, X, r):
        self.component_size = r.sum(axis=0)
        Xm = X.T.dot(r) / self.component_size
        d = X[:, :, None] - Xm
        S = np.einsum('nik,njk->ijk', d, r[:, None, :] * d) / self.component_size
        self.alpha = self.alpha0 + self.component_size
        self.beta = self.beta0 + self.component_size
        self.m = (self.beta0 * self.m0[:, None] + self.component_size * Xm) / self.beta
        d = Xm - self.m0[:, None]
        self.W = np.linalg.inv(
            np.linalg.inv(self.W0)
            + (self.component_size * S).T
            + (self.beta0 * self.component_size * np.einsum('ik,jk->ijk', d, d) / (self.beta0 + self.component_size)).T).T
        self.nu = self.nu0 + self.component_size

    def predict_proba(self, X):
        covs = self.nu * self.W
        precisions = np.linalg.inv(covs.T).T
        d = X[:, :, None] - self.m
        exponents = np.sum(np.einsum('nik,ijk->njk', d, precisions) * d, axis=1)
        gausses = np.exp(-0.5 * exponents) / np.sqrt(np.linalg.det(covs.T).T * (2 * np.pi) ** self.ndim)
        gausses *= (self.alpha0 + self.component_size) / (self.n_component * self.alpha0 + self.sample_size)
        return np.sum(gausses, axis=-1)

    def classify(self, X):
        return np.argmax(self.e_like_step(X), 1)

    def student_t(self, X):
        nu = self.nu + 1 - self.ndim
        L = nu * self.beta * self.W / (1 + self.beta)
        d = X[:, :, None] - self.m
        maha_sq = np.sum(np.einsum('nik,ijk->njk', d, L) * d, axis=1)
        return (
            gamma(0.5 * (nu + self.ndim))
            * np.sqrt(np.linalg.det(L.T))
            * (1 + maha_sq / nu) ** (-0.5 * (nu + self.ndim))
            / (gamma(0.5 * nu) * (nu * np.pi) ** (0.5 * self.ndim)))

    def predict_dist(self, X):
        return (self.alpha * self.student_t(X)).sum(axis=-1) / self.alpha.sum()


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 = VariationalGaussianMixture(n_component=10, alpha0=0.01)
    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_dist(X_test)
    plt.scatter(X[:, 0], X[:, 1], c=labels, cmap=cm.get_cmap())
    plt.contour(x_test, y_test, probs.reshape(100, 100))
    plt.xlim(-10, 10)
    plt.ylim(-10, 10)
    plt.show()


if __name__ == '__main__':
    main()

Ergebnis

Das folgende Video zeigt das Clustering mit der Bayes-Variante. Die Farbe der Punkte gibt an, zu welchem Cluster sie gehören, und mit fortschreitendem Lernen nimmt die Anzahl der gültigen gemischten Elemente ab, und schließlich werden es drei. (Der obige Code zeigt keinen solchen Fortschritt und nur das Endergebnis.) animation.gif

Am Ende

In Kapitel 10 von PRML werden Beispiele für die Verwendung der varianten Bayes-Methode für lineare Regression und logistische Regression sowie für gemischte Gauß vorgestellt. Wir werden sie daher implementieren, wenn wir die Möglichkeit dazu haben.

Recommended Posts

PRML Kapitel 10 Variante Mixed Gaussian Distribution Python-Implementierung
PRML Kapitel 9 Mixed Gaussian Distribution Python-Implementierung
PRML Kapitel 6 Gaussian Return Python-Implementierung
PRML Kapitel 2 Python-Implementierung von Student t-Distribution
PRML Kapitel 5 Python-Implementierung für neuronale Netze
PRML Kapitel 3 Evidence Ungefähre Python-Implementierung
PRML Kapitel 8 Summe der Produkte Algorithmus Python-Implementierung
PRML Kapitel 4 Bayesianische logistische Regression Python-Implementierung
PRML Kapitel 5 Python-Implementierung eines Netzwerks mit gemischter Dichte
PRML Kapitel 14 Bedingte gemischte Modell-Python-Implementierung
PRML Kapitel 1 Bayesian Curve Fitting Python-Implementierung
PRML Kapitel 11 Implementierung der Markov-Kette Monte Carlo Python
PRML Kapitel 12 Bayesianische Hauptanalyse Python-Implementierung
Python-Implementierung gemischte Bernoulli-Verteilung
PRML Kapitel 7 Verwandte Vector Machine Python-Implementierung für Regressionsprobleme
Erläuterung und Implementierung von PRML Kapitel 4
Implementierung einer gemischten Normalverteilung in Python
PRML Kapitel 2 Wahrscheinlichkeitsverteilung Nichtparametrische Methode
PRML Kapitel 13 Wahrscheinlichste Schätzung Python-Implementierung des Hidden-Markov-Modells
Schätzung der gemischten Gaußschen Verteilung nach der varianten Bayes'schen Methode
PRML-Implementierung Kapitel 3 Lineares Basisfunktionsmodell
Implementiert in Python PRML Kapitel 7 Nichtlineare SVM
Implementiert in Python PRML Kapitel 5 Neuronales Netzwerk
Implementiert in Python PRML Kapitel 1 Bayesianische Schätzung
Implementiert in Python PRML Kapitel 3 Bayesianische lineare Regression
Implementiert in Python PRML Kapitel 1 Polygonkurvenanpassung
[Python] Implementierung von Clustering mit einem gemischten Gaußschen Modell
Implementiert in Python PRML Kapitel 4 Klassifizierung nach Perceptron-Algorithmus
Implementierung der Bayes'schen Varianzschätzung des Themenmodells in Python
RNN-Implementierung in Python
ValueObject-Implementierung in Python
[Python] Kapitel 01-01 Über Python (Erster Python)
SVM-Implementierung in Python