PRML Kapitel 12 Bayesianische Hauptanalyse Python-Implementierung

In diesem Artikel werden wir die Bayes'sche Hauptkomponentenanalyse implementieren. Ich denke, die Hauptanwendung der Hauptkomponentenanalyse (PCA) besteht darin, die Projektion vom Beobachtungsraum (hohe Dimension), in dem die Zieldaten existieren, zum latenten Raum (niedrige Dimension) zu finden. Wenn der Zweck die Visualisierung ist, wird der latente Raum in zwei (oder drei) Dimensionen unterteilt. Wenn es sich jedoch um eine Datenvorverarbeitung handelt, ist es selten bekannt, wie viele Dimensionen des latenten Raums festgelegt werden sollen. Es gibt auch eine Methode zur Berechnung des Beitragssatzes, aber am Ende müssen wir den Schwellenwert zu diesem Zeitpunkt festlegen. Daher wird bei der Bayes'schen Hauptkomponentenanalyse die Dimension des latenten Raums automatisch durch die automatische Bestimmung des Relevanzgrades bestimmt.

Probabilistische Hauptkomponentenanalyse

Durch probabilistische Interpretation der Hauptkomponentenanalyse wird es möglich sein, sie später wie Bayes zu behandeln. In der probabilistischen Hauptkomponentenanalyse beobachteten die von uns beobachteten Daten $ x $ (D-Dimension) Projekte $ z $ (M-Dimension), die aus dem latenten Raum mit der Matrix $ W $ ((D, M) -Matrix) abgetastet wurden. Dann wird es so interpretiert, dass es sich parallel bewegt ($ + \ mu $ (D-Dimension)) und Rauschen hinzufügt ($ + \ epsilon $ (D-Dimension)).

{\bf x = Wz + \mu + \epsilon}

Nehmen wir an, dass $ z $ und $ \ epsilon $ einer Gaußschen Verteilung folgen. Es sind drei Parameter zu schätzen: $ W, \ mu, \ sigma ^ 2 $ (die Varianz der Gaußschen Verteilung, gefolgt von $ \ epsilon $), und ihre Wahrscheinlichkeitsfunktion ist wie folgt.

p({\bf x|W,\mu},\sigma^2) = \int p({\bf x|Wz+\mu},\sigma^2)p({\bf z}){\rm d}{\bf z}

In PRML werden zwei Methoden zur Maximierung dieser Wahrscheinlichkeitsfunktion eingeführt. Das erste ist eine Methode, die einfach die Singularitätszerlegung verwendet, und das zweite ist, wenn die Aktualisierung der posterioren Verteilung (E-Schritt) für $ z $ und die vollständigen Daten $ x, z $ unter Verwendung des EM-Algorithmus angegeben werden. Es wiederholt die Maximierung (M-Schritt) der Wahrscheinlichkeitsfunktion von.

Bayesianische Hauptkomponentenanalyse

Im obigen Beispiel wurde die Dimension des latenten variablen Raums auf M festgelegt. Bei der Bayes'schen Hauptkomponentenanalyse werden zusätzliche Dimensionen mithilfe der automatischen Relevanzbestimmung beschnitten. (Der Wert von M nimmt jedoch nicht tatsächlich ab.) Daher wird für den Parameter $ W $ die folgende vorherige Verteilung bereitgestellt.

p({\bf W}|{\bf \alpha}) = \prod_{i=1}^M\left({\alpha_i\over2\pi}\right)^{D/2}\exp\left\{-{1\over2}\alpha_i{\bf w}_i^\top {\bf w}_i\right\}

Wobei $ {\ bf w} \ _i $ der $ i $ -te Spaltenvektor von $ {\ bf W} $ ist. Und $ \ alpha \ _i $ fungiert als Genauigkeitsparameter für die einzelnen Gaußschen Verteilungen. Bei der Schätzung von $ \ alpha $ haben einige Komponenten sehr große Werte. In diesem Fall ist die Genauigkeit sehr hoch, sodass die Komponente des entsprechenden Spaltenvektors von $ {\ bf} W $ nur 0 ist, dh die Dimension wird beschnitten.

