PRML Chapitre 1 Implémentation de Python pour l'ajustement de courbe bayésienne

Pour ma propre étude, j'ai choisi certaines des méthodes introduites dans PRML et j'ai décidé de les implémenter en Python. Fondamentalement, nous coderons avec la règle selon laquelle seul numpy peut être utilisé autre que la bibliothèque Python standard (dont scipy peut être utilisé). Cependant, matplotlib etc. sont des exceptions car ils n'ont rien à voir avec la partie algorithme.

Ajustement de courbe de Bayes

L'ajustement de courbe de Bayes est une méthode d'ajustement d'une courbe en la traitant complètement comme Bayes. Alors que la distribution prédite par l'estimation la plus probable a la même variance en tous points, l'ajustement de la courbe bayésienne utilise la distribution a posteriori des paramètres pour obtenir la distribution prédite de manière bayésienne (en utilisant le théorème d'addition et de multiplication de la probabilité). Vous pouvez calculer la variance pour chaque point avec. Il calculera si vous êtes relativement confiant sur chaque point.

la mise en oeuvre

Comme indiqué au tout début, l'implémentation de l'algorithme utilise numpy.

import matplotlib.pyplot as plt
import numpy as np

Calcul de la matrice de conception

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

Estimation et prédiction de l'ajustement de la courbe bayésienne

D'abord, estimez la distribution a posteriori du paramètre $ {\ bf w} . $ p({\bf w}|{\bf x}, {\bf t}, \alpha, \beta) \propto p({\bf t}|{\bf\Phi}{\bf w}, \beta)p({\bf w}|\alpha) $$ Ici, $ {\ bf x} et {\ bf t} $ sont des ensembles d'entraînement. De plus, $ p $ sur le côté droit suppose une fonction gaussienne, et les paramètres de précision sont respectivement $ \ alpha $ et $ \ beta . Ensuite, la distribution de probabilité sur le côté gauche devient également une distribution gaussienne. $ p({\bf w}|{\bf x}, {\bf t}, \alpha, \beta) = \mathcal{N}({\bf m}, {\bf S}) $ Cependant $ {\bf m} = \beta{\bf S}{\bf\Phi}^{\rm T}{\bf t},~{\bf S}^{-1} = \alpha{\bf I} + \beta{\bf\Phi}^{\rm T}{\bf\Phi} $$ En utilisant cette distribution postérieure, la distribution de probabilité de $ t $ pour le nouveau point $ 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} $ C'est aussi une distribution gaussienne lorsqu'elle est calculée. $ 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

Code entier

Le code ci-dessous est ce que j'ai implémenté, y compris ce qui précède. Le flux général est

1. Créez des données d'entraînement (x, t = create_toy_data (func, ...))
2. Calculez la distribution postérieure du paramètre de poids $ {\ bf w} $ (regression.fit (X, t))
3. Calculez la distribution prédite en utilisant la distribution postérieure (y, y_std = regression.predict (X_test))
4. Illustrer les résultats Il est devenu.

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

résultat

Il a été prédit en effectuant un ajustement de courbe bayésienne à partir des points de données (points bleus) générés en ajoutant du bruit à la fonction $ {\ rm sin} (2 \ pi x) $. La valeur moyenne de chaque $ x $ de la distribution prévue est indiquée en rouge et l'écart type est indiqué en rose. La largeur de la bande rose correspond à l'incertitude de la prédiction. Si vous calculez la distribution de prédiction par l'estimation la plus probable, la largeur de la bande rose sera constante quel que soit $ x , et la prédiction de tout point sera également incertaine. Dans ce résultat, je pense qu'il n'y a pas beaucoup de points ( x \ le0.1, 0.7 \ le x \ le0.9 $) dans la zone où je suis relativement incertain de la prédiction, mais $ 0.4 \ le x \ le0.6 Dans la zone $, l'écart type n'augmente pas même s'il n'y a pas de points de données. De plus, ce comportement dépend du type de vecteur de caractéristiques utilisé. Cette fois, nous avons utilisé les caractéristiques des polynômes, mais d'autres caractéristiques utilisant des fonctions gaussiennes sont également souvent utilisées.

finalement

Cette fois, j'ai choisi les super paramètres $ \ alpha et \ beta $ qui semblaient fonctionner, mais il semble que les super paramètres ne puissent être estimés qu'à partir des points de données en utilisant la méthode appelée approximation des preuves au chapitre 3 de PRML. .. Je pense mettre cela en œuvre aussi.