Логистическая регрессия

Заметка о линейных классификаторах закончилась таблицей того, что оценивает каждая суррогатная функция потерь. Кусочно-линейная (hinge) восстанавливает только сторону границы; логистическая восстанавливает логарифм шансов (log-odds) $\log\frac{\eta}{1-\eta}$, из которого следует и сама вероятность $\eta$. Многим решениям нужна именно эта вероятность, а не только сторона: рекламодатель сравнивает шанс клика с его ценой, клиника соотносит риск пациента со стоимостью анализа, а следующая заметка выбирает пороги исходя из стоимости ошибок. Логистическая регрессия (logistic regression) — линейная модель, построенная для этого. В этой заметке мы выведем её из вероятностной модели метки, обучим методом максимального правдоподобия, найдём её градиент и разберём её поведение на разделимых данных, а затем выясним, когда её вероятностям можно верить буквально.

Обозначения. В этой заметке метки $y_i \in \{0, 1\}$, что удобно для распределения Бернулли; форма через отступ из предыдущей заметки использует метки $\pm 1$, то есть $2y_i - 1$. Признаки $x_i \in \mathbb R^d$ включают константу 1, поэтому свободный член — один из весов $w \in \mathbb R^d$. Оценка (score), или логит, — $z_i = \langle w, x_i\rangle$, сигмоида — $\sigma(z) = 1/(1 + e^{-z})$, а вероятность по модели — $p_i = \sigma(z_i)$. $X$ — матрица размера $n \times d$ со строками $x_i$, $y$ и $p$ — векторы меток и вероятностей, а $\eta(x) = P(y = 1 \mid x)$ — истинная вероятность.

Почему логарифм шансов?#

Линейная функция $\langle w, x\rangle$ принимает любые вещественные значения, а вероятность лежит в $[0, 1]$. Если моделировать линейно саму $p$, получатся «вероятности» $1.3$ и $-0.2$. Шансы (odds) $p/(1 - p)$ лежат в $(0, \infty)$, а логарифм шансов — в $(-\infty, \infty)$, поэтому логарифм шансов — величина, которую линейная модель может описывать, не выходя за пределы её области значений:

$$\log\frac{P(y = 1 \mid x)}{1 - P(y = 1 \mid x)} = \langle w, x\rangle \quad\Longleftrightarrow\quad P(y = 1 \mid x) = \frac{1}{1 + e^{-\langle w, x\rangle}} = \sigma(\langle w, x\rangle) .$$

Сигмоида — просто функция, обратная к логарифму шансов. Ниже используются два её свойства: $\sigma(-z) = 1 - \sigma(z)$ и $\sigma'(z) = \sigma(z)\,(1 - \sigma(z))$.

Почему линейный логарифм шансов естественен. По формуле Байеса в форме шансов апостериорный логарифм шансов равен априорному плюс логарифм отношения правдоподобий $\log\frac{p(x \mid y = 1)}{p(x \mid y = 0)}$. Когда оба класса гауссовские с общей ковариационной матрицей $\Sigma$ и средними $\mu_1, \mu_0$, квадратичные члены двух логарифмов плотностей сокращаются, и

$$\log\frac{p(x \mid y = 1)}{p(x \mid y = 0)} = (\mu_1 - \mu_0)^\top \Sigma^{-1} x + \text{const},$$

что линейно по $x$ (Bishop, 2006, Section 4.2.1; само гауссовское распределение разобрано в заметке о многомерном нормальном распределении). То же происходит и для многих других семейств условных распределений классов. Логистическая регрессия сохраняет вывод — линейный логарифм шансов — и отбрасывает предположение: она моделирует $P(y \mid x)$ напрямую и никак не описывает распределение $x$.

Разобранный пример. Классы $N(1, 1)$ и $N(-1, 1)$ на прямой, равновероятные. Логарифм отношения правдоподобий равен $\big[(x + 1)^2 - (x - 1)^2\big]/2 = 2x$, так что логарифм шансов ровно $2x$, и в точке $x = 0.5$ вероятность класса 1 равна $\sigma(1) = 0.731$. Если класс 1 реже, с априорной вероятностью $0.2$, добавляется априорный логарифм шансов $\log(0.2/0.8) = -1.386$: $\sigma(1 - 1.386) = 0.405$. Изменение базовой частоты сдвигает каждый логарифм шансов на одну и ту же константу — это правило курс теории вероятностей вывел для применения модели при другой распространённости.

