Средняя квадратичная ошибка прогноза

From Wikipedia, the free encyclopedia

In statistics the mean squared prediction error (MSPE), also known as mean squared error of the predictions, of a smoothing, curve fitting, or regression procedure is the expected value of the squared prediction errors (PE), the square difference between the fitted values implied by the predictive function widehat {g} and the values of the (unobservable) true value g. It is an inverse measure of the explanatory power of {displaystyle {widehat {g}},} and can be used in the process of cross-validation of an estimated model.
Knowledge of g would be required in order to calculate the MSPE exactly; in practice, MSPE is estimated.[1]

Formulation[edit]

If the smoothing or fitting procedure has projection matrix (i.e., hat matrix) L, which maps the observed values vector y to predicted values vector {displaystyle {hat {y}}=Ly,} then PE and MSPE are formulated as:

{displaystyle operatorname {PE_{i}} =g(x_{i})-{widehat {g}}(x_{i}),}
{displaystyle operatorname {MSPE} =operatorname {E} left[operatorname {PE} _{i}^{2}right]=sum _{i=1}^{n}operatorname {PE} _{i}^{2}/n.}

The MSPE can be decomposed into two terms: the squared bias (mean error) of the fitted values and the variance of the fitted values:

{displaystyle operatorname {MSPE} =operatorname {ME} ^{2}+operatorname {VAR} ,}
{displaystyle operatorname {ME} =operatorname {E} left[{widehat {g}}(x_{i})-g(x_{i})right]}
{displaystyle operatorname {VAR} =operatorname {E} left[left({widehat {g}}(x_{i})-operatorname {E} left[{g}(x_{i})right]right)^{2}right].}

The quantity SSPE=nMSPE is called sum squared prediction error.
The root mean squared prediction error is the square root of MSPE: RMSPE=MSPE.

Computation of MSPE over out-of-sample data[edit]

The mean squared prediction error can be computed exactly in two contexts. First, with a data sample of length n, the data analyst may run the regression over only q of the data points (with q < n), holding back the other n – q data points with the specific purpose of using them to compute the estimated model’s MSPE out of sample (i.e., not using data that were used in the model estimation process). Since the regression process is tailored to the q in-sample points, normally the in-sample MSPE will be smaller than the out-of-sample one computed over the n – q held-back points. If the increase in the MSPE out of sample compared to in sample is relatively slight, that results in the model being viewed favorably. And if two models are to be compared, the one with the lower MSPE over the n – q out-of-sample data points is viewed more favorably, regardless of the models’ relative in-sample performances. The out-of-sample MSPE in this context is exact for the out-of-sample data points that it was computed over, but is merely an estimate of the model’s MSPE for the mostly unobserved population from which the data were drawn.

Second, as time goes on more data may become available to the data analyst, and then the MSPE can be computed over these new data.

Estimation of MSPE over the population[edit]

When the model has been estimated over all available data with none held back, the MSPE of the model over the entire population of mostly unobserved data can be estimated as follows.

For the model y_{i}=g(x_{i})+sigma varepsilon _{i} where varepsilon _{i}sim {mathcal  {N}}(0,1), one may write

{displaystyle ncdot operatorname {MSPE} (L)=g^{text{T}}(I-L)^{text{T}}(I-L)g+sigma ^{2}operatorname {tr} left[L^{text{T}}Lright].}

Using in-sample data values, the first term on the right side is equivalent to

{displaystyle sum _{i=1}^{n}left(operatorname {E} left[g(x_{i})-{widehat {g}}(x_{i})right]right)^{2}=operatorname {E} left[sum _{i=1}^{n}left(y_{i}-{widehat {g}}(x_{i})right)^{2}right]-sigma ^{2}operatorname {tr} left[left(I-Lright)^{T}left(I-Lright)right].}

Thus,

{displaystyle ncdot operatorname {MSPE} (L)=operatorname {E} left[sum _{i=1}^{n}left(y_{i}-{widehat {g}}(x_{i})right)^{2}right]-sigma ^{2}left(n-operatorname {tr} left[Lright]right).}

If sigma ^{2} is known or well-estimated by widehat {sigma }^{2}, it becomes possible to estimate MSPE by

{displaystyle ncdot operatorname {widehat {MSPE}} (L)=sum _{i=1}^{n}left(y_{i}-{widehat {g}}(x_{i})right)^{2}-{widehat {sigma }}^{2}left(n-operatorname {tr} left[Lright]right).}

Colin Mallows advocated this method in the construction of his model selection statistic Cp, which is a normalized version of the estimated MSPE:

{displaystyle C_{p}={frac {sum _{i=1}^{n}left(y_{i}-{widehat {g}}(x_{i})right)^{2}}{{widehat {sigma }}^{2}}}-n+2p.}

where p the number of estimated parameters p and widehat {sigma }^{2} is computed from the version of the model that includes all possible regressors.
That concludes this proof.

See also[edit]

  • Akaike information criterion
  • Bias-variance tradeoff
  • Mean squared error
  • Errors and residuals in statistics
  • Law of total variance
  • Mallows’s Cp
  • Model selection

References[edit]

  1. ^ Pindyck, Robert S.; Rubinfeld, Daniel L. (1991). «Forecasting with Time-Series Models». Econometric Models & Economic Forecasts (3rd ed.). New York: McGraw-Hill. pp. 516–535. ISBN 0-07-050098-3.

Среднеквадратичная ошибка (Mean Squared Error) – Среднее арифметическое (Mean) квадратов разностей между предсказанными и реальными значениями Модели (Model) Машинного обучения (ML):

MSE как среднее дистанций между предсказаниями и реальными наблюдениями

Рассчитывается с помощью формулы, которая будет пояснена в примере ниже:

$$MSE = frac{1}{n} × sum_{i=1}^n (y_i — widetilde{y}_i)^2$$
$$MSEspace{}{–}space{Среднеквадратическая}space{ошибка,}$$
$$nspace{}{–}space{количество}space{наблюдений,}$$
$$y_ispace{}{–}space{фактическая}space{координата}space{наблюдения,}$$
$$widetilde{y}_ispace{}{–}space{предсказанная}space{координата}space{наблюдения,}$$

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

Пример. Исследуем линейную регрессию, изображенную на графике выше, и установим величину среднеквадратической Ошибки (Error). Фактические координаты точек-Наблюдений (Observation) выглядят следующим образом:

Мы имеем дело с Линейной регрессией (Linear Regression), потому уравнение, предсказывающее положение записей, можно представить с помощью формулы:

$$y = M * x + b$$
$$yspace{–}space{значение}space{координаты}space{оси}space{y,}$$
$$Mspace{–}space{уклон}space{прямой}$$
$$xspace{–}space{значение}space{координаты}space{оси}space{x,}$$
$$bspace{–}space{смещение}space{прямой}space{относительно}space{начала}space{координат}$$

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

$$y = 0,5252 * x + 17,306$$

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

Рассчитаем квадрат разницы между Y и Ỹ:

Сумма таких квадратов равна 4 445. Осталось только разделить это число на количество наблюдений (9):

$$MSE = frac{1}{9} × 4445 = 493$$

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

MSE и Scikit-learn

Среднеквадратическую ошибку можно вычислить с помощью SkLearn. Для начала импортируем функцию:

import sklearn
from sklearn.metrics import mean_squared_error

Инициализируем крошечные списки, содержащие реальные и предсказанные координаты y:

y_true = [5, 41, 70, 77, 134, 68, 138, 101, 131]
y_pred = [23, 35, 55, 90, 93, 103, 118, 121, 129]

Инициируем функцию mean_squared_error(), которая рассчитает MSE тем же способом, что и формула выше:

mean_squared_error(y_true, y_pred)

Интересно, что конечный результат на 3 отличается от расчетов с помощью Apple Numbers:

496.0

Ноутбук, не требующий дополнительной настройки на момент написания статьи, можно скачать здесь.

Автор оригинальной статьи: @mmoshikoo

Фото: @tobyelliott

Электронный учебник Statsoft

Анализ временных рядов


  • Общее введение
  • Две основные цели
  • Идентификация модели
    временных рядов

    • Систематическая
      составляющая и случайный шум
    • Два общих типа компонент
      временных рядов
    • Анализ тренда
    • Анализ сезонности
  • АРПСС (Бокс и Дженкинс) и
    автокорреляции

    • Общее введение
    • Два основных процесса
    • Модель АРПСС
    • Идентификация
    • Оценивание параметров
    • Оценивание модели
  • Прерванные временные ряды
  • Экспоненциальное
    сглаживание

    • Общее введение
    • Простое экспоненциальное
      сглаживание
    • Выбор лучшего значения
      параметра a (альфа)
    • Индексы качества подгонки
    • Сезонная и несезонная модели
      с трендом или без тренда
  • Сезонная декомпозиция (метод
    Census I)

    • Общее введение
    • Вычисления
  • Сезонная корректировка X-11
    (метод Census II)

    • Сезонная корректировка:
      основные идеи и термины
    • Метод Census II
    • Таблицы результатов
      корректировки X-11
    • Подробное описание всех
      таблиц результатов, вычисляемых в методе X-11
  • Анализ распределенных лагов
    • Общая цель
    • Общая модель
    • Распределенный лаг Алмона
  • Одномерный анализ Фурье
  • Кросс-спектральный анализ
    • Общее введение
    • Основные понятия и принципы
    • Результаты для каждой
      переменной
    • Кросс-периодограмма,
      кросс-плотность, квадратурная плотность и
      кросс-амплитуда
    • Квадрат когерентности,
      усиление и фазовый сдвиг
    • Как создавались данные для
      примера
  • Спектральный анализ —
    Основные понятия и принципы

    • Частота и период
    • Общая структура модели
    • Простой пример
    • Периодограмма
    • Проблема рассеяния
    • Добавление констант во
      временной ряд (пэддинг)
    • Косинус-сглаживание
    • Окна данных и оценки
      спектральной плотности
    • Подготовка данных к анализу
    • Результаты для случая, когда в
      ряде отсутствует периодичность
  • Быстрое преобразование Фурье
    • Общее введение
    • Вычисление БПФ во временных
      рядах

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

Общее введение

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

Подробное обсуждение этих методов можно найти
в следующих работах: Anderson (1976), Бокс и Дженкинс
(1976), Kendall (1984), Kendall and Ord (1990), Montgomery, Johnson, and Gardiner (1990),
Pankratz (1983), Shumway (1988), Vandaele (1983), Walker (1991), Wei (1989).

Две основные цели

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


Идентификация модели временных
рядов

  • Систематическая
    составляющая и случайный шум
  • Два общих типа компонент
    временных рядов
  • Анализ тренда
  • Анализ сезонности

За более полной информацией о простых
автокорреляциях (обсуждаемых в этом разделе) и
других автокорреляциях, см. Anderson (1976), Box and Jenkins
(1976), Kendall (1984), Pankratz (1983), and Vandaele (1983). См. также:

  • АРПСС (Бокс и Дженкинс) и
    автокорреляции
  • Прерванные временные ряды
  • Экспоненциальное
    сглаживание
  • Сезонная декомпозиция (метод
    Census I)
  • Сезонная корректировка X-11
    (метод Census II)
  • Таблицы результатов
    корректировки X-11
  • Анализ распределенных лагов
  • Одномерный анализ Фурье
  • Кросс-спектральный анализ
  • Основные понятия и принципы
  • Быстрое преобразование Фурье

Систематическая составляющая и
случайный шум

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

Два общих типа компонент
временных рядов

Большинство регулярных составляющих временных
рядов принадлежит к двум классам: они являются
либо трендом, либо сезонной составляющей. Тренд
представляет собой общую систематическую
линейную или нелинейную компоненту, которая
может изменяться во времени. Сезонная
составляющая — это периодически повторяющаяся
компонента. Оба эти вида регулярных компонент
часто присутствуют в ряде одновременно.
Например, продажи компании могут возрастать из
года в год, но они также содержат сезонную
составляющую (как правило, 25% годовых продаж
приходится на декабрь и только 4% на август).

График

Эту общую модель можно понять на
«классическом» ряде — Ряд G (Бокс и
Дженкинс, 1976, стр. 531), представляющем месячные
международные авиаперевозки (в тысячах) в
течение 12 лет с 1949 по 1960 (см. файл Series_g.sta).
График месячных перевозок ясно показывает почти
линейный тренд, т.е. имеется устойчивый рост
перевозок из года в год (примерно в 4 раза больше
пассажиров перевезено в 1960 году, чем в 1949). В то же
время характер месячных перевозок повторяется,
они имеют почти один и тот же характер в каждом
годовом периоде (например, перевозок больше в
отпускные периоды, чем в другие месяцы). Этот
пример показывает довольно определенный тип
модели временного ряда, в которой амплитуда
сезонных изменений увеличивается вместе с
трендом. Такого рода модели называются моделями
с мультипликативной сезонностью.

Анализ тренда

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

Сглаживание. Сглаживание всегда включает
некоторый способ локального усреднения данных,
при котором несистематические компоненты
взаимно погашают друг друга. Самый общий метод
сглаживания — скользящее среднее, в котором
каждый член ряда заменяется простым или
взвешенным средним n соседних членов, где n
— ширина «окна» (см. Бокс и Дженкинс, 1976; Velleman
and Hoaglin, 1981). Вместо среднего можно использовать
медиану значений, попавших в окно. Основное
преимущество медианного сглаживания, в
сравнении со сглаживанием скользящим средним,
состоит в том, что результаты становятся более
устойчивыми к выбросам (имеющимся внутри окна).
Таким образом, если в данных имеются выбросы
(связанные, например, с ошибками измерений), то
сглаживание медианой обычно приводит к более
гладким или, по крайней мере, более
«надежным» кривым, по сравнению со
скользящим средним с тем же самым окном. Основной
недостаток медианного сглаживания в том, что при
отсутствии явных выбросов, он приводит к более
«зубчатым» кривым (чем сглаживание
скользящим средним) и не позволяет использовать
веса.

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

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

Анализ сезонности

Периодическая и сезонная зависимость
(сезонность) представляет собой другой общий тип
компонент временного ряда. Это понятие было
проиллюстрировано ранее на примере
авиаперевозок пассажиров. Можно легко видеть,
что каждое наблюдение очень похоже на соседнее;
дополнительно, имеется повторяющаяся сезонная
составляющая, это означает, что каждое
наблюдение также похоже на наблюдение, имевшееся
в том же самом месяце год назад. В общем,
периодическая зависимость может быть формально
определена как корреляционная зависимость
порядка k между каждым i-м элементом ряда и
(i-k)-м элементом (Kendall, 1976). Ее можно измерить с
помощью автокорреляции (т.е. корреляции между
самими членами ряда); k обычно называют лагом
(иногда используют эквивалентные термины:
сдвиг, запаздывание). Если ошибка измерения не
слишком большая, то сезонность можно определить
визуально, рассматривая поведение членов ряда
через каждые k временных единиц.

Автокорреляционная коррелограмма. Сезонные
составляющие временного ряда могут быть найдены
с помощью коррелограммы. Коррелограмма
(автокоррелограмма) показывает численно и
графически автокорреляционную функцию (AКФ),
иными словами коэффициенты автокорреляции (и их
стандартные ошибки) для последовательности
лагов из определенного диапазона (например, от 1
до 30). На коррелограмме обычно отмечается
диапазон в размере двух стандартных ошибок на
каждом лаге, однако обычно величина
автокорреляции более интересна, чем ее
надежность, потому что интерес в основном
представляют очень сильные (а, следовательно,
высоко значимые) автокорреляции (см. Элементарные
понятия статистики
).

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

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

Частные автокорреляции. Другой полезный
метод исследования периодичности состоит в
исследовании частной автокорреляционной
функции (ЧАКФ), представляющей собой
углубление понятия обычной автокорреляционной
функции. В ЧАКФ устраняется зависимость между
промежуточными наблюдениями (наблюдениями внутри
лага). Другими словами, частная автокорреляция на
данном лаге аналогична обычной автокорреляции,
за исключением того, что при вычислении из нее
удаляется влияние автокорреляций с меньшими
лагами (см. Бокс и Дженкинс, 1976; см. также McDowall,
McCleary, Meidinger, and Hay, 1980). На лаге 1 (когда нет
промежуточных элементов внутри лага), частная
автокорреляция равна, очевидно, обычной
автокорреляции. На самом деле, частная
автокорреляция дает более «чистую» картину
периодических зависимостей.

Удаление периодической зависимости. Как
отмечалось выше, периодическая составляющая для
данного лага k может быть удалена взятием
разности соответствующего порядка. Это означает,
что из каждого i-го элемента ряда вычитается (i-k)
элемент. Имеются два довода в пользу таких
преобразований.

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

Во-вторых, удаление сезонных составляющих
делает ряд стационарным,
что необходимо для применения АРПСС
и других методов, например, спектрального
анализа.


АРПСС

  • Общее введение
  • Два основных процесса
  • Модель АРПСС
  • Идентификация
  • Оценивание параметров
  • Оценивание модели

Дополнительная информация о методах Анализа
временных рядов
дана также в следующих
разделах:

  • Идентификация модели
    временных рядов
  • Прерванные временные ряды
  • Экспоненциальное
    сглаживание
  • Сезонная декомпозиция (метод
    Census I)
  • Сезонная корректировка X-11
    (метод Census II)
  • Таблицы результатов
    корректировки X-11
  • Анализ распределенных лагов
  • Одномерный анализ Фурье
  • Кросс-спектральный анализ
  • Основные понятия и принципы
  • Быстрое преобразование Фурье

Общее введение

Процедуры оценки параметров и прогнозирования,
описанные в разделе Идентификация
модели временных рядов
, предполагают, что
математическая модель процесса известна. В
реальных данных часто нет отчетливо выраженных
регулярных составляющих. Отдельные наблюдения
содержат значительную ошибку, тогда как вы
хотите не только выделить регулярные компоненты,
но также построить прогноз. Методология АРПСС,
разработанная Боксом и Дженкинсом (1976), позволяет
это сделать. Данный метод чрезвычайно популярен
во многих приложениях, и практика подтвердила
его мощность и гибкость (Hoff, 1983; Pankratz, 1983; Vandaele, 1983).
Однако из-за мощности и гибкости, АРПСС — сложный
метод. Его не так просто использовать, и
требуется большая практика, чтобы овладеть им.
Хотя часто он дает удовлетворительные
результаты, они зависят от квалификации
пользователя (Bails and Peppers, 1982). Следующие разделы
познакомят вас с его основными идеями. Для
интересующихся кратким, рассчитанным на
применение, (нематематическим) введением в АРПСС,
рекомендуем книгу McCleary, Meidinger, and Hay (1980).

Два основных процесса

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

xt =
+ 1*x(t-1) + 2*x(t-2) + 3*x(t-3) + … +

Здесь:
                 —
константа (свободный член),
 1,
2,
3  
— параметры авторегрессии.

Вы видите, что каждое наблюдение есть сумма
случайной компоненты (случайное воздействие, errorblu.gif (835 bytes)) и линейной
комбинации предыдущих наблюдений.

Требование стационарности. Заметим, что
процесс авторегрессии будет стационарным
только, если его параметры лежат в определенном
диапазоне. Например, если имеется только один
параметр, то он должен находиться в интервале -1<<+1. В противном случае,
предыдущие значения будут накапливаться и
значения последующих xt могут быть
неограниченными, следовательно, ряд не будет стационарным.
Если имеется несколько параметров
авторегрессии, то можно определить аналогичные
условия, обеспечивающие стационарность (см.
например, Бокс и Дженкинс, 1976; Montgomery, 1990).

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

