Откуда берётся loss function: likelihood, а не соглашение

Вот показания и три линии, оценённые тремя способами: squared error, к которому тянутся все; absolute error, который мог бы выбрать статистик; и worst error, который выбрал бы станочник, потому что инспектора не волнует ваше среднее — он бракует одну конкретную деталь, вышедшую за допуск.

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 function в двадцати восьми главах после неё — это некоторый способ сжать список 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

Читайте столбцы, а не строки. Squared error выбирает B, absolute error выбирает A, worst error выбирает C: три правила, три победителя на одних и тех же двадцати точках.

Я специально выбрал эти три линии так, чтобы они не совпали, и должен сказать это прямо. Смысл в том, насколько легко это оказалось: несколько минут поиска по правдоподобным intercept и slopes дают сотни таких троек. Ранжирование — свойство выбранного вами правила, а не факт о линиях, поэтому правило не является деталью реализации: оно и есть определение задачи. Отсюда вопрос, ради которого существует эта глава: на каком основании вы его выбираете?

Сначала вопрос поменьше, потому что линий не три, а бесконечно много. Пока возьмём squared error, раз уж его берут все, и сведём задачу к одному числу с помощью трюка, который сэкономил perceptron одиннадцать тысяч эпох в главе 1: вычтем среднее из обоих столбцов. Когда облако точек центрировано в начале координат, лучшая линия по 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, один победитель: 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 дала бы плоские участки, где сдвиг линии вообще ничего не меняет. Возведение в квадрат, бесспорно, удобно — и удобство примерно и есть причина, которую дают большинство курсов, наряженная четырьмя способами: оно делает ошибки положительными (absolute value тоже); оно сильнее наказывает большие ошибки (а почему должно?); оно дифференцируемо (четвёртая степень тоже); так делают все (да, но это не аргумент).

Вот честная позиция. Squared error выбрала линию B, а absolute error выбрала линию A. Одна из них правильна для этого завода, другая нет, и ничто сказанное до сих пор не позволяет понять, какая именно. Чтобы выбрать правило, нужно знать что-то о том, как показания начали отличаться от линии, а это вопрос о мире, а не о математике. Для ответа нужна одна небольшая деталь механики.

Вот утверждение, которое превращает «какая линия лучше» в вопрос с ответом.

Предположим, что ширина детали — это линия плюс случайная ошибка, и предположим, что эта ошибка взята из Gaussian — колоколообразной кривой — со средним ноль и standard deviation $\sigma$:

Плотность Gaussian равна

Теперь сделаем то, чего perceptron не мог. Для данного кандидатного slope $\theta$ у каждого показания есть residual, и формула выше превращает этот residual в число: насколько правдоподобна ошибка ровно такого размера, если этот slope — истина? Показание на линии получает большое число, показание в полмиллиметра от линии — маленькое.

Показания независимы — штангенциркуль не помнит предыдущую деталь, — поэтому правило произведения говорит, что правдоподобие всего блокнота равно произведению отдельных плотностей. Это произведение и есть likelihood $\theta$.¹ Обратите внимание на направление, потому что именно о нём правило Байеса: данные фиксированы и известны, меняется parameter. Это не «вероятность slope». Это вероятность, которую model назначает данным, которые вы фактически получили, прочитанная как функция от 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, и в сто двадцать семь миллионов раз правдоподобнее, чем 0,35. Maximum likelihood — это принцип, по которому вы выбираете parameter, делающий то, что вы фактически наблюдали, как можно менее удивительным. Это не теорема, а предложение о том, что должно означать «лучший», — предложение содержательное, потому что оно заставляет вас сформулировать предположение о шуме до того, как вам позволено что-либо оценивать.

Запустите те же три строки кода на месяце смен вместо одной, и метод падает.

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, который можно поймать, а второй даже не печатает warning.

С математикой всё в порядке. Likelihood при этих настройках — вполне определённое конечное число: его natural logarithm равен 1400,91, так что само число примерно $10^{608}$. Проблема в том, что у вашего компьютера такого числа нет, и стоит точно понять, какие числа у него есть, потому что это не последний раз, когда он решит исход.

Обычное исправление для взрывающегося произведения — взять логарифмы. Логарифм превращает произведения в суммы, он строго возрастает, поэтому не может сдвинуть положение максимума, а сумма двух тысяч умеренных чисел — то, с чем float64 справляется без жалоб. По соглашению мы берём отрицательный log-likelihood, чтобы «лучше» означало «меньше». Теперь подставьте плотность Gaussian и посмотрите, что произойдёт.

