⚠️ Dieser Text wurde aus dem Englischen Original maschinell übersetzt.

Das Verständnis der Einzelheiten einer logistischen Regression ist nicht trivial. Viele Quellen behandeln entweder nur die theoretische Seite oder die praktische Umsetzung. In diesem Beitrag möchte ich eine zentrale Anlaufstelle für die theoretischen Grundlagen und die praktische Anwendung der logistischen Regression schaffen.

Importierte Packages Link to heading

from functools import partial
from IPython import get_ipython

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pandas.io.formats.style

from sklearn.datasets import load_iris
from sklearn.linear_model import LogisticRegression

import warnings

%matplotlib inline

pd.set_option("display.precision", 2)
warnings.filterwarnings('ignore')

Ein bisschen Formatierung Link to heading

def add_raw_tag(df):
    return "
\n" + df.to_html(classes="nb-table") + "\n\n"

html_formatter = get_ipython().display_formatter.formatters['text/html']

html_formatter.for_type(
    pd.DataFrame,
    lambda df: add_raw_tag(df)
)

html_formatter.for_type(
    pd.io.formats.style.Styler,
    lambda df: add_raw_tag(df)
)

Der Datensatz Link to heading

Wir werden den bekannten Iris-Datensatz als Grundlage für unsere Diskussion verwenden. Er ist einfach genug, um noch etwas Intuition einzusetzen, wenn wir versuchen, die logistische Regression zu verstehen. Wenn wir den Datensatz von sklearn laden, erhalten wir 150 Beobachtungen von Messungen einiger Merkmale der Pflanzen. Dann haben Botaniker die Iris in drei Unterarten mit schicken lateinischen Namen eingeteilt, aber da ich kein Botaniker bin, werden wir einfach so tun, als könnten wir - basierend auf diesen Messungen - bestimmen, ob die Blume rot, grün oder blau blüht. Die Beobachtungen werden in $X$ gespeichert, wobei jede Beobachtung aus jeweils 4 Messungen besteht, und die Klassifikation landet im Vektor $y$, wobei die Einträge 0 = rot, 1 = grün, 2 = blau sind.

X, y = load_iris(return_X_y=True)
columns = ('Measurement 1', 'Measurement 2', 'Measurement 3', 'Measurement 4')
dataset = pd.DataFrame(X, columns=columns)
dataset["True Classification"] = y
dataset

Verkürzter Output:

Measurement 1 Measurement 2 Measurement 3 Measurement 4 True Classification
0 5.1 3.5 1.4 0.2 0
1 4.9 3.0 1.4 0.2 0
2 4.7 3.2 1.3 0.2 0
3 4.6 3.1 1.5 0.2 0
4 5.0 3.6 1.4 0.2 0
5 5.4 3.9 1.7 0.4 0
6 4.6 3.4 1.4 0.3 0
7 5.0 3.4 1.5 0.2 0
8 4.4 2.9 1.4 0.2 0
9 4.9 3.1 1.5 0.1 0
10 5.4 3.7 1.5 0.2 0
11 4.8 3.4 1.6 0.2 0
12 4.8 3.0 1.4 0.1 0
96 5.7 2.9 4.2 1.3 1
97 6.2 2.9 4.3 1.3 1
98 5.1 2.5 3.0 1.1 1
99 5.7 2.8 4.1 1.3 1
100 6.3 3.3 6.0 2.5 2
101 5.8 2.7 5.1 1.9 2
102 7.1 3.0 5.9 2.1 2
colors = ['red', 'green', 'blue']

for X_el, y_el in zip(X, y):
    plt.plot(X_el, colors[y_el], alpha=0.1)

plt.xticks(np.arange(4), columns)
plt.show()

png

Wir können deutlich erkennen, dass Messung 3 viele Informationen enthält. Sie trennt die roten Blumen deutlich von den anderen und scheint sogar eine gute Möglichkeit zu sein, die grünen von den blauen Blumen zu unterscheiden. Das Gleiche gilt für Messung 4, während Messung 1 und 2 anscheinend keine oder nur sehr wenig Hinweise darauf geben, welcher Klasse ein bestimmtes Exemplar angehört.

⚠️ scikit-learn benutzt ab Version 0.22 hier einen standardmäßig einen anderen Algorithmus. Die folgenden Berechnungen gelten unter Verwendung des Wertes ‘ovr’ für den Parameter multi_class. Die ‘multinomiale’ Variante erzeugt jedoch bessere Ergebnisse.

clf = LogisticRegression(random_state=42, multi_class="ovr").fit(X, y)

Nachdem das Modell trainiert wurde, ermöglicht uns die Funktion “predict”, eine Vorhersage für jedes Eingabe-Tupel zu erhalten. Die Klassen sind nummeriert und nicht benannt, können aber natürlich mithilfe der entsprechenden Zuordnung in Namen umgewandelt werden.