xt = µ + t1*(t-1)2*(t-2)3*(t-3) — …

Здесь:
 µ                —
константа,
 1,
2,
3  —
параметры скользящего среднего.

Другими словами, текущее наблюдение ряда
представляет собой сумму случайной компоненты
  (случайное воздействие, errorblu.gif (835 bytes)) в данный момент и линейной
комбинации случайных воздействий в предыдущие
моменты времени.

Обратимость. Не вдаваясь в детали, отметим,
что существует «двойственность» между
процессами скользящего среднего и авторегрессии
(см. например, Бокс и Дженкинс, 1976; Montgomery, Johnson, and
Gardiner, 1990). Это означает, что приведенное выше
уравнение скользящего среднего можно переписать
(обратить) в виде уравнения авторегрессии
(неограниченного порядка), и наоборот. Это так
называемое свойство обратимости. Имеются
условия, аналогичные приведенным выше условиям стационарности,
обеспечивающие обратимость модели.

Модель АРПСС

Модель авторегрессии и скользящего среднего. Общая
модель, предложенная Боксом и Дженкинсом (1976)
включает как параметры авторегрессии, так и
параметры скользящего среднего. Именно, имеется
три типа параметров модели: параметры
авторегрессии (p), порядок разности (d), параметры
скользящего среднего (q). В обозначениях Бокса
и Дженкинса модель записывается как АРПСС (p, d, q).
Например, модель (0, 1, 2) содержит 0
(нуль) параметров авторегрессии (p) и 2
параметра скользящего среднего (q), которые
вычисляются для ряда после взятия разности с
лагом 1.

Идентификация. Как отмечено ранее, для
модели АРПСС необходимо, чтобы ряд был стационарным,
это означает, что его среднее постоянно, а
выборочные дисперсия и автокорреляция не
меняются во времени. Поэтому обычно необходимо
брать разности ряда до тех пор, пока он не станет
стационарным
(часто также применяют логарифмическое
преобразование для стабилизации дисперсии).
Число разностей, которые были взяты, чтобы
достичь стационарности, определяются параметром
d (см. предыдущий раздел). Для того чтобы
определить необходимый порядок разности, нужно
исследовать график ряда и автокоррелограмму.
Сильные изменения уровня (сильные скачки вверх
или вниз) обычно требуют взятия несезонной
разности первого порядка (лаг=1). Сильные
изменения наклона требуют взятия разности
второго порядка. Сезонная составляющая требует
взятия соответствующей сезонной разности (см.
ниже). Если имеется медленное убывание
выборочных коэффициентов автокорреляции в
зависимости от лага, обычно берут разность
первого порядка. Однако следует помнить, что для
некоторых временных рядов нужно брать разности
небольшого порядка или вовсе не брать их.
Заметим, что чрезмерное количество взятых
разностей
приводит к менее стабильным оценкам
коэффициентов.

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

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

Константа в моделях АРПСС. Дополнительно
модели АРПСС могут содержать константу,
интерпретация которой зависит от подгоняемой
модели. Именно, если (1) в модели нет параметров
авторегрессии, то константа есть среднее значение ряда, если (2)
параметры авторегрессии имеются, то константа
представляет собой свободный член. Если бралась
разность ряда, то константа представляет собой
среднее или свободный член преобразованного
ряда. Например, если бралась первая разность
(разность первого порядка), а параметров
авторегрессии в модели нет, то константа
представляет собой среднее значение
преобразованного ряда и, следовательно, коэффициент
наклона линейного тренда
исходного.

Идентификация

Число оцениваемых параметров. Конечно, до
того, как начать оценивание, вам необходимо
решить, какой тип модели будет подбираться к
данным, и какое количество параметров
присутствует в модели, иными словами, нужно
идентифицировать модель АРПСС. Основными
инструментами идентификации порядка модели
являются графики, автокорреляционная функция
(АКФ), частная автокорреляционная функция (ЧАКФ).
Это решение не является простым и требуется
основательно поэкспериментировать с
альтернативными моделями. Тем не менее,
большинство встречающихся на практике временных
рядов можно с достаточной степенью точности
аппроксимировать одной из 5 основных моделей (см.
ниже), которые можно идентифицировать по виду
автокорреляционной (АКФ) и частной
автокорреляционной функции (ЧАКФ). Ниже дается
список этих моделей, основанный на рекомендациях
Pankratz (1983); дополнительные практические советы
даны в Hoff (1983), McCleary and Hay (1980), McDowall, McCleary, Meidinger, and Hay
(1980), and Vandaele (1983). Отметим, что число параметров
каждого вида невелико (меньше 2), поэтому нетрудно
проверить альтернативные модели.

  1. Один параметр (p): АКФ — экспоненциально
    убывает; ЧАКФ — имеет резко выделяющееся значение
    для лага 1, нет корреляций на других лагах.
  2. Два параметра авторегрессии (p): АКФ имеет
    форму синусоиды или экспоненциально убывает;
    ЧАКФ имеет резко выделяющиеся значения на лагах 1,
    2, нет корреляций на других лагах.
  3. Один параметр скользящего среднего (q): АКФ
    имеет резко выделяющееся значение на лаге 1,
    нет корреляций на других лагах. ЧАКФ
    экспоненциально убывает.
  4. Два параметра скользящего среднего (q): АКФ
    имеет резко выделяющиеся значения на лагах 1, 2,
    нет корреляций на других лагах. ЧАКФ имеет форму
    синусоиды или экспоненциально убывает.
  5. Один параметр авторегрессии (p) и один параметр
    скользящего среднего (q)
    : АКФ экспоненциально
    убывает с лага 1; ЧАКФ — экспоненциально
    убывает с лага 1.

Сезонные модели. Мультипликативная сезонная
АРПСС представляет естественное развитие и
обобщение обычной модели АРПСС на ряды, в которых
имеется периодическая сезонная компонента. В
дополнении к несезонным параметрам, в модель
вводятся сезонные параметры для определенного
лага (устанавливаемого на этапе идентификации
порядка модели). Аналогично параметрам простой
модели АРПСС, эти параметры называются: сезонная
авторегрессия (ps), сезонная разность (ds) и
сезонное скользящее среднее (qs). Таким
образом, полная сезонная АРПСС может быть
записана как АРПСС (p,d,q)(ps,ds,qs).
Например, модель (0,1,2)(0,1,1) включает 0
регулярных параметров авторегрессии, 2
регулярных параметра скользящего среднего и 1
параметр сезонного скользящего среднего. Эти
параметры вычисляются для рядов, получаемых
после взятия одной разности с лагом 1 и далее
сезонной разности. Сезонный лаг, используемый
для сезонных параметров, определяется на этапе
идентификации порядка модели.

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

Оценивание параметров

Существуют различные методы оценивания
параметров, которые дают очень похожие оценки, но
для данной модели одни оценки могут быть более
эффективны, а другие менее эффективны. В общем, во
время оценивания порядка модели используется
так называемый квазиньютоновский алгоритм
максимизации правдоподобия (вероятности)
наблюдения значений ряда по значениям
параметров (см. Нелинейное
оценивание
). Практически это требует
вычисления (условных) сумм квадратов (SS)
остатков модели. Имеются различные способы
вычисления суммы квадратов остатков SS; вы
можете выбрать: (1) приближенный метод
максимального правдоподобия МакЛеода и Сейлза
(1983), (2) приближенный метод максимального
правдоподобия с итерациями назад, (3)точный метод
максимального правдоподобия по Meларду (1984).

Сравнение методов. В общем, все методы дают
очень похожие результаты. Также все методы
показали примерно одинаковую эффективность на
реальных данных. Однако метод 1 (см. выше) —
самый быстрый, и им можно пользоваться для
исследования очень длинных рядов (например,
содержащих более 30,000 наблюдений). Метод Меларда
(номер 3) может оказаться неэффективным, если
оцениваются параметры сезонной модели с большим
сезонным лагом (например, 365 дней). С другой
стороны, вы можете использовать вначале
приближенный метод максимального правдоподобия
(для того, чтобы найти прикидочные оценки
параметров), а затем точный метод; обычно
требуется только несколько итераций точного
метода (номер 3, выше), чтобы получить
окончательные оценки.

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

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

Оценивание модели

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

Другой критерий качества. Другой обычной
мерой надежности модели является сравнение
прогноза, построенного по урезанному ряду с
«известными (исходными) данными».

График прогноза

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

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

Ограничения. Следует напомнить, что модель
АРПСС является подходящей только для рядов,
которые являются стационарными
(среднее, дисперсия и автокорреляция примерно
постоянны во времени); для нестационарных рядов
следует брать разности. Рекомендуется иметь, как
минимум, 50 наблюдений в файле исходных данных.
Также предполагается, что параметры модели
постоянны, т.е. не меняются во времени.

Прерванные временные ряды

Обычный вопрос, возникающий при анализе
временных рядов, состоит в следующем,
воздействует или нет внешнее событие на
последовательность наблюдений. Например,
привела ли новая экономическая политика к росту
экономики, как обещалось; изменил ли новый закон
интенсивность преступлений и т.д. В общем, нужно
оценивать воздействия одного или нескольких
дискретных событий на значения ряда. Этот вид
анализа прерванных временных рядов подробно
описан в книге McDowall, McCleary, Meidinger, and Hay (1980).
Различают следующие три типа воздействий: (1)
устойчивое скачкообразное, (2) устойчивое
постепенное, (3) скачкообразное временное. См.
также следующие разделы:

  • Идентификация модели
    временных рядов
  • АРПСС
  • Экспоненциальное
    сглаживание
  • Сезонная декомпозиция (метод
    Census I)
  • Сезонная корректировка X-11
    (метод Census II)
  • Таблицы результатов
    корректировки X-11
  • Анализ распределенных лагов
  • Одномерный анализ Фурье
  • Кросс-спектральный анализ
  • Основные понятия и принципы
  • Быстрое преобразование Фурье

Экспоненциальное сглаживание

  • Общее введение
  • Простое экспоненциальное
    сглаживание
  • Выбор лучшего значения
    параметра a (альфа)
  • Индексы качества подгонки
  • Сезонная и несезонная модели
    с трендом или без тренда

См. также:

  • Идентификация модели
    временных рядов
  • АРПСС (Бокс и Дженкинс) и
    автокорреляции
  • Прерванные временные ряды
  • Сезонная декомпозиция (метод
    Census I)
  • Сезонная корректировка X-11
    (метод Census II)
  • Таблицы результатов
    корректировки X-11
  • Анализ распределенных лагов
  • Одномерный анализ Фурье
  • Кросс-спектральный анализ
  • Основные понятия и принципы
  • Быстрое преобразование Фурье

Общее введение

Экспоненциальное сглаживание — это очень
популярный метод прогнозирования многих
временных рядов. Исторически метод был
независимо открыт Броуном и Холтом. Броун служил
на флоте США во время второй мировой войны, где
занимался обнаружением подводных лодок и
системами наведения. Позже он применил открытый
им метод для прогнозирования спроса на запасные
части. Свои идеи он описал в книге, вышедшей в
свет в 1959 году. Исследования Холта были
поддержаны Департаментом военно-морского флота
США. Независимо друг от друга, Броун и Холт
открыли экспоненциальное сглаживание для
процессов с постоянным трендом, с линейным
трендом и для рядов с сезонной составляющей.

Gardner (1985), предложил «единую» классификацию
методов экспоненциального сглаживания.
Превосходное введение в эти методы можно найти в
книгах Makridakis, Wheelwright, and McGee (1983), Makridakis and Wheelwright (1989),
Montgomery, Johnson, and Gardiner (1990).

Простое экспоненциальное
сглаживание

Простая и прагматически ясная модель
временного ряда имеет следующий вид: Xt = b + t, где b
константа и
(эпсилон) — случайная ошибка. Константа b относительно
стабильна на каждом временном интервале, но
может также медленно изменяться со временем.
Один из интуитивно ясных способов выделения b
состоит в том, чтобы использовать сглаживание
скользящим средним, в котором последним
наблюдениям приписываются большие веса, чем
предпоследним, предпоследним большие веса, чем
пред-предпоследним и т.д. Простое
экспоненциальное именно так и устроено. Здесь
более старым наблюдениям приписываются
экспоненциально убывающие веса, при этом, в
отличие от скользящего среднего, учитываются все
предшествующие наблюдения ряда, а не те, что
попали в определенное окно. Точная формула
простого экспоненциального сглаживания имеет
следующий вид:

St = *Xt + (1-)*St-1

Когда эта формула применяется рекурсивно, то
каждое новое сглаженное значение (которое
является также прогнозом) вычисляется как
взвешенное среднее текущего наблюдения и
сглаженного ряда. Очевидно, результат
сглаживания зависит от параметра (альфа). Если равно 1, то
предыдущие наблюдения полностью игнорируются.
Если равно 0, то
игнорируются текущие наблюдения. Значения между 0, 1 дают
промежуточные результаты.

Эмпирические исследования Makridakis и др. (1982;
Makridakis, 1983) показали, что весьма часто простое
экспоненциальное сглаживание дает достаточно
точный прогноз.

Выбор лучшего значения
параметра (альфа)

Gardner (1985) обсуждает различные теоретические и
эмпирические аргументы в пользу выбора
определенного параметра сглаживания. Очевидно,
из формулы, приведенной выше, следует, что должно попадать в интервал между 0
(нулем) и 1 (хотя Brenner et al., 1968, для дальнейшего
применения анализа АРПСС считают, что 0<<2). Gardner (1985)
сообщает, что на практике обычно рекомендуется
брать меньше .30.
Однако в исследовании Makridakis et al., (1982), большее .30, часто
дает лучший прогноз. После обзора литературы,
Gardner (1985) приходит к выводу, что лучше оценивать
оптимально по данным
(см. ниже), чем просто «гадать» или
использовать искусственные рекомендации.

Оценивание лучшего значения с помощью данных. На
практике параметр сглаживания часто ищется с поиском
на сетке
. Возможные значения параметра
разбиваются сеткой с определенным шагом.
Например, рассматривается сетка значений от = 0.1 до = 0.9, с шагом 0.1.
Затем выбирается ,
для которого сумма квадратов (или средних
квадратов) остатков (наблюдаемые значения минус
прогнозы на шаг вперед) является минимальной.

Индексы качества подгонки

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

График прогноза

Такая визуальная проверка точности прогноза
часто дает наилучшие результаты. Имеются также
другие меры ошибки, которые можно использовать
для определения оптимального параметра (см. Makridakis, Wheelwright, and McGee,
1983):

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

Средняя абсолютная ошибка. Средняя
абсолютная ошибка (САО) вычисляется как среднее абсолютных
ошибок. Если она равна 0 (нулю), то имеем
совершенную подгонку (прогноз). В сравнении со
средней квадратической ошибкой, эта мера
«не придает слишком большого значения»
выбросам.

Сумма квадратов ошибок (SSE),
среднеквадратическая ошибка.
Эти величины
вычисляются как сумма (или среднее) квадратов
ошибок. Это наиболее часто используемые индексы
качества подгонки.

Относительная ошибка (ОО). Во всех
предыдущих мерах использовались действительные
значения ошибок. Представляется естественным
выразить индексы качества подгонки в терминах относительных
ошибок. Например, при прогнозе месячных продаж,
которые могут сильно флуктуировать (например, по
сезонам) из месяца в месяц, вы можете быть вполне
удовлетворены прогнозом, если он имеет точность
?10%. Иными словами, при прогнозировании
абсолютная ошибка может быть не так интересна
как относительная. Чтобы учесть относительную
ошибку, было предложено несколько различных
индексов (см. Makridakis, Wheelwright, and McGee, 1983). В первом
относительная ошибка вычисляется как:

ООt = 100*(Xt — Ft )/Xt

где Xt — наблюдаемое
значение в момент времени t, и Ft — прогноз (сглаженное
значение).

Средняя относительная ошибка (СОО).
Это значение вычисляется как среднее
относительных ошибок.

Средняя абсолютная относительная ошибка
(САОО).
Как и в случае с обычной средней
ошибкой отрицательные и положительные
относительные ошибки будут подавлять друг друга.
Поэтому для оценки качества подгонки в целом (для
всего ряда) лучше использовать среднюю абсолютную
относительную ошибку. Часто эта мера более
выразительная, чем среднеквадратическая ошибка.
Например, знание того, что точность прогноза ±5%,
полезно само по себе, в то время как значение 30.8
для средней квадратической ошибки не может быть
так просто проинтерпретировано.

Автоматический поиск лучшего параметра.
Для минимизации средней квадратической ошибки,
средней абсолютной ошибки или средней
абсолютной относительной ошибки используется
квази-ньютоновская процедура (та же, что и в АРПСС). В большинстве случаев
эта процедура более эффективна, чем обычный
перебор на сетке (особенно, если параметров
сглаживания несколько), и оптимальное значение alphanav.gif (845 bytes) можно быстро
найти.

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

Сезонная и несезонная модели с
трендом или без тренда

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

Аддитивная и мультипликативная
сезонность.
Многие временные ряды имеют
сезонные компоненты. Например, продажи игрушек
имеют пики в ноябре, декабре и, возможно, летом,
когда дети находятся на отдыхе. Эта
периодичность имеет место каждый год. Однако
относительный размер продаж может слегка
изменяться из года в год. Таким образом, имеет
смысл независимо экспоненциально сгладить
сезонную компоненту с дополнительным
параметром, обычно обозначаемым как (дельта). Сезонные
компоненты, по природе своей, могут быть
аддитивными или мультипликативными. Например, в
течение декабря продажи определенного вида
игрушек увеличиваются на 1 миллион долларов
каждый год. Для того чтобы учесть сезонное
колебание, вы можете добавить в прогноз на каждый
декабрь 1 миллион долларов (сверх
соответствующего годового среднего). В этом
случае сезонность — аддитивная. Альтернативно,
пусть в декабре продажи увеличились на 40%, т.е. в 1.4
раза. Тогда, если общие продажи малы, то
абсолютное (в долларах) увеличение продаж в
декабре тоже относительно мало (процент роста
константа). Если в целом продажи большие, то
абсолютное (в долларах) увеличение продаж будет
пропорционально больше. Снова, в этом случае
продажи увеличатся в определенное число раз, и
сезонность будет мультипликативной (в данном
случае мультипликативная сезонная составляющая
была бы равна 1.4). На графике различие между двумя
видами сезонности состоит в том, что в аддитивной
модели сезонные флуктуации не зависят от
значений ряда, тогда как в мультипликативной
модели величина сезонных флуктуаций зависит от
значений временного ряда.

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

Аддитивная модель:

Прогнозt = St + It-p

Мультипликативная модель:

Прогнозt = St*It-p

В этой формуле St
обозначает (простое) экспоненциально сглаженное
значение ряда в момент t, и It-p обозначает сглаженный
сезонный фактор в момент t
минус p (p
длина сезона). Таким образом, в сравнении с
простым экспоненциальным сглаживанием, прогноз
«улучшается» добавлением или умножением
сезонной компоненты. Эта компонента оценивается
независимо с помощью простого экспоненциального
сглаживания следующим образом:

Аддитивная модель:

It = It-p + *(1-)*et

Мультипликативная модель:

It = It-p + *(1-)*et/St

