Logistische Regression

Im ersten Abschnitt basteln wir eine S-förmige Kurve aus dem $\tanh$, haben aber keine Interpretation dafür. Danach bauen wir die S-Kurve von der Wahrscheinlichkeitsrechnung her auf, und dann fällt die Interpretation von selbst ab.

Lineare Regression zur Klassifikation¶

Wenn wir an Regression denken, dann oft an eine Regressionsgerade in einem $x$-$y$-Plot, bei dem auf der $x$-Achse die unabhängige Variable (der Regressor) und auf der $y$-Achse die abhängige Variable (der Regressand) aufgetragen ist. Bei der Klassifikation ist der Regressand die Klasse, $y \in \{0, 1\}$.

Schauen wir uns das Schritt für Schritt an, für einen eindimensionalen Regressor $x \in \mathbb{R}$ mit zwei Klassen, deren $x$-Werte sich überlappen:

1. Im eindimensionalen Plot liegen alle Punkte der Klasse 0 auf der Höhe $y = 0$ und alle Punkte der Klasse 1 auf $y = 1$.

2. Die lineare Regression fittet eine Gerade $\hat{y} = w x + b$ durch diese Punktwolke. Die Gerade nimmt Werte in ganz $\mathbb{R}$ an, auch weit außerhalb von $[0, 1]$.

3. Harte Klassifikation: Wir entscheiden für Klasse 1, wenn $\hat{y} > \tfrac12$. Die Entscheidungsgrenze ist der Punkt $x^*$ mit $w x^* + b = \tfrac12$; die vorhergesagte Klasse als Funktion von $x$ ist eine Sprungfunktion.

4. Für einen zweidimensionalen Regressor $x \in \mathbb{R}^2$ wird aus der Geraden eine Regressionshyperebene über der $(x_1, x_2)$-Ebene, und die Entscheidungsgrenze $\{ w \cdot x + b = \tfrac12 \}$ ist eine Gerade in der $(x_1, x_2)$-Ebene – genau wie beim Perzeptron.

S-Förmigkeit herstellen¶

Die Sprungfunktion $[\hat{y} > \tfrac12]$ ist unbefriedigend: Sie ist nicht differenzierbar, und sie verrät nicht, wie sicher die Entscheidung ist. Anstelle von $[\hat{y} > \tfrac12]$ könnte man $[\tanh(\hat{y}) > 0]$ verwenden; $\tanh$ liegt zwischen -1 und +1. Um zwischen 0 und 1 zu landen, kann man z.B.

nutzen. Das liefert eine S-förmige Kurve, die glatt von 0 nach 1 übergeht. Aber: Die Interpretation einer Zahl zwischen 0 und 1 als Wahrscheinlichkeit ist damit noch unklar. Genau das klären die nächsten Abschnitte dieses Kapitels, über Odds und Logit.

Wahrscheinlichkeiten und Chancen (Odds)¶

Chancen für ein Bernoulli-Experiment¶

Angenommen, wir haben ein Bernoulli-Experiment (Münzwurf) mit Parameter $p \in (0, 1)$ – zur Vereinfachung mit $p \neq 0$ und $p \neq 1$. Die Wahrscheinlichkeit (probability) für Erfolg ist $P(\text{Erfolg}) = p$. Die Chancen (odds) für Erfolg sind $p : (1-p)$.

Beispiel: Einen Standardwürfel (D6) rollen und gucken, ob es eine 6 ist. Dann ist $P(D = 6) = \tfrac16 = p$ die Wahrscheinlichkeit, und $1 : 5 = \tfrac15$ sind die Chancen für eine 6; $5 : 1 = 5$ sind die Chancen für eine Nicht-6. Als Formel:

Allgemein rechnet man von Wahrscheinlichkeiten zu Odds und zurück über

Die Skala für Odds¶

Wahrscheinlichkeiten liegen immer in $[0, 1]$. Wie sieht die Skala für Odds aus?

• Wenn $p = \tfrac12$, so ist $1 - p = \tfrac12$, also $\operatorname{odds}(p) = 1$.

• Wenn $0 < p < \tfrac12$, so ist $1 > 1 - p > \tfrac12$, also $0 < \operatorname{odds}(p) = \frac{p}{1-p} < \frac{p}{1/2} = 2p < 1$. Damit ist $\operatorname{odds}(p) \in (0, 1)$.

• Wenn $1 > p > \tfrac12$, so ist $0 < 1 - p < \tfrac12$, also $1 < 2p < \operatorname{odds}(p) < \frac{1}{1-p} < \infty$. Damit ist $\operatorname{odds}(p) \in (1, \infty)$.

