Линейная регрессия

В Алматы с 1941 по 2025 год было 31 046 дней, и для каждого известна средняя температура. Теплеет ли город, и если да, то насколько быстро? На глаз в этом облаке ничего не разобрать: январь отличается от июля почти на тридцать градусов, а искомая величина — десятые доли градуса за десятилетие. Линейная регрессия отвечает одной системой из шести уравнений. Ниже она решается прямо в браузере: включайте и выключайте столбцы модели, тяните края окна лет и смотрите, как меняется ответ.

Выключите годовую волну. Модель «среднее плюс тренд» объясняет одну тысячную разброса, R² = 0,001, а её типичный промах — те же десять градусов, что и у одного только среднего. Включите волну обратно: R² подскакивает до 0,85, промах падает до 3,9 градуса. А тренд почти не меняется: +0,13 и +0,12 градуса за десять лет. Времена года и потепление не мешают друг другу, и регрессия находит их одновременно. Теперь нажмите «с 1976»: тренд становится вчетверо круче. Почему так и какому числу верить, разберём ниже.

Идея

Модель — сумма столбцов с весами: прогноз для строки i равен β₀ + β₁x_i1 + … + β_d−1x_i,d−1. Веса выбирают так, чтобы сумма квадратов промахов по всем строкам была наименьшей. В этой точке промахи не связаны линейно ни с одним столбцом, и из этого получается система из d уравнений с d неизвестными — сколько бы строк ни было в данных.

Начнём с прямой y ≈ a + bx. Промах на точке — остаток e_i = y_i − a − bx_i. Остатки бывают разного знака, поэтому складывают их квадраты: такая сумма гладкая, у неё есть производные, и её минимум находится формулой. Это и есть метод наименьших квадратов. Приравняв нулю производные по a и b, получаем

Прямая проходит через центр масс облака (x̄, ȳ), а её наклон — это то, насколько согласованно x и y отклоняются от своих средних, делённое на разброс x. Вывод через проекции и виджет, где квадраты остатков нарисованы буквально, есть в «Царице наук», в опыте о лучшей прямой. У квадратов есть и обратная сторона, она там тоже видна: промах в десять единиц весит как сто промахов в одну, и одна грубая ошибка в данных разворачивает прямую.

Любое число столбцов

Сложим данные в таблицу X: n строк, d столбцов, первый столбец — единицы (он даёт свободный член β₀). Прогнозы всех строк — вектор Xβ, остатки — y − Xβ. Производная суммы квадратов остатков по весу β_j равна −2Σ_i x_ij(y_i − ŷ_i). В минимуме все d производных равны нулю, то есть вектор остатков перпендикулярен каждому столбцу:

Это нормальные уравнения; «нормальные» — от того самого перпендикуляра. Матрица X^TX квадратная, размера d × d: её элемент (j, k) — это сумма x_ijx_ik по всем строкам. Миллион строк или тридцать одна тысяча — система всё равно шесть на шесть.

Если столбцы зависимы — скажем, температура записана дважды, в градусах Цельсия и Фаренгейта, — минимумов бесконечно много: один столбец можно заменить другим. Метод Гаусса тогда наткнётся на нулевой ведущий элемент. NumPy в таком случае возвращает из всех решений самое короткое, но разумнее не подавать в модель одно и то же дважды.

«Линейная» — по коэффициентам

Температура за год — не прямая, а волна. Как её подогнать линейной регрессией? Любая синусоида с периодом в год раскладывается на синус и косинус: A sin(ωt + φ) = A cos φ · sin ωt + A sin φ · cos ωt. Амплитуда и сдвиг входят в формулу нелинейно, но коэффициенты при sin ωt и cos ωt — обычные числа, на которые умножаются столбцы. Поэтому в таблицу кладём два столбца, sin ωt и cos ωt с ω = 2π / 365,2425 дня, и метод наименьших квадратов сам найдёт и размах волны, и день, когда она выше всего. Ещё два столбца с двойной частотой — полугодовая волна — позволяют годовому ходу быть несимметричным: теплеть медленнее, чем остывать, или наоборот.

Слово «линейная» в названии относится к коэффициентам: прогноз — линейная комбинация столбцов. Сами столбцы могут быть любыми функциями данных — синусами, квадратами, логарифмами. Вот что дают модели из сцены — от одного столбца до шести:

columns                      trend °C/10y  wave °C      R²  miss °C
mean                                    -        -  0.0000    10.13
mean + trend                       +0.130        -  0.0010    10.12
mean + wave                             -    13.21  0.8509     3.91
mean + trend + wave                +0.121    13.21  0.8518     3.90
+ half-year wave                   +0.122    13.21  0.8553     3.85

Здесь R² — коэффициент детерминации: доля разброса вокруг среднего, которую объясняет модель, а промах — корень из среднего квадрата остатка. Годовая волна в Алматы — ±13,2 градуса вокруг среднего 8,76. С полугодовой поправкой самый холодный день модельного года — 8 января, самый тёплый — 24 июля, и между ними 26,6 градуса. Температура модели поднимается 197 дней, а опускается 168: остывает город быстрее, чем прогревается.

Почему тренд почти не заметил волну? Синус и косинус за целое число лет в сумме дают почти ноль, и столбец тренда с ними почти не связан: X^TX близка к блочной, и уравнения для тренда и для волны почти не зацепляют друг друга. Будь столбцы попарно перпендикулярны, каждый коэффициент считался бы отдельно, как будто других столбцов нет. Для зависимых столбцов это не так: добавишь один — меняются все остальные.

Код

Первый листинг — всё необходимое в общем виде: прямая по формуле со средними, метод Гаусса с выбором ведущего элемента и подгонка с любым числом столбцов через нормальные уравнения. Проверяем на четырёх наборах точек, о которых ещё пойдёт речь: формула и система дают одно и то же.

