Pyro прогнозирование сверху вниз

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

Вы когда-нибудь задумывались, как учесть неопределенности в прогнозах временных рядов?
Вы когда-нибудь думали, что должен быть способ генерировать точки данных из ранее увиденных данных и делать выводы о достоверности? Я знаю, что есть.
Если вы хотите создавать модели, которые фиксируют вероятности и сохраняют уверенность, мы рекомендуем использовать фреймворк вероятностного программирования, такой как Pyro.
В предыдущей статье мы рассмотрели NGBoosting и применили его к задаче прогнозирования M5 на Kaggle. Вкратце: задача прогнозирования M5 требует от нас предсказать, как продажи товаров Walmart будут развиваться с течением времени. Он предоставляет данные за 4–5 лет для товаров разных категорий из разных магазинов в разных штатах и просит нас прогнозировать 28 дней, о которых у нас нет информации. В качестве обзора задачи и набора данных мы по-прежнему рекомендуем этот удивительный блокнот. В прошлый раз мы пришли к выводу, что с Pyro мы достигли лучших результатов, и вот простой пример того, как мы этого достигли.

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

Sidenote: One more elegant way to do forecasting is with hierarchical models. These models allow accounting for different distributions over different categories in your data. You can assign individual priors and make use of information that you have about your model. The elegance of such an approach will be covered in another post.

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

prediction = regressor + trend + seasonal + motion + bias

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

Составление модели

Фактическая часть модели фиксируется регрессором. Эта часть ищет правильные веса модели и сравнивает их с функциями, которые поступают непосредственно из ввода данных. Каждый вес лежит на классической гауссовой кривой. Это означает, что мы выбираем из нормального распределения с центром вокруг нуля со стандартным отклонением 1 (μ = 0, σ² = 1).

weight = pyro.sample("weight", dist.Normal(0, 1).expand(
        [feature.size(-1)]).to_event(1))
 regressor = (weight * feature).sum(-1)

Второй параметр - это тренд. Мы консультируемся с нашим вероятностным зеркалом и затем решаем, что оно должно исходить из логнормального распределения, такого что l∼LN с μ = -2 и σ² = 1. Ключевым элементом является то, что мы учитываем временную характеристику, которую мы предоставляем модели.

trend_coef = pyro.sample("trend", dist.LogNormal(-2, 1))
trend = trend_coef * time

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

with pyro.plate("day_of_week", 7, dim=-1):
    seasonal = pyro.sample("seasonal", dist.Normal(0, 5))
seasonal = periodic_repeat(seasonal, duration, dim=-1)

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

Наконец, каждой профессиональной модели машинного обучения нужен некоторый обучаемый параметр, который учитывает тот факт, что иногда все бывает не так хорошо, как вы ожидаете - мы называем это предвзятостью. Для нашего смещения (b) мы решаем, что оно происходит из нормального распределения, такого что b∼N, где μ = 0 и σ² = 10 или loc = 0 и scale = 10, если вы хотите говорить о Pyro.

bias = pyro.sample("bias", dist.Normal(0, 10))

Прежде чем продолжить, давайте посмотрим на распределение наших параметров.

На этом вы, внимательный читатель, можете сделать вывод. У вас есть все необходимое, и мы согласны. То, что у нас сейчас есть, достаточно, чтобы дать нам оценку и сделать прогнозы на пару дней в будущее.
Но мы можем сделать это лучше. То, что мы прогнозируем до сих пор, является действительно жестким выражением. Говорят, что модель соответствует входным данным, но мы не допускаем свободы действий.
Давайте предположим, ради аргументации, что все люди в Калифорнии решают, что им лучше идти за покупками в полдень в четверг, а не в субботу; или наступает кризис, и все идут готовиться к апокалипсису. Наша модель, возможно, правильно предсказала все предыдущие субботы, но тот четверг действительно сбивает ее с толку. Чтобы учесть это, мы вводим так называемую повторную параметризацию. Подробный обзор этой темы см. В [2].
Чтобы сделать его еще более надежным, мы проводим повторную параметризацию с помощью другой повторной параметризации - почему мы это делаем, объясняется здесь.
Собирая все вместе, мы получаем дрейф, который определяет движение нашей модели. Это переводится на Pyro как таковое:

drift_stability = pyro.sample("drift_stability", dist.Uniform(1, 2))
drift_scale = pyro.sample("drift_scale", dist.LogNormal(-17, 5))
with self.time_plate:
   with poutine.reparam(config={"drift": LocScaleReparam()}):
      with poutine.reparam(config={"drift": SymmetricStableReparam()}):
         drift = pyro.sample("drift", dist.Stable(drift_stability, 0, drift_scale))
motion = drift.cumsum(dim=-1)

Теперь у нас есть все необходимое для работы нашей модели.

Ваша модель на работе

Для фактической подгонки модели мы используем механизм стохастического вариационного вывода Pyro или SVI. Этот алгоритм позволяет нам эффективно вычислить апостериорное распределение за разумный промежуток времени. В нашем случае это сделано для максимизации нижней границы доказательств или ELBO. Один из упрощенных способов думать об ELBO - это наблюдать за двумя разными распределениями, текущим и предыдущим. Если распределения сильно отличаются, мы получаем низкое значение ELBO. С другой стороны, если я установил хороший дистрибутив, мой следующий может быть таким же хорошим, хорошо согласуется с уже наблюдаемым, и я получаю высокое значение ELBO. Некоторые могут даже сказать, что расхождение обоих распределений сведено к минимуму. Чтобы прочитать об этом подробнее, обратимся к исходной статье для SVI (см. [3]). А пока давайте предположим, что это просто работает, потому что умные исследователи реализовали это в достаточной степени.

We will go into the interplay of Stochastic Variational Inference and ELBO (or minimizing KL) in another article.

Теперь, когда мы собрали рабочую модель, выбрали наш алгоритм для вычисления наших апостериорных распределений, нам нужно установить гиперпараметры, загрузить данные и начать обучение. Что касается обработки данных, мы уделяем большое внимание реализации M5 Starter Kit, найденной на github, которую команда Pyro PPL щедро предоставляет. Это простой способ получить агрегаты, например данные о продажах. Не стесняйтесь реализовать это самостоятельно. Это верный способ научиться обходить тензоры. В нашей модели для начинающих нас интересуют только данные о продажах с течением времени.
Теперь мы устанавливаем наши параметры и создаем экземпляр прогнозиста следующим образом:

forecaster_opt = {
    "learning_rate": 0.1,
    "learning_rate_decay": 0.1,
    "clip_norm": 10,
    "num_steps": 3501,
    "log_every": 100,
}
forecaster = Forecaster(TopDownModel(), data, covariates[:-28], **forecaster_opt)
samples = forecaster(data, covariates, num_samples=1000).exp().squeeze(-1).cpu()
pred = samples.mean(0)

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

Мы видим, что наши потери уменьшаются. Это хорошо и говорит нам о том, что вычисленные распределения по нашему прогнозу со временем становятся лучше. Потери здесь являются частным от отрицательного ELBO относительно данных и части документации прогнозиста (см. [4]).

Настоящая неуверенная вещь - с уверенностью и всем остальным

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

Чтобы вычислить достоверность отдельных шагов, мы запрашиваем у нашей модели 100 выборок и вычисляем квантили (0,1 и 0,9) для этих выборок. Это красная область, которую вы видите на рис.2. Поскольку эта модель находится в вашем распоряжении, вы можете запросить дополнительные образцы, если хотите сделать ее более точной. Продажи на рисунке - это агрегированные зарегистрированные значения для лучшего соответствия модели. Чтобы получить значения в более реалистичном диапазоне, преобразуйте значения с помощью numpy's exp.

Обобщить

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

Полный код модели

