손실 함수는 어디서 오는가: 관습이 아니라 likelihood

다음은 측정값과 세 직선이다. 세 가지 방식으로 채점했다. 모두가 먼저 집어 드는 squared error, 통계학자라면 떠올릴 absolute error, 그리고 기계공이라면 쓸 worst error다. 검사관은 평균에 관심이 없기 때문이다. 허용오차를 벗어난 단 하나의 부품을 불합격 처리할 뿐이다.

NumPy는 pure-Python perceptron 다음 장인 여기에서 등장한다. 이유는 하나다. 이 장이 끝날 무렵 우리는 후보 직선 40만 개를 각각 20개 측정값에 대해 평가하게 되며, 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}")               

강조된 줄의 양은 residual이다. 직선이 말한 값에서 캘리퍼가 말한 값을 뺀 것, 측정값 하나당 숫자 하나다. 이 장의 모든 채점 규칙, 그리고 이후 28개 장의 모든 손실 함수는 residual 목록을 어떤 방식으로든 눌러 하나의 숫자로 만드는 방법이다. 차이는 어떻게 누르느냐뿐이다.

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라고 말한다. 같은 20개 점 위에서 세 규칙, 세 승자다.

나는 이 세 직선이 서로 의견이 갈리도록 골랐다. 그 점은 분명히 말해야 한다. 핵심은 그것이 얼마나 쉬웠는가다. 그럴듯해 보이는 intercept와 slope를 몇 분만 뒤져도 이런 세 쌍이 수백 개 나온다. 순위는 당신이 고른 규칙의 속성이지 직선들에 대한 사실이 아니다. 그러므로 규칙은 구현 세부사항이 아니다. 그것은 문제의 정의 자체다. 이 장이 답하려는 질문은 여기서 생긴다. 어떤 근거로 그 규칙을 고를 것인가?

먼저 더 작은 문제부터 보자. 직선은 세 개가 아니라 무한히 많기 때문이다. 지금은 모두가 선택하는 squared error를 택하고, 1장에서 perceptron의 11,000 epoch를 줄여 준 trick을 써서 문제를 숫자 하나로 줄이자. 양쪽 열에서 평균을 빼는 것이다. 점구름이 원점 중심으로 옮겨지면, squared error 아래에서 최선의 직선은 정확히 원점을 지난다. 따라서 intercept는 결정되고, 남는 것은 slope 하나뿐이다.

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

후보 slope 601개, 승자 하나: 시간당 0.293밀리미터. 진실은 0.300이다. noise가 섞인 측정값 20개와 for-loop 하나가 시간당 0.01밀리미터 안쪽, 즉 2.33퍼센트까지 접근했다.

흥미로운 부분은 승자가 아니라 탐색의 모양이다. 곡선 전체를 출력하되, 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장은 601개 점을 모두 방문하지 않고 이 골짜기를 내려가는 이야기이고, 골짜기의 바닥이 하나보다 많아질 때 무엇이 달라지는지에 관한 이야기다.

우리가 골짜기를 얻은 것은 제곱했기 때문이다. Absolute error였다면 바닥에 꺾임이 생겼을 것이다. Worst error였다면 직선을 움직여도 아무것도 바뀌지 않는 평평한 구간들이 생겼을 것이다. 제곱은 분명 편리하다. 그리고 대개의 강의가 드는 이유도 대체로 편의성이다. 다만 네 가지 옷을 입고 나온다. error를 양수로 만든다, absolute value도 그렇다. 큰 error를 더 벌준다, 왜 그래야 하는가? differentiable하다, 네제곱도 그렇다. 모두가 쓴다, 그렇긴 하지만 그것은 논증이 아니다.

정직한 입장은 이렇다. Squared error는 직선 B를 골랐고 absolute error는 직선 A를 골랐다. 이 공장에는 둘 중 하나가 맞고 다른 하나는 틀리다. 그러나 지금까지 말한 것으로는 어느 쪽인지 알 수 없다. 규칙을 고르려면 측정값이 어떻게 직선과 달라졌는지에 대해 무언가를 알아야 한다. 그것은 수학에 대한 질문이 아니라 세계에 대한 질문이다. 답하려면 작은 장치 하나가 필요하다.

