PRML Kapitel 4 Bayesianische logistische Regression Python-Implementierung

Bei der Kurvenanpassung behandeln wir es wie Bayes und berechnen die posteriore Vorhersageverteilung, aber es besteht der Eindruck, dass so etwas bei der Klassifizierung nicht so häufig gemacht wird. Deshalb werden wir diesmal die logistische Regression, die häufig bei der Klassifizierung verwendet wird, eher Bayes'sch machen. Ich möchte den Code implementieren, um die posteriore Vorhersageverteilung zu berechnen.

Bayes Logistic Return

Wie ich oben geschrieben habe, habe ich nicht viel oder keinen Code gesehen, der die vorhergesagte Verteilung in der Klassifizierung berechnet (vielleicht nur ein paar Codes, die ich gesehen habe). Um die vorhergesagte Verteilung auf Bayes'sche Weise zu berechnen, muss die posteriore Verteilung der Gewichtsparameter verwendet werden, um die Parameter zu integrieren und zu eliminieren. Da die logistische Regression jedoch die logistische Sigmoidfunktion verwendet, ist es nicht möglich, die Parameter analytisch zu integrieren. Daher wird die vorhergesagte Verteilung unter Verwendung der Laplace-Näherung ungefähr erhalten.

Logistische Rückgabe

Betrachten wir zunächst die logistische Regression, die zwei Klassen klassifiziert. Sei $ \ phi $ der Merkmalsvektor eines Punktes $ x $ und sei dieser Punkt zur Klasse $ C_1 $

p(C_1|\phi) = y(\phi) = \sigma({\bf w}^\intercal\phi)

Es wird ausgedrückt als. Jedoch,\sigma(\cdot)Ist eine logistische Sigmoidfunktion\sigma(x)={1\over1+\exp(-x)}Und. Da wir hier eine Zwei-Klassen-Klassifizierung in Betracht ziehen, KlassenC_1Wenn Sie die Wahrscheinlichkeit kennen, zur anderen Klasse zu gehörenC_2Die Wahrscheinlichkeit der Zugehörigkeit zup(C_2|\phi(x)) = 1 - p(C_1|\phi(x))Gegeben in. Daher danach die KlasseC_1Es wird nur die Wahrscheinlichkeit einer Zugehörigkeit berücksichtigt.

Trainingsdatensatz $ \ {\ phi_n, t_n \} _ {n = 1} ^ N $ als $ \ phi_n = \ phi (x_n) $, $ t_n \ in \ {0,1 \} $ Vorausgesetzt, die Wahrscheinlichkeitsfunktion basiert auf der Bernoulli-Verteilung

\begin{align}
p({\bf t}|{\bf\Phi},{\bf w}) &= \prod_{n=1}^N{\rm Bern}(t_n|y_n)\\
&= \prod_{n=1}^N y_n^{t_n}(1 - y_n)^{1 - t_n}
\end{align}

Wird sein. $ {\ Bf \ Phi} $ ist jedoch eine Planungsmatrix, in der die horizontalen Vektoren $ \ phi ^ {\ rm T} $ vertikal angeordnet sind. Wie immer erfolgt die Definition der Kostenfunktion als negative Log-Likelihood-Funktion in Form einer Kreuzentropie.

E({\bf w}) = -\ln p({\bf t}|{\bf\Phi},{\bf w}) = -\sum_{n=1}^N\{t_n\ln y_n + (1 - t_n)\ln(1 - y_n)\}

Die Ergebnisse der Differenzierung erster und zweiter Ordnung dieser Kostenfunktion für den Parameter $ {\ bf w} $ sind

\begin{align}
\nabla E({\bf w}) &= \sum_{n=1}^N (y_n - t_n) {\bf\phi}_n = {\bf\Phi}^\intercal({\bf y} - {\bf t})\\
{\bf H} = \nabla\nabla E({\bf w}) &= \sum_{n=1}^N y_n(1 - y_n)\phi_n\phi_n^\intercal = {\bf\Phi}^\intercal{\bf R}{\bf\Phi}
\end{align}

Die Matrix $ {\ bf R} $ ist jedoch eine diagonale Matrix, deren Elemente $ R_ {nn} = y_n (1-y_n) $ sind. Mit diesen Ergebnissen

{\bf w}^{new} = {\bf w}^{old} - {\bf H}^{-1}\nabla E({\bf w})

As wird die Aktualisierung wiederholt, bis der Parameter $ {\ bf w} $ konvergiert.