Это результат, ради которого существует глава, и его стоит сформулировать без оговорок: squared error — не соглашение. Это negative log-likelihood Gaussian с отброшенными константами. Минимизировать squared error — ровно то же действие, что утверждать, что ваши ошибки Gaussian, и спрашивать, какой parameter делает ваши данные наименее удивительными. Вы всё это время делали это утверждение; просто вам об этом не говорили.

Эквивалентность проверяема, так что проверьте её: просканируйте те же шестьсот один slopes с полной negative 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

Разные числа на вертикальной оси, и одно из них отрицательно, чего сумма квадратов никогда не делает: negative log-likelihood может уходить ниже нуля, потому что плотность может превышать 1. То же дно той же долины, до последней точки сетки.

Мы зафиксировали $\sigma$ на 0,12 волевым решением, а никто на заводе не знает разброс ошибки штангенциркуля. Рассмотрим его как вторую неизвестную и дадим maximum likelihood выбрать его тоже. Здесь только что отброшенный constant term возвращается, потому что только он стоит между model и заявлением о совершенной точности.

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 не просто похожа на дисперсию. В этой model она и есть maximum-likelihood estimate дисперсии шума — число, которое вы всё это время минимизировали, было оценкой того, насколько шумен ваш sensor.

Одна тонкость, которую легко сформулировать и дорого заново открыть позже: эта estimate смещена вниз, потому что residuals измерялись относительно fit, который сам был выбран так, чтобы сделать их маленькими. Смоделируйте это: двести тысяч блокнотов по двадцать показаний, взятых из распределения с истинной дисперсией ровно 1, причём один parameter fit оценивается по самим показаниям. Деление суммы квадратов на $N$ даёт среднее 0,9501; деление на $N-1$ даёт 1,0001; а $(N-1)/N$ — ровно 0,95. Каждый parameter, который вы подгоняете, стоит одной степени свободы, и это самый маленький видимый экземпляр гораздо более большой проблемы: model всегда выглядит лучше на данных, по которым её подгоняли. Глава 4 превращает это в дисциплину удержания данных в стороне, а глава 6 даёт эффекту имя.

Если squared error утверждает, что шум Gaussian, следующий вопрос — что происходит, когда утверждение ложно. Не слегка ложно, а ложно так, как бывают ложными реальные измерения.

В цеху большинство показаний штангенциркуля точны до одной десятой миллиметра, а раз или два за смену стружка попадает под губку, и показание уходит на несколько миллиметров. Такие ошибки heavy-tailed: чаще всего маленькие, иногда огромные, причём огромные намного чаще, чем допускает колоколообразная кривая. Cauchy distribution — стандартная чистая model такого поведения, и её плотность так же проста, как у 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

Sample variance Gaussian сходится к 0,0144, то есть к $0.12^2$, и остаётся там. У Cauchy она растёт и продолжает расти столько, сколько вы sampling, потому что ей не к чему сходиться: у Cauchy distribution нет variance и нет mean. Squared error, вся работа которой — минимизировать среднее квадратов, просят найти величину, которой не существует.

Итак, вот одна смена, где штангенциркуль обманули. Те же двадцать часов, то же лезвие, тот же drift 0,30 миллиметра в час — только шум теперь Cauchy. Подгоним дважды: один раз минимизируя squared residuals, второй раз минимизируя negative log-likelihood того шума, который действительно сгенерировал данные. Трюк с центрированием здесь не помогает — он фиксирует intercept только для squared error, — поэтому оба fits выполняются грубой силой по сетке intercepts и slopes, ведь у нас всё ещё нет способа найти дно долины иначе, чем посетив его.

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

Две выделенные строки — вся разница между fits. Возьмите логарифм плотности 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 сообщает drift 0,108 миллиметра в час, примерно треть реальной скорости, и заключает, что лезвие годится до часа 19,6. Истинный ответ — час 11,7. Действуя по этому fit, завод гонит пресс ещё восемь лишних часов, производя детали вне допуска, опираясь на самую стандартную loss function в области. Cauchy fit, используя те же двадцать показаний, ту же сетку и разницу в одну строку кода, попадает в час 11,8.