# 0 = red, 1 = green, 2 = blue
clf.predict(X)
array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
   0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1,
   1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 2, 1, 2, 1, 1,
   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2,
   2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
   2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2])

Neben einer einfachen Klassenprognose können wir auch die Wahrscheinlichkeit für jede Klasse für jedes Eingabedatensample erhalten.

clf.predict_proba(X)[0:5]  # abbreviated output
array([[8.96807569e-01, 1.03191359e-01, 1.07219602e-06],
   [7.78979389e-01, 2.21019299e-01, 1.31168933e-06],
   [8.34864184e-01, 1.65134802e-01, 1.01485082e-06],
   [7.90001986e-01, 2.09996107e-01, 1.90723705e-06],
   [9.12050403e-01, 8.79485212e-02, 1.07537143e-06]])

Für das erste oben genannte Beispiel weisen wir eine Wahrscheinlichkeit von 98% zu, dass es zur Klasse ‘0’ gehört, eine Wahrscheinlichkeit von 2%, dass es zur Klasse ‘1’ gehört, und eine Wahrscheinlichkeit von fast 0%, dass es zur Klasse ‘2’ gehört. Mit der Funktion ‘argmax’ könnten wir dann die Wahrscheinlichkeiten auf das genaue Ergebnis der oben genannten ‘predict’-Funktion abbilden.

Für den ersten Teil haben wir also den Datensatz untersucht und eine Vorstellung davon bekommen, wie gut er in der Lage sein wird, die Klassen auf der Grundlage der gegebenen Merkmale zu trennen. Außerdem haben wir das Ergebnis der logistischen Regression von sklearn generiert, die durch ein Training auf dem vollständigen Datensatz bereitgestellt wird, und gelernt, wie man es interpretiert.

Bewertung des Modells Link to heading

Ein gängiges Maß zur Bewertung der Qualität von Vorhersagen ist der ‘Score’, der gemäß der sklearn-Hilfe definiert ist als:

In der Mehrklassenklassifikation handelt es sich dabei um die Teilmengengen-Genauigkeit (subset accuracy), was ein strenges Maß ist, da für jedes Sample verlangt wird, dass jeder Label-Satz korrekt vorhergesagt wird.

clf.score(X, y)

0.9533333333333334

Wir zählen einfach alle korrekt klassifizierten Labels und teilen diese durch die Gesamtanzahl:

correct = []
for y_true, y_pred in zip(y, clf.predict(X)):
    if y_true == y_pred:
        correct.append(1)
    else:
        correct.append(0)

sum(correct) / len(correct)

0.9533333333333334

Da wir den vollständigen Datensatz für das Training verwendet haben, ist das in Ordnung. Ein detaillierterer Ansatz zur Bewertung der Qualität unseres Modells basiert auf dem Klassifikationsbericht (classification report).

from sklearn.metrics import classification_report

print(classification_report(y_pred=clf.predict(X), y_true=y))
          precision    recall  f1-score   support

       0       1.00      1.00      1.00        50
       1       0.96      0.90      0.93        50
       2       0.91      0.96      0.93        50

accuracy                           0.95       150

macro avg 0.95 0.95 0.95 150 weighted avg 0.95 0.95 0.95 150

Während die Klasse “0 = rot” perfekte ‘precision’ und ‘recall’ aufweist, da die Klassen “1 = grün” und “2 = blau” in ihren Merkmalen teilweise überlappen, verwechselt der Algorithmus sie. Dies wird durch den ‘recall’ von 0,9 für die Klasse “1” gezeigt, was bedeutet, dass nur 90% aller grünen Blumen als grün klassifiziert wurden, und durch die ‘precision’ der Klasse “2”, was bedeutet, dass nur 91% aller als blau klassifizierten Blumen tatsächlich blau waren.

y_true_with_jitter = y + np.random.rand(y.shape[0]) * 0.25
y_classified_with_jitter = clf.predict(X) + np.random.rand(y.shape[0]) * 0.25
plt.xticks(np.arange(3), ('true red', 'true green', 'true blue'))
plt.yticks(np.arange(3), ('class. red', 'class. green', 'class. blue'))
plt.scatter(y_true_with_jitter,
            y_classified_with_jitter,
            color=[colors[y_el] for y_el in y])
plt.show()

png

Ist es ein gutes Modell? Nun, einerseits sind die reinen Werte unserer Statistiken wirklich ziemlich gut, andererseits haben wir einen winzigen Datensatz, und wir haben auf dem gesamten Datensatz trainiert, da wir keinen separaten Datensatz für Tests beiseite gelegt haben. Da wir jedoch nicht allzu viel Wert auf das tatsächliche Ergebnis der korrekten Vorhersage der Blumenfarben legen und nur verstehen wollen, wie wir zu den Vorhersagen kommen, werden wir es akzeptieren. Machen Sie sich jedoch mental darauf gefasst, dass dies in einer realen Anwendung kein gutes Ergebnis wäre.

