Знакомьтесь, линейные модели10.03.2016 16:03

Машинное обучение шагает по планете. Искусственный интеллект, поскрипывая нейронными сетями, постепенно опережает людей в тех задачах, до которых успел дотянуться своими нейронами. Однако не стоит забывать и про простую модель линейной регрессии. Во-первых, потому что на ней построены многие сложные методы машинного обучения, включая нейронные сети. А, во-вторых, потому что зачастую прикладные бизнес-задачи легко, быстро и качественно решаются именно линейными моделями.
И для начала небольшой тест. Можно ли с помощью линейной модели описать:
— зависимость веса человека от его роста?
— длительность ожидания в очереди в магазине в разное время суток?
— посещаемость сайта в фазе экспоненциального роста?
— динамику во времени количества человек, ожидающих поезда на станции метро?
— вероятность, что клиент не оформит заказ на сайте в зависимости от его производительности?
Как вы догадываетесь, на все вопросы ответ будет «Да, можно». Так что линейные модели не так просты, как может показаться на первый взгляд. Поэтому давайте познакомимся с их богатым разнообразием.

Постановка задачи

Пусть дан набор из M пар исходных данных (Xi, Yi), где
ℝ — множество вещественных чисел
ℕ — множество натуральных чисел
M Є ℕ, i Є ℕ, 1 ≤ i ≤ M
Xi Є ℝd, d Є ℕ, т.е. каждый Xi — это последовательность вещественных чисел длиной d
Yi Є ℝk, k Є ℕ, т.е. и каждый Yi — тоже последовательность вещественных чисел, но длиной k.
Х называют независимыми переменными (и также факторами или регрессорами). А Y называют зависимыми или объясняемыми переменными.

В самом простом случае d = k = 1, то есть нам заданы пары чисел, например, (рост человека, вес человека) или (объем детали, вес детали) или (время ожидания в очереди, признак ухода клиента).
Чаще всего d > 1, а k = 1, то есть для каждого Y даны несколько Х, например, ((рост, возраст), вес) или ((продажи вчера, продажи позавчера, продажи позапозавчера), продажи сегодня) или ((количество поставленных задач, суммарная сложность задач), количество прогулянных работником дней под благовидным предлогом).
Бывает, что и Y тоже состоит из нескольких значений и тогда пары (Xi, Yi) могут выглядеть, например, вот так — ((длительность физической нагрузки, уровень нагрузки, индекс массы тела, уровень тренированности), (уровень лактата в крови, длительность периода восстановления)).

И, как вы уже возможно догадались, мы вдруг делаем неожиданное предположение: наши Y не просто зависят, а линейно зависят от Х, то есть Y = bTX, где b — это параметры модели в виде вещественной матрицы размерности d ₓ k.

Однако сразу же вспомним про два важных аспекта. Во-первых, нет никаких причин, чтобы bTX вот так вот строго равнялось Y. Ведь и исходные данные мы могли измерить (и скорее всего измерили) с погрешностью. Кроме того, возможно (и даже скорее всего почти точно) мы упустили из виду какие-то факторы X, которые тоже влияют на Y. Поэтому в модели надо учесть случайную ошибку.
Во-вторых, нет никакой гарантии, что Y линейно зависят сразу напрямую от X. Ведь возможно они зависят от чего-то другого, что в свою очередь зависит от Х. Например, Y может быть равен b ∙ X2 или b ∙ log X. А тогда почему, собственно, мы жаждем именно зависимости самого Y, а не какого-то значения, зависящего от Y, допустим, log Y = bTX.

Таким образом мы приходим к обобщенной постановке линейной задачи:
Y* = bTX* + ℇ
где X* = f (X), Y* = g (Y),
ℇ — случайная величина.

Несмотря на то, что f и g могут быть нелинейными функциями, и Y в результате может весьма нелинейно зависеть от Х, модель все равно остается линейной относительно параметров b. Именно поэтому она и называется линейной моделью.

Но это еще не вся постановка задачи
Пока мы знаем только X и Y. Правда, нашей заслуги в этом нет, поскольку они нам были даны. А как нам узнать b и ℇ?

Вспомним, что сутью модели как раз являются коэффициенты b. Ведь только ради них мы и затеяли всю эту суету. И если бы не было никаких случайных ошибок ℇ, то b можно было бы легко вычислить, решив элементарную систему из d уравнений, где d — количество факторов, то есть длина вектора Xi. Иными словами, в качестве исходных данных достаточно получить ровно d пар (Xi, Yi), где каждый Xi состоит ровно из d значений — и точная модель готова.
Например, если исходные данные имеют вид ((размер детали, цена материала), себестоимость детали), то достаточно данных всего лишь о двух деталях, чтобы построить точную формулу себестоимости.

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

И все же они существуют
Из-за случайностей мы не можем легко рассчитать искомые коэффициенты b нашей модели. Поэтому придется решить оптимизационную задачу вида b* = argmin F (b | X, Y), где F — некий функционал (например, ∑ (Yi — bTXi)2, что приводит нас к методу наименьших квадратов, но можно задать и другие функционалы).

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

