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

Пиша Расписывая первые строчки этой заметки на 100+ страниц я не могу не ответить на важный вопрос, во первых, самому себе, во вторых, читателю - “Зачем опять писать про линейную регрессию”. Во первых, потому что так хочется. Во вторых, потому что я считаю что могу по новому и непохоже на других авторов рассказать про машинное обучение на примере линейной регрессии. Чем же хочется выделится на этой сцене:

  • Предоставить материал просто насколько это возможно, чтобы статья подходила в том числе и для начинающих специалистов;

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

  • Я в пункте выше писал про “комикс” - забудьте. Компьютерная наука - это среда, раскрывающаяся в динамике, поэтому я считаю, что важно показать ключевые элементы в виде анимаций;

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

  • Фокус на практике по принципу “от боли к решению” - при продвижении по вниз по тексту, каждый последующий шаг совершается для решения возникшей проблемы на текущем шаге. Я надеюсь что такая последовательность создает необходимую связность между блоками.

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

Для кого эта статья

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

О чем речь

В данном посте три повествовательных акта:

  • Что такое линейная регрессия, зачем нужна и как её построить

  • Как измерить качество модели

  • Как улучшить качество модели, если мы недовольны 

Глобально, мы поговорим о:

  1. моделировании на основе данных;

  2. аналитическом решении линейной регрессии и почему иногда им воспользоваться не получится; 

  3. способах измерения качества моделирования - визуально и при помощи метрик;

  4. многомерной (многофакторной) линейной регрессии когда для предсказания используется много данных;

  5. о вероятностной природе линейной регрессии потому что предсказания не точны и важно уметь измерять степень неуверенности;

  6. методах улучшения качества моделирования: от усложнения модели до её упрощения (регуляризации).

В частности, мы рассмотрим: 

  1. метод наименьших квадратов для парной регрессии;

  2. регрессионные метрики R2, RMSE, MAE, MAPE, SMAPE, коэффициент корреляции и коэффициент детерминации а также методы визуальной оценки, например графики остатков;

  3. метод максимального правдоподобия и предсказательный интервал;

  4. разбиение на тренировочную и тестовую выборки: зачем нужно и как делать;

  5. методы фильтрации выбросов: Консенсус случайной выборки, Махаланобисово расстояние, Локальный фактор выброса, Расстояние Кука;

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

  7. линейную алгебру стоящую за методом наименьших квадратов и то как она работает для многомерной регрессии;

  8. численные методы оптимизации, в том числе градиентный спуск;

  9. L1 и L2 регуляризацию линейных моделей;

  10. кросс валидацию и оптимизацию гиперпараметров. 

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

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

Все анимации и картинки в статье оригинальные и выполнены автором.

Маленький литературный обзор

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

Отдельно хочу выделить “Курс Линейные модели, дисперсионный и регрессионный анализ с использованием R” (см. материалы курса в гитхаб репозитории): курс охватывает большое количество тем и при этом соблюдается (на мой взгляд) баланс между простотой и формальностью изложения. 

И конечно, не обходите стороной классические работы. В данном разделе они не приведены в виде списка литературы, однако ниже на них будут ссылки. Эти статьи приведены сразу после фрагмента который на них ссылается в квадратных скобках в формате [Авторы(ы). Название. Год публикации. Ссылка на первоисточник]

Хорошая модель начинается с данных

Предположим, что у нас есть табличные данные с двумя столбцами: 

  • Количество комнат в квартире

  • Стоимость апартаментов (для удобства, стоимость будет указана в долларах $, иначе бы много нулей мелькало на экране)

Рисунок 1. Визуализация исходного набора данных о стоимости квартир (ссылка на код для генерации картинки)
Рисунок 1. Визуализация исходного набора данных о стоимости квартир (ссылка на код для генерации картинки)

То есть на момент анализа и подготовки модели, у нас уже должны быть данные. Сбор и предварительную подготовку выборки оставим за рамками текущего рассказа, тем более что в зависимости от предметной области процесс будет значительно отличаться. Главный принцип, который имеет смысл держать в голове, это “мусор на входе - мусор на выходе”, что в целом справедливо для задач машинного обучения с учителем. Для хорошей модели нужна, в первую очередь, хорошая выборка данных.

Зачем нужна модель 

Как говорил британский статистик Джордж Бокс, “Все модели неточны, но некоторые из них полезны”. И полезны они для выявления закономерностей в данных. Как только закономерности описаны в виде математического выражения (модели), мы можем использовать это для, например, генерации предсказаний (Рисунок 2). 

Рисунок 2. Превращение таблицы с данными в модель и что можно назвать моделью
Рисунок 2. Превращение таблицы с данными в модель и что можно назвать моделью

Моделирование зависимостей между данными - нетривиальная задача. Для этого могут использоваться математические модели самой разной природы, от односложных до самых современных многоступенчатых, например, нейронных сетей. Нам на данном этапе важно зафиксировать то, что моделями можно называть самые разные преобразователи одних данных (столбцы с признаками) в другие (столбец с откликом) (Рисунок 3). 

Рисунок 3. Моделью может быть (почти) что угодно (ссылка на код для генерации картинки)
Рисунок 3. Моделью может быть (почти) что угодно (ссылка на код для генерации картинки)

В случае линейной регрессии моделируются линейные взаимосвязи между массивами данных. Для парной регрессии (когда есть один предиктор и одна зависимая переменная) уравнение имеет вид:

y = b_0 + b_1 \times x

