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

Из песочницы Moneyball и Формула-1 модель прогнозирования результатов квалификаций

07.07.2020 18:11:11 |

Автор: admin

Сразу скажу: я не IT-специалист, а энтузиаст в сфере статистики. Помимо этого, я на протяжении многих лет участвовал в различных конкурсах прогнозов по Формуле-1. Отсюда вытекают и задачи, стоявшие перед моей моделью: выдавать прогнозы, которые были бы не хуже тех, которые создаются на глаз. А в идеале модель, конечно, должна обыгрывать человеческих оппонентов.

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

Для создания модели я свел в одну таблицу все результаты практик и квалификаций за сезоны 2018 и 2019. 2018-й год служил в качестве обучающей выборки, а 2019-й в качестве тестовой. По этим данным мы построили линейную регрессию. Если максимально просто объяснять регрессию, то наши данные это совокупность точек на координатной плоскости. Мы провели прямую, которая меньше всего отклоняется от совокупности этих точек. И функция, графиком которой является эта прямая это и есть наша линейная регрессия.

От известной из школьной программы формулы $inline$y = kx + b $inline$ нашу функцию отличает только то, что переменных у нас две. Первая переменная (X1) это отставание в третьей практике, а вторая переменная (X2) среднее отставание по предыдущим квалификациям. Эти переменные не равнозначны, и одна из наших целей определить вес каждой переменной в диапазоне от 0 до 1. Чем дальше переменная от нуля, тем большее значение она имеет при объяснении зависимой переменной. В нашем случае в качестве зависимой переменной выступает время на круге, выраженное в отставании от лидера (или точнее, от некоего идеального круга, поскольку у всех пилотов эта величина была положительной).

Поклонники книги Moneyball (в фильме этот момент не объясняется) могут вспомнить, что там с помощью линейной регрессии определили, что процент занятия базы, aka OBP (on-base percentage), более тесно связан с заработанными ранами, чем другие статистические показатели. Мы преследуем примерно такую же цель: понять, какие именно факторы наиболее тесно связаны с результатами квалификаций. Один из больших плюсов регрессии в том, что она не требует продвинутого знания математики: мы просто задаем данные, а потом Excel или другой табличный редактор выдает нам готовые коэффициенты.

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

Тут надо отметить важный момент. Конечный результат складывался из двух независимых параметров, каждый из которых вытекал из двух независимых регрессий. Первый параметр сила команды на этом этапе, точнее, отставание лучшего пилота команды от лидера. Второй параметр распределение сил внутри команды.

Что это значит на примере? Допустим, мы берем Гран-при Венгрии сезона-2019. Модель показывает, что отставание Феррари от лидера составит 0,218 секунды. Но это отставание первого пилота, а кто им будет Феттель или Леклер и какой разрыв между ними будет, определяется другим параметром. В этом примере модель показала, что впереди будет Феттель, а Леклер проиграет ему 0,096 секунды.

К чему такие сложности? Не проще ли рассматривать каждого пилота по отдельности вместо этой разбивки на отставание команды и отставание первого пилота от второго внутри команды? Возможно, это так, но мои личные наблюдения показывают, что смотреть на результаты команды гораздо надежнее, чем на результаты каждого пилота. Один пилот может допустить ошибку, или вылететь с трассы, или у него будут технические проблемы все это будет вносить хаос в работу модели, если только не отслеживать вручную каждую форс-мажорную ситуацию, что требует слишком много времени. Влияние форс-мажоров на результаты команды гораздо меньше.

Но вернемся к моменту, где мы хотели оценить, насколько хорошо выбранные нами независимые переменные объясняют изменения функции. Это можно сделать с помощью коэффициента детерминации. Он продемонстрирует, в какой степени результаты квалификации объясняются результатами практик и предыдущих квалификаций.

Поскольку мы строили две регрессии, то и коэффициента детерминации у нас тоже два. Первая регрессия отвечает за уровень команды на этапе, вторая за противостояние между пилотами одной команды. В первом случае коэффициент детерминации равен 0,82, то есть 82% результатов квалификаций объясняются выбранными нами факторами, а еще 18% какими-то другими факторами, которые мы не учли. Это достаточно неплохой результат. Во втором случае коэффициент детерминации составил 0,13.

Эти показатели, по сути, означают, что модель достаточно хорошо предсказывает уровень команды, но испытывает проблемы с определением разрыва между партнерами по команде. Однако для итоговой цели нам не нужно знать разрыв, нам достаточно знать, кто из двух пилотов будет выше, и с этим модель в основном справляется. В 62% случаев модель ставила выше того пилота, который действительно был выше по итогам квалификации.

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

Итоговая формула выглядит так:

Первый пилот:

$$display$$Y1 = (0,618 * X1 + 0,445 * X2)$$display$$

Второй пилот:

$$display$$Y2 = Y1 + (0,313 * X1 + 0,511 * X2)$$display$$

Напоминаю, что X1 это отставание в третьей практике, а X2 среднее отставание по предыдущим квалификациям.

Что означают эти цифры. Они означают, что уровень команды в квалификации на 60% определяется результатами третьей практики и на 40% результатами квалификаций на предыдущих этапах. Соответственно, результаты третьей практики в полтора раза более значимый фактор, чем результаты предыдущих квалификаций.

Поклонники Формулы-1 наверняка знают ответ на этот вопрос, но для остальных следует прокомментировать, почему я брал результаты именно третьей практики. В Формуле-1 проводится три практики. Однако именно в последней из них команды традиционно тренируют квалификацию. Однако же в тех случаях, когда третья практика срывается из-за дождя или других форс-мажоров, я брал результаты второй практики. Насколько я помню, в 2019 году был только один такой случай на Гран-при Японии, когда из-за тайфуна этап прошел в укороченном формате.