Обратите внимание, что предсказанная сезонная
компонента в момент t
вычисляется, как соответствующая компонента на
последнем сезонном цикле плюс ошибка (et,
наблюдаемое минус прогнозируемое значение в
момент t). Ясно, что параметр принимает значения
между 0 и 1. Если он равен нулю, то сезонная
составляющая на следующем цикле та же, что и на
предыдущем. Если
равен 1, то сезонная составляющая
«максимально» меняется на каждом шаге из-за
соответствующей ошибки (множитель  (1-) не
рассматривается из-за краткости введения). В
большинстве случаев, когда сезонность
присутствует, оптимальное значение лежит между 0 и 1.

Линейный, экспоненциальный,
демпфированный тренд.
Возвращаясь к примеру
с игрушками, мы можем увидеть наличие линейного
тренда (например, каждый год продажи
увеличивались на 1 миллион), экспоненциального
(например, каждый год продажи возрастают в 1.3
раза) или демпфированного тренда (в первом году
продажи возросли на 1 миллион долларов; во втором
увеличение составило только 80% по сравнению с
предыдущим, т.е. на $800,000; в следующем году вновь
увеличение было только на 80%, т.е. на $800,000 * .8 = $640,000
и т.д.). Каждый тип тренда по-своему проявляется в
данных. В целом изменение тренда — медленное в
течение времени, и опять (как и сезонную
компоненту) имеет смысл экспоненциально
сгладить его с отдельным параметром
[обозначаемым (гамма)
— для линейного и экспоненциального тренда, (фи) — для
демпфированного тренда].

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


Сезонная декомпозиция (метод Census I)

  • Общее введение
  • Вычисления

См. также:

  • Идентификация модели
    временных рядов
  • АРПСС (Бокс и Дженкинс) и
    автокорреляции
  • Прерванные временные ряды
  • Экспоненциальное
    сглаживание
  • Сезонная корректировка X-11
    (метод Census II)
  • Таблицы результатов
    корректировки X-11
  • Анализ распределенных лагов
  • Одномерный анализ Фурье
  • Кросс-спектральный анализ
  • Основные понятия и принципы
  • Быстрое преобразование Фурье

Общее введение

Предположим, что у вас имеются ежемесячные
данные о пассажиропотоке на международных
авиалиниях за 12 лет (см. Бокс и Дженкинс, 1976). Если
изобразить эти данные на графике, то будет хорошо
видно, что (1) объем пассажиропотока имеет во
времени возрастающий линейный тренд, и (2) в ряде
имеется ежегодно повторяющаяся закономерность — сезонность
(большинство перевозок приходится на летние
месяцы, кроме того, имеется пик меньшей высоты в
районе декабрьских каникул). Цель сезонной
декомпозиции и корректировки как раз и состоит в
том, чтобы отделить эти компоненты, то есть
разложить ряд на составляющую тренда, сезонную
компоненту и оставшуюся нерегулярную
составляющую. «Классический» прием,
позволяющий выполнить такую декомпозицию,
известен как метод Census I. Этот метод
описывается и обсуждается в работах Makridakis,
Wheelwright, and McGee (1983) и Makridakis and Wheelwright (1989).

Общая модель. Основная идея сезонной
декомпозиции проста. В общем случае временной
ряд типа того, который описан выше, можно
представить себе состоящим из четырех различных
компонент: (1) сезонной компоненты (обозначается St,
где t обозначает момент времени), (2) тренда (Tt),
(3) циклической компоненты (Ct) и (4)
случайной, нерегулярной компоненты или
флуктуации (It). Разница между
циклической и сезонной компонентой состоит в
том, что последняя имеет регулярную (сезонную)
периодичность, тогда как циклические факторы
обычно имеют более длительный эффект, который к
тому же меняется от цикла к циклу. В методе Census I
тренд и циклическую компоненту обычно
объединяют в одну тренд-циклическую компоненту
(TCt). Конкретные функциональные
взаимосвязи между этими компонентами могут
иметь самый разный вид. Однако, можно выделить
два основных способа, с помощью которых они могут
взаимодействовать: аддитивно и мультипликативно:

Аддитивная модель:

Xt = TCt + St + It

Мультипликативная модель:

Xt = Tt*Ct*St*It

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

Аддитивная и мультипликативная
сезонность.
Рассмотрим на примере различие
между аддитивной и мультипликативной сезонными
компонентами. График объема продаж детских
игрушек, вероятно, будет иметь ежегодный пик в
ноябре-декабре, и другой — существенно меньший по
высоте — в летние месяцы, приходящийся на
каникулы. Такая сезонная закономерность будет
повторяться каждый год. По своей природе
сезонная компонента может быть аддитивной или
мультипликативной. Так, например, каждый год
объем продаж некоторой конкретной игрушки может
увеличиваться в декабре на 3 миллиона долларов.
Поэтому вы можете учесть эти сезонные изменения, прибавляя
к своему прогнозу на декабрь 3 миллиона. Здесь
мы имеем аддитивную сезонность. Может
получиться иначе. В декабре объем продаж
некоторой игрушки может увеличиваться на 40%, то
есть умножаться на множитель 1.4. Это значит,
например, что если средний объем продаж этой
игрушки невелик, то абсолютное (в денежном
выражении) увеличение этого объема в декабре
также будет относительно небольшим (но в
процентном исчислении оно будет постоянным);
если же игрушка продается хорошо, то и абсолютный
(в долларах) рост объема продаж будет
значительным. Здесь опять, объем продаж
возрастает в число раз, равное определенному множителю,
а сезонная компонента, по своей природе,
мультипликативная
компонента (в данном случае
равная 1.4). Если перейти к графикам временных
рядов, то различие между этими двумя видами
сезонности будет проявляться так: в аддитивном
случае ряд будет иметь постоянные сезонные
колебания, величина которых не зависит от общего
уровня значений ряда; в мультипликативном случае
величина сезонных колебаний будет меняться в
зависимости от общего уровня значений ряда.

Аддитивный и мультипликативный тренд-цикл. Рассмотренный
пример можно расширить, чтобы проиллюстрировать
понятия аддитивной и мультипликативной
тренд-циклических компонент. В случае с
игрушками, тренд «моды» может привести к
устойчивому росту продаж (например, это может
быть общий тренд в сторону игрушек
образовательной направленности). Как и сезонная
компонента, этот тренд может быть по своей
природе аддитивным (продажи ежегодно
увеличиваются на 3 миллиона долларов) или
мультипликативным (продажи ежегодно
увеличиваются на 30%, или возрастают в 1.3 раза).
Кроме того, объем продаж может содержать
циклические компоненты. Повторим еще раз, что
циклическая компонента отличается от сезонной
тем, что она обычно имеет большую временную
протяженность и проявляется через неравные
промежутки времени. Так, например, некоторая
игрушка может быть особенно «горячей» в
течение летнего сезона (например, кукла,
изображающая персонаж популярного мультфильма,
которая к тому же агрессивно рекламируется). Как
и в предыдущих случаях, такая циклическая
компонента может изменять объем продаж
аддитивно, либо мультипликативно.

Вычисления

В вычислительном отношении процедура метода Сезонной
декомпозиции (Census I)
следует стандартным
формулам, см. Makridakis, Wheelwright, and McGee (1983) или Makridakis and
Wheelwright (1989).

График

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

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

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

График

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

График

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

Тренд-циклическая компонента.
Напомним, что циклическая компонента отличается
от сезонной компоненты тем, что
продолжительность цикла, как правило, больше, чем
один сезонный период, и разные циклы могут иметь
разную продолжительность. Приближение для
объединенной тренд-циклической компоненты можно
получить, применяя к ряду с сезонной поправкой
процедуру 5-точечного (центрированного)
взвешенного скользящего среднего с весами 1, 2, 3, 2,
1.

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


Сезонная корректировка X-11 (метод Census II)

Общие идеи, лежащие в основе сезонной
декомпозиции и корректировки, изложены в
разделе, посвященном методу сезонной
корректировки Census I (см. Сезонная
декомпозиция (метод Census I)
). Метод Census II (2)
является развитием и уточнением обычного метода
корректировки. На протяжении многих лет
различные варианты метода Census II развивались в
Бюро Переписи США (US Census Bureau); один из вариантов
этого метода, получивший широкую известность и
наиболее часто применяемый в государственных
органах и сфере бизнеса, называется «вариант X-11
метода Census II» (см. Shiskin, Young, and Musgrave, 1967).
Впоследствии этот усовершенствованный вариант
метода Census II стал называться просто X-11.
Помимо документации, которую можно получить из
Census Bureau, подробное описание метода дано в работах
Makridakis, Wheelwright and McGee (1983), Makridakis and Wheelwright (1989).

За дополнительной информацией обратитесь к
следующим разделам:

  • Сезонная корректировка:
    основные идеи и термины
  • Метод Census II
  • Таблицы результатов
    корректировки X-11
  • Подробное описание всех
    таблиц результатов, вычисляемых в методе X-11

За дальнейшей информацией обратитесь к Анализу временных рядов и
следующим разделам:

  • Идентификация модели
    временных рядов
  • АРПСС (Бокс и Дженкинс) и
    автокорреляции
  • Прерванные временные ряды
  • Экспоненциальное
    сглаживание
  • Сезонная декомпозиция (метод
    Census I)
  • Анализ распределенных лагов
  • Одномерный анализ Фурье
  • Кросс-спектральный анализ
  • Основные понятия и принципы
  • Быстрое преобразование Фурье

Сезонная корректировка: основные идеи и
термины

Предположим, что у вас имеются ежемесячные
данные о пассажиропотоке на международных
авиалиниях за 12 лет (см. Бокс и Дженкинс, 1976). Если
изобразить эти данные на графике, то будет хорошо
видно, что (1) объем пассажиропотока имеет во
времени возрастающий линейный тренд, и что (2) в
ряде имеется ежегодно повторяющаяся
закономерность — сезонность (большинство
перевозок приходится на летние месяцы, кроме
того, имеется пик меньшей высоты в районе
декабрьских каникул). Цель сезонной декомпозиции
и корректировки как раз и состоит в том, чтобы
отделить эти компоненты, то есть разложить ряд на
составляющую тренда, сезонную компоненту и
оставшуюся нерегулярную составляющую.
«Классический» прием, позволяющий выполнить
такую декомпозицию, известен как метод Census I
(см. раздел Census I). Этот метод
описывается и обсуждается в работах Makridakis,
Wheelwright, and McGee (1983) и Makridakis and Wheelwright (1989).

Общая модель. Основная идея сезонной
декомпозиции проста. В общем случае временной
ряд типа того, который описан выше, можно
представить себе состоящим из четырех различных
компонент: (1) сезонной компоненты (обозначается St,
где t обозначает момент времени), (2) тренда (Tt),
(3) циклической компоненты (Ct) и (4)
случайной, нерегулярной компоненты или
флуктуации (It). Разница между
циклической и сезонной компонентой состоит в
том, что последняя имеет регулярную (сезонную)
периодичность, тогда как циклические факторы
обычно имеют более длительный эффект, который к
тому же меняется от цикла к циклу. В методе Census I
тренд и циклическую компоненту обычно
объединяют в одну тренд-циклическую компоненту
(TCt). Конкретные функциональные
взаимосвязи между этими компонентами могут
иметь самый разный вид. Однако, можно выделить
два основных способа, с помощью которых они могут
взаимодействовать: аддитивно и мультипликативно:

Аддитивная модель:

Xt = TCt + St + It

Мультипликативная модель:

Xt = Tt*Ct*St*It

Здесь Xt обозначает
значение временного ряда в момент времени t.

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

Аддитивная и мультипликативная
сезонность.
Рассмотрим на примере различие
между аддитивной и мультипликативной сезонными
компонентами. График объема продаж детских
игрушек, вероятно, будет иметь ежегодный пик в
ноябре-декабре, и другой — существенно меньший по
высоте — в летние месяцы, приходящийся на
каникулы. Такая сезонная закономерность будет
повторяться каждый год. По своей природе
сезонная компонента может быть аддитивной или
мультипликативной. Так, например, каждый год
объем продаж некоторой конкретной игрушки может
увеличиваться в декабре на 3 миллиона долларов.
Поэтому вы можете учесть эти сезонные изменения, прибавляя
к своему прогнозу на декабрь 3 миллиона. Здесь
мы имеем аддитивную сезонность. Может
получиться иначе. В декабре объем продаж
некоторой игрушки может увеличиваться на 40%, то
есть умножаться на множитель 1.4. Это значит,
например, что если средний объем продаж этой
игрушки невелик, то абсолютное (в денежном
выражении) увеличение этого объема в декабре
также будет относительно небольшим (но в
процентном исчислении оно будет постоянным);
если же игрушка продается хорошо, то и абсолютный
(в долларах) рост объема продаж будет
значительным. Здесь опять, объем продаж
возрастает в число раз, равное определенному множителю,
а сезонная компонента, по своей природе,
мультипликативная
компонента (в данном случае
равная 1.4). Если перейти к графикам временных
рядов, то различие между этими двумя видами
сезонности будет проявляться так: в аддитивном
случае ряд будет иметь постоянные сезонные
колебания, величина которых не зависит от общего
уровня значений ряда; в мультипликативном случае
величина сезонных колебаний будет меняться в
зависимости от общего уровня значений ряда.

Аддитивный и мультипликативный тренд-цикл.
Рассмотренный пример можно расширить, чтобы
проиллюстрировать понятия аддитивной и
мультипликативной тренд-циклических компонент.
В случае с игрушками, тренд «моды» может
привести к устойчивому росту продаж (например,
это может быть общий тренд в сторону игрушек
образовательной направленности). Как и сезонная
компонента, этот тренд может быть по своей
природе аддитивным (продажи ежегодно
увеличиваются на 3 миллиона долларов) или
мультипликативным (продажи ежегодно
увеличиваются на 30%, или возрастают в 1.3 раза).
Кроме того, объем продаж может содержать
циклические компоненты. Повторим еще раз, что
циклическая компонента отличается от сезонной
тем, что она обычно имеет большую временную
протяженность и проявляется через неравные
промежутки времени. Так, например, некоторая
игрушка может быть особенно «горячей» в
течение летнего сезона (например, кукла,
изображающая персонаж популярного мультфильма,
которая к тому же агрессивно рекламируется). Как
и в предыдущих случаях, такая циклическая
компонента может изменять объем продаж
аддитивно, либо мультипликативно.

Метод Census II

Основной метод сезонной декомпозиции и
корректировки, рассмотренный в разделе Сезонная корректировка:
основные идеи и термины, может быть
усовершенствован различными способами. На самом
деле, в отличие от многих методов моделирования
временных рядов (в частности, АРПСС),
которые основаны на определенной теоретической
модели, вариант X-11 метода Census II представляет
собой просто результат многочисленных специально
разработанных
приемов и усовершенствований,
которые доказали свою работоспособность в
многолетней практике решения реальных задач (см.
Burman, 1979, Kendall and Ord, 1990, Makridakis and Wheelwright, 1989; Wallis, 1974).
Некоторые из наиболее важных усовершенствований
перечислены ниже.

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

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

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

Критерии и итоговые статистики. Помимо
оценки основных компонент ряда, можно вычислить
различные сводные статистики. Например, можно
сформировать таблицы дисперсионного анализа для
проверки значимости фактора сезонной
изменчивости и ряда и фактора рабочих дней (см.
выше), процедура метода X-11 вычисляет также
ежемесячные относительные изменения в случайной
и тренд-циклической компонентах. С увеличением
продолжительности временного промежутка,
измеряемого в месяцах или, в случае квартального
варианта метода X-11 — в кварталах года,
изменения в тренд-циклической компоненте, вообще
говоря, будут нарастать, в то время как изменения
случайной составляющей должны оставаться
примерно на одном уровне. Средняя длина
временного интервала, на котором изменения
тренд-циклической компоненты становятся
примерно равными изменениям случайной
компоненты, называется месяцем (кварталом)
циклического доминирования
, или сокращенно
МЦД (соответственно КЦД). Например, если МЦД равно
двум, то на сроках более двух месяцев
тренд-циклическая компонента станет
доминировать над флуктуациями нерегулярной
(случайной) компоненты. Эти и другие результаты
более подробно будут обсуждаться далее.

Таблицы результатов
корректировки X-11

Вычисления, которые производятся в процедуре
  X-11, лучше всего обсуждать в контексте
таблиц результатов, которые при этом выдаются.
Процедура корректировки разбивается на семь
этапов, которые обычно обозначаются буквами A
G.

  1. Априорная корректировка (помесячная сезонная
    корректировка).
    Перед тем, как к временному
    ряду, содержащему ежемесячные значения, будет
    применяться какая-либо сезонная корректировка,
    могут быть произведены различные корректировки,
    заданные пользователем. Можно ввести еще один
    временной ряд, содержащий априорные
    корректирующие факторы; значения этого ряда
    будут вычитаться из исходного ряда (аддитивная
    модель), или же значения исходного ряда будут
    поделены на значения корректирующего ряда
    (мультипликативная модель). В случае
    мультипликативной модели пользователь может
    также определить свои собственные поправочные
    коэффициенты (веса) на число рабочих дней. Эти
    веса будут использоваться для корректировки
    ежемесячных наблюдений, так чтобы учитывалось
    число рабочих дней в этом месяце.
  2. Предварительное оценивание вариации числа
    рабочих дней (месячный вариант X-11) и весов.
    На
    следующем шаге вычисляются предварительные
    поправочные коэффициенты на число рабочих дней
    (только в месячном варианте X-11) и веса,
    позволяющие уменьшить эффект выбросов.
  3. Окончательное оценивание вариации числа
    рабочих дней и нерегулярных весов (месячный
    вариант X-11).
    Поправки и веса, вычисленные в
    пункте B, используются для построения
    улучшенных оценок тренд-циклической и сезонной
    компонент. Эти улучшенные оценки используются
    для окончательного вычисления факторов числа
    рабочих дней (в месячном варианте X-11) и весов.
  4. Окончательное оценивание сезонных факторов,
    тренд-циклической, нерегулярной и сезонно
    скорректированной компонент ряда.
    Окончательные
    значения факторов рабочих дней и весов,
    вычисленные в пункте C, используются для
    вычисления окончательных оценок для компонент
    ряда.
  5. Модифицированные ряды: исходный, сезонно
    скорректированный и нерегулярный.
    Исходный и
    окончательный сезонно скорректированный ряды, а
    также нерегулярная компонента модифицируются
    путем сглаживания выбросов. Полученные в
    результате этого, модифицированные ряды
    позволяют пользователю проверить устойчивость
    сезонной корректировки.
  6. Месяц (квартал) циклического доминирования
    (МЦД, КЦД), скользящее среднее и сводные
    показатели.
    IНа этом этапе вычислений
    рассчитываются различные сводные
    характеристики (см. далее), позволяющие
    пользователю исследовать относительную
    важность разных компонент, среднюю флуктуацию от
    месяца к месяцу (от квартала к кварталу), среднее
    число идущих подряд изменений в одну сторону и
    др.
  7. Графики. Наконец, вы можете построить
    различные графики итоговых результатов.
    Например, можно построить окончательно
    скорректированный ряд в хронологическом порядке
    или по месяцам (см. ниже).

Подробное описание всех таблиц
результатов, вычисляемых в методе X-11

На каждом из этапов AG (см. раздел Таблицы результатов
корректировки X-11) вычислялись различные
таблицы результатов. Обычно все они нумеруются, а
также им приписывается буква, соответствующая
этапу анализа. Например, таблица B 11 содержит
предварительно сезонно скорректированный ряд; C
11
— это более точно сезонно скорректированный
ряд, а D 11 — окончательный сезонно
скорректированный ряд. Далее приводится
перечень всех таблиц. Таблицы, помеченные
звездочкой (*), недоступны (или неприменимы) при
анализе квартальных показателей. Кроме того, в
случае квартальной корректировки некоторые из
описанных ниже вычислений несколько
видоизменяются. Так, например, для вычисления
сезонных факторов вместо 12-периодного (т.е.
12-месячного) скользящего среднего используется
4-периодное (4-квартальное) скользящее среднее;
предварительная тренд-циклическая компонента
вычисляется по центрированному 4-периодному
скользящему среднему, а окончательная оценка
тренд-циклической компоненты вычисляется по
5-точечному среднему Хендерсона.