“어느 직선이 더 나은가”를 답이 있는 질문으로 바꾸는 주장은 다음과 같다.

부품의 폭은 직선에 random error를 더한 것이라고 가정하자. 그리고 그 error가 평균 0, standard deviation $\sigma$인 Gaussian, 즉 bell curve에서 뽑힌다고 가정하자.

Gaussian의 density는 다음과 같다.

이제 perceptron이 할 수 없었던 일을 해 보자. 주어진 후보 slope $\theta$에 대해, 모든 측정값에는 residual이 있다. 위 formula는 그 residual을 숫자 하나로 바꾼다. 이 slope가 진실이라면 정확히 그 크기의 error가 얼마나 plausible한가? 직선 위에 있는 측정값은 큰 숫자를 얻고, 0.5밀리미터 벗어난 측정값은 작은 숫자를 얻는다.

측정값들은 independent하다. 캘리퍼는 직전 부품을 기억하지 않는다. 따라서 product rule에 따르면 노트 전체의 plausibility는 individual density들의 곱이다. 그 곱이 $\theta$의 likelihood다.¹ 방향에 주목하라. Bayes' rule이 다루는 방향이 바로 이것이다. data는 고정되어 있고 알려져 있으며, 변하는 것은 parameter다. 이것은 “slope의 probability”가 아니다. 실제로 얻은 data에 model이 부여하는 probability를 slope의 함수로 읽은 것이다.

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

slope 0.293은 이 노트를 0.25보다 46,000배 더 plausible하게 만들고, 0.35보다 1억 2,700만 배 더 plausible하게 만든다. Maximum likelihood는 실제로 관측한 것이 가능한 한 덜 놀랍게 보이도록 만드는 parameter를 고른다는 원칙이다. theorem은 아니고, “best”가 무엇을 뜻해야 하는지에 대한 제안이다. 하지만 내용이 있는 제안이다. 무엇이든 채점하기 전에 noise에 대한 가정을 말하도록 강제하기 때문이다.

하루가 아니라 한 달치 교대 근무에 같은 세 줄의 code를 실행하면, 방법은 무너진다.

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

2,000번 곱하면 답은 inf다. 상수 하나를 바꿔 보자. 더 엉성한 캘리퍼라서 density가 1보다 크게가 아니라 작게 나온다고 하면, 같은 code는 0.0을 반환한다. 두 답 모두 틀렸고, 방향은 정반대이며, 어느 쪽도 잡을 수 있는 exception을 일으키지 않는다. 두 번째는 warning조차 출력하지 않는다.

수학에는 아무 문제가 없다. 그 설정에서 likelihood는 완벽하게 잘 정의된 finite number다. 자연로그는 1400.91이므로 숫자 자체는 대략 $10^{608}$다. 문제는 당신의 computer에 그 숫자가 없다는 것이다. 그리고 computer가 정확히 어떤 숫자들을 가지고 있는지 이해할 가치가 있다. 이것이 결과를 결정하는 마지막 순간이 아니기 때문이다.

폭발하는 곱을 고치는 방법은 늘 그렇듯 logarithm을 취하는 것이다. logarithm은 곱을 합으로 바꾸고, strictly increasing이므로 maximum의 위치를 옮길 수 없으며, 2,000개의 적당한 숫자의 합은 float64가 불평 없이 처리한다. 관례상 우리는 negative log-likelihood를 취한다. 그러면 더 낫다는 것은 더 작다는 뜻이 된다. 이제 Gaussian density를 대입하고 무슨 일이 일어나는지 보자.