Laplace-Näherung

Vorherige Verteilung für Parameter $ {\ bf w} $ $ p ({\ bf w}) = \ mathcal {N} ({\ bf w} | {\ bf 0}, \ alpha ^ {-1} {\ bf I. }) Führen Sie $ ein und verwenden Sie die Laplace-Näherung, um eine ungefähre hintere Verteilung in Form einer Gaußschen Funktion zu finden. Aus dem Satz von Bayes

p({\bf w}|{\bf t},{\bf\Phi}) \propto p({\bf t}|{\bf\Phi},{\bf w})p({\bf w})

Also, wenn Sie diesen Logarithmus nehmen

\ln p({\bf w}|{\bf t},{\bf\Phi}) = -{\alpha\over2}{\bf w}^\intercal{\bf w} + \sum_{n=1}^N\{t_n\ln y_n + (1 - t_n)\ln(1 - y_n)\}

Wird sein. Um die Gaußsche Näherung daraus zu berechnen, maximieren Sie zuerst die obige Gleichung für $ {\ bf w} $, dann den Parameter $ {\ bf w} _ {MAP} $ und das Differential zweiter Ordnung $ S_N ^ {-1 } Mit $ die ungefähre hintere Verteilung

q({\bf w}) = \mathcal{N}({\bf w}|{\bf w}_{MAP}, S_N)

Ist erforderlich als.

Ex-post-Vorhersageverteilung

Die posteriore Vorhersageverteilung für den neuen Punkt $ \ phi $ kann durch Marginalisieren der posterioren Verteilung der Parameter erhalten werden.

\begin{align}
p(C_1|\phi,{\bf t},{\bf\Phi}) &= \int p(C_1|\phi,{\bf w})p({\bf w}|{\bf t},{\bf\Phi}){\rm d}{\bf w}\\
&\approx \int \sigma({\bf w}^\intercal\phi)\mathcal{N}({\bf w}|{\bf w}_{MAP},S_N){\rm d}{\bf w}\\
&\approx \sigma\left({\mu\over\sqrt{1 + \pi\sigma^2/8}}\right)
\end{align}

Jedoch,

\begin{align}
\mu &= {\bf w}_{MAP}^\intercal\phi\\
\sigma^2 &= \phi^\intercal S_N\phi
\end{align}

Und. Es werden zwei Näherungen verwendet: Laplace-Näherung der posterioren Verteilung des Parameters $ {\ bf w} $ und Probit-Funktionsnäherung der logistischen Sigmoidfunktion. Siehe PRML für eine detaillierte Formeltransformation.

Implementierung

Paket

Zusätzlich zu der üblichen Matplotlib und Numpy wird eine Standard-Python-Bibliothek namens itertools verwendet.

import itertools
import matplotlib.pyplot as plt
import numpy as np

Entwurfsmatrix: Polygon-Features

Konvertieren eines zweidimensionalen Punktes $ (x_1, x_2) $ in einen quadratischen polymorphen Merkmalsvektor

\phi(x_1,x_2) =
\begin{bmatrix}
1\\
x_1\\
x_2\\
x_1^2\\
x_1x_2\\
x_2^2
\end{bmatrix}

Es wird sein. Unten finden Sie eine Klasse, die mehrere zweidimensionale Punkte in eine Entwurfsmatrix konvertiert. Zum Beispiel die Planungsmatrix der quadratischen polymorphen Merkmale der zweidimensionalen Punkte $ (2,5) $ und $ (-3,1) $

{\bf\Phi}=
\begin{bmatrix}
1 & 2 & 5 & 2^2 & 2 \times 5 & 5^2\\
1 & 3 & -1 & 3^2 & 3 \times (-1) & (-1)^2
\end{bmatrix}

Es sieht aus wie.

class PolynomialFeatures(object):

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

    def transform(self, x):
        x_t = x.transpose()
        features = [np.ones(len(x))]
        for degree in xrange(1, self.degree + 1):
            for items in itertools.combinations_with_replacement(x_t, degree):
                features.append(reduce(lambda x, y: x * y, items))
        return np.array(features).transpose()

Logistische Rückgabe

Es ist nicht notwendig, die Klassen, die eine logistische Regression durchführen, und eine Bayes'sche logistische Regression zu trennen, aber dieses Mal haben wir sie getrennt, um den Unterschied zwischen den beiden zu machen.