Die Fälle $p < \tfrac12$ und $p > \tfrac12$ sind symmetrisch, d.h. Erfolg und Misserfolg sind äquivalent – aber $\operatorname{odds}(p) \in (0, \infty)$ ist nicht symmetrisch: $(0, 1)$ ist beschränkt, $(1, \infty)$ ist unbeschränkt. Die Odds für eine 6 sind $\tfrac15$, die Odds für eine Nicht-6 sind 5.

Beachte: Logarithmieren rettet die Symmetrie. Es ist

und damit $\ln(1 : 5) = -\ln(5 : 1)$. Die Werte liegen in $(-\infty, +\infty)$.

Logistic Unit¶

Schon mal irgendwo gesehen? Im Data-Science-Skript haben wir den Informationsgehalt (die Überraschung) eines Ereignisses mit Wahrscheinlichkeit $p$ als $-\log p$ definiert und die Entropie als erwartete Überraschung (Informationsgehalt, Entropie). Der Logarithmus steht dort aus einem Grund, der uns noch oft nützen wird: Für unabhängige Ereignisse multiplizieren sich die Wahrscheinlichkeiten, und $\log$ macht daraus eine Summe. Zwei Überraschungen addieren sich. Ein fairer Münzwurf ist 1 bit Überraschung, eine 6 beim Würfeln $\log_2 6 \approx 2{,}58$ bit, beides zusammen $3{,}58$ bit.

Für das Bernoulli-Experiment ist die Entropie

und ihre Ableitung nach $p$ ist

(im DS-Skript: Die binäre Entropiefunktion, dort mit $\log_2$ und kleinem $h$). Interpretation: Bei Erhöhen von $p$ ist die Abnahme an Informationsgewinn durch $\operatorname{logit}(p)$ gegeben. Entropie misst man in bit ($\log_2$) oder nat ($\ln$), die Ableitung in logit. Zum Beispiel ist $H(0{,}5) = 1 \text{ bit} = 0{,}69 \text{ nat}$.

Logistische Sigmoidfunktion¶

Allgemeine S-Förmigkeit¶

Sigma $\approx$ S, Sigmoid $\approx$ S-förmig.

$M$ ist der Maximalwert (die obere Asymptote), $k$ die Steilheit, $\mu$ die Position des Wendepunkts. Bei $\sigma(w \cdot x + b)$ übernimmt $\|w\|$ die Rolle von $k$, siehe Kapitel „Neuron und Perzeptron“: Die Länge von $w$ bestimmt, wie steil die Wahrscheinlichkeit mit dem Abstand zur Entscheidungsgrenze wächst. Das Standardsigmoid steigt monoton von 0 (für $x \to -\infty$) nach 1 (für $x \to +\infty$) und hat bei $x = 0$ den Wert $\sigma(0) = \tfrac12$.

Symmetrie¶

Es gilt

also $\sigma(-x) = 1 - \sigma(x)$: Der Graph ist punktsymmetrisch um $(0, \tfrac12)$.

Andere Sigmoidfunktionen: Tangens hyperbolicus¶

Die Parameter $M$ und $k$ sind Skalierungsfaktoren. Möchte man die Skala des Standardsigmoids von $[0, 1]$ auf z.B. $[-1, 1]$ bringen, so nimmt man $2\sigma(x) - 1$, und mit $k = 2$:

Der Tangens hyperbolicus ist also nur ein umskaliertes Standardsigmoid, $\tanh(x) = 2\sigma(2x) - 1$.

Ableitung¶

Wir benötigen später die Ableitung der Sigmoidfunktion. Eventuell bekannt ist $\tanh'(x) = 1 - \tanh^2(x)$. Ähnlich gilt

Der Beweis mit der Quotientenregel ist eine lohnende Übung. Bemerkenswert daran: Die Ableitung lässt sich aus dem Funktionswert berechnen, ohne $e^{-x}$ noch einmal auszuwerten – das ist beim Training praktisch.

Sigmoid und Logit¶

Wir haben zwei Funktionen:

Also lassen sich beide Verkettungen bilden, $\sigma \circ \operatorname{logit} \colon (0, 1) \to (0, 1)$ und $\operatorname{logit} \circ \sigma \colon \mathbb{R} \to \mathbb{R}$. Was haben die Funktionen miteinander zu tun?

Für $p \in (0, 1)$ ist

