Woher eine Verlustfunktion kommt: Likelihood, nicht Konvention

Hier sind die Messwerte und die drei Linien, auf drei Arten bewertet: quadratischer Fehler, zu dem alle greifen; absoluter Fehler, zu dem vielleicht ein Statistiker greifen würde; und größter Fehler, zu dem der Maschinist greifen würde, weil der Prüfer sich nicht für deinen Durchschnitt interessiert — er weist das eine Teil zurück, das außerhalb der Toleranz liegt.

NumPy kommt hier dazu, ein Kapitel nach dem reinen Python-Perzeptron, aus einem Grund: Am Ende dieses Kapitels bewerten wir vierhunderttausend Kandidatenlinien gegen jeweils zwanzig Messwerte, und eine Python-Schleife ist dafür das falsche Werkzeug. Es ist außerdem die Notation, in der jede unten zitierte Quelle geschrieben ist.

import numpy as np

# Hours since the blade was changed, and the width of the part measured then.
SHIFT = np.array([
    (0.5, 20.17), (1.0, 20.28), (1.5, 20.53), (2.0, 20.61), (2.5, 20.69),
    (3.0, 20.94), (3.5, 21.21), (4.0, 21.31), (4.5, 21.27), (5.0, 21.35),
    (5.5, 21.58), (6.0, 21.80), (6.5, 21.67), (7.0, 22.07), (7.5, 22.10),
    (8.0, 22.31), (8.5, 22.48), (9.0, 22.66), (9.5, 22.90), (10.0, 23.13),
])
h, y = SHIFT[:, 0], SHIFT[:, 1]

LINES = {"A": (20.10, 0.26), "B": (20.20, 0.28), "C": (20.30, 0.26)}

for name, (a, b) in LINES.items():
    r = y - (a + b * h)                                      
    print(f"{name}   mean square {np.mean(r**2):.5f}"
          f"   mean absolute {np.mean(np.abs(r)):.5f}"
          f"   worst {np.max(np.abs(r)):.3f}")               

Die Größe in den hervorgehobenen Zeilen ist das Residuum: was die Linie sagte minus dem, was die Schieblehre sagte, eine Zahl pro Messwert. Jede Bewertungsregel in diesem Kapitel und jede Verlustfunktion in den achtundzwanzig Kapiteln danach ist irgendeine Art, eine Liste von Residuen auf eine einzige Zahl zusammenzudrücken. Sie unterscheiden sich nur darin, wie sie drücken.

A   mean square 0.02699   mean absolute 0.12600   worst 0.430
B   mean square 0.02524   mean absolute 0.13700   worst 0.350
C   mean square 0.03179   mean absolute 0.15000   worst 0.320

Lies die Spalten, nicht die Zeilen. Quadratischer Fehler sagt B, absoluter Fehler sagt A, größter Fehler sagt C: drei Regeln, drei Sieger, auf denselben zwanzig Punkten.

Ich habe diese drei Linien so gewählt, dass sie widersprechen, und das sollte ich offen sagen. Der Punkt ist, wie leicht das war — ein paar Minuten Suche über vernünftig aussehende Achsenabschnitte und Steigungen liefern Hunderte solcher Tripel. Die Rangfolge ist eine Eigenschaft der Regel, die du gewählt hast, keine Tatsache über die Linien. Deshalb ist die Regel kein Implementierungsdetail: Sie ist die Definition des Problems. Damit stellt sich die Frage, für deren Antwort dieses Kapitel existiert: Auf welcher Grundlage wählst du sie?

Zuerst eine kleinere Sache, denn es gibt nicht drei Linien, sondern unendlich viele. Nimm vorerst den quadratischen Fehler, weil ihn alle nehmen, und schrumpfe das Problem auf eine einzige Zahl mit dem Trick, der dem Perzeptron in Kapitel 1 elftausend Epochen erspart hat: Ziehe von beiden Spalten den Mittelwert ab. Sobald die Punktwolke um den Ursprung zentriert ist, geht die beste Linie unter quadratischem Fehler exakt durch den Ursprung — der Achsenabschnitt ist also erledigt, und nur die Steigung bleibt zu wählen.

