PRML Chapter 5 Implémentation Python de réseau à densité mixte

La dernière fois, j'ai implémenté un réseau neuronal ordinaire sans utiliser de bibliothèques telles que tensorflow et chainer. Cependant, il existe déjà de nombreuses implémentations de ce type sur le net, et ce n'est pas très intéressant, j'ai donc implémenté un réseau à densité mixte ** cette fois en guise de préparation pour la dernière fois. Puisque le code du réseau neuronal implémenté la dernière fois est utilisé, veuillez vous référer à Article précédent selon le cas.

Réseau à densité mixte

Problèmes avec les modèles conventionnels

Lors de la résolution du problème de régression, il peut ne pas être possible de le traiter correctement en utilisant la fonction d'erreur de la somme des carrés. La figure ci-dessous (reproduction de PRML Figure 5.19) montre le résultat lorsque le réseau de neurones a été formé en utilisant la fonction d'erreur de somme des carrés en utilisant les points bleus comme données d'apprentissage. inverse_problem.png Lorsque $ x = 0,5 $, la valeur de $ y $ semble bonne à 0,2, 0,5 ou 0,8. Cependant, comme je dois en choisir un, il devient 0,5 et je ne rentre pas dans les deux autres. ** J'ai formé le réseau de neurones mais je ne m'intégrais pas très bien car la fonction de coût est modélisée sur une distribution gaussienne monomodale **. La relation entre l'entrée $ x $ et la cible $ t $ est exprimée dans une distribution gaussienne comme suit.

p(t|x) = \mathcal{N}(t|\mu(x),\sigma^2),

Sur cette base, l'ajustement et l'estimation de la fonction $ \ mu (x) $ ont été effectués. Et l'estimation $ \ mu (x) $ est rouge dans la figure ci-dessus, ce qui est un mauvais résultat.

Distribution gaussienne mixte

Le réseau à densité mixte résout le problème ci-dessus en utilisant une fonction de coût modélisée avec une distribution gaussienne mixte multimodale. Préparez un coefficient de mélange $ \ pi \ _k $ pour chaque K distribution gaussienne, et trouvez la relation entre l'entrée $ x $ et la cible $ t $.

p(t|x) = \sum_{k=1}^K\pi_k(x)\mathcal{N}(t|\mu_k(x),\sigma^2_k(x))

Modèle avec. Cependant, la somme des coefficients de mélange est 1 ($ \ sum_k \ pi_k = 1 $). Dans l'exemple de données d'entraînement ci-dessus, à $ x = 0,5 $, il y a trois ensembles de points de données autour de $ y = 0,2,0,5,0,8 $. Par conséquent, il semble bon de définir $ K = 3 $. En ajustant les trois distributions gaussiennes à différents ensembles de points de données, même les données multimodales peuvent être traitées. Par conséquent, étant donné l'ensemble de données d'entraînement $ \ {x \ _n, t \ _n \} _ {n = 1} ^ N $, la fonction de coût à optimiser est

\begin{align}
E(w) &= \sum_{n=1}^N E_n\\
&= -\sum_{n=1}^N \ln p(t_n|x_n)\\
&= -\sum_{n=1}^N \ln \left\{\sum_{k=1}^K\pi_k(x_n,w)\mathcal{N}\left(t_n|\mu_k(x_n,w),\sigma^2_k(x_n,w)\right)\right\}
\end{align}

Ce sera. Ensuite, nous estimons le paramètre $ w $ qui minimise l'équation ci-dessus. Un réseau de neurones est utilisé pour les fonctions $ \ pi_k, \ mu_k, \ sigma_k $. Si vous placez l'entrée $ x $ dans le réseau avec le paramètre $ w $, le coefficient de mélange, la moyenne et la variance sont sortis et la distribution de probabilité de $ t $ est calculée.

Structure du réseau

Cette fois, l'entrée $ {\ bf x} $ est unidimensionnelle car la régression est effectuée sur les données comme indiqué dans la figure ci-dessus. Le nombre de nœuds au milieu $ {\ bf z} $ est de 5 selon PRML, et le nombre de nœuds dans la sortie $ {\ bf y} $ est de 9. Les paramètres de la première couche sont $ W ^ {(1)}, {\ bf b} ^ {(1)} $, la fonction d'activation est $ f (\ cdot) $, et les paramètres de la deuxième couche sont $ W ^ { Comme (2)}, {\ bf b} ^ {(2)} $, la propagation vers l'avant est la suivante.

