Регрессия: предсказать число

Диспетчеру проката нужно число: сколько велосипедов возьмут завтра в восемь утра. Мы построим лесенку моделей на настоящих данных, научимся читать остатки как дневник, поймаем в них выросший на две трети город и выясним, почему хороший прогноз — это не число, а число с подушкой.

Вопрос, который стоит велосипедов

Продолжаем работать с таблицей Capital Bikeshare из урока о режимах обучения: 17 379 часов, для каждого — спрос cnt, погода и календарь. Тогда мы задавали таблице пять разных вопросов; теперь останемся с одним, регрессионным, и доведём его до конца. Ночная смена собирает велосипеды по станциям, и вопрос звучит буднично: сколько аренд случится завтра в каждый час?

Ответ — число, и этим всё сказано о постановке: перед нами регрессия, обучение с учителем, где учителя зовут «счётчик станции». С классификацией из прошлого урока они двоюродные сёстры: там непрерывный счёт резали порогом до класса, здесь число остаётся числом до самого конца, и решение будет резать не выход модели, а распределение будущего. Весь словарь переносится: признаки, честный тест, базовая линия, цена ошибки — поменяется только то, что ошибка теперь имеет величину, а не только знак. Честную проверку устроим по стреле времени, как в уроке о статистике: обучаемся на всём 2011 годе (8645 часов), проверяемся на всём 2012-м (8734 часа). Никакая строка будущего не участвует в настройке — и, как мы увидим к середине урока, именно эта честность вскроет самое интересное свойство наших данных.

Прежде чем строить модели, договоримся, как мерить промах. Остаток часа $t$ — разница факта и прогноза:

$$r_t = y_t - \hat y_t ,$$

и он измеряется в велосипедах: остаток $+120$ значит, что спрос оказался на 120 аренд выше прогноза и у станций стояла очередь. Средний модуль остатка по всем часам теста называется средней абсолютной ошибкой:

$$\mathrm{MAE} =\frac{1}{m}\sum_{t=1}^{m}\bigl|y_t-\hat y_t\bigr| ,$$

и это главная валюта урока: «модель ошибается в среднем на столько-то велосипедов в час».

Лесенка из четырёх моделей

Строить будем от глупого к разумному, и каждая ступень лесенки будет честно измерена на 2012 годе. Все модели лесенки устроены одинаково: разбить часы на ячейки по значениям признаков и отвечать средним своей ячейки по обучающим данным,

$$\hat y(x) =\operatorname{mean}\bigl\{\,y_i:\ \text{cell}(x_i)=\text{cell}(x)\,\bigr\},$$

меняется только дробность разбиения. Это самый честный старт: ни одного подобранного веса, вся сила — в признаках.

Ступень первая: константа. Предскажем для любого часа среднее по обучению, $\hat y=143{,}8$ аренды. MAE на тесте — 168,3 велосипеда. Это наш ноль отсчёта, аналог фильтра-бездельника из урока о классификации: любая модель, не побившая константу, не выучила ничего.

Ступень вторая: у каждого часа суток своё среднее. Мы видели в уроке 06, как город дышит: ночью единицы, в пики сотни. Двадцать четыре числа вместо одного — и MAE падает до 118,2.

Ступень третья: час и тип дня. Рабочий и выходной ритмы различаются формой, поэтому заведём отдельное среднее для каждой пары «час × рабочий день»: 48 чисел. MAE — 108,5.

Ступень четвёртая: добавим погоду. Разделим часы на «сухо» и «осадки» по столбцу weathersit и получим до 96 средних. MAE — 104,4.

$$168{,}3 \;\xrightarrow{\text{час}}\; 118{,}2 \;\xrightarrow{\text{тип дня}}\; 108{,}5 \;\xrightarrow{\text{погода}}\; 104{,}4 .$$

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

Лесенка учит двум вещам сразу. Во-первых, дешёвые табличные модели — «среднее по ячейке» — уже сильны: мы не обучили ни одного веса, а ошибка упала на треть. Во-вторых, вклад признака измеряется только так: разностью ошибок с ним и без него на честном тесте, а не правдоподобием рассуждения «погода наверняка важна». Погода важна — на четыре велосипеда в час.

Взгляд на один день

Средние числа скрывают, как модель живёт внутри дня. Возьмём одну настоящую среду теста, 12 сентября 2012 года, и наложим прогноз третьей модели на факт.

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

Такой график — обязательная часть работы с временными рядами: метрика сжимает 8734 часа в одно число, а два-три дня крупным планом показывают, из чего это число сделано. Форма угадана, уровень занижен — уже диагноз, которого MAE сама по себе не выдаст.