u, v = h - h.mean(), y - y.mean()        # 5.25 hours, 21.553 mm

def mse(theta):
    return np.mean((v - theta * u) ** 2)

grid = np.arange(0.0, 0.6001, 0.001)
curve = np.array([mse(t) for t in grid])
print(grid.size, "candidates ->", f"theta={grid[curve.argmin()]:.3f}", f"mse={curve.min():.6f}")

601 candidates -> theta=0.293 mse=0.010115

Sechshunderteins Kandidatensteigungen, ein Sieger: 0,293 Millimeter pro Stunde gegenüber einer Wahrheit von 0,300. Zwanzig verrauschte Messwerte und eine for-Schleife kamen bis auf ein Hundertstel Millimeter pro Stunde heran — zweiein Drittel Prozent.

Der interessante Teil ist nicht der Sieger, sondern die Form der Suche. Gib die ganze Kurve aus, gedreht, sodass der Verlust von links nach rechts läuft:

ts = np.arange(0.0, 0.6001, 0.04)
ls = np.array([mse(t) for t in ts])
for t, l in zip(ts, ls):
    col = round(l / ls.max() * 50)
    print(f"theta={t:.2f} |{' ' * col}*{' ' * (50 - col)}| mse={l:7.4f}")

theta=0.00 |                                              *    | mse= 0.7244
theta=0.04 |                                  *                | mse= 0.5428
theta=0.08 |                        *                          | mse= 0.3878
theta=0.12 |                *                                  | mse= 0.2593
theta=0.16 |          *                                        | mse= 0.1575
theta=0.20 |     *                                             | mse= 0.0822
theta=0.24 |  *                                                | mse= 0.0336
theta=0.28 | *                                                 | mse= 0.0116
theta=0.32 | *                                                 | mse= 0.0161
theta=0.36 |   *                                               | mse= 0.0473
theta=0.40 |       *                                           | mse= 0.1050
theta=0.44 |            *                                      | mse= 0.1894
theta=0.48 |                   *                               | mse= 0.3004
theta=0.52 |                            *                      | mse= 0.4379
theta=0.56 |                                      *            | mse= 0.6021
theta=0.60 |                                                  *| mse= 0.7928

Das ist ein Tal, von der Seite gesehen. Es hat einen Boden, die Wände steigen auf beiden Seiten glatt an, und — das ist der Teil, den die Treppe aus Kapitel 1 nicht bieten konnte — an jedem einzelnen Punkt darauf gibt es eine wohldefinierte Richtung für „bergab“. Merk dir diese Form. Kapitel 3 dreht sich vollständig darum, sie hinunterzugehen, ohne alle sechshunderteins Punkte zu besuchen, und darum, was sich ändert, wenn ein Tal mehr als einen Boden hat.

Wir haben ein Tal, weil wir quadriert haben. Absoluter Fehler hätte unten einen Knick erzeugt; größter Fehler hätte flache Strecken erzeugt, auf denen das Verschieben der Linie gar nichts ändert. Quadrieren ist unbestreitbar bequem — und Bequemlichkeit ist ungefähr der Grund, den die meisten Kurse geben, vierfach verkleidet: Es macht Fehler positiv (das tut der Betrag auch); es bestraft große Fehler stärker (warum sollte es?); es ist differenzierbar (die vierte Potenz auch); es ist das, was alle verwenden (ja, und das ist kein Argument).

Hier ist die ehrliche Position. Quadratischer Fehler wählte Linie B, absoluter Fehler wählte Linie A. Eine davon ist für diese Fabrik richtig und die andere falsch, und nichts, was bisher gesagt wurde, kann dir sagen, welche. Um die Regel zu wählen, musst du etwas darüber wissen, wie die Messwerte von der Linie abgewichen sind, und das ist eine Frage über die Welt, nicht über Mathematik. Um sie zu beantworten, braucht es ein kleines Stück Maschinerie.

Hier ist die Behauptung, die aus „welche Linie ist besser“ eine Frage mit Antwort macht.