Также кто-то наверняка заметил, что модель использует среднее отставание в предыдущих квалификациях. Но как быть на первом этапе сезона? Я использовал отставания с предыдущего года, но не оставлял их как есть, а вручную их корректировал, основываясь на здравом смысле. Например, в 2019 году Феррари в среднем была быстрее, чем Ред Булл на 0,3 секунды. Однако похоже, что у итальянской команды не будет такого преимущества в этом году, а может, они и вовсе будут позади. Поэтому для первого этапа сезона 2020, Гран-при Австрии, я вручную приблизил Ред Булл к Феррари.

Таким образом я получал отставание каждого пилота, ранжировал пилотов по отставанию и получил итоговый прогноз на квалификацию. При этом важно понимать, что первый и второй пилот это чистые условности. Возвращаясь к примеру с Феттелем и Леклером, на Гран-при Венгрии модель посчитала первым пилотом Себастьяна, но на многих других этапах она отдавала предпочтение Леклеру.

Результаты

Как я говорил, задачей было создать такую модель, которая позволит прогнозировать не хуже людей. За основу я брал свои прогнозы и прогнозы своих партнеров по команде, которые создавались на глаз, но при внимательном изучении результатов практик и совместном обсуждении.

Система оценки была следующая. Учитывалась только первая десятка пилотов. За точное попадание прогноз получал 9 баллов, за промах в 1 позицию 6 баллов, за промах в 2 позиции 4 балла, за промах в 3 позиции 2 балла и за промах в 4 позиции 1 балл. То есть если в прогнозе пилот стоит на 3-м месте, а в результате он взял поул-позишн, то прогноз получал 4 балла.

При такой системе максимальное количество баллов за 21 Гран-при 1890.
Человеческие участники набрали 1056, 1048 и 1034 балла соответственно.
Модель набрала 1031 балл, хотя при легкой манипуляции с коэффициентами я также получал 1045 и 1053 балла.

Лично я доволен итогами, поскольку это мой первый опыт в построении регрессий, и он привел к достаточно приемлемым результатам. Безусловно, хотелось бы их улучшить, поскольку я уверен, что с помощью построения моделей, даже таких простых, как эта, можно добиваться лучшего результата, чем просто оценивая данные на глаз. В рамках этой модели можно было бы, например, учесть тот фактор, что некоторые команды слабы в практиках, но выстреливают в квалификациях. Например, есть наблюдение, что Мерседес часто не был лучшей командой в ходе тренировок, но гораздо лучше выступал в квалификациях. Однако эти человеческие наблюдения не нашли отражения в модели. Поэтому в сезоне 2020, который начнется в июле (если не произойдет ничего непредвиденного), я хочу проверить эту модель в соревновании против живых прогнозистов, а также найти, как ее можно сделать лучше.

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

Подробнее..

Категории: Алгоритмы , Big data , Data mining , Моделирование , Статистика , Статистический анализ , Статистическое моделирование , Регрессия , Линейная регрессия , Спортивная статистика , Формула-1 , Формула 1

Простыми словами о простых линейных функциях

06.06.2021 14:21:06 |

Автор: admin

Случайный лес (в буквальном смысле, сфотографировал с телефона)

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

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

Для начала есть смысл взглянуть на сами данные. Если вы отобразите массив данных с несколькими миллионами объектов, то человек не сможет всё это осмысленно прочитать. Придётся придумывать способ понятно описать его. Частая ошибка не изучив набор данных использовать среднее значение или медиану. Подобный подход чреват серьёзным искажением реальной картины.

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

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

В итоге мы видим распределение случайной величины. Регулируя число интервалов можно добиться адекватного представления. Теперь нам нужен способ увидеть взаимосвязь этой случайной величины с другой. Отобразим наблюдения в виде точек на плоскости, где по оси X будет score, а по Y наш единственный предиктор (признак). И вот перед вами появляется график:

На нём показано, что обе переменных сильно коррелируют между собой, следовательно, делаем предположение о линейной зависимости (score=a+bx). А теперь самая суть, но очень кратко: простая линейная регрессия это произведение коэффициента и признака, к которым добавляется смещение. Вот и всё. Степень линейной взаимосвязи вычисляют следующим образом:

Не все точки выстроились на одну линию, что говорит нам о небольшом влиянии неизвестного фактора. Это я специально сделал, для красоты. Подобрать коэффициенты можно с помощью готовых решений, допустим, scikit-learn. Так, например, класс LinearRegression после обучения (fit) позволяет получить смещение (intercept) и массив коэффициентов (coef). Подставляем значения в формулу и проверяем результат с помощью метрик mean_absolute_error и median_absolute_error. Собственно, это и есть решение нашей задачи. В случае множества признаков суть не меняется: под капотом всё устроено аналогичным образом (смещение + скалярное произведение вектора признаков и соответствующих коэффициентов).

А теперь превратим эту функцию в классификатор, который называется логистическая регрессия. Под капотом это обычная линейная регрессия, только результат её работы передают в сигмоиду. Далее выполняется проверка если результат больше 0.5, то класс 1, иначе класс 0. По сути, это формула гиперплоскости, разделяющей классы.

Обе упомянутые линейные модели это очень быстрые и простые алгоритмы, которые достаточно легко внедрить в любой проект. Надеюсь, теперь стало чуть проще понять, как это устроено под капотом. Сталкиваясь с этими алгоритмами, начинающему разработчику будет проще понять их принципы работы.

Подробнее..

Категории: Python , Математика , Машинное обучение , Data mining , Data engineering , Линейная регрессия , Логистическая регрессия

Основы линейной регрессии

13.08.2020 00:05:01 |

Автор: admin

Здравствуй, Хабр!

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

! Осторожно, трафик! В статье присутствует заметное число изображений для иллюстраций, часть в формате gif.

Содержание

Введение
Метод наименьших квадратов
Мультилинейная регрессия
- Линейная алгебра
Произвольный базис
Заключительные замечания
- Проблема выбора размерности
- Численные методы
Реклама и заключение

Введение