Два возражения заслуживают ответа, потому что оба — первое, что скажет хороший инженер.

Выброс очевиден — просто удалите его. Можно, это помогает, и этого недостаточно. Удаление одного худшего показания сдвигает least-squares slope с 0,108 до 0,239, что всё ещё ставит замену лезвия на час 13,1, на полтора часа позже; удаление худшего, повторная подгонка и удаление того, что стало худшим теперь, доводит вас до 0,286 — и обратите внимание, это уже процедура, а не наблюдение: удалите два крупнейших residuals исходного fit вместо этого, и получите 0,223. Но теперь вы приняли judgment calls, которые нельзя записать или защитить, а автоматизация правила его не спасает: drop-the-largest-residual-then-refit, запущенный на тысяче simulated shifts, имеет median slope error 0,0177 против 0,0100 у likelihood fit и ошибается больше чем на 0,05 в 14,7% смен против 1,3%. Удаление — заплатка поверх неверного предположения. Likelihood не нужна заплатка, потому что он никогда не предполагал, что outlier невозможен.

Вы выбрали удачный dataset. Это возражение ровно верно, поэтому последний experiment simulates тысячу независимых смен и refits обоими способами на каждой.

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, а не mean, по той же причине, что и всё остальное в этом разделе: least-squares errors движимы Cauchy, поэтому их average — неустойчивая величина для отчёта. Least squares сильно ошибается в двух сменах из пяти; likelihood fit сильно ошибается в одной смене из семидесяти семи, а его худший сбой на тысяче смен меньше одной пятой худшего сбоя least squares.

Всё это не делает squared error плохой. Это делает её конкретной, и арифметика точно говорит почему. Возьмите residual 0,1 мм и residual 7 мм. В квадрате плохое показание вносит в итог в 4 900 раз больше, чем хорошее, поэтому линия целиком тянется к нему; в Cauchy log-likelihood те же два residuals вносят 0,527 и 8,133, отношение 15,4. Плохое показание всё ещё учитывается, просто оно не получает права решать. Это начало robust statistics, где loss Хьюбера 1964 года занимает промежуточную позицию: ведёт себя квадратично для малых residuals и линейно для больших,⁷ а Тьюки уже показал, насколько малого contamination достаточно, чтобы sample variance стала худшим инструментом, чем mean absolute deviation.⁸

И одна историческая заметка, слишком хорошая, чтобы её опустить. Least squares первым опубликовал Лежандр в 1805 году как удобный алгебраический приём без обоснования кроме того, что он работал.⁹ Четыре года спустя Гаусс провёл рассуждение в обратную сторону: он принял как данность, что arithmetic mean — правильный способ объединять повторные измерения, спросил, какое распределение ошибок делает mean наиболее вероятным значением, и показал, что по существу только одно — то, которое теперь названо его именем.¹⁰ Вывод в этой главе — его вывод; ему больше двух веков, и именно эту часть большинство курсов всё ещё пропускает.

Но взгляните снова на долину. Стоя в $\theta = 0.20$ с loss 0,0822, направление «вниз» не загадка — вы видите его на странице, кривая спускается вправо. Если бы можно было спросить loss function, куда она наклонена в точке, где вы стоите, не вычисляя её нигде больше, вы могли бы сделать шаг туда, спросить снова и повторять, пока земля не станет плоской.

У этого вопроса есть имя. Наклон функции в точке — её производная, а для функции многих parameters набор наклонов во всех направлениях сразу — gradient. Глава 1 не могла им воспользоваться, потому что ошибка perceptron была лестницей без наклона, о котором можно было спросить. Эта глава построила нечто лучшее: loss, который гладок везде и происходит из явно сформулированного предположения, а не из предпочтения.

Так что вопрос для главы 3 уже не в том, существует ли slope. Вопрос в том, как его вычислить, почему движение против него ведёт вниз, а не вверх — знак, который почти каждый курс просит принять на веру, — и насколько далеко шагнуть перед следующим вопросом; как выясняется, именно это одно число решает, сойдётся ли training run, будет ли вечно колебаться вокруг ответа или убежит в 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 также выводит KL divergence, нужную главе 4; 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 and Bayes' rule; короткую заметку Tom Mitchell из CMU Estimating Probabilities: MLE and MAP (2016); и §22.7 из Dive into Deep Learning, который приходит к тому же результату в исполняемом коде.