Skąd bierze się funkcja straty: z likelihood, nie z konwencji

Oto odczyty i trzy linie, ocenione na trzy sposoby: błędem kwadratowym, po który sięga każdy; błędem bezwzględnym, po który mógłby sięgnąć statystyk; oraz najgorszym błędem, po który sięgnąłby operator maszyny, bo inspektora nie obchodzi twoja średnia — odrzuca pojedynczą część poza tolerancją.

NumPy pojawia się tutaj, jeden rozdział po perceptronie napisanym w czystym Pythonie, z jednego powodu: pod koniec tego rozdziału oceniamy czterysta tysięcy kandydujących linii względem dwudziestu odczytów każda, a pętla w Pythonie jest do tego złym narzędziem. To także notacja, w której zapisane jest każde cytowane niżej źródło.

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}")               

Wielkość w wyróżnionych liniach to residuum: to, co powiedziała linia, minus to, co powiedziała suwmiarka, jedna liczba na odczyt. Każda reguła oceny w tym rozdziale i każda funkcja straty w kolejnych dwudziestu ośmiu rozdziałach jest jakimś sposobem ściśnięcia listy residuów do jednej liczby. Różnią się tylko tym, jak ściskają.

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

Czytaj kolumny, nie wiersze. Błąd kwadratowy mówi B, błąd bezwzględny mówi A, najgorszy błąd mówi C: trzy reguły, trzej zwycięzcy, na tych samych dwudziestu punktach.

Wybrałem te trzy linie tak, żeby się nie zgadzały, i powinienem powiedzieć to wprost. Chodzi o to, jak łatwe to było — kilka minut przeszukiwania rozsądnie wyglądających wyrazów wolnych i nachyleń daje setki takich trójek. Ranking jest własnością reguły, którą wybrałeś, nie faktem o liniach, więc reguła nie jest szczegółem implementacji: ona jest definicją problemu. A to prowadzi do pytania, na które istnieje ten rozdział: na jakiej podstawie ją wybierasz?

Najpierw mniejsza sprawa, bo nie ma trzech linii, tylko nieskończenie wiele. Weźmy na razie błąd kwadratowy, skoro wszyscy go biorą, i zmniejszmy problem do jednej liczby przy użyciu sztuczki, która oszczędziła perceptronowi jedenaście tysięcy epok w rozdziale 1: odejmij średnią od obu kolumn. Gdy chmura punktów jest wycentrowana wokół początku układu, najlepsza linia według błędu kwadratowego przechodzi dokładnie przez początek — więc wyraz wolny jest ustalony i zostaje do wyboru tylko nachylenie.

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

Sześćset jeden kandydujących nachyleń, jeden zwycięzca: 0,293 milimetra na godzinę wobec prawdy 0,300. Dwadzieścia zaszumionych odczytów i pętla for zbliżyły się na setną milimetra na godzinę — dwa i jedną trzecią procenta.

Interesujący nie jest zwycięzca, lecz kształt przeszukiwania. Wypisz całą krzywą, obróconą tak, by strata biegła od lewej do prawej:

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

To dolina widziana z boku. Ma jedno dno, ściany po obu stronach wznoszą się gładko i — to część, której schody z rozdziału 1 nie mogły zaoferować — w każdym pojedynczym punkcie istnieje dobrze określony kierunek „w dół”. Zapamiętaj ten kształt. Rozdział 3 w całości dotyczy schodzenia po nim bez odwiedzania wszystkich sześciuset jeden punktów oraz tego, co się zmienia, gdy dolina ma więcej niż jedno dno.

Mamy dolinę, bo podnieśliśmy do kwadratu. Błąd bezwzględny dałby załamanie na dnie; najgorszy błąd dałby płaskie odcinki, gdzie przesunięcie linii niczego nie zmienia. Podnoszenie do kwadratu jest bezsprzecznie wygodne — i wygoda jest mniej więcej powodem, który podaje większość kursów, ubranym na cztery sposoby: sprawia, że błędy są dodatnie (wartość bezwzględna też); mocniej karze duże błędy (a dlaczego miałaby?); jest różniczkowalne (czwarta potęga też); wszyscy tego używają (używają, ale to nie argument).