Nimm an, die Breite eines Teils ist die Linie plus ein zufälliger Fehler, und nimm an, dieser Fehler wird aus einer Gauß-Verteilung — der Glockenkurve — mit Mittelwert null und Standardabweichung $\sigma$ gezogen:

Die Dichte der Gauß-Verteilung ist

Jetzt tu etwas, was das Perzeptron nicht konnte. Für eine gegebene Kandidatensteigung $\theta$ hat jeder Messwert ein Residuum, und die Formel oben verwandelt dieses Residuum in eine Zahl: Wie plausibel ist ein Fehler genau dieser Größe, wenn diese Steigung die Wahrheit ist? Ein Messwert auf der Linie bekommt eine große Zahl, ein Messwert einen halben Millimeter daneben eine kleine.

Die Messwerte sind unabhängig — die Schieblehre erinnert sich nicht an das letzte Teil —, also sagt die Produktregel, dass die Plausibilität des ganzen Notizbuchs das Produkt der einzelnen Dichten ist. Dieses Produkt ist die Likelihood von $\theta$.¹ Beachte die Richtung, denn genau um diese Richtung geht es in Bayes' Regel: Die Daten sind fix und bekannt, und der Parameter variiert. Das ist nicht „die Wahrscheinlichkeit der Steigung“. Es ist die Wahrscheinlichkeit, die das Modell den Daten zuweist, die du tatsächlich bekommen hast, gelesen als Funktion der Steigung.

SIGMA = 0.12

def gaussian(r, sigma):
    return np.exp(-r ** 2 / (2 * sigma ** 2)) / (sigma * np.sqrt(2 * np.pi))

def likelihood(theta):
    return np.prod(gaussian(v - theta * u, SIGMA))          

for t in (0.25, 0.293, 0.35):
    print(f"theta={t}   likelihood = {likelihood(t):.6g}")

theta=0.25   likelihood = 521.952
theta=0.293   likelihood = 2.42028e+07
theta=0.35   likelihood = 0.190312

Eine Steigung von 0,293 macht dieses Notizbuch sechsundvierzigtausendmal plausibler als 0,25 und einhundertsiebenundzwanzig Millionen Mal plausibler als 0,35. Maximum Likelihood ist das Prinzip, den Parameter zu wählen, der das tatsächlich Beobachtete so wenig überraschend wie möglich macht. Es ist kein Theorem, sondern ein Vorschlag dafür, was „best“ bedeuten sollte — ein Vorschlag mit Inhalt, weil er dich zwingt, deine Annahme über das Rauschen zu nennen, bevor du irgendetwas bewerten darfst.

Führe dieselben drei Codezeilen auf einem Monat von Schichten statt auf einer aus, und die Methode fällt auseinander.

rng = np.random.default_rng(7)
u_big = rng.uniform(-5.25, 5.25, 2000)                     # 2000 readings, not 20
v_big = 0.30 * u_big + 0.12 * rng.standard_normal(2000)

print("2000 readings, sigma = 0.12 mm :", np.prod(gaussian(v_big - 0.30 * u_big, 0.12)))
noisy = 0.30 * u_big + 2.0 * rng.standard_normal(2000)
print("2000 readings, sigma = 2.00 mm :", np.prod(gaussian(noisy - 0.30 * u_big, 2.0)))
print("largest float64 :", np.finfo(np.float64).max)

RuntimeWarning: overflow encountered in reduce
2000 readings, sigma = 0.12 mm : inf
2000 readings, sigma = 2.00 mm : 0.0
largest float64 : 1.7976931348623157e+308

Zweitausend Multiplikationen, und die Antwort ist inf. Ändere eine Konstante — eine ungenauere Schieblehre, sodass die Dichten kleiner als 1 statt größer werden —, und derselbe Code gibt 0.0 zurück. Beide Antworten sind falsch, in entgegengesetzte Richtungen, keine wirft eine Exception, die du abfangen kannst, und die zweite gibt nicht einmal eine Warnung aus.

