Από πού προέρχεται μια συνάρτηση απώλειας: πιθανότητα, όχι σύμβαση

Εδώ είναι οι μετρήσεις και οι τρεις γραμμές, βαθμολογημένες με τρεις τρόπους: squared error, στο οποίο καταφεύγουν όλοι· absolute error, στο οποίο μπορεί να κατέφευγε ένας στατιστικός· και worst error, στο οποίο θα κατέφευγε ο μηχανουργός, επειδή ο ελεγκτής δεν νοιάζεται για τον μέσο όρο σας — απορρίπτει το ένα εξάρτημα που είναι εκτός ανοχής.

Το NumPy εμφανίζεται εδώ, ένα κεφάλαιο μετά το καθαρό-Python perceptron, για έναν λόγο: μέχρι το τέλος αυτού του κεφαλαίου αξιολογούμε τετρακόσιες χιλιάδες υποψήφιες γραμμές απέναντι σε είκοσι μετρήσεις η καθεμία, και ένα Python loop είναι το λάθος εργαλείο για αυτό. Είναι επίσης η σημειογραφία στην οποία είναι γραμμένη κάθε πηγή που παρατίθεται παρακάτω.

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

Η ποσότητα στις επισημασμένες γραμμές είναι το υπόλοιπο: αυτό που είπε η γραμμή μείον αυτό που είπε το παχύμετρο, ένας αριθμός ανά μέτρηση. Κάθε κανόνας βαθμολόγησης σε αυτό το κεφάλαιο, και κάθε loss function στα είκοσι οκτώ κεφάλαια μετά από αυτό, είναι ένας τρόπος να συμπιεστεί μια λίστα υπολοίπων σε έναν μόνο αριθμό. Διαφέρουν μόνο στο πώς τη συμπιέζουν.

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

Διαβάστε τις στήλες, όχι τις σειρές. Το squared error λέει B, το absolute error λέει A, το worst error λέει C: τρεις κανόνες, τρεις νικητές, στα ίδια είκοσι σημεία.

Επέλεξα αυτές τις τρεις γραμμές ώστε να διαφωνούν, και πρέπει να το πω καθαρά. Το ζήτημα είναι πόσο εύκολο ήταν — λίγα λεπτά αναζήτησης σε λογικές τεταγμένες τομής και κλίσεις βρίσκουν εκατοντάδες τέτοιες τριάδες. Η κατάταξη είναι ιδιότητα του κανόνα που επιλέξατε, όχι γεγονός για τις γραμμές, άρα ο κανόνας δεν είναι λεπτομέρεια υλοποίησης: είναι ο ορισμός του προβλήματος. Κάτι που θέτει το ερώτημα για το οποίο υπάρχει αυτό το κεφάλαιο: με ποια κριτήρια τον επιλέγετε;

Πρώτα ένα μικρότερο ζήτημα, επειδή δεν υπάρχουν τρεις γραμμές αλλά άπειρες. Πάρτε προς το παρόν το squared error, αφού αυτό παίρνουν όλοι, και μικρύνετε το πρόβλημα σε έναν μόνο αριθμό χρησιμοποιώντας το τέχνασμα που γλίτωσε το perceptron από έντεκα χιλιάδες epochs στο Κεφάλαιο 1: αφαιρέστε τον μέσο όρο και από τις δύο στήλες. Μόλις το νέφος σημείων κεντραριστεί στην αρχή των αξόνων, η καλύτερη γραμμή με βάση το squared error περνά ακριβώς από την αρχή — άρα η τεταγμένη τομής έχει τακτοποιηθεί και μένει να επιλεγεί μόνο η κλίση.

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

Εξακόσιες μία υποψήφιες κλίσεις, ένας νικητής: 0,293 χιλιοστά ανά ώρα απέναντι σε αλήθεια 0,300. Είκοσι θορυβώδεις μετρήσεις και ένα for-loop έφτασαν μέσα σε ένα εκατοστό του χιλιοστού ανά ώρα — δύο και ένα τρίτο τοις εκατό.

Το ενδιαφέρον δεν είναι ο νικητής αλλά το σχήμα της αναζήτησης. Εκτυπώστε ολόκληρη την καμπύλη, περιστραμμένη ώστε το loss να τρέχει από αριστερά προς τα δεξιά:

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

Αυτή είναι μια κοιλάδα, όπως φαίνεται από το πλάι. Έχει έναν πυθμένα, τα τοιχώματα ανεβαίνουν ομαλά και στις δύο πλευρές, και — αυτό είναι το κομμάτι που η σκάλα του Κεφαλαίου 1 δεν μπορούσε να προσφέρει — σε κάθε σημείο της υπάρχει μια καλά ορισμένη κατεύθυνση προς το «κατηφορικά». Θυμηθείτε αυτό το σχήμα. Το Κεφάλαιο 3 αφορά εξ ολοκλήρου το πώς κατεβαίνουμε σε αυτή χωρίς να επισκεφθούμε και τα εξακόσια ένα σημεία, και τι αλλάζει όταν μια κοιλάδα έχει περισσότερους από έναν πυθμένες.

Έχουμε κοιλάδα επειδή υψώσαμε στο τετράγωνο. Το absolute error θα της είχε δώσει μια γωνία στον πυθμένα· το worst error θα της είχε δώσει επίπεδα τμήματα όπου η μετακίνηση της γραμμής δεν αλλάζει απολύτως τίποτα. Το τετράγωνο είναι αναμφισβήτητα βολικό — και η ευκολία είναι περίπου ο λόγος που δίνουν τα περισσότερα μαθήματα, ντυμένος με τέσσερις τρόπους: κάνει τα σφάλματα θετικά (το ίδιο κάνει και η απόλυτη τιμή)· τιμωρεί περισσότερο τα μεγάλα σφάλματα (γιατί θα έπρεπε;)· είναι παραγωγίσιμο (το ίδιο και η τέταρτη δύναμη)· είναι αυτό που χρησιμοποιούν όλοι (είναι, και αυτό δεν είναι επιχείρημα).

Η ειλικρινής θέση είναι η εξής. Το squared error επέλεξε τη γραμμή B και το absolute error επέλεξε τη γραμμή A. Μία από αυτές είναι σωστή για αυτό το εργοστάσιο και η άλλη λάθος, και τίποτα από όσα ειπώθηκαν μέχρι τώρα δεν μπορεί να σας πει ποια. Για να επιλέξετε τον κανόνα πρέπει να ξέρετε κάτι για το πώς οι μετρήσεις κατέληξαν να διαφέρουν από τη γραμμή, και αυτό είναι ερώτημα για τον κόσμο, όχι για τα μαθηματικά. Η απάντησή του χρειάζεται ένα μικρό κομμάτι μηχανισμού.

Εδώ είναι ο ισχυρισμός που μετατρέπει το «ποια γραμμή είναι καλύτερη» σε ερώτηση με απάντηση.

Υποθέστε ότι το πλάτος ενός εξαρτήματος είναι η γραμμή συν ένα τυχαίο σφάλμα, και υποθέστε ότι αυτό το σφάλμα αντλείται από μια Gaussian — την καμπύλη καμπάνας — με μέσο όρο μηδέν και τυπική απόκλιση $\sigma$:

Η πυκνότητα της Gaussian είναι

Τώρα κάντε κάτι που το perceptron δεν μπορούσε. Για μια δεδομένη υποψήφια κλίση $\theta$, κάθε μέτρηση έχει ένα υπόλοιπο, και ο παραπάνω τύπος μετατρέπει αυτό το υπόλοιπο σε έναν αριθμό: πόσο εύλογο είναι ένα σφάλμα ακριβώς αυτού του μεγέθους, αν αυτή η κλίση είναι η αλήθεια; Μια μέτρηση πάνω στη γραμμή παίρνει μεγάλο αριθμό, μια μέτρηση μισό χιλιοστό μακριά παίρνει μικρό.

Οι μετρήσεις είναι ανεξάρτητες — το παχύμετρο δεν θυμάται το προηγούμενο εξάρτημα — οπότε ο κανόνας γινομένου λέει ότι η ευλογοφάνεια ολόκληρου του σημειωματάριου είναι το γινόμενο των επιμέρους πυκνοτήτων. Αυτό το γινόμενο είναι η likelihood του $\theta$.¹ Προσέξτε την κατεύθυνση, επειδή αυτή είναι η κατεύθυνση για την οποία μιλά ο κανόνας του Bayes: τα δεδομένα είναι σταθερά και γνωστά, και η παράμετρος είναι αυτή που μεταβάλλεται. Αυτό δεν είναι «η πιθανότητα της κλίσης». Είναι η πιθανότητα που το μοντέλο αναθέτει στα δεδομένα που όντως πήρατε, ιδωμένη ως συνάρτηση της κλίσης.

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

Μια κλίση 0,293 κάνει αυτό το σημειωματάριο σαράντα έξι χιλιάδες φορές πιο εύλογο από το 0,25, και εκατόν είκοσι επτά εκατομμύρια φορές πιο εύλογο από το 0,35. Maximum likelihood είναι η αρχή ότι επιλέγετε την παράμετρο που κάνει αυτό που πραγματικά παρατηρήσατε όσο το δυνατόν λιγότερο απροσδόκητο. Δεν είναι θεώρημα αλλά πρόταση για το τι θα έπρεπε να σημαίνει «καλύτερο» — μια πρόταση με περιεχόμενο, επειδή σας αναγκάζει να δηλώσετε την υπόθεσή σας για τον θόρυβο πριν σας επιτραπεί να βαθμολογήσετε οτιδήποτε.

Τρέξτε τις ίδιες τρεις γραμμές κώδικα σε έναν μήνα βαρδιών αντί για μία, και η μέθοδος καταρρέει.

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

Δύο χιλιάδες πολλαπλασιασμοί και η απάντηση είναι inf. Αλλάξτε μία σταθερά — ένα πιο αδέξιο παχύμετρο, ώστε οι πυκνότητες να βγαίνουν μικρότερες από 1 αντί για μεγαλύτερες — και ο ίδιος κώδικας επιστρέφει 0.0. Και οι δύο απαντήσεις είναι λάθος, σε αντίθετες κατευθύνσεις, καμία δεν σηκώνει exception που μπορείτε να πιάσετε, και η δεύτερη δεν τυπώνει καν προειδοποίηση.

Δεν υπάρχει τίποτα λάθος στα μαθηματικά. Η likelihood σε αυτές τις ρυθμίσεις είναι ένας απολύτως καλά ορισμένος πεπερασμένος αριθμός: ο φυσικός λογάριθμός της είναι 1400,91, άρα ο ίδιος ο αριθμός είναι περίπου $10^{608}$. Το πρόβλημα είναι ότι ο υπολογιστής σας δεν έχει αυτόν τον αριθμό, και αξίζει να καταλάβετε ακριβώς ποιους αριθμούς έχει, επειδή δεν είναι η τελευταία φορά που θα αποφασίσει το αποτέλεσμα.

Η διόρθωση για το γινόμενο που εκρήγνυται είναι η συνηθισμένη: παίρνουμε λογαρίθμους. Ο λογάριθμος μετατρέπει γινόμενα σε αθροίσματα, είναι αυστηρά αύξων άρα δεν μπορεί να μετακινήσει τη θέση του μέγιστου, και ένα άθροισμα δύο χιλιάδων μέτριων αριθμών είναι κάτι που το float64 χειρίζεται χωρίς παράπονο. Κατά σύμβαση παίρνουμε την αρνητική log-likelihood, ώστε καλύτερο να σημαίνει μικρότερο. Τώρα αντικαταστήστε την πυκνότητα της Gaussian και δείτε τι συμβαίνει.

Αυτό είναι το αποτέλεσμα για το οποίο υπάρχει το κεφάλαιο, και αξίζει να διατυπωθεί χωρίς επιφυλάξεις: το squared error δεν είναι σύμβαση. Είναι η αρνητική log-likelihood μιας Gaussian, με τις σταθερές αφαιρεμένες. Η ελαχιστοποίηση του squared error είναι ακριβώς η ίδια πράξη με το να ισχυρίζεστε ότι τα σφάλματά σας είναι Gaussian και να ρωτάτε ποια παράμετρος κάνει τα δεδομένα σας λιγότερο απροσδόκητα. Κάνατε αυτόν τον ισχυρισμό από την αρχή· απλώς δεν σας το έλεγαν.

Η ισοδυναμία είναι ελέγξιμη, οπότε ελέγξτε τη: σαρώστε τις ίδιες εξακόσιες μία κλίσεις με την πλήρη αρνητική log-likelihood, σταθερές και όλα, και με απλό squared error.

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

Διαφορετικοί αριθμοί στον κατακόρυφο άξονα, και ο ένας είναι αρνητικός, κάτι που ένα άθροισμα τετραγώνων δεν είναι ποτέ: μια αρνητική log-likelihood μπορεί να πέσει κάτω από το μηδέν, επειδή μια πυκνότητα μπορεί να ξεπεράσει το 1. Ίδιος πυθμένας της ίδιας κοιλάδας, μέχρι το τελευταίο σημείο του πλέγματος.

Καθορίσαμε το $\sigma$ στο 0,12 με αυθαίρετη απόφαση, και κανείς στο εργοστάσιο δεν ξέρει τη διασπορά του σφάλματος του παχύμετρου. Αντιμετωπίστε το ως δεύτερο άγνωστο και αφήστε τη maximum likelihood να το αποφασίσει κι αυτό. Εδώ ο σταθερός όρος που μόλις απορρίψαμε επιστρέφει, επειδή είναι το μόνο πράγμα που στέκεται ανάμεσα στο μοντέλο και σε έναν ισχυρισμό τέλειας ακρίβειας.

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

Τα δύο συμφωνούν σε τέσσερα δεκαδικά ψηφία, και όχι τυχαία: αν παραγωγίσετε αυτή την έκφραση και τη θέσετε ίση με μηδέν, παίρνετε ακριβώς $\hat{\sigma}^2 = \frac{1}{N}\sum r_i^2$. Άρα το mean squared error δεν είναι απλώς σαν διακύμανση. Κάτω από αυτό το μοντέλο είναι η maximum-likelihood εκτίμηση της διακύμανσης του θορύβου — ο αριθμός που ελαχιστοποιούσατε από την αρχή ήταν εκτίμηση του πόσο θορυβώδης είναι ο αισθητήρας σας.

Μια λεπτομέρεια, φτηνή στη διατύπωση και ακριβή αν την ξανανακαλύψετε αργότερα: αυτή η εκτίμηση είναι μεροληπτικά χαμηλή, επειδή τα υπόλοιπα μετρήθηκαν απέναντι σε μια προσαρμογή που είχε επιλεγεί η ίδια για να τα κάνει μικρά. Προσομοιώστε το — διακόσιες χιλιάδες σημειωματάρια των είκοσι μετρήσεων το καθένα, αντλημένα από μια κατανομή της οποίας η πραγματική διακύμανση είναι ακριβώς 1, με τη μία παράμετρο της προσαρμογής να εκτιμάται από τις ίδιες τις μετρήσεις. Η διαίρεση του αθροίσματος τετραγώνων με $N$ δίνει μέσο όρο 0,9501· η διαίρεση με $N-1$ δίνει 1,0001· και το $(N-1)/N$ είναι ακριβώς 0,95. Κάθε παράμετρος που προσαρμόζετε κοστίζει έναν βαθμό ελευθερίας, και αυτή είναι η μικρότερη ορατή περίπτωση ενός πολύ μεγαλύτερου προβλήματος: ένα μοντέλο φαίνεται πάντα καλύτερο στα δεδομένα στα οποία προσαρμόστηκε. Το Κεφάλαιο 4 το μετατρέπει αυτό στην πειθαρχία του να κρατάμε δεδομένα στην άκρη, και το Κεφάλαιο 6 δίνει στο φαινόμενο το όνομά του.

Αν το squared error ισχυρίζεται ότι ο θόρυβος είναι Gaussian, η επόμενη ερώτηση είναι τι συμβαίνει όταν ο ισχυρισμός είναι ψευδής. Όχι ελαφρώς ψευδής — ψευδής με τον τρόπο που είναι ψευδείς οι πραγματικές μετρήσεις.

Στο δάπεδο του εργοστασίου, οι περισσότερες μετρήσεις του παχύμετρου είναι καλές στο ένα δέκατο του χιλιοστού, και μία ή δύο φορές ανά βάρδια ένα ρινίσμα μπαίνει κάτω από τη σιαγόνα και η μέτρηση φεύγει κατά αρκετά χιλιοστά. Τέτοια σφάλματα είναι heavy-tailed: μικρά τον περισσότερο χρόνο, περιστασιακά τεράστια, και τεράστια πολύ συχνότερα απ’ όσο επιτρέπει μια καμπύλη καμπάνας. Η κατανομή Cauchy είναι το κλασικό καθαρό μοντέλο αυτής της συμπεριφοράς, και η πυκνότητά της είναι τόσο απλή όσο της Gaussian:

Η διαφορά είναι η ουρά: η Gaussian πέφτει σαν $e^{-\varepsilon^2}$, βίαια γρήγορα, και η Cauchy σαν $1/\varepsilon^2$, σχεδόν καθόλου. Η συνέπεια φαίνεται ευκολότερα απ’ όσο λέγεται:

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

Η δειγματική διακύμανση της Gaussian σταθεροποιείται στο 0,0144, που είναι $0.12^2$, και μένει εκεί. Της Cauchy ανεβαίνει, και συνεχίζει να ανεβαίνει όσο συνεχίζετε να δειγματοληπτείτε, επειδή δεν υπάρχει τίποτα στο οποίο να συγκλίνει: η κατανομή Cauchy δεν έχει διακύμανση, ούτε μέσο όρο. Το squared error, του οποίου όλη η δουλειά είναι να ελαχιστοποιεί έναν μέσο όρο τετραγώνων, καλείται να βρει μια ποσότητα που δεν υπάρχει.

Ορίστε λοιπόν μία βάρδια όπου το παχύμετρο ξεγελάστηκε. Ίδιες είκοσι ώρες, ίδια λεπίδα, ίδια μετατόπιση 0,30 χιλιοστά ανά ώρα — μόνο που ο θόρυβος είναι τώρα Cauchy. Προσαρμόστε το δύο φορές: μία ελαχιστοποιώντας τα squared residuals, μία ελαχιστοποιώντας την αρνητική log-likelihood του θορύβου που πραγματικά δημιούργησε τα δεδομένα. Το τέχνασμα του κεντραρίσματος δεν βοηθά εδώ — καθηλώνει την τεταγμένη τομής μόνο για το squared error — άρα και οι δύο προσαρμογές γίνονται με brute force πάνω σε πλέγμα τεταγμένων τομής και κλίσεων, αφού ακόμη δεν έχουμε τρόπο να βρούμε τον πυθμένα μιας κοιλάδας παρά μόνο επισκεπτόμενοι τον.

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

Οι δύο επισημασμένες γραμμές είναι ολόκληρη η διαφορά ανάμεσα στις προσαρμογές. Πάρτε τον λογάριθμο της πυκνότητας Cauchy, πετάξτε τις σταθερές ακριβώς όπως πριν, και το $\sum \log\left(1 + (r/s)^2\right)$ είναι αυτό που επιβιώνει. Ίδια συνταγή, διαφορετικός ισχυρισμός για τον θόρυβο.

     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 αναφέρουν μετατόπιση 0,108 χιλιοστά την ώρα, περίπου το ένα τρίτο του πραγματικού ρυθμού, και συμπεραίνουν ότι η λεπίδα είναι εντάξει μέχρι την ώρα 19,6. Η πραγματική απάντηση είναι η ώρα 11,7. Αν ενεργήσει με βάση αυτή την προσαρμογή, το εργοστάσιο λειτουργεί την πρέσα για οκτώ επιπλέον ώρες φτιάχνοντας εξαρτήματα εκτός ανοχής, με την εξουσία της πιο καθιερωμένης loss function του πεδίου. Η προσαρμογή Cauchy, χρησιμοποιώντας τις ίδιες είκοσι μετρήσεις, το ίδιο πλέγμα και διαφορά μίας γραμμής στον κώδικα, καταλήγει στην ώρα 11,8.

Δύο ενστάσεις αξίζουν απάντηση, επειδή είναι και οι δύο το πρώτο πράγμα που λέει ένας καλός μηχανικός.

Το outlier είναι προφανές — απλώς διαγράψτε το. Μπορείτε, και βοηθά, και δεν αρκεί. Η διαγραφή της μίας χειρότερης μέτρησης μετακινεί την κλίση least-squares από 0,108 σε 0,239, που εξακολουθεί να βάζει την αλλαγή λεπίδας στην ώρα 13,1, μιάμιση ώρα αργά· η διαγραφή της χειρότερης, η επαναπροσαρμογή και η διαγραφή όποιας είναι χειρότερη τώρα σας φέρνει στο 0,286 — και προσέξτε ότι αυτό είναι ήδη διαδικασία, όχι παρατήρηση: διαγράψτε αντί γι’ αυτό τα δύο μεγαλύτερα υπόλοιπα της αρχικής προσαρμογής και καταλήγετε στο 0,223. Αλλά τώρα έχετε κάνει κρίσεις που δεν μπορείτε να γράψετε ή να υπερασπιστείτε, και η αυτοματοποίηση του κανόνα δεν τον σώζει: το drop-the-largest-residual-then-refit, τρεγμένο σε χίλιες προσομοιωμένες βάρδιες, έχει διάμεσο σφάλμα κλίσης 0,0177 απέναντι στο 0,0100 της likelihood προσαρμογής, και αποκλίνει πάνω από 0,05 στο 14,7% των βαρδιών απέναντι στο 1,3%. Η διαγραφή είναι μπάλωμα πάνω σε λάθος υπόθεση. Η likelihood δεν χρειάζεται μπάλωμα, επειδή ποτέ δεν υπέθεσε ότι το outlier ήταν αδύνατο.

Διαλέξατε τυχερό dataset. Αυτή η ένσταση είναι ακριβώς σωστή, γι’ αυτό και το τελευταίο πείραμα προσομοιώνει χίλιες ανεξάρτητες βάρδιες και επαναπροσαρμόζει και με τους δύο τρόπους σε καθεμία.

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

Διάμεσος, όχι μέσος όρος, για τον ίδιο λόγο με όλα τα άλλα σε αυτή την ενότητα: τα σφάλματα least-squares οδηγούνται από μια Cauchy, άρα ο μέσος όρος τους δεν είναι σταθερό πράγμα για να αναφέρετε. Τα least squares είναι πολύ λάθος σε δύο βάρδιες στις πέντε· η likelihood προσαρμογή είναι πολύ λάθος σε μία βάρδια στις εβδομήντα επτά, και η χειρότερη αποτυχία της σε χίλιες βάρδιες είναι μικρότερη από το ένα πέμπτο της χειρότερης των least squares.

Τίποτα από αυτά δεν κάνει το squared error κακό. Το κάνει συγκεκριμένο, και η αριθμητική λέει ακριβώς γιατί. Πάρτε ένα υπόλοιπο 0,1 mm και ένα 7 mm. Στο τετράγωνο, η κακή μέτρηση συνεισφέρει 4.900 φορές όσο η καλή στο σύνολο, άρα η γραμμή σέρνεται ολόκληρη προς αυτήν· κάτω από τη Cauchy log-likelihood τα ίδια δύο υπόλοιπα συνεισφέρουν 0,527 και 8,133, λόγος 15,4. Η κακή μέτρηση εξακολουθεί να μετρά, απλώς δεν αποφασίζει. Αυτή είναι η αρχή της robust statistics, όπου το loss του Huber του 1964 χωρίζει τη διαφορά συμπεριφερόμενο τετραγωνικά για μικρά υπόλοιπα και γραμμικά για μεγάλα,⁷ και όπου ο Tukey είχε ήδη δείξει πόσο λίγη μόλυνση χρειάζεται για να γίνει η δειγματική διακύμανση χειρότερο εργαλείο από τη μέση απόλυτη απόκλιση.⁸

Και μια ιστορική σημείωση, πολύ καλή για να μείνει έξω. Τα least squares δημοσιεύτηκαν πρώτα, από τον Legendre το 1805, ως βολικό αλγεβρικό εργαλείο χωρίς δικαιολόγηση πέρα από το ότι δούλευε.⁹ Τέσσερα χρόνια αργότερα ο Gauss έτρεξε το επιχείρημα ανάποδα: πήρε ως δεδομένο ότι ο αριθμητικός μέσος είναι ο σωστός τρόπος να συνδυάζονται επαναλαμβανόμενες μετρήσεις, ρώτησε ποια κατανομή σφάλματος κάνει τον μέσο την πιο πιθανή τιμή, και έδειξε ότι ουσιαστικά μόνο μία το κάνει — αυτή που σήμερα φέρει το όνομά του.¹⁰ Η παραγωγή σε αυτό το κεφάλαιο είναι δική του, είναι πάνω από δύο αιώνες παλιά, και εξακολουθεί να είναι το κομμάτι που τα περισσότερα μαθήματα αφήνουν έξω.

Κοιτάξτε όμως ξανά την κοιλάδα. Στεκόμενοι στο $\theta = 0.20$ με loss 0,0822, η κατεύθυνση του «κατηφορικά» δεν είναι μυστήριο — τη βλέπετε στη σελίδα, η καμπύλη κατεβαίνει προς τα δεξιά. Αν μπορούσατε να ρωτήσετε τη loss function προς ποια κατεύθυνση γέρνει στο σημείο όπου στέκεστε, χωρίς να την αξιολογήσετε πουθενά αλλού, θα μπορούσατε να κάνετε ένα βήμα προς τα εκεί, να ρωτήσετε ξανά και να επαναλάβετε μέχρι το έδαφος να γίνει επίπεδο.

Αυτή η ερώτηση έχει όνομα. Η κλίση μιας συνάρτησης σε ένα σημείο είναι η παράγωγός της, και για μια συνάρτηση πολλών παραμέτρων η συλλογή των κλίσεων προς κάθε κατεύθυνση ταυτόχρονα είναι το gradient. Το Κεφάλαιο 1 δεν μπορούσε να χρησιμοποιήσει ένα, επειδή το σφάλμα του perceptron ήταν σκάλα χωρίς κλίση για να ρωτήσει. Αυτό το κεφάλαιο έχτισε κάτι καλύτερο: ένα loss που είναι ομαλό παντού και που προήλθε από δηλωμένη υπόθεση αντί για προτίμηση.

Άρα το ερώτημα για το Κεφάλαιο 3 δεν είναι πλέον αν υπάρχει κλίση. Είναι πώς να την υπολογίσετε, γιατί η κίνηση ενάντια σε αυτή κατεβαίνει αντί να ανεβαίνει — ένα πρόσημο που σχεδόν κάθε μάθημα σας ζητά να δεχτείτε με πίστη — και πόσο μεγάλο βήμα να κάνετε πριν ρωτήσετε ξανά, κάτι που αποδεικνύεται ότι είναι ο ένας αριθμός που αποφασίζει αν ένα training run συγκλίνει, ταλαντώνεται για πάντα γύρω από την απάντηση ή φεύγει στο άπειρο.

Αξίζει επίσης να διαβάσετε παράλληλα με αυτό το κεφάλαιο: Prince, Understanding Deep Learning §5.1–5.2 και Appendix C, που χτίζει κάθε loss του βιβλίου από maximum likelihood με τη σειρά που χρησιμοποιείται εδώ· Goodfellow, Bengio and Courville, Deep Learning §3.1–3.11 και §5.5, των οποίων η ενότητα maximum-likelihood παράγει επίσης την KL divergence που χρειάζεται το Κεφάλαιο 4· Murphy, Probabilistic Machine Learning: An Introduction chapter 2 και §4.2, για το τι εγγυάται και τι δεν εγγυάται η maximum likelihood· Deisenroth, Faisal and Ong, Mathematics for Machine Learning §6.1–6.4 για τον κανόνα αθροίσματος, τον κανόνα γινομένου και τον κανόνα του Bayes σωστά διατυπωμένους· το σύντομο σημείωμα CMU του Tom Mitchell Estimating Probabilities: MLE and MAP (2016)· και §22.7 του Dive into Deep Learning, που φτάνει στο ίδιο αποτέλεσμα με εκτελέσιμο κώδικα.