Погода как число: температурный горб

В лесенке погода участвовала грубой категорией «сухо или осадки». Но в таблице есть и температура — непрерывное число, и его связь со спросом поучительна: спрос растёт с теплом лишь до поры. Средний спрос по корзинам температуры поднимается от десятков в мороз к плато около 335 аренд при +28…+32 °C, а за +33 отступает: в тридцатишестиградусной жаре спрос на четверть ниже пикового. Велосипед в пекле нужен немногим.

Зависимость с горбом — приговор для «среднего по ячейке» с крупными корзинами и для прямой линии сразу: прямая обязана расти или убывать всюду, а правда сначала растёт, потом падает. Здесь впервые виден предел линейного мира и причина, по которой урок 50 будет гнуть прямые полиномиальными признаками, а нейросети второй главы — строить горбы из изломов.

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

MAE или RMSE: спор о цене хвостов

Вторая популярная метрика возводит остатки в квадрат:

$$\mathrm{RMSE} =\sqrt{\frac{1}{m}\sum_{t=1}^{m}\bigl(y_t-\hat y_t\bigr)^{2}} ,$$

и у нашей третьей модели она равна 157,1 велосипеда — в полтора раза больше её же MAE 108,5. Расхождение двух метрик само по себе диагностика: квадрат раздувает большие остатки, поэтому RMSE много больше MAE ровно тогда, когда в остатках есть тяжёлые хвосты — редкие, но огромные промахи. У симметричного узкого шума обе метрики почти совпадают.

Какая метрика правильная? Неправильный вопрос; правильный — какая соответствует цене ошибки. Если два промаха по 50 велосипедов стоят столько же, сколько один на 100, ваша метрика MAE. Если один промах на 100 — катастрофа, которую не искупают десять точных часов (сорванный час пик, потерянные клиенты), цена ошибки выпуклая, и RMSE ближе к правде. Выбор метрики — это выбор функции потерь из урока о режимах, продолженный в отчётность.

Среднее и медиана: маленькая теорема с большими последствиями

У спора MAE против RMSE есть точная математическая сердцевина. Спросим: каким одним числом $c$ лучше всего предсказывать выборку $y_1,\ldots,y_n$?

Следствие практично до грубости. Обучая модель на квадратичную потерю, вы просите её предсказывать условное среднее; обучая на модуль — условную медиану. На симметричном шуме они совпадают, но спрос на велосипеды асимметричен: изредка случаются часы-гиганты, и среднее подтягивается к ним вверх, а медиана — нет. Модель «под RMSE» будет систематически ставить прогноз выше типичного часа, страхуясь от гигантов; модель «под MAE» — целиться в типичный час, смиряясь с редкими крупными недолётами. Ни одна не «врёт» — они отвечают на разные вопросы, и заказчик обязан знать, какой из них задал.

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

Прямая по данным: минимум, который считается руками

Раз уж температура — число, проведём через данные первую настоящую обучаемую модель: прямую $\hat y=a+b\,x$, где $x$ — температура. Какие $a$ и $b$ лучшие? По квадратичной метрике ответ выводится так же, как выводилась оптимальная константа: приравнять производные суммы квадратов по $a$ и $b$ к нулю. Получаются формулы метода наименьших квадратов:

$$b=\frac{\sum_i (x_i-\bar x)(y_i-\bar y)} {\sum_i (x_i-\bar x)^{2}}, \qquad a=\bar y-b\,\bar x :$$

наклон — отношение сонаправленности $x$ и $y$ к разбросу $x$, свободный член — подгонка прямой под точку средних. На часах с температурой до +25 °C, где связь ещё монотонна, формулы дают $b=7{,}1$ и $a=73{,}7$: каждый градус тепла приносит примерно семь аренд в час. Это уже не таблица средних — это модель с параметрами, у которых есть смысл и единицы: «аренды в час на градус». МНК мы встречали в уроке о Галилее, где он вытаскивал $g$ из зашумлённых высот; полную линейную регрессию со многими признаками построит урок 49.

Куда исчезают последние сто: пол шума

Лесенка застряла около сотни велосипедов, и стоит спросить: а сколько вообще можно? Разложим ошибку любой «табличной» модели на две части: насколько её ячейковые средние не совпали с истинными и насколько сам час гуляет вокруг истинного среднего своей ячейки:

$$\mathbb E\bigl[(y-\hat\mu_{\text{cell}})^{2}\bigr] =\underbrace{(\mu_{\text{cell}}-\hat\mu_{\text{cell}})^{2}}_{\text{ошибка модели}} +\underbrace{\mathbb E\bigl[(y-\mu_{\text{cell}})^{2}\bigr]}_{\text{шум ячейки}} .$$

Второе слагаемое статичная таблица убрать не может: это внутренний разброс спроса при фиксированных часе и типе дня. Его размер измеряется оракулом: подсмотрим средние ячеек прямо в тестовом году. Такой «идеальный» статичный прогноз даёт MAE 69,1 велосипеда — пол шума для всего класса моделей «постоянное число на ячейку».

У пола есть мелкий шрифт, и он поучителен: пол принадлежит классу моделей, а не задаче. Поправка скользящим множителем из прошлого раздела формально пользуется той же сеткой признаков, но её уровень плывёт во времени — она уже не «постоянное число на ячейку», и её MAE, как вы увидите в лаборатории, проламывает статичный пол. Часть «неустранимого шума» ячейки была на самом деле ростом сервиса внутри года, и модель со временем в руках до неё дотягивается. Правильное чтение пола шума: «точнее нельзя, оставаясь в этом классе моделей»; смена класса двигает и пол. Требовать от подрядчика MAE 30 при статичном поле 69 можно — но только вместе с правом модели смотреть на свежие данные.

R²: метрика, которую спросят на любой защите

В отчётах о регрессии живёт ещё одно число — коэффициент детерминации:

$$R^{2} =1-\frac{\sum_t (y_t-\hat y_t)^2} {\sum_t (y_t-\bar y)^2} ,$$

доля квадратичной ошибки константы, которую модель сумела убрать. Присмотритесь к знаменателю: это ошибка нашей первой ступени. $R^2$ — не что иное, как сравнение с базовой линией, зашитое в формулу: ноль у константы, единица у оракула. Этим он и хорош — безразмерен, сравним между задачами, — и этим же опасен. Он наследует квадрату всю чувствительность к хвостам; он растёт от любого добавленного признака, даже шумового, если мерить на обучении; и он молчит о единицах: модель с $R^2=0{,}9$ может ошибаться на сотню велосипедов там, где решению нужна точность в двадцать. Правило простое: $R^2$ — для сравнения моделей между собой, ошибка в единицах — для разговора о пригодности к делу.

Насколько точна сама ошибка

MAE 104,4 против 108,5 — это победа погоды или случайность удачного года? Вопрос из урока о статистике, и ответ тем же аппаратом: MAE — среднее $m=8734$ модулей, у него есть стандартная ошибка

$$\operatorname{se}(\mathrm{MAE}) =\frac{s_{|r|}}{\sqrt{m}} ,$$

где $s_{|r|}$ — выборочный разброс модулей остатков; для наших моделей она около 1,2 велосипеда. Разность 4,1 велосипеда при таких se — различие реальное, не шум. Но помните оговорку того же урока 05: соседние часы зависимы, тяжёлый день портит модель много часов подряд, поэтому честная неопределённость больше формульной, и аккуратные команды сравнивают модели парно по дням:

$$d_j =\mathrm{MAE}^{(A)}_{\text{день }j} -\mathrm{MAE}^{(B)}_{\text{день }j}, \qquad \bar d\pm 2\,\frac{s_d}{\sqrt{366}} ,$$

той же техникой, что сравнивала прогнозы автобусов.

Остатки — дневник модели

Урок о модели и эксперименте приучил нас читать остатки как письмо от данных. Прочитаем письмо третьей модели. Возьмём шесть самых больших по модулю остатков 2012 года: все шесть положительные, все — недолёты на 570–610 велосипедов, все в вечерние часы пик, и все — в конце лета и сентябре 2012-го. Это не шум: шум не выбирает сезон и знак. Это подпись.

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

Разгадка проста и поучительна: сервис вырос. Средний час 2011 года — 143,8 аренды, 2012-го — 234,7: плюс 63% за год. Прокат открывал станции, город привыкал, спрос рос — а модель, обученная на 2011-м, честно воспроизводит уровень 2011-го. Мы снова встретили сдвиг мира из урока о статистике, только там он был обвалом, а здесь — ростом. Обвал модель замечает мгновенно и катастрофически; рост она недолетает тихо, теряя клиентов по чуть-чуть, и без графика остатков его можно не видеть месяцами.

Чинится сдвиг уровня одной строкой. Оценим коэффициент роста по свежему окну — скажем, по последним четырём неделям — и умножим прогноз:

$$\hat y_t^{\,\text{испр}} =g\cdot\hat y_t, \qquad g=\frac{\overline{y}_{\text{свежие 28 дней}}} {\overline{\hat y}_{\text{свежие 28 дней}}} ,$$