Как читать веса. Увеличение признака $j$ на единицу прибавляет $w_j$ к логарифму шансов, то есть умножает шансы на $e^{w_j}$, каковы бы ни были остальные признаки: вес $0.7$ удваивает шансы ($e^{0.7} = 2.01$). Это отношение шансов (odds ratio), которое приводят статьи по медицине и социальным наукам. Оно описывает модель, а причинным оно становится, только если данные допускают причинную интерпретацию.

Метод максимального правдоподобия даёт логистическую функцию потерь#

Каждая метка — бернуллиевская величина с вероятностью успеха $p_i = \sigma(\langle w, x_i\rangle)$, и при данных признаках метки независимы. Правдоподобие обучающих меток равно $\prod_i p_i^{y_i} (1 - p_i)^{1 - y_i}$, а его отрицательный логарифм —

$$\mathcal L(w) = -\sum_{i=1}^n \Big[ y_i \log p_i + (1 - y_i) \log(1 - p_i) \Big] ,$$

логистическая функция потерь (log-loss), или бинарная кросс-энтропия (binary cross-entropy, BCE). То, что кросс-энтропийные функции потерь — отрицательные логарифмы правдоподобия категориальной модели, выведено в заметке о методе максимального правдоподобия; теоретико-информационное прочтение — неожиданность наблюдаемых меток при убеждениях модели — дано в заметке об энтропии из курса генеративных моделей.

Это и есть логистическая суррогатная функция. Для $y_i = 1$ слагаемое равно $-\log\sigma(z_i) = \log(1 + e^{-z_i})$. Для $y_i = 0$ оно равно $-\log(1 - \sigma(z_i)) = -\log\sigma(-z_i) = \log(1 + e^{z_i})$. С отступом $M_i = (2y_i - 1)\, z_i$ оба случая записываются как

$$\mathcal L(w) = \sum_{i=1}^n \log\big(1 + e^{-M_i}\big) ,$$

то есть логистическая суррогатная функция из предыдущей заметки с натуральным логарифмом. Значит, суррогатная функция, выбранная там за свою форму, — это в точности метод максимального правдоподобия в модели линейного логарифма шансов.

Разобранный пример. Один признак, $w = 1$, без свободного члена, три точки: $x = (2, -1, 0.5)$ с метками $y = (1, 0, 0)$. Вероятности равны $p = \sigma(x) = (0.881, 0.269, 0.622)$. Потери: $-\log 0.881 = 0.127$ для уверенно верной точки, $-\log(1 - 0.269) = 0.313$ для умеренно верной и $-\log(1 - 0.622) = 0.974$ для третьей точки, которую модель помещает не на ту сторону; сумма равна $1.414$. Модель, которая сказала бы $p = 0.99$ для третьей точки, заплатила бы за неё $-\log 0.01 = 4.6$: логистическая функция потерь штрафует уверенные ошибки без ограничения.

Градиент: остатки, умноженные на входы#

Продифференцируем одно слагаемое по его оценке. Для $y_i = 1$ имеем $\frac{d}{dz}\big[-\log\sigma(z)\big] = -(1 - \sigma(z))$; для $y_i = 0$ имеем $\frac{d}{dz}\big[-\log(1 - \sigma(z))\big] = \sigma(z)$; в обоих случаях получается $\sigma(z_i) - y_i$. По правилу дифференцирования сложной функции, с $\partial z_i / \partial w = x_i$,

$$\nabla \mathcal L(w) = \sum_{i=1}^n \big(\sigma(z_i) - y_i\big)\, x_i = X^\top (p - y) .$$

Форма та же, что у метода наименьших квадратов, $X^\top(Xw - y)$: предсказание минус целевое значение, умноженное на вход. Точка тянет веса пропорционально тому, насколько неверна её вероятность, а верно и уверенно классифицированная точка, у которой $p_i \approx y_i$, почти не тянет. В разобранном примере $p - y = (-0.119, 0.269, 0.622)$, градиент по $w$ равен $-0.119 \cdot 2 + 0.269 \cdot (-1) + 0.622 \cdot 0.5 = -0.196$, а по свободному члену он был бы $\sum_i (p_i - y_i) = 0.772$: предсказания в среднем завышены (их среднее $0.591$, доля положительных $1/3$), и шаг градиентного спуска уменьшает свободный член.

Выпуклость. Ещё одно дифференцирование, с $\sigma' = \sigma(1 - \sigma)$, даёт гессиан