Mit der Mathematik ist nichts falsch. Die Likelihood bei diesen Einstellungen ist eine vollkommen wohldefinierte endliche Zahl: Ihr natürlicher Logarithmus ist 1400,91, also ist die Zahl selbst etwa $10^{608}$. Das Problem ist, dass dein Computer diese Zahl nicht hat, und es lohnt sich, genau zu verstehen, welche Zahlen er hat, denn das ist nicht das letzte Mal, dass er über das Ergebnis entscheidet.

Die Lösung für das explodierende Produkt ist die übliche: Logarithmen nehmen. Der Logarithmus verwandelt Produkte in Summen, er ist streng monoton steigend und kann deshalb den Ort des Maximums nicht verschieben, und eine Summe aus zweitausend moderaten Zahlen verarbeitet float64 ohne Beschwerde. Konventionsgemäß nehmen wir die negative Log-Likelihood, damit besser kleiner bedeutet. Setze nun die Gauß-Dichte ein und sieh, was passiert.

Das ist das Ergebnis, für das dieses Kapitel existiert, und es verdient, ohne Einschränkung formuliert zu werden: Der quadratische Fehler ist keine Konvention. Er ist die negative Log-Likelihood einer Gauß-Verteilung, mit entfernten Konstanten. Quadratischen Fehler zu minimieren ist exakt derselbe Akt wie zu behaupten, dass deine Fehler gaußverteilt sind, und zu fragen, welcher Parameter deine Daten am wenigsten überraschend macht. Du hast diese Behauptung die ganze Zeit aufgestellt; man hat es dir nur nicht gesagt.

Die Äquivalenz ist prüfbar, also prüfe sie: Scanne dieselben sechshunderteins Steigungen mit der vollständigen negativen Log-Likelihood, Konstanten und allem, und mit schlichtem quadratischem Fehler.

N = v.size

def nll(theta):
    r = v - theta * u
    return N * np.log(SIGMA * np.sqrt(2 * np.pi)) + np.sum(r ** 2) / (2 * SIGMA ** 2)

nlls = np.array([nll(t) for t in grid])
mses = np.array([mse(t) for t in grid])
print(f"argmin of the negative log-likelihood : theta={grid[nlls.argmin()]:.3f}  nll={nlls.min():.6f}")
print(f"argmin of the mean squared error      : theta={grid[mses.argmin()]:.3f}  mse={mses.min():.6f}")
print("same index:", nlls.argmin() == mses.argmin())

argmin of the negative log-likelihood : theta=0.293  nll=-17.001977
argmin of the mean squared error      : theta=0.293  mse=0.010115
same index: True

Verschiedene Zahlen auf der vertikalen Achse, und eine davon ist negativ, was eine Summe von Quadraten nie ist: Eine negative Log-Likelihood kann unter null gehen, weil eine Dichte größer als 1 sein kann. Derselbe Boden desselben Tals, bis zum letzten Gitterpunkt.

Wir haben $\sigma$ per Dekret auf 0,12 gesetzt, und niemand im Werk kennt die Streuung des Fehlers der Schieblehre. Behandle sie als zweite Unbekannte und lass Maximum Likelihood auch darüber entscheiden. Hier kommt der konstante Term, den wir gerade weggeworfen haben, zurück, weil er das Einzige ist, was zwischen dem Modell und einer Behauptung perfekter Präzision steht.

r = v - 0.293 * u
sigmas = np.arange(0.01, 1.0001, 0.0001)
nll_sigma = N * np.log(sigmas * np.sqrt(2 * np.pi)) + np.sum(r ** 2) / (2 * sigmas ** 2)

print("best sigma on the grid       :", round(float(sigmas[nll_sigma.argmin()]), 4))
print("sqrt(mean squared residual)  :", round(float(np.sqrt(np.mean(r ** 2))), 4))

best sigma on the grid       : 0.1006
sqrt(mean squared residual)  : 0.1006

