Einführung

Ich war besorgt, dass es nicht viele korrekte Implementierungen der Bayes'schen linearen Regression auf der Welt gibt, und es gibt nur wenige Implementierungen, die mehrdimensionale Eingaben unterstützen. Deshalb habe ich sie als benutzerfreundliche Klasse implementiert. Die Beschreibung folgt grundsätzlich PRML.

Klasse implementiert

Das heißt, es ist nicht so groß wie eine Implementierung, es sind ungefähr 50 Zeilen ... Ich werde die ganze Klasse hier schreiben. Die Funktionen sind wie folgt. Es ist sehr einfach.

Geben Sie eine Kombination aus $ \ phi $ und $ t $ an, um die hintere Verteilung zu berechnen
Zufällige Stichprobe von Mittelwert- und Kovarianzparametern für die Bayes'sche lineare Regression
Geben Sie $ \ phi $ an, um den Mittelwert und die Varianz der vorhergesagten Verteilung auszugeben

`Beyes_LR.py`


import numpy as np
from scipy.stats import multivariate_normal

class BeyesLinearRegression:
    def __init__(self, mu, S, beta):
        self.mu = mu
        self.S = S
        self.beta = beta

    def calc_posterior(self, phi, t):
        S_inv = np.linalg.inv(self.S)

        if len(phi.shape) == 1:
            phi = phi.reshape(1, -1)
            t = t.reshape(1, 1)
        self.S = np.linalg.inv(S_inv + self.beta * phi.T @ phi)
        self.mu = self.S @ (S_inv @ self.mu + np.squeeze(self.beta * phi.T @ t))

    def sampling_params(self, n=1, random_state=0):
        np.random.seed(random_state)
        return np.random.multivariate_normal(self.mu, self.S, n)

    def probability(self, x):
        dist = multivariate_normal(mean=self.mu, cov=self.S)
        return dist.logpdf(x)

    def predict(self, phi):
        if len(phi.shape) == 1:
            phi = phi.reshape(1, -1)
        pred = np.array([self.mu.T @ _phi for _phi in phi])
        S_pred = np.array([(1 / self.beta) + _phi.T @ self.S @ _phi for _phi in phi])

        # Above is a simple implementation.
        # This may be better if you want speed.
        # pred = self.mu @ phi.T
        # S_pred = (1 / self.beta) + np.diag(phi @ self.S @ phi.T)
        return pred, S_pred

Der gesamte Code ist auch in Git. (Obwohl es ungefähr 50 Zeilen sind)

GitHub

Über die Bayes'sche lineare Regression

Da im Folgenden die Ableitung der Formel detailliert beschrieben wird, werde ich die Details nicht schreiben.

Berechnung der posterioren Verteilung

Wichtig ist, die Distribution wie folgt zu aktualisieren. Einfach ausgedrückt ist $ \ phi $ die erklärende Variable und $ t $ die Antwort. Die Mittelwert- und Kovarianzmatrizen werden entsprechend aktualisiert.

M_N=S_N(S_0^{-1}m_0+\beta\Phi^Tt)

S_N^{-1}=S_0^{-1}+\beta\Phi^T\Phi

Voraussichtliche Verteilung

Die vorhergesagte Verteilung ist unten gezeigt. Ich werde die Details auch hier weglassen, aber der Punkt ist, dass die Verteilung durch den Mittelwert und die Varianz ausgedrückt wird, wie unten für den neuen Punkt $ x $ gezeigt.

N(m_N^T\phi(x), 1/\beta+\phi(x)^TS_N\phi(x))

Verwenden Sie die implementierte Klasse

Versuchen wir von hier aus die Bayes'sche lineare Regression mit der Klasse, die wir tatsächlich implementiert haben. Die Klasse, die ich erstellt habe, muss $ \ phi $ als Feature von $ x $ erhalten. In einer gängigen Implementierung ist der Generierungsteil von $ \ phi $ (z. B. ein Polynom) ebenfalls in der Klasse enthalten, und es ist schwierig zu sagen, ob eine Bayes'sche lineare Regression durchgeführt oder Features mithilfe von Polynomen entworfen werden. Aber hier ist es getrennt.

Wenn Sie also die Originaldaten $ x $, $ y $ eingeben und dann die Funktion zum Erstellen von $ \ phi $ implementieren, ist dies grundsätzlich in Ordnung.

Regression mit Sinuswellendaten

Versuchen wir es zuerst mit Toy-Daten.

Datengenerierung und Design von Funktionsfunktionen

Die Eingangsdaten sind die Beobachtungsdaten, indem der wahren Verteilung der Sinuswelle Rauschen hinzugefügt wird. Darüber hinaus ist die Merkmalsfunktion als zusammengesetzte Welle von Dreiecksfunktionen mit mehreren Frequenzen ausgelegt. Die Methode x_to_phi vektorisiert 10 Wellen und _phi repräsentiert eine zusammengesetzte Welle. Die Amplitude ist der Parameter, der durch die Bayes'sche lineare Regression erhalten wird. Das Bild unten ist mathematisch. (Wenn Sie darüber nachdenken, ist der erste Punkt Null ... ich brauche ihn nicht ...)

y=w_1sin(0)+w_2sin(2\pi x)+w_3sin(2\times2\pi x)+\cdots+w_9sin(9\times2\pi x)

def x_to_phi(x):
    if len(x.shape) == 1:
        x = x.reshape(-1, 1)
    return np.concatenate([np.sin(2 * np.pi * x * m) for m in range(0, 10)], axis=1)