Есть три сходных между собой понятия, три сестры: интерполяция, аппроксимация и регрессия.
У них общая цель: из семейства функций выбрать ту, которая обладает определенным свойством.

Интерполяция способ, как из семейства функций выбрать ту, которая проходит через заданные точки. Затем ее обычно используют для вычисления в точках, отличных от заданных. Например, мы вручную задаем цвет нескольким точкам и хотим чтобы цвета остальных точек образовали плавные переходы между заданными. Или задаем ключевые кадры анимации и хотим плавные переходы между ними. Классические примеры: интерполяция полиномами Лагранжа, сплайн-интерполяция.

$display$

Аппроксимация способ, как из семейства простых функций выбрать приближение для сложной функции на отрезке, при этом ошибка не должна превышать определенного предела. Аппроксимацию используют, когда нужно получить функцию, похожую на данную, но более удобную для вычислений и манипуляций (дифференцирования, интегрирования и т.п). При оптимизации критических участков кода часто используют аппроксимацию: если значение функции вычисляется много раз в секунду и не нужна абсолютная точность, то можно обойтись более простым аппроксимантом с меньшей ценой вычисления. Классические примеры включают ряд Тейлора на отрезке, аппроксимацию функции ортогональными многочленами, аппроксимацию Паде, аппроксимацию синуса Бхаскара и т.п.

Регрессия способ, как из семейства функций выбрать ту, которая минимизирует функцию потерь. Последняя характеризует насколько сильно пробная функция отклоняется от значений в заданных точках. Если точки получены в эксперименте, они неизбежно содержат ошибку измерений, шум, поэтому разумнее требовать, чтобы функция передавала общую тенденцию, а не точно проходила через все точки. В каком-то смысле регрессия это интерполирующая аппроксимация: мы хотим провести кривую как можно ближе к точкам и при этом сохранить ее максимально простой чтобы уловить общую тенденцию. За баланс между этими противоречивыми желаниями как-раз отвечает функция потерь (в английской литературе loss function или cost function).

В этой статье мы рассмотрим линейную регрессию. Это означает, что семейство функций, из которых мы выбираем, представляет собой линейную комбинацию наперед заданных базисных функций ${f_i}$

$f = \sum_i w_i f_i.$

Цель регрессии найти коэффициенты этой линейной комбинации, и тем самым определить регрессионную функцию $inline$ (которую также называют моделью). Отмечу, что линейную регрессию называют линейной именно из-за линейной комбинации базисных функций это не связано с самыми базисными функциями (они могут быть линейными или нет).

Регрессия с нами уже давно: впервые метод опубликовал Лежандр в 1805 году, хотя Гаусс пришел к нему раньше и успешно использовал для предсказания орбиты кометы (на самом деле карликовой планеты) Цереры. Существует множество вариантов и обобщений линейной регрессии: LAD, метод наименьших квадратов, Ridge регрессия, Lasso регрессия, ElasticNet и многие другие.

гифка

Точки генерируются случайно по распределению Гаусса с заданным средним и вариациями. Синяя линия регрессионная прямая.

Можно поиграться с демонстрацией в GoogleColab.
Много других материалов по классическому машинному обучению на соответствующей страничке на GitHub

Метод наименьших квадратов

Начнём с простейшего двумерного случая. Пусть нам даны точки на плоскости $\{(x_1,y_1),\cdots,(x_N,y_N)\}$ и мы ищем такую аффинную функцию

$f(x) = a + b \cdot x,$

чтобы ее график ближе всего находился к точкам. Таким образом, наш базис состоит из константной функции и линейной $inline$ .

Как видно из иллюстрации, расстояние от точки до прямой можно понимать по-разному, например геометрически это длина перпендикуляра. Однако в контексте нашей задачи нам нужно функциональное расстояние, а не геометрическое. Нас интересует разница между экспериментальным значением и предсказанием модели для каждого $inline$ поэтому измерять нужно вдоль оси $inline$ .

Первое, что приходит в голову, в качестве функции потерь попробовать выражение, зависящее от абсолютных значений разниц $inline$ . Простейший вариант сумма модулей отклонений $\sum_i |f(x_i) - y_i|$ приводит к Least Absolute Distance (LAD) регрессии.

Впрочем, более популярная функция потерь сумма квадратов отклонений регрессанта от модели. В англоязычной литературе она носит название Sum of Squared Errors (SSE)

$\text{SSE}(a,b)=\text{SS}_{res[iduals]}=\sum_{i=1}^N{\text{отклонение}_i}^2=\sum_{i=1}^N(y_i-f(x_i))^2=\sum_{i=1}^N(y_i-a-b\cdot x_i)^2,$

Метод наименьших квадратов (по англ. OLS) линейная регрессия c $\text{SSE}(a,b)$ в качестве функции потерь.

Такой выбор прежде всего удобен: производная квадратичной функции линейная функция, а линейные уравнения легко решаются. Впрочем, далее я укажу и другие соображения в пользу $\text{SSE}(a,b)$ .

гифка

Регрессионная прямая (синяя) и пробная прямая (зеленая). Справа показана функция потерь и точки соответствующие параметра пробной и регрессионной прямых.

Математический анализ

Простейший способ найти $\text{argmin}_{a,b} \, \text{SSE}(a,b)$ вычислить частные производные по $inline$ и $inline$ , приравнять их нулю и решить систему линейных уравнений

$\begin{aligned} \frac{\partial}{\partial a}\text{SSE}(a,b)&=-2\sum_{i=1}^N(y_i-a-bx_i), \\ \frac{\partial}{\partial b}\text{SSE}(a,b)&=-2\sum_{i=1}^N(y_i-a-bx_i)x_i. \end{aligned}$

Значения параметров, минимизирующие функцию потерь, удовлетворяют уравнениям