$$\nabla^2 \mathcal L(w) = \sum_{i=1}^n p_i (1 - p_i)\, x_i x_i^\top = X^\top S X, \qquad S = \operatorname{diag}\big(p_i (1 - p_i)\big) .$$

Для любого вектора $v$ выполнено $v^\top X^\top S X v = \sum_i p_i(1 - p_i) \langle x_i, v\rangle^2 \ge 0$, поэтому логистическая функция потерь выпукла, а при полном столбцовом ранге $X$ строго выпукла (при конечных $w$ каждое $p_i$ лежит строго между 0 и 1). Каждый локальный минимум глобален.

Замкнутой формулы нет, но есть метод Ньютона. Приравнивание градиента к нулю даёт $X^\top(\sigma(Xw) - y) = 0$ — уравнения, нелинейные по $w$, в отличие от нормальных уравнений метода наименьших квадратов. Их решают итеративно. Градиентный спуск работает; при нескольких сотнях признаков метод Ньютона (Newton’s method) сходится за считанные шаги:

$$w \leftarrow w - (X^\top S X)^{-1} X^\top (p - y) .$$

Каждый шаг — задача взвешенных наименьших квадратов с весами $S$ и рабочим откликом (working response) $Xw + S^{-1}(y - p)$, поэтому метод называется итеративно перевзвешенным методом наименьших квадратов (iteratively reweighted least squares, IRLS; Bishop, 2006, Section 4.3.3). На данных о раке молочной железы с двумя стандартизованными признаками 25 шагов Ньютона доводят градиент до $10^{-15}$ и дают те же веса, что и scikit-learn (раздел «Код»).

Что гарантирует оптимум. Свободный член — вес постоянного признака, поэтому уравнение градиента для него имеет вид $\sum_i (p_i - y_i) = 0$: в решении метода максимального правдоподобия средняя предсказанная вероятность равна доле положительных в обучающей выборке. То же верно внутри любой группы, выделенной бинарным признаком модели. В модели для рака молочной железы предсказанные вероятности злокачественности в сумме дают $212.0$ — ровно число злокачественных опухолей.

Рисунок

Логистическая регрессия: граница и линии уровня вероятности

Рисунок показывает геометрию. Вероятность зависит от $x$ только через оценку $z = \langle w, x\rangle$, поэтому её линии уровня — гиперплоскости $z = \operatorname{logit}(p)$, все параллельные границе. Сдвиг вдоль $w$ на расстояние $1/\lVert w\rVert$ меняет логарифм шансов на единицу; чем больше $\lVert w\rVert$, тем круче переход от «наверняка 0» к «наверняка 1». Ширина полосы между $p = 0.1$ и $p = 0.9$, $2\log 9/\lVert w\rVert$, наглядно показывает, насколько модель уверена.

Разделимые данные: веса уходят на бесконечность#

Если некоторый $w_0$ классифицирует обучающую выборку безошибочно, все отступы $M_i = (2y_i - 1)\langle w_0, x_i\rangle$ положительны, и умножение $w_0$ на $c > 1$ умножает их все на $c$. Каждое слагаемое $\log(1 + e^{-cM_i})$ убывает к нулю с ростом $c$, тогда как $\mathcal L(w) > 0$ при любом конечном $w$. Инфимум, ноль, не достигается: оценки максимального правдоподобия не существует. Она существует ровно тогда, когда классы перекрываются, то есть не разделены ни полностью, ни квазиполностью (complete / quasi-complete separation; Albert and Anderson, 1984).

Градиентный спуск на разделимых данных не останавливается. Его направление стабилизируется — к разделяющей гиперплоскости с максимальным отступом из заметки о методе опорных векторов (Soudry et al., 2018), — а норма растёт неограниченно, хотя лишь логарифмически по числу шагов, и все предсказанные вероятности уплывают к 0 и 1. На четырёх точках $x = -2, -1, 1, 2$, из которых две правые положительны, градиентный спуск с шагом $0.5$ из нуля даёт наклон $5.55$ после $10^3$ шагов, $7.83$ после $10^4$ и $10.13$ после $10^5$: примерно $\log(t/4)$, так что каждое десятикратное увеличение числа шагов добавляет к весу $\log 10 = 2.30$ и делит функцию потерь, примерно $2/t$, на десять.