Второй листинг читает все 31 046 дней Алматы. Обратите внимание, что он хранит: не строки, а для каждого года матрицу X^TX, вектор X^Ty, сумму квадратов y и число дней. Этого хватает на всё остальное. Окно лет — сумма матриц своих лет. Модель с частью столбцов — подматрица: выбросить столбец из X значит выбросить строку и столбец из X^TX. Даже R² не требует второго прохода: в точке минимума сумма квадратов остатков равна Σy² − β^TX^Ty — раскройте |y − Xβ|² и подставьте нормальные уравнения.

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

Сколько это стоит

Каждая строка добавляет к X^TX все d² попарных произведений своих чисел (из-за симметрии хватило бы половины), к X^Ty — ещё d. Проход по данным стоит O(n·d²), решение системы методом Гаусса — O(d³), память — O(d²): строки можно читать из файла или из сети и сразу забывать. Для Алматы это около миллиона умножений и система шесть на шесть.

Линии на логарифмических осях почти параллельны, и наклон у них близок к единице: время растёт пропорционально числу строк. Расстояние между линиями — цена столбцов: при втрое большем числе столбцов умножений на строку почти в девять раз больше, а на замере время выросло примерно в семь раз — часть его уходит на сам цикл по строкам.

Раз суммы складываются, их можно хранить префиксными суммами: тогда наклон по любому окну считается за O(d²) действий, без единого прохода по данным. На этом построена вторая задача статьи.

Где формула перестаёт работать? Если столбцов сотни тысяч — слова в текстах, пиксели, признаки рекламных показов, — матрица X^TX не помещается в память, а d³ не дождаться. Тогда минимум ищут градиентным спуском: шаг за шагом, по случайным пачкам строк, не составляя систему вовсе. Так учат и нейросети; первая тренировка в главе курса о машинном обучении — та же прямая, найденная спуском.

Насколько верить наклону

Вот тот же тренд по разным окнам лет. Первое число — по дням, моделью со всеми шестью столбцами, как в сцене. Второе — прямая через 85 (или меньше) средних годовых и её стандартная ошибка:

1941-2025: +0.12 °C per 10 years; through the yearly means +0.12 ± 0.04
1941-1975: -0.72 °C per 10 years; through the yearly means -0.72 ± 0.14
1976-2025: +0.49 °C per 10 years; through the yearly means +0.49 ± 0.07
2000-2025: +0.98 °C per 10 years; through the yearly means +0.98 ± 0.17

Стандартная ошибка наклона — √(Σe² / (m − 2) / Σ(x − x̄)²) для m точек — говорит, насколько гулял бы наклон, если бы те же годы прожить заново с другой погодой. Формула предполагает, что остатки разных точек независимы. Для лет это близко к правде, для дней — нет: холодная неделя — это семь похожих друг на друга остатков, а не семь независимых. Посчитанная по дням, ошибка вышла бы во много раз меньше и обманывала бы. Поэтому её считают по годовым средним.

Четыре окна — четыре разных ответа, и каждый честный для своего окна. С 1976 года тренд +0,49 ± 0,07 градуса за десять лет, с 2000-го — около градуса. А за 1941–1975 — минус 0,72: в этих данных самый тёплый год всего ряда — 1941-й, 11,72 °C, теплее любого из последних. Точки у краёв окна тянут прямую сильнее всех: в формуле наклона каждая точка входит с множителем x_i − x̄, у крайней он наибольший. Так тёплые 1940-е опускают тренд за все 85 лет до +0,12.

Верны ли тёплые 1940-е — регрессия не знает, она проводит прямую через то, что ей дали. Данные здесь — реанализ ERA5: модель атмосферы, в которую подмешивают наблюдения. По данным ECMWF, в 1940 году их было около 17 тысяч в сутки на всю планету, в 2022 году — 25 миллионов, и ранние десятилетия реанализа поэтому надёжны меньше поздних. Отсюда правило: число «столько-то градусов за десятилетие» без окна и без источника данных мало что значит, и окно лучше выбирать до того, как посмотрели на ответ.

И последнее. Наведите лупу на годы за краем окна: там кривая модели пунктирная. Это уже прогноз — продолжение формулы туда, где её никто не проверял. Прямая, которая хорошо описывает 1976–2025 годы, ничего не обещает про 2100-й: линейная модель — описание данных в окне, а не закон природы.

Что прячет одно число

В 1973 году статистик Фрэнк Анскомб опубликовал четыре набора по одиннадцать точек. У всех четырёх одинаковые средние и разброс x и y, одна и та же прямая y = 3 + 0,5x и R² = 0,67 — первый листинг это и печатает. Посмотрите на них.

Первый набор — обычное облако вокруг прямой. Второй лежит на гладкой дуге, и прямая для него — неверная модель; переключите на остатки — они выстроятся аркой. Лечится это тем же приёмом, что и времена года: добавить столбец x², и R² станет 0,999999 (последняя строка листинга). В третьем наборе десять точек лежат точно на прямой, а одиннадцатая её перекосила: нажмите на неё, и прямая встанет на место. В четвёртом всю прямую задаёт одна точка: десять остальных имеют один и тот же x = 8. Уберите её — и наклон не определён вовсе: знаменатель формулы обратится в ноль.

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

История

Метод наименьших квадратов первым опубликовал Адриен Мари Лежандр в 1805 году, в приложении к книге о вычислении орбит комет. В 1809 году Карл Фридрих Гаусс изложил его подробнее, связал с законом распределения ошибок и заметил, что пользуется им с 1795 года; спор о первенстве так и не разрешился. Эту историю подробно рассказывает «Царица наук».

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

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

Задачи

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