В соответствии со стандартом метода X-11,
принятым Бюро переписи США, предусмотрены три
степени подробности вывода: Стандартный (17 — 27
таблиц), Длинный (27 — 39 таблиц) и Полный (44 —
59 таблиц). Имеется также возможность выводить
только таблицы результатов, выбранные
пользователем. В следующих далее описаниях
таблиц, буквы С, Д и П рядом с
названием таблицы указывают, какие таблицы
выводятся и/или распечатываются в
соответствующем варианте вывода. (Для графиков
предусмотрены два уровня подробности вывода: Стандартный
и Все.)

Щелкните на имени таблицы для
получения информации о ней.

* A 1.
Исходный ряд (С)
*
A 2. Априорные месячные поправки (С)
*
A 3. Исходный ряд, скорректированный с помощью
априорных месячных поправок (С)
* A
4. Априорные поправки на рабочие дни (С)
B
1. Ряд после априорной корректировки либо
исходный ряд (С)
B 2.
Тренд-цикл (Д)
B
3. Немодифицированные S-I разности или отношения (П)
B
4. Значения для замены выбросов S-I разностей
(отношений) (П)
B 5.
Сезонная составляющая (П)
B 6.
Сезонная корректировка ряда (П)
B 7.
Тренд-цикл (Д)
B 8.
Немодифицированные S-I разности (отношения) (П)
B
9. Значения для замены выбросов S-I разностей
(отношений) (П)
B 10.
Сезонная составляющая (Д)
B 11.
Сезонная корректировка ряда (П)
B 12. (не используется)
B 13.
Нерегулярная составляющая ряда (Д)
Таблицы B 14 — B 16, B 18 и B 19: Поправка на число
рабочих дней.
Эти таблицы доступны только при
анализе ежемесячных данных. Число разных дней
недели (понедельников, вторников и т.д.)
колеблется от месяца к месяцу. Бывают ряды, в
которых различия в числе рабочих дней в месяце
могут давать заметный разброс ежемесячных
показателей (например, месячный доход парка
аттракционов сильно зависит от того, сколько в
этом месяце было выходных дней). Пользователь
имеет возможность определить начальные веса для
каждого дня недели (см. A
4
), и/или эти веса могут быть оценены по данным
(пользователь также может сделать использование
этих весов условным, т.е. только в тех случаях,
когда они объясняют значительную часть
дисперсии).
*
B 14. Выбросы нерегулярной составляющей,
исключенные из регрессии рабочих дней (Д)
* B
15. Предварительная регрессия рабочих дней (Д)
*
B 16. Поправки на число рабочих дней, полученные из
коэффициентов регрессии (П)
B
17. Предварительные веса нерегулярной компоненты (Д)
*
B 18. Поправки на число рабочих дней, полученные из
комбинированных весов дней недели (П)
*
B 19. Исходный ряд с поправками на рабочие дни и
априорную вариацию (П)
C
1. Исходный ряд, модифицированный с помощью
предварительных весов, с поправкой на рабочие
дни и априорную вариацию (Д)
C 2.
Тренд-цикл (П)
C 3. (не используется)
C
4. Модифицированные S-I разности (отношения) (П)
C 5.
Сезонная составляющая (П)
C 6.
Сезонная корректировка ряда (П)
C 7.
Тренд-цикл (Д)
C 8. (не используется)
C
9. Модифицированные S-I разности (отношения) (П)
C 10.
Сезонная составляющая (Д)
C 11.
Сезонная корректировка ряда (П)
C 12. (не используется)
C 13.
Нерегулярная составляющая (С)
Таблицы C 14 — C 16, C 18 и C 19: Поправка на число
рабочих дней.
Эти таблицы доступны только при
анализе ежемесячных данных и если при этом
требуется поправка на различное число рабочих
дней. В этом случае поправки на число рабочих
дней вычисляются по уточненным значениям
сезонно скорректированных рядов аналогично
тому, как это делалось в пункте B (B 14 B 16,
B 18, B 19).
*
C 14. Выбросы нерегулярной составляющей,
исключенные из регрессии рабочих дней (С)
* C
15. Регрессия рабочих дней — окончательный вариант
(С)
*
C 16. Поправки на число рабочих дней, полученные из
коэффициентов регрессии, — окончательный вариант
(С)
C
17. Окончательные веса нерегулярной компоненты (С)
*
C 18. Поправки на число рабочих дней, полученные из
комбинированных весов дней недели —
окончательный вариант (С)
*
C 19. Исходный ряд с поправками на рабочие дни и
априорную вариацию (С)
D
1. Исходный ряд, модифицированный с помощью
окончательных весов, с поправкой на рабочие дни и
априорную вариацию (Д)
D 2.
Тренд-цикл (П)
D 3. (не используется)
D
4. Модифицированные S-I разности (отношения) (П)
D 5.
Сезонная составляющая (П)
D 6.
Сезонная корректировка ряда (П)
D 7.
Тренд-цикл (Д)
D
8. Немодифицированные S-I разности (отношения) —
окончательный вариант (С)
D
9. Окончательные значения для замены выбросов S-I
разностей (отношений) (С)
D 10.
Сезонная составляющая — окончательный вариант (С)
D
11. Сезонная корректировка ряда — окончательный
вариант (С)
D 12.
Тренд-циклическая компонента — окончательный
вариант (С)
D 13.
Нерегулярная составляющая — окончательный
вариант (С)
E 1.
Модифицированный исходный ряд (С)
E 2.
Модифицированный ряд с сезонной поправкой (С)
E 3.
Модифицированная нерегулярная составляющая (С)
E
4. Разности (отношения) годовых сумм (С)
E
5. Разности (относительные изменения) исходного
ряда (С)
E
6. Разности (относительные изменения)
окончательного варианта ряда с сезонной
поправкой (С)
F 1.
МЦД (КЦД) скользящее среднее (С)
F 2.
Сводные показатели (С)
G 1. График (С)
G 2. График (С)
G 3. График (В)
G 4. График (В)

Анализ распределенных лагов

  • Общая цель
  • Общая модель
  • Распределенный лаг Алмона

За дальнейшей информацией обратитесь к Анализу временных рядов и
следующим разделам:

  • Идентификация модели
    временных рядов
  • АРПСС (Бокс и Дженкинс) и
    автокорреляции Вводный обзор АРПСС
  • Прерванные временные ряды
  • Экспоненциальное
    сглаживание
  • Сезонная декомпозиция (метод
    Census I)
  • Сезонная корректировка X-11
    (метод Census II)
  • Таблицы результатов
    корректировки X-11
  • Одномерный анализ Фурье
  • Кросс-спектральный анализ
  • Основные понятия и принципы
  • Быстрое преобразование Фурье

Общая цель

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

временной сдвиг (лаг) (см. также автокорреляции и
кросскорреляции).

Такого рода зависимости с запаздыванием
особенно часто возникают в эконометрике.
Например, доход от инвестиций в новое
оборудование отчетливо проявится не сразу, а
только через определенное время. Более высокий
доход изменяет выбор жилья людьми; однако эта
зависимость, очевидно, тоже проявляется с
запаздыванием. [Подобные задачи возникают в
страховании, где временной ряд клиентов и ряд
денежных поступлений сдвинуты друг относительно
друга].

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

Подробные обсуждения зависимостей с
распределенными лагами имеются в
эконометрических учебниках, например, в Judge, Griffith,
Hill, Luetkepohl, and Lee (1985), Maddala (1977), and Fomby, Hill, and Johnson (1984).
Ниже дается краткое описание этих методов.
Предполагается, что вы знакомы с понятием
корреляции (см. Основные
статистики и таблицы
), кросскорреляции и
основными идеями множественной регрессии (см. Множественная регрессия).

Общая модель

Пусть y — зависимая переменная, a независимая
или объясняющая x. Эти переменные измеряются
несколько раз в течение определенного отрезка
времени. В некоторых учебниках по эконометрике
зависимая переменная называется также эндогенной
переменной, a зависимая или объясняемая
переменная экзогенной переменной.
Простейший способ описать зависимость между
этими двумя переменными дает следующее линейное
уравнение:

Yt = i*xt-i

В этом уравнении значение зависимой переменной
в момент времени  t является
линейной функцией переменной x,
измеренной в моменты t, t-1,
t-2
и т.д. Таким образом, зависимая
переменная представляет собой линейные функции x и x, сдвинутых на
1, 2, и т.д. временные периоды. Бета коэффициенты
(i) могут
рассматриваться как параметры наклона в этом
уравнении. Будем рассматривать это уравнение как
специальный случай уравнения линейной регрессии
(см. раздел Множественная
регрессия
). Если коэффициент переменной с
определенным запаздыванием (лагом) значим, то
можно заключить, что переменная y
предсказывается (или объясняется) с
запаздыванием.

Распределенный лаг Алмона

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

Алмон (1965) предложил специальную процедуру,
которая в данном случае уменьшает
мультиколлинеарность. Именно, пусть каждый
неизвестный коэффициент записан в виде:

i =
0 + 1*i + … + q*iq

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

Неправильная спецификация. Общая
проблема полиномиальной аппроксимации, состоит
в том, что длина лага и степень полинома
неизвестны заранее. Последствия
неправильного определения (спецификации) этих
параметров потенциально серьезны (в силу
смещения, возникающего в оценках при
неправильном задании параметров). Этот вопрос
подробно обсуждается в книгах Frost (1975), Schmidt and Waud
(1973), Schmidt and Sickles (1975) и Trivedi and Pagan (1979).


Одномерный анализ Фурье

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

Наиболее известный пример применения
спектрального анализа — циклическая природа
солнечных пятен (например, см. Блумфилд, 1976 или
Шамвэй, 1988). Оказывается, что активность
солнечных пятен имеет 11-ти летний цикл. Другие
примеры небесных явлений, изменения погоды,
колебания в товарных ценах, экономическая
активность и т.д. также часто используются в
литературе для демонстрации этого метода. В
отличие от АРПСС или метода экспоненциального
сглаживания (см. разделы АРПСС
и Экспоненциальное
сглаживание), цель спектрального анализа —
распознать сезонные колебания различной длины, в
то время как в предшествующих типах анализа,
длина сезонных компонент обычно известна (или
предполагается) заранее и затем включается в
некоторые теоретические модели скользящего
среднего или автокорреляции.

Классический текст по спектральному анализу —
Bloomfield (1976); однако другие подробные обсуждения
могут быть найдены в Jenkins and Watts (1968), Brillinger (1975), Brigham
(1974), Elliott and Rao (1982), Priestley (1981), Shumway (1988) или Wei (1989).

За дальнейшей информацией обратитесь к Анализу временных рядов и
следующим разделам:

  • Основные понятия и принципы
  • Быстрое преобразование Фурье
  • Идентификация модели
    временных рядов
  • АРПСС (Бокс и Дженкинс) и
    автокорреляции Вводный обзор АРПСС
  • Прерванные временные ряды
  • Анализ распределенных лагов
  • Сезонная декомпозиция (метод
    Census I)
  • Экспоненциальное
    сглаживание
  • Кросс-спектральный анализ

Кросс-спектральный анализ

  • Общее введение
  • Основные понятия и принципы
  • Результаты для каждой
    переменной
  • Кросс-периодограмма,
    кросс-плотность, квадратурная плотность и
    кросс-амплитуда
  • Квадрат когерентности,
    усиление и фазовый сдвиг
  • Как создавались данные для
    примера

За дальнейшей информацией обратитесь к Анализу временных рядов и
следующим разделам:

  • Идентификация модели
    временных рядов
  • АРПСС (Бокс и Дженкинс) и
    автокорреляции Вводный обзор АРПСС
  • Прерванные временные ряды
  • Экспоненциальное
    сглаживание
  • Сезонная декомпозиция (метод
    Census I)
  • Сезонная корректировка X-11
    (метод Census II)
  • Таблицы результатов
    корректировки X-11
  • Анализ распределенных лагов
  • Одномерный анализ Фурье
  • Основные понятия и принципы
  • Быстрое преобразование Фурье

Общее введение

Кросс-спектральный анализ развивает Одномерный анализ Фурье и
позволяет анализировать одновременно два ряда.
Мы предполагаем, что вы уже прочитали введение к
разделу одномерного спектрального анализа.
Подробное обсуждение кросс-спектрального
анализа можно найти в книгах Bloomfield (1976), Jenkins and Watts
(1968), Brillinger (1975), Brigham (1974), Elliott and Rao (1982), Priestley (1981),
Shumway (1988), or Wei (1989).

Периодичность ряда на определенных
частотах.
Наиболее известный пример
применения спектрального анализа — циклическая
природа солнечных пятен (например, см. Блумфилд,
1976 или Шамвэй, 1988). Оказывается, что активность
солнечных пятен имеет 11-ти летний цикл. Другие
примеры небесных явлений, изменения погоды,
колебания в товарных ценах, экономическая
активность и т.д. также часто используются в
литературе для демонстрации этого метода.

Основные понятия и принципы

Простой пример. Рассмотрим следующие два
ряда с 16 наблюдениями:

  ПЕРЕМ1 ПЕРЕМ2
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
1.000
1.637
1.148
-.058
-.713
-.383
.006
-.483
-1.441
-1.637
-.707
.331
.441
-.058
-.006
.924
-.058
-.713
-.383
.006
-.483
-1.441
-1.637
-.707
.331
.441
-.058
-.006
.924
1.713
1.365
.266

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

Незавмсимая (X):
ПЕРЕМ1
Зависимая (Y): ПЕРЕМ2
 
Частота
 
Период
X
плотность
Y
плотность
Кросс
плотность
Кросс
квадр.
Кросс
амплит.
0.000000
.062500
.125000
.187500
.250000
.312500
.375000
.437500
.500000
 
16.00000
8.00000
5.33333
4.00000
3.20000
2.66667
2.28571
2.00000
.000000
8.094709
.058771
3.617294
.333005
.091897
.052575
.040248
.037115
.024292
7.798284
.100936
3.845154
.278685
.067630
.036056
.026633
0.000000
-.00000
2.35583
-.04755
-2.92645
-.26941
-.07435
-.04253
-.03256
0.00000
0.00000
-7.58781
.06059
2.31191
.14221
.02622
.00930
.00342
0.00000
.000000
7.945114
.077020
3.729484
.304637
.078835
.043539
.032740
0.000000

Результаты для каждой
переменной

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

Кросс-периодограмма,
кросс-плотность, квадратурная плотность и
кросс-амплитуда

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

Квадрат когерентности, усиление
и фазовый сдвиг

Существуют дополнительные статистики, которые
будут показаны в полной итоговой таблице
результатов.

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

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

Фазовый сдвиг. В заключение, оценки
фазового сдвига вычисляются как арктангенс (tan**-1)
коэффициента пропорциональности оценки
квадр-плотности и оценки кросс-плотности. Оценки
фазового сдвига (обычно обозначаемые греческой
буквой y) измеряют, насколько каждая частотная
компонента одного ряда опережает частотные
компоненты другого.

Как создавались данные для
примера

Теперь вернемся к примеру данных, приведенному
выше. Большие оценки спектральной плотности для
обоих рядов и значения кросс-амплитуды для
частот = 0.0625 и = .1875 предполагают две
существенных синхронных периодичности с этими
частотами в обоих рядах. Фактически, два ряда
создавались как:

v1 = cos(2**.0625*(v0-1))
+ .75*sin(2**.2*(v0-1))

v2 = cos(2**.0625*(v0+2)) +
.75*sin(2**.2*(v0+2))

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


Спектральный анализ — Основные понятия и
принципы

  • Частота и период
  • Общая структура модели
  • Простой пример
  • Периодограмма
  • Проблема рассеяния
  • Добавление констант во
    временной ряд (пэддинг)
  • Косинус-сглаживание
  • Окна данных и оценки
    спектральной плотности
  • Подготовка данных к анализу
  • Результаты для случая, когда в
    ряде отсутствует периодичность

За дальнейшей информацией обратитесь к Анализу временных рядов и
следующим разделам:

  • Идентификация модели
    временных рядов
  • АРПСС (Бокс и Дженкинс) и
    автокорреляции Вводный обзор АРПСС
  • Прерванные временные ряды
  • Экспоненциальное
    сглаживание
  • Сезонная декомпозиция (метод
    Census I)
  • Сезонная корректировка X-11
    (метод Census II)
  • Таблицы результатов
    корректировки X-11
  • Анализ распределенных лагов
  • Одномерный анализ Фурье
  • Кросс-спектральный анализ
  • Быстрое преобразование Фурье

Частота и период

Длина волны функций синуса или косинуса, как
правило, выражается числом циклов (периодов) в
единицу времени (Частота), часто обозначается
греческой буквой ню (; в некоторых учебниках также
используют f). Например, временной ряд,
состоящий из количества писем, обрабатываемых
почтой, может иметь 12 циклов в году. Первого числа
каждого месяца отправляется большое количество
корреспонденции (много счетов приходит именно
первого числа каждого месяца); затем, к середине
месяца, количество корреспонденции уменьшается;
и затем вновь возрастает к концу месяца. Поэтому
каждый месяц колебания в количестве
корреспонденции, обрабатываемой почтовым
отделением, будут проходить полный цикл. Таким
образом, если единица анализа — один год, то будет равно 12 (поскольку
имеется 12 циклов в году). Конечно, могут быть и
другие циклы с различными частотами. Например,
годичные циклы (=1)
и, возможно, недельные циклы (=52 недели в год).

Период Т функций синуса или косинуса
определяется как продолжительность по времени
полного цикла. Таким образом, это обратная
величина к частоте: T = 1/. Возвратимся к примеру с почтой из
предыдущего абзаца, здесь месячный цикл будет
равен 1/12 = 0.0833 года. Другими словами, это период
составляет 0.0833 года.

Общая структура модели

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

xt = a0 + [ak*cos(k*t) + bk*sin(k*t)]    (для k = 1 до q)

Следующее общее понятие классического
гармонического анализа в этом уравнении — (лямбда) -это
круговая частота, выраженная в радианах в
единицу времени, т.е. = 2**k, где константа пи =
3.1416 и k = k/q. Здесь
важно осознать, что вычислительная задача
подгонки функций синусов и косинусов разных длин
к данным может быть решена с помощью
множественной линейной регрессии. Заметим, что
коэффициенты ak при
косинусах и коэффициенты bk
при синусах — это коэффициенты регрессии,
показывающие степень, с которой соответствующие
функции коррелируют с данными [заметим, что сами
синусы и косинусы на различных частотах не
коррелированы или, другим языком, ортогональны.
Таким образом, мы имеем дело с частным случаем
разложения по ортогональным полиномам.] Всего
существует q различных синусов
и косинусов (см. также Множественная регрессия);
интуитивно ясно, что число функций синусов и
косинусов не может быть больше числа данных в
ряде. Не вдаваясь в подробности, отметим, если n —
количество данных, то будет n/2+1 функций
косинусов и n/2-1 функций синусов. Другими
словами, различных синусоидальных волн будет
столько же, сколько данных, и вы сможете
полностью воспроизвести ряд по основным
функциям. (Заметим, если количество данных в ряде
нечетно, то последнее наблюдение обычно
опускается. Для определения синусоидальной
функции нужно иметь, по крайней мере, две точки:
высокого и низкого пика.)

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