class LogisticRegression(object):

    def __init__(self, iter_max, alpha=0):
        self.iter_max = iter_max
        self.alpha = alpha

    def _sigmoid(self, a):
        return np.divide(1, 1 + np.exp(-a))

    def fit(self, X, t):
        self.w = np.zeros(np.size(X, 1))
        for i in xrange(self.iter_max):
            w = np.copy(self.w)
            y = self.predict_proba(X)
            grad = X.T.dot(y - t) + self.alpha * w
            hessian = (X.T.dot(np.diag(y * (1 - y))).dot(X)
                       + self.alpha * np.identity(len(w)))
            try:
                self.w -= np.linalg.solve(hessian, grad)
            except np.linalg.LinAlgError:
                break
            if np.allclose(w, self.w):
                break
            if i == self.iter_max - 1:
                print "weight parameter w may not have converged"

    def predict(self, X):
        return (self.predict_proba(X) > 0.5).astype(np.int)

    def predict_proba(self, X):
        return self._sigmoid(X.dot(self.w))

Bayes Logistic Return

Bei der Schätzung der maximalen posterioren Wahrscheinlichkeit eines Parameters wird nicht nur der häufigste Wert, sondern auch eine ungefähre Varianz erhalten. Die obige logistische Regression suchte auch nach der Hessen-Matrix (der Präzisionsmatrix der ungefähren Parameter), aber ich habe diese Matrix in keiner Variablen gespeichert, da sie diese Matrix nur beim Schätzen verwendet. Die Bayes'sche logistische Regression verwendet jedoch auch die Varianz der Parameter, um die vorhergesagte Verteilung zu berechnen und zu speichern. Dann berechnet die Predict_Dist-Methode die vorhergesagte Verteilung unter Verwendung der Verteilungsmatrix dieses Parameters.

class BayesianLogisticRegression(LogisticRegression):

    def __init__(self, iter_max, alpha=0.1):
        super(BayesianLogisticRegression, self).__init__(iter_max, alpha)

    def fit(self, X, t):
        super(BayesianLogisticRegression, self).fit(X, t)
        y = self.predict_proba(X)
        hessian = (X.T.dot(np.diag(y * (1 - y))).dot(X)
                   + self.alpha * np.identity(len(self.w)))
        self.w_var = np.linalg.inv(hessian)

    def predict_dist(self, X):
        mu_a = X.dot(self.w)
        var_a = np.sum(X.dot(self.w_var) * X, axis=1)
        return self._sigmoid(mu_a / np.sqrt(1 + np.pi * var_a / 8))

Ganzer Code

import itertools
import matplotlib.pyplot as plt
import numpy as np


class PolynomialFeatures(object):

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

    def transform(self, x):
        x_t = x.transpose()
        features = [np.ones(len(x))]
        for degree in xrange(1, self.degree + 1):
            for items in itertools.combinations_with_replacement(x_t, degree):
                features.append(reduce(lambda x, y: x * y, items))
        return np.array(features).transpose()


class LogisticRegression(object):

    def __init__(self, iter_max, alpha=0):
        self.iter_max = iter_max
        self.alpha = alpha

    def _sigmoid(self, a):
        return np.divide(1, 1 + np.exp(-a))

    def fit(self, X, t):
        self.w = np.zeros(np.size(X, 1))
        for i in xrange(self.iter_max):
            w = np.copy(self.w)
            y = self.predict_proba(X)
            grad = X.T.dot(y - t) + self.alpha * w
            hessian = (X.T.dot(np.diag(y * (1 - y))).dot(X)
                       + self.alpha * np.identity(len(w)))
            try:
                self.w -= np.linalg.inv(hessian).dot(grad)
            except np.linalg.LinAlgError:
                break
            if np.allclose(w, self.w):
                break
            if i == self.iter_max - 1:
                print "weight parameter w may not have converged"

    def predict(self, X):
        return (self.predict_proba(X) > 0.5).astype(np.int)

    def predict_proba(self, X):
        return self._sigmoid(X.dot(self.w))


class BayesianLogisticRegression(LogisticRegression):

    def __init__(self, iter_max, alpha=0.1):
        super(BayesianLogisticRegression, self).__init__(iter_max, alpha)

    def fit(self, X, t):
        super(BayesianLogisticRegression, self).fit(X, t)
        y = self.predict_proba(X)
        hessian = (X.T.dot(np.diag(y * (1 - y))).dot(X)
                   + self.alpha * np.identity(len(self.w)))
        self.w_var = np.linalg.inv(hessian)

    def predict_dist(self, X):
        mu_a = X.dot(self.w)
        var_a = np.sum(X.dot(self.w_var) * X, axis=1)
        return self._sigmoid(mu_a / np.sqrt(1 + np.pi * var_a / 8))