Разделимость — не экзотика. Как показала предыдущая заметка, случайные метки на $n$ точках в размерности $d$ разделимы с вероятностью $1/2$ при $n = 2(d + 1)$, а реальные данные далеко не случайны. Данные о раке молочной железы со всеми 30 стандартизованными признаками линейно разделимы: при ослаблении штрафа норма весов растёт с $3.8$ до $1\,891$, а логистическая функция потерь на обучении падает с $0.053$ до $0.004$ (раздел «Код»).

Штраф это исправляет. Минимизация

$$\mathcal L(w) + \frac{\lambda}{2}\lVert w\rVert^2$$

добавляет к гессиану $\lambda I$, поэтому целевая функция строго выпукла и неограниченно растёт по всем направлениям: единственная конечная точка минимума существует при любых данных. Штраф — это гауссовское априорное распределение $w \sim N(0, \lambda^{-1} I)$, а решение со штрафом — его MAP-оценка; свободный член обычно не штрафуют. Чем больше $\lambda$, тем меньше $\lVert w\rVert$ и тем дальше друг от друга линии уровня вероятности. LogisticRegression из scikit-learn по умолчанию применяет этот штраф, где C = 1 — величина, обратная его силе; C = np.inf даёт чистый метод максимального правдоподобия.

Калибровка: значат ли вероятности то, что говорят?#

Модель откалибрована (calibrated), если среди примеров, которым она приписывает вероятность $p$, доля положительных равна $p$:

$$P\big(y = 1 \mid \hat p(x) = p\big) = p \qquad \text{для всех } p .$$

Синоптик, который говорит «вероятность дождя 70%», откалиброван, если дождь идёт в 70% таких дней. Калибровку проверяют диаграммой надёжности (reliability diagram): тестовые предсказания разбивают на корзины (часто на десять корзин равной ширины) и для каждой корзины откладывают среднее предсказание против доли положительных в ней. У откалиброванной модели точки лежат на диагонали. Ожидаемая ошибка калибровки (expected calibration error, ECE) сводит расстояние к одному числу:

$$\mathrm{ECE} = \sum_{b} \frac{n_b}{n} \,\big\lvert \bar y_b - \bar p_b \big\rvert ,$$

где в корзине $b$ лежат $n_b$ примеров со средним предсказанием $\bar p_b$ и долей положительных $\bar y_b$. Это число зависит от корзин и от размера выборки. Доля положительных в каждой корзине шумна, поэтому даже у идеально откалиброванных предсказаний, с метками $y \sim \mathrm{Bernoulli}(p)$, ECE положительна: на 10 000 примеров она в среднем равна $0.010$ при десяти корзинах и $0.022$ при пятидесяти. Поэтому диаграмма надёжности должна показывать число примеров в каждой корзине, а ECE стоит читать рядом с логистической потерей и мерой Брайера ниже, которым корзины не нужны. Одна калибровка — ещё не качество: предсказание базовой частоты для всех идеально откалибровано и бесполезно. Хороший вероятностный классификатор откалиброван и разделяет классы; последнее измеряет следующая заметка.

Почему логистическая функция потерь нацелена на откалиброванные вероятности. В отдельной точке с истинной вероятностью $\eta$ предсказание $p$ имеет ожидаемую логистическую потерю $-\eta\log p - (1 - \eta)\log(1 - p)$. Её производная $-\eta/p + (1 - \eta)/(1 - p)$ обращается в ноль при $p = \eta$: ожидаемая потеря минимальна, когда говорят правду. Функции потерь с этим свойством называются правильными оценочными правилами (proper scoring rules; Gneiting and Raftery, 2007). Квадратичная ошибка $(y - p)^2$, мера Брайера (Brier score; Brier, 1950), тоже правильная; абсолютная ошибка $\lvert y - p\rvert$ — нет, поскольку её математическое ожидание $\eta(1 - p) + (1 - \eta)p$ линейно по $p$ и минимально при $0$ или $1$. Правильность — утверждение о наилучшем возможном предсказателе. Обученная логистическая регрессия откалибрована в целом (in the large) на своей обучающей выборке — по уравнению для свободного члена выше — и откалибрована всюду лишь настолько, насколько истинный логарифм шансов линеен, а выборка велика. Переобучение делает её переуверенной, сильный штраф — недостаточно уверенной.