Комплексные числа (действительные и мнимые
числа).
Во многих учебниках по спектральному
анализу структурная модель, показанная выше,
представлена в комплексных числах; т.е. параметры
оцениваемого процесса описаны с помощью
действительной и мнимой части преобразования
Фурье. Комплексное число состоит из
действительного и мнимого числа. Мнимые числа, по
определению, — это числа, умноженные на константу i,
где i определяется как квадратный корень из -1.
Очевидно, корень квадратный из -1 не существует в
обычном сознании (отсюда термин мнимое число);
однако арифметические операции над мнимыми
числами могут производиться естественным
образом [например, (i*2)**2= -4]. Полезно представление
действительных и мнимых чисел, образующих
двумерную координатную плоскость, где
горизонтальная или X-ось представляет все
действительные числа, а вертикальная или Y-ось
представляет все мнимые числа. Комплексные числа
могут быть представлены точками на двумерной
плоскости. Например, комплексное число 3+i*2 может
быть представлено точкой с координатами {3,2} на
этой плоскости. Можно также представить
комплексные числа как углы; например, можно
соединить точку, соответствующую комплексному
числу на плоскости с началом координат
(комплексное число 0+i*0), и измерить угол наклона
этого вектора к горизонтальной оси. Таким
образом интуитивно ясно, каким образом формула
спектрального разложения, показанная выше, может
быть переписана в комплексной области. В таком
виде математические вычисления часто более
изящны и проще в выполнении, поэтому многие
учебники предпочитают представление
спектрального анализа в комплексных числах.

Простой пример

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

x = 1*cos(2**.0625*(v0-1))
+ .75*sin(2**.2*(v0-1))

Эта переменная состоит из двух основных
периодичностей — первая с частотой =.0625 (или периодом 1/=16; одно наблюдение
составляет 1/16-ю длины полного цикла, или весь
цикл содержит каждые 16 наблюдений) и вторая с
частотой =.2 (или
периодом 5). Коэффициент при косинусе (1.0) больше
чем коэффициент при синусе (.75). Итоговая таблица
результатов спектрального анализа показана
ниже.

  Спектральный
анализ: ПЕРЕМ1 (shumex.sta)
Число наблюдений: 16
 
t
Час-
тота
 
Период
Косинус
корэфф.
Синус
корэфф.
Периодо-
грамма
0
1
2
3
4
5
6
7
8
.0000
.0625
.1250
.1875
.2500
.3125
.3750
.4375
.5000
 
16.00
8.00
5.33
4.00
3.20
2.67
2.29
2.00
.000
1.006
.033
.374
-.144
-.089
-.075
-.070
-.068
0.000
.028
.079
.559
-.144
-.060
-.031
-.014
0.000
.000
8.095
.059
3.617
.333
.092
.053
.040
.037

Теперь рассмотрим столбцы таблицы результатов.
Ясно, что наибольший коэффициент при косинусах
расположен напротив частоты .0625. Наибольший
коэффициент при синусах соответствует частоте
.1875. Таким образом, эти две частоты, которые были
«внесены» в данные, отчетливо проявились.

Периодограмма

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

Pk = синус-коэффициентk2
+ косинус-коэффициентk2 * N/2

где Pk — значения
периодограммы на частоте  k , и N — общая
длина ряда. Значения периодограммы можно
интерпретировать как дисперсию (вариацию) данных
на соответствующей частоте. Обычно значения
периодограммы изображаются в зависимости от
частот или периодов.

График периодограммы

Проблема рассеяния

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

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

Ниже смотрите описание каждого из этих
подходов.

Добавление констант во
временной ряд (пэддинг)

Так как частотные величины вычисляются как N/t,
можно просто добавить в ряд константы (например,
нули), и таким образом получить увеличение
частот. Фактически, если вы добавите в файл
данных, описанный в примере выше, десять нулей,
результаты не изменятся; т.е. наибольшие пики
периодограммы будут находиться по-прежнему на
частотах близких к .0625 и .2. (Добавление констант
во временной ряд также часто желательно для
увеличения вычислительной эффективности; см.
ниже.)

Косинус-сглаживание

Так называемый процесс косинус-сглаживания
рекомендуемое преобразование ряда,
предшествующее спектральному анализу. Оно
обычно приводит к уменьшению рассеяния в
периодограмме. Логическое обоснование этого
преобразования подробно объясняется в книге
Bloomfield (1976, стр. 80-94). По существу, количественное
отношение (p) данных в начале и в конце ряда
преобразуется при помощи умножения на веса:

wt = 0.5*{1-cos[*(t — 0.5)/m]}     (для t=0 до m-1)
wt = 0.5*{1-cos[*(N — t +
0.5)/m]}     (для t=N-m до N-1)

где m выбирается так, чтобы 2*m/N было равно коэффициенту
пропорциональности сглаживаемых данных (p).

Окна данных и оценки
спектральной плотности

На практике, при анализе данных обычно не очень
важно точно определить частоты основных функций
синусов или косинусов. Скорее, т.к. значения
периодограммы — объект существенного случайного
колебания, можно столкнуться с проблемой многих
хаотических пиков периодограммы. В этом случае
хотелось бы найти частоты с большими спектральными
плотностями
, т.е. частотные области, состоящие
из многих близких частот, которые вносят
наибольший вклад в периодическое поведение
всего ряда. Это может быть достигнуто путем
сглаживания значений периодограммы с помощью
преобразования взвешенного скользящего
среднего. Предположим, ширина окна скользящего
среднего равна m (должно быть нечетным
числом); следующие наиболее часто используемые
преобразования (заметим: p = (m-1)/2).

Окно Даниэля (равные веса). Окно
Даниэля (Daniell, 1946) означает простое (с равными
весами) сглаживание скользящим средним значений
периодограммы; т.е. каждая оценка спектральной
плотности вычисляется как среднее m/2
предыдущих и последующих значений
периодограммы.

Окно Тьюки. В окне Тьюки (Blackman and Tukey, 1958)
или Тьюки-Ханна (Hanning) (названное в честь Julius Von Hann),
для каждой частоты веса для взвешенного
скользящего среднего значений периодограммы
вычисляются как:

wj = 0.5 + 0.5*cos(*j/p)    (для j=0 до p)
w-j = wj    (для j 0)

Окно Хемминга. В окне Хемминга
(названного в честь R. W. Hamming) или Тьюки-Хемминга
(Blackman and Tukey, 1958), для каждой частоты, веса для
взвешенного скользящего среднего значений
периодограммы вычисляются как:

wj = 0.54 + 0.46*cos(*j/p)    (для j=0 до p)
w-j = wj    (для j 0)

Окно Парзена. В окне Парзена (Parzen, 1961),
для каждой частоты, веса для взвешенного
скользящего среднего значений периодограммы
вычисляются как:

wj = 1-6*(j/p)2 + 6*(j/p)3    (для
j = 0 до p/2)
wj = 2*(1-j/p)3    (для j = p/2 + 1 до p)
w-j = wj    (для j 0)

Окно Бартлетта. В окне Бартлетта (Bartlett,
1950) веса вычисляются как:

wj = 1-(j/p)    (для j = 0 до p)
w-j = wj    (для j 0)

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

Подготовка данных к анализу

Теперь рассмотрим несколько других
практических моментов спектрального анализа.
Обычно, полезно вычесть среднее из значений ряда
и удалить тренд (чтобы добиться стационарности)
перед анализом. Иначе периодограмма и
спектральная плотность «забьются» очень
большим значением первого коэффициента при
косинусе (с частотой 0.0). По существу, среднее — это
цикл частоты 0 (нуль) в единицу времени; т.е.
константа. Аналогично, тренд также не
представляет интереса, когда нужно выделить
периодичность в ряде. Фактически оба этих
эффекта могут заслонить более интересные
периодичности в данных, поэтому и среднее, и
(линейный) тренд следует удалить из ряда перед
анализом. Иногда также полезно сгладить данные
перед анализом, чтобы убрать случайный шум,
который может засорять существенные
периодические циклы в периодограмме.

Результаты для случая, когда в
ряде отсутствует периодичность

В заключение, зададим вопрос: что, если
повторяющихся циклов в данных нет, т.е. если
каждое наблюдение совершенно независимо от всех
других наблюдений? Если распределение
наблюдений соответствует нормальному, такой
временной ряд может быть белым шумом (подобный
белый шум можно услышать, настраивая радио). Если
исходный ряд — белый шум, то значения
периодограммы будут иметь экспоненциальное
распределение. Таким образом, проверкой на
экспоненциальность значений периодограммы
можно узнать, отличается ли исходный ряд от
белого шума. Пользователь может также построить
одновыборочную статистику d статистику
Колмогорова-Смирнова (cм. также раздел Непараметрическая статистика и
распределения
).

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


Быстрое преобразование Фурье (БПФ)

  • Общее введение
  • Вычисление БПФ во временных
    рядах

За дальнейшей информацией обратитесь к Анализу временных рядов и
следующим разделам:

  • Идентификация модели
    временных рядов
  • АРПСС (Бокс и Дженкинс) и
    автокорреляции Вводный обзор АРПСС
  • Прерванные временные ряды
  • Экспоненциальное
    сглаживание
  • Сезонная декомпозиция (метод
    Census I)
  • Сезонная корректировка X-11
    (метод Census II)
  • Таблицы результатов
    корректировки X-11
  • Анализ распределенных лагов
  • Одномерный анализ Фурье
  • Кросс-спектральный анализ
  • Основные понятия и принципы

Общее введение

Интерпретация результатов спектрального
анализа обсуждается в разделе Основные
понятия и принципы
, однако там мы не
обсуждали вычислительные проблемы, которые в
действительности очень важны. До середины 1960-х
для представления спектрального разложения
использовались точные формулы для нахождения
параметров синусов и косинусов. Соответствующие
вычисления требовали как минимум N**2 (комплексных)
умножений. Таким образом, даже сегодня
высокоскоростному компьютеру потребовалось бы
очень много времени для анализа даже небольшого
временного ряда (для 8,000 наблюдений
потребовалось бы по меньшей мере 64 миллиона
умножений).

Ситуация кардинально изменилась с открытием
так называемого алгоритма
быстрого преобразования Фурье, или БПФ для
краткости. Достаточно сказать, что при
применении алгоритма БПФ время выполнения
спектрального анализа ряда длины N стало
пропорционально N*log2(N) что конечно
является огромным прогрессом.

Однако недостаток стандартного алгоритма БПФ
состоит в том, что число данных ряда должно быть
равным степени 2 (т.е. 16, 64, 128, 256, …). Обычно это
приводит к необходимости добавлять нули во
временной ряд, который, как описано выше, в
большинстве случаев не меняет характерные пики
периодограммы или оценки спектральной
плотности. Тем не менее, в некоторых случаях,
когда единица времени значительна, добавление
констант во временной ряд может сделать
результаты более громоздкими.

Вычисление БПФ во временных
рядах

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

Как упоминалось ранее, для применения
стандартного (и наиболее эффективного) алгоритма
БПФ требуется, чтобы длина исходного ряда была
равна степени 2. Если это не так, должны быть
проведены дополнительные вычисления. Будут
использоваться простые точные вычислительные
формулы, пока исходный ряд относительно мал, и
вычисления можно выполнить за относительно
короткое время. Для длинных временных рядов,
чтобы применить алгоритм БПФ, используется
основной подход, описанный Monro и Branch (1976). Этот
метод требует значительно больше памяти; однако
ряд рассматриваемой длины может анализироваться
все еще очень быстро, даже если число наблюдений
не является степенью 2.

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