$\begin{aligned} 0 &= -2\sum_{i=1}^N(y_i-\hat{a}-\hat{b}x_i), \\ 0 &= -2\sum_{i=1}^N(y_i-\hat{a}-\hat{b}x_i)x_i, \end{aligned}$

которые легко решить

$\begin{aligned} \hat{a}&=\frac{\sum_i y_i}{N}-\hat{b}\frac{\sum_i x_i}{N},\\ \hat{b}&=\frac{\frac{\sum_i x_i y_i}{N}-\frac{\sum_i x_i\sum_i y_i}{N^2}}{\frac{\sum_i x_i^2}{N}-\left(\frac{\sum_i x_i^2}{N}\right)^2}. \end{aligned}$

Мы получили громоздкие и неструктурированные выражения. Сейчас мы их облагородим и вдохнем в них смысл.

Статистика

Полученные формулы можно компактно записать с помощью статистических эстиматоров: среднего $\langle{\cdot}\rangle$ , вариации $\sigma_{\cdot}$ (стандартного отклонения), ковариации $\sigma({\cdot},{\cdot})$ и корреляции $\rho({\cdot},{\cdot})$

$\begin{aligned} \hat{a}&=\langle{y}\rangle-\hat{b}\langle{x}\rangle, \\ \hat{b}&=\frac{\langle{xy}\rangle-\langle{x}\rangle\langle{y}\rangle}{\langle{x^2}\rangle-\langle{x}\rangle^2}. \end{aligned}$

Перепишем $\hat{b}$ как

$\hat{b} = \frac{\sigma(x,y)}{\sigma_x^2},$

где $\sigma_x$ это нескорректированное (смещенное) стандартное выборочное отклонение, а $\sigma(x,y)$ ковариация. Теперь вспомним, что коэффициент корреляции (коэффициент корреляции Пирсона)

$\rho(x,y)=\frac{\sigma(x,y)}{\sigma_x \sigma_y}$

и запишем

$\hat{b}=\rho(x,y)\frac{\sigma_y}{\sigma_x}.$

Теперь мы можем оценить все изящество дескриптивной статистики и записать уравнение регрессионной прямой так

$\boxed{y-\langle {y} \rangle = \rho(x,y)\frac{\sigma_y}{\sigma_x}(x-\langle {x} \rangle)}.$

Во-первых, это уравнение сразу указывает на два свойства регрессионной прямой:

прямая проходит через центр масс $(\langle{x}\rangle, \langle{y}\rangle)$ ;
если по оси $inline$ за единицу длины выбрать $\sigma_x$ , а по оси $inline$ $\sigma_y$ , то угол наклона прямой будет от $-45^\circ$ до $45^\circ$ . Это связано с тем, что $-1 \leq\rho(x,y)\leq 1$ .

Во-вторых, теперь становится понятно, почему метод регрессии называется именно так. В единицах стандартного отклонения $inline$ отклоняется от своего среднего значения меньше чем $inline$ , потому что $|\rho(x,y)|\leq1$ . Это называется регрессией(от лат. regressus возвращение) по отношению к среднему. Это явление было описано сэром Фрэнсисом Гальтоном в конце XIX века в его статье Регрессия к посредственности при наследовании роста. В статье показано, что черты (такие как рост), сильно отклоняющиеся от средних, редко передаются по наследству. Характеристики потомства как-бы стремятся к среднему на детях гениев природа отдыхает.

Возведя коэффициент корреляции в квадрат, получим коэффициент детерминации $R = \rho^2$ . Квадрат этой статистической меры показывает насколько хорошо регрессионная модель описывает данные. $inline$ , равный $inline$ , означает что функция идеально ложится на все точки данные идеально скоррелированны. Можно доказать, что $inline$ показывает какая доля вариативности в данных объясняется лучшей из линейных моделей. Чтобы понять, что это значит, введем определения

$\begin{aligned} \text{Var}_{data} &= \frac{1}{N}\sum_i (y_i-\langle y \rangle)^2, \\ \text{Var}_{res} &= \frac{1}{N} \sum_i (y_i-\text{модель}(x_i))^2, \\ \text{Var}_{reg} &= \frac{1}{N} \sum_i (\text{модель}(x_i)-\langle y \rangle)^2. \end{aligned}$

$\text{Var}_{data}$ вариация исходных данных (вариация точек $inline$ ).

$\text{Var}_{res}$ вариация остатков, то есть вариация отклонений от регрессионной модели от $inline$ нужно отнять предсказание модели и найти вариацию.

$\text{Var}_{reg}$ вариация регрессии, то есть вариация предсказаний регрессионной модели в точках $inline$ (обратите внимание, что среднее предсказаний модели совпадает с $\langle y \rangle$ ).

Дело в том, что вариация исходных данных разлагается в сумму двух других вариаций: вариации, которая объясняется моделью, и вариации случайного шума (остатков)

$\boxed{{\color{red}{\text{Var}_{data}}} ={\color{green}{\text{Var}_{res}}}+ {\color{blue}{\text{Var}_{reg}}}.}$

или

$\sigma^2_{data} =\sigma^2_{res}+ \sigma^2_{reg}.$

Как видим, стандартные отклонения образуют прямоугольный треугольник.

Мы стремимся избавиться от вариативности связанной с шумом и оставить лишь вариативность, которая объясняется моделью, хотим отделить зерна от плевел. О том, насколько это удалось лучшей из линейных моделей, свидетельствует $inline$ , равный единице минус доля вариации ошибок в суммарной вариации

$R^2=\frac{\text{Var}_{data}-\text{Var}_{res}}{\text{Var}_{data}}=1-\frac{\color{green}{\text{Var}_{res}}}{\color{red}{\text{Var}_{data}}}$

или доле объясненной вариации (доля вариации регрессии в полной вариации)

$R^2=\frac{\color{blue}{\text{Var}_{reg}}}{\color{red}{\text{Var}_{data}}}.$