Code

Bibliothek

Verwenden Sie wie gewohnt nur Matplotlib und Numpy.

import matplotlib.pyplot as plt
import numpy as np

Hauptkomponentenanalyse nach der wahrscheinlichsten Methode

Methode unter Verwendung der Eigenwertzerlegung wie gewohnt

#Klasse für die Hauptkomponentenanalyse
class PCA(object):

    def __init__(self, n_component):
        #Geben Sie die Dimension des latenten Raums an
        self.n_component = n_component

    #Führen Sie die Hauptkomponentenanalyse nach der wahrscheinlichsten Methode durch
    def fit(self, X):
        #PRML-Formel(12.1)Berechnen Sie die wahrscheinlichste Schätzung von mu
        self.mean = np.mean(X, axis=0)

        #PRML-Formel(12.2)Datenkovarianzmatrix
        cov = np.cov(X, rowvar=False)

        #Einzigartige Wertzerlegung
        values, vectors = np.linalg.eigh(cov)
        index = np.size(X, 1) - self.n_component

        #PRML-Formel(12.46) sigma^Höchstwahrscheinlich Schätzung von 2
        if index == 0:
            self.var = 0
        else:
            self.var = np.mean(values[:index])

        #PRML-Formel(12.45)Höchstwahrscheinlich Schätzung der Projektionsmatrix W.
        self.W = vectors[:, index:].dot(np.sqrt(np.diag(values[index:]) - self.var * np.eye(self.n_component)))

Bayesianische Hauptkomponentenanalyse

Die folgenden Methoden sind alle Methoden in der PCA-Klasse. Zum Vergleich werden diesmal zum ersten Mal die Parameter höchstwahrscheinlich geschätzt, und diese werden als Anfangswerte für das Beschneiden durch automatische Bestimmung der Relevanz verwendet.

    #Führen Sie eine Bayes'sche Hauptkomponentenanalyse durch
    def fit_bayesian(self, X, iter_max=100):
        #Datenraumdimensionen
        self.ndim = np.size(X, 1)

        #Schätzen Sie den Parameter anhand der wahrscheinlichsten Schätzung und verwenden Sie ihn als Anfangswert
        self.fit(X)

        #Initialisierung von Präzisionsparametern(Erste Schätzung)
        self.alpha = self.ndim / np.sum(self.W ** 2, axis=0)

        #0 Mittelung von Daten
        D = X - self.mean

        #Wiederholen Sie den EM-Algorithmus eine bestimmte Anzahl von Malen
        for i in xrange(iter_max):
            #Ausreichende Statistik für Schritt z
            Ez, Ezz = self.expectation(D)

            #M Schritt W.,sigma^2 Schätzungen
            self.maximize(D, Ez, Ezz)

            #PRML-Formel(12.62)Super Parameter Update
            self.alpha = self.ndim / np.sum(self.W ** 2, axis=0).clip(min=1e-10)

    #Ausreichende Statistik für Schritt z E.[z]、E[zz^T]Berechnung von
    def expectation(self, D):
        #PRML-Formel(12.41)
        M = self.W.T.dot(self.W) + self.var * np.eye(self.n_component)
        Minv = np.linalg.inv(M)

        #PRML-Formel(12.54) E[z]
        Ez = D.dot(self.W).dot(Minv)

        #PRML-Formel(12.55) E[zz^T]
        Ezz = self.var * Minv + np.einsum('ni,nj->nij', Ez, Ez)
        return Ez, Ezz

    #M Schritt W.,sigma^2 Schätzungen
    def maximize(self, D, Ez, Ezz):
        #PRML-Formel(12.63)Schätzung von W.
        self.W = D.T.dot(Ez).dot(np.linalg.inv(np.sum(Ezz, axis=0) + self.var * np.diag(self.alpha)))

        #PRML-Formel(12.57) sigma^2 Schätzungen
        self.var = np.mean(
            np.mean(D ** 2, axis=-1)
            - 2 * np.mean(Ez.dot(self.W.T) * D, axis=-1)
            + np.trace(Ezz.dot(self.W.T).dot(self.W).T) / self.ndim)