Für $x \in \mathbb{R}$ ist – erst $\sigma(x)$ mit $e^x$ erweitern, dann $\ln$ anwenden und $\ln(e^x) = x$ nutzen –

Wir haben gesehen: $\sigma$ und $\operatorname{logit}$ sind Inverse zueinander,

Darum nennt man $\sigma$ als Inverse von $\operatorname{logit}$ auch expit.

Wie interpretieren wir das?¶

Gegeben einen Wert $x \in \mathbb{R}$ gibt $\sigma(x)$ einen Wert $p := \sigma(x) \in (0, 1)$, sodass $x = \operatorname{logit}(p)$. Jeder Wert $x$ lässt sich also als Logit auffassen – als logarithmierte Chance –, und $\sigma(x)$ gibt die Wahrscheinlichkeit des Bernoulli-Experiments an, zu der $x$ der Logit-Wert ist. Das ist die Interpretation, die uns im ersten Abschnitt gefehlt hat: Eine Zahl zwischen 0 und 1, die aus $\sigma$ kommt, ist eine Wahrscheinlichkeit, wenn die Eingabe ein Logit war.

Logistische Regression¶

Alle Zutaten sind bereit, jetzt wird gekocht. Die lineare Regression spuckt $\hat{y}$-Werte aus, die in ganz $\mathbb{R}$ liegen. Wenn $y = 0$ und $y = 1$ zu den zwei Klassen korrespondieren, so fassen wir $\hat{y} = w \cdot x + b$ als Logit auf, und die Werte $\sigma(\hat{y})$ geben die Wahrscheinlichkeit an, in Klasse 1 zu liegen. Wenn $\sigma(\hat{y}) > \tfrac12$, so entscheiden wir für Klasse 1 – und weil $\sigma$ monoton mit $\sigma(0) = \tfrac12$ ist, ist das dasselbe wie $\hat{y} > 0$.

In Scikit-Learn wird in der logistischen Regression unterschieden zwischen predict und predict_proba. Den Wert $\hat{y} = w \cdot x + b$ selbst liefert decision_function, und damit ist

predict_proba gibt also eine Spalte pro Klasse zurück, in der Reihenfolge von classes_: bei Labels $0/1$ zuerst $P(Y = 0)$, dann $P(Y = 1)$. Die Wahrscheinlichkeit für Klasse 1 steht in der zweiten Spalte. Beim Softmax weiter unten begegnet uns dieselbe Form für $K$ Klassen.

Lineare Regression in PyTorch¶

In PyTorch macht man eine lineare Regression so:

import torch
import torch.nn as nn

class LinearRegression(nn.Module):
    def __init__(self, input_dim, output_dim):
        super().__init__()
        self.linear = nn.Linear(input_dim, output_dim)

    def forward(self, x):
        return self.linear(x)

Ein Modell ist eine Unterklasse von nn.Module; im Konstruktor werden die Bausteine mit Parametern angelegt (hier eine einzige lineare Schicht nn.Linear, das ist $x \mapsto Wx + b$), und forward beschreibt die Inferenz. Das Training:

import numpy as np

x = np.linspace(0, 10, 50)                      # irgendein eindimensionaler Datensatz
y = 2.5 * x + 1.0 + np.random.normal(0, 1, 50)  # verrauschte Gerade

model = LinearRegression(1, 1)
criterion = nn.MSELoss()
optimizer = torch.optim.SGD(model.parameters(), lr=0.01)

inputs = torch.from_numpy(x.reshape(-1, 1).astype("float32"))
labels = torch.from_numpy(y.reshape(-1, 1).astype("float32"))

number_of_epochs = 100
losses = np.zeros(number_of_epochs, dtype="float32")

for epoch in range(number_of_epochs):
    optimizer.zero_grad()
    outputs = model(inputs)
    loss = criterion(outputs, labels)
    loss.backward()
    optimizer.step()
    losses[epoch] = loss.item()

y_pred = model(inputs).detach().numpy()

Die fünf Zeilen in der Schleife sind die Trainingsschleife von PyTorch, die uns in allen weiteren Kapiteln begleitet: Gradienten auf null setzen, Vorwärtsrechnung, Loss berechnen, Gradienten rückwärts berechnen (backward), Parameter anpassen (step). Was criterion (die Loss-Funktion, hier der mittlere quadratische Fehler MSE) und optimizer (hier stochastischer Gradientenabstieg SGD mit Lernrate $0{,}01$) genau tun, klären die nächsten beiden Kapitel.

Logistische Regression in PyTorch¶