где x - предиктор, y - целевая переменная [James, G., et al. Linear Regression. An Introduction to Statistical Learning, 2021. Открытая версия книги https://www.statlearning.com/]

То есть выражение y = 1 + 10x это модель линейной регрессии. И y = 15 - 21x - тоже. Разница только в коэффициентах. Поскольку коэффициенты в уравнении являются важнейшими параметрами, они имеют свои названия:

  • b_0 - свободный член, или параметр (коэффициент) смещения

  • b_1 - коэффициент наклона

Итак, когда мы строим модель линейной регрессии мы делаем следующее предположение: 

Предположение 1. Зависимость между признаками (независимыми переменными) и откликом (зависимой переменной) линейна [Kim, Hae-Young. Statistical notes for clinical researchers: simple linear regression 1 – basic concepts, 2018. https://www.rde.ac/upload/pdf/rde-43-e21.pdf]

Пример линейной модели с уже подобранными коэффициентами смещения и наклона (почему названия именно такие обсудим чуть ниже) можно увидеть на Рисунке 4.

Рисунок 4. Модель линейной регрессии для получения предсказаний (ссылка на код для генерации картинки)
Рисунок 4. Модель линейной регрессии для получения предсказаний (ссылка на код для генерации картинки)

Таким образом на приведенном наборе данных из Рисунка 1 для того чтобы узнать цену апартаментов в долларах нужно умножить количество комнат на 10000. 

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

Немного про интерполяцию, экстраполяцию и аппроксимацию
Дополнительный рисунок 1. Разница между терминами интерполяции, экстраполяции и аппроксимации
Дополнительный рисунок 1. Разница между терминами интерполяции, экстраполяции и аппроксимации

Как построить простую модель

Итак, нам требуется подобрать коэффициенты b_0 и b_1 для уравнения ниже таким образом, чтобы прямая линия как можно точнее описывала эмпирические наблюдения (реальные данные): y = b_0 + b_1 \cdot x, где x – количество комнат в квартире, y - цена апартаментов, $

Почему именно такое уравнение и почему два коэффициента

Несмотря на кажущаяся простоту, уравнение линейной регрессии воспроизводит самые разные линейные зависимости, что видно на Рисунке 5. Для каждого набора данных оптимальной будет своя прямая.

Рисунок 5. Примеры уравнений с разными оптимальными значениями коэффициентов (ссылка на код для генерации картинки)
Рисунок 5. Примеры уравнений с разными оптимальными значениями коэффициентов (ссылка на код для генерации картинки)

Аналитическое решение

Для нахождения оптимальных значений коэффициентов воспользуемся аналитическим решением: подставим эмпирические данные из предыдущего раздела в заранее известную и выведенную умными людьми (Карл Гаусс, Адриен Лежандр) формулу. Аналитическое решение можно разбить на четыре шага (Рисунок 6). [Hastie, T., et al. Linear Methods for Regression (Chapter 3 in The Elements of Statistical Learning: Data Mining, Inference, and Prediction). 2009. https://hastie.su.domains/ElemStatLearn]

Рисунок 6. Аналитическое решение для уравнения парной линейной регрессии. В шаге два представлен Python-подобный псевдокод для расчета углового коэффициента (см. код на Python для произведения таких расчетов)
Рисунок 6. Аналитическое решение для уравнения парной линейной регрессии. В шаге два представлен Python-подобный псевдокод для расчета углового коэффициента (см. код на Python для произведения таких расчетов)

Ошибка - тоже часть модели (когда двух коэффициентов недостаточно)

Чуть выше мы отметили, что модель линейной регрессии - аппроксимационный алгоритм. Это значит, что не накладывается ограничений на то чтобы линия строго проходила через набор известных значений. То есть мы уже на этом этапе допускаем возможность несоответствия предсказанных (произведенных моделью) и наблюдаемых значений цены квартиры. И такое расхождение, важно подчеркнуть, - абсолютно нормально - в реальном мире очень сложно найти процесс порождающий данные, которые идеально ровной прямой ложатся на график  (Рисунок 7).

Рисунок 7. Реальные данные иногда невозможно описать моделью без невязки. Поэтому в уравнении линейной регрессии есть компонента ошибки (ссылка на код для генерации картинки)
Рисунок 7. Реальные данные иногда невозможно описать моделью без невязки. Поэтому в уравнении линейной регрессии есть компонента ошибки (ссылка на код для генерации картинки)

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

Как измерить качество модели

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

  • Визуальная оценка

  • Расчет метрик

Перед тем как предметно рассуждать о каждом пункте, предлагаю определить “качественность” как нашу удовлетворенность моделью. Удовлетворены же мы будем (в рамках данного поста) в том случае, если компонента ошибки будет минимально возможной.

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

Таблица 1. Сравнение ошибки на одном наблюдаемом значении (при количестве комнат равным двум) при различных значениях коэффициентов b_0 и b_1

Количество комнат

Модель (b_0 + b_1x, где x количество комнат)

Предсказанное значение

Реальное значение

Ошибка (реальное - предсказанное)

2

0 + 10 000 \cdot 2

20 000

20 000

0

2

0 + 5 000 \cdot 2

10 000

20 000

10 000

2

500 + 1 000 \cdot 2

2 500

20 000

17 500

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

Но сначала усложним немного нам самим задачу: увеличим объем выборки в трех разных вариантах, один простой и два сложных для моделирования, и ниже попробуем проанализировать модель уже для этих данных.

Рисунок 8. Три набора данных. Пример расширенных выборок A, B, C  с ценами квартир для оценки качества моделирования (ссылка на код для генерации картинки)
Рисунок 8. Три набора данных. Пример расширенных выборок A, B, C  с ценами квартир для оценки качества моделирования (ссылка на код для генерации картинки)

Рисунок 8 больше похож на реальность - квартиры бывают разные, и несмотря на то что площадь и количество комнат совпадают, совсем необязательно, чтобы цена для разных объектов получилась идентичной.

Визуальная оценка 

Воспользуемся формулой из раздела про Аналитическое решение (Рисунок 6). Подставим данные и получим следующие модели для каждого из датасетов: 

  1. A: 0 + 10000 \cdot x

  2. B: 0 + 10000 \cdot x

  3. C: 6800 + 6600 \cdot x

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

Первый график, который есть смысл привести здесь, - это классическая точечная диаграмма (диаграмма рассеяния), где на ось x наносят значения признака, а по оси y приводятся показанные разным цветом предсказанные значения и реальные наблюдения. Анализировать такой рисунок довольно просто, - чем ближе модельная линия к реальным данным, тем качественнее модель. Ещё на графике видно зависимости между откликом и независимой переменной, поскольку перед глазами имеется и признак. [Piñeiro, G., et al. How to evaluate models: Observed vs. predicted or predicted vs. observed? 2008. https://doi.org/10.1016/j.ecolmodel.2008.05.006]

Рисунок 9. Визуальная оценка качества моделирования: предсказанные значения вместе с реальными на диаграмме рассеяния (ссылка на код для генерации картинки)
Рисунок 9. Визуальная оценка качества моделирования: предсказанные значения вместе с реальными на диаграмме рассеяния (ссылка на код для генерации картинки)

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

Решить проблему сложности визуализации, когда признаков становится много, можно с помощью диаграммы рассеяния предсказания vs наблюдения. График строится так: по оси x откладывают реальные значения, а по оси y - предсказанные (Рисунок 10) [Moriasi, D. N., et al. Hydrologic and Water Quality Models: Performance Measures and Evaluation Criteria. 2015. https://web.ics.purdue.edu/~mgitau/pdf/Moriasi%20et%20al%202015.pdf]. Я встречал вариации, когда поступают наоборот: по оси x откладывают предсказанные значения. Тем не менее, свою функцию такие графики выполняют и так, и так - выбирайте на свой вкус.

Рисунок 10. Визуальная оценка качества моделирования: диаграмма рассеяния наблюдения vs предсказания (ссылка на код для генерации картинки)
Рисунок 10. Визуальная оценка качества моделирования: диаграмма рассеяния наблюдения vs предсказания (ссылка на код для генерации картинки)

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

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

Первый вариант Q-Q графика, с порядковыми статистиками: значения предсказаний упорядочиваются по возрастанию, также поступаем с выборкой отклика, и визуализируем как на Рисунке 10. Тогда будет показано столько же точек, сколько объектов в выборке.

Второй вариант, - классический Q-Q график двух выборок. В нём визуализируются не сырые упорядоченные значения а квантили, - данные разбиваются таким образом на конечное (я обычно использую около 100) число уровней. Такой график полезен, когда важно сравнить не отдельные пары “предсказание - отклик”, а общий рисунок данных (форму распределений): где находится медиана и насколько часто встречаются большие или очень маленькие значения.

Напоминание про то что такое квантиль и перцентиль

По определению из википедии: квантиль - это значение, которое заданная случайная величина не превышает с фиксированной вероятностью.

Временно вынося за скобки вероятность из определения выше, можно сказать, что квантиль - это такое значение для выборки, которое делит данные на части. Например, квантиль уровня 0.25 - это число, ниже которого (меньше) лежит 25% всех значений выборки. А квантиль уровня 0.9 - это значение, ниже которого 90% данных.

Для выборки [1, 3, 5, 7, 9] квантиль 0.5 (медиана) равен 5. Больше 5 есть только два значения: 7 и 9, и меньше тоже только два: 1 и 3.  

Квантиль уровня 0.25 ≈ 3. Квантиль уровня 0.75 ≈ 7. Смотрите пояснения на Рисунке ниже.

Дополнительный рисунок 2. Немного про квантили и перцентили
Дополнительный рисунок 2. Немного про квантили и перцентили

25й перцентиль (персентиль) также называется первым квартилем, 50й перцентиль - медиана или второй квартиль, 75й перцентиль - третий квартиль.

Рисунок  11. Визуальная оценка качества моделирования: Q-Q график. На графике цифрами и черной окантовкой показаны 25й, 50й и 75й перцентили (квантили уровней 0.25, 0.50 и 0.75 соответственно) (ссылка на код для генерации картинки)
Рисунок  11. Визуальная оценка качества моделирования: Q-Q график. На графике цифрами и черной окантовкой показаны 25й, 50й и 75й перцентили (квантили уровней 0.25, 0.50 и 0.75 соответственно) (ссылка на код для генерации картинки)

Во втором варианте, независимо от того насколько много объектов в датасете, количество отображаемых на графике точек всегда равно 99, а значит такой подход отлично масштабируется на сколь угодно большие выборки. На Рисунке 11 видно, что для датасета A квантили реальных и предсказанных значений лежат близко к диагонали - модель хорошая. В то время как  для датасета B правый хвост распределений (справа вверху) реальных и предсказанных значений расходится. Это значит, что модель хуже предсказывает цену дорогих квартир.

Для датасета C: 

  • Ниже 25го перцентиля модельные данные лежат выше реальных

  • В интерквартильном размахе (от 25го до 75го перцентиля) модельные значения лежат ниже 

  • Выше 75го перцентиля опять хвост распределения предсказанных значений лежит выше реальных 

Если вам хочется углубиться в тему Q-Q графиков, то могу порекомендовать хабр статью Q-Q Plots. От чайника до профессионала за один гайд.

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

Рисунок 12. Визуальная оценка качества моделирования: график остатков (ссылка на код для генерации картинки)
Рисунок 12. Визуальная оценка качества моделирования: график остатков (ссылка на код для генерации картинки)

Именно по графику остатков удобнее всего анализировать следующие предположения для инициализации модели линейной регрессии (Предположение 1 было выше по тексту - линейность): 

  • Предположение 2. Нормальность остатков - ошибки модели (разности между фактическими и предсказанными значениями) должны быть нормально распределены (на вопрос “Что значит нормально распределены?” - ответ будет следующий. Большинство ошибок - маленькие, около нуля. А большие ошибки - редкость, и они встречаются примерно одинаково часто в обе стороны (в плюс и в минус)

  • Предположение 3. Гомоскедастичность (постоянная дисперсия ошибок) - модель предсказывает с примерно одинаковой ошибкой и дешевые квартиры, и средний ценовой сегмент, и дорогие объекты

  • Предположение 4. Отсутствие автокорреляции остатков - ошибки должны быть независимыми друг от друга

Из Рисунка 12 видно, что для датасета B не подтверждается второе предположение - с ростом количества комнат модель ошибается все больше и больше - дисперсия остатков увеличивается (разброс значений растет слева направо). Это значит, что ошибка для модели непостоянна и зависит от значения предиктора, - модель не учитывает какую-то скрытую (от нас и модели) закономерность и потому нестабильна.

Для датасета C не соблюдается условие нормальности остатков - модель иногда систематически завышает или систематически занижает значения уходя то выше, то ниже нуля в среднем. Кроме того на графике остатков вырисовываются “фигуры” (паттерны на графике) - это может быть признаком того что ошибки не независимы друг от друга (справедливости ради замечу, что может и не быть XD Но в любом случае, это сигнал, что что-то с моделью не ладно).

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

Рисунок 13. Визуальная оценка качества моделирования: графики остатков вместе с частотным распределением (ссылка на код для генерации картинки)
Рисунок 13. Визуальная оценка качества моделирования: графики остатков вместе с частотным распределением (ссылка на код для генерации картинки)
Я не знаю что такое частотное распределение

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

Частотное распределение - это упорядоченное представление, показывающее, сколько раз значения случайной величины попадают в определенные интервалы.

Для его построения:

  1. Весь диапазон значений разбивают на k интервалов (интервальные частоты).

  2. Подсчитывают, сколько наблюдений попало в каждый интервал - это называется абсолютной частотой.

  3. Делят абсолютную частоту на общий объем выборки n - получают относительную частоту.

На рисунке такая последовательность действий для переменной V выглядит следующим образом: 

Дополнительный рисунок 3. Визуализация частотного распределения V в виде гистограммы: методика расчета (ссылка на код для генерации картинки)
Дополнительный рисунок 3. Визуализация частотного распределения V в виде гистограммы: методика расчета (ссылка на код для генерации картинки)

Такую же визуализацию можно сделать и для переменной U, но в данном разделе сконцентрируемся только на V, для примера. Далее по тексту я график переверну на бок, чтобы удобнее было соотносить исходные данные и вертикальную нотацию визуализаций распределений.

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

Дополнительный рисунок 4. Визуализации частотного распределения с разным количеством интервалов k: 5, 10, 20. Вертикальная ось намеренно не подписана, чтобы не возникло соблазна интерпретировать положение точек на графике по оси ординат: значения этой размерности могут быть какими угодно, на распределение V это не влияет (ссылка на код для генерации картинки)
Дополнительный рисунок 4. Визуализации частотного распределения с разным количеством интервалов k: 5, 10, 20. Вертикальная ось намеренно не подписана, чтобы не возникло соблазна интерпретировать положение точек на графике по оси ординат: значения этой размерности могут быть какими угодно, на распределение V это не влияет (ссылка на код для генерации картинки)

Существуют эмпирические формулы, которые позволяют подобрать адекватное значение количества интервалов в зависимости от объема выборки. Например формула Стерджеса или формула Райса (см. Рисунок ниже) [Sturges. The Choice of a Class Interval. 1926. DOI: 10.1080/01621459.1926.10502161], [Lane, David M., et. al. Histograms. https://onlinestatbook.com/2/graphing_distributions/histograms.html].

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

Можно использовать альтернативу - визуализировать с помощью ядерной оценки плотности. Ядерная оценка плотности - это сглаженный аналог гистограммы. Вместо прямоугольных «столбиков» применяется непрерывная кривая, построенная на основе суммы «ядер» гладких функций - как правило, нормальных распределений.

Дополнительный рисунок 6. Ядерная оценка плотности распределения для выборки V (ссылка на код для генерации картинки)
Дополнительный рисунок 6. Ядерная оценка плотности распределения для выборки V (ссылка на код для генерации картинки)

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

Дополнительная анимация 1. Представление ядерной оценки плотности как выпадение насыпи (ссылка на код для генерации анимации)
Дополнительная анимация 1. Представление ядерной оценки плотности как выпадение насыпи (ссылка на код для генерации анимации)

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

Еще одним популярным способом представлять информацию о распределениях являются ящики с усами (Boxplot). Ящик с усами - это график, обобщённо описывающий распределение по квартилям. Он отображает:

  • Медиану (второй квартиль)

  • Первый Q1 и третий Q3 квартиль (25й и 75й перцентили соответственно) - границы «ящика»

  • «Усы» - границы диапазона данных без выбросов

  • Отдельные точки - выбросы

Дополнительный рисунок 7. Визуализация распределения V. Ящики с усами (ссылка на код для генерации картинки)
Дополнительный рисунок 7. Визуализация распределения V. Ящики с усами (ссылка на код для генерации картинки)

Чтобы подытожить все вышесказанное, визуализируем выборки разного размера и формы всеми рассмотренными способами. Для этого будем рассматривать разные теоретические распределения из которых будем семплировать (извлекать) выборки разного размера: по 30 и 500 элементов.

Дополнительный рисунок 8. Визуализация выборок различными способами (частотные гистограммы, ядерная оценка плотности, ящики с усами) для двух теоретических распределений: нормальное и бимодальное (смесь двух нормальных) (ссылка на код для генерации картинки)
Дополнительный рисунок 8. Визуализация выборок различными способами (частотные гистограммы, ядерная оценка плотности, ящики с усами) для двух теоретических распределений: нормальное и бимодальное (смесь двух нормальных) (ссылка на код для генерации картинки)

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

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

Продолжая тему работы с распределениями: распределение остатков, как мы заметили выше должно быть похоже на нормальное. И потому имеет смысл подготовить визуализации сравнения распределений - теоретического нормального и наблюдаемого у остатков модели. Удобно это делать через графики плотности распределения и квантильные диаграммы (квантили остатков Vs квантили нормального распределения). Оценка параметров нормального распределения производится по выборке остатков.  Поскольку такие графики работают лучше всего на больших объемах данных, то я предлагаю для наглядности искусственно увеличить размер выборок у остатков до 500 элементов в каждой, сохранив при этом основные черты поведения остатков модели для каждого датасета (Рисунок 14). 

Рисунок 14. Квантильная диаграмма с нормальным распределением и остатками модели (ряд снизу). Для наглядности размер выборок остатков для датасетов A, B и C был искусственно увеличен (ссылка на код для генерации картинки)
Рисунок 14. Квантильная диаграмма с нормальным распределением и остатками модели (ряд снизу). Для наглядности размер выборок остатков для датасетов A, B и C был искусственно увеличен (ссылка на код для генерации картинки)

Как видно из рисунка 14, для датасета A* и B* распределение остатков хорошо аппроксимируется нормальным, хотя и для B* хвосты немного расходятся - модель допускает большие ошибки немного чаще, чем нам этого хотелось бы. В случае же с бимодальным C* различия особенно заметны - распределение остатков совсем не похоже на нормальное. 

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

Подытожим. Модель редко бывает идеальной - она будет ошибаться. А анализ ошибок при помощи графиков - удобный способ диагностики модели:

  • Если предиктор всего один, то удобно построить график где показаны предсказанные значения вместе с реальными на оси y, а на оси x - признак. Сразу видно форму взаимосвязи между признаком и откликом;

  • В дополнение строим график отклик (ось x) vs предсказания (ось y), на котором чем ближе данные к прямой линии, выходящей из левого нижнего угла, тем лучше. Еще эта визуализация удобна тем, что не зависит от количества признаков в модели;

  • Хочется рассмотреть распределения предсказаний и откликов целиком, а не отдельными парами - выбираем Q-Q график;

  • Если выборки очень большие, то можно снизить когнитивную нагрузку путем группировки  данных по уровням квантилей на Q-Q графике и тогда можно например, обойтись всего 100 точками;

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

  • Для более удобного сравнения распределения остатков с теоретическим нормальным - используем квантильную диаграмму.

Метрики

Дисклеймер про обозначения X и Y

В визуализациях этого раздела вы можете встретить не совсем привычные (в сравнении со смежной по теме литературой) обозначения. Так например, предсказанные значения будут обозначаться как X в то время как отклик - Y. Сделано это намеренно чтобы у читателя несмотря на привязку к конкретной теме анализа качества модели не сложилось впечатление будто метрики можно рассчитывать только для пары “предсказания - наблюдения”. Вместо абстрактных X и Y могут быть абсолютно любые массивы - все зависит от решаемой вами задачи.

Другой причиной использовать именно эту пару является заметное отличие в написании между x и y, на графиках и анимациях легче отличить именно такие пары, чем например U и V, или классическую y и ŷ. Кроме того на анимациях ниже переменные X и Y размещены на соответствующих осях графиков.

Несмотря на эффектность визуальных методов диагностики, по настоящему эффективно измерить качество моделирования можно только в связке с метриками (численными показателями качества модели).  Расчет подходящей метрики привлекателен тем, что вместо производства дополнительной когнитивной нагрузки на инженера / исследователя через генерацию графиков для анализа, мы ее снимаем - сводим задачу анализа к оценке одного числа (Рисунок 15).

Рисунок 15. Зачем нужны метрики - они помогают одним числом (иногда несколькими числами) оценить насколько хорошей получилась модель. На графике показана метрика Средняя абсолютная процентная ошибка (Mean absolute percentage error - MAPE) (ссылка на код для генерации картинки)
Рисунок 15. Зачем нужны метрики - они помогают одним числом (иногда несколькими числами) оценить насколько хорошей получилась модель. На графике показана метрика Средняя абсолютная процентная ошибка (Mean absolute percentage error - MAPE) (ссылка на код для генерации картинки)

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

В раздел “Метрики” также входят тесты, которые помогают оценить статистическую значимость полученных коэффициентов в частности и модели в целом (об этом тоже поговорим ниже по тексту). 

Начнем с неисчерпывающего перечня: 

  • Коэффициент детерминации R2 - [Kvalseth, Tarald O. Cautionary Note about R². 1985. https://www.tandfonline.com/doi/abs/10.1080/00031305.1985.10479448];

  • Смещение (Bias);

  • Средняя абсолютная ошибка (Mean absolute error - MAE);

  • Корень из средней квадратичной ошибки (Root mean square error - RMSE);

  • Средняя абсолютная процентная ошибка (Mean absolute percentage error - MAPE);

  • Симметричная средняя абсолютная процентная ошибка (Symmetric mean absolute percentage error - SMAPE);

  • F тест (критерий Фишера, или F критерий) для проверки значима ли модель целиком;

  • t критерий для проверки значимости предикторов и отклика;

  • Критерий Дарбина - Уотсона для анализа остатков (Durbin-Watson test).

На Рисунке 16 показаны метрики рассчитанные путем сопоставления действительных значений цены квартир с предсказанными

Рисунок 16. Некоторые метрики моделей на датасетах A, B и C. Внимание - шкала по оси Y для каждого столбца в нижних трех подграфиках поделена по цветам. Поэтому смотреть имеет смысл на высоту столбцов только при сравнении значений метрик для разных датасетов, например значение коэффициента детерминации для датасетов A, B и C, а не MAE и коэффициент корреляции для одного датасета (ссылка на код для генерации картинки)
Рисунок 16. Некоторые метрики моделей на датасетах A, B и C. Внимание - шкала по оси Y для каждого столбца в нижних трех подграфиках поделена по цветам. Поэтому смотреть имеет смысл на высоту столбцов только при сравнении значений метрик для разных датасетов, например значение коэффициента детерминации для датасетов A, B и C, а не MAE и коэффициент корреляции для одного датасета (ссылка на код для генерации картинки)

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

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

Вторая группа - зеленая. Эти метрики показывают ошибку в тех же величинах, что и отклик, то есть в данном случае в долларах $. Для всех трех метрик интерпретация звучит так - чем ближе значение метрики к нулю, тем лучше (Анимация 2).

Анимация 2. Как рассчитать смещение, среднюю абсолютную ошибку (MAE), корень из средней квадратичной ошибки (RMSE). Обозначения: X - предсказанные значения, Y - реальные (ссылка на код для генерации анимации)
Анимация 2. Как рассчитать смещение, среднюю абсолютную ошибку (MAE), корень из средней квадратичной ошибки (RMSE). Обозначения: X - предсказанные значения, Y - реальные (ссылка на код для генерации анимации)

Из любопытного: на Рисунке 16 метрика смещения равна нулю во всех случаях, что буквально значит, что ошибки модели не смещены ни в меньшую ни в большую сторону. Потому вопрос читателю: почему это будет справедливо для модели линейной регрессии построенной по любым данным (можете попробовать поиграться с датасетами и изменять значения в исходных данных)? 

Из Анимации 2 и Рисунка 16 можно заметить, что при увеличении расхождения между массивами X и Y, RMSE сильнее реагирует на большие ошибки, чем MAE. Это происходит потому что в RMSE ошибки возводятся в квадрат.

Третья группа - синяя. Здесь показаны метрики ошибки, которые измеряются в процентах. Чем меньше, тем лучше. MAPE здесь чувствительна к ошибкам на маленьких значениях, потому что в формуле расчета разница наблюдения и предсказания делится на само значение наблюдения. Если фактическое значение маленькое, даже небольшая ошибка даст большой процент, и сильно повлияет на итоговый результат (Рисунок 17).

Анимация 3. Как рассчитать Среднюю абсолютную процентную ошибку, симметричную среднюю абсолютную процентную ошибку. Обозначения: X - предсказанные значения, Y - реальные (ссылка на код для генерации анимации)
Анимация 3. Как рассчитать Среднюю абсолютную процентную ошибку, симметричную среднюю абсолютную процентную ошибку. Обозначения: X - предсказанные значения, Y - реальные (ссылка на код для генерации анимации)
Рисунок 17. Особенности расчетов метрики MAPE и SMAPE на двух наборах данных, где значения наблюдений близки к нулю. Показаны изменения в метриках при увеличении значений наблюдений и предсказаний второго набора данных на 10 единиц
Рисунок 17. Особенности расчетов метрики MAPE и SMAPE на двух наборах данных, где значения наблюдений близки к нулю. Показаны изменения в метриках при увеличении значений наблюдений и предсказаний второго набора данных на 10 единиц

На Рисунке 17 видно, что разница, исчисляемая в единицах измерения величины между действительными и предсказанными значениями (абсолютное отклонение), не меняется в обоих случаях: она будет равна 0 в первой паре, 8 во второй и 47 для третьей. Для метрик, измеряемых в процентах, значения ошибок по понятным причинам уменьшатся поскольку абсолютные значения наблюдений растут. При этом для MAPE уменьшение при сдвиге более заметно потому что величины нормируются на значения выборки, а не на среднее суммы модулей как это реализовано в SMAPE. Эффект от различий в формулах в таком случае будет особенно заметен на малых значениях наблюдений и нивелироваться по мере удаления от нуля, что мы и видим (если хотите изучить этот момент подробнее, раскрывайте секцию “Особенности некоторых метрик”).

Особенности некоторых метрик

Нюансы в расчетах метрик крайне важно проговорить. На примере MAPE и SMAPE (и немного MAE) мы посмотрим на то как по разному ведут себя метрики на игрушечных датасетах. Надеюсь этот раздел преподнесет читателю важное правило в удобоваримой форме: перед тем как начинать решать задачу машинного обучения, хорошо подумайте какую метрику, или метрики, имеет смысл использовать для измерения качества: не все они одинаково хороши для вашей текущей задачи / данных. 

Начнем с небольшого эксперимента: используем данные из Рисунка 17, и исходный датасет наблюдения [1, 2, 3] - предсказания [1, 10, 50] будем двигать дальше и дальше от нуля путем прибавления 10 ко всем значениям на протяжении 10 итераций. И на каждом шаге будем рассчитывать три метрики: MAPE, SMAPE, MAE. Результаты можно увидеть на графике ниже: 

Дополнительный рисунок 9. Исследуем ассиметричность MAPE. Значения метрики MAPE, SMAPE (левая шкала) и MAE (правая шкала) при различной удаленности значений наблюдений и предсказаний модели от нуля. Абсолютное отклонение между наблюдениями и предсказаниями в данном эксперименте неизменно на каждом сдвиге (ссылка на код для генерации картинки)
Дополнительный рисунок 9. Исследуем ассиметричность MAPE. Значения метрики MAPE, SMAPE (левая шкала) и MAE (правая шкала) при различной удаленности значений наблюдений и предсказаний модели от нуля. Абсолютное отклонение между наблюдениями и предсказаниями в данном эксперименте неизменно на каждом сдвиге (ссылка на код для генерации картинки)

Как видно из Рисунка выше, чем большие значения включены в датасет, тем меньше различия между MAPE и SMAPE, и тем меньше ошибки, измеряемые в процентах. Выравнивание MAPE и SMAPЕ объясняется особенностями расчета, которые позволяют устранить эффект несимметричности MAPE особенно заметный именно на малых значениях наблюдений. MAE при этом неизменна, как и ожидалось.

Теперь разберем почему это называется “ассиметричностью”. Лучше всего показать на примере. Представим, что модель предсказала значение 110 когда наблюдение было равно 100. Тогда MAPE для данного наблюдения составит 10%. В случае же когда наблюдение равно 110, а предсказание 100, ошибка MAPE составит 9.1% несмотря на то что абсолютное отклонение в обоих случаях равно 10. Таким образом, метрика MAPE ассиметрична, потому что одинаковое абсолютное отклонение воспринимается по-разному, в зависимости от того, превышает ли прогноз истинное значение или наоборот.

К другим недостаткам MAPE можно отнести невозможность расчета когда есть наблюдения, равные нулю. В таком случае обычно принято заменять значения в выборке на очень малые, например 0.000001 во время вычисления метрик. Тем не менее, мы знаем, что в таком случае MAPE будет завышена.

Свои особенности есть и у других метрик. Так, например, мы уже упоминали что RMSE по сравнению с MAE, более чувствительна к большим ошибкам. Перед данным разделом не стоит задача обозреть все такие нюансы, однако основная мысль, надеюсь понятна - будьте рациональны в выборе метрик. Используйте рекомендованные метрики из вашей предметной области или, если таковых нет, переходите к общеупотребимым. 

Чтобы подытожить, - единицы измерения метрик и также диапазоны возможных значений сведены в справочную Таблицу 2.

Таблица 2. Некоторые метрики для задачи регрессии

Метрика

В чем измеряется

Диапазон

Значение

Коэффициент корреляции Пирсона (между предсказаниями и реальными данными)

Безразмерный

От -1 до 1

Чем ближе к 1, тем модель лучше

Коэффициент детерминации R2

Безразмерный

От −∞ до 1

Чем ближе к 1, тем модель лучше

Смещение (Bias)

Та же единица, что и целевая переменная

От −∞ до∞

Чем ближе к нулю, тем модель лучше

Средняя абсолютная ошибка (MAE)

Та же единица, что и целевая переменная

От 0 до ∞

Чем ближе к нулю, тем модель лучше

Корень из средней квадратичной ошибки (RMSE)

Та же единица, что и целевая переменная

От 0 до ∞

Чем ближе к нулю, тем модель лучше

Средняя абсолютная процентная ошибка (MAPE)

Проценты (%)

От 0 до ∞

Чем ближе к нулю, тем модель лучше

Симметричная средняя абсолютная процентная ошибка (SMAPE)

Проценты (%)

От 0 до 200

Чем ближе к нулю, тем модель лучше

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

Подведем промежуточный итог. Для оценки качества моделирования мы начинали с того, что составляли таблицу с предсказанными и реальными значениями (Таблица 1). Большие таблицы анализировать было бы тяжело, поэтому мы привели информацию в более удобоваримый вид при помощи графиков, то есть перешли к визуальной оценке (Рисунки 9-14). Затем упростили задачу еще немного - вместо проведения экспертной оценки по графикам мы рассчитали метрики (Рисунки 15, 16, 17 и Анимации 1-3). Но и тут есть подвох - пусть число(а) мы получили, все ещё наша ответственность - решить хорошая метрика или нет. Для Рисунка 15 использовался порог в 5% для MAPE. Однако эту эвристику никак нельзя обобщить на все линейные регрессии в мире: данные разные, (бизнес-)задачи разные и т.д. Для некоторых данных хорошей будет считаться такая модель, ошибка которой не превысила порог 7.5%, для других - 11.2%. 

F критерий

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

Если вы до этого никогда не встречались со статистическими тестами, имеет смысл начать с упрощенного определения. Статистический тест - это способ проверить, случайность ли мы наблюдаем, или существует закономерность. Воспринимайте такой тест как черный ящик, который принимает на вход некоторые данные, и по определенным формулам рассчитывает ответ - несколько промежуточных значений (статистик и p-value) и вердикт (Рисунок 18) [Sureiman, Onchiri, et al. F-test of overall significance in regression analysis simplified. 2020. https://www.tandfonline.com/doi/full/10.1080/00031305.2016.1154108].

Рисунок 18. Статистические тесты и проверка гипотез. Приведен пример расчета для F критерия проверки значимости модели. Входные данные показаны оранжевым, то, что получается из расчетов теста - желтым (ссылка на код для расчетов)
Рисунок 18. Статистические тесты и проверка гипотез. Приведен пример расчета для F критерия проверки значимости модели. Входные данные показаны оранжевым, то, что получается из расчетов теста - желтым (ссылка на код для расчетов)

Как видно из Рисунка 18, важно перед тестированием определить порог (да, самое время поговорить про нюанс, нам и тут с порогом нужно повозиться, но тут значительно проще, - есть общепринятые наборы значений из которых можно выбрать), значение которого сильно влияет на заключительный ответ. Такой порог называется уровнем значимости. Значение 0.05 означает, что мы допускаем 5% шанс ошибочного отклонения нулевой гипотезы (гипотеза может звучать так: модель не лучше, чем наивное предсказание средним). Этот порог мы можем варьировать, например, в некоторых научных дисциплинах используют 0.01 или даже 0.001 (более строго), а в других 0.10 (менее строго). Если вам не до конца понятен “физический смысл” уровня значимости на этом этапе, - ничего страшного. В конце этого раздела есть подробное пояснение. На данный же момент будет достаточно зафиксировать следующее - статистические тесты которые мы будем обсуждать далее имеют параметр α, который мы, как исследователи-инженеры определяем под задачу. В нашем случае он равен 0.05.

Итак, статистический критерий позволяет взять данные и некоторые параметры, на основе чего рассчитывает статистики, которые используются для сравнения (больше ли статистика порогового значения или меньше). На основе этого сравнения делается вывод о значимости модели. Я рекомендую не переизобретать велосипед здесь и пользоваться статистическими пакетами для расчета критериев (отчасти поэтому я здесь не привожу формулы для расчета), - так надежнее. Что же касается вопроса какие значения сравнивать - F статистику с критическим значением F или p-value с уровнем значимости, лично я тяготею ко второму варианту в силу привычки. 

Мы можем использовать F тест для ответа на вопрос “Значима ли модель?”. Поскольку статистика - это математическая наука, предлагаю ниже формально описать две точки зрения на получившуюся модель. Статистический тест поможет нам выбрать наиболее вероятную из этих гипотез. 

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

Альтернативная гипотеза (H_1) тогда будет звучать так: Хотя бы один коэффициент не равен нулю. Модель значима ведь она объясняет какую-то часть вариации целевой переменной.

Теперь проведем тесты на наших трех датасетах: A, B и C (Рисунок 19).

Рисунок 19. Расчет F критерия проверки значимости модели. Показаны результаты тестирования с применением пакета Python statsmodels для уровня значимости 0.05. x - признак модели (количество комнат). Модели разные, но каждая из них в целом статистически значима (ссылка на код для генерации картинки)
Рисунок 19. Расчет F критерия проверки значимости модели. Показаны результаты тестирования с применением пакета Python statsmodels для уровня значимости 0.05. x - признак модели (количество комнат). Модели разные, но каждая из них в целом статистически значима (ссылка на код для генерации картинки)

Как мы видим из Рисунка 19, во всех трех случаях p value меньше 0.05 (уровень значимости, мы его выбрали именно таким потому что это общепринятый порог по умолчанию и поскольку принятие неверной гипотезы в случае предсказания цены квартиры не влечет страшных последствий, как например в медицинских данных, предлагаю не делать его строже), значит мы отвергаем нулевую гипотезу H₀ для моделей A, B и C. После такой проверки мы можем говорить о статистической значимости моделей в целом - хотя бы один из признаков статистически связан с откликом.

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

К недостаткам такого метода оценки качества можно отнести его узкоспециализированность: F-критерий - это параметрический тест, разработанный специально для линейных моделей, то есть его не применить к, скажем, алгоритму случайного леса (это такой алгоритм в машинном обучении), в отличие от MAPE или MAE. Да и для линейных моделей этот статистический тест тоже требует выполнения стандартных предположений (см. Предположения 2-4 выше по тексту: независимость наблюдений, нормальность остатков и гомоскедастичность). Тем не менее, если вы заинтересовались темой, можете самостоятельно изучить материалы, например, по t-критерию для каждого предиктора, когда гипотезы проверяется для каждого коэффициента модели, а также критерий Дарбина - Уотсона. Или выберите любой другой критерий для изучения самостоятельно, мы здесь рассмотрели только лишь основной принцип. P.S. можете уделить особенное внимание формулам расчета статистик и той математике, которая за ними стоит.

Если не до конца понятно про уровень значимости a, пожалуйста посмотрите этот раздел

Каждый раз когда я пытался разобраться в том что значит уровень значимости, я сталкивался с невидимой стеной. Примеры посложнее приводили выкладки, которые я не понимал. Источники же “попроще” доносили - вот смотри пример, где все интуитивно понятно: 

  • H₀ (нулевая гипотеза): У пациента нет рака

  • Ошибка первого рода: Тест говорит "рак есть", хотя его на самом деле нет

  • Если уровень значимости α принимается равным 0.05 - В 5% случаев тест может ошибочно напугать здорового человека, сообщив, что у него рак.

  • Поэтому в медицине часто выбирают низкий α (например, 0.01), чтобы свести к минимуму ложные тревоги

Но у нас то здесь данные и коэффициенты модели - всё фиксировано. Мы применяем F-критерий и получаем p-value < 0.05. Можно хоть 100 раз провести этот тест, - результат будет таким же, ведь модель одна и та же, и коэффициенты одни и те же. Вот пожалуйста - 100 раз получим подтверждение, что модель значима. Причем тут вообще порог в 5 процентов? Откуда берётся эта “вероятность”?

Давайте вместе разберемся. Начнем с фразы - "Модель значима на уровне 0.05". Она, несмотря на формулировку, не относится к конкретной модели, а скорее определяет уровень нашего доверия к данным, на которых модель обучалась. То есть если мы из большого внешнего мира соберем данные и обучим модель, затем соберем новую порцию и обучим еще модель, и сделаем так много раз, то в n случаев мы получим статистически значимую связь, даже если в реальном мире никакой связи на самом деле между переменными нет. И уровень значимости помогает обработать такие случаи.

Подытожим, 0.05 p-value: даже если связи на самом деле нет, примерно в 5 случаях из 100 тест всё равно скажет “связь есть” из-за случайности в данных.

Чтобы разбавить текст, позвольте продемонстрировать это в виде анимации - мы сгенерируем 100 точек случайным образом и попытаемся из этой совокупности извлекать наборы данных по 30 объектов и строить по ним модель линейной регрессии. Повторим извлечение 20 раз. В таком случае при уровне значимости 5% мы допускаем, что 1 раз из 20 мы при помощи F критерия посчитаем что модель значима, хотя на самом деле взаимосвязи между величинами нет.

Дополнительная анимация 2. Разбираемся в том какой смысл в уровне значимости при тестировании моделей линейной регрессии. Генеральная совокупность сгенерирована случайным образом. Показаны результаты при уровне значимости 0.05 (ссылка на код для генерации анимации)
Дополнительная анимация 2. Разбираемся в том какой смысл в уровне значимости при тестировании моделей линейной регрессии. Генеральная совокупность сгенерирована случайным образом. Показаны результаты при уровне значимости 0.05 (ссылка на код для генерации анимации)

Правда, в случае когда никакой связи между x и y в действительности не было в 1 случае из 20 мы получили результат p-value при тестировании меньше 0.05. Выбери мы порог уровня значимости строже, например, 0.01 и мы бы избежали ошибки первого рода (случай когда мы отвергаем гипотезу H₀ (связи между x и y нет) и принимаем альтернативную (связь есть), хотя H₀ на самом деле верна).

И для сравнения, - сгенерируем генеральную совокупность, где линейная зависимость явно присутствует и повторим эксперимент: 20 выборок и те же 20 попыток построить модель линейной регрессии.

Дополнительная анимация 3. Разбираемся в том какой смысл в уровне значимости при тестировании моделей линейной регрессии. Генеральная совокупность имеет линейную взаимосвязь. Показаны результаты при уровне значимости 0.05 (ссылка на код для генерации анимации)
Дополнительная анимация 3. Разбираемся в том какой смысл в уровне значимости при тестировании моделей линейной регрессии. Генеральная совокупность имеет линейную взаимосвязь. Показаны результаты при уровне значимости 0.05 (ссылка на код для генерации анимации)

Итогом данной обзорной главы про регрессионные метрики и F критерий обозначим:

  • Помимо визуальных методов оценки ошибки прогноза можно пользоваться метриками. Они удобны тем, что одним числом позволяют охарактеризовать модель - хорошая получилась или не очень;

  • Метрики используются при оптимизации модели и потому важно понимать их свойства, например:

    • Метрики “зеленой группы” (RMSE, MAE, смещение) удобны тем что предоставляют число в исходных единицах измерения

    • Корень из средней квадратичной ошибки (RMSE) сильнее реагирует на большие ошибки и выбросы, чем средняя абсолютная ошибка (MAE)

    • “Синяя группа” (MAPE и SMAPE) - вычисляются в процентах и поэтому их часто удобно обсуждать в бизнес-контексте. При этом, когда значения отклика близки к нулю, эти метрики могут становиться нестабильными и давать вводящие в заблуждение оценки.

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

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

Неопределенность прогноза. Предсказательный интервал

На мой взгляд, интересной комбинацией визуальной оценки и формальных метрик является построение предсказательного интервала. Предсказательный интервал - это диапазон значений, внутри которого с заданной вероятностью окажется новое наблюдение. Он позволяет показать неопределённость прогноза, объединяя статистические метрики и наглядность визуализаций (Рисунок 20).

Рисунок 20. Точечное предсказание и предсказательный интервал
Рисунок 20. Точечное предсказание и предсказательный интервал

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

Вспомним. Для получения точечного прогноза в модель подставляют значения признаков (в случае парной регрессии - одного) и вычисляют предсказание. Однако прогноз крайне редко бывает абсолютно точным: как правило, присутствует случайная ошибка модели. 

При инициализации модели мы предполагаем, что делаем несерьезные (маленькие) ошибки чаще больших и при этом ошибки одинаково вероятны и в большую и в меньшую стороны. На основе этих двух предположений выстраивается статистико-вероятностный взгляд на тему инициализации модели линейной регрессии как неразрывной связки значений коэффициентов модели и распределения ошибок (Рисунок 21) [Fisher, R. A. On the Mathematical Foundations of Theoretical Statistics. 1922. https://doi.org/10.1098/rsta.1922.0009]. 

Рисунок 21. Метод максимального правдоподобия как способ оценки коэффициентов модели линейной регрессии. На примере упрощенной модели только со свободным членом (ссылка на код для генерации основных элементов картинки)
Рисунок 21. Метод максимального правдоподобия как способ оценки коэффициентов модели линейной регрессии. На примере упрощенной модели только со свободным членом (ссылка на код для генерации основных элементов картинки)

Как видно из Рисунка 21, разброс ошибок модели можно оценить путем расчета стандартного отклонения ошибок “сигма” (\sigma). Тут же можно говорить и о дисперсии ошибок - тоже подходящая метрика для измерения разброса, стандартное отклонение это квадратный корень из дисперсии.  Чем больше значение стандартного отклонения \sigma, тем выше неопределенность прогноза (см. секцию 2 на Рисунке 21).

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

  • Шум в данных: чем больше шума, тем больше неопределенность

  • Объем выборки: чем больше данных модель “видела” при обучении, тем надёжнее оценены коэффициенты и уже интервал

  • Удаленность точки от центра данных: чем дальше новое значение признака от среднего, тем выше неопределённость.

Упрощенно, алгоритм построения предсказательного интервала выглядит так: 

  1. Инициализируем модель (по формуле из раздела выше - Рисунок 6)

  2. Считаем компоненту ошибки - остатки

  3. Из остатков рассчитываем типичный размер ошибки s

  4. Получаем  точечный прогноз

  5. Умножаем s на поправки “сколько данных было в обучении”, “как далеко значения признака​ от центра данных”, “какой уровень доверия модели” (это вероятность того, что интервал накрывает значение параметра. Этот параметр мы выберем сами исходя из особенностей задачи, принцип здесь такой же как и уровня значимости при статистическом тестировании выше по тексту, классическое значение по умолчанию - 0.95)

Рассмотрим простейший пример, сгенерируем набор данных в 30 элементов с идеальной линейной зависимостью между признаком и откликом, построим модель и рассчитаем предсказательный интервал. А затем 1) добавим шум в данные, 2) увеличим объем выборки, 3) увеличим уровень доверия (от 90% до 99%, где при 99% ширина предсказательного интервала будет максимальной) - см. Анимация 4.

Анимация 4. Предсказательный интервал и его зависимость от особенностей обучающей выборки (размера и зашумленности) и уровня доверия (ссылка на код для генерации анимации)
Анимация 4. Предсказательный интервал и его зависимость от особенностей обучающей выборки (размера и зашумленности) и уровня доверия (ссылка на код для генерации анимации)

И отдельно рассмотрим как выглядит предсказательный интервал для датасетов A, B, C (Рисунок 22).

Рисунок 22. Предсказательные интервалы с разным уровнем доверия для моделей построенных на датасетах A, B, C (ссылка на код для генерации картинки)
Рисунок 22. Предсказательные интервалы с разным уровнем доверия для моделей построенных на датасетах A, B, C (ссылка на код для генерации картинки)

На Рисунке 22 хорошо видно, что несмотря на то, что коэффициенты у модели A и B одинаковые, предсказательные интервалы отличаются по ширине - для датасета B он намного более широкий. А если говорить про абсолютную ширину предсказательного интервала, самая большая, ожидаемо, получена для модели на датасете C. 

Разбиение на обучение и тест и оценка метрик 

Все вышеперечисленные оценки качества моделирования затрагивали изучение поведения модели на объектах, по которым проводилось обучение модели. Однако на практике мы хотим быть уверенными, что модель будет показывать достойный результат и на новых данных, которые она “до этого не видела”. 

Поэтому общепринятой практикой в машинном обучении является разбиение исходной выборки на части (подвыборки), на одной из которых модель “обучается” - обучающая, а на второй проверяется её обобщающая способность - тестовая (Рисунок 23).

Рисунок 23. Концепт разбиения выборки на обучающую и тестовую. В большом количестве случаев рекомендовано производить разбиение случайно, а не, например, так: первые 70 процентов датасета в обучение и последующие 30 в тест потому что данные могут быть упорядочены
Рисунок 23. Концепт разбиения выборки на обучающую и тестовую. В большом количестве случаев рекомендовано производить разбиение случайно, а не, например, так: первые 70 процентов датасета в обучение и последующие 30 в тест потому что данные могут быть упорядочены

Если скомбинировать методы диагностики качества моделирования в одну большую визуализацию, то получится следующая картина:

Рисунок 24. Оценка качества моделирования на обучающей и тестовой выборках с графиками остатков и метриками (осторожно, это насыщенный информацией график - его нужно рассматривать по чуть чуть, если есть затруднения - переходите в текст ниже). Предсказательные интервалы показаны для уровня доверия 95% и рассчитаны по обучающей выборке (ссылка на код для генерации картинки)
Рисунок 24. Оценка качества моделирования на обучающей и тестовой выборках с графиками остатков и метриками (осторожно, это насыщенный информацией график - его нужно рассматривать по чуть чуть, если есть затруднения - переходите в текст ниже). Предсказательные интервалы показаны для уровня доверия 95% и рассчитаны по обучающей выборке (ссылка на код для генерации картинки)

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

  • наконец то “ожила” метрика смещения - на тестовых данных она теперь не равна нулю, как для обучения, а изменяется и в большую (датасеты A и B) и в меньшую стороны (датасет C); 

  • чем сложнее датасет (A самый простой для моделирования линейной зависимости, B сложнее и C самый сложный), тем заметнее растут метрики (или уменьшаются, в случае коэффициентов корреляции и детерминации) и “расползаются” остатки на графиках от обучения к тесту.

В текущем разделе важно сказать, что способ разбиения на обучение и тест может сильно повлиять на то, как будет выглядеть наша модель (Анимация 5).

Анимация 5. “Данные одни, а коэффициенты разные”. Визуализация влияния способов разбиения на обучение и тест для датасета B на коэффициенты линейной регрессии и метрики. Пропорции: 60% - обучение, 40% тест. x - признак модели (количество комнат) (ссылка на код для генерации анимации)
Анимация 5. “Данные одни, а коэффициенты разные”. Визуализация влияния способов разбиения на обучение и тест для датасета B на коэффициенты линейной регрессии и метрики. Пропорции: 60% - обучение, 40% тест. x - признак модели (количество комнат) (ссылка на код для генерации анимации)

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

  • Географическая или пространственная зависимость. Когда данные имеют пространственную привязку, например, измерения температуры, уровень загрязнения воздуха, урожайность на разных участках. Объекты, расположенные близко друг к другу, часто оказываются сильно коррелированными. В таких случаях тестовую выборку имеет смысл формировать из пространственно удаленных регионов, чтобы избежать переоценки качества модели;

  • Сценарное тестирование. В ряде бизнес-задач важно заранее оценить, как модель поведёт себя в определённых критичных или редких ситуациях (например, при высоких или экстремальных значениях признаков). Такие случаи можно специально включить в тестовую выборку, даже если они не представлены или слабо представлены в тренировочной.

Представим что в мире всего 45 квартир…

Чтобы нам было проще идти дальше по тексту, примем важную условность для этой статьи - представим, что наш гипотетический мир в рамках которого мы строим модели очень маленький и в нем построили всего 45 квартир. Тогда все предыдущие попытки построить модель на выборках A, B и C были ничем иным как частными шагами к воспроизведению вот этой самой “первородной” зависимости из всех доступных данных. 

Тогда выборки A, B и C - это не совсем отдельные датасеты (представим, что это данные собранные из разных городов - A, B и C соответственно), а части генеральной совокупности D. Предположим, что мы можем объединить эти выборки и работать с ними как с одним целым (Рисунок 25).

Рисунок 25. Объединение датасетов A, B и C в общий D. Представим что это все данные которые есть вообще (ссылка на код для генерации картинки)
Рисунок 25. Объединение датасетов A, B и C в общий D. Представим что это все данные которые есть вообще (ссылка на код для генерации картинки)

Важно держать в голове что все что мы делаем: делим выборки на части для обучения и тестирования, предобрабатываем данные, считаем метрики, проводим статистические тесты и все остальное - всё для того, чтобы полученная модель в конечном итоге хорошо моделировала данные из генеральной совокупности. Задача статистики (для задач машинного обучения с учителем это тоже справедливо) - сделать вывод о всей совокупности, имея данные только по выборке. 

То есть если мы вдруг построим модель которая будет идеально предсказывать цену для этих 45 объектов, мы получим инструмент всегда выдающий правильный ответ, потому что других данных, на которых модель могла бы ошибиться, в этой реальности не существует (опять же, тут все упирается в “если бы”). Вернемся к реальности и попробуем описать все данные при помощи одной модели линейной регрессии (Рисунок 26).

Рисунок 26. Модель построенная на всех данных из генеральной совокупности (“эталонная линейная модель”) - Значение метрик из картинки будем считать эталоном (с некоторыми допущениями) к которому далее по тексту мы будем стремиться (ссылка на код для генерации картинки)
Рисунок 26. Модель построенная на всех данных из генеральной совокупности (“эталонная линейная модель”) - Значение метрик из картинки будем считать эталоном (с некоторыми допущениями) к которому далее по тексту мы будем стремиться (ссылка на код для генерации картинки)

В реальном мире  собрать данные по всем квартирам физически невозможно (слишком долго, дорого и трудоёмко), поэтому мы всегда берем только часть. Так и здесь мы собирали выборки и пытались оценить зависимость между переменными таким образом, чтобы быть максимально приближенными к зависимости из D. 

Очень очень важная оговорка: Далее по тексту мы воспользуемся условностями нашего мира и иногда будем периодически подглядывать как построенная модель ведет себя на (всех данных) генеральной совокупности. Это поможет нам понять были ли наши модификации успешны, когда метрика ошибки падает, или нет, когда метрика ошибки растет. Тем не менее, пожалуйста помните, в реальном мире такого не бывает и проверить модель на всех объектах невозможно!

Улучшение качества модели

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

  • Расширение выборки - увеличение количества объектов в выборке;

  • Уменьшение выборки - исключение выбросов и других нежелательных строк из таблицы с данными;

  • Усложнение модели - включение новых признаков, как из наблюдений, так и сгенерированных;

  • Упрощение модели - сокращение числа предикторов (иногда и это помогает улучшить метрики!);

  • Настройка модели - поиск оптимальных гиперпараметров (такие параметры которые не оптимизируются во время обучения модели) для выбранной модели.

Рассмотрим каждый из методов по порядку, начиная с расширения выборки. Для демонстрации поставим эксперимент. 

Увеличение размеров выборки 

Помним, что значения из генеральной совокупности нам недоступны и мы можем извлекать их только частями. В рамках эксперимента мы будем формировать выборки по 10 и 20 квартир  случайным образом. Эксперименты повторим по 30 раз для каждого размера. Метрики будем замерять по 1) обучающей выборке 2) тестовой и 3) по всей совокупности значений - по 45 объектам. Попробуем оценить как увеличение объема выборок помогает построить более надежную модель для генеральной совокупности (Анимация 6).

Анимация 6. Анализ взаимосвязи между количеством объектов в выборке и метриками рассчитанными по генеральной совокупности. Визуализация проведения экспериментов для первых пяти попыток из 30 для каждого способа семплирования (по 10 и по 20 объектов) (ссылка на код для генерации анимации)
Анимация 6. Анализ взаимосвязи между количеством объектов в выборке и метриками рассчитанными по генеральной совокупности. Визуализация проведения экспериментов для первых пяти попыток из 30 для каждого способа семплирования (по 10 и по 20 объектов) (ссылка на код для генерации анимации)

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

Рисунок 27. Результаты экспериментов по увеличению размеров семплов: чем больше объектов в семпле, тем ближе метрики на обучении и тесте к тому что модель показывает на генеральной совокупности. И тем, качественнее модель (ссылка на код для генерации картинки)
Рисунок 27. Результаты экспериментов по увеличению размеров семплов: чем больше объектов в семпле, тем ближе метрики на обучении и тесте к тому что модель показывает на генеральной совокупности. И тем, качественнее модель (ссылка на код для генерации картинки)
Если вам более привычны ящики с усами (Boxplot) - раскрывайте раздел
Дополнительный рисунок 10. Результаты экспериментов по увеличению размеров семплов в виде ящиков с усами: чем больше объектов в семпле, тем ближе метрики на обучении и тесте к тому что модель показывает на генеральной совокупности. И тем, качественнее модель (ссылка на код для генерации картинки)
Дополнительный рисунок 10. Результаты экспериментов по увеличению размеров семплов в виде ящиков с усами: чем больше объектов в семпле, тем ближе метрики на обучении и тесте к тому что модель показывает на генеральной совокупности. И тем, качественнее модель (ссылка на код для генерации картинки)

Несмотря на то что мы здесь оперировали очень маленькими выборками (это отчасти сделано для наглядности), все же некоторые выводы, справедливые и для бОльших объемов, сделать мы можем на основе Анимации 6 и Рисунка 27. А именно:

  • Среднее значение метрики ошибки RMSE на генеральной совокупности меньше в случае когда размер сэмплов составлял 20 элементов - 4088 против 4419 при 10 объектах. Это значит, что модель, построенная на большем объеме данных, имеет меньшую ошибку на генеральной совокупности.

  • Оценки метрик более устойчивые для больших выборок - разница между метрикой RMSE на обучающей выборке, тестовой и на генеральной совокупности меньше для 20 объектов.  

Как видно, мы добились прироста в метрике на генеральной совокупности в случае когда использовали большие выборки (по 20 элементов). Так и Вы, уважаемый читатель, не забывайте измерять метрики после каких-либо изменений в данных или модели. Если изменение приводит к улучшению метрики - оставляем, наоборот - откатываем назад, - полагаемся на инженерный подход, а не на удачу! Хоть и в реальном имре замерять метрики по всем объектам у нас нет возможности, метрики на обучающих и тестовых выборках помогут выбрать правильный путь.

Уменьшение выборки (фильтрация выбросов)

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

Итак, мы попытались собрать побольше данных в выборку для обучения модели. Однако, представим, что нам не повезло: полученный семпл несмотря на свой размер в 20 элементов, не дает возможность построить модель которая была бы похожа на эталонную (Рисунок 28).

Рисунок 28. “Неудачное” семплирование из генеральной совокупности. Эталонная модель показана черной линией (ссылка на код для генерации картинки) 
Рисунок 28. “Неудачное” семплирование из генеральной совокупности. Эталонная модель показана черной линией (ссылка на код для генерации картинки

В дополнение к не самому удачному сэмплированию могут появляться и другие усложняющие нам работу факторы. Такие искажения довольно распространены в реальных данных из-за множества причин: неточности измерений, различные технические ошибки при хранении и передачи данных, человеческие ошибки. В нашем случае, представим, что для некоторых наблюдений риелторы, у которых мы просили данные, допускали ошибки при ручном вводе информации с бумажных носителей: ставили соседнюю цифру вместо 4 - 3, добавляли или убирали нолики (Рисунок 29).

Рисунок 29. Некоторые выборки содержат испорченные данные (ссылка на код для генерации картинки)
Рисунок 29. Некоторые выборки содержат испорченные данные (ссылка на код для генерации картинки)

Если по таким сырым данным построить модель, то результат будет далек от эталонной модели, и опять, мы недовольны качеством моделирования. 

На этот раз попробуем решить проблему путем исключения некоторых объектов, наиболее непохожих на остальные (выбросы). Для этого существует большое количество методов, так или иначе оперирующих простым принципом отделения похожих объектов от непохожих по некоторому порогу (Рисунок 30) [Mandic-Rajcevic, et al. Methods for the Identification of Outliers and Their Influence on Exposure Assessment in Agricultural Pesticide Applicators: A Proposed Approach and Validation Using Biological Monitoring. 2019. https://doi.org/10.3390/toxics7030037]:

  • Метод межквартильного размаха (IQR method) - непараметрический метод*

  • Метод 3-сигм - параметрический* (предполагается распределение, чаще всего нормальное)

  • Z-оценка (Z-score method) - параметрический*

  • Модифицированная Z-оценка (на основе медианы) - параметрический*

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

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

В одномерных методах (Рисунок 30) не используется информация о предикторах - только об одной переменной - отклике. И потому, в частности, видно, что такие разбиения не учитывают тренд в данных. Другим ограничением является необходимость задавать порог, будь это 1.5 перед значениями разности квартилей, 3 в формуле трех сигм или порог для Z-оценки.

Следующая особенность: три из четырех приведенных методов фильтрации выбросов полагаются на предположение о форме распределения отклика. Если данные распределены нормально, или хотя бы с одной модой и без сильной асимметрии, то метод 3-сигм, Z оценка и модифицированная Z-оценка будут выдавать адекватный результат, однако если форма распределения в какой-то степени необычна, выбросы могут оказаться не выбросами вовсе. Поскольку на Рисунке 30 форма распределения не сильно отличается от нормального купола, в данном случае можно использовать именно эти стандартные методы.

Еще из интересностей здесь можно отметить, что метод трёх сигм - это частный случай Z-оценки при пороге 3.0, просто выраженный в исходных координатах y, а не в стандартизированном виде (в пространстве z-оценки). Можете обратить внимание на график: как там соотносятся линии 2σ для метода трех сигм и линии для z-оценки при пороге 2.0.

Если применить все приведенные выше способы фильтрации уже к нашим данным, то у нас получится построить следующие модели - Рисунок 31.

Рисунок 31. Построенные модели по данным, отфильтрованным одномерными методами (ссылка на код для генерации картинки)
Рисунок 31. Построенные модели по данным, отфильтрованным одномерными методами (ссылка на код для генерации картинки)

Рассматривая Рисунок 31 заметим, что “худшая” модель с точки зрения метрики RMSE на генеральной совокупности - именно та, которая была инициализирована на данных с выбросами. Лучшая по метрике получилась модель, которая обучалась на данных очищенных по Z оценке с порогом в 1.5.

По Рисунку 31 довольно легко делать выводы про эффективность того или иного метода исключения выбросов. Однако такое впечатление обманчиво по той причине, что мы сверяемся с метриками на генеральной совокупности D, что при реальной разработке недоступно. На закономерный вопрос “А что же тогда делать” - ответ один, экспериментировать. Иногда удобнее и быстрее всего “вычистить” тестовую выборку и замерять метрику по ней. В других случаях можно считать исключение выбросов успешным, если разница ошибок на обучении и тесте уменьшается. Серебрянной пули подходящей для всех случаев тут нет.

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

Рисунок 32. Фильтрация выбросов как способ нахождения непохожих наблюдений, - разбираемся как работают методы “многомерной” фильтрации (ссылка на код для генерации картинки)
Рисунок 32. Фильтрация выбросов как способ нахождения непохожих наблюдений, - разбираемся как работают методы “многомерной” фильтрации (ссылка на код для генерации картинки)

Каждый метод из Рисунка 32 можно обсуждать отдельно поскольку они уже значительно более продвинутые по сравнению с одномерными. Я же позволю здесь ограничиться визуализацией и не идти вглубь. Мы рассмотрим эти методы сугубо с практической точки зрения - посмотрим как их использование влияет на коэффициенты и метрики модели линейной регрессии (Рисунок 33).

Рисунок 33. Построенные модели по данным, отфильтрованным многомерными методами (ссылка на код для генерации картинки)
Рисунок 33. Построенные модели по данным, отфильтрованным многомерными методами (ссылка на код для генерации картинки)

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

Рассмотрим метод расчета рычага и Критерий Кука. Рычаг - это характеристика, показывающая “насколько объект странный по оси x”. Далеко ли расположен x​_i от центра. Если далеко, то рычаг у этого объекта будет большим. Подходящая метафора: качели - чем дальше сидишь от центра, тем сильнее влияешь на процесс качания. Критерий Кука определяет насколько выбранная точка способна изменить модель, если её убрать (зависит от величины рычага и остатка).

Анимация 7. Принцип работы метода расчета рычага и расстояния Кука. Формулы расчета приведены для одной точки, p - количество параметров в модели. При исключении объекта, измеряется метрика ошибки новой модели, если метрика лучше, то выбирается новая модель, если нет - рассматривается альтернативный вариант (ссылка на код для генерации анимации)
Анимация 7. Принцип работы метода расчета рычага и расстояния Кука. Формулы расчета приведены для одной точки, p - количество параметров в модели. При исключении объекта, измеряется метрика ошибки новой модели, если метрика лучше, то выбирается новая модель, если нет - рассматривается альтернативный вариант (ссылка на код для генерации анимации)

В примере который я привел выше для наглядности расчеты производятся итеративно. Однако в Python библиотеке sklearn например, расчет построен таким образом, что расстояние Кука определяется без реальной переоценки модели n раз.

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

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

  • Выбор случайным образом подвыборки размера n объектов

  • Обучение модели по выборке из выбранных n объектов

  • Отсечение выбросов - исключение объектов, для которых ошибка модели превышает заданный порог

  • Опциональный шаг: Обучение модели ещё раз на оставшихся не-выбросах и отсечение выбросов еще раз

  • Подсчет количества не-выбросов (согласованных точек) - m

  • Повторение первых пяти шагов несколько (мы сами определяем количество итераций) раз и последующий выбор той модели, количество не-выбросов m для которой максимально 

Анимация 8. Принцип работы алгоритма RANSAC (ссылка на код для генерации анимации)
Анимация 8. Принцип работы алгоритма RANSAC (ссылка на код для генерации анимации)

Результаты применения алгоритма RANSAC и метода расчета расстояния Кука показаны на Рисунке 34.

Рисунок 34. Модель линейной регрессии построенная при помощи методов фильтрации выбросов RANSAC и расстояния Кука. RMSE эталонной модели на генеральной совокупности - 3873 (ссылка на код для генерации картинки)
Рисунок 34. Модель линейной регрессии построенная при помощи методов фильтрации выбросов RANSAC и расстояния Кука. RMSE эталонной модели на генеральной совокупности - 3873 (ссылка на код для генерации картинки)

Исходя из результатов представленных на Рисунке 34, наиболее перспективной в дуэли оказалась модель построенная при помощи метода RANSAC.

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

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

Усложнение модели. Многомерная (многофакторная) линейная регрессия

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

Рисунок 35. Многомерная линейная регрессия (ссылка на код для генерации основных элементов картинки)
Рисунок 35. Многомерная линейная регрессия (ссылка на код для генерации основных элементов картинки)

Генерация новых признаков

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

Рисунок 36. Полиномиальное преобразование признаков степени 2
Рисунок 36. Полиномиальное преобразование признаков степени 2

Обратите внимание на то, что полученное уравнение теперь - модель полиномиальной регрессии, что позволяет моделировать нелинейные взаимосвязи в данных. Чем выше степень полинома, тем больше степеней свободы модель имеет (Рисунок 37).

Рисунок 37. Примеры полиномов построенных по выборке. Когда становится возможным моделировать нелинейность (ссылка на код для генерации картинки)
Рисунок 37. Примеры полиномов построенных по выборке. Когда становится возможным моделировать нелинейность (ссылка на код для генерации картинки)

Существует большое количество самых разнообразных преобразований, которые можно прикрутить к исходным данным. Однако в случае их применения, модель, получится не совсем линейная (что и можно заметить по форме получившихся модельных кривых на Рисунке 37). Поэтому я не буду останавливаться на них подробно в рамках этой статьи. Если стало любопытно, то можете прочитать подробнее про следующие преобразования которые можно делать над исходными признаками (и как пример литературы на которую можно опереться - Trevor Hastie, Robert Tibshirani, Jerome Friedman - The Elements of Statistical Learning): 

  • Функциональные преобразования

    • Логарифмы: log(x + ε)

    • Обратные значения: 1/x, 1/(x + ε)

    • Корни: \sqrt{x},\ x^{1/3}

    • Экспоненты: \exp(x),\ \exp(-x)

    • Тригонометрические: sin(x), cos(x), tan(x) особенно если признак имеет периодичность

    • Сигмоиды: 1 / (1 + exp(-x))

  • Бинаризация и дискретизация

    • Биннинг: разбить признак X на интервалы, например, [x < 10], [10 ≤ x < 20], [x ≥ 20]

    • Квантильный биннинг: разделение на равные по числу наблюдений группы

  • Пороговые функции (привет, нейросети)

  • Сплайны

  • Вейвлет-преобразования

  • Фурье-признаки

  • и многие другие

Ставьте лайки - если наберется больше 500 плюсов, то обязательно выпущу отдельный пост посвященный генерации признаков, естественно, с большим количеством визуализаций XD

Сбор новых признаков

Если генерация новых признаков не помогает улучшить метрику, то можно подойти к более “тяжелому” подходу - собрать больше данных, но теперь не новых объектов (строк) как мы делали выше, а новых характеристик (столбцов). 

Предположим, что у нас есть возможность собрать некоторое количество дополнительных потенциальных предикторов. В случае со стоимостью квартир имеет смысл рассмотреть: 

  • Площадь квартиры, квадратные метры

  • Удаленность от ближайшей станции метро, м

  • Город

  • Есть ли кондиционер

Обновленный датасет будет выглядеть следующим образом:

Рисунок 38. Датасет D с новыми признаками: площадью, удаленностью от метро, названием города, наличие кондиционера в квартире
Рисунок 38. Датасет D с новыми признаками: площадью, удаленностью от метро, названием города, наличие кондиционера в квартире

К вопросу визуализации 

Оглядываясь назад к Рисунку 1 (или к большинству рисунков встречавшихся в статье раньше) несложно заметить, что мы уже не охватываем все признаки двухмерной картинкой. Поэтому пора сделать новые визуализации, чтобы посмотреть на данные под новым углом (Рисунок 39 и Анимация 9). 

Рисунок 39. Визуализация взаимосвязи нескольких предикторов между собой и откликом. Столбцы и строки - предикторы. По главной диагонали (где предиктор пересекается сам с собой) показаны двумерные диаграммы “признак (ось x) - отклик (ось y)”. Верхний треугольник (то что выше главной диагонали) состоит из 3d графиков “два признака (оси x и y) и отклик (ось z)”. Нижний треугольник, - это немного другая форма представления трехмерных данных в виде контурных карт, где оси - признаки (ссылка на код для генерации картинки) 
Рисунок 39. Визуализация взаимосвязи нескольких предикторов между собой и откликом. Столбцы и строки - предикторы. По главной диагонали (где предиктор пересекается сам с собой) показаны двумерные диаграммы “признак (ось x) - отклик (ось y)”. Верхний треугольник (то что выше главной диагонали) состоит из 3d графиков “два признака (оси x и y) и отклик (ось z)”. Нижний треугольник, - это немного другая форма представления трехмерных данных в виде контурных карт, где оси - признаки (ссылка на код для генерации картинки

Рисунок лучше рассмотреть во всех деталях (Рисунок 40).

Рисунок 40. Предыдущая визуализация (см. Рисунок 39) многомерных данных с комментариями
Рисунок 40. Предыдущая визуализация (см. Рисунок 39) многомерных данных с комментариями
Анимация 9. Трехмерные точечные диаграммы для двух наборов признаков: количество комнат и удаленность от станции метро, площадь квартиры и наличие кондиционера (ссылка на код для генерации анимации)
Анимация 9. Трехмерные точечные диаграммы для двух наборов признаков: количество комнат и удаленность от станции метро, площадь квартиры и наличие кондиционера (ссылка на код для генерации анимации)

Из Анимации 9 видны две выразительных особенности набора данных:

  • Чем ближе квартира к метро, тем выше стоимость; квартиры у метро обычно меньшей площади (наблюдение номер 2 из Рисунка 40);

  • Наличие кондиционера - переменная, по которой заметно разделяется отклик (цена квартиры): когда в квартире установлен кондиционер, квартира обычно стоит дороже (наблюдение номер 6 из Рисунка 40).

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

Вспоминаем про Рисунок 5: почему цена все таки убывала

Вернемся к одному из первых рисунков в заметке (кликайте сюда чтобы перенестись -> figure_5) - там, где объяснялась концепция описания данных прямой линией, был пример из трех объектов, где цена убывала несмотря на увеличение количества комнат. Однако все становится очевидным если мы визуализируем данные с дополнительным признаком:

Дополнительная анимация 4. Объяснение убывание цены квартир несмотря на увеличение количества комнат. Цена растет не из-за уменьшения количества комнат, а потому что квартиры ближе к метро! (ссылка на код для генерации анимации)
Дополнительная анимация 4. Объяснение убывание цены квартир несмотря на увеличение количества комнат. Цена растет не из-за уменьшения количества комнат, а потому что квартиры ближе к метро! (ссылка на код для генерации анимации)

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

Когда количество признаков увеличивается становится сложно строить такие попарные визуализации, которые были приведены на Рисунках 39 и 40. И если в вашем датасете много численных признаков, то довольно популярно отрисовывать матрицы корреляций (Рисунок 41) - их, я уверен, вы будете встречать очень часто, если продолжите интересоваться темой обработки данных.

Рисунок 41. Матрица с численными признаками и рассчитанными коэффициентами корреляции (ссылка на код для генерации картинки)
Рисунок 41. Матрица с численными признаками и рассчитанными коэффициентами корреляции (ссылка на код для генерации картинки)

Тут тот же принцип, что и с оценкой моделирования, - когнитивная нагрузка на инженера меньше при интерпретации чисел (по одному на пару) вместо разглядывания подграфиков. Итак, на Рисунке 41 видно, что цена положительно скоррелирована с признаками “Количество комнат” и “Площадь” и отрицательно с “Расстоянием до метро”, что логично, ведь обычно чем квартира ближе к метро и чем больше её площадь, тем она дороже. 

Отдельно хочется отметить почему часто визуализируется коэффициент корреляции, - всегда полезно смотреть есть ли в датасете сильно скоррелированные друг с другом предикторы (явление мультиколлинеарности). Именно это мы и видем для пары “Количество комнат” -  “Площадь”, где коэффициент пары равен единице. В таких случаях имеет смысл исключить один из признаков, поскольку он не несет полезной информации в модель в то время как ресурсы (например, вычислительные при подготовке данных и при оптимизации модели) забирает. Мультиколлинеарность приводит и к другим неприятным последствиям но поговорим об этом чуть позже.

Про важность предобработки (категориальных) предикторов

Как видно из Рисунка 39, теперь в таблице не только ровные числа количества комнат, но и не очень аккуратные дистанции до метро, так и совсем непонятные названия городов или ответ на вопрос “есть ли что-то в квартире” в виде текста.  

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

Начнем с признака попроще - “Есть ли кондиционер?”, поскольку он принимает только два значения - “да” или “нет”. Такие признаки принято кодировать (переводить из строкового представления) двумя числами, например (Рисунок 42): 

  • да - 1

  • нет - 0

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

Обратим внимание на то, что на Рисунке 42 не две разные модели, построенные каждая для своей подвыборки, а одна. Здесь коэффициент наклона b1 фиксирован а сдвиг модельной линии будет отличаться для объектов когда бинарный признак равен 0 и когда равен 1. Так происходит потому что когда признак равен нулю, будет зануляться и произведение значения предиктора и коэффициента. Это работает хорошо, когда зависимость между признаками и откликом линейна и одинаково направлена для всех объектов выборки, но бинарный признак не поможет при моделировании более сложных разнонаправленных зависимостей (Рисунок 43).

Рисунок 43. Различия зависимостей между признаками и откликом в разных подвыборках приводит к тому, что одна модель с бинарным признаком не способно адекватно моделировать ни одну из частей датасета (ссылка на код для генерации основных элементов картинки)
Рисунок 43. Различия зависимостей между признаками и откликом в разных подвыборках приводит к тому, что одна модель с бинарным признаком не способно адекватно моделировать ни одну из частей датасета (ссылка на код для генерации основных элементов картинки)

Как видно из рисунка 43, в худшем случае модель с бинарным признаком сходится к модели с одним численным признаком. Чтобы эту “проблему” разрешить, можно вспомнить приемы из предыдущей главы, и сгенерировать новый признак - взаимодействие, или мы можем построить две разные модели для разных частей датасета (Рисунок 44).

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

Теперь когда с бинарным признаком разобрались, имеет смысл поговорить про более сложный случай, когда в столбце закодировано более двух уникальных значений. Существует большое количество методов кодировки категориальных значений, - некоторые из них показаны на Рисунке 45 [Unidata Blog, Encoding Categorical Variables: One-Hot vs. Label Encoding and Beyond. 2026. https://unidata.pro/blog/encoding-categorical-variables-one-hot-vs-label/], [Hancock, J. T., et al. Survey on categorical data for neural networks. 2020. https://doi.org/10.1186/s40537-020-00305-w]. Однако я не буду приводить все из них здесь поскольку из своего личного опыта OneHotEncoding было достаточно для решения задач. Просто будьте в курсе что существуют разные механизмы кодировки.

Рисунок 45. Методы кодировки категориальных признаков (см. код на Python для произведения трансформаций над данными)
Рисунок 45. Методы кодировки категориальных признаков (см. код на Python для произведения трансформаций над данными)

Измерение важности признаков 

Теперь когда мы научились усложнять модель путем добавления новых признаков, имеет смысл поговорить о том, как комбинировать независимые переменные, оставлять только полезные. Безусловно, при раздувании признакового пространства путем генерации признаков или через сбор новых данных срабатывают такие ограничители как “здравый смысл” и “модель долго обучается”. Однако мы вполне можем положиться и на более эффективные эвристики, чтобы определить какие признаки имеет смысл в модели оставить. Начнем с самого простого, - рассмотрим поближе коэффициенты модели многомерной регрессии (Рисунок 46).

Рисунок 46. Размер коэффициента как индикатор важности признака (ссылка на код для генерации основных элементов картинки)
Рисунок 46. Размер коэффициента как индикатор важности признака (ссылка на код для генерации основных элементов картинки)

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

Нормализация и стандартизация предикторов

Нормализация, это трансформация данных, которая приводит значения массивов к одинаковому диапазону (Рисунок 47).

Рисунок 47. Демонстрация результатов применения алгоритмов нормализации данных на примере двух признаков: количество комнат и расстояние до метро (ссылка на код для генерации картинки) 
Рисунок 47. Демонстрация результатов применения алгоритмов нормализации данных на примере двух признаков: количество комнат и расстояние до метро (ссылка на код для генерации картинки

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

Какие именно формулы используются для нормализации и стандартизации, показаны на Рисунке 48.

Рисунок 48. Способы нормализации. Приведены экстремальные случаи с выбросами, тем не менее в случае репрезентативной обучающей выборки появление подобных выбросов будет сведено к минимуму (см. код на Python для произведения нормализации)
Рисунок 48. Способы нормализации. Приведены экстремальные случаи с выбросами, тем не менее в случае репрезентативной обучающей выборки появление подобных выбросов будет сведено к минимуму (см. код на Python для произведения нормализации)

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

Пространство коэффициентов модели когда признаки стандартизированы

После того как исходные признаки стандартизированы, а значит коэффициенты b_0, b_1, b_2 и т.д. “выровнены” (так удобнее варьировать коэффициенты), самое время более пристально посмотреть на то как значения этих коэффициентов влияют на ошибку модели. Для измерения ошибки воспользуемся метриками MAE и MAPE для парной регрессии и RMSE для многомерной.

Анимация 10. Взаимосвязь коэффициентов ,  парной линейной регрессии и метрики MAE. Предиктор в модели - количество комнат. Пояснение про изменяющийся коэффициент сдвига в исходных единицах: мы меняем наклон в стандартизованных данных, а при переводе обратно сдвиг автоматически пересчитывается (ссылка на код для генерации анимации)
Анимация 10. Взаимосвязь коэффициентов b_0, b_1 парной линейной регрессии и метрики MAE. Предиктор в модели - количество комнат. Пояснение про изменяющийся коэффициент сдвига в исходных единицах: мы меняем наклон в стандартизованных данных, а при переводе обратно сдвиг автоматически пересчитывается (ссылка на код для генерации анимации)

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

Для наглядности, давайте посмотрим на то как по разному будут выглядеть ландшафты метрик, например для средней абсолютной процентной ошибки мы увидим следующее - поскольку MAPE чувствительна к ошибкам с небольшим значением отклика (“дешевые квартиры”), то минимум “растягивается” в вытянутую долину. И потому многие комбинации дают похожий MAPE, если линия хорошо попадает в область маленьких y, даже если на дорогих квартирах она заметно ошибается. (Анимация 11).

Анимация 11. Взаимосвязь коэффициентов ,  парной линейной регрессии и метрики MAPE. Предиктор в модели - количество комнат (ссылка на код для генерации анимации)
Анимация 11. Взаимосвязь коэффициентов b_0, b_1 парной линейной регрессии и метрики MAPE. Предиктор в модели - количество комнат (ссылка на код для генерации анимации)

Далее увеличим количество признаков в модели, чтобы было необходимо найти оптимальное сочетание уже не двух, а трех коэффициентов (Анимации 12 и 13). 

Анимация 12. Взаимосвязь коэффициентов , ,  и метрики RMSE. Предикторы в модели - количество комнат (), удаленность от метро () (ссылка на код для генерации анимации)
Анимация 12. Взаимосвязь коэффициентов b_0, b_1, b_2 и метрики RMSE. Предикторы в модели - количество комнат (x_1), удаленность от метро (x_2) (ссылка на код для генерации анимации)
Анимация 13. Взаимосвязь коэффициентов , ,  и метрики RMSE. Предикторы в модели - количество комнат (), площадь квартиры () (ссылка на код для генерации анимации)
Анимация 13. Взаимосвязь коэффициентов b_0, b_1, b_2 и метрики RMSE. Предикторы в модели - количество комнат (x_1), площадь квартиры (x_2) (ссылка на код для генерации анимации)

На анимациях выше видно, что признаки имеют сильную линейную связь, например в Анимации 12, на проекции коэффициентов b1 vs b2 (плоскость слева на нижнем левом графике) виден линейный паттерн. Это говорит о том что 1) присутствует сильная обратная корреляция между признаками “количество комнат” и “расстояние до метро” и 2) несмотря на то что коэффициенты “плавают вдоль долины низких значений RMSE”, предсказания модели стабильны, - ошибка модели практически не изменяется. Значит, в том числе и о то что признаки дают похожую информацию. Такая же ситуация и с Анимацией 13, но здесь линейная связь между признаками еще сильнее, и уже прямая, а не обратная.

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

Обобщение аналитического решения на многомерный случай

Выше по тексту когда мы исследовали ландшафт функционала ошибки, мы визуально наблюдали локации где значение ошибки модели было минимальным. Модель такой информации не имеет и находит оптимум (оптимальное сочетание коэффициентов b_0, b_1, b_2, и т.д.) по формуле. Для парной регрессии, когда признак один, формулу мы привели (Рисунок 6). Однако теперь признаков несколько и, когда они предобработаны, самое время задаться вопросом, - как найти оптимальные коэффициенты для многомерной линейной регрессии, как можно обобщить решение для данных большей размерности.

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

Когда мы рассматривали формулу аналитического решения ранее, мы приводили расчеты в скалярном виде. Однако намного эффективнее перейти к векторной записи, для этого визуализируем исходные данные не в пространстве признаков, а в пространстве наблюдений (Рисунок 49).

Рисунок 49. Игрушечный исходный датасет и данные в пространстве наблюдений (объектов)
Рисунок 49. Игрушечный исходный датасет и данные в пространстве наблюдений (объектов)

Несмотря на то что поначалу такой взгляд на данные может казаться контринтуитивным, за ним не стоит никакой магии, данные все те же, изменилась лишь форма. А мы продолжаем, - из школьной программы нас учили (по крайней мере меня) называть вектором направленные отрезки. Эти “направленные отрезки” можно умножать на число и складывать друг с другом. Задача алгоритма линейной регрессии в векторном пространстве - найти такое преобразование вектора x, чтобы получившийся вектор предсказаний \hat{y} (обозначется как y с крышечкой наверху) был как можно больше похож на вектор откликов y. Давайте попробуем перебрать такие преобразования, начиная с умножения на число (Рисунок 50). 

Рисунок 50. Строим простейшую линейную регрессию: только угловой коэффициент , умножение вектора x на число
Рисунок 50. Строим простейшую линейную регрессию: только угловой коэффициент b_1, умножение вектора x на число

Изучаем Рисунок 50, начиная с левого верхнего угла: модель никак признак x не преобразовывает (коэффициент b_1 равен единице) и потому значение предсказания равно значению признака, вектор x полностью совпадает с вектором предсказаний. Если же коэффициент b_1 больше 1, умножение вектора x на этот коэффициент приводит к тому, что длина вектора предсказаний пропорционально увеличивается. Также вектор признака можно сжимать (коэффициент b_1 от 0 до 1) и направлять в противоположную сторону (коэффициент b_1 меньше 0).

Рисунок 51. Что делать если умножения на коэффициент  недостаточно
Рисунок 51. Что делать если умножения на коэффициент b_1 недостаточно

Если из Рисунка 50 хорошо видно как производится умножение вектора на число, то вот на Рисунке 51 были замечены еще две операции над вектором, напомнить про которые имеет смысл отдельно (Рисунок 52).

Рисунок 52. Маленькое но важное напоминание: параллельный перенос и сложение векторов
Рисунок 52. Маленькое но важное напоминание: параллельный перенос и сложение векторов

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

Анимация 14. Увеличиваем размер выборки до трех объектов (наблюдений). Попробуйте мысленно провести прямую на графике слева так, чтобы она пересекала все три точки (ссылка на код для генерации анимации)
Анимация 14. Увеличиваем размер выборки до трех объектов (наблюдений). Попробуйте мысленно провести прямую на графике слева так, чтобы она пересекала все три точки (ссылка на код для генерации анимации)

Чем больше объектов, тем больше размерность, тем больше осей будет на графике. Нам, как людям, такое представить трудно, и потому в более высокие размерности далее по тексту мы не пойдем, - нет необходимости. Тем более что принципы, которые мы здесь обсуждаем будут работать и там. В частности, задача всё та же - нам нужно найти такую комбинацию векторов v (состоит из единиц) и x (вектор-признак из датасета) чтобы получившийся вектор предсказаний был как можно ближе к вектору откликов y. Все что мы можем варьировать здесь, это изменять коэффициенты, на которые мы умножаем вектор v (b_0) и вектор x (b_1). Давайте попробуем перебрать разные комбинации и посмотрим на то как будет выглядеть решение в пространстве признаков и в векторном пространстве (Анимация 15).

Анимация 15. Исследуем коэффициенты модели парной линейной регрессии для трех объектов: визуализация векторов отклика и предсказаний, составленного из векторов признаков  и . Визуализация подпространства  (ссылка на код для генерации анимации)
Анимация 15. Исследуем коэффициенты модели парной линейной регрессии для трех объектов: визуализация векторов отклика и предсказаний, составленного из векторов признаков v и x. Визуализация подпространства Col(X) (ссылка на код для генерации анимации)

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

Мы собрали это подпространство в визуализации выше путем сложения векторов v и x с разными коэффициентами перед ними. Однако существует и другая запись того же самого выражения, но более компактная - через матричное умножение. Для этого составим еще один вектор, но на этот раз, из коэффициентов b_0 и b_1. Обозначим его как вектор B. Его можно преобразовать путем умножения на матрицу, чтобы повернуть вектор, растянуть или сжать его, а также отобразить в другое (под)пространство. И если взять матрицу X составленную из векторов столбцов v и x и умножить на вектор \vec{b}, составленный из коэффициентов, то мы получим отображение y на подпространство предсказаний Col(X) (Рисунок 53).

Рисунок 53. Преобразование вектора откликов  в вектор предсказаний
Рисунок 53. Преобразование вектора откликов y в вектор предсказаний

Обратим внимание на то, что в соответствии с нашими предположениями, вектор отклика не лежит в подпространстве предсказаний: если между двумя точками можно провести прямую, точно через них проходящую, то в случае с тремя и более возрастает вероятность того, что идеальной модели с нулевой ошибкой не существует. Именно поэтому вектор отклика не лежит на гиперплоскости даже для оптимальной модели (см. Черный вектор модели C на Рисунке 54).

Рисунок 54. Визуализация двух неудачных моделей (A и B) и одной оптимальной (C) (ссылка на код для генерации картинки)
Рисунок 54. Визуализация двух неудачных моделей (A и B) и одной оптимальной (C) (ссылка на код для генерации картинки)

Присмотримся к картинке и попробуем найти отличия между векторами предсказаний разных моделей A, B и C: вектор модели C выглядит как тень вектора отклика на плоскости. Тогда процесс решения задачи линейной регрессии можно интерпретировать как проецирование вектора y на подпространство Col(X). Лучший прогноз из всех возможных - вектор заканчивающийся в ближайшей к отклику точка на плоскости. А ближайшая точка на плоскости, как можно вспомнить из базовых геометрических знаний, получается в основании перпендикуляра. Получившаяся вертикаль, это тоже вектор и его называет вектором остатков e, потому что он получается путем вычитания предсказаний из отклика (вспоминаем формулу остатков из главы про оценку качества модели с помощью визуализаций).

Итак, нам известны вектор откликов y и вектор признаков x. Мы хотим найти такой вектор коэффициентов \vec{b} чтобы получившийся вектор предсказаний был как можно ближе к y. Мы не знаем вектор остатков e, но знаем что этот вектор ортогонален пространству Col(X). А это в свою очередь значит, что вектор e будет ортогонален каждому направлению в плоскости - значит, в частности, перпендикулярен каждому столбцу матрицы X (вектору v и вектору x).

Рисунок 55. Используем свойство ортогональности для вывода формулы. Чтобы найти вектор коэффициентов нам надо транспонировать, перемножить и обратить матрицу признаков. Метод наименьших квадратов
Рисунок 55. Используем свойство ортогональности для вывода формулы. Чтобы найти вектор коэффициентов нам надо транспонировать, перемножить и обратить матрицу признаков. Метод наименьших квадратов

Название аналитического метода который мы только что рассмотрели - метод наименьших квадратов (МНК, или Ordinary Least Squares, OLS, по английски). Название такое, потому что мы подбирали такие коэффициенты, чтобы минимизировать сумму квадратов остатков модели (Рисунок 6, кликайте чтобы перенестись -> figure_6). В векторном пространстве размер остатков - это квадрат евклидова расстояния от точки отлика до подпространства Col(X) (Рисунок 55). Наименьшие квадраты равно “наименьшее расстояние в квадрате”.

Вспомним цель текущего раздела - мы продирались через формулы и визуализации выше чтобы обобщить аналитическое решение на многомерный случай. И теперь самое время проверить как будет работать формула для случая когда признаков не один, а два! Рассмотрим на примере датасета с тремя объектами в который добавим один признак (Анимация 16).

Анимация 16. Что если количество признаков возрастает, - многомерная регрессия в векторном виде. Формула все та же, только в матрицу  добавился новый вектор (). Подпространство  для наглядности ограничено многоугольником (ссылка на код для генерации анимации)
Анимация 16. Что если количество признаков возрастает, - многомерная регрессия в векторном виде. Формула все та же, только в матрицу X добавился новый вектор (x_2). Подпространство Col(X) для наглядности ограничено многоугольником (ссылка на код для генерации анимации)

На анимации 16 прокомментируем три важных наблюдения:

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

  2. Смотрим направо - количество векторов не изменилось, поскольку в выборке как и прежде, три объекта

  3. Теперь подпространство Col(X) не плоскость на графике, а все пространство. Для наглядности значения были ограничены объемной фигурой, - параллелепипедом. И поскольку это подпространство полностью покрывает собой вектор отклик y, проекция отклика будет тривиальна. На Анимации это выглядит так: вектор отклика и вектор предсказаний совпадают. Остаток нулевой.

Когда аналитическое решение встречает трудности 

Представим, что нам не повезло, и новый признак x_2 новой информации не привносит. Пусть новый признак линейно выражается через два другие (сдвиг, и признак x_1) и тогда многоугольник Col(X) опять сжимается до состояния плоскости - Анимация 17.

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

И несмотря на то что проекцию в таком подпространствы мы ранее находили без всяких проблем, теперь вектор предсказаний составляется не из двух векторов (сдвиг и x_1), а из трех (сдвиг, x_1 и x_2). И поскольку степеней свободы больше, существует больше одного решения: на графике слева это визуализировано через две отдельные модельные поверхности которые с точки зрения метода наименьших квадратов “одинаково хорошо” описывают данные. На картинке справа показаны векторы признаков каждой из моделей, складывающиеся в один вектор предсказаний.

При таких исходных данных в расчетах проблема возникнет на этапе расчета обратной матрицы (Рисунок 56).

Рисунок 56. Невозможность применить формулу для аналитического решения, которую использовали ранее. Точно такая же проблема возникнет и на нашем основном датасете с ценами квартир!
Рисунок 56. Невозможность применить формулу для аналитического решения, которую использовали ранее. Точно такая же проблема возникнет и на нашем основном датасете с ценами квартир!

Как видно из Рисунка 56, матрица вырождена, а значит формула с обратной матрицей неприменима и единственного решения нет. Стоить отметить, что даже если точной линейной зависимости нет, но признаки сильно коррелируют друг с другом (например, площадь и количество комнат), проблема остается: матрица становится плохо обусловленной, и решение получается численно неустойчивым. Могут присутствовать и другие проблемы, например с one-hot закодированными признаками, но уже и этого достаточно чтобы задуматься об альтернативных методах решения.

Помимо приведенного выше, решение линейной регрессии в аналитическом виде невозможно, если:

  1. Используется не-квадратичная или негладкая функция потерь (L1, квантильная и т.д.) - тогда задача уже не сводится к методу наименьших квадратов;

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

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

Численные решения

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

Оптимизация полным перебором

Самый простой способ, - исследовать все комбинации коэффициентов с некоторым шагом. Метод полного перебора (exhaustive search) в данном случае - это перебор всех пар коэффициентов из заранее заданной дискретной сетки значений и выбор пары с минимальной ошибкой. Обычно выбирают метрику MSE (как RMSE но без квадратного корня) для измерения ошибки.

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

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

Анимация 18. Метод полного перебора для нахождения решения парной регрессии (ссылка на код для генерации анимации)
Анимация 18. Метод полного перебора для нахождения решения парной регрессии (ссылка на код для генерации анимации)

Оптимизация случайным поиском

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

Анимация 19. Случайный поиск для нахождения оптимального набора коэффициентов парной регрессии (ссылка на код для генерации анимации)
Анимация 19. Случайный поиск для нахождения оптимального набора коэффициентов парной регрессии (ссылка на код для генерации анимации)

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

Использование информации о направлении

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

Для наглядности упростим задачу - возьмем срез по одному параметру, одномерный массив, и рассмотрим на его примере. При движении по массиву воспользуемся тем, что мы на предыдущем шаге уже рассчитали значение ошибки (возьмем MSE для данного примера) и сравнивая текущее и предыдущее сделаем вывод куда имеет смысл сделать следующий шаг - Рисунок 57.

Рисунок 57. Спуск с использованием попарных сравнений. Оптимизируем значения коэффициента в срезе по параметру сдвига  (ссылка на код для генерации картинки)
Рисунок 57. Спуск с использованием попарных сравнений. Оптимизируем значения коэффициента в срезе по параметру сдвига b_0 (ссылка на код для генерации картинки)

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

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

Анимация 20. Анимация спуска с использованием попарных сравнений по параболе. Показаны примеры для двух начальных приближений (желтого и зеленого) (ссылка на код для генерации анимации)
Анимация 20. Анимация спуска с использованием попарных сравнений по параболе. Показаны примеры для двух начальных приближений (желтого и зеленого) (ссылка на код для генерации анимации)

Явно отметитим, что в Анимации 20 шаг всегда равен одному интервалу (шагу сетки) и мы пока не используем никаких производных (предвкушая алгоритм градиентного спуска) - просто попарно сравниваем метрики.  

У описанного выше подхода есть один существенный недостаток, - он сильно зависит от размеров сетки. Так например, если сетка будет подробная, то придется сделать очень много шагов в направлении оптимума. И наоборот, если сетка слишком грубая - пропустим оптимум (Анимация 21).

Анимация 21. Спуск с использованием попарных сравнений. Скорость сходимости и размер шага (ссылка на код для генерации анимации)
Анимация 21. Спуск с использованием попарных сравнений. Скорость сходимости и размер шага (ссылка на код для генерации анимации)

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

Оптимизация методом градиентного спуска 

При уменьшении шага сетки при попарных сравнениях мы приходим к определению производной через предел (Рисунок 58).

Рисунок 58. Градиент на срезе функционала ошибки: в одномерном случае он совпадает с производной и задаёт направление изменения MSE (ссылка на код для генерации основных элементов картинки)
Рисунок 58. Градиент на срезе функционала ошибки: в одномерном случае он совпадает с производной и задаёт направление изменения MSE (ссылка на код для генерации основных элементов картинки)

Предлагаю ниже посерфить по ландшафту ошибки - рассмотрим анимацию где показаны величины градиента и антиградиента (Анимация 22). Как мы видим, шаг мы теперь можем выбирать произвольный - регулярной сеткой теперь не ограничены [Goh, Gabriel. Why Momentum Really Works. 2017. https://distill.pub/2017/momentum/].

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

При обобщении на многомерные пространства, например когда одновременно оптимизируем коэффициенты сдвига и наклона, градиент будет состоять из частных производных (Рисунок 59).

Рисунок 59. Выбор направления когда градиент рассчитывается с учетом двух коэффициентов (ссылка на код для генерации картинки)
Рисунок 59. Выбор направления когда градиент рассчитывается с учетом двух коэффициентов (ссылка на код для генерации картинки)

Самое время посмотреть на градиентный спуск в действии (Анимация 23).

Анимация 23. Градиентный спуск для нахождения оптимального набора коэффициентов парной регрессии. Как правило, начальное приближение выбирается в координатах 0, 0 или близких к этому, однако далее по тексту чтобы добавить вариативности в визуализации я буду задавать начальное приближение в разных местах (ссылка на код для генерации анимации)
Анимация 23. Градиентный спуск для нахождения оптимального набора коэффициентов парной регрессии. Как правило, начальное приближение выбирается в координатах 0, 0 или близких к этому, однако далее по тексту чтобы добавить вариативности в визуализации я буду задавать начальное приближение в разных местах (ссылка на код для генерации анимации)
Посмотрите как сходится градиентный спуск при разных значениях темпа обучения (learning rate)
Дополнительная анимация 5. Медленно продвигаемся в сторону оптимума с темпом обучения 0.06. Максимально разрешенное количество итераций - 25 (ссылка на код для генерации анимации)
Дополнительная анимация 5. Медленно продвигаемся в сторону оптимума с темпом обучения 0.06. Максимально разрешенное количество итераций - 25 (ссылка на код для генерации анимации)
Дополнительная анимация 6. Перепрыгиваем оптимум с темпом обучения 3.0 (ссылка на код для генерации анимации)
Дополнительная анимация 6. Перепрыгиваем оптимум с темпом обучения 3.0 (ссылка на код для генерации анимации)

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

Рисунок 60. Модель можно оптимизировать по разному. Функционал ошибки Тьюки как способ борьбы с выбросами (ссылка на код для генерации основных элементов картинки)
Рисунок 60. Модель можно оптимизировать по разному. Функционал ошибки Тьюки как способ борьбы с выбросами (ссылка на код для генерации основных элементов картинки)

В случае использования функционал ошибки Тьюки, процесс оптимизации будет выглядеть следующим образом - Анимация 24.

Анимация 24. Замена функционала ошибки MSE на функцию потерь Тьюки (ссылка на код для генерации анимации)
Анимация 24. Замена функционала ошибки MSE на функцию потерь Тьюки (ссылка на код для генерации анимации)

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

Анимация 25. Градиентный спуск - метод локальной оптимизации, где важно начальное приближение. Показано на примере функции потерь Тьюки (ссылка на код для генерации анимации)
Анимация 25. Градиентный спуск - метод локальной оптимизации, где важно начальное приближение. Показано на примере функции потерь Тьюки (ссылка на код для генерации анимации)

Переходим к многомерной регрессии. Если рассмотреть историю сходимости решения к оптимальным коэффициентам, то можно увидеть как постепенно коэффициенты перед “важными” признаками растут и как ошибка также постепенно падает (Рисунок 61).

Рисунок 61. Процесс схождения к оптимальному решению модели множественной линейной регрессии (ссылка на код для генерации картинки)
Рисунок 61. Процесс схождения к оптимальному решению модели множественной линейной регрессии (ссылка на код для генерации картинки)

Регуляризация 

Вспомним эффект из Анимации 5, где разные обучающие выборки приводили к разным оцененным коэффициентам несмотря на то что мы пытались восстановить одну и ту же исходную зависимость между признаком и откликом. Модель получалась нестабильна - чувствительна к разбиению. 

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

Поэтому в текущем разделе будем рассматривать задачу оценки коэффициентов в двух плоскостях:

  • Как регуляризация помогает при неоднородных разбиениях

  • Как регуляризация помогает модели хорошо обобщать данные

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

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

Один из способов обеспечения численной стабилизации, это наложение ограничений на коэффициенты - регуляризация (Рисунок 63).

Рисунок 63. Наложение ограничений на значения коэффициентов перед признаками модели линейной регрессии. Lasso и Ridge регрессии. Разбиение номер 2
Рисунок 63. Наложение ограничений на значения коэффициентов перед признаками модели линейной регрессии. Lasso и Ridge регрессии. Разбиение номер 2

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

Рисунок 64. Процесс сходимости коэффициентов при регуляризации L1 (Lasso) и L2 (Ridge). Разбиение номер 2 (ссылка на код для генерации картинки)
Рисунок 64. Процесс сходимости коэффициентов при регуляризации L1 (Lasso) и L2 (Ridge). Разбиение номер 2 (ссылка на код для генерации картинки)

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

А теперь изучим как выглядят графики с оцененными коэффициентами, если к моделям на разных разбиениях была применена регуляризация на Рисунке 65.

Рисунок 65. Коэффициенты модели многомерной линейной регрессии полученные при помощи гребневой линейной регрессии, в сравнении с коэффициентами без регуляризации (ссылка на код для генерации картинки)
Рисунок 65. Коэффициенты модели многомерной линейной регрессии полученные при помощи гребневой линейной регрессии, в сравнении с коэффициентами без регуляризации (ссылка на код для генерации картинки)

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

Переобучение

И раз мы заговорили о регуляризации, нельзя не упомянуть что степень регуляризации можно варьировать (Анимация 26).

Анимация 26. Диаграмма рассеяния предсказания vs действительные значения и значения метрик моделей, полученных при помощи регуляризации моделей (ссылка на код для генерации анимации)
Анимация 26. Диаграмма рассеяния предсказания vs действительные значения и значения метрик моделей, полученных при помощи регуляризации моделей (ссылка на код для генерации анимации)

Из Анимации 26 видно:

  • Линия 1: Показаны значения коэффициентов перед признаками, метрики на обучении и тесте, а также график сравнения предсказаний и фактических значения для модели без регуляризации;

  • Линия 2: Как себя ведет Lasso регрессия при разной силе регуляризации, ошибка сначала уменьшается на тестовых данных, но затем модель деградирует к предсказанию средним из-за чрезмерной регуляризации и значения коэффициентов перед признаками зануляются;

  • Линия 3: Гребневая регрессия при увеличении силы регуляризации показывает всё меньшие и меньшие значения ошибки на тестовой выборке несмотря на то, что ошибка на обучении постепенно увеличивается.

Итог из анимации 26, - при слабой регуляризации, модель показывает себя отлично на обучающей выборке, но на тесте качество значительно проседает. Это хороший пример переобучения (Рисунок 66).

Рисунок 66. Явление переобучения - когда модель плохо работает на новых данных
Рисунок 66. Явление переобучения - когда модель плохо работает на новых данных

Приведем искусственный но очень наглядный пример на сгенерированных признаках для полиномиальной регрессии.

Анимация 27. Регуляризация с полиномиальными признаками. Когда модель учится находить важные закономерности, а не пытается выучить шум в данных. Данные сгенерированы синтетически, где исходная зависимость линейная, при этом на обучающую выборку наложен шум, а на тестовую - нет (ссылка на код для генерации анимации)
Анимация 27. Регуляризация с полиномиальными признаками. Когда модель учится находить важные закономерности, а не пытается выучить шум в данных. Данные сгенерированы синтетически, где исходная зависимость линейная, при этом на обучающую выборку наложен шум, а на тестовую - нет (ссылка на код для генерации анимации)

Настройка гиперпараметров

Выше мы затронули очень важную тему, - как понять какое значение гиперпараметра альфа подходит нашему датасету. Мы можем разбить выборку на обучение и тест и попробовать обучить n моделей на обучающей выборке с последующим измерением метрики на тестовых для каждой модели. И выберем ту, где ошибка на тесте будет минимальной (Рисунок 67). 

Рисунок 67. Концепт подбора гиперпараметров поиском по сетке с измерением метрик на тестовой выборке для того чтобы найти оптимальное значение (ссылка на код для генерации основных элементов картинки) 
Рисунок 67. Концепт подбора гиперпараметров поиском по сетке с измерением метрик на тестовой выборке для того чтобы найти оптимальное значение (ссылка на код для генерации основных элементов картинки

Однако так мы рискуем настроить модель под конкретную тестовую выборку, поэтому в машинном обучении принято использовать кросс валидацию (Рисунок 68).

Рисунок 68. Разделение на обучение, тест и валидацию, обучение модели на данных (ссылка на код для генерации основных элементов картинки)
Рисунок 68. Разделение на обучение, тест и валидацию, обучение модели на данных (ссылка на код для генерации основных элементов картинки)

Линейная регрессия - это целый мир

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

В статистике и теории вероятностей, это оценка параметров, распределения остатков,  предсказательные интервалы, статистическое тестирование.

В линейной алгебре это вектора, матричные операции, проекции в подпространства признаков и многое другое.

Анимация 28. Спасибо за внимание!
Анимация 28. Спасибо за внимание!

Заключение

Комплимент дня Спасибо всем тем кто дочитал до этого момента.

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

Если вам понравились визуализации и примеры, и вы хотите использовать их в своих лекциях или выступлениях, пожалуйста делайте это. Все материалы и исходный код для их генерации есть в репозитории на GitHub - https://github.com/Dreamlone/linear-regression 

Искренне Ваш, Михаил Сарафанов

Благодарности

Выражаю искреннюю благодарность другу, Клюшникову Александру, за рецензирование данной статьи :)