Hinton-Diagramm

Dieses Mal werden wir PRML Abbildung 12.14 reproduzieren. Dazu benötigen wir eine Funktion zum Zeichnen eines Hinton-Diagramms, das jedes Element der Matrix als Quadrat darstellt. Seite, die von Google mit matplotlib hinton veröffentlicht wurde wurde geringfügig geändert.

def hinton(matrix, max_weight=None, ax=None):
    """Draw Hinton diagram for visualizing a weight matrix."""
    ax = ax if ax is not None else plt.gca()

    if not max_weight:
        max_weight = 2 ** np.ceil(np.log(np.abs(matrix).max()) / np.log(2))

    ax.patch.set_facecolor('gray')
    ax.set_aspect('equal', 'box')
    ax.xaxis.set_major_locator(plt.NullLocator())
    ax.yaxis.set_major_locator(plt.NullLocator())

    for (x, y), w in np.ndenumerate(matrix):
        color = 'white' if w > 0 else 'black'
        size = np.sqrt(np.abs(w) / max_weight)
        rect = plt.Rectangle([y - size / 2, x - size / 2], size, size,
                             facecolor=color, edgecolor=color)
        ax.add_patch(rect)

    ax.autoscale_view()
    ax.invert_yaxis()
    plt.xlim(-0.5, np.size(matrix, 1) - 0.5)
    plt.ylim(-0.5, len(matrix) - 0.5)
    plt.show()

Hauptfunktion

def create_toy_data(sample_size=100, ndim_hidden=1, ndim_observe=2, std=1.):
    Z = np.random.normal(size=(sample_size, ndim_hidden))
    mu = np.random.uniform(-5, 5, size=(ndim_observe))
    W = np.random.uniform(-5, 5, (ndim_hidden, ndim_observe))

    #PRML-Formel(12.33)
    X = Z.dot(W) + mu + np.random.normal(scale=std, size=(sample_size, ndim_observe))
    return X


def main():
    #Erstellen Sie 100 Datenpunkte, die durch Projektion von einem dreidimensionalen latenten Raum in einen 10-dimensionalen Raum erstellt wurden
    X = create_toy_data(sample_size=100, ndim_hidden=3, ndim_observe=10, std=1.)

    #Führen Sie die PCA nach der wahrscheinlichsten Methode mit dem latenten Raum als 9 Dimensionen durch
    pca = PCA(9)
    pca.fit(X)
    hinton(pca.W)

    #Beschneiden aus bis zu 9 Dimensionen zur Durchführung der Bayes'schen PCA
    pca.fit_bayesian(X)
    hinton(pca.W)

Ganzer Code

pca.py

import matplotlib.pyplot as plt
import numpy as np


class PCA(object):

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

    def fit(self, X):
        self.mean = np.mean(X, axis=0)
        cov = np.cov(X, rowvar=False)
        values, vectors = np.linalg.eigh(cov)
        index = np.size(X, 1) - self.n_component
        if index == 0:
            self.var = 0
        else:
            self.var = np.mean(values[:index])
        self.W = vectors[:, index:].dot(np.sqrt(np.diag(values[index:]) - self.var * np.eye(self.n_component)))

    def fit_bayesian(self, X, iter_max=100):
        self.ndim = np.size(X, 1)
        self.fit(X)
        self.alpha = self.ndim / np.sum(self.W ** 2, axis=0)
        D = X - self.mean
        for i in xrange(iter_max):
            Ez, Ezz = self.expectation(D)
            self.maximize(D, Ez, Ezz)
            self.alpha = self.ndim / np.sum(self.W ** 2, axis=0).clip(min=1e-10)

    def expectation(self, D):
        M = self.W.T.dot(self.W) + self.var * np.eye(self.n_component)
        Minv = np.linalg.inv(M)
        Ez = D.dot(self.W).dot(Minv)
        Ezz = self.var * Minv + np.einsum('ni,nj->nij', Ez, Ez)
        return Ez, Ezz

    def maximize(self, D, Ez, Ezz):
        self.W = D.T.dot(Ez).dot(np.linalg.inv(np.sum(Ezz, axis=0) + self.var * np.diag(self.alpha)))
        self.var = np.mean(
            np.mean(D ** 2, axis=-1)
            - 2 * np.mean(Ez.dot(self.W.T) * D, axis=-1)
            + np.trace(Ezz.dot(self.W.T).dot(self.W).T) / self.ndim)