class TopDownModel(ForecastingModel):
  """
  Top-Down Hierarchical Forecasting Model
  """
  def model(self, zero_data, covariates):
    # check univariate data
    assert zero_data.size(-1) == 1 
    duration = zero_data.size(-2)
    time, feature = covariates[..., 0], covariates[..., 1:]
    bias = pyro.sample("bias", dist.Normal(0, 10))
    trend_coef = pyro.sample("trend", dist.LogNormal(-2, 1))
    trend = trend_coef * time
    weight = pyro.sample("weight", dist.Normal(0, 1).expand(
        [feature.size(-1)]).to_event(1))
    regressor = (weight * feature).sum(-1)
    # weekly seasonality as independent events
    with pyro.plate("day_of_week", 7, dim=-1):
      seasonal = pyro.sample("seasonal", dist.Normal(0, 5))
    seasonal = periodic_repeat(seasonal, duration, dim=-1)
    drift_stability = pyro.sample("drift_stability", dist.Uniform(1, 2))
    drift_scale = pyro.sample("drift_scale", dist.LogNormal(-17, 5))
    # introduce drift
    with self.time_plate:
          # We combine two different reparameterizers: the inner SymmetricStableReparam
          # is needed for the Stable site, and the outer LocScaleReparam improves inference.
          with poutine.reparam(config={"drift": LocScaleReparam()}):
              with poutine.reparam(config={"drift": SymmetricStableReparam()}):
                  drift = pyro.sample("drift",
                                      dist.Stable(drift_stability, 0, drift_scale))
    motion = drift.cumsum(dim=-1)
    # predict
    prediction = regressor + trend + seasonal + motion + bias
    # Pyro Forecast is multivariate - univariate timeseries is needed
    prediction = prediction.unsqueeze(-1)
    # heavy tail nose to account for outliers
    stability = pyro.sample("noise_stability", dist.Uniform(1, 2).expand([1]).to_event(1))
    skew = pyro.sample("noise_skew", dist.Uniform(-1, 1).expand([1]).to_event(1))
    scale = pyro.sample("noise_scale", dist.LogNormal(-5, 5).expand([1]).to_event(1))
    noise_dist = dist.Stable(stability, skew, scale)
    with poutine.reparam(config={"residual": StableReparam()}):
      self.predict(noise_dist, prediction)

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

# sample from the fitted model to forecast data
samples = forecaster(data[T0:T1], covariates[T0:T2], num_samples=100)
samples.clamp_(min=0)
p10, p50, p90 = quantile(samples[:, :], (0.1, 0.5, 0.9)).squeeze(-1)
fig, ax = plt.subplots()
ax.set_title("aggregated (log1p) sales over 42 days")
ax.plot(data[T1:T2].squeeze(-1).cpu().numpy(), "r-",lw=4, label="truth")
ax.fill_between(np.arange(0, (T2-T1)), p10.cpu().numpy(), p90.cpu().numpy(), alpha=0.3, color="red", label="confidence")
ax.plot(samples.squeeze(-1).T.cpu().numpy()[:, 10:13], color="green", alpha=0.35, label="sample")
ax.plot(np.mean(samples.squeeze(-1).T.cpu().numpy()[:, :], axis=1),color="green", lw=2, alpha=0.75, label="mean")
ax.set_xlabel("t in days")
ax.set_ylabel("aggregated sales")
labels = np.arange(len(covariates)-42, len(covariates), step=8)
ax.set_xticklabels(labels)
plt.ylim((9.5, 11.2))
plt.legend(loc='lower right')
plt.show()

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

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

использованная литература

Pyro Documentation - Учебное пособие по прогнозированию I, II и III
М. Горинова, Д. Мур, М. Хоффман. Автоматическая повторная параметризация вероятностных программ. 2019 опубликовано на ArXiv
М. Хоффман и др. в Стохастический вариационный вывод. 2013 г. в JMLR
Uber Technologies. 2018. Pyro Forecaster см. Документация

Pyro прогнозирование сверху вниз | Случай применения