PRML Chapitre 12 Mise en œuvre de l'analyse principale bayésienne Python

Dans cet article, nous allons implémenter l'analyse en composantes principales bayésiennes. Je pense que l'utilisation principale de l'analyse en composantes principales (ACP) est de trouver la projection de l'espace d'observation (haute dimension) où les données cibles existent vers l'espace latent (faible dimension). Si le but est la visualisation, l'espace latent est transformé en 2 (ou 3) dimensions, mais s'il s'agit d'un prétraitement de données, il est rare de savoir combien de dimensions de l'espace latent doivent être définies. Il existe également une méthode de calcul du taux de cotisation, mais en fin de compte, il faut fixer le seuil à ce moment-là. Par conséquent, dans l'analyse bayésienne en composantes principales, la dimension de l'espace latent est automatiquement déterminée par la détermination automatique du degré de pertinence.

Analyse probabiliste en composantes principales

En interprétant l'analyse en composantes principales de manière probabiliste, il sera possible de la gérer comme Bayes plus tard. Dans l'analyse probabiliste en composantes principales, les données $ x $ (dimension D) que nous avons observées projettent $ z $ (dimension M) échantillonnées à partir de l'espace latent avec la matrice $ W $ ((D, M) matrice). Ensuite, il est interprété comme se déplaçant en parallèle ($ + \ mu $ (dimension D)) et ajoutant du bruit ($ + \ epsilon $ (dimension D)).

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

Supposons maintenant que $ z $ et $ \ epsilon $ suivent une distribution gaussienne. Il y a trois paramètres à estimer: $ W, \ mu, \ sigma ^ 2 $ (la variance de la distribution gaussienne suivie de $ \ epsilon $), et sa fonction de vraisemblance est la suivante.

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

Deux méthodes pour maximiser cette fonction de vraisemblance sont introduites dans PRML. La première est une méthode qui utilise simplement la décomposition de singularité, et la seconde est lorsque la mise à jour de la distribution postérieure (étape E) pour $ z $ et les données complètes $ x, z $ sont données en utilisant l'algorithme EM. Il répète la maximisation (étape M) de la fonction de vraisemblance de.

Analyse en composantes principales bayésiennes

Dans l'exemple ci-dessus, la dimension de l'espace variable latent a été fixée à M. Dans l'analyse bayésienne en composantes principales, les dimensions supplémentaires sont élaguées à l'aide de la détermination automatique de la pertinence. (Cependant, la valeur de M ne diminue pas réellement.) Par conséquent, la distribution précédente suivante est fournie pour le paramètre $ W $.

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\}

Où $ {\ bf w} \ _i $ est le vecteur de colonne $ i $ ème de $ {\ bf W} $. Et $ \ alpha \ _i $ agit comme un paramètre de précision pour les distributions gaussiennes individuelles. En estimant $ \ alpha $, certains composants ont des valeurs très élevées. Dans ce cas, la précision est très élevée, donc la composante du vecteur colonne correspondant de $ {\ bf} W $ n'est que de 0, c'est-à-dire que la dimension est élaguée.

code

Bibliothèque

N'utilisez que matplotlib et numpy comme d'habitude.

import matplotlib.pyplot as plt
import numpy as np

Analyse des composants principaux par la méthode la plus probable

Méthode utilisant la décomposition des valeurs propres comme d'habitude

#Classe d'analyse en composantes principales
class PCA(object):

    def __init__(self, n_component):
        #Spécifiez la dimension de l'espace latent
        self.n_component = n_component

    #Effectuer l'analyse des composants principaux par la méthode la plus probable
    def fit(self, X):
        #Formule PRML(12.1)Calculer l'estimation la plus probable de mu
        self.mean = np.mean(X, axis=0)

        #Formule PRML(12.2)Matrice de covariance des données
        cov = np.cov(X, rowvar=False)

        #Décomposition de valeur unique
        values, vectors = np.linalg.eigh(cov)
        index = np.size(X, 1) - self.n_component

        #Formule PRML(12.46) sigma^Estimation la plus probable de 2
        if index == 0:
            self.var = 0
        else:
            self.var = np.mean(values[:index])

        #Formule PRML(12.45)Estimation la plus probable de la matrice de projection W
        self.W = vectors[:, index:].dot(np.sqrt(np.diag(values[index:]) - self.var * np.eye(self.n_component)))

