PRML Kapitel 1 Bayesian Curve Fitting Python-Implementierung

Für meine eigene Studie habe ich einige der in PRML eingeführten Methoden aufgegriffen und beschlossen, sie in Python zu implementieren. Grundsätzlich werden wir nach der Regel codieren, dass nur numpy verwendet werden kann, außer der Standard-Python-Bibliothek (von der scipy verwendet werden kann). Matplotlib usw. sind jedoch Ausnahmen, da sie nichts mit dem Algorithmus-Teil zu tun haben.

Bayes-Kurvenanpassung

Die Bayes'sche Kurvenanpassung ist eine Methode zum Anpassen einer Kurve, indem sie vollständig wie die Bayes'sche Kurve behandelt wird. Während die vorhergesagte Verteilung durch die wahrscheinlichste Schätzung an allen Punkten die gleiche Varianz aufweist, verwendet die Bayes'sche Kurvenanpassung die posteriore Verteilung von Parametern, um die vorhergesagte Verteilung auf Bayes'sche Weise zu erhalten (unter Verwendung des Additions- und Multiplikationssatzes der Wahrscheinlichkeit). Sie können die Varianz für jeden Punkt mit berechnen. Es wird berechnet, ob Sie in Bezug auf jeden Punkt relativ sicher sind.

Implementierung

Wie eingangs erwähnt, verwendet die Algorithmusimplementierung numpy.

import matplotlib.pyplot as plt
import numpy as np

Berechnung der Entwurfsmatrix

class PolynomialFeatures(object):

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

    def transform(self, x):
        features = [x ** i for i in xrange(self.degree + 1)]
        return np.array(features).transpose()

Schätzung und Vorhersage der Bayes'schen Kurvenanpassung

Schätzen Sie zunächst die hintere Verteilung des Parameters $ {\ bf w} . $ p({\bf w}|{\bf x}, {\bf t}, \alpha, \beta) \propto p({\bf t}|{\bf\Phi}{\bf w}, \beta)p({\bf w}|\alpha) $$ Hier sind $ {\ bf x} und {\ bf t} $ Trainingssätze. Außerdem nimmt $ p $ auf der rechten Seite eine Gaußsche Funktion an, und die Genauigkeitsparameter sind $ \ alpha $ bzw. $ \ beta . Dann wird die Wahrscheinlichkeitsverteilung auf der linken Seite auch eine Gaußsche Verteilung. $ p({\bf w}|{\bf x}, {\bf t}, \alpha, \beta) = \mathcal{N}({\bf m}, {\bf S}) $ Allerdings $ {\bf m} = \beta{\bf S}{\bf\Phi}^{\rm T}{\bf t},~{\bf S}^{-1} = \alpha{\bf I} + \beta{\bf\Phi}^{\rm T}{\bf\Phi} $$ Unter Verwendung dieser posterioren Verteilung die Wahrscheinlichkeitsverteilung von $ t $ für den neuen Punkt $ x $ $ p(t|x,{\bf x}, {\bf t}, \alpha, \beta) = \int p(t|x,{\bf w},\beta)p({\bf w}|{\bf x}, {\bf t}, \alpha, \beta){\rm d}{\bf w} $ Dies ist auch eine Gaußsche Verteilung, wenn sie berechnet wird. $ p(t|x,{\bf x}, {\bf t}, \alpha, \beta) = \mathcal{N}(t|{\bf m}^{\rm T}{\bf\phi}(x), {1\over\beta}+{\bf\phi}(x)^{\rm T}{\bf S}{\bf\phi}(x)) $

class BayesianRegression(object):

    def __init__(self, alpha=0.1, beta=0.25):
        self.alpha = alpha
        self.beta = beta

    def fit(self, X, t):
        self.w_var = np.linalg.inv(
            self.alpha * np.identity(np.size(X, 1))
            + self.beta * X.T.dot(X))
        self.w_mean = self.beta * self.w_var.dot(X.T.dot(t))

    def predict(self, X):
        y = X.dot(self.w_mean)
        y_var = 1 / self.beta + np.sum(X.dot(self.w_var) * X, axis=1)
        y_std = np.sqrt(y_var)
        return y, y_std