\begin{align}
{\bf z} &= f(W^{(1)}{\bf x} + {\bf b}^{(1)})\\
{\bf y} &= W^{(2)}{\bf z} + {\bf b}^{(2)}
\end{align}

Ces neuf sorties sont les activités des paramètres de distribution gaussienne mixte $ \ pi_k, \ mu_k, \ sigma_k $. Premièrement, la sortie $ {\ bf y} $ est l'activité correspondant à $ \ pi, \ mu, \ sigma $ trois à partir du haut.

{\bf y} =
\begin{bmatrix}
a_{\pi_1}\\
a_{\pi_2}\\
a_{\pi_3}\\
a_{\mu_1}\\
a_{\mu_2}\\
a_{\mu_3}\\
a_{\sigma_1}\\
a_{\sigma_2}\\
a_{\sigma_3}
\end{bmatrix}

Puisque le coefficient de mélange $ \ pi_k $ a une contrainte de $ \ sum_k \ pi_k = 1 $, le coefficient de mélange est émis par la fonction softmax.

\pi_k = {\exp(a_{\pi_k})\over\sum_{l=1}^K\exp(a_{\pi_l})}

La moyenne est $ \ mu_k = a_ {\ mu_k} $ sans la transformation non linéaire. Puisque l'écart type $ \ sigma $ est égal ou supérieur à 0, utilisez la fonction exponentielle $ \ sigma_k = \ exp (a_ {\ sigma_k}) $. Vous avez maintenant les fonctions $ \ pi_k, \ mu_k, \ sigma_k $ nécessaires pour calculer la distribution gaussienne mixte.

Pente

Comme je l'ai écrit dans le précédent Implementation of Neural Network, pour former un réseau de neurones, un gradient obtenu en différenciant la fonction de coût avec la sortie du réseau de neurones est nécessaire. Avec cela, vous pouvez également utiliser la rétropropagation d'erreur pour calculer le différentiel pour le paramètre de fonction de coût $ w $. Le gradient en sortie du réseau neuronal de cette fonction de coût est

\gamma_{nk}(t_n|x_n) = {\pi_k\mathcal{N}(t_n|\mu_k(x_n,w),\sigma_k^2(x_n,w))\over\sum_l\pi_l\mathcal{N}(t_n|\mu_l(x_n,w),\sigma_l^2(x_n,w))}

En utilisant,

\begin{align}
{\partial E_n\over\partial a_{\pi_k}} &= \pi_k - \gamma_{nk}\\
{\partial E_n\over\partial a_{\mu_k}} &= \gamma_{nk}{\mu_k - t_n\over\sigma^2_k}\\
{\partial E_n\over\partial a_{\sigma_k}} &= \gamma_{nk}\left(1 - {||t_n - \mu_k||^2\over \sigma^2_k}\right)
\end{align}

Ce sera. La dérivation de formules est un exercice dans PRML, veuillez donc vous référer à la réponse ou à l'article sur le réseau à densité mixte.

la mise en oeuvre

Fonction de coût

Le code de la formule écrite ci-dessus est le suivant.