Die beiden stimmen auf vier Dezimalstellen überein, und nicht zufällig: Differenzierst du diesen Ausdruck und setzt ihn auf null, erhältst du exakt $\hat{\sigma}^2 = \frac{1}{N}\sum r_i^2$. Der mittlere quadratische Fehler ist also nicht nur wie eine Varianz. Unter diesem Modell ist er die Maximum-Likelihood-Schätzung der Varianz des Rauschens — die Zahl, die du die ganze Zeit minimiert hast, war eine Schätzung dafür, wie verrauscht dein Sensor ist.

Eine Feinheit, billig auszusprechen und teuer später neu zu entdecken: Diese Schätzung ist nach unten verzerrt, weil die Residuen gegen eine Anpassung gemessen wurden, die selbst so gewählt wurde, dass sie klein werden. Simuliere es — zweihunderttausend Notizbücher mit je zwanzig Messwerten, gezogen aus einer Verteilung, deren wahre Varianz exakt 1 ist, wobei der eine Parameter der Anpassung aus den Messwerten selbst geschätzt wird. Die Summe der Quadrate durch $N$ zu teilen ergibt einen Durchschnitt von 0,9501; durch $N-1$ zu teilen ergibt 1,0001; und $(N-1)/N$ ist punktgenau 0,95. Jeder Parameter, den du anpasst, kostet einen Freiheitsgrad, und das ist die kleinste sichtbare Instanz eines viel größeren Problems: Ein Modell sieht auf den Daten, an die es angepasst wurde, immer besser aus. Kapitel 4 macht daraus die Disziplin, Daten zurückzuhalten, und Kapitel 6 gibt dem Effekt seinen Namen.

Wenn quadratischer Fehler behauptet, dass das Rauschen gaußverteilt ist, lautet die nächste Frage, was passiert, wenn die Behauptung falsch ist. Nicht leicht falsch — falsch auf die Art, wie echte Messungen falsch sind.

In der Werkhalle sind die meisten Schieblehrenmessungen bis auf ein Zehntel Millimeter gut, und ein- oder zweimal pro Schicht gerät ein Span unter die Backe, und die Messung liegt um mehrere Millimeter daneben. Solche Fehler sind heavy-tailed: meistens klein, gelegentlich enorm, und enorm viel häufiger, als eine Glockenkurve erlaubt. Die Cauchy-Verteilung ist das klassische saubere Modell für dieses Verhalten, und ihre Dichte ist so einfach wie die der Gauß-Verteilung:

Der Unterschied liegt im Schwanz: Die Gauß-Verteilung fällt wie $e^{-\varepsilon^2}$ ab, brutal schnell, und die Cauchy-Verteilung wie $1/\varepsilon^2$, kaum überhaupt. Die Konsequenz ist leichter zu sehen als zu sagen:

rng = np.random.default_rng(3)
g = 0.12 * rng.standard_normal(10 ** 6)          # Gaussian noise
c = 0.12 * rng.standard_cauchy(10 ** 6)          # Cauchy noise, same scale
for k in (10 ** 2, 10 ** 3, 10 ** 4, 10 ** 5, 10 ** 6):
    print(f"{k:>9,} samples   gaussian var {g[:k].var():.4f}   cauchy var {c[:k].var():10.2f}")

      100 samples   gaussian var 0.0164   cauchy var       0.26
    1,000 samples   gaussian var 0.0146   cauchy var      59.88
   10,000 samples   gaussian var 0.0145   cauchy var     358.17
  100,000 samples   gaussian var 0.0144   cauchy var    3097.98
1,000,000 samples   gaussian var 0.0144   cauchy var   32886.10

Die Stichprobenvarianz der Gauß-Verteilung pendelt sich bei 0,0144 ein, also $0.12^2$, und bleibt dort. Die der Cauchy-Verteilung steigt, und steigt weiter, solange du sampelst, weil es nichts gibt, wohin sie konvergieren könnte: Die Cauchy-Verteilung hat keine Varianz und auch keinen Mittelwert. Quadratischer Fehler, dessen ganzes Geschäft darin besteht, einen Durchschnitt von Quadraten zu minimieren, wird nach einer Größe gefragt, die nicht existiert.