이것이 이 장이 존재하는 이유인 결과이며, 애매하게 말할 필요가 없다. squared error는 관습이 아니다. 상수를 제거한 Gaussian의 negative log-likelihood다. squared error를 minimise하는 것은 error가 Gaussian이라고 주장하고, 어떤 parameter가 data를 가장 덜 놀랍게 만드는지 묻는 행위와 정확히 같다. 당신은 줄곧 그 주장을 하고 있었다. 다만 아무도 말해 주지 않았을 뿐이다.

이 equivalence는 확인할 수 있으니 확인하자. 같은 601개 slope를 full negative log-likelihood, 즉 constants까지 모두 포함한 것과 plain squared error로 scan한다.

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

vertical axis의 숫자는 다르고, 그중 하나는 negative다. sum of squares는 결코 negative가 아니지만 negative log-likelihood는 0 아래로 내려갈 수 있다. density가 1을 넘을 수 있기 때문이다. 그러나 같은 골짜기의 바닥은 마지막 grid point까지 동일하다.

우리는 $\sigma$를 0.12로 제멋대로 고정했다. 그리고 공장 안 누구도 캘리퍼 error의 spread를 모른다. 그것을 두 번째 unknown으로 취급하고 maximum likelihood가 그것도 결정하게 하자. 여기서는 방금 버린 constant term이 돌아온다. model이 perfect precision을 주장하는 것을 막는 유일한 것이기 때문이다.

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

둘은 소수 네 자리까지 일치하며, 우연이 아니다. 그 expression을 differentiate하고 0으로 두면 정확히 $\hat{\sigma}^2 = \frac{1}{N}\sum r_i^2$가 나온다. 따라서 mean squared error는 variance와 단지 비슷한 것이 아니다. 이 model 아래에서는 noise variance의 maximum-likelihood estimate 그 자체다. 당신이 줄곧 minimise하던 숫자는 sensor가 얼마나 noisy한지에 대한 estimate였다.

말하기는 쉽지만 나중에 다시 발견하려면 비싼 주름 하나가 있다. 그 estimate는 낮게 biased되어 있다. residual이 자기 자신을 작게 만들도록 선택된 fit을 기준으로 측정되었기 때문이다. simulate해 보라. true variance가 정확히 1인 distribution에서 각각 20개 측정값을 뽑아 만든 notebook 200,000개를 만들고, fit의 parameter 하나는 측정값 자체에서 estimate한다. sum of squares를 $N$로 나누면 평균은 0.9501이다. $N-1$로 나누면 1.0001이다. 그리고 $(N-1)/N$은 정확히 0.95다. fitting하는 parameter 하나마다 degree of freedom 하나가 든다. 이것은 훨씬 더 큰 문제의 가장 작은 가시적 사례다. model은 언제나 자신이 fitted된 data 위에서 더 좋아 보인다. 4장은 이것을 data를 따로 떼어 두는 discipline으로 바꾸고, 6장은 그 효과에 이름을 붙인다.

squared error가 noise가 Gaussian이라고 주장한다면, 다음 질문은 그 주장이 false일 때 무슨 일이 일어나는가다. 살짝 false가 아니라, 실제 측정값들이 false인 바로 그 방식으로 false일 때다.

작업 현장에서 대부분의 캘리퍼 측정값은 0.1밀리미터 정도까지 좋다. 그리고 교대 근무마다 한두 번은 절삭 부스러기 조각이 jaw 아래에 끼어, 측정값이 몇 밀리미터씩 빗나간다. 이런 error는 heavy-tailed다. 대부분의 시간에는 작지만 가끔 거대하고, bell curve가 허용하는 것보다 훨씬 더 자주 거대하다. Cauchy distribution은 그런 행동의 표준적인 clean model이며, density는 Gaussian만큼 단순하다.

차이는 tail이다. 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의 sample variance는 0.0144, 즉 $0.12^2$에 자리 잡고 그대로 머문다. Cauchy의 것은 올라가고, sample을 계속하는 한 계속 올라간다. converge할 대상이 없기 때문이다. Cauchy distribution에는 variance가 없고 mean도 없다. squared error는 squares의 average를 minimise하는 것이 전부인데, 존재하지 않는 양을 요구받고 있다.