Другие модели даже не целятся. Оценка, обученная с кусочно-линейной функцией потерь, вообще не вероятность. Наивный байесовский классификатор перемножает свидетельства коррелированных признаков, как будто они независимы, учитывает один и тот же сигнал несколько раз и обычно толкает свои вероятности к 0 и 1. Методы с максимальным отступом и бустинг над деревьями обычно отталкивают их от 0 и 1, давая S-образные кривые надёжности (Niculescu-Mizil and Caruana, 2005). Современные глубокие сети, обученные до нулевой ошибки на обучении, склонны к переуверенности (Guo et al., 2017). Это наблюдаемые тенденции семейств моделей, а не гарантии. В разделе «Код» наивный байесовский классификатор на гауссовских классах с коррелированными признаками имеет ECE $0.236$ и логистическую потерю $1.133$, а логистическая регрессия, которая там является верной моделью, — $0.012$ и $0.633$. На этих 10 000 тестовых примерах значения ECE около $0.01$ находятся на уровне идеально откалиброванных предсказаний, так что хорошо откалиброванные модели по ним между собой не упорядочить.

Перекалибровка. Оценку, которая хорошо ранжирует, но плохо откалибрована, можно исправить монотонным отображением, обученным на данных, не использованных при обучении модели (на отложенной выборке или кросс-валидацией, как в заметке о валидации):

• Масштабирование Платта (Platt scaling; Platt, 1999) обучает методом максимального правдоподобия логистическую регрессию с одним признаком — оценкой: $\hat p = \sigma(a s + b)$, где $s$ — оценка модели или логит её вероятности. Чтобы ограничить переобучение на малых выборках, Платт заменил целевые значения $1$ и $0$ на $(N_+ + 1)/(N_+ + 2)$ и $1/(N_- + 2)$, где $N_\pm$ — размеры классов. Температурное масштабирование (temperature scaling) нейронных сетей — частный случай $b = 0$, один множитель при логитах (Guo et al., 2017).
• Изотоническая регрессия (isotonic regression; Zadrozny and Elkan, 2002) подбирает наилучшую по квадратичной ошибке неубывающую ступенчатую функцию от оценки алгоритмом объединения соседних нарушителей (pool-adjacent-violators; Ayer et al., 1955): идём по примерам, упорядоченным по оценке, и всякий раз, когда среднее блока больше среднего следующего блока, объединяем их. Для шести меток в порядке оценки $0, 1, 0, 0, 1, 1$ вторая метка больше третьей, объединённый блок $(1, 0)$ со средним $0.5$ больше четвёртой метки, и результат — $0, \tfrac13, \tfrac13, \tfrac13, 1, 1$. Она исправляет любое монотонное искажение, но ей нужно больше данных: если примеров для калибровки меньше, чем от нескольких сотен до тысячи, масштабирование Платта обычно работает лучше (Niculescu-Mizil and Caruana, 2005).

Отображение Платта строго возрастает (при $a > 0$) и не меняет ранжирование, а значит, и ROC-кривую из следующей заметки. Изотоническое отображение лишь неубывающее: оно может слить соседние оценки в одну, что мало меняет AUC. Оба утверждения относятся к одной фиксированной оценке. CalibratedClassifierCV с cv=5, как в разделе «Код», обучает модель и её отображение на каждом из пяти фолдов и усредняет пять откалиброванных выходов; среднее разных моделей может упорядочить примеры иначе, чем любая из них. В разделе «Код» AUC наивного байесовского классификатора равна $0.695$ до и после любой из двух перекалибровок.

Рисунок

Диаграмма надёжности: переуверенная модель до и после калибровки

Несколько классов: softmax#

При $K$ классах у каждого класса свой вектор весов $w_k$ и своя оценка $z_k = \langle w_k, x\rangle$, а оценки превращаются в вероятности функцией softmax:

$$P(y = k \mid x) = \frac{e^{z_k}}{\sum_{j=1}^K e^{z_j}} .$$

Метод максимального правдоподобия даёт кросс-энтропийную функцию потерь $-\sum_i \log P(y_i \mid x_i)$, и её градиент имеет ту же форму через остатки, $X^\top(P - Y)$, где $P$ — матрица размера $n \times K$ предсказанных вероятностей, а $Y$ — метки в виде one-hot-векторов. Важны только разности оценок: прибавление одного и того же вектора ко всем $w_k$ ничего не меняет, поэтому один класс можно зафиксировать, $w_K = 0$. При $K = 2$ вероятность класса 1 по softmax равна $e^{z_1}/(e^{z_1} + e^{z_2}) = \sigma(z_1 - z_2)$, то есть это бинарная логистическая регрессия с $w = w_1 - w_2$. Заметка о многоклассовой классификации развивает эту модель и её альтернативы.