def hinton(matrix, max_weight=None, ax=None):
    """Draw Hinton diagram for visualizing a weight matrix."""
    ax = ax if ax is not None else plt.gca()

    if not max_weight:
        max_weight = 2 ** np.ceil(np.log(np.abs(matrix).max()) / np.log(2))

    ax.patch.set_facecolor('gray')
    ax.set_aspect('equal', 'box')
    ax.xaxis.set_major_locator(plt.NullLocator())
    ax.yaxis.set_major_locator(plt.NullLocator())

    for (x, y), w in np.ndenumerate(matrix):
        color = 'white' if w > 0 else 'black'
        size = np.sqrt(np.abs(w) / max_weight)
        rect = plt.Rectangle([y - size / 2, x - size / 2], size, size,
                             facecolor=color, edgecolor=color)
        ax.add_patch(rect)

    ax.autoscale_view()
    ax.invert_yaxis()
    plt.xlim(-0.5, np.size(matrix, 1) - 0.5)
    plt.ylim(-0.5, len(matrix) - 0.5)
    plt.show()


def create_toy_data(sample_size=100, ndim_hidden=1, ndim_observe=2, std=1.):
    Z = np.random.normal(size=(sample_size, ndim_hidden))
    mu = np.random.uniform(-5, 5, size=(ndim_observe))
    W = np.random.uniform(-5, 5, (ndim_hidden, ndim_observe))
    X = Z.dot(W) + mu + np.random.normal(scale=std, size=(sample_size, ndim_observe))
    return X


def main():
    X = create_toy_data(sample_size=100, ndim_hidden=3, ndim_observe=10, std=1.)

    pca = PCA(9)
    pca.fit(X)
    hinton(pca.W)

    pca.fit_bayesian(X)
    hinton(pca.W)


if __name__ == '__main__':
    main()

Ergebnis

Das Ergebnis der Schätzung der Projektionsmatrix W durch PCA nach dem wahrscheinlichsten Verfahren ist wie folgt. Wenn dann die Dimension des latenten Raums unter Verwendung der Bayes'schen Hauptkomponentenanalyse beschnitten wird, wird die Projektionsmatrix W wie in der folgenden Abbildung gezeigt. Die 6 Spalten links sind verschwunden und die entsprechenden Abmessungen wurden beschnitten. Es ist möglich, die Tatsache zu erfassen, dass es aus drei Dimensionen projiziert wurde. PRML Die in Abbildung 12.14 gezeigten Ergebnisse wurden erhalten.

Am Ende

Als PRML an diesen Punkt kam, wurde es üblicher, das bisher Gelernte zu kombinieren und auf neue Modelle anzuwenden. Dieses Mal werden durch Beobachtung der automatischen Bestimmung des Relevanzgrads in Kapitel 6 und des EM-Algorithmus in Kapitel 9 die Beobachtungsdaten tatsächlich auf das Modell angewendet, das aus einem Raum mit niedrigeren Dimensionen projiziert wird. Es scheint, dass es auf Daten mit fehlenden Werten angewendet werden kann, also würde ich das auch gerne versuchen.