그래서 캘리퍼가 속은 교대 근무 하나를 보자. 같은 20시간, 같은 칼날, 같은 시간당 0.30밀리미터 drift다. 다만 noise만 이제 Cauchy다. 두 번 fit하자. 한 번은 squared residuals를 minimise하고, 한 번은 실제로 data를 생성한 noise의 negative log-likelihood를 minimise한다. 여기서는 centring trick이 도움이 되지 않는다. 그것은 squared error에 대해서만 intercept를 고정해 준다. 그러므로 두 fit 모두 intercept 와 slope의 grid를 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")

강조된 두 줄이 두 fit의 차이 전체다. Cauchy density의 log를 취하고, 이전과 정확히 같은 방식으로 constants를 버리면, $\sum \log\left(1 + (r/s)^2\right)$가 남는다. 같은 recipe, noise에 대한 다른 claim이다.

     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는 drift가 시간당 0.108밀리미터라고 보고한다. 실제 rate의 대략 3분의 1이며, 칼날은 19.6시간까지 괜찮다고 결론낸다. 진짜 답은 11.7시간이다. 그 fit을 믿고 행동하면 공장은 이 분야에서 가장 표준적인 손실 함수의 권위 아래, 허용오차를 벗어난 부품을 만들며 press를 8시간 더 돌린다. Cauchy fit은 같은 20개 측정값, 같은 grid, code 한 줄 차이만으로 11.8시간에 도착한다.

두 가지 반론에는 답할 필요가 있다. 둘 다 좋은 engineer가 가장 먼저 할 말이기 때문이다.

Outlier는 뻔하니 그냥 지워라. 그럴 수 있고, 도움이 되며, 충분하지는 않다. 단 하나의 최악 측정값을 지우면 least-squares slope는 0.108에서 0.239로 이동한다. 그래도 칼날 교체 시점은 13.1시간이고, 한 시간 반 늦다. 최악을 지우고, refit한 다음, 이제 최악인 것을 지우면 0.286에 도달한다. 하지만 이것은 이미 관찰이 아니라 procedure라는 점에 주목하라. 원래 fit의 가장 큰 residual 두 개를 대신 지우면 0.223에 도착한다. 이제 당신은 적어 둘 수도, 방어할 수도 없는 judgement call을 한 것이다. 규칙을 자동화해도 구해지지 않는다. drop-the-largest-residual-then-refit을 simulated shift 1,000개에 실행하면 median slope error는 likelihood fit의 0.0100에 비해 0.0177이고, error가 0.05를 넘는 비율은 14.7%로 likelihood fit의 1.3%에 비해 훨씬 높다. deletion은 잘못된 assumption 위에 얹은 patch다. likelihood에는 patch가 필요 없다. outlier가 impossible하다고 가정한 적이 없기 때문이다.

운 좋은 dataset을 골랐다. 이 반론은 정확히 맞다. 그래서 마지막 experiment는 independent shift 1,000개를 simulate하고 각각에서 두 방식으로 refit한다.

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

Mean이 아니라 median이다. 이 절의 모든 것과 같은 이유다. least-squares error는 Cauchy에 의해 driven되므로 그 average는 보고하기에 stable한 것이 아니다. Least squares는 5번 중 2번 크게 틀린다. likelihood fit은 77번 중 1번 크게 틀리고, 1,000개 shift 전체에서의 최악 실패도 least squares의 최악 실패의 5분의 1보다 작다.