#Classe de fonction de coût
class GaussianMixture(object):
    """Negative Log Likelihood of Gaussian Mixture model"""
    def __init__(self, n_components):
        #Nombre de distributions gaussiennes, 3 dans les exemples précédents
        self.n_components = n_components

    #Méthode de calcul de la valeur de la fonction de coût
    def __call__(self, X, targets):

        #Convertir la sortie réseau X avec la fonction d'activation pour calculer l'écart type, le coefficient de mélange, la moyenne
        sigma, weight, mu = self.activate(X)

        #Fonction gaussienne N(t|mu,sigma^2)Calculez la valeur de
        gauss = self.gauss(mu, sigma, targets)

        #Log de vraisemblance négatif E pour une distribution gaussienne mixte(w):PRML(Équation 5.153)
        return -np.sum(np.log(np.sum(weight * gauss, axis=1)))

    #Convertir avec la fonction d'activation
    def activate(self, X):
        assert np.size(X, 1) == 3 * self.n_components

        #Divisez X où il correspond à l'écart type, au coefficient de mélange et à la moyenne
        X_sigma, X_weight, X_mu = np.split(X, [self.n_components, 2 * self.n_components], axis=1)

        #Convertir l'écart type avec la fonction d'activation
        sigma = np.exp(X_sigma)

        #Convertissez le coefficient de mélange avec la fonction d'activation et soustrayez la valeur maximale pour que les chiffres ne débordent pas
        weight = np.exp(X_weight - np.max(X_weight, 1, keepdims=True))
        weight /= np.sum(weight, axis=1, keepdims=True)

        return sigma, weight, X_mu

    #Fonction gaussienne N(target|mu,sigma^2)Calculer
    def gauss(self, mu, sigma, targets):
        return np.exp(-0.5 * (mu - targets) ** 2 / np.square(sigma)) / np.sqrt(2 * np.pi * np.square(sigma))

    #Différencier la fonction de coût par activité
    def delta(self, X, targets):
        sigma, weight, mu = self.activate(X)
        var = np.square(sigma)
        gamma = weight * self.gauss(mu, sigma, targets)
        gamma /= np.sum(gamma, axis=1, keepdims=True)

        #Calculez chaque différentiel
        delta_mu = gamma * (mu - targets) / var
        delta_sigma = gamma * (1 - (mu - targets) ** 2 / var)
        delta_weight = weight - gamma

        #Connectez-vous puis revenez
        delta = np.hstack([delta_sigma, delta_weight, delta_mu])
        return delta

Création de données d'entraînement

#Créer des données d'entraînement qui suivent la fonction inverse de la fonction func
def create_toy_dataset(func, n=300):
    t = np.random.uniform(size=(n, 1))
    x = func(t) + np.random.uniform(-0.05, 0.05, size=(n, 1))
    return x, t

Traitement par lots mini

L'une des méthodes que j'ai utilisées à la hâte car il y avait des moments où l'apprentissage ne se passait pas bien. Au lieu d'utiliser toutes les données d'entraînement pour calculer le gradient, nous en prenons certaines des données d'entraînement et effectuons le calcul du gradient. En faisant cela, vous pourrez peut-être sortir de la solution locale.

#Ensemble de données d'entraînement x,Sélectionnez aléatoirement n paires de t
def sample(x, t, n=None):
    assert len(x) == len(t)
    N = len(x)
    if n is None:
        n = N
    indices = np.random.choice(N, n, replace=False)
    return x[indices], t[indices]

Apprentissage des réseaux à densité mixte

C'est la fonction principale qui apprend le réseau à densité mixte et illustre le résultat.

def main():

    def func(x):
        return x + 0.3 * np.sin(2 * np.pi * x)

    #Générer des points de données selon la fonction inverse de la fonction func
    x, t = create_toy_dataset(func)

    #Déterminez la structure du réseau neuronal.
    layers = [TanhLayer(1, 5, std=0.1), LinearLayer(5, 9, std=0.1)]
    #Déterminer la fonction de coût à optimiser
    cost_function = GaussianMixture(3)
    #Créer un réseau de neurones
    nn = NeuralNetwork(layers, cost_function)
    #Décommentez la ligne ci-dessous pour vérifier si l'implémentation fonctionne.
    # nn._gradient_check(np.array([[0.5]]), np.array([[0.5]]))

    #Définir le premier coefficient d'apprentissage
    learning_rate = 1e-4
    for i in xrange(500000):
        #Calculez la valeur de la fonction de coût une fois toutes les 10000 fois et définissez le coefficient d'apprentissage sur 0.9 fois
        if i % 10000 == 0:
            print "step %6d, cost %f" % (i, nn.cost(x, t))
            learning_rate *= 0.9
        #Traitement par lots mini
        batch = sample(x, t, n=100)
        nn.fit(*batch, learning_rate=learning_rate)

    #Créer des données de test
    x_test = np.linspace(x.min(), x.max(), 100)
    y_test = np.linspace(t.min(), t.max(), 100)
    X_test, Y_test = np.meshgrid(x_test, y_test)
    test = np.array([X_test, Y_test]).transpose(1, 2, 0).reshape(-1, 2)

    #Calculer la probabilité des données de test
    sigma, weight, mu = nn(test[:, 0].reshape(-1, 1))
    probs = cost_function.gauss(mu, sigma, test[:, 1].reshape(-1, 1))
    probs = np.sum(weight * probs, axis=1)
    Probs = probs.reshape(100, 100)

    #PRML Figure 5.21(a)Reproduction de
    plt.plot(x_test, weight[:100, 0], color="blue")
    plt.plot(x_test, weight[:100, 1], color="red")
    plt.plot(x_test, weight[:100, 2], color="green")
    plt.title("weights")
    plt.show()

    #PRML Figure 5.21(b)Reproduction de
    plt.plot(x_test, mu[:100, 0], color="blue")
    plt.plot(x_test, mu[:100, 1], color="red")
    plt.plot(x_test, mu[:100, 2], color="green")
    plt.title("means")
    plt.show()

    #PRML Figure 5.21(c)Reproduction de
    plt.scatter(x, t, alpha=0.5, label="observation")
    levels_log = np.linspace(0, np.log(probs.max()), 21)
    levels = np.exp(levels_log)
    levels[0] = 0
    plt.contourf(X_test, Y_test, Probs, levels, alpha=0.5)
    plt.colorbar()
    plt.xlim(x.min(), x.max())
    plt.ylim(t.min(), t.max())
    plt.show()