Дополнительная информация по методам анализа данных, добычи данных,
визуализации и прогнозированию содержится на
Портале StatSoft (http://www.statsoft.ru/home/portal/default.asp)
и в Углубленном Учебнике StatSoft (Учебник с формулами).


Все права на материалы электронного учебника принадлежат компании StatSoft


Часть 1 – Энтропия

Часть 2 – Mutual Information

В этой 3-й части мы поговорим про Machine Learning, а именно, про задачу прогноза, в контексте теории информации.

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

Что такое Machine Learning?

Чтобы двигаться дальше, нам нужны базовые понятия из ML, такие как: факторы, target, Loss-функция, обучающий и тестовый пулы, overfitting и underfitting и их варианты, регуляризация, разные виды утечки данных.

Факторы (features) – это то, что вам даётся на вход, а сигнал (target) – это то, что нужно спрогнозировать, используя факторы. Например, вам нужно спрогнозировать температуру завтра в 12:00 в определённом городе – это target, а дано вам множество чисел про сегодня и предыдущие дни: температура, давление, влажность, направление и сила ветра в этом и соседних городах в разное время дня – это факторы.

Обучающие Данные – это множество примеров (aka сэмплов) с известными правильными ответами, то есть строчки в таблице, в которой есть и поля с факторами features=(f1, f2, …, fn), и поле target. Данные принято делить на две части – обучающий пул (train set) и тестовый пул (test set). Выглядит это примерно так.

Обучающий пул:

id

f1

f2

target

predict

1

1.234

3.678

1.23

?

2

2.345

6.123

2.34

?

18987

1.432

3.444

….

5.67

?

Тестовый пул:

id

f1

f2

target

predict

18988

6.321

6.545

4.987

?

18989

4.123

2.348

3.765

?

….

30756

2.678

3.187

2.593

?

В общих чертах задачу прогноза можно сформулировать в духе соревнования на kaggle.com:

Задача прогноза (задача ML). Вам дан обучающий пул. Реализуйте в виде кода функцию mathcal{P}(f_1, ldots, f_k), которая по заданным факторам возвращает значение, максимально близкое к target. Мера близости задана с помощью Loss-функции, и значение этой функции называется ошибкой прогноза: error=Loss(predict, target), где predict=mathcal{P}(f_1,ldots,f_k). Качество прогноза определяется средним значением ошибки во время применения этого прогноза в жизни, но по факту для оценки этой средней ошибки будет использоваться некоторый тестовый пул, который от вас скрыт.

У величины xi = predict - target есть специальное название – невязка. Два популярных варианта Loss-функции для задач прогноза вещественной величины – это:

  • средний квадрат невязки, что также называют среднеквадратической ошибкой или mse (mean square error);

  • средний модуль невязки, что также называют L1-ошибкой.

На практике задача ML более общая и высокоуровневая. А именно, нужно разработать в меру универсальную ML-модель – способ получения функций mathcal{P}(f_1, ldots, f_n) по данному обучающему пулу и заданной Loss-функции. А также нужно уметь делать model evaluation: мониторить качество прогноза в работающей системе, уметь обновлять обученные модели (пошагово, делая новые версии, или модифицируя внутренние веса модели в режиме online), улучшать качество модели, контролировать чистоту и качество факторов.

Процесс получения функции mathcal{P}(f_1, ldots, f_n)по обучающему пулу называется обучением. Популярны такие варианты ML-моделей (классы моделей):

  • Линейная модель, в которой mathcal{P}(f_1,ldots,f_n)=sum_{i} w_icdot f_i, и процесс обучения сводится к подбору параметров (w_1,w_2,ldots, w_n), что делается обычно методом градиентного спуска.

  • Gradient boosted trees (GBT) – модель, в которой функция выглядит как сумма нескольких слагаемых (сотен, тысяч), где каждое слагаемое – это решающее дерево, в узлах которого находятся условия на факторы, а в листьях – конкретные числа; каждое слагаемое можно представлять как какую-то систему вложенных if-конструкций с условиями на значения факторов, где в конечных узлах находятся просто числа. GBT – это не только про то, что решение есть сумма деревьев, но и вполне конкретный алгоритм получения этих слагаемых. Для обучения GBT есть много готовых программ: CatBoost, xgboost. См. статью на хабре.

  • Neural Networks – в простейшем базовом варианте это модель вида mathcal{P}(vec{f}) = W_kodot W_{k-1}odot ldots W_1 odot vec{f}, где W_j– матрицы, размеры которых определяет разработчик модели, факторы представлены в виде вектора vec{f} = |f_1,ldots,f_n|^{T}, оператор odot соответствует умножению вектора на матрицу с последующим обнулением всех отрицательных чисел в результирующем векторе. Операторы выполняются в порядке справа налево – это важно в случае наличия такого обнуления. Матрицы называются слоями нейросети, а количество матриц – глубиной нейросети. Вместо обнуления отрицательных можно брать произвольные другие нелинейные преобразования. Если бы не было нелинейных преобразований после умножения вектора на матрицу, все матрицы можно было бы «схлопнуть» в одну и пространство возможных функций не отличалось бы того, что задаётся линейной моделью. Я описал линейную архитектуру нейросети, но возможны и более сложные, например,

     mathcal{P}(vec{f})=W_6odot(W_1odot W_2 cdot vec{f}  + W_3odot W_4 odot W_5odot W_2odot vec{f})+W_7odot vec{f}.

    Кроме самых разных поэлементных нелинейных преобразования и умножения на матрицу в нейросетях могут использоваться операторы скалярного произведения векторов и оператор поэлементного максимума для двух векторов одной размерности, объединения векторов в один более длинный и др. Можно думать об общей архитектуре для функции predict, когда на вход поступает вектор, а дальше с помощью операторов {+,-,cdot, max }, нелинейных функций {tanh, mathrm{abs}, max(0, cdot), ldots}и весов {w_1, w_2, ldots} конструируется ответ. В этом смысле нейросети могут представлять самой функции довольно общего вида. По факту, архитектуру функции mathcal{P} называют нейросетью тогда, когда в ней присутствует что-то похожее на цепочку вида W_kodot W_{k-1}odot ldots W_1 odot vec{f}. А так, мы по сути имеем задачу регрессии – подобрать параметры (веса) в параметрически заданной функции, чтобы минимизировать ошибку. Есть множество методов обучения нейросетей, они в своем большинстве итерационные, и за шаг изменения весов отвечает модуль, который программисты называют оптимизатором, см. например AdamOptimizer.

Терминология ML

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

Overfitting

Overfitting (переобучение) – это когда выбранная вами модель сложнее, чем реальность, стоящая за target И/ИЛИ данных слишком мало, чтобы позволить себе обучать такую сложную модель. Есть две основные причины overfitting:

  • too complex model: Устройство модели сильно сложнее, чем реальность или в своей сложности ей не соответствует. Проще это визуализировать на примере однофакторной модели.
    Пусть в обучающем пуле 11 точек (x,y)= (f_i, target_i), а модель – суть многочлен 10-й степени. Можно подобрать такие коэффициенты в многочлене, что он «заглянет» в каждую точку обучающего пула, но при этом не будет обладать хорошим качеством прогноза:

    Голубая линия – многочлен 10 степени, 
который смог в точности воспроизвести обучающий пул из 11 точек. 
Но правильная модель скорее линейная (чёрная линия), а отклонения от неё – это или шум, или что-то, объясняемое факторами, которых у нас нет

    Голубая линия – многочлен 10 степени,
    который смог в точности воспроизвести обучающий пул из 11 точек.
    Но правильная модель скорее линейная (чёрная линия), а отклонения от неё – это или шум, или что-то, объясняемое факторами, которых у нас нет

    Многочлен 10-й степени – это очевидный и часто используемый пример переобучения. Многочлены невысокой степени, но от многих переменных тоже могут давать overfitting. Например, вы можете подбирать модель predict = многочлен степени 3 от 100 факторов (кстати, сколько в нём коэффициентов?), в то время как реальность соответствует predict = многочлен степени 2 от факторов + случайный шум. При этом, если у вас данных довольно много, то классические методы регрессии могут дать приемлемый результат, и коэффициенты при членах степени 3 будут очень маленькими. Эти члены степени 3 будут вносить маленький вклад в ответ при прогнозе для типичных значений факторов (как в случае с интерполяцией). Но в зоне крайних значений факторов, где прогноз уже скорее экстраполяция, а не интерполяция, высокие степени будут вносить заметный вклад и портить качество прогноза. Ну и, конечно, когда данных мало, шум может начать восприниматься за чистую монету, и прогноз будет стараться заглянуть в каждую точку обучающего пула.

  • insufficient train data: Модель может более менее соответствовать реальности, но при этом вы можете не иметь достаточное количество обучающих данных. Приведём опять пример с многочленом: пусть и реальность и ваша модель есть многочлен 3-й степени от двух факторов. Этот многочлен задаётся 10 коэффициентами, то есть пространство возможных прогнозаторов у вас 10-мерное. Если в обучающем пуле у вас 9 примеров, то, записав равенство «многочлен(featuresi) = targeti» для каждого примера, вы получите 9 уравнений на эти коэффициенты, и их недостаточно, чтобы однозначно определить 10 коэффициентов. В пространстве возможных прогнозаторов вы получите прямую и каждая точка на этой прямой есть прогнозатор, который идеально повторяет то, что у вас в обучающем пуле. Вы можете случайно выбирать один из них, и это с высокой вероятностью будет плохой прогнозатор.

    • Важное замечание про обучение многопараметрических моделей. Пример выше может показаться искусственным, но правда в том, что современные нейросетевые модели могут содержать миллионы параметров или даже больше. Например, языковая модель GPT3 содержит 175 млрд. параметров.

      Если в вашей модели N = 175 млрд. параметров, а размер обучающего пула M = 1 млрд, то в общем случае вы имеете бесконечное множество моделей идеально подходящих под обучающий лог, и это множество – суть многообразие размерности NM = 174 млрд.

      В случае с глубокими многопараметрическими нейросетями эффект «insufficient train data» проявляется в полный рост, если действовать неправильно. А именно, если искать строгий минимум ошибки на обучающем пуле, то получается плохой, overfitted, прогнозатор. Поэтому на практике так не делают. Выработаны интересные техники и интуиция про то, как обучать модель с N параметрами на M примерах, где M заметно меньше N. Используются регуляризационные добавки к Loss-функции (L2, L1), Stochastic Gradient Descent, Dropout Layers, Early Stopping, prunning и другие техники.
      ВАЖНО: Знание этих техник и умение их применять во многом и определяет экспертизу в ML.

Underfitting

Underfitting – это когда модель прогноза проще, чем реальность. Также как и в случае с overfitting здесь есть две основных причины:

  • too simplistic model: выбрана слишком простая модель, сильно упрощающая то, что в реальности стоит за target. Такие случаи по-прежнему встречаются в продакшене, и обычно это линейные модели. Линейные модели привлекают своей простотой и наличием математических теорем, обосновывающих и описывающих их предиктивные способности. Но можно с уверенностью сказать, что если вы решите с помощью линейных моделей прогнозировать курс валют, погоду, вероятность покупки или вероятность возврата кредита, то вы получите слабую модель – under-fitted weak model.

  • too early stopping: обычно процессы обучения моделей итеративны, и если проделать слишком мало итераций, то получается недообученная модель;

Data leakage

Data leakage – это когда вы во время обучения использовали то, что есть в тестовом пуле или данные, которые недоступны или отличаются от тех, что будут доступны на практике при применении модели в жизни. Модели, которые сумеют воспользоваться этой утечкой, будут иметь неоправданно низкую ошибку на тестовом пуле и выбраны для использования в жизни. Есть несколько видов data leakage:

  • утечка таргета в факторы: Если, в задаче с прогнозом температуры вы добавите фактор – средняя температура за K дней, и в эти K дней при обучении случайно попадёт в том числе и целевой день, на который нужно сделать прогноз, то вы получите испорченный или нерабочий прогнозатор. Могут быть более сложные утечки через факторы. В случае прогноза какого-то будущего события вы должны многократно убедиться, что при вычислении факторов для обучающего пула не используются данные, которые получены во время этого события или после.

  • простая утечка тестового пула: Иногда всё просто – в обучающий пул попадают строчки из тестового пула.

  • утечка через подбор гиперпараметров: Правило, что для обучения разрешено использовать только обучающий пул, обманчиво просто. Здесь можно легко обмануться. Типичный пример неявной утечки – это подбор гиперапараметров при обучении. У всякого алгоритма обучения есть гиперпараметры. Например, у градиентного спуска есть параметр «скорость обучения» и параметр, задающий критерий останова. Ещё вы к Loss-функции можете добавить регуляризационные компоненты R_2=lambda_2 sum w_i^2и R_1=lambda_1 sum |w_i|, и вот вам ещё два гиперапараметраlambda_1и lambda_2. Приставка «гипер» используется, чтобы отличать внутренние параметры модели (веса), от параметров влияющих на процесс обучения.
    Итак, предположим, вы решили в цикле проверять разные гиперпараметры обучения и выбрать те гиперпараметры, на которых достигается минимум ошибки на тестовом пуле. Надо понимать, что здесь происходит утечка тестового пула в модель. «Засаду», которая здесь скрывается, легко понять, именно с точки зрения теории информации. Когда вы задаете вопрос, чему равна ошибка на тестовом пуле, вы получаете информацию о тестовом пуле. А когда вы, на основании этой информации, принимаете какие-то решение, влияющие в конечном итоге на веса вашей модели, то значит биты этой информации оказываются в весах модели.
    Ещё можно проиллюстрировать концепцию утечки тестового пула доведя до абсурда понятие гиперапараметров. Ведь разница между весами (параметрами) и гиперпараметрами (параметрами процесса обучения) чисто формальная. Давайте веса модели, например коэффициенты многочлена 10 степени все назовоём гиперапараметрами, а процесса обучения у нас не будет (множество внутренних весов у нас будет пустым). И давайте запустим безумный цикл по перебору всех возможных комбинаций значений гиперпараметров (каждый параметр будем перебирать от -1000 до +1000 с шагом 10^{-6}) и вычислению ошибки на тестовом пуле. Понятно, что через условные миллиард лет мы найдём такую комбинацию параметров, которая даёт очень маленькую ошибку на тестовом пуле и картинку выше мы сможем запостить снова, но уже с другой подписью:

    Голубая линия – многочлен 10 степени, 
который смог в точности воспроизвести тестовый пул из 11 точек.

    Голубая линия – многочлен 10 степени,
    который смог в точности воспроизвести тестовый пул из 11 точек.

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

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

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

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

  • утечка данных из будущего: например, когда вы обучаете прогнозаторы вероятности клика на рекламное объявление, то в идеале обучающий пул должен содержать данные, доступные исключительно на даты меньше даты X, где X – минимальная дата событий в тестовом пуле. Иначе возможны хитрые утечки типа следующей. Пусть мы взяли множества событий показы рекламы размеченной target = 1 для кликнутых событий показа, и target = 0 для некликнутой. И пусть мы разбили множество событий на тестовый и обучающий пул не по границе даты X, а случайным образом, например, в пропорции 50:50. Известно, что пользователи любят кликать сразу на несколько объявлений одной тематики подряд в течение нескольких минут, отбирая нужный товар или услугу. И тогда такие серии нескольких кликов будут случайно делиться на две части – одна пойдёт в обучающий, другая – в тестовый пул. Можно искусственным образом сделать модель с утечкой: возьмём самую хорошую правильную модель без утечки и дополнительно запомним факты «пользователь U кликал на тематику T» из обучающего лога. Затем при использовании модели будем совсем ненамного повышать вероятность клика для случаев (U, T) из этого множества. Это улучшит модель с точки зрения ошибки на тестовом пуле, но ухудшит её с точки применения на практике на новых данных. Мы описали искусственную модель, но несложно представить естественный механизм такого запоминания пар (U, T) и завышения вероятности для них в нейросетях и других популярных алгоритмах ML. Наличие таких утечек из будущего усугубляет проблему утечки через гиперпараметры. На таком тестовом пуле с утечкой будут выигрывать модели перекрученные в сторону запоминания данных. Более общим образом эту проблему можно описать так: модель может переобучиться под test set, если MI между примером в train set и примером в test set больше, чем между примером в train set и реальным примером при использовании модели в жизни.

Задачи

Задача 3.1. Пусть реальность такова, что target = w_0 + w_1cdot f + w_2cdot f^2+ nu, где nu– случайное число из mathcal{N}(0,alpha^2), alpha=0.2. Факторf сэмплируется из mathcal{N}(0,1). Сколько в среднем нужно обучающих данных, чтобы, имея соответствующую реальности модель, получить оценки весов равные истинным с точностью до среднеквадратичной ошибки varepsilon=0.01? Считайте, что априорное распределение весов – это нормальное распределение mathcal{N}(0,sigma^2). Запишите ответ как функцию от alpha, varepsilon, sigma. Как ошибка прогноза зависит от размера обучающего лога?

Ответ

Есть такой способ получения ответа «от физиков». Пусть размер обучающего лога равен n. Из задачи 1.16 мы имеем интуицию, что varepsilon = c / sqrt{n}. Константа c в этой задаче является функцией от alpha и sigma. Но есть ли зависимость c от sigma? На самом деле нет. Действительно, шум alphacdot nu как бы стирает цифры в десятичной записи значения target, начиная с какого-то места после десятичной точки. При увеличении sigma остается больше значащих цифр и ровно настолько больше значащих цифр мы узнаем про коэффициенты {w_0, w_1, w_2}. И тем самым при увеличении sigma позиция в десятичной записи {w_0, w_1, w_2}, начиная с которой есть неопределённость остается той же. Из того, что sigma не участвует в формуле и из соображений размерности, формула должна выглядеть так

varepsilon =c_0 cdot {alpha over sqrt{n}}

Где константа c_0 вычисляется эмпирически равна примерно 0.9. Соответственно получаем ответ n = 0.8cdot left({alpha over varepsilon} right)^2.

Ошибка самого значения targetравна

  mse=sqrt{varepsilon^2 + alpha^2}=alpha cdot sqrt{1+0.8/n}

Это следствие того, что при сложении двух независимых случайных величин w'_0 + w'_1cdot f + w'_2cdot f^2 и  nu дисперсии складываются. Они должны быть независимы при правильном прогнозе, потому то если они зависимы, то модель можно улучшить, сделав некую калибровку прогноза.

Таким образом, в этой задаче прогноза и многих практических задачах есть неустранимая ошибка alpha. Она определяет две вещи:
(1) то, какую ошибку вы получите для случая очень большого обучающего пула (при nto inftyимеем mse=alpha sqrt{1  + beta /n}to alpha)
(2) а также размер пула, необходимый для достижения какой-либо заданной ошибки внутренних параметров, квадратично зависит от этой неустранимой ошибки: n=beta cdot (alpha/varepsilon)^2

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

n={beta over (mse/alpha)^2 -1 }approx  beta cdot (alpha / mse)^2 ; ;mathrm{при;маленьких} ; alpha/mse

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

К тому же факторы часто зависимы и сами по себе содержат в себе шумную компоненту. Контроль пользы от новых факторов и отбор факторов (feature selection) – это отдельная сложная тема для разговора.

Задача 3.2. Пусть реальность такова, что target есть многочлен 2-й степени от 10 факторов, при этом коэффициенты в многочлене – это числа разово сэмплированные из mathcal{N}(0,1), а факторы – независимые случайные числа из mathcal{N}(0,1). Сколько нужно данных в обучающем логе, чтобы получить приемлемое качество прогноза, когда все факторы лежат на отрезке [-1, 1], для случая, когда модель верна, то есть является многочленом 2-й степени?
А если модель есть многочлен 3-й степени?

Задача 3.3. Постройте модель пользователей, которые в каждый момент времени склонны больше кликать на объявления какой-то одной тематики, заинтересовавшей их в данный момент, и убедитесь, что деление пула на обучающий и тестовый пул не по времени, а случайно, приводит к утечке данных из будущего.
Можно взять, например, такую модель. У каждого пользователя есть 10 любимых тематик, мы их знаем и храним в профиле пользователя. Пользовательская активность разделена на сессии, каждая сессия длится 5 минут, и во время сессии каждый пользователь видит ровно 10 объявлений по одной из его 10 тематик, тематика из любимых выбирается случайно. В каждую сессию пользователь интересуется одной из этих 10 тематик особенно сильно, и нет никаких факторов, позволяющих угадать, какой именно. Вероятность клика для такой «горячей» тематики в 2 раза больше, чем обычно. Для каждого пользователя по истории его кликов мы хорошо знаем вероятности клика для его 10 тематик, и пусть, если их упорядочить по интересности для конкретного пользователя, вектор этих вероятностей равен {0.10, 0.11, 0.12, …, 0.19} .
Идентификатор пользователя и номер категории и есть доступные нам факторы. Насколько можно увеличить правдоподобие прогноза вероятности на тестовом пуле, если предположить, что каждая сессия разделилась ровно пополам между тестовым и обучающим пулами (5 событий пошли в один пул и 5 – во второй), и статистику по 5 событиям из обучающего пула можно использовать для прогноза вероятности клика на 5 других?
Предлагается экспериментально получить overfitting из-за утечки данных из будущего , используя какую-либо программу для ML, например, CatBoost.

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

Задача 3.4. Какая средняя varepsilon-ИС функции вида

p(f_1,ldots,f_n)=sum_{i=1}^n w_icdot f_i

от n факторов, в которой веса сэмплированы из распределения mathcal{N}(0,sigma^2), а факторы сэмплируются из mathcal{N}(0,1).

Ответ

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

delta p^2 = sum_i delta w_i ^2 f_i^2

В среднем по всем возможным значениям факторов получаем delta p^2 = {sum}_i delta w_i ^2. Если мы хотим delta p^2 le varepsilon^2, то веса нужно передавать с точностью delta w_i le  varepsilon^2 /n.

Таким образом, согласно ответу в задаче 1.9, в среднем нужно будет передать ncdot log(sqrt{n}cdot sigma/varepsilon)бит. Здесь видно важное свойство очень сложных моделей – для определения сложности сложной модели не так важно какое и насколько хорошее у вас априорное знание о весах, и не так важно, какая вам требуется точность или какая точность по факту достигается. Основная компонента в ИС модели это 0.5 cdot ncdot log n, где n – некое эффективное количество весов, равноправно и с пользой участвующих в модели.

Задача 3.5. Функция t(id) задаётся таблицей как функция от одного фактора – id категории. Есть N = 1 млн. категорий, их доли в данных различны, и можно считать, что их доли получены сэмплированием 1 млн. чисел из экспоненциального распределения с последующей нормализацией, чтобы сумма долей была 1 (такое сэмплирование соответствует разовому сэмплированию из распределения Дирихле с параметрами {1,1,…,1} – 1 млн. единичек). Значения функции t для этих категорий сэмплированы из бета-распределенияB(alpha, beta), alpha=2,;beta=10 независимо от их доли.
Какова средняяvarepsilon-ИС таких функции?

Ответ

Во-первых, давайте опишем поведение этой функции в районе очень маленьких varepsilon. Когда varepsilon ll 1/10^6, мы вынуждены хранить табличную функцию – 1 млн. значений. Можно численно или по формуле из википедии почитать энтропию бета-распределения B(alpha, beta). Далее можно положить, что мы хотим передать другому программисту, что значение t для какой-то конкретной категории i лежит в районе некоторого фиксированного значения t_iс разрешённой погрешностью varepsilon. Это все равно что сузить «колпак» распределения B(alpha, beta) до, скажем, нормального распределения mathcal{N}(t_i, varepsilon^2). Разница энтропий B(alpha, beta) и mathcal{N}(t_i, varepsilon^2) равна примерно -2.38 -log varepsilon. Это надо сделать для 1 миллиона категорий, поэтому для очень маленьких varepsilon ответ Ncdot (-2.38 + log 1/varepsilon).

Для больших varepsilon (но всё ещё заметно меньше 1) работает другая стратегия. Давайте разобьем промежуток [0,1] на одинаковые отрезки длины delta, и будем для каждой категории передавать номер отрезка, в который попало значение функции для этой категория. Эта информация имеет объём Ncdot (H_{B(2,10)}+ log(1/delta)) (см. задачу 1.7). В результате в каждом ответе наша ошибка будет ограничена длиной отрезка, в который мы попали. Если ответом считать середину отрезка и предполагать равномерность распределения (что допустимо, когда отрезки маленькие), то ошибка в рамках отрезка длины deltaравна delta/sqrt{12} . Таким образом, чтобы получить заданную среднюю ошибку varepsilon, нужно брать отрезки длиной delta=sqrt{12}cdot varepsilon. Итого ответ равен Ncdot (-2.20487 + log 1/varepsilon).

Оба подхода дают ответ вида Ncdot (c + log 1/varepsilon). Идея использования фильтра Блюма для хранения множеств категорий для каждого отрезка даёт тот же результат.

Задача 3.6. Пусть в предыдущей задаче смысл функции – вероятность клика на рекламное объявление. Перенумеруем все категории в порядке убывания их истинной доли и обозначим порядковый номер как order_id, а исходный случайный идентификатор как id. Предлагается рассмотреть разные варианты ML с огрублением номера (идентификатора) категории до натурального числа от 1 до 1000. Это ограничение может быть связано, например, с тем, что вы хотите хранить статистику по историческим данным, и у вас есть возможность хранить только 1000 пар (показы, клики). Варианты огрубления идентификатора категории:
(a) огрубленный идентификатор равный id’ = id // 1000 (целочисленное деление на 1000);
(b) огрубленный идентификатор равный id’ = order_id для топовых 999 категорий по доле, а всем остальным назначить идентификатор 1000;
(с) какой-то другой способ огрубления id’ = G(id), где G – детерминированная функция, которую вы можете сконструировать на базе статистики кликов на 100 млн показов, где каждый показ имеет какую-то категорию с вероятностью ровной доле категорий, а клик происходит согласно вероятности клика для этой категории.
Оцените значение ошибки error в этих трёх вариантах, где error = (LL_0 - LL)/n, а LL=-{sum}_{i=1}^n t_i log(p_i) + (1-t_i)log(1-p_i), а LL_0– это значение LL для идеального прогноза p_i = t_i.

Задача 3.7. Пусть реальность такова, что target ={sum}_{i=1}^{45} w_icdot f_i +0.02 cdot nu, где nu– это шум, случайное число из mathcal{N}(0,1). А ваша модельpredict = {sum}_{i=1}^{50} w_icdot f_i +0.02 cdot nu,
где веса w_1,w_2, ldots,w_{70} вам неизвестны и имеют априорное распределение лапласа со средним 0 и дисперсией 1. Последние 5 факторов модели по факту бесполезны для прогноза. Какого размера должен быть обучающий пул, чтобы получить среднеквадратическую ошибку прогноза

mse={sum}_{i=1}^n(predict_i - target_i)^2 / n  le 0.04^2.

Сгенерируйте три пула №1, №2, №3 c размерами 50, 50, 100000 и узнайте как лучше действовать – (а), (б) или (в):
(а) соединить два пула в один и на основании объединённого пула найти наилучшие веса, минимизирующие mseна этом пуле
(б) для различных пар (lambda_1, lambda_2), найти наилучшие веса, минимизирующие mse  + lambda_1  sum |w_i| + lambda_2 cdot sum w_i^2 на пуле №1; из всех пар выбрать такую пару lambda_1, lambda_2 , на которых достигается минимум mseна пуле №2; пары (lambda_1, lambda_2) предлагается взять из множества A times B, где A и B геометрические прогрессии с q=1.2 от1/10^3до 1. Нарисуйте график ошибки как функцию от lambda_1, положивlambda_2 = 0.
(в) соединить два пула в один и на основании объединённого пула найти наилучшие веса, минимизирующие mse  + lambda_1  sum |w_i| + lambda_2 cdot sum w_i^2на этом пуле; значения lambda_1, lambda_2 взять из пункт (б).
Метод тем лучше, чем меньше ошибка на пуле №3, который соответствует применению модели в реальности. Стабилен ли победитель, когда вы генерируете новые пулы? Нарисуйте таблицу 3×3 mse ошибок на трёх пулах в этих трёх методах.
Как падает ошибка на пуле №3 с ростом размеров пулов №1 и №2?

Задача 3.8. Пусть реальность такова, что target =l_1(vec{f})+ l_2(vec{f}) +l_3(vec{f})+nu, где

Какая будет mse ошибка прогноза в среднем на достаточно больших тестовых и обучающих пулах по различным вариантам функции target(vec{f})(то есть для различных вариантов сэмплирования коэффициентов w_i, ; w_{ij},; w_{ijk}) для случая, когда ваша модель есть

  • (а0) mathcal{P}(vec{f})=l_1(vec{f}) с истинными значениями коэффициентов.

  • (б0) mathcal{P}(vec{f})=l_1(vec{f})+l_2(vec{f}) с истинными значениями коэффициентов.

  • (в0) mathcal{P}(vec{f})=l_1(vec{f})+l_2(vec{f})+l_3(vec{f}) с истинными значениями коэффициентов.

  • (а1) mathcal{P}(vec{f})=l_1(vec{f}) с обучением (то есть с коэффициентами, полученными регрессией).

  • (б1) mathcal{P}(vec{f})=l_1(vec{f})+l_2(vec{f}) с обучением, без знания, какие коэффициенты равны 0.

  • (в1) mathcal{P}(vec{f})=l_1(vec{f})+l_2(vec{f})+l_3(vec{f}) с обучением, без знания, какие коэффициенты равны 0.

  • (в2) mathcal{P}(vec{f})=l_1(vec{f})+l_2(vec{f})+l_3(vec{f}) с обучением и знанием, какие именно коэффициенты равны 0.

  • (г1) mathcal{P}(vec{f}) нейросеть глубины d = 5, размер внутренних слоёв h=5 (внутренние матрицы имеют размер h x h).

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

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

Пронаблюдайте явление переобучения в этих моделях (кроме первых двух), уменьшив размер обучающего пула до некоторого значения. Для каждой модели это будет свой размер. Попробуйте избавится от явления переобучения
1) добавляя к Loss-функции регуляризационные члены R_1=lambda_1cdot sum |w|, R_2=lambda_2 cdot sum w^2;
2) используя early stopping;
3) используя SGD;
4) меняя архитектуру нейросети;
5) добавляя Dropout;
6) комбинацию перечисленных методов.

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

