PRML Chapitre 6 Implémentation Python Gaussian Return

Prévoyez d'implémenter l'algorithme PRML presque exclusivement avec Numpy, mais la version japonaise de PRML est déjà entrée dans le deuxième volume. Dans la méthode du noyau du chapitre 6, le produit interne des points de données dans l'espace des caractéristiques est calculé à l'aide de la fonction noyau (astuce du noyau), et une régression linéaire ou une séparation linéaire est effectuée dans cet espace. Je ne connais pas grand-chose au processus gaussien, donc je pense qu'il y a beaucoup d'erreurs. Le code Ichiou a une implémentation qui fonctionne comme ça.

Processus gaussien

Dans la régression linéaire sur laquelle nous avons travaillé au chapitre 3, etc., le paramètre $ {\ bf w} $ suivait souvent la distribution gaussienne. Donc $ y = {\ bf w} ^ {\ rm T} \ phi (x) $ suit également une distribution gaussienne unidimensionnelle, et $ {\ bf y} = {\ bf \ Phi} {\ bf w} $ Il suit une distribution gaussienne multidimensionnelle. Cependant, $ {\ bf \ Phi} $ est une matrice de planification de $ \ {x_1, \ dots, x_N \} $. Dans de tels cas, on dit que $ p ({\ bf y}) $ suit le processus gaussien. ** Le processus gaussien suit une distribution gaussienne non seulement en une dimension ($ y $ est un scalaire), en plusieurs dimensions ($ y $ est un nombre fini de vecteurs), mais aussi en dimensions infinies ($ y $ est un nombre infini de vecteurs). Elle peut être interprétée comme une distribution gaussienne de dimension infinie **.

La distribution simultanée de $ {\ bf y} $ est une distribution gaussienne avec la moyenne de $ {\ bf 0} $ et la covariance de la matrice gramme $ K $.

p({\bf y}) = \mathcal{N}({\bf y}|{\bf 0},K)