Code entier

Le neural_network sur la troisième ligne est le neural_network.py implémenté la dernière fois. Mettez-le dans le même répertoire que ce code.

mixture_density_network.py


import matplotlib.pyplot as plt
import numpy as np
from neural_network import TanhLayer, LinearLayer, NeuralNetwork


class GaussianMixture(object):
    """Negative Log Likelihood of Gaussian Mixture model"""
    def __init__(self, n_components):
        self.n_components = n_components

    def __call__(self, X, targets):
        sigma, weight, mu = self.activate(X)
        gauss = self.gauss(mu, sigma, targets)
        return -np.sum(np.log(np.sum(weight * gauss, axis=1)))

    def activate(self, X):
        assert np.size(X, 1) == 3 * self.n_components
        X_sigma, X_weight, X_mu = np.split(X, [self.n_components, 2 * self.n_components], axis=1)
        sigma = np.exp(X_sigma)
        weight = np.exp(X_weight - np.max(X_weight, 1, keepdims=True))
        weight /= np.sum(weight, axis=1, keepdims=True)
        return sigma, weight, X_mu

    def gauss(self, mu, sigma, targets):
        return np.exp(-0.5 * (mu - targets) ** 2 / np.square(sigma)) / np.sqrt(2 * np.pi * np.square(sigma))

    def delta(self, X, targets):
        sigma, weight, mu = self.activate(X)
        var = np.square(sigma)
        gamma = weight * self.gauss(mu, sigma, targets)
        gamma /= np.sum(gamma, axis=1, keepdims=True)

        delta_mu = gamma * (mu - targets) / var
        delta_sigma = gamma * (1 - (mu - targets) ** 2 / var)
        delta_weight = weight - gamma
        delta = np.hstack([delta_sigma, delta_weight, delta_mu])
        return delta


def create_toy_dataset(func, n=300):
    t = np.random.uniform(size=(n, 1))
    x = func(t) + np.random.uniform(-0.05, 0.05, size=(n, 1))
    return x, t


def sample(x, t, n=None):
    assert len(x) == len(t)
    N = len(x)
    if n is None:
        n = N
    indices = np.random.choice(N, n, replace=False)
    return x[indices], t[indices]


def main():

    def func(x):
        return x + 0.3 * np.sin(2 * np.pi * x)

    x, t = create_toy_dataset(func)

    layers = [TanhLayer(1, 5, std=0.1), LinearLayer(5, 9, std=0.1)]
    cost_function = GaussianMixture(3)
    nn = NeuralNetwork(layers, cost_function)
    # nn._gradient_check(np.array([[0.5]]), np.array([[0.5]]))
    learning_rate = 1e-4
    for i in xrange(500000):
        if i % 10000 == 0:
            print "step %6d, cost %f" % (i, nn.cost(x, t))
            learning_rate *= 0.9
        batch = sample(x, t, n=100)
        nn.fit(*batch, learning_rate=learning_rate)

    x_test = np.linspace(x.min(), x.max(), 100)
    y_test = np.linspace(t.min(), t.max(), 100)
    X_test, Y_test = np.meshgrid(x_test, y_test)
    test = np.array([X_test, Y_test]).transpose(1, 2, 0).reshape(-1, 2)

    sigma, weight, mu = nn(test[:, 0].reshape(-1, 1))
    probs = cost_function.gauss(mu, sigma, test[:, 1].reshape(-1, 1))
    probs = np.sum(weight * probs, axis=1)
    Probs = probs.reshape(100, 100)

    plt.plot(x_test, weight[:100, 0], color="blue")
    plt.plot(x_test, weight[:100, 1], color="red")
    plt.plot(x_test, weight[:100, 2], color="green")
    plt.title("weights")
    plt.show()

    plt.plot(x_test, mu[:100, 0], color="blue")
    plt.plot(x_test, mu[:100, 1], color="red")
    plt.plot(x_test, mu[:100, 2], color="green")
    plt.title("means")
    plt.show()

    plt.scatter(x, t, alpha=0.5, label="observation")
    levels_log = np.linspace(0, np.log(probs.max()), 21)
    levels = np.exp(levels_log)
    levels[0] = 0
    plt.contourf(X_test, Y_test, Probs, levels, alpha=0.5)
    plt.colorbar()
    plt.xlim(x.min(), x.max())
    plt.ylim(t.min(), t.max())
    plt.show()