Analyse en composantes principales bayésiennes

Les méthodes ci-dessous sont toutes les méthodes de la classe PCA. Cette fois, à des fins de comparaison, les paramètres sont d'abord estimés les plus probables, et ils sont utilisés comme valeurs initiales pour l'élagage par détermination automatique de la pertinence.

    #Effectuer une analyse en composantes principales bayésiennes
    def fit_bayesian(self, X, iter_max=100):
        #Dimensions de l'espace de données
        self.ndim = np.size(X, 1)

        #Estimer le paramètre par l'estimation la plus probable et l'utiliser comme valeur initiale
        self.fit(X)

        #Initialisation des paramètres de précision(Première estimation)
        self.alpha = self.ndim / np.sum(self.W ** 2, axis=0)

        #Moyenne zéro des données
        D = X - self.mean

        #Répétez l'algorithme EM un nombre spécifié de fois
        for i in xrange(iter_max):
            #Statistiques suffisantes pour l'étape z
            Ez, Ezz = self.expectation(D)

            #M étape W,sigma^2 estimations
            self.maximize(D, Ez, Ezz)

            #Formule PRML(12.62)Mise à jour des super paramètres
            self.alpha = self.ndim / np.sum(self.W ** 2, axis=0).clip(min=1e-10)

    #Statistiques suffisantes pour l'étape z E[z]、E[zz^T]Calculs de
    def expectation(self, D):
        #Formule PRML(12.41)
        M = self.W.T.dot(self.W) + self.var * np.eye(self.n_component)
        Minv = np.linalg.inv(M)

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

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

    #M étape W,sigma^2 estimations
    def maximize(self, D, Ez, Ezz):
        #Formule PRML(12.63)Estimation de W
        self.W = D.T.dot(Ez).dot(np.linalg.inv(np.sum(Ezz, axis=0) + self.var * np.diag(self.alpha)))

        #Formule PRML(12.57) sigma^2 estimations
        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)

Diagramme de Hinton

Cette fois, nous reproduirons la figure 12.14 du PRML. Pour ce faire, nous avons besoin d'une fonction pour dessiner un diagramme de Hinton qui représente chaque élément de la matrice sous forme de carré. La page publiée par google avec matplotlib hinton est légèrement modifiée.

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()

Fonction principale

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))

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


def main():
    #Créez 100 points de données créés en projetant d'un espace latent en 3 dimensions vers un espace en 10 dimensions
    X = create_toy_data(sample_size=100, ndim_hidden=3, ndim_observe=10, std=1.)

    #Effectuer l'ACP par la méthode la plus probable avec l'espace latent en 9 dimensions
    pca = PCA(9)
    pca.fit(X)
    hinton(pca.W)

    #Taille jusqu'à 9 dimensions pour réaliser l'ACP bayésienne
    pca.fit_bayesian(X)
    hinton(pca.W)

Code entier

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()

résultat

Le résultat de l'estimation de la matrice de projection W par PCA par la méthode la plus probable est le suivant. Ensuite, lorsque les dimensions de l'espace latent sont élaguées à l'aide de l'analyse bayésienne en composantes principales, la matrice de projection W devient comme le montre la figure ci-dessous. Les six colonnes de gauche ont disparu et les dimensions correspondantes ont été élaguées. Il est possible de saisir le fait qu'il a été projeté à partir de trois dimensions. PRML Les résultats présentés sur la figure 12.14 ont été obtenus.

À la fin

Lorsque PRML est arrivé à ce point, il est devenu plus courant de combiner ce que nous avions appris jusqu'à présent et de l'appliquer à de nouveaux modèles. Cette fois, en combinant la détermination automatique du degré de pertinence au chapitre 6 et l'algorithme EM au chapitre 9, les données d'observation sont effectivement appliquées au modèle qui est projeté à partir d'un espace dimensionnel inférieur. Il semble que cela puisse être appliqué aux données avec des valeurs manquantes, je voudrais donc essayer cela également.