Ein Blick auf die inneren Abläufe der sklearn Logistic Regression Link to heading

Wenn wir uns die Variablen des Log-Regression-Klassifikators nach dem Training ansehen, finden wir drei Sätze von Koeffizienten und drei verschiedene Intercepts. Das liegt daran, dass die Log-Regression im Wesentlichen binär ist, d.h. nur eine Ja/Nein- oder 1/0-Klassifikation durchführt. Wenn wir $n > 2$ Klassen haben, müssen wir dieses Problem in $n$ separate “1 gegen Rest”-Klassifikationsprobleme aufteilen. Jeder Satz von Koeffizienten und jeder Intercept gehört zu einer dieser Teil-Klassifikationen.

clf.__dict__["coef_"] # just running clf.__dict__ spits out all the info about the trained model

array([[-0.44501376, 0.89999242, -2.32353827, -0.97345836], [-0.1792787 , -2.12866718, 0.69665417, -1.27480129], [-0.39444787, -0.5133412 , 2.93087523, 2.41709879]])

clf.__dict__["intercept_"]
array([  6.69040651,   5.58615272, -14.43121671])

Wenn wir nun einen Satz von Merkmalen $x^i$ in den trainierten Klassifikator einspeisen, können wir die Wahrscheinlichkeiten berechnen, dass $x^i$ zu einer Klasse gehört oder nicht, und zwar über:

$$p(x_i)=\frac{1}{1+e^{-(\beta_{0} + \beta_{1}x_1^{i} + \beta_{2}x_2^{i} + \beta_{3}x_3^{i} + \beta_{4}x_4^{i})}}$$

wobei $\beta_{0}$ ein Intercept ist und $\beta_{1}..\beta_{4}$ die Koeffizienten für jede Eintrag im Merkmalsvektor $x^i = (x^{i}_1,x^{i}_2,x^{i}_3,x^{i}_4)$ sind. Wir werden später erkunden, warum dieser Ausdruck der richtige ist. Lassen Sie uns das oben für unsere allererste Beobachtung im Datensatz berechnen:

# @ is shorthand for matrix multiplication
p_0 = 1 / (1 + np.exp(-(clf.intercept_[0] + X[0] @ clf.coef_[0])))
p_1 = 1 / (1 + np.exp(-(clf.intercept_[1] + X[0] @ clf.coef_[1])))
p_2 = 1 / (1 + np.exp(-(clf.intercept_[2] + X[0] @ clf.coef_[2])))
print('p_0 =', p_0)
print('p_1 =', p_1)
print('p_2 =', p_2)

p_0 = 0.9840648521799207 p_1 = 0.11323163760159996 p_2 = 1.1765181958016876e-06

Mit dieser Berechnung haben wir nun festgestellt, dass unser $x^0$ eine Wahrscheinlichkeit von $0.98$ hat, zur Klasse “0 = rot” zu gehören, im Vergleich zu jeder der anderen Klassen - entweder grün oder blau. Wie wir jedoch sehen können, addieren sich diese drei Wahrscheinlichkeiten nicht zu 100%, und warum sollten sie das auch? Diese Wahrscheinlichkeiten gehören zu drei mathematisch unabhängigen Problemen:

Gehört $x^i$ zur Klasse 0 oder nicht? Gehört $x^i$ zur Klasse 1 oder nicht? Gehört $x^i$ zur Klasse 2 oder nicht?

Was passiert jedoch, wenn wir diese Wahrscheinlichkeiten linear skalieren, damit sie sich zu 1 addieren?

p_sum = p_0 + p_1 + p_2
print("p_0_scaled =", p_0 / p_sum)
print("p_1_scaled =", p_1 / p_sum)
print("p_2_scaled =", p_2 / p_sum)

p_0_scaled = 0.896807568632095 p_1_scaled = 0.10319135917188017 p_2_scaled = 1.0721960247752252e-06

Wir haben diese genauen Zahlen bereits gesehen und können unsere Vorhersage mit argmax treffen:

clf.predict_proba(X)[0]

array([8.96807569e-01, 1.03191359e-01, 1.07219602e-06])

np.argmax(clf.predict_proba(X)[0]) # 0 = red, 1 = green, 2 = blue
0

Jetzt haben wir ein Verständnis dafür, wie die von sklearn generierten Daten des Trainingsprozesses interpretiert werden können, und wir haben über die clf.predict-Funktion hinausgeschaut, um zu verstehen, wie die Vorhersagen im trainierten Modell ausgewählt werden.