Hier ist also eine Schicht, in der die Schieblehre getäuscht wurde. Dieselben zwanzig Stunden, dieselbe Klinge, derselbe Drift von 0,30 Millimetern pro Stunde — nur ist das Rauschen jetzt Cauchy. Passe zweimal an: einmal durch Minimieren quadrierter Residuen, einmal durch Minimieren der negativen Log-Likelihood des Rauschens, das die Daten tatsächlich erzeugt hat. Der Zentrierungstrick hilft hier nicht — er fixiert den Achsenabschnitt nur für quadratischen Fehler —, also laufen beide Fits per brute force über ein Gitter von Achsenabschnitten und Steigungen, da wir immer noch keine Möglichkeit haben, den Boden eines Tals zu finden, außer es zu besuchen.

SWARF = np.array([
    (0.5, 20.08), (1.0, 21.95), (1.5, 20.86), (2.0, 27.51), (2.5, 20.64),
    (3.0, 20.75), (3.5, 21.01), (4.0, 21.03), (4.5, 21.37), (5.0, 20.60),
    (5.5, 22.03), (6.0, 21.95), (6.5, 21.98), (7.0, 22.01), (7.5, 21.73),
    (8.0, 22.97), (8.5, 22.60), (9.0, 22.66), (9.5, 22.44), (10.0, 22.78),
])
hs, ys = SWARF[:, 0], SWARF[:, 1]

A = np.arange(18.0, 22.001, 0.005)      # 801 intercepts
B = np.arange(-0.20, 0.8001, 0.002)     # 501 slopes
R = ys - (A[:, None, None] + B[None, :, None] * hs)      # every line against every point

SCALE = 0.12
square = np.sum(R ** 2, axis=2)                          # least squares          
cauchy = np.sum(np.log(1 + (R / SCALE) ** 2), axis=2)    # Cauchy likelihood      

for name, surface in (("least squares", square), ("Cauchy likelihood", cauchy)):
    i, j = np.unravel_index(surface.argmin(), surface.shape)
    print(f"{name:>18}:  width = {A[i]:.3f} + {B[j]:.4f} * hours"
          f"   -> 23.5 mm at hour {(23.5 - A[i]) / B[j]:.2f}")
print(f"{'the truth':>18}:  width = 20.000 + 0.3000 * hours"
      f"   -> 23.5 mm at hour {(23.5 - 20.0) / 0.30:.2f}")
print(f"{A.size * B.size:,} candidate lines evaluated")

Die beiden hervorgehobenen Zeilen sind der gesamte Unterschied zwischen den Fits. Nimm den Logarithmus der Cauchy-Dichte, streiche die Konstanten genau wie zuvor, und $\sum \log\left(1 + (r/s)^2\right)$ bleibt übrig. Dasselbe Rezept, eine andere Behauptung über das Rauschen.

     least squares:  width = 21.380 + 0.1080 * hours   -> 23.5 mm at hour 19.63
 Cauchy likelihood:  width = 19.935 + 0.3020 * hours   -> 23.5 mm at hour 11.80
         the truth:  width = 20.000 + 0.3000 * hours   -> 23.5 mm at hour 11.67
401,301 candidate lines evaluated

Least Squares meldet einen Drift von 0,108 Millimetern pro Stunde, ungefähr ein Drittel der echten Rate, und schließt, dass die Klinge bis Stunde 19,6 gut ist. Die wahre Antwort ist Stunde 11,7. Handelt das Werk nach diesem Fit, läuft die Presse acht zusätzliche Stunden und produziert Teile außerhalb der Toleranz, mit Autorität der standardmäßigsten Verlustfunktion des Feldes. Der Cauchy-Fit, mit denselben zwanzig Messwerten, demselben Gitter und einer einzigen Codezeile Unterschied, landet bei Stunde 11,8.

Zwei Einwände verdienen Antworten, weil beide das Erste sind, was ein guter Engineer sagt.

