Простые методы оценки параметров моделей

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

где \(y_t\) — объём продаж мороженого, \(x_t\) — температура воздуха, а \(e_t\) — случайная ошибка модели.

Вроде бы ничего особенного и ничего сложного, но откуда же нам взять коэффициенты \(a\) и \(b\)? Надо их как-то рассчитать (или, как это говорят в статистике, «оценить»). Для этого существует множество различных методов, но все они преследуют одну и ту же цель — оптимизировать какую-нибудь функцию (минимизировать или максимизировать) путём подбора этих самых значений \(a\) и \(b\). По сути вся задача сводится к тому, чтобы провести некоторым образом линию через «облако» значений \(x\) и \(y\) так, чтобы минимизировать целевую функцию. На рисунке ниже показан пример такого облака и нескольких линий, проходящих через него.

Мороженое и палочки

Каждая из линий на рисунке соответствует одной и той же линейной функции \eqref, но с разными коэффициентами. Фактически любая такая линия — это последовательность некоторых расчётных значений \(\hat_t\), позволяющая оценить, какими будут продажи мороженного в среднем если температура будет принимать значения \(x_t\) (отложенные по оси абсцисс). Очевидно, что таких прямых на плоскости можно провести сколько угодно, и каждая из них будет минимизировать либо максимизировать какой-нибудь свой критерий. Обычно это происходит путём подбора коэффициентов с использованием различных оптимизационных методов, путём решения задач программирования (линейного, нелинейного, целочисленного) либо с использованием эвристических методов (генетические и эволюционные алгоритмы). В достаточно редких случаях для расчёта коэффициентов в нашем распоряжении могут быть формулы, гарантирующие получения оптимума (такие формулы даёт, например, Метод Наименьших Квадратов, который мы обсудим в главе про регрессии).

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

Рассмотрим подробней самые простые и популярные в прогнозировании методы оценки коэффициентов моделей.

Первой строкой мы создаём нормально распределённую величину, состоящую из 100 наблюдений. Второй строкой мы создаём вторую переменную, состоящую из значений "NA" - "Np (то есть никак не заданных). Третьей и четвёртой строками мы заполняем вторую переменную так, что первые 95 наблюдений сгенерированы исходя из одной зависимости от x, а последние 5 - исходя из другой. Следующими четырьмя строками мы создаём из этих двух переменных объект с данными. Строкой за ними мы создаём вектор остатков. Пока что он пустой. Последней строкой мы создаём вектор для двух коэффициентов нашей будущей модели. Они тоже пока пустые.

MSE расшифровывается как «Mean Squared Error» и переводится как «Средняя квадратическая ошибка». Суть метода заключается по сути в том, чтобы минимизировать сумму квадратов отклонений фактических значений от расчётных (SSE - «Sum of Squared Errors»). Если полученную сумму разделить на число наблюдений, то получится та самая MSE. Формула целевой функции в этом случае выглядит следующим образом:

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

Стоит отметить, что, как мы помним из «Статистического анализа», возведение в квадрат существенно увеличивает те значения, которые лежат далеко от всех остальных. Например, если продажи мороженого колеблются в основном в районе ста штук, но есть наблюдение в 200, то это наблюдение будет тянуть на себя одеяло — влиять существенно на оценки коэффициентов. То есть MSE является не робастной величиной и вообще позволяет получить некоторую среднюю оценку (в среднем фактическое значение \(y_t\) будет соответствовать \(\hat_t\) при данном значении \(x_t\)).

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

Самой первой строкой мы определяем, что "CF" это объект типа "функция". Далее мы указываем, как рассчитать остатки модели: это фактические значения (второй столбец наших данных) минус расчётные, полученные по формуле \eqref (здесь \(C[1] = a \text C[2] = b \text data[,1] = x\)). Третьей строкой мы определяем нашу целевую функцию. В данном случае это средняя, взятая по остаткам, возведённым в квадрат (MSE).

Теперь проведём оптимизацию и запишем полученные коэффициенты:

Функция "nlminb" позволяет подбирать оптимальные параметры путём минимизации некой функции. c(1,1) в данном случае задаёт функции вектор стартовых параметров для поиска. Вместо этих значений можно было бы задать и что-нибудь другое. Второй строкой мы записываем найденные параметры в наш вектор C.

Построим теперь точечную диаграмму по ряду данных и проведём на том же графике полученную линию:

С первой функцией мы уже сталкивались ранее, а вот со второй - пока нет. Она позволяет провести прямую линию с заданными параметрами на уже имеющемся графике. a и b соответствуют тем самым коэффициентам той самой функции \eqref.

MAE расшифровывается как «Mean Absolute Error», что с зарубежного на отечественный переводится как «Средняя абсолютная ошибка». Рассчитывается она так:

Модули в формуле \eqref всё так же позволяют избавиться от знаков и получить некоторую оценку расстояния от фактических до расчётных значений, которое нужно будет потом минимизировать. Несомненным преимуществом MAE является то, что модули не увеличивают в разы отклонения, считающиеся выбросами. Поэтому эта оценка является более робастной, чем MSE и фактически соответствует медиане (в 50% случаев при данном значении \(x_t\) зависимая величина \(y_t\) будет не меньше полученной \(\hat_t\)).

Разница между оценками коэффициентов, полученными при использовании MSE и MAE, наглядно показана на следующем рисунке:

Модели, оценённые с использованием MSE и MAE.

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

Теперь функция возвращает средние абсолютные ошибки (функция "abs" возвращает значения по модулю), а не квадратические.

Ещё раз проведём оптимизацию, запишем коэффициенты и построим график:

Если вы проделали всё это, то теперь можете сравнить полученные линии в случае с MSE и в случае с MAE. Здесь я специально опустил функцию "plot" для того, чтобы новая линия для наглядности была нанесена на тот же график.

Прозорливые читатели уже, наверно, догадались, что помимо получения средней и медианной величин, можно так же получить и квантильные. В статистическом анализе регрессии, оценённые таким методом, называются «квантильными регрессиями».

Целевая функция, соответствующая таким квантильным оценкам, имеет вид:

\begin \label CF_\alpha = (1 - \alpha) \sum_T\), то мы получим следующую функцию правдоподобия (здесь нужно иметь в виду, что при этом предполагается независимость наблюдений):

где \(Y\) - это вектор всех фактических значений: \(Y = \beginy_1 \\ y_2 \\ \vdots \\ y_T \end\)

Теперь мы можем подставлять в \eqref разные значения параметров модели \(\theta \text < и >\sigma^2 \) и получать разные вероятности. Естественно, обычно мы заинтересованы в том, чтобы вероятность была максимальной, так что функцию правдоподобия \eqref можно максимизировать, задавая разные значения параметров, и получить какие-нибудь оптимальные оценки коэффициентов.

Всё. Счастье есть! Можно откинуться на спинку стула и потягивать чай / мохито / виски / ром - в зависимости от времени года и предпочтений прогнозиста.

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

где \(e_t = y_t - \hat_t\).

Подставляя \eqref в \eqref, мы получим следующую функцию правдоподобия на основе нормального распределения:

Очевидно, что функцию \eqref неудобно максимизировать в таком её виде, так как приходится сталкиваться с произведениями и экспонированием. Однако оптимальное значение коэффициентов не изменится, если мы линеаризуем эту функцию (например, прологарифмируем левую и правую части, используя натуральный логарифм). В этом случае мы получим следующую функцию (которая у них там на западе называется "log-likelihood"):

\begin \label \ell(\theta, \sigma^2 | Y) = \ln ( L(\theta, \sigma^2 | Y) ) = -\frac \ln \left( 2 \pi \sigma^2 \right) -\frac \sum_^T \frac . \end

В правой части \eqref дисперсию можно вынести за знак суммы, так как она представлена константой, тогда получим:

Теперь можно попытаться математически вывести при каком значении параметров функция будет максимизироваться. Проще всего это сделать вначале с дисперсией. Для этого возьмём первую производную \eqref по \(\sigma^2\) и приравняем её к нулю:

Если теперь в \eqref сократить то, что сокращается и перенести правое слагаемое в правую часть, то получим:

Умножим теперь левую и правую части \eqref на \(\sigma^3\) и разделим на \(T\):

Полученная формула позволяет нам оценить дисперсию ошибок по выборке с помощью статистически выверенного и корректного метода - путём максимизации функции правдоподобия на основе нормального закона распределения вероятностей. Ничего нового нам эта формула, вроде бы, не даёт, потому что мы итак знаем, как рассчитать дисперсию ошибок, но нам она полезна как раз с точки зрения статистического обоснования. Однако, чтобы не перепутать, что мы вообще-то имеем дело не с истинной дисперсией (существующей в генеральной совокупности и воображении статистиков), а с оценённой по выборке, мы будем полученную с помощью \eqref дисперсию обозначать \(\hat^2\). Кроме того, отметим, что мы только что показали, что дисперсия, оценённая по формуле \eqref является некой оптимальной оценкой только в случае, если мы имеем дело с нормальным распределение ошибок!

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

Заметим, пока не поздно, что функция правдоподобия, в которой используется оценённая дисперсия, называется "концентрированной" ("concentrated log-likelihood").

В \eqref суммы в правой части легко сокращаются, а произведение под знаком логарифма легко преобразуется в сумму логарифмов таким вот образом:

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

\begin \label \ell(\theta | Y) = -\frac \left( \ln \left( 2 \pi e \right) + \ln \left( \frac \sum_^T e_t^2 \right) \right) . \end

Максимизировать функцию \eqref теперь очень просто. На самом деле ещё легче - минимизировать \(-\frac \ell(\theta | Y)\), так как размер выборки у нас фиксирован. \eqref тогда преобразуется в:

Ну, и в качестве финального аккорда, можно избавиться от константы \(\ln \left( 2 \pi e \right)\), так как она никак не влияет на получение оптимальных коэффициентов. А полученное в итоге число можно вообще проэкспонировать, для того, чтобы прийти к следующей целевой функции, минимизируя которую мы будем получать оценки максимального правдоподобия (не забываем про допущение о нормальном распределении остатков модели):

Удивительным образом эта целевая функция целиком и полностью соответствует MSE (формула \eqref).

В этом месте может возникнуть вопрос: зачем же все эти выводы и сложные процедуры, если мы в итоге получили всего лишь известный нам критерий MSE? Ответов на этот вопрос несколько:

1. Путём простых математических выкладок мы доказали, что MSE является эффективным критерием оценки коэффициентов моделей в том случае, если остатки модели распределены нормально. Например, для наших примеров в предыдущих параграфах (см. рисунок с MSE и MAE) эта предпосылка нарушается, поэтому MSE не является эффективным критерием.
2. Используя логарифмированную концентрированную функцию правдоподобия \eqref можно осуществлять выбор наилучшей модели из ряда имеющихся в распоряжении исследователя. Эту тему мы обсудим подробней в одном из следующих параграфов.
3. Кроме того, всё та же функция позволяет нам получать оценки дисперсий коэффициентов любой модели. Они нам могут понадобиться для проверки всяких статистических гипотез или построения доверительных интервалов. Для получения этих дисперсий нужно обратиться к некой вещи под названием "Информация Фишера", но здесь мы вдаваться в детали пока не будем.
4. Если нам не нравится предположение о нормальности распределения остатков, мы можем использовать любую другую функцию распределения в \eqref. Конечно же, придётся произвести все последующие выводы аналогично тем, что мы произвели только что (а в некоторых случаях ничего такого сделать и не получится, функцию придётся максимизировать как есть), но зато это даст нам красивые, статистически выверенные оценки.
5. Используя этот подход, можно показать, что минимизация MAE является оптимальным методом оценки в случае, если остатки имеют распределение Лапласа, а квантиальные оценки дают максимум функции правдоподобия для асимметричного распределения Лапласа.

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