Der mathematische Hintergrund Link to heading

Wir haben oben die Formel verwendet:

$$p(x^i)=\frac{1}{1+e^{-(\beta_{0} + \beta_{1}x_{1}^{i} + \beta_{2}x_{2}^{i} + \beta_{3}x_{3}^{i} + \beta_{4}x_{4}^{i})}}$$

Diese Formel ist technisch gesehen eine Wahl, keine mathematische Vorgabe (tatsächlich gibt es andere Formeln, die ebenfalls funktionieren), aber warum ergibt sie Sinn?

Um das zu verstehen, müssen wir zunächst die Wahrscheinlichkeiten und ihre Beziehung zu den Odds verstehen. Wenn ein Ereignis eine 50%ige Chance hat, einzutreten, sind die Odds dafür 1:1 (1 Mal tritt es ein, 1 Mal tritt es nicht ein). Wenn ein Ereignis eine 33,3%ige Chance hat, einzutreten, sind die Odds 1:2 (1 Mal tritt es ein, 2 Mal tritt es nicht ein), 25% repräsentieren Odds von 1:3 und 20% repräsentieren Odds von 1:4. Oder als allgemeine Formel:

$$ Odds = \frac{p}{1-p} $$

z.B. für $p = 0,25$ ergibt das $ Odds = \frac{0,25}{1-0,25} = 0,333… $, d.h. für jede 3-mal, dass das betrachtete Ereignis nicht eintritt, wird es 1-mal eintreten oder 1 Erfolg : 3 Misserfolge.

Wenn wir das $p$ in der Odds-Formel durch das $p(x^i)$ von oben ersetzen, erhalten wir:

$$ Odds = e^{\beta_{0} + \beta_{1}x_{1}^{i} + \beta_{2}x_{2}^{i} + \beta_{3}x_{3}^{i} + \beta_{4}x_{4}^{i}} $$

oder

$$ \log(Odds) = \beta_{0} + \beta_{1}x_{1}^{i} + \beta_{2}x_{2}^{i} + \beta_{3}x_{3}^{i} + \beta_{4}x_{4}^{i} $$

Damit haben wir eine Verbindung zu unserem Merkmalsraum $X$ hergestellt und können jede Beobachtung von Merkmalen auf eine Wahrscheinlichkeit $ p \in (0,1) $ abbilden. Alles, was wir dann brauchen, ist eine mehr oder weniger willkürliche Grenzregel, normalerweise $ p>0,5 $.

Außerdem erhalten wir die Interpretierbarkeit quasi kostenlos: Der Koeffizient $\beta_i$ beschreibt die Änderung der Odds, wenn wir $x_i$ um eine Einheit erhöhen. Schauen Sie sich noch einmal das Spaghetti-Diagramm mit den Farben oben an. Je größer der Wert der Messung 3, desto geringer ist die Wahrscheinlichkeit, dass wir ein rotes Exemplar in der Hand haben. Tatsächlich beträgt der $\beta$-Koeffizient der Messung 3 in unserem “Rot vs. Rest”-Problem -2,26, was bedeutet, dass eine Erhöhung der Messung 3 um eine Einheit die Odds des Exemplars, rot zu sein, um $exp(-2.26) \approx 0,104$ verringert, was sehr grob 1 zu 9,5 entspricht.

Lassen Sie uns diese Funktion $p(x)$ zeichnen, um zu sehen, wie sie aussieht:

import matplotlib.pyplot as plt


def map_to_p(log_odds):
    odds = np.exp(log_odds)
    p = odds / (1 + odds)
    return p


lots_of_log_odds_for_drawing = np.linspace(-10, 10, num=1000)
mapped_to_p = list(map(map_to_p, lots_of_log_odds_for_drawing))

plt.xlabel("log(Odds)")
plt.ylabel("p")
plt.plot(lots_of_log_odds_for_drawing, mapped_to_p);

png

Diese Funktion wird als Sigmoid-Funktion bezeichnet, die - in Bezug auf unser logistisches Regressionsmodell - die sogenannte Verbindungsfunktion ist, da sie einen Prädiktor, die lineare Kombination $\log(Odds)$, mit einer Antwort in $p \in (0,1)$ verknüpft.

Anpassung der Parameter Link to heading

Die Anpassung der Parameter ist etwas knifflig, da wir keine Methode der kleinsten Quadrate wie im linearen Fall verwenden können. Wir müssen numerische Methoden wie die Methode der maximalen Likelihood (MLE) verwenden. Intuitiv möchten wir, dass unsere $\beta$-Werte so gewählt werden, dass die linearen Kombinationen mit den Merkmalsvektoren $x$ so weit wie möglich von der Mitte des Sigmoids und so weit wie möglich auf einer Seite für “Erfolge” (normalerweise die positiven) und auf der anderen Seite für “Misserfolge” liegen. “Erfolg” bedeutet, Mitglied einer bestimmten Klasse $i$ zu sein, und Misserfolg bedeutet eine höhere Wahrscheinlichkeit für jede andere Klasse.