Der Ausreißer ist offensichtlich — lösch ihn einfach. Kannst du tun, und es hilft, und es reicht nicht. Das Löschen des einzelnen schlimmsten Messwerts verschiebt die Least-Squares-Steigung von 0,108 auf 0,239, was den Klingenwechsel immer noch auf Stunde 13,1 setzt, anderthalb Stunden zu spät; den schlimmsten zu löschen, neu zu fitten und dann zu löschen, was jetzt am schlimmsten ist, bringt dich auf 0,286 — und beachte, dass das bereits ein Verfahren ist, keine Beobachtung: Lösche stattdessen die zwei größten Residuen des ursprünglichen Fits, und du landest bei 0,223. Aber du hast jetzt Ermessensentscheidungen getroffen, die du nicht aufschreiben oder verteidigen kannst, und die Regel zu automatisieren rettet sie nicht: Größtes-Residuum-löschen-dann-neu-fitten, über tausend simulierte Schichten ausgeführt, hat einen medianen Steigungsfehler von 0,0177 gegenüber 0,0100 beim Likelihood-Fit und liegt in 14,7 % der Schichten um mehr als 0,05 daneben, gegenüber 1,3 %. Löschen ist ein Flicken auf einer falschen Annahme. Die Likelihood braucht keinen Flicken, weil sie nie angenommen hat, der Ausreißer sei unmöglich.

Du hast einen glücklichen Datensatz ausgewählt. Dieser Einwand ist genau richtig, und deshalb simuliert das letzte Experiment tausend unabhängige Schichten und fittet jede auf beide Arten neu.

A = np.arange(18.0, 22.001, 0.02)        # a coarser grid: a thousand fits to do
B = np.arange(-0.20, 0.8001, 0.005)
lines = A[:, None, None] + B[None, :, None] * hs
rng = np.random.default_rng(2026)
err_sq, err_ca = [], []

for _ in range(1000):                                        # 1000 independent shifts
    ys = 20.00 + 0.30 * hs + SCALE * rng.standard_cauchy(hs.size)
    R = ys - lines
    _, j = np.unravel_index(np.sum(R ** 2, axis=2).argmin(), (A.size, B.size))
    _, q = np.unravel_index(np.sum(np.log1p((R / SCALE) ** 2), axis=2).argmin(), (A.size, B.size))
    err_sq.append(abs(B[j] - 0.30))
    err_ca.append(abs(B[q] - 0.30))

err_sq, err_ca = np.array(err_sq), np.array(err_ca)
for name, e in (("least squares", err_sq), ("Cauchy likelihood", err_ca)):
    print(f"{name:>18}: median slope error {np.median(e):.4f} mm/h"
          f"   off by more than 0.05 in {100 * np.mean(e > 0.05):4.1f}% of shifts"
          f"   worst {e.max():.3f}")
print(f"the likelihood fit is the closer of the two in {100 * np.mean(err_ca < err_sq):.1f}% of shifts")

     least squares: median slope error 0.0350 mm/h   off by more than 0.05 in 40.4% of shifts   worst 0.500
 Cauchy likelihood: median slope error 0.0100 mm/h   off by more than 0.05 in  1.3% of shifts   worst 0.090
the likelihood fit is the closer of the two in 75.6% of shifts

Median, nicht Mittelwert, aus demselben Grund wie alles andere in diesem Abschnitt: Die Least-Squares-Fehler werden von einer Cauchy-Verteilung getrieben, also ist ihr Durchschnitt keine stabile Größe, die man berichten sollte. Least Squares liegt in zwei von fünf Schichten grob falsch; der Likelihood-Fit liegt in einer von siebenundsiebzig Schichten grob falsch, und sein schlimmster Fehlschlag über tausend Schichten ist kleiner als ein Fünftel des schlimmsten von Least Squares.