Форму суток по-прежнему помнит 2011 год, уровень — последний месяц. Это та же логика короткой и длинной памяти, что в экспоненциальном сглаживании урока 05: форма меняется медленно, уровень — быстро, и хранить их лучше отдельно.

Лаборатория: соберите прогноз сами

Порядок опытов. Поднимайтесь по лесенке моделей и следите за MAE на 2012 годе и за гистограммой остатков: с каждым признаком она сжимается, но остаётся смещённой вправо. Затем включите поправку роста со скользящим окном и не пропустите главное событие урока: MAE падает со 108 до 51 велосипеда. Один плывущий коэффициент оказался сильнее всех признаков вместе взятых и пробил пол статичной таблицы — свежесть данных здесь стоит больше их ширины. В конце переключитесь с прогноза-числа на прогноз-квантиль и подберите подушку для вечернего пика; сводка покажет, в какую долю часов подушки не хватило.

Число с подушкой: квантильный прогноз

Диспетчеру мало числа 320. Ему нужно решение: сколько велосипедов привезти, чтобы почти наверняка хватило, но без горы лишних. «Почти наверняка» — это не среднее и не медиана, это квантиль уровня $q$: значение $Q_q$, которое спрос не превысит с вероятностью $q$,

$$\Pr\bigl(y\le Q_q\bigr)=q .$$

Подвезти 90-й процентиль спроса, $Q_{0{,}9}$, значит остаться без велосипедов лишь в один час из десяти. Медиана — это $Q_{0{,}5}$, и вся шкала квантилей — это способ говорить о будущем не точкой, а профилем рисков.

Оптимальный уровень подушки выводится из цен, как порог в уроке о классификации. Пусть лишний привезённый велосипед стоит $c_o$ (перегон и парковка), а нехватка одного — $c_u$ (потерянная поездка и злой клиент). Тогда ожидаемые потери минимизирует квантиль уровня

$$q^\star=\frac{c_u}{c_u+c_o} :$$

при $c_u=40$ и $c_o=10$ нужен 80-й процентиль. Формула — близнец порога $\tau^\star$ из прошлого урока, и это не совпадение: классификация и регрессия сшиваются на уровне решений одной и той же арифметикой цен. Задачу о запасе под неопределённый спрос экономисты зовут задачей газетчика; ей больше века, и каждая пекарня решает её ежеутренне.

Научить модель целиться в квантиль можно напрямую, заменив модуль в функции потерь на его перекошенную версию, pinball loss:

$$\ell_q(y,\hat y)= \begin{cases} q\,(y-\hat y), & y\ge \hat y,\\[2pt] (1-q)\,(\hat y-y), & y<\hat y, \end{cases}$$

при $q=0{,}5$ это половина модуля и цель — медиана; при $q=0{,}9$ недолёт штрафуется в девять раз дороже перелёта, и модель выучивается завышать ровно до 90-го процентиля. Три модели с $q=0{,}1$, $0{,}5$ и $0{,}9$ дают не прогноз, а коридор, и коридор для решений полезнее точки.

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

MAPE: ловушка процентов

Есть метрика, которую менеджеры любят больше всех, — средняя абсолютная процентная ошибка:

$$\mathrm{MAPE} =\frac{100\%}{m}\sum_{t} \frac{\bigl|y_t-\hat y_t\bigr|}{y_t} ,$$

«в среднем ошибаемся на столько-то процентов». Звучит универсально, а на наших данных взрывается: в четыре утра спрос равен двум-трём арендам, и промах на пять велосипедов — это ошибка в двести процентов. Ночные часы, самые лёгкие в велосипедах, становятся самыми тяжёлыми в процентах и захватывают метрику; при нулевом спросе формула и вовсе делит на ноль. MAPE честна только вдали от нуля и при сопоставимых масштабах — для спроса с ночными нулями она источник вечных недоразумений. Если проценты необходимы, делить лучше на средний уровень, а не на каждое значение; а внутри команды надёжнее всего говорить в велосипедах.

Логарифм вместо процентов: метрика соревнований

У этой таблицы есть спортивное прошлое: в 2014 году она была задачей соревнования Kaggle Bike Sharing Demand, и тысячи команд мерились на ней силами. Организаторы выбрали метрикой не MAE и не MAPE, а логарифмическую ошибку:

$$\mathrm{RMSLE} =\sqrt{\frac{1}{m}\sum_t \Bigl(\ln(1+\hat y_t)-\ln(1+y_t)\Bigr)^{2}} ,$$