Dieses Beta, sobald wir es gefunden haben, werden wir $\hat\beta$ nennen.

Schauen wir uns das Sigmoid noch einmal an und verwenden die Parameter für die Klasse “0 = Rot” aus der logistischen Regression von sklearn, um unsere $x$-Werte zu verteilen:

colors = ['red', 'green', 'blue']

# this is the result of b_0 + b_1 * x_1 + b_2 * x_2 + ...
# clf is still our trained classifier
log_odds_from_our_dataset = [clf.intercept_[0] + x @ clf.coef_[0] for x in X]

plt.xlabel("log(Odds)")
plt.ylabel("p")
plt.plot(lots_of_log_odds_for_drawing, mapped_to_p)

for x, y_cl in zip(log_odds_from_our_dataset, y):
    #plt.scatter(x, 0, s=10, c=colors[y_cl])
    plt.scatter(x, map_to_p(x), s=10, c=colors[y_cl])
plt.show()

png

Jetzt sehen wir, dass alle roten Punkte, die die Mitglieder der Klasse $0$ repräsentieren, auf der rechten Seite von 0 liegen (und daher $p > .5$ haben), und die anderen beiden Klassen auf der linken Seite liegen. Kurz gesagt müssen wir $\hat\beta$ so wählen, dass die Sigmoid-Funktion für einige $x$-Werte nahe bei 0 und für andere nahe bei 1 liegt. Dies ist ein sehr umständliches Problem. Glücklicherweise arbeiten wir mit dem Intervall $p \in (0,1)$, was bedeutet, dass wir wissen, dass das Maximum 1 ist. Daher können wir die $x$-Werte, die zu einem Wert von $p$ nahe 0 führen sollen, umdrehen, indem wir $(1-p)$ berechnen. Jetzt haben wir ein Maximierungsproblem für alle $x$ in unserem Bereich.

Auch ein Wort zur Konstanten (Intercept): Während $\beta_1 .. \beta_4$ im Wesentlichen unsere $x$-Werte so “dehnen”, dass es so wenig Überlappung zwischen den Gruppen wie möglich gibt, ändert $\beta_0$ die Position des gesamten Punktesatzes, sodass er schön um 0 zentriert ist. Tatsächlich müssen wir dem Intercept keine besondere Behandlung geben, da wir unseren Merkmalsvektor $x$ einfach um eine konstante 1 erweitern können, wie z.B. $x^i = (1, x_{1}^{i}, x_{2}^{i},x_{3}^{i},x_{4}^{i})$. Die Verwendung dieser erweiterten Vektoren in den folgenden Schritten wird die Dinge erheblich vereinfachen.

Darüber hinaus, da wir alle unsere individuellen $p$ und $1-p$ maximieren möchten, können wir auch versuchen, das Produkt all dieser Werte zu maximieren. Da wir möchten, dass alle Terme des Produkts möglichst nahe bei 1 liegen, möchten wir, dass das gesamte Produkt möglichst nahe bei 1 liegt.

In unseren Eingabedaten wird die Klasse, zu der eine bestimmte Probe gehört, durch $y_i \in {0,1,2}$ bezeichnet, wobei diese Zahlen Rot, Grün und Blau repräsentieren. Wenn wir das in ein binäres Problem übersetzen möchten, benötigen wir ein $y_{binary, i} \in {1,0}$, wobei z.B. für die erste Klassifikation “rot vs. nicht rot” ein Erfolg mit “1” (Blume ist rot) und ein Misserfolg mit “0” (Blume ist nicht rot) bezeichnet wird. Ich werde die Bezeichnung “binär” der Kürze halber weglassen, aber bitte merken Sie sich, dass die $y_i$s, mit denen wir von nun an arbeiten, nicht mehr dieselben wie oben sind.

Mathematisch ausgedrückt möchten wir unseren Schätzer $\hat\beta$ finden, der die Likelihood-Funktion maximiert:

$$ l(\beta) = \prod_{x_i; y_i = 1} p(x_i) \times \prod_{x_i; y_i = 0} (1-p(x_i)) $$

oder:

$$ \hat\beta = \arg\max_{\beta} l(\beta) $$

was folgendermaßen vereinfacht werden kann:

$$ l(\beta) = \prod_{i} p(x_i)^{y_i} (1-p(x_i))^{1-y_i}$$

Wir müssen jetzt das nächste Konzept einführen - um dieses Maximierungsproblem anzugehen, nutzen wir die Tatsache, dass $log(a)$ monoton mit $a$ steigt, daher ist die Maximierung von $log(a)$ äquivalent zur Maximierung von $a$, daher:

$$ \hat\beta = \arg\max_{\beta} l(\beta) \iff \hat\beta = \arg\max_{\beta} \log l(\beta) $$

unter Verwendung dieser Eigenschaft können wir die Multiplikation in $l(\beta)$ zu einer Summ in $\log l(\beta)$ (oder kurz $ll(\beta)$) umwandeln.

$$ ll(\beta) = \sum_{i} y_i \log(p(x^i)) + (1-y_i)\log(1-p(x^i))$$

Wir werden diese Gleichung jetzt weiter vereinfachen, beginnend mit der Log-Wahrscheinlichkeit im ersten Term (die blaue Hervorhebung wird weiter unten deutlicher werden):

$$ y_i\log(p(x^i)) = y_i\log\frac{1}{1+e^{-\beta x^i}} = \color{blue}{-y_i\log(1+e^{-\beta x^i})} $$

Während das ziemlich einfach war, ist der zweite Term etwas anspruchsvoller. Für den Moment werden wir den Term $(1-y_i)$ weglassen und uns auf die logarithmierte inverse Wahrscheinlichkeit im zweiten Term konzentrieren:

$$ \log(1-p(x^i)) = \log(1-\frac{1}{1+e^{-\beta x^i}}) = \log(\frac{e^{-\beta x^i}}{1+e^{-\beta x^i}})$$

Um fortzufahren, müssen wir uns entscheiden, was wir mit dem Term in den Klammern tun wollen: Wir können entweder das $e^{-\beta x_i}$ im Zähler nach unten in den Nenner bringen oder die Logarithmusregel für Brüche verwenden, um den Bruch in eine Subtraktion von zwei logarithmusfreien Logarithmen aufzuteilen. Es stellt sich heraus, dass wir beides tun müssen, um zur einfachsten möglichen Form der kombinierten Gleichung zu gelangen. Nachdem wir die Terme aus dem $(1-y_i)$ verteilt haben, müssen wir jeden mit einer anderen Strategie behandeln. Für den $1$-Term (beachten Sie, dass das Minuszeichen im Nenner verschwunden ist):

$$ 1 \log(\frac{e^{-\beta x^i}}{1+e^{-\beta x^i}}) = \log(\frac{1}{e^{\beta x^i}(1+e^{-\beta x^i})}) = \log(\frac{1}{e^{\beta x^i}+1}) = \color{green}{ -\log(e^{\beta x^i}+1)} $$

und für den $-y_i$ Term:

$$ -y_i \log(\frac{e^{-\beta x^i}}{1+e^{-\beta x^i}}) = -y_i (\log(e^{-\beta x^i}) - \log(1+e^{-\beta x^i})) $$

$$ = \color{green}{ -y_i(-\beta x^i)} \color{blue}{ + y_i \log(1+e^{-\beta x^i})} $$

Wenn wir die Puzzlestücke nun wieder zusammenbauen, eliminieren sich die blauen Terme und die grünen bleiben übrig:

$$ ll(\beta) = \sum_{i} \color{blue}{-y_i\log(1+e^{-\beta x^i})} \color{green}{ -\log(e^{\beta x^i}+1)} \color{green}{ -y_i(-\beta x^i)} \color{blue}{ + y_i \log(1+e^{-\beta x^i})} $$

$$ = \sum_{i} y_i\beta x^i-\log(e^{\beta x^i}+1) $$

Um unser Gedächtnis aufzufrischen:

  • $y_i \in {1,0}$, was bedeutet, dass ein $x^i$ zu einer bestimmten Klasse (“1”) gehört oder nicht (“0”).
  • $x^i \in \mathbb{R}^5$, wobei das erste Element auf 1 festgelegt ist, und den Feature-Vektor einer Beobachtung repräsentiert.
  • $\beta \in \mathbb{R}^5$, der Vektor der Koeffizienten, wobei der erste Eintrag das Intercept repräsentiert.

Wir haben es jetzt geschafft, unser Optimierungsproblem in vergleichsweise einfachen Begriffen darzustellen. Das Einzige, was jetzt noch fehlt, ist das $\beta$, das den obigen Ausdruck maximieren wird, aber alle anderen Variablen sind klar definiert. Wir können das optimale $\beta$ jedoch nicht algebraisch berechnen und müssen auf numerische Methoden zurückgreifen.

Ein einfaches Beispiel Link to heading

Um uns auf die nächsten Schritte vorzubereiten, in denen wir tatsächlich die $\beta$-Koeffizienten anpassen, müssen wir die theoretische Mathematik in Python-Code übersetzen. Außerdem ist es notwendig, das Problem mit drei Klassen von roten, grünen und blauen Blumen in mehrere binäre Probleme wie “Blume ist rot vs. Blume ist nicht rot” zu übersetzen. Zuerst die Übersetzung der Log-Likelihood-Funktion.

