Planen Sie die Implementierung des PRML-Algorithmus fast ausschließlich mit Numpy, aber die japanische Version von PRML hat bereits den zweiten Band eingegeben. Bei der Kernelmethode von Kapitel 6 wird das innere Produkt von Datenpunkten im Merkmalsraum unter Verwendung der Kernelfunktion (Kerneltrick) berechnet, und in diesem Raum wird eine lineare Regression oder lineare Trennung durchgeführt. Ich weiß nicht viel über den Gaußschen Prozess, daher denke ich, dass es viele Fehler gibt. Ichiou-Code hat eine Implementierung, die so funktioniert.

Gaußscher Prozess

In der linearen Regression, an der wir in Kapitel 3 usw. gearbeitet haben, folgte der Parameter $ {\ bf w} $ häufig der Gaußschen Verteilung. Also folgt $ y = {\ bf w} ^ {\ rm T} \ phi (x) $ auch einer eindimensionalen Gaußschen Verteilung und $ {\ bf y} = {\ bf \ Phi} {\ bf w} $ Es folgt eine mehrdimensionale Gaußsche Verteilung. $ {\ Bf \ Phi} $ ist jedoch eine Planungsmatrix von $ \ {x_1, \ dots, x_N \} $. In solchen Fällen soll $ p ({\ bf y}) $ dem Gaußschen Prozess folgen. ** Der Gaußsche Prozess folgt einer Gaußschen Verteilung nicht nur in einer Dimension ($ y $ ist ein Skalar), in mehreren Dimensionen ($ y $ ist eine endliche Anzahl von Vektoren), sondern auch in unendlichen Dimensionen ($ y $ ist eine unendliche Anzahl von Vektoren). Es kann als unendlich dimensionale Gaußsche Verteilung ** interpretiert werden.

Die gleichzeitige Verteilung von $ {\ bf y} $ ist eine Gaußsche Verteilung mit dem Mittelwert von $ {\ bf 0} $ und der Kovarianz der Gramm-Matrix $ K $.

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