Für die logistische Regression ändern sich genau zwei Zeilen: Die Aktivierungsfunktion in forward

class LogisticRegression(nn.Module):
    def __init__(self, input_dim, output_dim):
        super().__init__()
        self.linear = nn.Linear(input_dim, output_dim)

    def forward(self, x):
        return torch.sigmoid(self.linear(x))  # <---

und die Loss-Funktion:

model = LogisticRegression(1, 1)
criterion = nn.BCELoss()  # <---
optimizer = torch.optim.SGD(model.parameters(), lr=0.01)

Der Rest der Trainingsschleife bleibt identisch. BCELoss ist die binäre Kreuzentropie: die erwartete Überraschung, wenn man mit den vorhergesagten Wahrscheinlichkeiten $(\sigma(\hat{y}), 1 - \sigma(\hat{y}))$ rechnet, die Wahrheit aber $(y, 1 - y)$ ist. Sie ist bis auf eine Konstante die Kullback-Leibler-Divergenz aus dem DS-Skript (dort $H(p \Vert r)$ geschrieben); den genauen Zusammenhang und woher diese Loss-Funktion kommt, klärt das nächste Kapitel.

Weg von der 2-Klassen-Gesellschaft: Softmax¶

Von Bernoulli zu Multinoulli¶

Von binomial (2 Klassen) zu multinomial ($K$ Klassen): Anstelle eines einzelnen $\hat{y}$-Werts der linearen Regression, der uns mit $\sigma$ die Wahrscheinlichkeit für Klasse 1 angibt, benötigen wir nun einen Vektor der Länge $K$, der sich auf 1 aufsummiert. Für $K = 2$ erwarten wir den Vektor $(\sigma(\hat{y}), 1 - \sigma(\hat{y}))$.

Wenn wir die Klasse one-hot kodieren, also Vektoren der Länge $K$ mit einer 1 an Stelle $k$ und sonst Nullen als Labels verwenden, dann liefert multinomiale lineare Regression – eine Matrix $W \in \mathbb{R}^{K \times n}$ statt eines Vektors $w$ – Vorhersagevektoren $v = Wx + b \in \mathbb{R}^K$. Die Entscheidung für Klasse $k \in \{1, \dots, K\}$ fällt, wenn der Eintrag des Vektors an Index $k$ am größten ist, $\operatorname{argmax}_j v_j = k$. Gesucht ist, analog zu $\sigma$, ein „soft argmax“.

Man darf sich das zweistufig vorstellen: Im ersten Schritt wird der Vektor $v$ komponentenweise exponenziert (alle Einträge werden positiv, und die Reihenfolge bleibt erhalten), im zweiten Schritt wird der Vektor normiert, sodass die Summe 1 ist. Für den Fall $K = 2$ mit $v = (x, y)$ können wir leicht einsehen, dass

Softmax verallgemeinert also tatsächlich das Sigmoid, wobei nur die Differenz der Logits zählt. Obacht beim Vergleich mit $\sigma$: In der Softmax-Formel steht kein Minus im Exponenten. Das Minus in $\sigma(x) = 1/(1 + e^{-x})$ entsteht erst beim Kürzen für $K = 2$, wie eben gerechnet. Die Probe, dass $\operatorname{sm}(v)$ wirklich eine Wahrscheinlichkeitsverteilung ist: Alle Einträge sind positiv, und die Summe ist $\sum_j e^{v_j} \big/ \sum_k e^{v_k} = 1$.

Übung¶

1. Zeigen Sie $\sigma(-x) = 1 - \sigma(x)$. Tipp: Erweitern Sie $\sigma(-x) = \frac{1}{1+e^{x}}$ mit $e^{-x}$.

2. Zeigen Sie $\sigma'(x) = \sigma(x)\,\sigma(-x)$. Tipp: Quotientenregel, dann mit $1/\sigma(-x)$ multiplizieren und so lange umformen, bis die Definition von $\sigma(x)$ dasteht. Mit Teil 1 ist das dasselbe wie die aus der Literatur bekannte Form $\sigma(x)(1-\sigma(x))$.

3. Wie groß ist $\sigma'(0)$? Und wie groß sind die Odds bei $p = 0{,}5$?

Antwortprobe: Die Odds sind 1. Die andere Zahl ist nach Teil 2 ein Produkt aus zwei gleichen Faktoren, also eine Quadratzahl; steht bei Ihnen $\tfrac12$, haben Sie $\sigma(0)$ statt $\sigma'(0)$ ausgerechnet.