def log_likelihood(x, y, beta):
    ll = 0
    for x_el, y_el in zip(x, y):
        ll += y_el * (beta @ x_el) - np.log(np.exp(beta @ x_el) + 1)

    return ll

Wir teilen das Problem mit drei Klassen in 3 binäre Teilprobleme auf. Daher müssen wir den Klassenvektor so ändern, dass $y$ nur “1”-Einträge für die einzelne Klasse hat, die wir testen, und “0”-Einträge für die anderen beiden Klassen:

y_binary_red = [1 if y_el == True else 0 for y_el in y == 0]
y_binary_green = [1 if y_el == True else 0 for y_el in y == 1]
y_binary_blue = [1 if y_el == True else 0 for y_el in y == 2]

Außerdem müssen wir eine “1” am Anfang jedes Merkmalsvektors $x$ hinzufügen, um den Intercept zu berücksichtigen.

# make a vector with just "1"s and glue it to the left side of X
X_with_intercept = np.hstack((np.ones((X.shape[0], 1)), X))

Jetzt verlassen wir die Komfortzone der algebraischen Gewissheit und müssen numerische Methoden verwenden, um unsere Schätzung $\hat\beta$ zu erhalten. Bevor wir das tun, wollen wir sehen, was wir im Prinzip mit einem Brute-Force-Algorithmus erreichen möchten. Wir beginnen mit zufälligen $\beta$-Werten - hier habe ich ein wenig geschummelt, da ich bereits grob weiß, in welchem Bereich die Variablen liegen sollten. Mit diesem halb zufälligen $\beta$ berechnen wir die Log-Likelihood-Funktion, die Werte zwischen $-\infty$ und $0$ annehmen kann, da es sich um einen $log$ einer Wahrscheinlichkeit $p \in (0,1)$ handelt. Wenn das zufällige $\beta$ unsere Wahrscheinlichkeit erhöht, behalten wir es, andernfalls werfen wir es weg und wählen ein anderes zufälliges $\beta$. Um die Ergebnisse zu visualisieren, programmieren wir zuerst eine kleine Hilfsfunktion.

# define a plot function to visualize the result
def plot_separation(x, beta, y, color='red'):
    color = ['grey', color]

    for x_el, y_el in zip(x, y):
        log_odds = beta @ x_el
        plt.scatter(log_odds, map_to_p(beta @ x_el), c=color[y_el])
    plt.plot(lots_of_log_odds_for_drawing, mapped_to_p)
    plt.show()
# choose a random, but very small likelihood as basis
ll_hat = -1000000

for step in range(10001):
    # choose "random" beta vector 10.000 times 
    # each entry will be between -3 and 3
    beta_random = 6 * np.random.random(5) - 3

    ll = log_likelihood(X_with_intercept, y_binary_red, beta_random)

    # if our log-likelihood has improved, overwrite old beta, save likelihood for futher iterations
    if ll > ll_hat:
        beta_hat = beta_random
        ll_hat = ll

    # draw the result every 5000 steps
    if step % 5000 == 0:
        print("Step:", step, ", beta_hat:", beta_hat, ", ll_hat:", ll_hat)
        plot_separation(X_with_intercept, beta_hat, y_binary_red)

        print()
Step: 0 , beta_hat: [-0.06781811  0.82358851 -0.58925447  2.8207669  -1.35588074] , ll_hat: -1496.509430016339

png

Step: 5000 , beta_hat: [-2.85756301  0.74906527  2.64267461 -2.84641646 -2.83370022] , ll_hat: -0.6742826828096209

png

Step: 10000 , beta_hat: [-2.85756301  0.74906527  2.64267461 -2.84641646 -2.83370022] , ll_hat: -0.6742826828096209

png

Wie wir sehen können, trennt dies die roten von den nicht-roten Punkten sauber auf die “1” und “0” Seite der Sigmoid-Kurve. Um ehrlich zu sein, scheint das nur deshalb zu funktionieren, weil die roten Exemplare von Anfang an ziemlich sauber vom Rest getrennt sind. Wenn Sie die anderen Farben ausprobieren, werden die Ergebnisse nicht so gut sein.

Aber wir haben das allgemeine Prinzip gesehen. Indem wir wild zufällige $\beta$-Werte wählen und diejenigen behalten, die die Wahrscheinlichkeit erhöhen, verschieben wir die Log-Odds unserer “Erfolgs”-Klasse so weit wie möglich nach rechts, während wir die “Fehlschläge” auf der linken Seite behalten. Nun bleibt nur noch der Übergang von zufälligem Raten zu einem anspruchsvolleren Prozess.