Sera. Cependant, les éléments de la matrice gramme sont $ K_ {nm} = k (x_n, x_m) $ avec $ k (x, x ') $ comme fonction de noyau. Les fonctions gaussiennes $ k (x, x ') = a \ exp \ left (-b (x --x') ^ 2 \ right) $ sont souvent utilisées comme fonctions du noyau, avec $ a et b $ comme constantes. Je vais. Pour $ x_n, x_m $ qui sont similaires les uns aux autres, la corrélation entre $ y (x_n) et y (x_m) $ est définie comme élevée. Cette «similitude» dépend du type de fonction du noyau que vous utilisez.

Retour par procédé gaussien

Soit la cible observée $ t_n = y_n + \ epsilon_n $. Ici, $ y_n = {\ bf w} ^ {\ rm T} \ phi (x_n) $ et le bruit $ \ epsilon_n $ sont des bruits qui suivent la distribution gaussienne ajoutée à la nième valeur d'observation. Avec le paramètre de précision $ \ beta $

p(t_n|y_n) = \mathcal{N}(t_n|y_n,\beta^{(-1)})

Ce sera. Par conséquent, comme $ {\ bf t} = (t_1, \ dots, t_N) ^ {\ rm T} $,

p({\bf t}|{\bf y}) = \mathcal{N}({\bf t}|{\bf y},\beta^{(-1)}{\bf I}_N)

Ce sera. Maintenant que la probabilité simultanée de $ {\ bf y} $ a déjà été définie ci-dessus, la distribution simultanée de $ {\ bf t} $ peut être obtenue.

\begin{align}
p({\bf t}) &= \int p({\bf t}|{\bf y})p({\bf y})d{\bf y}\\
&= \int \mathcal{N}({\bf t}|{\bf y},\beta^{(-1)}{\bf I}_N)\mathcal{N}({\bf y}|{\bf 0},K)d{\bf y}\\
&= \mathcal{N}({\bf t}|{\bf 0},{\bf C}_N)
\end{align}

Cependant, $ {\ bf C} \ _N = K + \ beta ^ {(-1)} {\ bf I} \ _N $.

La probabilité simultanée de $ {\ bf t}, t_ {N + 1} $ est supérieure à $ t_ {N + 1} $ d'être nouvellement prédite en plus de N observations $ {\ bf t} $. De la discussion

p({\bf t},t_{N+1}) = \mathcal{N}({\bf t},t_{N+1}|{\bf 0}, {\bf C}_{N+1})

Cependant, la matrice de covariance est $ {\ bf k} = \ left (k (x_1, x_ {N + 1}), \ dots, k (x_N, x_ {N + 1}) \ right), c = k (x_ {N + 1}, x_ {N + 1}) comme $

{\bf C}_{N+1} =
\begin{bmatrix}
{\bf C}_N & {\bf k}\\
{\bf k}^{\rm T} & c
\end{bmatrix}

Ce sera. Maintenant que nous connaissons la distribution simultanée de $ {\ bf t}, t_ {N + 1} $, nous pouvons trouver la distribution conditionnelle $ p (t_ {N + 1} | {\ bf t}) $.

p(t_{N+1}|{\bf t}) = \mathcal{N}(t_{N+1}|{\bf k}^{\rm T}{\bf C}_N^{-1}{\bf t},c-{\bf k}^{\rm T}{\bf C}_N^{-1}{\bf k})

Apprentissage des super paramètres

La discussion ci-dessus a corrigé la fonction du noyau $ k (x, x ') $. Par exemple, dans la fonction noyau, $ k (x, x ') = a \ exp (-b (x-x') ^ 2) $, $ a, b $ a été défini comme une constante. Ensuite, que devons-nous faire des constantes utilisées pour la fonction noyau? Cette estimation de paramètre correspond à l'estimation de super paramètre implémentée dans PRML Chapitre 3. À ce moment-là, la fonction de preuve $ p ({\ bf t} | \ theta) $ a été maximisée. $ \ Theta $ est un paramètre de fonction du noyau, par exemple $ \ theta = (a, b) ^ {\ rm T} $. La fonction de preuve logarithmique est

\ln p({\bf t}|\theta) = -{1\over2}\ln |{\bf C}_N| - {1\over2}{\bf t}^{\rm T}{\bf C}_N^{-1}{\bf t} - {N\over2}\ln(2\pi)

Et quand cela est différencié par rapport aux paramètres,

{\partial\over\partial\theta_i}\ln p({\bf t}|\theta) = -{1\over2}{\rm Tr}({\bf C}_N^{-1}{\partial{\bf C}_N\over\partial\theta_i}) + {1\over2}{\bf t}^{\rm T}{\bf C}_N^{-1}{\partial{\bf C}_N\over\partial\theta_i}{\bf C}_N^{-1}{\bf t}

Ce sera. J'ai mis à jour les paramètres en utilisant la méthode du gradient.

la mise en oeuvre

Bibliothèque

import matplotlib.pyplot as plt
import numpy as np

Il existe deux bibliothèques à importer: matplotlib, qui est une bibliothèque pour dessiner des graphiques, et numpy, qui effectue des calculs matriciels.

Fonction noyau

Cette fois, j'ai utilisé $ a \ exp (- {b \ over2} (x-x ') ^ 2) $ comme fonction du noyau.

#Classe de fonction du noyau
class GaussianKernel(object):
    #Paramètres de la fonction noyau a,Initialiser b
    def __init__(self, params):
        assert np.shape(params) == (2,)
        self.__params = params

    #Paramètres de la fonction noyau a,Renvoie b
    def get_params(self):
        return np.copy(self.__params)

     # x,Calculer la valeur de la fonction noyau avec y comme expression PRML d'entrée(6.63)
    def __call__(self, x, y):
        return self.__params[0] * np.exp(-0.5 * self.__params[1] * (x - y) ** 2)

    # x,Calculer le différentiel avec les paramètres de la fonction noyau lorsque y est entré
    def derivatives(self, x, y):
        sq_diff = (x - y) ** 2
        #Différenciation avec le paramètre a
        delta_0 = np.exp(-0.5 * self.__params[1] * sq_diff)
        #Différenciation avec le paramètre b
        delta_1 = -0.5 * sq_diff * delta_0 * self.__params[0]
        return (delta_0, delta_1)

    #Paramètres de fonction du noyau mis à jour
    def update_parameters(self, updates):
        assert np.shape(updates) == (2,)
        self.__params += updates

Retour par procédé gaussien

#Classe qui effectue une régression par processus gaussien
class GaussianProcessRegression(object):
    #Initialisation des fonctions du noyau et des paramètres de précision du bruit
    def __init__(self, kernel, beta=1.):
        self.kernel = kernel
        self.beta = beta

    #Régression sans estimation des paramètres de la fonction noyau
    def fit(self, x, t):
        self.x = x
        self.t = t
        #Calcul de la matrice de Gram Formule PRML(6.54)
        Gram = self.kernel(*np.meshgrid(x, x))
        #Calcul de la formule PRML de la matrice de covariance(6.62)
        self.covariance = Gram + np.identity(len(x)) / self.beta
        #Calcul de la matrice de précision
        self.precision = np.linalg.inv(self.covariance)

    #Régression pour estimer les paramètres de la fonction noyau
    def fit_kernel(self, x, t, learning_rate=0.1, iter_max=10000):
        for i in xrange(iter_max):
            params = self.kernel.get_params()
            #Retour avec les paramètres courants de la fonction noyau
            self.fit(x, t)
            #Différencier la fonction de preuve logarithmique avec des paramètres
            gradients = self.kernel.derivatives(*np.meshgrid(x, x))
            #Calculer la quantité de mises à jour des paramètres Formule PRML(6.70)
            updates = np.array(
                [-np.trace(self.precision.dot(grad)) + t.dot(self.precision.dot(grad).dot(self.precision).dot(t)) for grad in gradients])
            #Mettre à jour les paramètres
            self.kernel.update_parameters(learning_rate * updates)
            #Arrêter la mise à jour si le nombre de mises à jour des paramètres est faible
            if np.allclose(params, self.kernel.get_params()):
                break
        else:
            #Si le nombre de mises à jour de paramètres n'est pas faible même après la mise à jour du nombre de mises à jour par défaut, l'instruction suivante est générée.
            print "parameters may not have converged"

    #Distribution prévue en sortie
    def predict_dist(self, x):
        K = self.kernel(*np.meshgrid(x, self.x, indexing='ij'))
        #Calculer la moyenne de la formule PRML de distribution prédite(6.66)
        mean = K.dot(self.precision).dot(self.t)
        #Calculer la variance de la formule PRML de distribution prédite(6.67)
        var = self.kernel(x, x) + 1 / self.beta - np.sum(K.dot(self.precision) * K, axis=1)
        return mean.ravel(), np.sqrt(var.ravel())

Fonction principale

Le remplacement de regression.fit_kernel (...) par regression.fit (x, t) provoquera la régression sans estimer les paramètres.

def main():
    #Définissez la fonction que suivent les données d'entraînement
    def func(x):
        return np.sin(2 * np.pi * x)
    #Générer des données d'entraînement
    x, t = create_toy_data(func, high=0.7, std=0.1)

    #Les réglages et paramètres de la fonction du noyau sont enfin définis
    kernel = GaussianKernel(params=np.array([1., 1.]))
    #Fonctions du noyau et paramètres de précision utilisés pour la régression de processus gaussien(Vraie valeur)paramètres de
    regression = GaussianProcessRegression(kernel=kernel, beta=100)
    #Régression, y compris l'estimation des paramètres des fonctions du noyau.fit(x,t)Si vous passez à, revenez sans estimer les paramètres
    regression.fit_kernel(x, t, learning_rate=0.1, iter_max=10000)

    #Distribution de prédiction de sortie pour les données de test
    x_test = np.linspace(0, 1, 100)
    y, y_std = regression.predict_dist(x_test)

    #Tracer les résultats de la régression
    plt.scatter(x, t, alpha=0.5, color="blue", label="observation")
    plt.plot(x_test, func(x_test), color="blue", label="sin$(2\pi x)$")
    plt.plot(x_test, y, color="red", label="predict_mean")
    plt.fill_between(x_test, y - y_std, y + y_std, color="pink", alpha=0.5, label="predict_std")
    plt.legend(loc="lower left")
    plt.show()

Code entier

gaussian_process_regression.py

import matplotlib.pyplot as plt
import numpy as np


class GaussianKernel(object):

    def __init__(self, params):
        assert np.shape(params) == (2,)
        self.__params = params

    def get_params(self):
        return np.copy(self.__params)

    def __call__(self, x, y):
        return self.__params[0] * np.exp(-0.5 * self.__params[1] * (x - y) ** 2)

    def derivatives(self, x, y):
        sq_diff = (x - y) ** 2
        delta_0 = np.exp(-0.5 * self.__params[1] * sq_diff)
        delta_1 = -0.5 * sq_diff * delta_0 * self.__params[0]
        return (delta_0, delta_1)

    def update_parameters(self, updates):
        assert np.shape(updates) == (2,)
        self.__params += updates


class GaussianProcessRegression(object):

    def __init__(self, kernel, beta=1.):
        self.kernel = kernel
        self.beta = beta

    def fit(self, x, t):
        self.x = x
        self.t = t
        Gram = self.kernel(*np.meshgrid(x, x))
        self.covariance = Gram + np.identity(len(x)) / self.beta
        self.precision = np.linalg.inv(self.covariance)

    def fit_kernel(self, x, t, learning_rate=0.1, iter_max=10000):
        for i in xrange(iter_max):
            params = self.kernel.get_params()
            self.fit(x, t)
            gradients = self.kernel.derivatives(*np.meshgrid(x, x))
            updates = np.array(
                [-np.trace(self.precision.dot(grad)) + t.dot(self.precision.dot(grad).dot(self.precision).dot(t)) for grad in gradients])
            self.kernel.update_parameters(learning_rate * updates)
            if np.allclose(params, self.kernel.get_params()):
                break
        else:
            print "parameters may not have converged"

    def predict_dist(self, x):
        K = self.kernel(*np.meshgrid(x, self.x, indexing='ij'))
        mean = K.dot(self.precision).dot(self.t)
        var = self.kernel(x, x) + 1 / self.beta - np.sum(K.dot(self.precision) * K, axis=1)
        return mean.ravel(), np.sqrt(var.ravel())


def create_toy_data(func, low=0, high=1., n=10, std=1.):
    x = np.random.uniform(low, high, n)
    t = func(x) + np.random.normal(scale=std, size=n)
    return x, t


def main():

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

    x, t = create_toy_data(func, high=0.7, std=0.1)

    kernel = GaussianKernel(params=np.array([1., 1.]))
    regression = GaussianProcessRegression(kernel=kernel, beta=100)
    regression.fit_kernel(x, t, learning_rate=0.1, iter_max=10000)

    x_test = np.linspace(0, 1, 100)
    y, y_std = regression.predict_dist(x_test)

    plt.scatter(x, t, alpha=0.5, color="blue", label="observation")
    plt.plot(x_test, func(x_test), color="blue", label="sin$(2\pi x)$")
    plt.plot(x_test, y, color="red", label="predict_mean")
    plt.fill_between(x_test, y - y_std, y + y_std, color="pink", alpha=0.5, label="predict_std")
    plt.legend(loc="lower left")
    plt.show()

if __name__ == '__main__':
    main()

résultat

Vous avez maintenant un diagramme comme la figure 6.8 dans PRML. La variance de la distribution prédite est grande dans la région sans données d'apprentissage et petite dans la région avec données d'apprentissage.

À propos, le résultat lorsque l'estimation la plus probable du super paramètre n'est pas effectuée avec les mêmes données d'apprentissage que ci-dessus est comme indiqué sur la figure ci-dessous. La distribution prédite ci-dessus est plus appropriée pour les données d'entraînement (également pour les données de test).

À la fin

Après tout, ce que vous faites est similaire à la régression linéaire bayésienne. Cependant, au lieu d'utiliser explicitement le vecteur de caractéristiques, j'ai utilisé une astuce du noyau pour retourner en utilisant uniquement le produit interne du vecteur de caractéristiques. Bien que non mentionné dans le PRML, le processus gaussien est utilisé dans l'optimisation bayésienne, une technique qui semble attirer l'attention récemment. Il semble que l'optimisation bayésienne soit utilisée pour déterminer des paramètres tels que le coefficient d'apprentissage du réseau neuronal. J'aimerais mettre en œuvre l'optimisation bayésienne si j'en ai l'occasion.