Nichts davon macht quadratischen Fehler schlecht. Es macht ihn spezifisch, und die Arithmetik sagt genau warum. Nimm ein Residuum von 0,1 mm und eines von 7 mm. Quadriert trägt die schlechte Messung 4.900-mal so viel zur Summe bei wie die gute, also wird die Linie körperlich zu ihr hingezogen; unter der Cauchy-Log-Likelihood tragen dieselben beiden Residuen 0,527 und 8,133 bei, ein Verhältnis von 15,4. Die schlechte Messung zählt immer noch, sie darf nur nicht entscheiden. Das ist der Anfang robuster Statistik, wo Hubers Loss von 1964 den Unterschied teilt, indem er sich für kleine Residuen quadratisch und für große linear verhält,⁷ und wo Tukey bereits gezeigt hatte, wie wenig Kontamination nötig ist, damit die Stichprobenvarianz ein schlechteres Werkzeug ist als die mittlere absolute Abweichung.⁸

Eine historische Notiz, zu gut, um sie wegzulassen. Least Squares wurde zuerst veröffentlicht, von Legendre im Jahr 1805, als bequemes algebraisches Hilfsmittel ohne Rechtfertigung über die Tatsache hinaus, dass es funktionierte.⁹ Vier Jahre später drehte Gauss das Argument um: Er nahm als gegeben, dass das arithmetische Mittel die richtige Art ist, wiederholte Messungen zu kombinieren, fragte, welche Fehlerverteilung den Mittelwert zum wahrscheinlichsten Wert macht, und zeigte, dass im Wesentlichen nur eine das tut — die, die heute nach ihm benannt ist.¹⁰ Die Herleitung in diesem Kapitel ist seine, sie ist mehr als zwei Jahrhunderte alt, und sie ist immer noch der Teil, den die meisten Kurse weglassen.

Schau aber noch einmal auf das Tal. Wenn du bei $\theta = 0.20$ mit einem Verlust von 0,0822 stehst, ist die Richtung von „bergab“ kein Geheimnis — du kannst sie auf der Seite sehen, die Kurve fällt nach rechts. Wenn du die Verlustfunktion fragen könntest, in welche Richtung sie an dem Punkt, an dem du stehst, geneigt ist, ohne sie irgendwo anders auszuwerten, könntest du einen Schritt in diese Richtung machen, wieder fragen und wiederholen, bis der Boden flach ist.

Diese Frage hat einen Namen. Die Steigung einer Funktion an einem Punkt ist ihre Ableitung, und bei einer Funktion vieler Parameter ist die Sammlung der Steigungen in jede Richtung zugleich der gradient. Kapitel 1 konnte keinen verwenden, weil der Fehler des Perzeptrons eine Treppe ohne Steigung war, nach der man hätte fragen können. Dieses Kapitel hat etwas Besseres gebaut: einen Loss, der überall glatt ist und aus einer benannten Annahme stammt, nicht aus einer Vorliebe.

Die Frage für Kapitel 3 ist also nicht mehr, ob eine Steigung existiert. Sie lautet, wie man sie berechnet, warum eine Bewegung gegen sie bergab statt bergauf führt — ein Vorzeichen, das fast jeder Kurs dich auf Glauben akzeptieren lässt — und wie weit man geht, bevor man wieder fragt. Das stellt sich als die eine Zahl heraus, die entscheidet, ob ein Trainingslauf konvergiert, für immer um die Antwort oszilliert oder ins Unendliche davonläuft.

Ebenfalls lesenswert neben diesem Kapitel: Prince, Understanding Deep Learning §5.1–5.2 und Anhang C, der jeden Loss im Buch aus Maximum Likelihood in der hier verwendeten Reihenfolge aufbaut; Goodfellow, Bengio und Courville, Deep Learning §3.1–3.11 und §5.5, dessen Maximum-Likelihood-Abschnitt auch die KL-Divergenz herleitet, die Kapitel 4 braucht; Murphy, Probabilistic Machine Learning: An Introduction Kapitel 2 und §4.2, dazu, was Maximum Likelihood garantiert und was nicht; Deisenroth, Faisal und Ong, Mathematics for Machine Learning §6.1–6.4 für Summenregel, Produktregel und Bayes' Regel sauber durchgeführt; Tom Mitchells kurze CMU-Notiz Estimating Probabilities: MLE and MAP (2016); und §22.7 von Dive into Deep Learning, das dasselbe Ergebnis in ausführbarem Code erreicht.