def _phi(x, params):
    return np.array([p * np.sin(2 * np.pi * x * i) for i, p in enumerate(params)]).sum(axis=0)

Klicken Sie hier für den eigentlichen Datengenerierungsteil.

x = np.arange(0, 1, 0.01)
phi = x_to_phi(x)

e = np.random.randn(len(x)) / 10
y_true = np.sin(2 * np.pi * x)
y = y_true + e

Wenn nur ein Punkt beobachtet wird

Betrachten Sie zunächst den Fall, in dem nur ein Punkt der 50. Daten beobachtet wird.

train_idx = 50
x_train = x[train_idx]
phi_train = phi[train_idx]
y_train = y[train_idx]
plt.scatter(x_train, y_train, c='crimson', marker='o', label='observation')
plt.plot(x, y_true, label='true')

Berechnen wir sofort die posteriore Verteilung für diese Daten. Wenn Sie nur lernen möchten, ist es eine Zeile wie diese:

#Anfangswert der Bayes'schen linearen Regression
input_dim = phi.shape[1]
mu = np.zeros(input_dim)
S = np.identity(input_dim)
sigma = 0.1
beta = 1.0 / (sigma ** 2)

#Modelldefinition und Training
beyes_linear_model = BeyesLR.BeyesLinearRegression(mu, S, beta)
beyes_linear_model.calc_posterior(phi_train, y_train)

Einige zufällig abgetastete Wellenformen aus der posterioren Verteilung nach dem Training werden mit grün gepunkteten Linien angezeigt, und die vorhergesagte Verteilung wird hellblau angezeigt. Die meisten von ihnen sind blau, weil ich nur einen Punkt gelernt habe.

Bei der Beobachtung von 50 Punkten

Machen wir genau das Gleiche aus 50 Beobachtungsdaten. Der einzige Unterschied im Code besteht darin, zufällig 50 train_idx auszuwählen. Die vorhergesagte Verteilung der Ergebnisse ist in der folgenden Abbildung dargestellt.

Zurück im Werbedatensatz

Als nächstes kommt ein echtes Problem, das auch eine mehrdimensionale lineare Regression löst. Wenn Sie Features in mehreren Dimensionen extrahieren, werden die Dimensionen zu stark vergrößert und sind schwer zu verstehen. Daher ist $ \ phi $ eine lineare Funktion.

Eingabedaten

Es ist ein Werbedatensatz, der der ISLR vertraut ist. Dieses Mal werden wir TV- und Radio-Werbekosten als Input und Verkäufe als Antwort verwenden.

Dieses Mal ist $ \ phi $ linear, fügen Sie also einfach den Abschnittsbegriff hinzu. Die Formel für die Regression lautet wie folgt. $Sales = w_0+w_1TV+w_2Radio$

def x_to_phi(x, typ='linear', degree=3):
    if len(x.shape) == 1:
        x = x.reshape(-1, 1)
    return np.concatenate([np.ones(x.shape[0]).reshape(-1, 1), x], axis=1)


df = pd.read_csv(ADVERTISING_DATASET)
x = df[['TV', 'Radio']].values
y = df['Sales'].values

phi = x_to_phi(x)
x_train, x_test, phi_train, phi_test, y_train, y_test = \
    train_test_split(x, phi, y, train_size=0.05, random_state=0)

Sie müssen lediglich wie im vorherigen Beispiel lernen.

input_dim = phi.shape[1]
mu = np.zeros(input_dim)
S = np.identity(input_dim)
sigma = 10
beta = 1.0 / (sigma ** 2)

beyes_linear_model = BeyesLR.BeyesLinearRegression(mu, S, beta)
beyes_linear_model.calc_posterior(phi_train, y_train)

Im Code wird train_size auf 0,05 gesetzt, aber die Regressionsebene, die beim Ändern gezeichnet wird, lautet wie folgt. Es wird durch Bayes'sche lineare Regression gelernt, und 5 Ebenen werden durch Zufallsstichprobe extrahiert und gezeichnet. Zufällig, wenn die Anzahl der Lernmuster gering ist Eine Ebene wird gezeichnet, konvergiert jedoch mit zunehmender Anzahl von Datenpunkten.

abschließend

Zum Schluss noch eine kleine Förderung der Bayes'schen linearen Regression. Obwohl es noch einige Teile gibt, die nicht vollständig verstanden werden, wird das Design der Merkmalsextraktionsfunktion $ \ phi $ der Bayes'schen linearen Regression wichtig. Ich erkenne, dass die Gaußsche Prozessregression die verteilte, gemeinsam verteilte Matrix als Planungsmatrix unter Verwendung von Kernelfunktionen behandelt, ohne dies explizit aufzuschreiben. Die Erfahrung zeigt jedoch, dass die lineare Bayes'sche Regression unter dem Gesichtspunkt der Erklärung ausreichend ist, wenn Vorkenntnisse vorliegen, die gute Merkmale extrahieren können. Aufgrund der Berechnungsmethode der Bayes'schen linearen Regression kann das sequentielle Lernen so wie es ist durchgeführt werden. Sie müssen lediglich die berechnete posteriore Verteilung als vorherige Verteilung lernen. Es ist nicht erforderlich, die Gramm-Matrix wie die Gaußsche Prozessregression neu zu berechnen. Es kann eine Online-Version geben, aber ...

Übrigens, haben Sie ein komfortables Bayes'sches lineares Rücklaufleben!

Python-Implementierung der Bayes'schen linearen Regressionsklasse