Разберемся с ошибками
С самими случайными ошибками ℇ все еще сложнее. Вроде мы их не знаем, и вычислить нам их не из чего, но все же они (точнее информация о них) нужна, как минимум, для точного формулирования вышеупомянутой оптимизационной задачи (а иначе как ее решать?!). Поэтому придется делать какие-то параметрические предположения.
Допустим, исходя из характера процесса, который генерирует исходные данные, мы априорно заявляем, что ℇ ∽ ℥(θ), то есть наша случайная величина ℇ описывается известным распределением ℥ с набором параметров θ, например, нормальным распределением.
Или, не зная распределения, мы можем высказать предположения о его важных свойствах. Например, что E[ℇ] = 0, то есть математическое ожидание ошибок равно нулю. Или что D[ℇi] = σ2 = const, то есть дисперсии всех ошибок одинаковы и конечны.

Но зачем это всё? Ведь можно просто взять argmin от выбранного функционала и получить значение b.
Действительно можно. Вопрос лишь в качестве этого значения. Если непонятно что взять и посчитать, то получится непонятно что, а не модель.
Например, если ℇ имеет распределение Коши, то решение методом МНК даст хаотический и бессмысленный результат.
А вот, например, если одновременно выполняются условия
— E[ℇ] = 0
— E[ℇ | X] = 0, то есть ошибки не зависят от X
— D[ℇi] = σ2 = const
— cov (ℇi, ℇj) = 0, где i ≠ j
то в результате расчета по МНК мы получим не просто оценку искомых параметров b, а наиболее эффективную, несмещенную и состоятельную оценку. И это все очень четко обосновано, доказано и обвязано теоремами со всех сторон.

Заметим важный факт. Поскольку произведение bTX* абсолютно детерминировано, то можно сказать, что Y* тоже является случайной величиной, которая имеет такую же форму распределения, что и ℇ. И тогда можно сделать удобный вывод, переписав модель в виде E[Y*] = bTX*, что означает, что наша линейная модель предсказывает не само значение Y*, а его математическое ожидание. А дальше мы подключаем для решения нашей задачи богатейший статистический матаппарат.
В частности, имея исходные данные и задав параметрическое предположение относительно формы распределения ℇ (а значит и Y*), можно построить функцию правдоподобия ℒ(b) и затем максимизировать ее, что будет равнозначно построению нашей линейной модели. Иными словами, подобрав параметры распределения (а это все те же b), мы научимся генерировать такие случайные величины, которые будут максимально похожи на Y. Чего мы в общем-то и добивались.

Многообразие линейных моделей

Повторим описание нашей модели:
Y* = bTX* + ℇ
b* = argmin F (b | X*, Y*).
Меняя формат и устройство компонентов (Y*, X*, b*, ℇ, F), будем получать разные модели, которые обладают отличающимися свойствами и применимы в различных задачах.

По размерности независимой переменной (X) можно выделить однофакторные (univariable) и многофакторые (multivariable) модели.
Если зависимая переменная (Y) является скалярным значением, то имеем одинарную (univariate) модель. Когда зависимая переменная является многомерной, то есть представлена вектором, получаем множественную или общую (multivariate, general) модель.

Также не будем забывать, что независимые переменные могут содержать как исходные данные, так и их преобразования, в том числе неоднократные. Например, пусть на входе имеем единственную переменную Х, тогда из нее можно сделать несколько факторов — [X, X2, X3] — и получим тем самым, как бы это странно вместе ни звучало, полиномиальную линейную модель.
Еще одним отличные примером является преобразования категорийных переменных. Например, одна из переменных в исходных данных принимает значения «метро», «автомобиль», «велосипед». Из нее создается сразу три фактора для модели: Х*метро, Х*автомобиль, Х*велосипед, таких что Х*метро=1, если Хвид транспорта=метро и так далее.
Благодаря этому вместо весьма неоднозначной модели вида Y = bХвид транспорта мы переходим к удобной и гибкой модели вида Y = b1 Х*метро + b2 Х*автомобиль + b3Х*велосипед.

Зависимая переменная также может содержать как исходные данные — и это будет простой моделью — так и их преобразование. Причем когда преобразованная зависимая переменная принадлежит к экспоненциальному семейству распределений, то речь уже идет о так называемой обобщенной линейной модели (GLM, generalized linear model), к которым в частности относятся нормальная, логистическая, Пуассоновская, экспоненциальная, биномиальная и многие другие модели. Обобщенные модели очень важны и удобны в использовании, поскольку для них доказаны и параметры сходимости, и качества получаемых оценок, и влияние функционалов разных видов. В идеале старайтесь свести вашу задачу к какой-нибудь GLM-модели.

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