Ganzer Code

Der folgende Code ist das, was ich implementiert habe, einschließlich des obigen. Der allgemeine Fluss ist

1. Trainingsdaten erstellen (x, t = create_toy_data (func, ...))
2. Berechnen Sie die posteriore Verteilung des Gewichtsparameters $ {\ bf w} $ (regression.fit (X, t))
3. Berechnen Sie die vorhergesagte Verteilung mit der posterioren Verteilung (y, y_std = Regression.predict (X_test))
4. Illustrieren Sie die Ergebnisse Es ist geworden.

bayesian_curve_fitting.py

import matplotlib.pyplot as plt
import numpy as np


class PolynomialFeatures(object):

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

    def transform(self, x):
        features = [x ** i for i in xrange(self.degree + 1)]
        return np.array(features).transpose()


class BayesianRegression(object):

    def __init__(self, alpha=0.1, beta=0.25):
        self.alpha = alpha
        self.beta = beta

    def fit(self, X, t):
        self.w_var = np.linalg.inv(
            self.alpha * np.identity(np.size(X, 1))
            + self.beta * X.T.dot(X))
        self.w_mean = self.beta * self.w_var.dot(X.T.dot(t))

    def predict(self, X):
        y = X.dot(self.w_mean)
        y_var = 1 / self.beta + np.sum(X.dot(self.w_var) * X, axis=1)
        y_std = np.sqrt(y_var)
        return y, y_std


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


def main():

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

    x, t = create_toy_data(func, low=0, high=1, size=10, sigma=0.25)
    plt.scatter(x, t, s=50, marker='o', alpha=0.5, label="observation")

    features = PolynomialFeatures(degree=5)

    regression = BayesianRegression(alpha=1e-3, beta=2)
    X = features.transform(x)
    regression.fit(X, t)

    x_test = np.linspace(0, 1, 100)
    X_test = features.transform(x_test)
    y, y_std = regression.predict(X_test)

    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()
    plt.title("Predictive distribution")
    plt.xlabel("x")
    plt.ylabel("t")
    plt.show()


if __name__ == '__main__':
    main()

Ergebnis

Es wurde vorhergesagt, indem eine Bayes'sche Kurvenanpassung aus den Datenpunkten (blauen Punkten) durchgeführt wurde, die durch Hinzufügen von Rauschen zur Funktion $ {\ rm sin} (2 \ pi x) $ erzeugt wurden. Der Durchschnittswert jedes $ x $ der vorhergesagten Verteilung wird in Rot angezeigt, und die Standardabweichung wird in Rosa angezeigt. Die Breite des rosa Bandes entspricht der Unsicherheit der Vorhersage. Wenn Sie die Vorhersageverteilung anhand der wahrscheinlichsten Schätzung berechnen, ist die Breite des rosa Bandes unabhängig von $ x $ konstant, und die Vorhersage eines beliebigen Punkts ist ebenfalls ungewiss. In diesem Ergebnis denke ich, dass es nicht viele Punkte in dem Bereich gibt, in denen ich über die Vorhersage relativ unsicher bin ($ x \ le0.1, 0.7 \ le x \ le0.9 $), sondern $ 0.4 \ le x \ le0.6 Im Bereich $ erhöht sich die Standardabweichung nicht, selbst wenn keine Datenpunkte vorhanden sind. Darüber hinaus hängt ein solches Verhalten davon ab, welche Art von Merkmalsvektor verwendet wird. Dieses Mal haben wir die Merkmale von Polynomen verwendet, aber auch andere Merkmale, die Gaußsche Funktionen verwenden, werden häufig verwendet.

Schließlich

Dieses Mal habe ich die Superparameter $ \ alpha und \ beta $ ausgewählt, die zu funktionieren schienen, aber es scheint, dass die Superparameter nur aus den Datenpunkten mithilfe der in Kapitel 3 von PRML als Evidenznäherung bezeichneten Methode geschätzt werden können. .. Ich denke auch daran, das umzusetzen.