$inline$ равен косинусу угла в прямоугольном треугольнике $(\sigma_{data}, \sigma_{reg}, \sigma_{res})$ . Кстати, иногда вводят долю необъясненной вариации $inline$ и она равна квадрату синуса в этом треугольнике. Если коэффициент детерминации мал, возможно мы выбрали неудачные базисные функции, линейная регрессия неприменима вовсе и т.п.

Теория вероятностей

Ранее мы пришли к функции потерь $\text{SSE}(a,b)$ из соображений удобства, но к ней же можно прийти с помощью теории вероятностей и метода максимального правдоподобия (ММП). Напомню вкратце его суть. Предположим, у нас есть $inline$ независимых одинаково распределенных случайных величин (в нашем случае результатов измерений). Мы знаем вид функции распределения (напр. нормальное распределение), но хотим определить параметры, которые в нее входят (например $\mu$ и $\sigma$ ). Для этого нужно вычислить вероятность получить $inline$ датапоинтов в предположении постоянных, но пока неизвестных параметров. Благодаря независимости измерений, мы получим произведение вероятностей реализации каждого измерения. Если мыслить полученную величину как функцию параметров (функция правдоподобия) и найти её максимум, мы получим оценку параметров. Зачастую вместо функции правдоподобия используют ее логарифм дифференцировать его проще, а результат тот же.

Вернемся к задаче простой регрессии. Допустим, что значения $inline$ нам известны точно, а в измерении $inline$ присутствует случайный шум (свойство слабой экзогенности). Более того, положим, что все отклонения от прямой (свойство линейности) вызваны шумом с постоянным распределением (постоянство распределения). Тогда

$y = a + bx + \epsilon,$

где $\epsilon$ нормально распределенная случайная величина

$\epsilon \sim \mathcal{N}(0,\,\sigma^{2}), \qquad p(\epsilon) = \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-\frac{\epsilon^2}{2\sigma^2}}.$

Исходя из предположений выше, запишем функцию правдоподобия

$\begin{aligned} L(a,b|\mathbf{y})&=P(\mathbf{y}|a,b)=\prod_i P(y_i|a,b)=\prod_i p(y_i-a-bx|a,b)=\\ &= \prod_i \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-\frac{(y_i-a-bx)^2}{2\sigma^2}}= \frac{1}{\sqrt{2 \pi \sigma^2}}e^{-\frac{\sum_i (y_i-a-bx)^2}{2 \sigma^2}}=\\ &= \frac{1}{\sqrt{2 \pi \sigma^2}}e^{-\frac{\text{SSE}(a,b)}{2 \sigma^2}} \end{aligned}$

и ее логарифм

$l(a,b|\mathbf{y})=\log{L(a,b|\mathbf{y})}=-\text{SSE}(a,b)+const.$

Таким образом, максимум правдоподобия достигается при минимуме $\text{SSE}$

$(\hat{a},\hat{b})=\text{argmax}_{a,b} \, l(a,b|\mathbf{y}) = \text{argmin}_{a,b} \, \text{SSE}(a,b),$

что дает основание принять ее в качестве функции потерь. Кстати, если

$\begin{aligned} \epsilon \sim \text{Laplace}(0, \alpha), \qquad p_{L}(\epsilon; \mu, \alpha) =\frac{\alpha}{2}e^{-\alpha |\epsilon-\mu|} \end{aligned}$

мы получим функцию потерь LAD регрессии

$E_{LAD}(a,b)=\sum_i |y_i-a-bx_i|,$

которую мы упоминали ранее.

Подход, который мы использовали в этом разделе один из возможных. Можно прийти к такому же результату, используя более общие свойства. В частности, свойство постоянства распределения можно ослабить, заменив на свойства независимости, постоянства вариации (гомоскедастичность) и отсутствия мультиколлинеарности. Также вместо ММП эстимации можно воспользоваться другими методами, например линейной MMSE эстимацией.

Мультилинейная регрессия

До сих пор мы рассматривали задачу регрессии для одного скалярного признака $inline$ , однако обычно регрессор это $inline$ -мерный вектор $\mathbf{x}$ . Другими словами, для каждого измерения мы регистрируем $inline$ фич, объединяя их в вектор. В этом случае логично принять модель с $inline$ независимыми базисными функциями векторного аргумента $inline$ степеней свободы соответствуют $inline$ фичам и еще одна регрессанту $inline$ . Простейший выбор линейные базисные функции $(1, x_1, \cdots, x_n)$ . При $inline$ получим уже знакомый нам базис $inline$ .

Итак, мы хотим найти такой вектор (набор коэффициентов) $\mathbf{w}$ , что

$\sum_{j=0}^n w_j x_j^{(i)}= \mathbf{w}^{\top}\mathbf{x}^{(i)} \simeq y_i, \qquad \qquad \qquad \qquad i=1\dots N.$

Знак " $\simeq$ " означает, что мы ищем решение, которое минимизирует сумму квадратов ошибок

$\hat{\mathbf{w}}=\text{argmin}_\mathbf{w} \, \sum_{i=1}^N \left({y_i - \mathbf{w}^{\top}\mathbf{x}^{(i)}}\right)^2$

Последнее уравнение можно переписать более удобным образом. Для этого расположим $\mathbf{x}^{(i)}$ в строках матрицы (матрицы информации)

$X= \begin{pmatrix} - & \mathbf{x}^{(1)\top} & - \\ \cdots & \cdots & \cdots\\ - & \mathbf{x}^{(N)\top} & - \end{pmatrix} = \begin{pmatrix} | & | & & | \\ \mathbf{x}_0 & \mathbf{x}_1 & \cdots & \mathbf{x}_n \\ | & | & & | \end{pmatrix} = \begin{pmatrix} 1 & x^{(1)}_{1} & \cdots & x^{(1)}_{n} \\ \cdots & \cdots & \cdots & \cdots\\ 1 & x^{(N)}_{1} & \cdots & x^{(N)}_{n} \end{pmatrix}.$