이 중 어느 것도 squared error가 나쁘다는 뜻은 아니다. 그것은 squared error가 specific하다는 뜻이며, 산술은 정확히 그 이유를 말한다. residual 0.1 mm 하나와 7 mm 하나를 생각하자. 제곱하면 나쁜 측정값은 좋은 측정값보다 total에 4,900배 더 많이 기여하므로, 직선 전체가 그것 쪽으로 끌려간다. Cauchy log-likelihood 아래에서는 같은 두 residual이 0.527과 8.133을 기여하며, ratio는 15.4다. 나쁜 측정값은 여전히 반영되지만, 결정권을 갖지는 않는다. 이것이 robust statistics의 시작이다. Huber의 1964년 loss는 작은 residual에는 quadratic하게, 큰 residual에는 linear하게 행동하여 양쪽의 중간을 취했고,⁷ Tukey는 이미 아주 적은 contamination만으로도 sample variance가 mean absolute deviation보다 나쁜 도구가 될 수 있음을 보였다.⁸

역사적 메모 하나도 빼기에는 너무 좋다. Least squares는 Legendre가 1805년에 먼저 발표했다. 그것은 잘 작동한다는 것 외에는 정당화가 없는 편리한 algebraic device였다.⁹ 4년 뒤 Gauss는 논증을 거꾸로 달렸다. 그는 arithmetic mean이 repeated measurement를 결합하는 올바른 방법이라고 받아들이고, 어떤 error distribution이 mean을 가장 probable한 값으로 만드는지 물었으며, 본질적으로 하나만 그렇다는 것을 보였다. 지금 그의 이름이 붙은 바로 그 distribution이다.¹⁰ 이 장의 derivation은 그의 것이다. 200년이 넘었고, 여전히 대부분의 강의가 빠뜨리는 부분이다.

하지만 골짜기를 다시 보라. $\theta = 0.20$에 서 있고 loss가 0.0822라면, “내리막”의 방향은 수수께끼가 아니다. 페이지 위에서 볼 수 있듯 곡선은 오른쪽으로 내려간다. 다른 곳을 전혀 evaluate하지 않고도, loss function에게 당신이 서 있는 그 점에서 어느 방향으로 기울어지는지 물어볼 수 있다면, 그쪽으로 한 걸음 내딛고, 다시 묻고, 땅이 평평해질 때까지 반복할 수 있을 것이다.

그 질문에는 이름이 있다. 한 점에서 함수의 slope는 derivative이고, 많은 parameter를 가진 함수에서는 모든 방향의 slope를 한꺼번에 모은 것이 gradient다. 1장은 그것을 사용할 수 없었다. perceptron의 error가 물어볼 slope가 없는 계단이었기 때문이다. 이 장은 더 나은 것을 만들었다. 어디에서나 smooth하고, preference가 아니라 명시된 assumption에서 나온 loss다.

따라서 3장의 질문은 더 이상 slope가 존재하는지가 아니다. 그것을 어떻게 계산할 것인가, 왜 그 방향 반대로 움직이면 uphill이 아니라 downhill로 가는가, 거의 모든 강의가 믿으라고 요구하는 그 sign의 이유는 무엇인가, 그리고 다시 묻기 전에 얼마나 멀리 step해야 하는가다. 마지막 것은 training run이 converge할지, 답 주변에서 영원히 oscillate할지, 아니면 infinity로 달아날지를 결정하는 단 하나의 숫자로 드러난다.

이 장과 함께 읽을 만한 자료도 있다. Prince, Understanding Deep Learning §5.1–5.2 and Appendix C는 책의 모든 loss를 여기서 쓴 순서대로 maximum likelihood에서 만든다. Goodfellow, Bengio and Courville, Deep Learning §3.1–3.11 and §5.5의 maximum-likelihood section은 4장에 필요한 KL divergence도 유도한다. Murphy, Probabilistic Machine Learning: An Introduction chapter 2 and §4.2는 maximum likelihood가 무엇을 보장하고 무엇을 보장하지 않는지 다룬다. Deisenroth, Faisal and Ong, Mathematics for Machine Learning §6.1–6.4는 sum rule, product rule, Bayes' rule을 제대로 다룬다. Tom Mitchell의 짧은 CMU note Estimating Probabilities: MLE and MAP (2016), 그리고 같은 결과를 runnable code로 도달하는 Dive into Deep Learning §22.7도 함께 볼 만하다.