Uczciwe stanowisko jest takie. Błąd kwadratowy wybrał linię B, a błąd bezwzględny wybrał linię A. Jedna z nich jest właściwa dla tej fabryki, a druga błędna, i nic z tego, co dotąd powiedziano, nie mówi ci która. Aby wybrać regułę, musisz wiedzieć coś o tym, jak odczyty zaczęły różnić się od linii, a to pytanie o świat, nie o matematykę. Odpowiedź wymaga jednego małego mechanizmu.

Oto twierdzenie, które zamienia „która linia jest lepsza” w pytanie z odpowiedzią.

Załóżmy, że szerokość części to linia plus błąd losowy, i załóżmy, że ten błąd jest losowany z rozkładu Gaussa — krzywej dzwonowej — o średniej zero i odchyleniu standardowym $\sigma$:

Gęstość rozkładu Gaussa to

Teraz zrób coś, czego perceptron nie potrafił. Dla danego kandydującego nachylenia $\theta$ każdy odczyt ma residuum, a powyższy wzór zamienia to residuum w liczbę: jak wiarygodny jest błąd dokładnie tej wielkości, jeśli to nachylenie jest prawdą? Odczyt na linii dostaje dużą liczbę, odczyt oddalony o pół milimetra — małą.

Odczyty są niezależne — suwmiarka nie pamięta poprzedniej części — więc reguła iloczynu mówi, że wiarygodność całego notesu jest iloczynem pojedynczych gęstości. Ten iloczyn to likelihood $\theta$.¹ Zwróć uwagę na kierunek, bo o nim mówi reguła Bayesa: dane są stałe i znane, a zmienia się parametr. To nie jest „prawdopodobieństwo nachylenia”. To prawdopodobieństwo, które model przypisuje danym, które faktycznie dostałeś, odczytane jako funkcja nachylenia.

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

Nachylenie 0,293 czyni ten notes czterdzieści sześć tysięcy razy bardziej wiarygodnym niż 0,25 i sto dwadzieścia siedem milionów razy bardziej wiarygodnym niż 0,35. Maximum likelihood to zasada, według której wybierasz parametr czyniący to, co faktycznie zaobserwowałeś, możliwie najmniej zaskakującym. To nie twierdzenie, lecz propozycja, co „najlepsze” powinno znaczyć — propozycja z treścią, bo zmusza cię do podania założenia o szumie, zanim wolno ci cokolwiek ocenić.

Uruchom te same trzy linie kodu na miesiącu zmian zamiast na jednej, a metoda się wykłada.

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

Dwa tysiące mnożeń i odpowiedź to inf. Zmień jedną stałą — mniej dokładną suwmiarkę, przez co gęstości wychodzą mniejsze niż 1 zamiast większe — a ten sam kod zwraca 0.0. Obie odpowiedzi są błędne, w przeciwnych kierunkach, żadna nie zgłasza wyjątku, który możesz złapać, a druga nawet nie wypisuje ostrzeżenia.

Z matematyką nie jest nic nie tak. Likelihood przy tych ustawieniach jest doskonale określoną skończoną liczbą: jej logarytm naturalny wynosi 1400,91, więc sama liczba to około $10^{608}$. Problem w tym, że twój komputer nie ma tej liczby, i warto dokładnie zrozumieć, które liczby ma, bo to nie ostatni raz, kiedy zdecyduje o wyniku.

Naprawa eksplodującego iloczynu jest standardowa: weź logarytmy. Logarytm zamienia iloczyny w sumy, jest ściśle rosnący, więc nie może przesunąć położenia maksimum, a suma dwóch tysięcy umiarkowanych liczb to coś, z czym float64 radzi sobie bez narzekania. Z konwencji bierzemy ujemny logarytm likelihood, tak aby lepsze znaczyło mniejsze. Teraz podstaw gęstość Gaussa i zobacz, co się stanie.

To wynik, dla którego istnieje ten rozdział, i zasługuje na stwierdzenie bez asekuracji: błąd kwadratowy nie jest konwencją. Jest ujemnym logarytmem likelihood rozkładu Gaussa, po usunięciu stałych. Minimalizowanie błędu kwadratowego jest dokładnie tym samym aktem co stwierdzenie, że twoje błędy są gaussowskie, i zapytanie, który parametr czyni twoje dane najmniej zaskakującymi. Cały czas składałeś to stwierdzenie; po prostu nikt ci o tym nie mówił.

Równoważność da się sprawdzić, więc ją sprawdź: przeskanuj te same sześćset jeden nachyleń pełnym ujemnym logarytmem likelihood, ze stałymi i wszystkim, oraz zwykłym błędem kwadratowym.

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

Inne liczby na osi pionowej, a jedna z nich jest ujemna, czym suma kwadratów nigdy nie jest: ujemny logarytm likelihood może spaść poniżej zera, bo gęstość może przekraczać 1. To samo dno tej samej doliny, do ostatniego punktu siatki.

Ustaliliśmy $\sigma$ na 0,12 arbitralnie, a nikt w zakładzie nie zna rozrzutu błędu swojej suwmiarki. Potraktuj go jako drugą niewiadomą i pozwól, by maximum likelihood zdecydowało również o nim. Tutaj składnik stały, który właśnie odrzuciliśmy, wraca, bo jest jedyną rzeczą stojącą między modelem a twierdzeniem o idealnej precyzji.

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

Oba wyniki zgadzają się do czterech miejsc po przecinku i nie przez przypadek: zróżniczkowanie tego wyrażenia i przyrównanie do zera daje dokładnie $\hat{\sigma}^2 = \frac{1}{N}\sum r_i^2$. Zatem średni błąd kwadratowy nie jest jedynie podobny do wariancji. W tym modelu jest estymatą maximum likelihood wariancji szumu — liczba, którą przez cały czas minimalizowałeś, była estymatą tego, jak zaszumiony jest twój czujnik.

Jedna komplikacja, tania do wypowiedzenia i droga do ponownego odkrycia później: ta estymata jest obciążona w dół, bo residua mierzono względem dopasowania, które samo zostało wybrane po to, by były małe. Zasymuluj to — dwieście tysięcy notesów po dwadzieścia odczytów każdy, wylosowanych z rozkładu, którego prawdziwa wariancja wynosi dokładnie 1, z jednym parametrem dopasowania estymowanym z samych odczytów. Dzielenie sumy kwadratów przez $N$ daje średnio 0,9501; dzielenie przez $N-1$ daje 1,0001; a $(N-1)/N$ to równo 0,95. Każdy parametr, który dopasowujesz, kosztuje jeden stopień swobody, a to najmniejszy widoczny przypadek znacznie większego problemu: model zawsze wygląda lepiej na danych, do których został dopasowany. Rozdział 4 zamienia to w dyscyplinę odkładania danych na bok, a rozdział 6 nadaje temu efektowi nazwę.

Jeśli błąd kwadratowy stwierdza, że szum jest gaussowski, następne pytanie brzmi: co się dzieje, gdy to stwierdzenie jest fałszywe. Nie trochę fałszywe — fałszywe tak, jak fałszywe bywają prawdziwe pomiary.

Na hali produkcyjnej większość odczytów suwmiarki jest dobra do jednej dziesiątej milimetra, a raz czy dwa na zmianę pod szczękę dostaje się wiór i odczyt myli się o kilka milimetrów. Takie błędy są ciężkoogonowe: zwykle małe, czasem ogromne, i ogromne znacznie częściej, niż pozwala na to krzywa dzwonowa. Rozkład Cauchy’ego jest standardowym czystym modelem takiego zachowania, a jego gęstość jest tak prosta jak gaussowska:

Różnica tkwi w ogonie: rozkład Gaussa spada jak $e^{-\varepsilon^2}$, brutalnie szybko, a Cauchy’ego jak $1/\varepsilon^2$, prawie wcale. Konsekwencję łatwiej zobaczyć niż opisać:

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

Wariancja próby z rozkładu Gaussa stabilizuje się na 0,0144, czyli $0.12^2$, i tam zostaje. Dla rozkładu Cauchy’ego rośnie i rośnie tak długo, jak próbkujesz, bo nie ma niczego, do czego mogłaby zbiec: rozkład Cauchy’ego nie ma wariancji ani nawet średniej. Błąd kwadratowy, którego całym zajęciem jest minimalizowanie średniej kwadratów, jest proszony o wielkość, która nie istnieje.

Oto więc jedna zmiana, podczas której suwmiarka dała się oszukać. Te same dwadzieścia godzin, to samo ostrze, ten sam dryf 0,30 milimetra na godzinę — tylko szum jest teraz Cauchy’ego. Dopasuj go dwa razy: raz minimalizując kwadraty residuów, raz minimalizując ujemny logarytm likelihood szumu, który faktycznie wygenerował dane. Sztuczka z centrowaniem tutaj nie pomaga — unieruchamia wyraz wolny tylko dla błędu kwadratowego — więc oba dopasowania idą brute force po siatce wyrazów wolnych i nachyleń, skoro wciąż nie mamy sposobu, by znaleźć dno doliny inaczej niż przez odwiedzenie go.

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")

Dwie wyróżnione linie to cała różnica między dopasowaniami. Weź logarytm gęstości Cauchy’ego, odrzuć stałe dokładnie tak jak wcześniej, a $\sum \log\left(1 + (r/s)^2\right)$ jest tym, co zostaje. Ten sam przepis, inne twierdzenie o szumie.

     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

Metoda najmniejszych kwadratów zgłasza dryf 0,108 milimetra na godzinę, mniej więcej jedną trzecią rzeczywistego tempa, i wnioskuje, że ostrze jest dobre do godziny 19,6. Prawdziwa odpowiedź to godzina 11,7. Działając według tego dopasowania, zakład prowadzi prasę przez osiem dodatkowych godzin, produkując części poza tolerancją, z autorytetu najbardziej standardowej funkcji straty w tej dziedzinie. Dopasowanie Cauchy’ego, używając tych samych dwudziestu odczytów, tej samej siatki i różnicy jednej linii w kodzie, trafia w godzinę 11,8.

Dwa zastrzeżenia zasługują na odpowiedź, bo oba są pierwszą rzeczą, którą powie dobry inżynier.

Odstający punkt jest oczywisty — po prostu go usuń. Możesz, pomaga, i to za mało. Usunięcie pojedynczego najgorszego odczytu przesuwa nachylenie najmniejszych kwadratów z 0,108 do 0,239, co nadal ustawia wymianę ostrza na godzinę 13,1, półtorej godziny za późno; usunięcie najgorszego, ponowne dopasowanie i usunięcie tego, co jest najgorsze teraz, prowadzi do 0,286 — i zauważ, że to już procedura, nie obserwacja: usuń dwa największe residua pierwotnego dopasowania, a wylądujesz na 0,223. Podjąłeś jednak decyzje uznaniowe, których nie możesz zapisać ani obronić, a automatyzacja reguły jej nie ratuje: procedura „usuń największe residuum, potem dopasuj ponownie”, uruchomiona na tysiącu symulowanych zmian, ma medianę błędu nachylenia 0,0177 wobec 0,0100 dla dopasowania likelihood i myli się o więcej niż 0,05 w 14,7% zmian wobec 1,3%. Usuwanie jest łatą na błędnym założeniu. Likelihood nie potrzebuje łaty, bo nigdy nie zakładało, że punkt odstający jest niemożliwy.

Wybrałeś szczęśliwy zbiór danych. To zastrzeżenie jest dokładnie trafne, dlatego ostatni eksperyment symuluje tysiąc niezależnych zmian i dopasowuje oba sposoby na każdej.

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

Mediana, nie średnia, z tego samego powodu co wszystko inne w tej sekcji: błędy najmniejszych kwadratów są napędzane przez rozkład Cauchy’ego, więc ich średnia nie jest stabilną rzeczą do raportowania. Najmniejsze kwadraty są bardzo błędne w dwóch zmianach na pięć; dopasowanie likelihood jest bardzo błędne w jednej zmianie na siedemdziesiąt siedem, a jego najgorsza porażka w tysiącu zmian jest mniejsza niż jedna piąta najgorszej porażki najmniejszych kwadratów.

Nic z tego nie czyni błędu kwadratowego złym. Czyni go konkretnym, a arytmetyka dokładnie mówi dlaczego. Weź residuum 0,1 mm i residuum 7 mm. Po podniesieniu do kwadratu zły odczyt dokłada do sumy 4 900 razy więcej niż dobry, więc linia zostaje fizycznie przeciągnięta w jego stronę; przy log-likelihood Cauchy’ego te same dwa residua wnoszą 0,527 i 8,133, stosunek 15,4. Zły odczyt nadal się liczy, po prostu nie może decydować. To początek statystyki odpornej, gdzie funkcja straty Hubera z 1964 roku dzieli różnicę, zachowując się kwadratowo dla małych residuów i liniowo dla dużych,⁷ oraz gdzie Tukey już wcześniej pokazał, jak niewiele zanieczyszczenia wystarcza, by wariancja z próby stała się gorszym narzędziem niż średnie odchylenie bezwzględne.⁸

Także jedna uwaga historyczna, zbyt dobra, by ją pominąć. Metoda najmniejszych kwadratów została opublikowana jako pierwsza, przez Legendre’a w 1805 roku, jako wygodne narzędzie algebraiczne bez uzasadnienia poza tym, że działało.⁹ Cztery lata później Gauss przeprowadził argument w drugą stronę: przyjął, że średnia arytmetyczna jest właściwym sposobem łączenia powtarzanych pomiarów, zapytał, który rozkład błędu czyni średnią najbardziej prawdopodobną wartością, i pokazał, że zasadniczo robi to tylko jeden — ten nazwany dziś jego nazwiskiem.¹⁰ Wyprowadzenie w tym rozdziale jest jego, ma ponad dwa stulecia i nadal jest częścią, którą większość kursów pomija.

Spójrz jednak z powrotem na dolinę. Stojąc przy $\theta = 0.20$ ze stratą 0,0822, kierunek „w dół” nie jest tajemnicą — widać go na stronie, krzywa opada w prawo. Gdybyś mógł zapytać funkcję straty, w którą stronę nachyla się w punkcie, w którym stoisz, bez obliczania jej nigdzie indziej, mógłbyś zrobić krok w tę stronę, zapytać ponownie i powtarzać, aż grunt będzie płaski.

To pytanie ma nazwę. Nachylenie funkcji w punkcie to jej pochodna, a dla funkcji wielu parametrów zbiór nachyleń we wszystkich kierunkach naraz to gradient. Rozdział 1 nie mógł z niego skorzystać, bo błąd perceptronu był schodami bez nachylenia, o które można pytać. Ten rozdział zbudował coś lepszego: loss, która jest gładka wszędzie i pochodzi z podanego założenia, nie z preferencji.

Pytanie na rozdział 3 nie brzmi więc już, czy nachylenie istnieje. Brzmi: jak je obliczyć, dlaczego ruch przeciw niemu prowadzi w dół, a nie w górę — znak, który niemal każdy kurs każe przyjąć na wiarę — oraz jak daleki krok zrobić przed kolejnym pytaniem, co okazuje się jedną liczbą decydującą, czy trening zbiegnie, będzie wiecznie oscylował wokół odpowiedzi, czy ucieknie do nieskończoności.

Warto też czytać równolegle z tym rozdziałem: Prince, Understanding Deep Learning §5.1–5.2 i dodatek C, gdzie każda loss w książce budowana jest z maximum likelihood w kolejności użytej tutaj; Goodfellow, Bengio i Courville, Deep Learning §3.1–3.11 i §5.5, których sekcja o maximum likelihood wyprowadza także dywergencję KL potrzebną w rozdziale 4; Murphy, Probabilistic Machine Learning: An Introduction rozdział 2 i §4.2, o tym, co maximum likelihood gwarantuje, a czego nie; Deisenroth, Faisal i Ong, Mathematics for Machine Learning §6.1–6.4 dla porządnego omówienia reguły sumy, reguły iloczynu i reguły Bayesa; krótka notatka Toma Mitchella z CMU Estimating Probabilities: MLE and MAP (2016); oraz §22.7 z Dive into Deep Learning, który dochodzi do tego samego wyniku w uruchamialnym kodzie.