Тогда столбцы матрицы $\mathbf{x}_{i}$ отвечают измерениям $inline$ -ой фичи. Здесь важно не запутаться: $inline$ количество измерений, $inline$ количество признаков (фич), которые мы регистрируем. Систему можно записать как

$X \, \mathbf{w} \simeq \mathbf{y}.$

Квадрат нормы разности векторов в правой и левой частях уравнения образует функцию потерь

$\text{SSE}(\mathbf{w}) = {\|\mathbf{y}-X \mathbf{w}\|}^2, \qquad \qquad \mathbf{w} \in \mathbb{R}^{n+1}; \, \mathbf{y} \in \mathbb{R}^{N},$

которую мы намерены минимизировать

$\begin{aligned} \hat{\mathbf{w}}&=\text{argmin}_\mathbf{w} \, \text{SSE}(\mathbf{w}) = \text{argmin}_\mathbf{w} \, (\mathbf{y}-X \mathbf{w})^{\top}(\mathbf{y}-X \mathbf{w})=\\ &= \text{argmin}_\mathbf{w} \,(\mathbf{y}^{\top}\mathbf{y}-2\mathbf{w}^{\top}X^{\top}\mathbf{y}+\mathbf{w}^{\top}X^{\top}X\mathbf{w}). \end{aligned}$

Продифференцируем финальное выражение по $\mathbf{w}$ (если забыли как это делается загляните в Matrix cookbook)

$\frac{\partial \, \text{SSE}(\mathbf{w})}{\partial \mathbf{w}}=-2 X^{\top}\mathbf{y}+2 X^{\top}X\mathbf{w},$

приравняем производную к $\mathbf{0}$ и получим т.н. нормальные уравнения

$X^{\top}X \, \hat{\mathbf{w}}=X^{\top}\mathbf{y}.$

Если столбцы матрицы информации $inline$ линейно независимы (нет идеально скоррелированных фич), то матрица $X^{\top}X$ имеет обратную (доказательство можно посмотреть, например, в видео академии Хана). Тогда можно записать

$\boxed{\hat{\mathbf{w}} = (X^{\top}X)^{-1}X^{\top}\mathbf{y}=X^{+}\mathbf{y}},$

где

$X^{+}=(X^{\top}X)^{-1}X^{\top}$

псевдообратная к $inline$ . Понятие псевдообратной матрицы введено в 1903 году Фредгольмом, она сыграла важную роль в работах Мура и Пенроуза.

Напомню, что обратить $X^{\top}X$ и найти $X^{+}$ можно только если столбцы $inline$ линейно независимы. Впрочем, если столбцы $inline$ близки к линейной зависимости, вычисление $(X^{\top}X)^{-1}$ уже становится численно нестабильным. Степень линейной зависимости признаков в $inline$ или, как говорят, мультиколлинеарности матрицы $X^{\top}X$ , можно измерить числом обусловленности отношением максимального собственного значения к минимальному. Чем оно больше, тем ближе $X^{\top}X$ к вырожденной и неустойчивее вычисление псевдообратной.

Линейная алгебра

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

$X \, \mathbf{w} \simeq \mathbf{y}.$

Если количество переменных равно количеству неизвестных и уравнения линейно независимы, то система имеет единственное решение. Однако, если число измерений превосходит число признаков, то есть уравнений больше чем неизвестных система становится несовместной, переопределенной. В этом случае лучшее, что мы можем сделать выбрать вектор $\mathbf{w}$ , образ которого $X\mathbf{w}$ ближе остальных к $\mathbf{y}$ . Напомню, что множество образов или колоночное пространство $\mathcal{C}(X)$ это линейная комбинация вектор-столбцов матрицы $inline$

$\begin{pmatrix} | & | & & | \\ \mathbf{x}_0 & \mathbf{x}_1 & \cdots & \mathbf{x}_n \\ | & | & & | \end{pmatrix} \mathbf{w} = w_0 \mathbf{x}_0 + w_1 \mathbf{x}_1 + \cdots w_n \mathbf{x}_n .$

$\mathcal{C}(X)$ $inline$ -мерное линейное подпространство (мы считаем фичи линейно независимыми), линейная оболочка вектор-столбцов $inline$ . Итак, если $\mathbf{y}$ принадлежит $\mathcal{C}(X)$ , то мы можем найти решение, если нет будем искать, так сказать, лучшее из нерешений.

Если в дополнение к векторам $\mathcal{C}(X)$ мы рассмотрим все вектора им перпендикулярные, то получим еще одно подпространство и сможем любой вектор из $\mathbb{R}^{N}$ разложить на две компоненты, каждая из которых живет в своем подпространстве. Второе, перпендикулярное пространство, можно характеризовать следующим образом (нам это понадобится в дальнейшем). Пускай $\mathbf{v} \in \mathbb{R}^{N}$ , тогда

$X^\top \mathbf{v} = \begin{pmatrix} - & \mathbf{x}_0^{\top} & - \\ \cdots & \cdots & \cdots\\ - & \mathbf{x}_n^{\top} & - \end{pmatrix} \mathbf{v} = \begin{pmatrix} \mathbf{x}_0^{\top} \cdot \mathbf{v} \\ \cdots \\ \mathbf{x}_n^{\top} \cdot \mathbf{v} \\ \end{pmatrix}$

равен нулю в том и только в том случае, если $\mathbf{v}$ перпендикулярен всем $\mathbf{x}_i$ , а значит и целому $\mathcal{C}(X)$ . Таким образом, мы нашли два перпендикулярных линейных подпространства, линейные комбинации векторов из которых полностью, без дыр, покрывают все $\mathbb{R}^N$ . Иногда это обозначают c помощью символа ортогональной прямой суммы

где $\text{ker}(X^{\top})=\{\mathbf{v}|X^{\top}\mathbf{v}=\mathbf{0}\}$ . В каждое из подпространств можно попасть с помощью соответствующего оператора проекции, но об этом ниже.