Wird sein. Die Elemente der Gramm-Matrix sind jedoch $ K_ {nm} = k (x_n, x_m) $ mit $ k (x, x ') $ als Kernelfunktion. Gaußsche Funktionen $ k (x, x ') = a \ exp \ left (-b (x - x') ^ 2 \ right) $ werden häufig als Kernfunktionen verwendet, wobei $ a und b $ als Konstanten dienen. Ich werde. Für $ x_n, x_m $, die einander ähnlich sind, wird die Korrelation zwischen $ y (x_n) und y (x_m) $ hoch gesetzt. Diese "Ähnlichkeit" hängt davon ab, welche Art von Kernelfunktion Sie verwenden.

Rückkehr nach dem Gaußschen Verfahren

Das beobachtete Ziel sei $ t_n = y_n + \ epsilon_n $. Hier sind $ y_n = {\ bf w} ^ {\ rm T} \ phi (x_n) $ und Rauschen $ \ epsilon_n $ Rauschen, das der Gaußschen Verteilung folgt, die zum n-ten Beobachtungswert addiert wird. Mit dem Präzisionsparameter $ \ beta $

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

Es wird sein. Daher ist als $ {\ 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)

Es wird sein. Nachdem wir oben bereits die gleichzeitige Wahrscheinlichkeit von $ {\ bf y} $ definiert haben, können wir die gleichzeitige Verteilung von $ {\ bf t} $ finden.

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

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

Die gleichzeitige Wahrscheinlichkeit von $ {\ bf t}, t_ {N + 1} $ liegt über $ t_ {N + 1} $, die zusätzlich zu N Beobachtungen $ {\ bf t} $ neu vorhergesagt werden soll. Aus der Diskussion

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

Die Kovarianzmatrix ist jedoch $ {\ 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}) als $

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

Es wird sein. Nachdem wir nun die gleichzeitige Verteilung von $ {\ bf t}, t_ {N + 1} $ kennen, können wir die bedingte Verteilung $ p (t_ {N + 1} | {\ bf t}) $ finden.

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

Super Parameter lernen

Die obige Diskussion hat die Kernelfunktion $ k (x, x ') $ behoben. Zum Beispiel wurde in der Kernelfunktion $ k (x, x ') = a \ exp (-b (x-x') ^ 2) $, $ a, b $ als Konstante gesetzt. Was sollen wir dann mit den Konstanten tun, die für die Kernelfunktion verwendet werden? Diese Parameterschätzung entspricht der in PRML Kapitel 3 implementierten Superparameterschätzung. Zu diesem Zeitpunkt wurde die Evidenzfunktion $ p ({\ bf t} | \ theta) $ maximiert. $ \ Theta $ ist ein Kernelfunktionsparameter, zum Beispiel $ \ theta = (a, b) ^ {\ rm T} $. Die logarithmische Beweisfunktion ist

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

Und wenn dies in Bezug auf die Parameter unterschieden wird,

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

Es wird sein. Ich habe die Parameter mit der Gradientenmethode aktualisiert.

Implementierung

Bibliothek

import matplotlib.pyplot as plt
import numpy as np

Es müssen zwei Bibliotheken importiert werden: matplotlib, eine Bibliothek zum Zeichnen von Diagrammen, und numpy, die Matrixberechnungen durchführt.

Kernelfunktion

Diesmal habe ich $ a \ exp (- {b \ over2} (x-x ') ^ 2) $ als Kernelfunktion verwendet.

#Kernelfunktionsklasse
class GaussianKernel(object):
    #Kernelfunktionsparameter a,Initialisieren b
    def __init__(self, params):
        assert np.shape(params) == (2,)
        self.__params = params

    #Kernelfunktionsparameter a,Rückgabe b
    def get_params(self):
        return np.copy(self.__params)

     # x,Berechnen Sie den Wert der Kernelfunktion mit y als Eingabe-PRML-Ausdruck(6.63)
    def __call__(self, x, y):
        return self.__params[0] * np.exp(-0.5 * self.__params[1] * (x - y) ** 2)

    # x,Berechnen Sie das Differential mit den Parametern der Kernelfunktion, wenn y eingegeben wird
    def derivatives(self, x, y):
        sq_diff = (x - y) ** 2
        #Differenzierung mit Parameter a
        delta_0 = np.exp(-0.5 * self.__params[1] * sq_diff)
        #Differenzierung mit Parameter b
        delta_1 = -0.5 * sq_diff * delta_0 * self.__params[0]
        return (delta_0, delta_1)

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

Rückkehr nach dem Gaußschen Verfahren

#Klasse, die eine Regression nach dem Gaußschen Prozess durchführt
class GaussianProcessRegression(object):
    #Initialisierung von Kernelfunktionen und Rauschgenauigkeitsparametern
    def __init__(self, kernel, beta=1.):
        self.kernel = kernel
        self.beta = beta

    #Regression ohne Parameterschätzung der Kernelfunktion
    def fit(self, x, t):
        self.x = x
        self.t = t
        #Gramm Matrix Berechnung PRML Formel(6.54)
        Gram = self.kernel(*np.meshgrid(x, x))
        #Berechnung der PRML-Formel der Kovarianzmatrix(6.62)
        self.covariance = Gram + np.identity(len(x)) / self.beta
        #Berechnung der Präzisionsmatrix
        self.precision = np.linalg.inv(self.covariance)

    #Regression zum Schätzen von Kernelfunktionsparametern
    def fit_kernel(self, x, t, learning_rate=0.1, iter_max=10000):
        for i in xrange(iter_max):
            params = self.kernel.get_params()
            #Kehren Sie mit den aktuellen Parametern der Kernelfunktion zurück
            self.fit(x, t)
            #Unterscheiden Sie die logarithmische Evidenzfunktion mit Parametern
            gradients = self.kernel.derivatives(*np.meshgrid(x, x))
            #Berechnen Sie die Anzahl der Parameteraktualisierungen der PRML-Formel(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])
            #Parameter aktualisieren
            self.kernel.update_parameters(learning_rate * updates)
            #Beenden Sie die Aktualisierung, wenn die Anzahl der Parameteraktualisierungen gering ist
            if np.allclose(params, self.kernel.get_params()):
                break
        else:
            #Wenn der Parameteraktualisierungsbetrag auch nach dem Aktualisieren der Standardanzahl von Aktualisierungen nicht klein ist, wird die folgende Anweisung ausgegeben.
            print "parameters may not have converged"

    #Prognostizierte Verteilung ausgeben
    def predict_dist(self, x):
        K = self.kernel(*np.meshgrid(x, self.x, indexing='ij'))
        #Berechnen Sie den Mittelwert der PRML-Formel für die vorhergesagte Verteilung(6.66)
        mean = K.dot(self.precision).dot(self.t)
        #Berechnen Sie die Varianz der PRML-Formel für die vorhergesagte Verteilung(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())

Hauptfunktion

Das Ersetzen von "Regression.fit_kernel (...)" durch "Regression.fit (x, t)" bewirkt die Regression, ohne die Parameter zu schätzen.

def main():
    #Stellen Sie die Funktion ein, der die Trainingsdaten folgen
    def func(x):
        return np.sin(2 * np.pi * x)
    #Trainingsdaten generieren
    x, t = create_toy_data(func, high=0.7, std=0.1)

    #Die Einstellungen und Parameter der Kernelfunktion werden endgültig festgelegt
    kernel = GaussianKernel(params=np.array([1., 1.]))
    #Kernfunktionen und Genauigkeitsparameter, die für die Gaußsche Prozessregression verwendet werden(Wahrer Wert)Einstellungen von
    regression = GaussianProcessRegression(kernel=kernel, beta=100)
    #Regression, einschließlich Schätzung der Kernelfunktionsparameter.fit(x,t)Wenn Sie zu wechseln, kehren Sie zurück, ohne die Parameter zu schätzen
    regression.fit_kernel(x, t, learning_rate=0.1, iter_max=10000)

    #Prognostizierte Verteilung für Testdaten ausgeben
    x_test = np.linspace(0, 1, 100)
    y, y_std = regression.predict_dist(x_test)

    #Zeichnen Sie die Ergebnisse der Regression
    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()

Ganzer Code

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

Ergebnis

Sie haben jetzt ein Diagramm wie in Abbildung 6.8 in PRML. Die Varianz der vorhergesagten Verteilung ist in der Region ohne Trainingsdaten groß und in der Region mit Trainingsdaten klein.

Übrigens ist das Ergebnis, wenn die wahrscheinlichste Schätzung des Superparameters nicht mit den gleichen Trainingsdaten wie oben durchgeführt wird, wie in der folgenden Abbildung gezeigt. Die oben vorhergesagte Verteilung eignet sich besser für Trainingsdaten (auch für Testdaten).

Am Ende

Schließlich ähnelt das, was Sie tun, der Bayes'schen linearen Regression. Anstatt den Feature-Vektor explizit zu verwenden, habe ich einen Kernel-Trick verwendet, um nur das innere Produkt des Feature-Vektors zu verwenden. Obwohl in der PRML nicht erwähnt, wird der Gaußsche Prozess bei der Bayes'schen Optimierung verwendet, eine Technik, die in letzter Zeit an Aufmerksamkeit zu gewinnen scheint. Es scheint, dass die Bayes'sche Optimierung verwendet wird, um Parameter wie den Lernkoeffizienten des neuronalen Netzwerks zu bestimmen. Ich möchte die Bayes'sche Optimierung implementieren, wenn ich die Gelegenheit dazu habe.

PRML Kapitel 6 Gaussian Return Python-Implementierung

Gaußscher Prozess

Rückkehr nach dem Gaußschen Verfahren

Super Parameter lernen

Implementierung

Bibliothek

Kernelfunktion

Rückkehr nach dem Gaußschen Verfahren

Hauptfunktion

Ganzer Code

gaussian_process_regression.py

Ergebnis

Am Ende

`gaussian_process_regression.py`