if __name__ == '__main__':
    main()

résultat

À la suite de l'exécution du code ci-dessus pour former le réseau à densité mixte en utilisant le point bleu comme point de données d'apprentissage, le graphique illustré à la figure 5.21 de PRML a été reproduit. Cependant, la couleur de la ligne, etc. peut être différente. Il est ajusté aux points de données d'apprentissage par trois distributions gaussiennes, bleu, rouge et vert. result_weights.png result_means.png result_distributions.png

État de l'apprentissage

Comme la dernière fois, j'ai fait une vidéo du processus d'apprentissage. anime_gaussian_mixture.gif

À la fin

Je n'avais pas entendu parler du réseau à densité mixte avant de le lire sur PRML, donc quand j'ai cherché sur le net, l'article écrit par l'auteur de PRML CM Bishop en 1994 était un succès, donc cette personne lui-même Peut avoir été préconisé par. Le réseau à densité mixte semble être plus efficace que le modèle conventionnel s'il est multimodal, je voudrais donc l'appliquer à certaines données réelles.

Recommended Posts

PRML Chapter 5 Implémentation Python de réseau à densité mixte
PRML Chapitre 5 Implémentation Python du réseau neuronal
PRML Chapitre 9 Implémentation Python de distribution gaussienne mixte
PRML Chapitre 14 Implémentation Python de modèle mixte conditionnel
PRML Chapitre 3 Preuve Implémentation approximative de Python
PRML Chapitre 8 Algorithme Somme des produits Implémentation Python
PRML Chapitre 4 Implémentation Python de la régression logistique bayésienne
PRML Chapitre 10 Implémentation Python de distribution gaussienne mixte
PRML Chapitre 6 Implémentation Python Gaussian Return
PRML Chapter 2 Student t-Distribution Python Implementation
PRML Chapitre 1 Implémentation de Python pour l'ajustement de courbe bayésienne
PRML Chapitre 11 Implémentation Python Monte Carlo Chaîne de Markov
PRML Chapitre 12 Mise en œuvre de l'analyse principale bayésienne Python
Implémentation Python Distribution mixte de Bernoulli
Implémentation de réseau neuronal en python
PRML Chapter 7 Implémentation de Python Vector Machine associée pour les problèmes de régression
Explication et mise en œuvre de PRML Chapitre 4
Implémentation de distribution normale mixte en python
PRML Chapter 13 Estimation la plus probable Implémentation Python du modèle de Markov caché
Implémentation PRML Chapitre 3 Modèle de fonction de base linéaire
Implémenté en Python PRML Chapitre 7 SVM non linéaire
Implémenté dans Python PRML Chapter 5 Neural Network
Implémenté en Python PRML Chapitre 1 Estimation bayésienne
Implémenté en Python PRML Chapitre 3 Régression linéaire bayésienne
Python vs Ruby "Deep Learning from scratch" Chapitre 3 Implémentation d'un réseau neuronal à 3 couches
Implémenté dans Python PRML Chapitre 1 Ajustement de courbe polygonale
[Python] Implémentation du clustering à l'aide d'un modèle gaussien mixte
Implémenté en Python PRML Chapitre 4 Classification par algorithme Perceptron
Implémentation RNN en python
Implémentation ValueObject en Python
[Python] Chapitre 01-01 À propos de Python (First Python)
Implémentation SVM en python