Не будем забывать и про случайные ошибки, которые могут иметь разные распределения, оказывая тем самым сильнейшее влияние на модель и метод ее построения. Например, logit и probit модель внешне устроены совершенно одинаково, принимают одни и те же данные и предсказывают вероятность некоторого события Y при заданных X. Вот только в probit модели ошибки распределены нормально, а в logit модели имеют логистическое распределение. Естественно и результат эти модели дают разный, поэтому не стоит их путать.

Функция потерь
С формулой модели вроде разобрались, переходим к оптимизируемому функционалу. И он может быть простым, когда учитывает только функцию потерь, то есть отличие предсказанных моделью значений от фактических. Типичными примерами простых функционалов являются:
— наименьшие квадраты: min ∑ (Yi* — bTXi*)2
— взвешенные наименьшие квадраты: min ∑ Wi (Yi* — bTXi*)2, так мы можем, например, недавним данным придать больший вес, и тем самым понизить значимость данных за прошлые годы.
— обобщенные наименьшие квадраты по расстоянию Махаланобиса: min ∑ (Yi* — bTXi*)T Ω-1 (Yi* — bTXi*)
— функция Хубера, которая интересна тем, что рядом с минимумом она ведет себя квадратичным образом, а в остальных местах линейно.
— обратная функция Хубера, которая, наоборот, везде квадратична, а в окрестности минимума линейна.
Этим, конечно, варианты возможных функционалов совсем не ограничиваются. Пожалуй, самым широким классом являются функционалы максимального правдоподобия. Введя параметрическое предположение о распределении случайных ошибок ℇ (а значит и Y*), можно в явном виде написать функцию правдоподобия, а затем, построив максимизирующий функционал, рассчитать искомые параметры.
Кстати, любопытный факт: если Y* распределено нормально, то функционал максимального правдоподобия фактически эквивалентен функционалу наименьших квадратов.Регуляризация
Сложный функционал содержит регуляризацию, которая обычно представлена в виде дополнительного регуляризационного слагаемого: min ℒ (b) + λ ℱ (b), где ℒ (b) — функция потерь, ℱ (b) — регуляризационная функция, λ — параметр, задающий степень влияния регуляризации.
Регуляризация предназначена для регулирования сложности модели и ее целью является упрощение модели. Это, в частности, помогает бороться с переобучением и позволяет увеличить обобщающую способность модели.

Типичные примеры регуляризационных функций:
1. L1 = ∑ |b|
Известная как LASSO-регуляризация (Least Absolute Shrinkage and Selection Operator), и, как несложно догадаться из названия, она позволяет снижать размерность коэффициентов, обращая некоторые из них в нули. И это весьма удобно, когда исходные данные сильно коррелированы.

2. L2 = ∑ |b|2
Иногда ее называют ridge-регуляризацией, и она позволяет минимизировать значения коэффициентов модели, а заодно сделать ее робастной к незначительным изменениям исходных данных. А еще она хорошо дифференцируется, а значит модель можно рассчитать аналитически.

3. LEN = α L1 + (1 — α) L2
Совмещая LASSO и ridge, получаем ElasticNet, которая объединяет два мира со всеми их плюсами и минусами.

4. LN = ∑ E [A (b, Ž)] — A (b, X), где А — log partition функция
Введя новую переменную вида Ž = X + ℥, где ℥ — случайная величина, мы фактически добавляем в исходные данные случайный шум, что очевидным образом помогает бороться с переобучением.
Для самой простой линейной регрессии введение аддитивного шума идентично L2-регуляризации, но для других моделей аддитивный шум может давать очень интересный результат. Например, в логистической регресии он по сути штрафует за предсказания близкие к ½ (проще говоря, поощряет категоричность предсказаний и наказывает за неопределенность).

5. Dropout
Еще один хитрый подход, активно применяемый в нейронных сетях. Введем новую переменную вида Ž = X * ℥, где ℥ — вектор Бернуллевских случайных величин длиной d. Проще говоря, мы случайным образом выбираем некоторое подмножество факторов Х и строим модель по ним, а потом выбираем такую модель, которая меньше всех зависит от этой случайности.
Для самой простой линейной регрессии dropout снова аналогичен L2-регуляризации. А вот, например, в логистической регрессии он позволяет учитывает влияние редких, но весьма характерных факторов (проще говоря, для некоторых очень маленьких Xij он будет подбирать большие коэффициенты bj, тем самым повышая их влияние на результат.

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

А теперь решение
После того как мы построили функционал можно приступать к его решению. И тут есть два основных пути:
— решение в аналитическом виде
— численное решение.

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

Обычно же аналитическое решение вообще недоступно, и приходится прибегать к численным методам и тут, конечно, мы сталкиваемся к гигантским разнообразием алгоритмов:
— стохастический градиентный спуск
— стохастический средний градиент
— метод сопряженных градиентов
— Алгоритм Бройдена — Флетчера — Гольдфарба — Шанно, а также его модификация с ограниченной памятью L-BFGS.

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