Теперь представим $\mathbf{y}$ в виде разложения

$\mathbf{y} = \mathbf{y}_{\text{proj}} + \mathbf{y}_{\perp}, \qquad \mathbf{y}_{\text{proj}} \in \mathcal{C}(X), \qquad \mathbf{y}_{\perp} \in \text{ker}(X^{\top}).$

Если мы ищем решение $\hat{\mathbf{w}}$ , то естественно потребовать, чтобы $|| \mathbf{y} - X\mathbf{w} ||$ была минимальна, ведь это длина вектора-остатка. Учитывая перпендикулярность подпространств и теорему Пифагора

$\text{argmin}_\mathbf{w} || \mathbf{y} - X\mathbf{w} || = \text{argmin}_\mathbf{w} || \mathbf{y}_{\perp} + \mathbf{y}_{\text{proj}} - X\mathbf{w} || = \text{argmin}_\mathbf{w} \sqrt{|| \mathbf{y}_{\perp} ||^2 + || \mathbf{y}_{\text{proj}} - X\mathbf{w} ||^2},$

но поскольку, выбрав подходящий $\mathbf{w}$ , я могу получить любой вектор колоночного пространства, то задача сводится к

$X\hat{\mathbf{w}} = \mathbf{y}_{\text{proj}},$

а $\mathbf{y}_{\perp}$ останется в качестве неустранимой ошибки. Любой другой выбор $\hat{\mathbf{w}}$ сделает ошибку только больше.

Если теперь вспомнить, что $X^{\top} \mathbf{y}_{\perp} = \mathbf{0}$ , то легко видеть

$X^\top X \mathbf{w} = X^{\top} \mathbf{y}_{\text{proj}} = X^{\top} \mathbf{y}_{\text{proj}} + X^{\top} \mathbf{y}_{\perp} = X^{\top} \mathbf{y},$

что очень удобно, так как $\mathbf{y}_{\text{proj}}$ у нас нет, а вот $\mathbf{y}$ есть. Вспомним из предыдущего параграфа, что $X^{\top} X$ имеет обратную при условии линейной независимости признаков и запишем решение

$\mathbf{w} = (X^\top X)^{-1} X^\top \mathbf{y} = X^{+} \mathbf{y},$

где $X^{+}$ уже знакомая нам псевдообратная матрица. Если нам интересна проекция $\mathbf{y}_{\text{proj}}$ , то можно записать

$\mathbf{y}_{\text{proj}} = X \mathbf{w} = X X^{+} \mathbf{y} = \text{Proj}_X \mathbf{y},$

где $\text{Proj}_X$ оператор проекции на колоночное пространство.

Выясним геометрический смысл коэффициента детерминации.

Заметьте, что фиолетовый вектор $\bar{y} \cdot \boldsymbol{1}=\bar{y} \cdot (1,1,\dots,1)^{\top}$ пропорционален первому столбцу матрицы информации $inline$ , который состоит из одних единиц согласно нашему выбору базисных функций. В RGB треугольнике

${\color{red}{\mathbf{y}-\hat{y} \cdot \boldsymbol{1}}}={\color{green}{\mathbf{y}-\bar{\mathbf{y}}}}+{\color{blue}{\hat{\mathbf{y}}-\bar{y} \cdot \boldsymbol{1}}}.$

Так как этот треугольник прямоугольный, то по теореме Пифагора

${\color{red}{\|\mathbf{y}-\hat{y} \cdot \boldsymbol{1}\|^2}}={\color{green}{\|\mathbf{y}-\bar{\mathbf{y}}\|^2}}+{\color{blue}{\|\hat{\mathbf{y}}-\bar{y} \cdot \boldsymbol{1}\|^2}}.$

Это геометрическая интерпретация уже известного нам факта, что

${\color{red}{\text{Var}_{data}}} = {\color{green}{\text{Var}_{res}}}+{\color{blue}{\text{Var}_{reg}}}.$

Мы знаем, что

$R^2=\frac{\color{blue}{\text{Var}_{reg}}}{\color{red}{\text{Var}_{data}}},$

а значит

$R=\cos{\theta}.$

Красиво, не правда ли?

Произвольный базис

Как мы знаем, регрессия выполняется на базисных функциях $inline$ и её результатом есть модель

$f = \sum_i w_i f_i,$

но до сих пор мы использовали простейшие $inline$ , которые просто ретранслировали изначальные признаки без изменений, ну разве что дополняли их постоянной фичей $f_0(\mathbf{x}) = 1$ . Как можно было заметить, на самом деле ни вид $inline$ , ни их количество ничем не ограничены главное чтобы функции в базисе были линейно независимы. Обычно, выбор делается исходя из предположений о природе процесса, который мы моделируем. Если у нас есть основания полагать, что точки $\{(x_1,y_1),\cdots,(x_N,y_N)\}$ ложатся на параболу, а не на прямую, то стоит выбрать базис $inline$ . Количество базисных функций может быть как меньшим, так и большим, чем количество изначальных фич.

гифка

Регрессия в полиномиальном базисе. Выделенная часть кода демонстрирует использование стандартных функций scikit-learn для выполнения регрессии полиномами разной степени, снизу визуализация результата работы.

Если мы определились с базисом, то дальше действуем следующим образом. Мы формируем матрицу информации

$\Phi = \begin{pmatrix} - & \boldsymbol{f}^{(1)\top} & - \\ \cdots & \cdots & \cdots\\ - & \boldsymbol{f}^{(N)\top} & - \end{pmatrix} = \begin{pmatrix} {f}_{0}\left(\mathbf{x}^{(1)}\right) & {f}_{1}\left(\mathbf{x}^{(1)}\right) & \cdots & {f}_{n}\left(\mathbf{x}^{(1)}\right) \\ \cdots & \cdots & \cdots & \cdots\\ {f}_{0}\left(\mathbf{x}^{(1)}\right) & {f}_{1}\left(\mathbf{x}^{(N)}\right) & \cdots & {f}_{n}\left(\mathbf{x}^{(N)}\right) \end{pmatrix},$