Где логистическая регрессия появляется снова#

• Первый текстовый классификатор курса NLP — многоклассовая логистическая регрессия на числах вхождений слов, модель мешка слов; она достигает доли правильных ответов (accuracy) 0.918 на AG News и остаётся базовой моделью в заметке о классификации текстов.
• Негативное сэмплирование в word2vec обучает логистическую регрессию, которая отличает наблюдаемые пары (слово, контекст) от случайных, а выходной слой языковой модели — softmax по словарю, обучаемый с кросс-энтропией, как в заметке об n-граммах.
• Дискриминатор GAN — классификатор, обучаемый ровно с этой логистической функцией потерь (заметка о GAN), а шумоконтрастное оценивание выучивает ненормированную плотность логистической регрессией данных против шума (заметка о моделях, основанных на энергии).

Что изменилось#

	Линейный классификатор с суррогатной функцией потерь	Логистическая регрессия
Выход	сторона гиперплоскости	вероятность $\sigma(\langle w, x\rangle)$
Функция потерь	выбрана за форму	выведена: отрицательный логарифм правдоподобия Бернулли
Смысл оценки	зависит от функции потерь	логарифм шансов; $e^{w_j}$ — отношение шансов
Обучение	градиентные методы	то же, плюс метод Ньютона (IRLS) с $\nabla = X^\top(p - y)$
Разделимые данные	много решений	оценки максимального правдоподобия нет; штраф (гауссовское априорное распределение) её возвращает
Вероятности	недоступны	откалиброваны, если модель верна; проверяются диаграммой надёжности, исправляются отображениями Платта или изотоническими

Ограничения#

• Линейный логарифм шансов. Взаимодействия и нелинейные эффекты нужно подавать в виде признаков; иначе модель плохо откалибрована там, где истинный логарифм шансов изгибается.
• Разделимость и почти разделимость. Решения без штрафа расходятся или дают огромные неустойчивые коэффициенты; тогда важна сила штрафа, и её нужно подбирать на валидации, а признаки приводить к сопоставимым масштабам.
• Коэффициенты ненадёжно интерпретировать. Коррелированные признаки делят вес и перекидывают его друг другу, а отношение шансов условно относительно всех остальных признаков модели.
• Калибровка — не разделение. Откалиброванная модель может всё равно плохо разделять классы, а калибровка, измеренная на одной генеральной совокупности, ломается при сдвиге базовой частоты или распределения признаков.
• Редкие классы. Вероятности могут быть верными, а порог 0.5 по умолчанию почти никогда не срабатывает; выбор порога — отдельное решение.

Что нужно запомнить#

• Логистическая регрессия моделирует логарифм шансов линейно: $P(y = 1 \mid x) = \sigma(\langle w, x\rangle)$. Гауссовские классы с общей ковариационной матрицей дают ровно эту форму.
• Метод максимального правдоподобия в модели Бернулли даёт логистическую функцию потерь, которая совпадает с логистической суррогатной функцией $\sum_i \log(1 + e^{-M_i})$.
• Градиент равен $X^\top(p - y)$, гессиан — $X^\top S X$ с $S = \operatorname{diag}(p_i(1 - p_i))$; функция потерь выпукла, точки минимума в замкнутом виде нет, а метод Ньютона — это IRLS.
• В оптимуме при наличии свободного члена средняя предсказанная вероятность равна доле положительных.
• На разделимых данных веса растут неограниченно, а вероятности уходят к 0 и 1; L2-штраф, то есть гауссовское априорное распределение, даёт единственное конечное решение.
• Логистическая функция потерь — правильное оценочное правило. Калибровку проверяют диаграммой надёжности с числом примеров в корзинах и ECE рядом с логистической потерей или мерой Брайера, исправляют масштабированием Платта или изотонической регрессией, обученными на отложенных данных. Масштабирование Платта сохраняет ранжирование; изотоническая регрессия может лишь добавить совпадения оценок.
• Softmax с кросс-энтропией — версия для $K$ классов; при $K = 2$ она сводится к сигмоиде.

Код#

Гауссовский логарифм шансов, метод Ньютона, разделимость, калибровка наивного байесовского классификатора, объединение соседних нарушителей и softmax для двух классов:

import numpy as np
from scipy.special import expit as sigmoid, softmax
from scipy.stats import norm
from sklearn.datasets import load_breast_cancer
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.naive_bayes import GaussianNB
from sklearn.calibration import CalibratedClassifierCV
from sklearn.isotonic import IsotonicRegression
from sklearn.metrics import log_loss, roc_auc_score

rng = np.random.default_rng(0)

# 1. гауссовские классы N(+1, 1) и N(-1, 1): логарифм шансов ровно 2x; априорная вероятность 0.2 добавляет log(0.25)
x0 = 0.5
log_odds = norm.logpdf(x0, 1) - norm.logpdf(x0, -1)
print(log_odds, round(sigmoid(log_odds), 4), round(sigmoid(log_odds + np.log(0.2 / 0.8)), 4))   # 1.0 0.7311 0.4046

# 2. максимальное правдоподобие методом Ньютона (IRLS) на двух стандартизованных признаках breast cancer
data = load_breast_cancer()
X = StandardScaler().fit_transform(data.data[:, :2])        # средний радиус, средняя текстура
y = 1.0 - data.target                                       # 1 = злокачественная
Xb = np.c_[np.ones(len(X)), X]
w = np.zeros(3)
for _ in range(25):
    p = sigmoid(Xb @ w)
    grad = Xb.T @ (p - y)                                   # градиент X^T (sigma - y)
    H = Xb.T @ (Xb * (p * (1 - p))[:, None])                # гессиан X^T S X
    w -= np.linalg.solve(H, grad)
sk = LogisticRegression(C=np.inf, solver='newton-cholesky').fit(X, y)
print(w.round(4), np.r_[sk.intercept_, sk.coef_[0]].round(4))   # [-0.7076  3.722   0.9374] дважды
print(f'{np.abs(grad).max():.1e}', round(p.sum(), 4), y.sum())  # 1.2e-15 212.0 212.0: среднее предсказание = доля положительных

# 3. все 30 признаков: обучающая выборка разделима, и функция потерь без штрафа уходит к нулю
X30 = StandardScaler().fit_transform(data.data)
for C in (1.0, 1e2, 1e4, 1e6):
    m = LogisticRegression(C=C, max_iter=100_000).fit(X30, y)
    p = m.predict_proba(X30)[:, 1]
    print(C, round(np.linalg.norm(m.coef_), 1), f'{log_loss(y, p):.2g}', int((m.predict(X30) != y).sum()))
# C = 1: |w| 3.8, потеря 0.053, 7 ошибок | 1e2: 21.4, 0.03, 4 | 1e4: 113.8, 0.022, 4 | 1e6: 1891.1, 0.0042, 0

# 4. градиентный спуск на четырёх разделимых точках: вес растёт как log(t/4), потеря падает как 2/t
Xs, ys = np.array([[1., -2], [1, -1], [1, 1], [1, 2]]), np.array([0, 0, 1, 1.])
w = np.zeros(2)
for t in range(1, 100_001):
    w -= 0.5 * Xs.T @ (sigmoid(Xs @ w) - ys) / 4
    if t in (1_000, 10_000, 100_000):
        print(t, round(w[1], 3), f'{np.mean(np.logaddexp(0, -(2 * ys - 1) * (Xs @ w))):.2e}')
# 1000 5.548 1.95e-03 | 10000 7.828 1.99e-04 | 100000 10.127 2.00e-05

# 5. калибровка: гауссовские классы с коррелированными признаками; наивный Байес учитывает коррелированные свидетельства дважды
d, rho = 8, 0.7
L = np.linalg.cholesky(rho * np.ones((d, d)) + (1 - rho) * np.eye(d))
yc = rng.integers(0, 2, 20_000)
Xc = rng.normal(size=(20_000, d)) @ L.T + np.where(yc[:, None] == 1, 0.3, -0.3)
tr, te = slice(0, 10_000), slice(10_000, None)
def ece(p, y, bins=10):
    k = np.minimum((p * bins).astype(int), bins - 1)
    return sum((k == b).mean() * abs(p[k == b].mean() - y[k == b].mean()) for b in range(bins) if (k == b).any())
models = {'naive Bayes': GaussianNB(), '+ Platt': CalibratedClassifierCV(GaussianNB(), method='sigmoid', cv=5),
          '+ isotonic': CalibratedClassifierCV(GaussianNB(), method='isotonic', cv=5), 'logistic regression': LogisticRegression()}