Задача 3.9. Пусть реальность такова, что

target =w_0+sum_{i=1}^{m} w_icdot f_i+nu,

где

Как зависит ошибка определения весов w_i от размера обучающего лога n и m, ; delta?

Ответ

В первом приближении ответ равен

mse(w_i) = {delta over sqrt{n}}.;;;;;(1)

И это хорошее приближение для случая n gg m (n заметно больше m).

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

Средний квадрат ошибки весов, нормированный на δ^2 для m=50.
Ось X нормирована на m, то есть X=1  соответствует n = m.
Красная линия соответствует функции 1/n.

Средний квадрат ошибки весов, нормированный на δ^2 для m=50.
Ось X нормирована на m, то есть X=1 соответствует n = m.
Красная линия соответствует функции 1/n.

Если у кого-то есть хороший ответ про аналитическое приближение в окрестности n ~ m, пишите в личку.

Задача 3.10. Пусть реальность такова, что

target =sum_{i=1}^{11} w_icdot g_i+nu,

где

Как будет выглядеть наилучший прогноз, если вам известны точные значения весовw_i?
Как лучше всего реализовать обучение в этой задаче (когда веса неизвестны и их надо «обучить»)?

Ответ

Рассмотрим частную задачу про добавление нового фактора f_1 к имеющемуся прогнозу F_0. Пусть

target =  F_0 + w_1cdot g_1 + nu_1,; \ nu_1 sim mathcal{N}(0, delta_1^2), ; g_1sim mathcal{N}(0,1), \ predict_0 = F_0, \ f_1 = g_1 + zeta_1, ; ;  zeta_1sim  mathcal{N}(0,z_1^2)

Можно ли уменьшить ошибку mse, взяв новый прогноз

predict_1 = predict_0 + w_1cdot f_1?

Текущая ошибка равна дисперсии D(target - predict_0) =  w_1^2 + delta_1^2. Ошибка нового прогноза будет равна D(target - predict_1) = w_1^2 cdot z_1^2+delta_1^2, она меньше старой ошибки, если z_1 < 1.

Давайте построим последовательность улучшений прогноза

predict_0 = 0,\ predict_1 = predict_0 + t_1 cdot w_1 cdot f_1, \ predict_2=predict_1 + t_2cdot w_2cdot f_2, \ ldots \ predict_{m} = predict_{m-1} + t_{m}cdot w_{m}cdot f_{m}

Числа t_i равны 0 или 1 – это индикаторы того, берём мы фактор f_iв линейный прогноз или нет. D(target - predict_0) =  sigma_0^2 = w_1^2 +  delta_1^2, где

delta_1^2 = sum_{i=2}^{m} w_i^2 + delta^2

Потом, если мы возьмём первый фактор, ошибка будет равна

D(target - predict_1) =sigma_1^2=  w_1cdot z_1^2 + delta_1^2

Логично брать новый фактор f_i в линейный прогноз тогда, когда z_i <1.

В итоге ошибка конечного прогноза равна

  D(target - predict_{m}) = sigma_m^2=sum_{i=1}^m min(1, z_i^2)cdot w_i^2+delta^2

Итак в случае линейного прогноза и известных весов w_iможно брать или не брать слагаемые, и логично быть слагаемые w_icdot f_i с теми факторами f_i, в которых дисперсия шума меньше, чем дисперсия самого незашумлённого фактора g_i. Но ведь кроме вариантов «брать» и «не брать» есть ещё вариант брать с другим, меньшим по модулю весом. Кроме того, модель прогноза не обязательно должна быть линейной. Правильное решение задачи выглядит сложнее, в действительности мы всегда можем извлечь пользу из фактора, если в нём есть новая информация, даже если он очень зашумлён.

Пусть у нас есть два несмещённых прогноза predict_0 и predict_1и их ошибки равны

 ; mse_0 = D(target - predict_0) = sigma_0^2,\ mse_1=D(target - predict_1) = sigma_1^2.

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

predict_c = {m_0 cdot predict_0 + m_1 cdot predict_1 over  m_0 + m_1},\  m_0 = 1/sigma_0^2,; m_1 = 1/sigma_1^2.

И погрешность этого прогноза была бы равна

sigma_c^2 = {1 over 1/sigma_0^2 + 1/sigma_1^2}

То есть m_c = 1/sigma_c^2 = m_0 + m_1.

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

В нашем случае оценки зависимы, в ошибках

target - predict_0 = -w_1cdot g_1 + nu_1, \ target - predict_1= - w_1cdot  zeta_i + nu_1

есть общая часть nu_1– общая неустранимая ошибка для этих двух прогнозов. При комбинировании мы как бы комбинируем только две независимые части -w_1cdot g_1 и -w_1 cdot zeta_1 и получаем ответ

predict_c= F_0 + {m_0cdot 0 + m_1cdot (w_1 cdot f_i) over m_0 + m_1}=\=F_0 + {1over m_0/m_1+1}cdot w_1cdot f_1

где

m_0 = 1/w_1^2,; m_1= 1/(w_1^2cdot z_1^2).

Итого:

predict_c = F_0 + {1 over z_1^2+1}cdot w_1cdot f_1

И ошибка у такого прогноза будет

sigma_c^2 = w_1^2 {z_1^2over z_1^2 + 1} +   delta_1^2

что меньше погрешностей sigma_0^2=w_1^2 + delta_1^2 и sigma_1^2=w_1cdot z_1^2 + delta_1^2 обоих прогнозов predict_0 и predict_1.

Итоговая формула наилучшего линейного прогноза выглядит так

predict = sum_{i=1}^m {1over z_i^2 + 1} cdot w_i cdot f_i

А ошибка этого прогноза равна

  D(target - predict) = sigma^2=sum_{i=1}^m {z_i^2over z_i^2+1} cdot w_i^2+delta^2

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

Веса можно оценить по формуле

w'_i= mathrm{avg}(target cdot f_i)/mathrm{avg}(f_i^2)

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

M(target cdot f_i) =\= M((ldots + w_icdot g_i)cdot (g_i + zeta_i)) =\= M(w_icdot g_i^2) = w_i

А средний квадрат фактора равен

M((g_i + zeta_i)^2) = M(g_i^2 + zeta_i^2) = 1 + z_i^2

Таким образом выражение w'_i = mathrm{avg}(target cdot f_i)/mathrm{avg}(f_i^2)в пределе равно w_i / (1+z_i^2), то есть в точности тому, что нужно подставить в формулу predict. Величина w'_iбудет вычислена с погрешностью, то есть будет иметь отклонение от w_i / (1+z_i^2) – из-за конечности обучающего пула (нормировка будет неидеальной, и средние величины по выборке будут отличаться от истинных мат. ожиданий). Пусть нам дана оценка среднего квадрата относительной погрешности:

epsilon_i=left( { w'_i over  w_i / (1+z_i^2)}-1right),;; d^2_i=M(epsilon_i^2).

Если мы не берём в прогноз слагаемое w'_icdot f_i, то дополнительное слагаемое к квадрату ошибки (от потери истинного слагаемого w_icdot g_i) равно w_i^2. Если же мы вместо w_icdot  g_i берём w'_icdot f_i, то дополнительное слагаемое к квадрату ошибки равно

M((w'_icdot f_i - w_icdot g_i)^2) =\=Mleft(left({w_i over 1+z_i^2}cdot (1 + epsilon_i)cdot(g_i + zeta_i) - w_icdot g_iright)^2right) = \ = w_i^2cdot Mleft(left(left({1 + epsilon_iover 1+z_i^2} - 1right)cdot g_i + left({1 + epsilon_iover 1+z_i^2}right)cdot zeta_i right)^2 right)

Здесь мы просто

  • вынесли за скобки w_i

  • разделили выражение в скобках на две части, которые являются независимыми случайными величинами со средним равным 0.

Если это выражение больше w_i^2,то фактор лучше отбросить. Выражение в итоге упрощается до

w_i^2 cdotleft({1  + d_i^2 over (1 + z_i^2)^2} - {2 over 1 + z_i^2}+1 + {1+d_i^2 over (1+z_i^2)^2} z_i^2   right) = \ =  w_i^2 cdotleft({d_i^2 + z_i^2  over 1 + z_i^2}   right)

Таким образом, если относительная погрешность d_i^2оценки w'_i больше 1, то фактор лучше отбросить. Погрешность w'_i можно оценить, например, методом бутстрэпа.

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

Квадратичная Loss-функция и Mutual Information

Обычно задача прогноза формулируется как задача минимизации ошибки:

ML-task №3.1: По данному обучающему пулу обучите модель predict=mathcal{P}(f_1, ldots), которая минимизирует среднюю ошибку Loss(predict, target). Качество прогноза будет измеряться по средней ошибке на скрытом тестовом пуле.

Но можно попробовать сформулировать задачу прогноза как задачу максимизации MI.

ML-task №3.2: По данному обучающему пулу обучите модель mathcal{P}(f_1, ldots)такую, что

Это короткая, но не совсем правильная формулировка, требующая ряд уточнений. А именно, уточнения требуют следующие моменты:

  • В каком смысле пара(predict, target)может интерпретироваться как пара зависимых случайных величин?

  • Требование несмещённости выглядит недостижимым, если множество данных, на котором мы обучаемся и тестируемся, конечно и фиксированно.

Эти моменты, в принципе, решаются просто. При тестировании мы как бы сэмплируем пример из потенциально бесконечного пула и получаем пару случайных величин (predict, target). А несмещённость нужно понимать так, как этот обычно делают в мат. статистике в регрессионных задачах – несмещённость в среднем по всем возможным обучающим и тестовым пулам, а не на данных конкретных пулах.

Подробнее

Детерминированная функцияmathcal{P}(f_1, ldots, f_k)интерпретируется как случайная величина следующим образом: мы сэмплируем строчку из тестового пула, берём из неё набор факторов (f_1, ldots, f_k) и target, подставляем факторы в функцию и получаем значение predict и пару зависимых случайных величин (predict, target).

Требование несмещённости следует воспринимать как несмещённость в среднем при рассмотрении обучающего и тестового пулов как случайных величин. То есть конкретный конечный обучающий пул однозначно даст модель со смещённым прогнозом – среднее значение разницы predict - targetне будет равно нулю в среднем на случайном примере из потенциально бесконечного тестового пула. Но можно рассмотреть разницу avg(predict) - target, где среднее берётся по всевозможным моделям, которые можно было бы получить, обучаясь на гипотетических обучающих пулах и тестируя на гипотетических тестовых пулах. Другой вариант легализации несмещённости – это требование несмещённости в пределе, при стремлении размера обучающего пула в бесконечность. Такие интерпретации несмещённости дают шанс на существование решения задачи ML-task №3.2.2. Но честно говоря, ни инженеры, ни теоретики ML не уделяют вопросу смещённости большого внимания, как это обычно делается в методах мат. статистики. ML-щики это признают, и открыто предпочитают тестировать и оценивать свои методы на искусственных или реальных задачах по метрикам ошибки, не анализируя их смещённость.

Но мне требование несмещённости потребовалось, чтобы переформулировать задачу прогноза в терминах максимизации MI.

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

Определение 3.2: Loss-функцию равную мат. ожиданию квадратичной ошибки M_{xi^2} будем называть L_2– ошибкой или квадратичной ошибкой или mse. А Loss-функцию M_{|xi|^p} будем называть L_p– ошибкой.

Утверждение 3.1: Пусть ваша модель в задаче ML-task №3.2.2 или ML-task №3.2.1 такова, что невязка xi = predict - target обладает двумя свойствами:

  1. Является нормальной величиной с нулевым средним (то есть прогноз не смещён) при любом фиксированное значении predict.

  2. Имеет дисперсию, которая не зависит от predict.

Тогда задачи ML-task №3.1 и ML-task №3.2 эквивалентны для Loss=mse, то есть задача «maximize MI» даёт тот же ответ, что и задача «minimize mse-ошибку».

Конечно, описанные в утверждении свойства редко встречаются на практике, но тем не менее, этот факт интересен.
Во-первых, эти свойства достижимы тогда, когда вектор из n факторов и target(f_1, f_2, ldots, f_n, target)является измерением многомерной нормальной величины, а именно, так и будет, если target есть линейная комбинация нескольких независимых гауссовских величин, некоторые из которых нам известны и даны как факторы. Тогда прогноз естественно тоже искать как линейную комбинацию данных факторов.
В этом случае

mathrm{MI}(target, predict) = log({sigma_t^2 / sigma_e^2}),

где sigma_e^2 и есть среднее значение квадратичной ошибки прогноза, то есть M_{xi^2} = M_{{(target - predict)}^2}, а sigma_t^2 — дисперсия target-а (см. задачу 2.8 про MI двух зависимых нормальных величин). Понятно, что в этом случае задача максимизации mathrm{MI} эквивалентна задаче минимизации M_{xi^2}=sigma_e^2

Надо понимать, что хотя чистый и не зависящий отpredict гаусс для ошибки xi на практике не встречается, тем не менее, распределение target при фиксированном predict часто выглядит как гауссовский «колокол» с центром в predict, и утверждение про эквивалентность тоже почти верно, а именно, задача максимизации MI эквивалентна задаче минимизации ошибки в L_2-подобной метрике, где усреднение делается с весом, неким образом зависящим от predict.

Точнее так — если колокол распределения rho_xi одномодальный, симметричный, с квадратичной вершинкой, то и найдётся и эквивалентная задача про максимизацию MI. Это мы обсудим ниже.

А сейчас попробуем доказать утверждение про эквивалентность задач «maximize MI» и «minimize L2» для случая, описанного в утверждении 3.1.

Введём короткие обозначения t = target, p = predict.

Доказательство основано на том, что мы просто записываем значение mathrm{MI} в таком виде:

mathrm{MI}(p, t) = H(t) - int H(t | p = x) rho_{p}(x) dx

Величина H(t) не зависит от прогноза и равна C + 0.5cdot log(D(t)). Величина H(t | p = x) равна C + 0.5cdot log(D(t|p=x)), где как D(t|p=x) мы обозначили дисперсию прогнозируемой величины при заданном значении прогноза. Согласно второму условию утверждения она константа. В случае несмещённого прогноза этоM_{xi^2} (то есть как раз квадратичная ошибка) при условии p=x.

Выражение int ldots rho_{p}(x) dx — суть просто усреднение по пулу.

Если M_{xi^2} не зависит от значения p=x, то мы получим нужное утверждение. Действительно,

mathrm{MI}(p, t) = С +  0.5cdot log(D(t)) - int (C +  0.5 log(mathrm{mse})) rho_{p}(x) dx = \ = 0.5 cdot log(D(t) / mathrm{mse})

Максимизация этого выражения равносильна минимизации mse. Конец доказательства.

Для большего понимания распишем более подробно выражение

mathrm{MI}(p, t) =\= H(t) - int H(t | p = x) rho_{p}(x) dx =\=  H(t) - int rho_{t}(y | p = x) log(1/rho_{t}(y | p = x)) rho_{p}(x) dx dy =\  = H(t) - int log(1/rho_{t}(y | p = x)) cdot rho_{t,p}(x, y) dx dy

Интегрирование

int ldots (rho_{t,p}(x, y) dx dy

снова соответствует просто усреднению по пулу, а выражение log(1/rho_{t}(y | p = x)) и есть то, что нужно интерпретировать как значение ошибки. Если перейти к дискретизированному случаю, и вспомнить про кодирование Хаффмана, то log(1/P(t= y | p = x)) есть количество бит, которое требуется для записи кода значения t=y(с точностью до varepsilon) при условии, что известно значение p=x. Если прогноз не смещен, и он нам дан, то, чтобы закодировать t с некоторой заданной точностью, проще записать код для невязки xi=t- p, которая есть случайная величина с нулевым средним и меньшим значением дисперсии, нежели дисперсия t. Эти рассуждения наводят ещё на одну формулировку задачи прогноза как задачи максимизации MI:

ML-task №3.3: По данному обучающему пулу обучите модель predict(f_1, ldots, f_k)такую, чтобы задача кодирования невязки target-predict с точностью до varepsilon требовала минимального количества бит.

В этой задаче имеется в виду наивное кодирование чисел, когда записываются все значащие цифры. Например, пусть фиксирована точность varepsilon = 0.000001, тогда число -0.000123456789pm 0.000001 логично кодировать как «-123». Здесь мы тратим один бит на знак (+ или — в начале), далее нас не интересуют все цифры начиная с седьмого места после десятичной точки, и мы записываем только значащие цифры, не перечисляя лидирующие нули, и мы понимаем, что последняя цифра в коде стоит на шестом места после десятичной точки. В таком кодировании чем меньше ошибка в среднем, тем меньше бит потребуется для кодирования всех ошибок на тестовом пуле. Значение varepsilon должно быть достаточно мало. Но поверх этого наивного кодирования важно применить ещё кодирование Хаффмана, чтобы различие в вероятностях разных ошибок позволило дополнительно сжать данные.

В этой формулировке мы как бы невольны назначать свою Loss-функцию. Некоторые склонны рассматривать это как критерий правильного выбора Loss-функции: если вы минимизируйте mse, то хорошо, когда log(1/P(t= y | p = x)) с ростом обучающего и тестового пулов стремится к A cdot (y  - x)^2 + B, то есть распределение невязки близко к нормальному, а когда вы минимизируйте L1— ошибку, то хорошо, если log(1/P(t= y | p = x)) стремится к Acdot |y  - x| + B, то есть распределение невязки близко к распределению Лапласа. В общем случае ожидается приближённое равенствоP(t = y| p=x) approx e^{-A cdot Loss(y,x) + B} и если это не так, то либо вы «недоделали ML», либо ваша задача не относится к классу хороших.

Видимо, при определённых условиях, верны рассуждения и в другую сторону. Если у вас есть идеальное решение задачи ML-task №3.3, то изучив распределение невязки для её решения, вы узнаёте какую Loss-функцию нужно использовать, чтобы извлечь максимум информации о сигнале.

Log-Likelihood и Mutual Information

Пример задачи из жизни: Реализуйте детерминированную функцию
mathcal{P}(пол, возраст, последние 10 поисковых запросов, ldots)
прогнозирующую вероятность p того, что пользователь кликнет на данное рекламное объявление.

Задачи такого типа, когда нужно спрогнозировать булеву величину (величину, принимающую значения «Да» или «Нет»), называются задачами бинарной классификации (Binary Classification Problems).
Для их решения обычно используется метод максимизации правдоподобия, а именно, предполагается, что прогноз predict должен возвращать вещественное число от 0 до 1, соответствующее вероятности того, что target = 1(пользователь кликнет на рекламу). При обучении прогнозатора на обучающем множестве можно подбирать такие значения внутренних параметров модели (aka весов), на которых достигается максимальная вероятность видеть то, что находится в обучающем множестве.

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

  • во время обучения алгоритм не «видит» данных из тестового пула;

  • прогнозаторы сравниваются по оценке вероятностей событий на тестовом пуле.

Вероятность видеть то, что мы видим в пуле, в применении к модели называется правдоподобием модели (model likelihood). То есть мы говорим ‘вероятность события’ и ‘правдоподобие модели’, но не говорим ‘вероятность модели’ или ‘правдоподобие событий’.

Вот пример тестового пула из трёх строчек:

i

f1

f1

f1

target

predict=
P(target=1)

P(target=0)

1

0

p_1

1 - p_1

2

1

p_2

1-p_2

3

0

p_3

1-p_3

mathrm{Likelihood}=(1-p_1)cdot p_2cdot(1-p_3)

Удобнее оперировать значением LogLikelihood:

mathrm{LogLikelihood} = log((1 - p_1)cdot p_2cdot (1 - p_3)) = log(1 - p_1) + log (p_2) + log(1 - p_3)

Строчки с target=1 входят в суммарный LogLikelihood как log(predict), а строчки с target = 0— как log(1 - predict).

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

ML-task №3.4: Дан обучающий лог. Найдите детерминированную функцию predict(f_1,f_2,ldots), чтобы значение mathrm{LogLikelihood}(predict, target) на тестовом пуле было максимально.

Сформулируем аналогичную задачу в терминах максимизации MI.

ML-task №3.5: Найдите детерминированную функцию p=predict(f_1,f_2,ldots),чтобы случайная дискретная величина xi принимающая значения 0 и 1 с вероятностями {p,; 1 - p}, имела максимально возможное значение взаимной информации mathrm{MI}(xi, target) с величиной target.

Утверждение 3.2: Если ваша модель такова, что xi = predict - target является случайной величиной с нулевым средним (то есть прогноз не смещён), тогда задачи ML-task №3.2.2 и ML-task №3.2.1 эквивалентны, то есть задача «maximize MI(target, predict)» даёт тот же ответ, что и задача «maximize LogLikilihood».

То, что прогноз не смещён или, другими словами, не требует калибровки не является сложным требованием. В задачах классификации, если вы используете Gradient Boosted Trees или нейросети с правильными гиперпараметрами (скорость обучения, число итераций) прогноз получается несмещённым. А именно, если вы берёте события сpredict in [x - varepsilon, x + varepsilon], то получаете множество событий с

p_{fact} = { mathrm{number_of_lines_with_target_1} over  mathrm{number_of_lines}} = {mathrm{clicks} over  mathrm{impressions}}approx x

В задачах классификации строчки, у которых target = 1 мне привычно называть кликами (clicks), а строчки с target = 0 – некликами (notclicks). Клики + неклики дают множество всех событий, которые называются показами (impressions). Фактическое отношение кликов к показам называется CTR — Click Through Rate.

Чтобы доказать утверждение 3.2, давайте воспользуемся решением задачи:

Задача 3.10. Как по N измерениям двух случайных величин — дискретной величины xi, принимающей значения от 1 до M, и булевой случайной величины  mu, оценить значение mathrm{MI}(xi,; mu)?

Можно думать про эту задачу так: xi — это некоторая категориальная величина про рекламное объявление, её значения интерпретировать как идентификатор класса, а mu = mathrm{IsClick} — это индикатор того, был ли клик. Данные про измерение пар (xi, mu) — это просто лог кликов и некликов.

Для решения этой задачи удобно воспользоваться формулой

mathrm{MI}(xi,mu) = H(mu) - H(mu | xi)=H(mu)-{sum}_i P(xi=i)cdot H(mu|xi=i)

Значение H(mu)=H(mathrm{ctr}_{0}), где mathrm{ctr}_{0} = mathrm{ctr}_{total} = mathrm{clicks}_{total} / mathrm{impressions}_{total} — средний CTR по всему логу.

Значение H(mu | xi = i) = H(mathrm{ctr}_i), где mathrm{ctr}_{i} = mathrm{clicks}_{i} / mathrm{impressions}_{i} — CTR по событиям, в которых xi = i, что естественно называть CTR в классе i. Подставляем это в формулу и получаем:

.

mathrm{MI}(mu, xi) = H(mu) - H(mu | xi) = \ = H(mathrm{ctr}_{0}) - displaystylesum_{iin classes} { mathrm{impressions}_i over mathrm{impressions}_{total} } cdot H(mathrm{ctr}_i) = \  = - (mathrm{ctr}_{0}cdot log(mathrm{ctr}_{0}) + (1 - mathrm{ctr}_{0})cdot log(1 - mathrm{ctr}_{0})) +\  +  displaystylesum_{iin classes} {mathrm{impressions}_i over mathrm{impressions}_{total}} cdot ( mathrm{ctr}_{i}cdot log(mathrm{ctr}_{i}) + (1 - mathrm{ctr}_{i})cdot log(1 - mathrm{ctr}_{i}))

Давайте займёмся первым слагаемым и запишем его в виде суммы по классам:

H(mu) =  - left(displaystyle{mathrm{clicks}_{total} over mathrm{impressions}_{total}} log(mathrm{ctr}_{0}) + displaystyle{mathrm{notclicks}_{total} over mathrm{impressions}_{total}}cdot log(1 - mathrm{ctr}_{0})right)  = \ = - displaystylesum_{i in classes}left({mathrm{clicks}_i over mathrm{impressions}_{total}} log(mathrm{ctr}_{0}) + {mathrm{notclicks}_i over mathrm{impressions}_{total}}cdot log(1 - mathrm{ctr}_{0})right) = - displaystylesum_{i in classes}left({mathrm{ctr}_i cdot mathrm{impressions}_i over mathrm{impressions}_{total}} log(mathrm{ctr}_{0}) + {(1 - mathrm{ctr}_i)cdot mathrm{impressions}_i over mathrm{impressions}_{total}} log(1 - mathrm{ctr}_{0})right).

Подставив это в выражение для MI, и поместив обе части в одну сумму по классам, получим итоговое выражение для MI:

displaystylesum_{i in classes} { mathrm{impressions}_i over mathrm{impressions}_{total}} cdot left(mathrm{ctr}_{i}cdot logleft({mathrm{ctr}_{i} over mathrm{ctr}_{0}}right) + (1 - mathrm{ctr}_{i})cdot logleft({1 - mathrm{ctr}_{i} over 1 - mathrm{ctr}_{0}}right)right) = \ = displaystyle{1 over mathrm{impressions}_{total}} displaystylesum_{e in mathrm{impressions}} mathrm{ctr}_{i(e)}cdot log{left({mathrm{ctr}_{i(e)} over mathrm{ctr}_{0}}right)} + (1 - mathrm{ctr}_{i(e)})cdot log{left({1 - mathrm{ctr}_{i(e)} over 1 - mathrm{ctr}_{0}}right)}

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

В случае, когда классы есть просто корзины прогноза, и прогноз не смещён, то есть средний прогноз в корзине совпадает с реальным mathrm{ctr} в корзине, имеем:

mathrm{LogLikelihood(predict)}=\ = sum_{ein mathrm{impressions}}mathrm{ctr}_{i(e)}cdot log(mathrm{ctr}_{i(e)}) + (1-mathrm{ctr}_{i(e)})cdotlog(1-mathrm{ctr}_{i(e)})=\=

А если прогноз равен константе mathrm{ctr}_0 , то выражение для правдоподобия равно

mathrm{LogLikelihood}(mathrm{predict}=mathrm{ctr}_0)=\ = sum_{ein mathrm{impressions}}mathrm{ctr}_{0}cdot log(mathrm{ctr}_{0}) + (1-mathrm{ctr}_{0})cdotlog(1-mathrm{ctr}_{0})=\= sum_{ein mathrm{impressions}}mathrm{ctr}_{i(e)}cdot log(mathrm{ctr}_{0}) + (1-mathrm{ctr}_{i(0)})cdotlog(1-mathrm{ctr}_{0})

Замена mathrm{ctr}_0 на mathrm{ctr}_{i(e)} валидна, так как выражения c логарифмами константны и среднее значение mathrm{ctr}_{i(e)} по всем impressions равно mathrm{ctr}_0.

B здесь мы видим, что итоговое выражение для MI есть просто нормированная разница двух выражений для LogLikelihood:

mathrm{MI}(xi, mu) =  mathrm{MI}(predict, mu) =\  {1 over mathrm{impressions}_{total}}cdot (mathrm{LogLikelihood}(predict) - mathrm{LogLikelihood}(predict = mathrm{ctr}_{0})

То есть mathrm{MI}(predict, target) в случае несмещённого прогноза есть линейная функция от LogLikelihood.

Кстати, величина mathrm{LogLikelihood}(predict_1) - mathrm{LogLikelihood}(predict_2) называется Log Likelihood Ratio двух прогнозаторов, и её естественно нормировать на количество событий (impressions). В качестве predict_2часто берут какой-нибудь простейший прогноз, в нашем случае, это константный прогноз. Часто имеет смысл мониторить график величины LogLikelihoodRatio / impressions, а не LogLikelihood / impressions, беря в качестве бейслайнового прогноза predict_2 какой-либо робастный (простой, надёжный неломающийся) прогноз, основанный на небольшом числе факторов. Также иногда можно брать  mathrm{LogLikelihoodRatio} / mathrm{impressions} cdot mathrm{ctr}_0^{gamma}, для некоторого gamma, чтобы устранить корреляцию или антикорреляцию с mathrm{ctr}_0 и лучше видеть точки поломок прогноза.

Итак, для несмещённого прогноза величина MI между прогнозом и сигналом равна удельному (нормированному на число строк) Log Likelihood Ratio вашего прогноза и наилучшего константного прогноза.

Оценка Mutual Information = Machine Learning

Два утверждения — 3.1 и 3.2 — говорят о том, что Mutual Information является метрикой качества, которая в некоторых допущениях соответствует двум метрикам в задачах прогноза — квадратичной ошибке в задаче прогноза вещественной величины с нормальным распределением и LogLikelihood в задаче бинарной классификации.

Сама величина Mutual Information не может использоваться как Loss-функция, потому что она не является метрикой на данных, то есть не представляется как сумма значений Loss-функции по элементам пула. Смысл «соответствия» можно пояснить так: практически все Loss-функции в конечном итоге занимаются максимизацией Mutual Information predict-а и target-а, но с разными дополнительными добавками.

Возможно, свет на это прольют следующие два утверждения (опять без доказательства):

Утверждение 3.3: Если у вас есть прогноз на факторах f_1,ldots,f_k, и есть новый фактор f_{k+1}, такой что mathrm{MI}({f_{k+1}, predict}, target) > mathrm{MI}(predict, target), то при достаточно большом обучающем логе этот фактор уменьшит вашу Loss-функцию, если она нормальная. Loss-функция называется нормальной если она падает, когда вы в любом элементе пула значение predict приближаете к target. Квадратичная ошибка mse, L_p-ошибка, LogLikelihood — это всё нормальные Loss-функции. Далее будем считать, что Loss-функция нормальная.

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

Утверждение 3.4 (требует уточнения условий): Если вы нашли такое изменение весов (ака внутренних параметров) вашей модели, которое увеличивает mathrm{MI}(predict, target), то после надлежащей калибровки прогноза вырастет и ваша Loss-функция.

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

Задача 3.11: Докажите, что для случая, когда target=mathcal{P}(f_1,f_2,ldots,w_1,w_2,ldots) + nu, где nu случайный шум, и модель соответствует реальности, при увеличении обучающего лога веса стремятся к правильным значениям весов для обоих Loss-функций L_1иL_2. При этом в случае, если шум имеет симметричное распределение, обе Loss-функций дают несмещённый прогноз.

Задача 3.12: Приведите пример модели и реального target, когда Loss-функции L_1 и L_2дают разные прогнозы даже на очень большом обучающем пуле, так что один нельзя перевести во второй никаким монотонным преобразованием (один прогноз нельзя получить калибровкой второго).

Утверждение 3.5: Задача оценки mi = mathrm{MI}({f_1,ldots, f_k}, target) эквивалентна задаче построения прогноза predict=mathcal{P}(f_1,ldots, f_k) в какой-то Loss-функции.

Это утверждение говорит, что истинное значение mi примерно равно supremum mathrm{MI}(predict, target)всем парам (ML-метод, Loss-функция), где
ML-метод ::= («структура модели», «метод получения весов модели»),
«метод получения весов модели» ::= («алгоритм», «гиперпараметры алгоритма»),
а равенство тем точнее, чем больше обучающий пул.

Собственно та пара (ML-метод, Loss-функция), которая даст самое большое значение mathrm{MI}(predict, target) близкое к mi, и есть модель, самая близкая к реальности, то есть тому, как на самом деле связаны f_1,ldots,f_k и target.

Итак: Мы сформулировали две глубокие связи между задачами прогноза и MI:

  • Во-первых, максимизация MI в некотором смысле соответствует задаче минимизации какой-то Loss-функции

  • Во-вторых, «оценка MI» = «ML», а именно, задача оценки mathrm{MI}({f_1,ldots, f_k}, target) эквивалентна задаче построения прогноза mathcal{P}(f_1,ldots, f_k) для какой-то Loss-функции. И собственно, ML-метод и Loss-функция, которые дают максимум mathrm{MI}(predict, target) и является наиболее правдоподобной моделью.


  Перевод


  Ссылка на автора

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

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

В этом руководстве вы узнаете показатели производительности для оценки прогнозов временных рядов с помощью Python.

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

После завершения этого урока вы узнаете:

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

Давайте начнем.

Ошибка прогноза (или остаточная ошибка прогноза)

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

Это называется остаточной ошибкой прогноза.

forecast_error = expected_value - predicted_value

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

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

expected = [0.0, 0.5, 0.0, 0.5, 0.0]
predictions = [0.2, 0.4, 0.1, 0.6, 0.2]
forecast_errors = [expected[i]-predictions[i] for i in range(len(expected))]
print('Forecast Errors: %s' % forecast_errors)

При выполнении примера вычисляется ошибка прогноза для каждого из 5 прогнозов. Список ошибок прогноза затем печатается.

Forecast Errors: [-0.2, 0.09999999999999998, -0.1, -0.09999999999999998, -0.2]

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

Средняя ошибка прогноза (или ошибка прогноза)

Средняя ошибка прогноза рассчитывается как среднее значение ошибки прогноза.

mean_forecast_error = mean(forecast_error)

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

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

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

expected = [0.0, 0.5, 0.0, 0.5, 0.0]
predictions = [0.2, 0.4, 0.1, 0.6, 0.2]
forecast_errors = [expected[i]-predictions[i] for i in range(len(expected))]
bias = sum(forecast_errors) * 1.0/len(expected)
print('Bias: %f' % bias)

При выполнении примера выводится средняя ошибка прогноза, также известная как смещение прогноза.

Bias: -0.100000

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

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

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

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

mean_absolute_error = mean( abs(forecast_error) )

кудаабс ()делает ценности позитивными,forecast_errorодна или последовательность ошибок прогноза, иимею в виду()рассчитывает среднее значение.

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

from sklearn.metrics import mean_absolute_error
expected = [0.0, 0.5, 0.0, 0.5, 0.0]
predictions = [0.2, 0.4, 0.1, 0.6, 0.2]
mae = mean_absolute_error(expected, predictions)
print('MAE: %f' % mae)

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

MAE: 0.140000

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

Средняя квадратическая ошибка

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

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

mean_squared_error = mean(forecast_error^2)

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

from sklearn.metrics import mean_squared_error
expected = [0.0, 0.5, 0.0, 0.5, 0.0]
predictions = [0.2, 0.4, 0.1, 0.6, 0.2]
mse = mean_squared_error(expected, predictions)
print('MSE: %f' % mse)

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

MSE: 0.022000

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

Среднеквадратическая ошибка

Средняя квадратичная ошибка, описанная выше, выражается в квадратах единиц прогнозов.

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

rmse = sqrt(mean_squared_error)

Это можно рассчитать с помощьюSQRT ()математическая функция среднего квадрата ошибки, рассчитанная с использованиемmean_squared_error ()функция scikit-learn.

from sklearn.metrics import mean_squared_error
from math import sqrt
expected = [0.0, 0.5, 0.0, 0.5, 0.0]
predictions = [0.2, 0.4, 0.1, 0.6, 0.2]
mse = mean_squared_error(expected, predictions)
rmse = sqrt(mse)
print('RMSE: %f' % rmse)

При выполнении примера вычисляется среднеквадратичная ошибка.

RMSE: 0.148324

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

Дальнейшее чтение

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

  • Раздел 3.3 Измерение прогнозирующей точности, Практическое прогнозирование временных рядов с помощью R: практическое руководство,
  • Раздел 2.5 Оценка точности прогноза, Прогнозирование: принципы и практика
  • scikit-Learn Metrics API
  • Раздел 3.3.4. Метрики регрессии, scikit-learn API Guide

Резюме

В этом руководстве вы обнаружили набор из 5 стандартных показателей производительности временных рядов в Python.

В частности, вы узнали:

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

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

Понравилась статья? Поделить с друзьями:
  • Средняя квадратическая ошибка как найти
  • Среднеквадратичная ошибка excel
  • Среднеквадратичная ошибка мнк
  • Средняя квадратичная ошибка пример
  • Средняя квадратическая ошибка измерения угла