записываем функцию потерь

$E(\mathbf{w})={\|{\boldsymbol{\epsilon}}(\mathbf{w})\|}^2={\|\mathbf{y}-\Phi \, \mathbf{w}\|}^2$

и находим её минимум, например с помощью псевдообратной матрицы

$\hat{\mathbf{w}} = \text{argmin}_\mathbf{w} \,E(\mathbf{w}) = (\Phi^{\top}\Phi)^{-1}\Phi^{\top}\mathbf{y}=\Phi^{+}\mathbf{y},$

или другим методом.

Заключительные замечания

Проблема выбора размерности

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

Есть два способа выйти из ситуации. Первый: последовательно наращивать количество базисных функций, проверять качество регрессии и вовремя остановиться. Или же второй: выбрать функцию потерь, которая определит число степеней свободы автоматически. В качестве критерия успешности регрессии можно использовать коэффициент детерминации, о котором уже упоминалось выше, однако, проблема в том, что $inline$ монотонно растет с ростом размерности базиса. Поэтому вводят скорректированный коэффициент

$\bar{R}^2=1-(1-R^2)\left[\frac{N-1}{N-(n+1)}\right],$

где $inline$ размер выборки, $inline$ количество независимых переменных. Следя за $\bar{R}^2$ , мы можем вовремя остановиться и перестать добавлять дополнительные степени свободы.

Вторая группа подходов регуляризации, самые известные из которых Ridge( $inline$ /гребневая/Тихоновская регуляризация), Lasso( $inline$ регуляризация) и Elastic Net(Ridge+Lasso). Главная идея этих методов: модифицировать функцию потерь дополнительными слагаемыми, которые не позволят вектору коэффициентов $\mathbf{w}$ неограниченно расти и тем самым воспрепятствуют переобучению

$\begin{aligned} E_{\text{Ridge}}(\mathbf{w})&=\text{SSE}(\mathbf{w})+\alpha \sum_i |w_i|^2 = \text{SSE}(\mathbf{w})+\alpha \| \mathbf{w}\|_{L_2}^2,\\ E_{Lasso}(\mathbf{w})&=\text{SSE}(\mathbf{w})+\beta \sum_i |w_i| =\text{SSE}(\mathbf{w})+\beta \| \mathbf{w}\|_{L_1}, \\ E_{EN}(\mathbf{w})&=\text{SSE}(\mathbf{w})+\alpha \| \mathbf{w}\|_{L_2}^2+\beta \| \mathbf{w}\|_{L_1}. \\ \end{aligned}$

где $\alpha$ и $\beta$ параметры, которые регулируют силу регуляризации. Это обширная тема с красивой геометрией, которая заслуживает отдельного рассмотрения. Упомяну кстати, что для случая двух переменных при помощи вероятностной интерпретации можно получить Ridge и Lasso регрессии, удачно выбрав априорное распределения для коэффициента $inline$

$y = a + bx + \epsilon,\qquad \epsilon \sim \mathcal{N}(0,\,\sigma^{2}),\qquad \left\{\begin{aligned} &b \sim \mathcal{N}(0,\,\tau^{2})&\leftarrow\text{Ridge},\\ &b \sim \text{Laplace} (0,\,\tau^{2})&\leftarrow\text{Lasso}. \end{aligned}\right.$

Численные методы

Скажу пару слов, как минимизировать функцию потерь на практике. SSE это обычная квадратичная функция, которая параметризируется входными данными, так что принципиально ее можно минимизировать методом скорейшего спуска или другими методами оптимизации. Так, реализация Lasso регрессии в scikit-learn использует метод координатного спуска.

Также можно решить нормальные уравнения с помощью численных методов линейной алгебры. Эффективный метод, который используется в scikit-learn для МНК нахождение псевдообратной матрицы с помощью сингулярного разложения. Поля этой статьи слишком узки, чтобы касаться этой темы, за подробностями советую обратиться к курсу лекций К.В.Воронцова.

Реклама и заключение

Эта статья сокращенный пересказ одной из глав курса классического машинного обучения в Киевском академическом университете (преемник Киевского отделения Московского физико-технического института, КО МФТИ). Автор статьи помогал в создании этого курса. Технически курс выполнен на платформе Google Colab, что позволяет совмещать формулы, форматированные LaTeX, исполняемый код Python и интерактивные демонстрации на Python+JavaScript, так что студенты могут работать с материалами курса и запускать код с любого компьютера, на котором есть браузер. На главной странице собраны ссылки на конспекты, рабочие тетради для практик и дополнительные ресурсы. В основу курса положены следующие принципы:

все материалы должны быть доступны студентам с первой пары;
лекция нужны для понимания, а не для конспектирования (конспекты уже готовы, нет смысла их писать, если не хочется);
конспект больше чем лекция (материала в конспектах больше, чем было озвучено на лекции, фактически конспекты представляют собой полноценный учебник);
наглядность и интерактивность (иллюстрации, фото, демки, гифки, код, видео с youtube).

Если хотите посмотреть на результат загляните на страничку курса на GitHub.

Надеюсь Вам было интересно, спасибо за внимание.

Подробнее..

Категории: Математика , Машинное обучение , Линейная регрессия , Метод наименьших квадратов , Линейная алгебра

	Русский
	English

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

Из песочницы Moneyball и Формула-1 модель прогнозирования результатов квалификаций

Результаты

Простыми словами о простых линейных функциях

Основы линейной регрессии

Содержание

Введение

Метод наименьших квадратов

Математический анализ

Статистика

Теория вероятностей

Мультилинейная регрессия

Линейная алгебра

Произвольный базис

Заключительные замечания

Проблема выбора размерности

Численные методы

Реклама и заключение

Категории

Последние комментарии