Anpassung der Parameter mit Gradient Descent Link to heading

Als nächster Schritt werden wir die Brute-Force-Methode durch eine numerische Optimierungsmethode wie den Gradient Descent oder den Newton-Raphson-Algorithmus ersetzen. Heute verwenden wir die Methode des Gradient Descents. Vereinfacht ausgedrückt werden wir $\beta$ in kleinen Schritten in Richtung der Minimierung unserer Fehlerfunktion bewegen, die den tatsächlichen $y$-Wert minus unser berechnetes Ergebnis für $y$ unter unserer Vermutung für $\beta$ ist. Diese Richtung ist zufällig der negative Gradient. Das tatsächliche $y$ ist einfach unsere Eingabedaten für $y$.

def gradient_descent(X, y, steps, learning_rate):

    beta = np.zeros(X.shape[1])

    for _ in range(1, steps):

        # calculate log odds for all x in X at once
        log_odds = np.dot(X, beta)

        # calculate result based on current beta
        tentative_y = list(map(map_to_p, log_odds))

        # calculate difference between current estimate and truth
        error = np.subtract(y, tentative_y)

        # see below for explanation
        gradient = np.dot(error, X)

        # move beta in opposite direction from error
        beta += learning_rate * gradient

    return beta

Also, was ist der Gradient von $ll(\beta)$ in Bezug auf $\beta$?

$$\nabla_{\beta} ll(\beta) = \sum_{i} \nabla_{\beta} y_i\beta x_i-\nabla_{\beta}\log(e^{\beta x_i}+1) $$ $$= \sum_{i} y_i x_i-\nabla_{\beta}\log(e^{\beta x_i}+1) $$ $$= \sum_{i} y_i x_i - x_i e^{\beta x_i} \frac{1}{e^{\beta x_i} +1} $$ $$= \sum_{i} y_i x_i - x_i \frac{1}{1 + e^{-\beta x_i}} $$ $$= \sum_{i} y_i x_i - x_i p(x_i) = \sum_{i} (y_i -p(x_i)) x_i $$ Dies ist nichts anderes als das wahre $y_i$ minus die berechnete Näherung für $y_i$, die $p(x_i)$ mal den Merkmalsvektor für jede Beobachtung ist, oder in Matrixform:

$$\nabla_{\beta} ll(\beta) = (y_{true} - y_{estimate}(\beta))X $$ Um die Dinge ein wenig zu ändern, versuchen wir, mit der Methode des Gradientenabstiegs die blauen Exemplare anstelle der roten zu identifizieren:

beta_hat = gradient_descent(X_with_intercept,
                            y_binary_blue,
                            steps=100001,
                            learning_rate=10e-5)
plot_separation(X_with_intercept, beta_hat, y_binary_blue, color='blue')

png

Und da haben wir eine Trennung der blauen Punkte in Richtung der positiven reellen Zahlen und der Rest in Richtung der negativen Zahlen und ihrer jeweiligen Wahrscheinlichkeiten, die sich 1 und 0 nähern. Wie wir jedoch sehen können, funktioniert die Trennung nicht so gut wie bei der “Rot gegen Rest” -Aufgabe.

Wenn wir jedoch dasselbe mit den grünen Exemplaren versuchen, funktioniert es überhaupt nicht gut. Das war jedoch irgendwie zu erwarten, wie wir am Anfang gesehen haben. Wir sollten uns damit trösten, dass auch die sklearn-Implementierung nicht sehr gut abschneidet, wie aus der Verwirrungsmatrix hervorgeht. Wenn wir das Ergebnis für die Teilaufgabe “Grün gegen Rest” plotten, können wir kaum zwischen den grünen und den grauen Punkten unterscheiden.

# parameters for green from sklearn log-reg
beta_0 = np.hstack((clf.__dict__["intercept_"][0], clf.__dict__["coef_"][0]))  # red
beta_1 = np.hstack((clf.__dict__["intercept_"][1], clf.__dict__["coef_"][1]))  # green
beta_2 = np.hstack((clf.__dict__["intercept_"][2], clf.__dict__["coef_"][2]))  # blue
plot_separation(X_with_intercept, beta_0, y_binary_red, color='red')

png

plot_separation(X_with_intercept, beta_1, y_binary_green, color='green')

png

plot_separation(X_with_intercept, beta_2, y_binary_blue, color='blue')

png

Abschließend können wir feststellen, dass unsere Lösung zwar nicht so gut ist wie die in sklearn implementierte Version, aber bereits sehr nahe an die Ergebnisse heranreicht. Ein Unterschied zwischen den beiden Algorithmen besteht z.B. darin, dass sklearn Lösungen mit großen Koeffizienten in seinem Optimierer bestraft.

Sources: