מאיפה מגיעה פונקציית loss: likelihood, לא מוסכמה

הנה הקריאות ושלושת הקווים, מנוקדים בשלוש דרכים: שגיאה ריבועית, שאליה כולם ניגשים; שגיאה מוחלטת, שסטטיסטיקאי אולי היה בוחר; והשגיאה הגרועה ביותר, שהמכונאי היה בוחר, כי המפקח לא מתעניין בממוצע שלכם — הוא פוסל את החלק היחיד שחורג מהסבילות.

NumPy נכנסת כאן, פרק אחד אחרי ה-perceptron ב-Python טהור, מסיבה אחת: עד סוף הפרק הזה נעריך ארבע מאות אלף קווים מועמדים מול עשרים קריאות כל אחד, ולולאת Python היא הכלי הלא נכון לזה. זו גם הנוטציה שבה כתוב כל מקור שמצוטט למטה.

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

הכמות בשורות המודגשות היא ה-residual: מה שהקו אמר פחות מה שהקליבר אמר, מספר אחד לכל קריאה. כל כלל ניקוד בפרק הזה, וכל פונקציית loss בעשרים ושמונה הפרקים שאחריו, היא דרך כלשהי לדחוס רשימת residuals למספר יחיד. הם נבדלים רק באופן הדחיסה.

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

קראו את העמודות, לא את השורות. שגיאה ריבועית אומרת B, שגיאה מוחלטת אומרת A, השגיאה הגרועה ביותר אומרת C: שלושה כללים, שלושה מנצחים, על אותן עשרים נקודות.

בחרתי את שלושת הקווים האלה כך שהם יחלקו זה על זה, וראוי לומר זאת במפורש. הנקודה היא כמה קל היה לעשות זאת — כמה דקות של חיפוש על פני חותכים ושיפועים שנראים סבירים מגלות מאות שלשות כאלה. הדירוג הוא תכונה של הכלל שבחרתם, לא עובדה על הקווים, ולכן הכלל אינו פרט מימוש: הוא ההגדרה של הבעיה. וזה מעלה את השאלה שלשמה הפרק הזה קיים: על סמך מה בוחרים אותו?

קודם עניין קטן יותר, כי אין שלושה קווים אלא אינסוף. קחו בינתיים שגיאה ריבועית, כי זה מה שכולם לוקחים, וצמצמו את הבעיה למספר יחיד בעזרת הטריק שחסך ל-perceptron אחד-עשר אלף epochs בפרק 1: מחסרים את הממוצע משתי העמודות. ברגע שענן הנקודות ממורכז סביב הראשית, הקו הטוב ביותר תחת שגיאה ריבועית עובר בדיוק דרך הראשית — כך שהחותך נקבע ונשאר לבחור רק את השיפוע.

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 הגיעו למרחק של מאית מילימטר לשעה — שניים ושליש אחוזים.

החלק המעניין אינו המנצח אלא צורת החיפוש. הדפיסו את כל העקומה, מסובבת כך שה-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 עוסק כולו בהליכה במורד שלה בלי לבקר בכל שש מאות ואחת הנקודות, ובמה שמשתנה כשלעמק יש יותר מתחתית אחת.

יש לנו עמק כי העלינו בריבוע. שגיאה מוחלטת הייתה נותנת לו שפיץ בתחתית; השגיאה הגרועה ביותר הייתה נותנת לו מקטעים שטוחים שבהם הזזת הקו לא משנה דבר. העלאה בריבוע נוחה בלי ספק — ונוחות היא בערך הסיבה שרוב הקורסים נותנים, בארבע עטיפות שונות: היא הופכת שגיאות לחיוביות (גם ערך מוחלט עושה זאת); היא מענישה שגיאות גדולות יותר (למה בעצם?); היא גזירה (גם חזקה רביעית); זה מה שכולם משתמשים בו (נכון, וזה לא נימוק).

זו העמדה הכנה. שגיאה ריבועית בחרה את קו B ושגיאה מוחלטת בחרה את קו A. אחד מהם נכון למפעל הזה והשני שגוי, ושום דבר שנאמר עד כה לא יכול לומר לכם איזה. כדי לבחור את הכלל צריך לדעת משהו על איך הקריאות הגיעו להיות שונות מהקו, וזו שאלה על העולם, לא על מתמטיקה. כדי לענות עליה צריך מכונה קטנה אחת.

הנה הטענה שהופכת את ״איזה קו טוב יותר״ לשאלה שיש לה תשובה.

נניח שרוחב של חלק הוא הקו בתוספת שגיאה אקראית, ונניח שהשגיאה הזאת נמשכת מגאוסיאן — עקומת הפעמון — עם ממוצע אפס וסטיית תקן $\sigma$:

הצפיפות של הגאוסיאן היא

עכשיו עשו משהו שה-perceptron לא יכול היה לעשות. עבור שיפוע מועמד נתון $\theta$, לכל קריאה יש residual, והנוסחה למעלה הופכת את ה-residual הזה למספר: עד כמה סבירה שגיאה בדיוק בגודל הזה, אם השיפוע הזה הוא האמת? קריאה על הקו מקבלת מספר גדול, קריאה במרחק חצי מילימטר מקבלת מספר קטן.

הקריאות עצמאיות — הקליבר לא זוכר את החלק הקודם — ולכן כלל המכפלה אומר שהסבירות של כל המחברת היא מכפלת הצפיפויות הבודדות. המכפלה הזאת היא ה-likelihood של $\theta$.¹ שימו לב לכיוון, כי זה הכיוון שעליו מדבר כלל בייס: הנתונים קבועים וידועים, והפרמטר הוא זה שמשתנה. זו אינה ״ההסתברות של השיפוע״. זו ההסתברות שהמודל מקצה לנתונים שקיבלתם בפועל, כשהיא נקראת כפונקציה של השיפוע.

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 השלילי, כך שטוב יותר פירושו קטן יותר. עכשיו הציבו את הצפיפות הגאוסיאנית וצפו במה שקורה.

זו התוצאה שלשמה הפרק קיים, וראוי לומר אותה בלי הסתייגות: שגיאה ריבועית אינה מוסכמה. היא ה-negative log-likelihood של גאוסיאן, לאחר הסרת הקבועים. מזעור שגיאה ריבועית הוא בדיוק אותו מעשה כמו לטעון שהשגיאות שלכם גאוסיאניות ולשאול איזה פרמטר הופך את הנתונים שלכם לפחות מפתיעים. טענתם זאת כל הזמן; פשוט לא אמרו לכם.

אפשר לבדוק את השקילות, אז בדקו אותה: סרקו את אותם שש מאות ואחד שיפועים עם ה-negative log-likelihood המלא, קבועים והכול, ועם שגיאה ריבועית רגילה.

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

מספרים שונים על הציר האנכי, ואחד מהם שלילי, דבר שסכום ריבועים לעולם אינו: negative 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 של שונות הרעש — המספר שמזערתם כל הזמן היה אומדן לכמה רועש החיישן שלכם.

קמט אחד, זול לניסוח ויקר לגילוי מחדש מאוחר יותר: האומדן הזה מוטה כלפי מטה, כי ה-residuals נמדדו מול התאמה שנבחרה בעצמה כדי להקטין אותם. סימולציה מראה זאת — מאתיים אלף מחברות של עשרים קריאות כל אחת, נמשכות מהתפלגות שהשונות האמיתית שלה בדיוק 1, עם הפרמטר היחיד של ההתאמה שנאמד מתוך הקריאות עצמן. חלוקת סכום הריבועים ב-$N$ נותנת ממוצע של 0.9501; חלוקה ב-$N-1$ נותנת 1.0001; ו-$(N-1)/N$ הוא 0.95 בדיוק. כל פרמטר שמותאם עולה דרגת חופש אחת, וזה המקרה הגלוי הקטן ביותר של בעיה גדולה בהרבה: מודל תמיד נראה טוב יותר על הנתונים שעליהם הותאם. פרק 4 הופך זאת למשמעת של שמירת נתונים בצד, ופרק 6 נותן לאפקט את שמו.

אם שגיאה ריבועית טוענת שהרעש גאוסיאני, השאלה הבאה היא מה קורה כשהטענה שגויה. לא שגויה מעט — שגויה באופן שבו מדידות אמיתיות שגויות.

ברצפת הייצור, רוב קריאות הקליבר טובות עד עשירית מילימטר, ופעם או פעמיים במשמרת שבב מתכת קטן נכנס מתחת ללסת והקריאה סוטה בכמה מילימטרים. שגיאות כאלה הן בעלות זנבות כבדים: קטנות רוב הזמן, עצומות מדי פעם, ועצומות הרבה יותר לעיתים קרובות מכפי שעקומת פעמון מאפשרת. התפלגות קושי היא המודל הנקי הסטנדרטי להתנהגות הזאת, והצפיפות שלה פשוטה כמו זו של הגאוסיאן:

ההבדל הוא הזנב: הגאוסיאן דועך כמו $e^{-\varepsilon^2}$, מהר באכזריות, וקושי כמו $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

שונות המדגם של הגאוסיאן מתייצבת על 0.0144, שהוא $0.12^2$, ונשארת שם. זו של קושי מטפסת, וממשיכה לטפס כל עוד דוגמים, כי אין לה לאן להתכנס: להתפלגות קושי אין שונות, וגם אין ממוצע. שגיאה ריבועית, שכל עניינה הוא מזעור ממוצע של ריבועים, מתבקשת למצוא כמות שאינה קיימת.

אז הנה משמרת אחת שבה הקליבר הוטעה. אותן עשרים שעות, אותו להב, אותו סחף של 0.30 מילימטר לשעה — רק שהרעש עכשיו הוא קושי. התאימו אותה פעמיים: פעם אחת על ידי מזעור residuals ריבועיים, פעם אחת על ידי מזעור ה-negative log-likelihood של הרעש שבאמת יצר את הנתונים. טריק המרכז אינו עוזר כאן — הוא מקבע את החותך רק עבור שגיאה ריבועית — ולכן שתי ההתאמות נעשות בכוח גס על פני סריג של חותכים ושיפועים, כי עדיין אין לנו דרך למצוא את תחתית העמק מלבד לבקר בה.

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

שתי השורות המודגשות הן כל ההבדל בין ההתאמות. קחו את ה-log של צפיפות קושי, זרקו את הקבועים בדיוק כמו קודם, ו-$\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 הסטנדרטית ביותר בתחום. התאמת קושי, עם אותן עשרים קריאות, אותו סריג, והבדל של שורה אחת בקוד, נוחתת על שעה 11.8.

שתי התנגדויות ראויות לתשובות, כי שתיהן הדבר הראשון שמהנדס טוב אומר.

ה-outlier ברור — פשוט תמחקו אותו. אפשר, וזה עוזר, וזה לא מספיק. מחיקת הקריאה הגרועה היחידה מזיזה את שיפוע least-squares מ-0.108 ל-0.239, מה שעדיין שם את החלפת הלהב בשעה 13.1, באיחור של שעה וחצי; מחיקת הגרועה ביותר, התאמה מחדש, ומחיקת מה שהכי גרוע עכשיו מביאה אתכם ל-0.286 — ושימו לב שזה כבר נוהל, לא תצפית: מחקו במקום זאת את שני ה-residuals הגדולים ביותר של ההתאמה המקורית ותנחתו על 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 מונעות על ידי קושי, ולכן הממוצע שלהן אינו דבר יציב לדיווח. least squares שוגה מאוד בשתי משמרות מתוך חמש; התאמת ה-likelihood שוגה מאוד במשמרת אחת מתוך שבעים ושבע, והכישלון הגרוע ביותר שלה על פני אלף משמרות קטן מפחות מחמישית מהכישלון הגרוע ביותר של least squares.

כל זה לא הופך שגיאה ריבועית לרעה. זה הופך אותה לספציפית, והאריתמטיקה אומרת בדיוק למה. קחו residual של 0.1 מ״מ ואחד של 7 מ״מ. בריבוע, הקריאה הרעה תורמת לסך הכול פי 4,900 מהטובה, ולכן הקו נגרר בגופו לעברה; תחת log-likelihood של קושי אותם שני residuals תורמים 0.527 ו-8.133, יחס של 15.4. הקריאה הרעה עדיין נחשבת, היא פשוט לא זוכה להחליט. זו תחילתה של סטטיסטיקה רובסטית, שבה פונקציית loss של Huber מ-1964 מפצלת את ההבדל על ידי התנהגות ריבועית עבור residuals קטנים וליניארית עבור גדולים,⁷ ושבה Tukey כבר הראה כמה מעט זיהום צריך כדי להפוך את שונות המדגם לכלי גרוע יותר מהסטייה המוחלטת הממוצעת.⁸

גם הערה היסטורית אחת, טובה מכדי להשמיט. least squares פורסמה ראשונה, על ידי Legendre ב-1805, כמכשיר אלגברי נוח ללא הצדקה מעבר לכך שהיא עבדה.⁹ ארבע שנים מאוחר יותר Gauss הריץ את הטיעון לאחור: הוא קיבל כנתון שהממוצע האריתמטי הוא הדרך הנכונה לשלב מדידות חוזרות, שאל איזו התפלגות שגיאה הופכת את הממוצע לערך הסביר ביותר, והראה שבעצם רק אחת עושה זאת — זו שנקראת כיום על שמו.¹⁰ הגזירה בפרק הזה היא שלו, היא בת יותר ממאתיים שנה, והיא עדיין החלק שרוב הקורסים משאירים בחוץ.

אבל הביטו שוב בעמק. כשעומדים ב-$\theta = 0.20$ עם loss של 0.0822, כיוון ה״ירידה״ אינו מסתורי — אפשר לראות אותו על הדף, העקומה יורדת ימינה. אם הייתם יכולים לשאול את פונקציית loss לאיזה כיוון היא נוטה בנקודה שבה אתם עומדים, בלי להעריך אותה בשום מקום אחר, הייתם יכולים לצעוד בכיוון הזה, לשאול שוב, ולחזור עד שהקרקע שטוחה.

לשאלה הזאת יש שם. השיפוע של פונקציה בנקודה הוא הנגזרת שלה, ועבור פונקציה של פרמטרים רבים אוסף השיפועים בכל הכיוונים בבת אחת הוא הגרדיאנט. פרק 1 לא יכול היה להשתמש באחד, כי השגיאה של ה-perceptron הייתה מדרגות בלי שיפוע לשאול עליו. הפרק הזה בנה משהו טוב יותר: loss חלקה בכל מקום, שמגיעה מהנחה מוצהרת ולא מהעדפה.

לכן השאלה לפרק 3 כבר אינה אם שיפוע קיים. היא איך לחשב אותו, למה תנועה נגדו יורדת למטה ולא עולה למעלה — סימן שכמעט כל קורס מבקש מכם לקבל באמונה — וכמה רחוק לצעוד לפני ששואלים שוב, שמתברר כמספר האחד שמחליט אם ריצת אימון מתכנסת, מתנדנדת סביב התשובה לנצח, או בורחת לאינסוף.

כדאי לקרוא לצד הפרק הזה גם את Prince, Understanding Deep Learning §5.1–5.2 ונספח 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 פרק 2 ו-§4.2, על מה maximum likelihood כן ולא מבטיח; Deisenroth, Faisal and Ong, Mathematics for Machine Learning §6.1–6.4 עבור כלל הסכום, כלל המכפלה וכלל בייס בצורה מסודרת; ההערה הקצרה של Tom Mitchell מ-CMU, Estimating Probabilities: MLE and MAP (2016); ו-§22.7 של Dive into Deep Learning, שמגיע לאותה תוצאה בקוד שניתן להרצה.