for name, m in models.items():
    p = m.fit(Xc[tr], yc[tr]).predict_proba(Xc[te])[:, 1]
    print(name, round(ece(p, yc[te]), 3), round(log_loss(yc[te], p), 3), round(roc_auc_score(yc[te], p), 3))
# ECE, log-loss, AUC: naive Bayes 0.236 1.133 0.695 | + Platt 0.008 0.64 0.695 | + isotonic 0.01 0.634 0.695
#                     logistic regression 0.012 0.633 0.695

# 6. изотоническая регрессия (объединение соседних нарушителей) на шести метках, упорядоченных по оценке
print(IsotonicRegression().fit_transform(np.arange(6), [0, 1, 0, 0, 1, 1]).round(3))   # [0. 0.333 0.333 0.333 1. 1.]

# 7. softmax для двух классов = сигмоида от разности оценок
z = np.array([1.3, -0.4])
print(softmax(z)[0].round(4), sigmoid(z[0] - z[1]).round(4))    # 0.8455 0.8455

# 8. ECE идеально откалиброванных предсказаний (y ~ Bernoulli(p)) на 10 000 примеров, 200 повторений, 10 и 50 корзин
p0 = rng.random((200, 10_000)); y0 = (rng.random((200, 10_000)) < p0).astype(float)
print([round(float(np.mean([ece(p, y, b) for p, y in zip(p0, y0)])), 3) for b in (10, 50)])   # [0.01, 0.022]

Дальше#

Логистическая регрессия возвращает вероятность, а вероятность — ещё не решение. Где-то порог превращает её в «пометить» или «пропустить», и правильный порог зависит от того, насколько редки положительные примеры и сколько стоит каждый вид ошибки. Чтобы судить о модели, тоже нужно число, а доля правильных ответов (accuracy), очевидный кандидат, почти ничего не говорит, когда положительных мало: модель, которая никогда ничего не помечает, набирает 99%, когда положительных 1%. Следующая заметка строит словарь для этих вопросов: матрица ошибок, точность и полнота, ROC-кривая и кривая точность–полнота, пороги, выбранные по стоимости ошибок.

Литература#

• D. R. Cox. The regression analysis of binary sequences. Journal of the Royal Statistical Society, Series B, 20(2), 215–232, 1958.
• J. Berkson. Application of the logistic function to bio-assay. Journal of the American Statistical Association, 39(227), 357–365, 1944.
• A. Albert, J. A. Anderson. On the existence of maximum likelihood estimates in logistic regression models. Biometrika, 71(1), 1–10, 1984.
• D. Soudry, E. Hoffer, M. S. Nacson, S. Gunasekar, N. Srebro. The implicit bias of gradient descent on separable data. Journal of Machine Learning Research, 19(70), 1–57, 2018.
• J. C. Platt. Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods. In A. Smola, P. Bartlett, B. Schölkopf, D. Schuurmans (eds.), Advances in Large Margin Classifiers, 61–74. MIT Press, 1999.
• B. Zadrozny, C. Elkan. Transforming classifier scores into accurate multiclass probability estimates. KDD, 694–699, 2002.
• M. Ayer, H. D. Brunk, G. M. Ewing, W. T. Reid, E. Silverman. An empirical distribution function for sampling with incomplete information. The Annals of Mathematical Statistics, 26(4), 641–647, 1955.
• A. Niculescu-Mizil, R. Caruana. Predicting good probabilities with supervised learning. ICML, 625–632, 2005.
• C. Guo, G. Pleiss, Y. Sun, K. Q. Weinberger. On calibration of modern neural networks. ICML, 2017.
• T. Gneiting, A. E. Raftery. Strictly proper scoring rules, prediction, and estimation. Journal of the American Statistical Association, 102(477), 359–378, 2007.
• G. W. Brier. Verification of forecasts expressed in terms of probability. Monthly Weather Review, 78(1), 1–3, 1950.
• C. M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006. Раздел 4.2 (генеративные модели с линейным логарифмом шансов), разделы 4.3.2–4.3.4 (логистическая регрессия, IRLS, многоклассовый случай).
• T. Hastie, R. Tibshirani, J. Friedman. The Elements of Statistical Learning, 2nd ed. Springer, 2009. Раздел 4.4 (логистическая регрессия).
• E. Sokolov et al. Машинное обучение (HSE): лекция 5 (логистическая регрессия и логистическая функция потерь), лекция 6 (многоклассовая логистическая регрессия), семинар 6: калибровка.