def create_data_set():
    x = np.random.normal(size=50).reshape(-1, 2)
    y = np.random.normal(size=50).reshape(-1, 2)
    y += np.array([2., 2.])
    return (np.concatenate([x, y]), np.concatenate([np.zeros(25), np.ones(25)]))


def main():
    x, labels = create_data_set()
    colors = ['blue', 'red']
    plt.scatter(x[:, 0], x[:, 1], c=[colors[int(label)] for label in labels])

    features = PolynomialFeatures(degree=3)

    classifier = BayesianLogisticRegression(iter_max=100, alpha=0.1)
    classifier.fit(features.transform(x), labels)

    X_test, Y_test = np.meshgrid(np.linspace(-2, 4, 100), np.linspace(-2, 4, 100))
    x_test = np.array([X_test, Y_test]).transpose(1, 2, 0).reshape(-1, 2)
    probs = classifier.predict_proba(features.transform(x_test))
    Probs = probs.reshape(100, 100)
    dists = classifier.predict_dist(features.transform(x_test))
    Dists = dists.reshape(100, 100)
    levels = np.linspace(0, 1, 5)
    cp = plt.contour(X_test, Y_test, Probs, levels, colors='k', linestyles='dashed')
    plt.clabel(cp, inline=True, fontsize=10)
    plt.contourf(X_test, Y_test, Dists, levels, alpha=0.5)
    plt.colorbar()
    plt.xlim(-2, 4)
    plt.ylim(-2, 4)
    plt.show()


if __name__ == '__main__':
    main()

Ergebnis

• Lerndaten: Blaue (Klasse 0), rote (Klasse 1) Punkte
• ** Schwarze gepunktete Linie: Logistische Regressionsausgabe **, eine Linie mit einer Wahrscheinlichkeit der Zugehörigkeit zur Klasse 1 von oben 0,75,0,5,0,25, obwohl sie schwer zu erkennen ist
• ** Farbcodierung: Bayesianische logistische Regressionsausgabe **, jeweils mit einer Wahrscheinlichkeit der Zugehörigkeit zur Klasse 1
• Rot: 0,75 ~
• Gelb: 0,5 bis 0,75

--Blau: ~ 0,25 Betrachtet man dies, so ist die Linie mit einer Wahrscheinlichkeit von 0,5 (die gepunktete Linie in der Mitte und die Grenze zwischen hellblau und gelb) in beiden Fällen üblich, während die Linien von 0,25 und 0,75 wie Bayes behandelt werden. Der eine ist weiter von der 0,5-Linie entfernt und sagt uns, dass unklar ist, zu welcher Klasse sie gehören. Wenn oben links in der Abbildung keine Datenpunkte vorhanden sind, wird die vorhergesagte Verteilung mehrdeutiger als in anderen Regionen. Andererseits ist die Breite von 0,25 bis 0,75 schmal, da es viele Datenpunkte um die Mitte gibt. Nach der Berechnung der Wahrscheinlichkeit der Zugehörigkeit zu Klasse 1 nach einer der beiden Methoden wird sie nach einem bestimmten Kriterium einer der beiden Klassen zugewiesen. Wenn das Kriterium jedoch eine Minimierung der Fehlklassifizierung ist, ändert sich nicht, welches Kriterium verwendet wird. Weil die Linie mit einer Wahrscheinlichkeit von 0,5 üblich ist. Wenn Sie jedoch kompliziertere Kriterien verwenden möchten, können Sie diese wie Bayesian behandeln. Schließlich implementierte ich diesmal eine Methode zur Durchführung der Bayes'schen logistischen Regression unter Verwendung der Laplace-Näherung. Der Vorteil einer Bayes'schen Behandlung besteht darin, dass sie sich in Bereichen, in denen es viele Datenpunkte gibt und in denen es keine gibt, anders verhält. Zusätzlich zur Verwendung der Laplace-Näherung scheint es eine Methode zu geben, die vorhergesagte Verteilung unter Verwendung der Variante Bayes abzuleiten. Kapitel 10 von PRML verwendet auch die Variante Bayes, um den Superparameter $ \ alpha $ zu schätzen. Ich denke, ich werde das auch umsetzen.