и выбор глубже, чем кажется. Разность логарифмов — это логарифм отношения: метрика штрафует за «во сколько раз», а не за «на сколько». Промах 10 велосипедов ночью при факте 5 весит как промах 200 в пик при факте 500. Это лечит и болезнь MAE, для которой ночь невидима, и болезнь MAPE с делением на ноль: единица под логарифмом страхует нулевые часы. Платит метрика привычным: она безразмерна, её значение не перевести в велосипеды, и модель под RMSLE целится в геометрию, а не в арифметику спроса. Одна и та же лесенка моделей может ранжироваться по MAE и RMSLE по-разному — какая расстановка «правильная», решает не математика, а то, чьи часы дороже: диспетчера пиков или планировщика ночных развозок.

Откуда слово «регрессия»

Название дисциплине подарил курьёз. В 1886 году Фрэнсис Гальтон, сравнивая рост родителей и взрослых детей, обнаружил: у очень высоких отцов сыновья в среднем высоки, но ниже отцов; у очень низких — низки, но выше. Он назвал это «регрессией к посредственности», возвращением к среднему. Никакой мистики: рост — сумма наследственности и случайности, и экстремальное значение обычно означает, что случайность подыграла; в следующем поколении она подыгрывает заново, без памяти о прошлом успехе.

Эффект живёт всюду, где есть отбор по шумной величине. Выберите десять лучших дней проката месяца — следующая неделя у «лучших» окажется скромнее, хотя ничего не сломалось. Ученик после блестящей контрольной чаще пишет следующую чуть хуже, после провальной — чуть лучше, и похвала с руганью тут ни при чём. Не узнав регрессию к среднему в лицо, легко приписать её действию лекарства, реформы или наказания; формальный аппарат для таких ловушек даст урок о нормальном распределении.

Паспорт регрессии

Соберём чек-лист отчёта, парный к паспорту классификатора из прошлого урока. Ошибка в единицах задачи, посчитанная на честном отрезке будущего, рядом — объём теста. Базовая линия: константа и лучшее «среднее по ячейке»; вклад модели — разность с ними. Пара MAE и RMSE вместе: их отношение сообщает о хвостах. График остатков по времени: смещение остатков — сигнал сдвига мира. И форма ответа: точка, коридор или квантиль под конкретное решение, с ценами $c_u$ и $c_o$, если решение о запасе. К пяти строкам — две сноски мелким шрифтом: пол шума, посчитанный оракулом, чтобы обещания не выходили за физику задачи, и парное сравнение по дням для каждого «модель A лучше модели B». Семь пунктов помещаются на страницу, и почти все провалы регрессионных проектов ловятся каким-то из них.

Что осталось за лесенкой

Наши модели были таблицами средних; у них есть потолок. Ячейка «8 утра × рабочий день × дождь» усредняет весь год и не видит ни температуры как числа, ни тренда внутри года, ни взаимодействий тоньше своей сетки; а добавить пятый признак — значит снова разрезать данные, и в ячейках станет пусто. Линейная регрессия $\hat y=\langle w,x\rangle$ обходит потолок, обучая веса непрерывных признаков совместно, и урок 49 выведет её из принципа максимума правдоподобия, а урок 50 научит гнуть прямую полиномиальными признаками. Там же нас ждёт обратная сторона гибкости: модель, которая может выучить что угодно, с удовольствием выучит и шум — это переобучение, главный дракон второй половины курса. Честное деление данных, которым мы пользовались сегодня как рефлексом, станет отдельным уроком.

Сборка: что диспетчер знает теперь

Регрессия отвечает числом и меряется в единицах задачи на честном будущем. Лесенка моделей от константы вверх заменяет спор о важности признаков измерением: час дал пятьдесят велосипедов, погода — четыре, а плывущий уровень — пятьдесят семь, и это главный сюжет урока: свежесть данных победила ширину признаков. Пара MAE/RMSE вскрывает хвосты, а выбор между ними — это выбор между медианой и средним, между вопросами «сколько типично» и «сколько в среднем». Остатки читаются как дневник: их смещение поймало рост сервиса на 63% раньше любого дашборда. Пол шума отделяет «можно улучшить» от «нельзя в этом классе моделей» — и напоминает, что смена класса двигает сам пол. Наконец, решение о запасе требует не точки, а квантиля, и уровень подушки выводится из цен нехватки и излишка той же формулой, что порог классификатора. Прогноз, который нельзя превратить в решение, не закончен; прогноз с коридором и ценами — рабочий инструмент.

Задачи