Предсказываем будущее с помощью библиотеки Facebook Prophet

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


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


    Сегодня же мы познакомимся с библиотекой для прогнозирования временных рядов Facebook Prophet (в переводе с английского, "пророк", выпущена в open-source 23-го февраля 2017 года), а также попробуем в жизненной задаче – прогнозировании числа постов на Хабрехабре.



    Библиотека


    Согласно статье Facebook Prophet, был разработан для прогнозирования большого числа различных бизнес-показателей и строит достаточно хорошие default'ные прогнозы. Кроме того, библиотека дает возможность, изменяя человеко-понятные параметры, улучшать прогноз и не требует от аналитиков глубоких знаний устройства предсказательных моделей.


    Давайте немного обсудим, как же работает библиотека Prophet. По сути, это additive regression model, состоящая из следующих компонент:


    $y(t) = g(t) + s(t) + h(t) + \epsilon_{t}$


    1. Сезонные компоненты $s(t)$ отвечают за моделирование периодических изменений, связанных с недельной и годовой сезонностью. Недельная сезонность моделируется с помощью dummy variables. Добавляются 6 дополнительных признаков, например, [monday, tuesday, wednesday, thursday, friday, saturday], которые принимают значения 0 и 1 в зависимости от даты. Признак sunday, соответствующий седьмому дню недели, не добавляют, потому что он будут линейно зависеть от других дней недели и это будет влиять на модель.
      Годовая же сезонность моделируется рядами Фурье.
    2. Тренд $g(t)$ — это кусочно-линейная или логистическая функция. С линейной функцией все понятно. Логистическая же функция вида $g(t) = \frac{C}{1+exp(-k(t-b))}$ позволяет моделировать рост с насыщением, когда при увеличении показателя снижается темп его роста. Типичный пример — это рост аудитории приложения или сайта.
      Кроме всего прочего, библиотека умеет по историческим данным выбирать оптимальные точки изменения тренда. Но их также можно задать и вручную (например, если известны даты релизов новой функциональности, которые сильно повлияли на ключевые показатели).
    3. Компонента $h(t)$ отвечает за заданные пользователем аномальные дни, в том числе и нерегулярные, такие как, например, Black Fridays.
    4. Ошибка $\epsilon_{t}$ содержит информацию, которая не учтена моделью.

    Подробнее про алгоритмы можно прочитать в публикации Sean J. Taylor, Benjamin Letham "Forecasting at scale".


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



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


    MAPE (mean absolute percentage error) — это средняя абсолютная ошибка нашего прогноза. Пусть $y_{i}$ — это показатель, а $\hat{y_{i}}$ — это соответствущий этой величине прогноз нашей модели. Тогда $e_{i} = y_{i} - \hat{y_{i}}$ — это ошибка прогноза, a $p_{i} =\frac{e_{i}}{y_{i}}$ — это относительная ошибка прогноза.


    $MAPE = mean(|p_{i}|)$


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


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


    $MAE = mean(|e_{i}|)$


    Cтоит сказать пару слов о тех алгоритмах, с которыми сравнивали Prophet в публикации, тем более, большинство из них очень простые и их часто используют как baseline:


    • naive — наивный прогноз, когда мы прогнозируем все последующие значения последней точкой;
    • snaive - seasonal naive — такой прогноз подходит для данных с явно выраженной сезонностью. Например, если мы говорим о показатели с недельной сезонностью, то для каждого последующего понедельника мы будем брать значение за последний понедельник, для вторника — за последний вторник и так далее;
    • mean — в качестве прогноза берется среднее значение показателя;
    • arima - autoregressive integrated moving average — подробности на wiki;
    • ets - exponential smoothing — подробности на wiki.

    Практика


    Установка


    Для начала необходимо установить библиотеку. Библиотека Prophet доступна для python и R. Я предпочитаю python, поэтому использовала именно его. Для python библиотека ставится с помощью PyPi следующим образом:


    pip install fbprophet

    Под R у библиотеки есть CRAN package. Подробные инструкции по установке можно найти в документации.


    Данные


    В качестве показателя для предсказания я выбрала количество постов, опубликованных на Хабрахабре. Данные я взяла из учебного конкурса на Kaggle "Прогноз популярности статьи на Хабре". Тут рассказано подробнее о соревновании и курсе машинного обучения, в рамках которого оно проводится.


    Код для предобработки данных
    import pandas as pd
    
    habr_df = pd.read_csv('howpop_train.csv')
    habr_df['published'] = pd.to_datetime(habr_df.published)
    habr_df = habr_df[['published', 'url']]
    
    aggr_habr_df = habr_df.groupby('published')[['url']].count() 
    aggr_habr_df.columns = ['posts']
    
    aggr_habr_df = aggr_habr_df.resample('D').apply(sum)

    Для начала посмотрим на данные и построим time series plot за весь период. На таком длинном периоде удобнее смотреть на недельные точки.


    Для визуализации я как обычно буду использовать библиотеку plot.ly, которая позволяет строить в python интерактивные графики. Подробнее про нее и визуализацию в целом можно почитать в статье Открытый курс машинного обучения. Тема 2: Визуализация данных c Python.


    Код для визуализации
    from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
    from plotly import graph_objs as go
    
    # инициализируем plotly
    init_notebook_mode(connected = True)
    
    # опишем функцию, которая будет визуализировать все колонки dataframe в виде line plot
    def plotly_df(df, title = ''):
        data = []
    
        for column in df.columns:
            trace = go.Scatter(
                x = df.index,
                y = df[column],
                mode = 'lines',
                name = column
            )
            data.append(trace)
    
        layout = dict(title = title)
        fig = dict(data = data, layout = layout)
        iplot(fig, show_link=False)
    
    plotly_df(aggr_habr_df.resample('W').apply(sum), title = 'Опубликованные посты на Хабрахабре')


    Построение прогноза


    Библиотека Prophet имеет интерфейс похожий на sklearn, сначала мы создаем модель, затем вызываем у нее метод fit и затем получаем прогноз. На вход методу fit библиотека принимает dataframe с двумя колонками:


    • ds — время, поле должно быть типа date или datetime,
    • y — числовой показатель, который мы хотим предсказывать.

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


    #импортируем библиотеку
    from fbprophet import Prophet
    
    predictions = 30
    
    # приводим dataframe к нужному формату
    df = aggr_habr_df.reset_index()
    df.columns = ['ds', 'y']
    
    # отрезаем из обучающей выборки последние 30 точек, чтобы измерить на них качество
    train_df = df[:-predictions] 

    Далее создаем объект класса Prophet (все параметры модели задаются в конструкторе класса, для начала возьмем default'ные параметры) и обучаем его.


    m = Prophet()
    m.fit(train_df)

    С помощью вспомогательной функции Prophet.make_future_dataframe создаем dataframe, который содержит все исторические временные точки и еще 30 дней, для которых мы хотели построить прогноз.


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


    future = m.make_future_dataframe(periods=predictions)
    forecast = m.predict(future)

    В библиотеке Prophet есть встроенные средства визуализации, которые позволяют оценить результат построенной модели.


    Во-первых, метод Prophet.plot отображает прогноз. Честно говоря, в данном случае такая визуализация не очень показательна. Основной вывод из этого графика, который я сделала — в данных много outlier'ов. Однако если при прогнозировании будет меньше исторических точек, то по ней можно будет что-нибудь понять.


    m.plot(forecast)


    Вторая функция Prophet.plot_components, на мой взгляд, гораздо более полезная. Она позволяет посмотреть отдельно на компоненты: тренд, годовую и недельную сезонность. Если при построении модели были заданы аномальные дни/праздники, то они также будут отображаться на этом графике.


    m.plot_components(forecast)


    На графике видно, что Prophet хорошо хорошо подстроился под рост числа постов "ступенькой" в начале 2015-го. По недельной сезонности можно сделать вывод, что меньше постов приходится на воскресенье и понедельник. На графике годовой сезонности ярче всего выделяется провал активности в Новогодние каникулы, также виден спад и на майских праздниках.


    Оценка качества модели


    Давайте оценим качество алгоритма и посчитаем MAPE для последних 30 дней, которые мы предсказывали. Для расчета нам нужны наблюдения $y_{i}$ и прогнозы для них $ \hat{y_{i}}$.


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


    print(', '.join(forecast.columns))

    ds, t, trend, seasonal_lower, seasonal_upper, trend_lower, trend_upper, yhat_lower, yhat_upper, weekly, weekly_lower, weekly_upper, yearly, yearly_lower, yearly_upper, seasonal, yhat

    Прежде чем продолжить, нам нужно объединить forecast с нашими исходными наблюдениями.


    cmp_df = forecast.set_index('ds')[['yhat', 'yhat_lower', 'yhat_upper']].join(df.set_index('ds'))

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


    import numpy as np
    cmp_df['e'] = cmp_df['y'] - cmp_df['yhat']
    cmp_df['p'] = 100*cmp_df['e']/cmp_df['y']
    print 'MAPE', np.mean(abs(cmp_df[-predictions:]['p']))
    print 'MAE', np.mean(abs(cmp_df[-predictions:]['e']))

    В результате мы получили качество около 37.35%, а в среднем в прогнозе модель ошибается на 10.62 поста.


    Визуализация


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


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


    Код для визуализации
    
    # функция для визуализации построенного прогноза
    def show_forecast(cmp_df, num_predictions, num_values):
        # верхняя граница доверительного интервала прогноза
        upper_bound = go.Scatter(
            name='Upper Bound',
            x=cmp_df.tail(num_predictions).index,
            y=cmp_df.tail(num_predictions).yhat_upper,
            mode='lines',
            marker=dict(color="444"),
            line=dict(width=0),
            fillcolor='rgba(68, 68, 68, 0.3)',
            fill='tonexty')
    
        # прогноз
        forecast = go.Scatter(
            name='Prediction',
            x=cmp_df.tail(predictions).index,
            y=cmp_df.tail(predictions).yhat,
            mode='lines',
            line=dict(color='rgb(31, 119, 180)'),
        )
    
        # нижняя граница доверительного интервала
        lower_bound = go.Scatter(
            name='Lower Bound',
            x=cmp_df.tail(num_predictions).index,
            y=cmp_df.tail(num_predictions).yhat_lower,
            marker=dict(color="444"),
            line=dict(width=0),
            mode='lines')
    
        # фактические значения
        fact = go.Scatter(
            name='Fact',
            x=cmp_df.tail(num_values).index,
            y=cmp_df.tail(num_values).y,
            marker=dict(color="red"),
            mode='lines',
        )
    
        # последовательность рядов в данном случае важна из-за применения заливки
        data = [lower_bound, upper_bound, forecast, fact]
    
        layout = go.Layout(
            yaxis=dict(title='Посты'),
            title='Опубликованные посты на Хабрахабре',
            showlegend = False)
    
        fig = go.Figure(data=data, layout=layout)
        iplot(fig, show_link=False)
    
    show_forecast(cmp_df, predictions, 200)

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



    Визуально прогноз модели, кажется достаточно хорошим и разумным. Скорее всего такая низкая оценка качества объясняется аномальным высоким количеством постов 13 и 17 октября и снижением активности с 7 октября.


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


    Сравнение с ARIMA моделью


    На глаз прогноз получился вполне разумным, но давайте сравним его с классической моделью SARIMA - Seasonal autoregressive integrated moving average с недельным периодом.


    На Хабрахабре уже есть несколько статей про ARIMA-модели, всем интересующимся советую почитать их: Построение модели SARIMA с помощью Python+R и Анализ временных рядов с помощью python.


    Для построения прогноза я также вдохновлялась учебными материалами курса Прикладные задачи анализа данных на Сoursera, в котором подробно описаны ARIMA-модели и как их строить на python.


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


    Но в данном случае усилия были не напрасны и предсказание SARIMA получилось более точным: MAPE=16.54%, MAE=7.28 поста. Лучшая модель с параметрами: D=1, d=1, Q=1, q=4, P=1, p=3.



    Но и Prophet, конечно же, можно еще потюнить. Например, если предсказывать в этой библиотеке не исходный ряд, а после преобразования Бокса-Кокса, нормализующего дисперсию ряда, то мы получим прирост качества: MAPE=26.79%, MAE=8.49 поста.


    Код для прогноза ряда с преобразованием Бокса-Кокса
    from scipy import stats
    import statsmodels.api as sm
    
    def invboxcox(y,lmbda):
        if lmbda == 0:
            return(np.exp(y))
        else:
            return(np.exp(np.log(lmbda*y+1)/lmbda))
    
    train_df2 = train_df.copy().fillna(14)
    train_df2 = train_df2.set_index('ds')
    train_df2['y'], lmbda_prophet = stats.boxcox(train_df2['y'])
    
    train_df2.reset_index(inplace=True)
    
    m2 = Prophet()
    m2.fit(train_df2)
    future2 = m2.make_future_dataframe(periods=30)
    
    forecast2 = m2.predict(future2)
    forecast2['yhat'] = invboxcox(forecast2.yhat, lmbda_prophet)
    forecast2['yhat_lower'] = invboxcox(forecast2.yhat_lower, lmbda_prophet)
    forecast2['yhat_upper'] = invboxcox(forecast2.yhat_upper, lmbda_prophet)
    
    cmp_df2 = forecast2.set_index('ds')[['yhat', 'yhat_lower', 'yhat_upper']].join(df.set_index('ds'))
    
    cmp_df2['e'] = cmp_df2['y'] - cmp_df2['yhat']
    cmp_df2['p'] = 100*cmp_df2['e']/cmp_df2['y']
    np.mean(abs(cmp_df2[-predictions:]['p'])), np.mean(abs(cmp_df2[-predictions:]['e']))

    Итог


    Мы познакомились с open-source библиотекой Prophet и ее использованием для предсказания временных рядов на практике.
    Я бы не стала говорить, что эта библиотека творит чудеса и идеально предсказывает будущее. В нашем случае прогноз получился хуже стандартной SARIMA. Однако, библиотека Prophet достаточно удобная, легко кастомизируется (чего только стоит возможность добавления заранее известных аномальных дней), поэтому ее полезно иметь в своем аналитическом toolbox'e.


    Полезные ссылки


    Немного материалов для более глубокого изучения библиотеки Prophet и предсказаний временных рядов, в общем:


    Open Data Science 232,53
    Крупнейшее русскоязычное Data Science сообщество
    Поделиться публикацией
    AdBlock похитил этот баннер, но баннеры не зубы — отрастут

    Подробнее
    Реклама
    Комментарии 14
    • +5
      А в Prophet есть возможность включать дополнительные признаки?
      • +1

        В документации такой возможности не встречала.

      • +4
        Надо проверить на USDRUB
        • +1
          У меня тоже зачесались руки проверить валютную пару какую-то.
          Кстати, видел дипломную по предсказанию USD/RUB
          Там получилось предсказать дельту, но не ± направления движения
        • –3
          Стоит отметить, что построение ARIMA модели требует гораздо больших затрат по сравнению с Prophet: нужно исследовать исходный ряд, привести его к стационарному, подобрать начальные приближения и потратить немало времени на подбор гипер-параметров алгоритма


          Все то, что вы описали — не является ни проблемой, ни «большими затратами». К слову, привести временной ряд к стационарному можно c помощью невероятно сложной команды того же Eviews: D(x1). Занимает 0 целых, и сколько-то там сотых секунд времени.

          Очень сомневаюсь, что Prophet где-то бьет Ариму(нормально выполненную Ариму)
          • +2

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


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

            • 0
              Eviews — это то, что создано непосредственно для прогнозирования временных рядов. Но даже если это Питон, в чем проблема взять first difference, чтобы привести к стационарности?

              Это просто банальная разница между сегодняшним и вчерашним значением (в случае, когда речь идет о чем-то простом, вроде анализа количества постов), это сложно реализовать в питоне?)

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

              Зачем вообще нужна стационарность? Если посмотреть на любой временной ряд (практически), то его среднее — не является адекватной, ибо пересекает значения ряда всего несколько раз. Гораздо лучше работать с growth rate (темпами роста), тогда решается проблема средней.

              Арима, при условии, что мы выровняли временной ряд, и очистили его от сезонности — лучше справится с задачей, чем Prophet
              • +3
                В Eviews не так просто загрузить данные и сделать что-то гибкое и кастомное))
                И суть статьи наверное в том, что бы быстро получить результат при минимальных усилиях. Для Аримы все-таки нужно немного подумать, придумать интерпретацию и зафитить параметры. Prophet еще можно использовать как бейслайн.
                Например, приходишь к начальнику и говоришь: «Я делаю лучше прогноз, чем библиотека от самого фейсбука! » И сразу поднимаешься в глазах руководства :)
                • +1
                  Да, с подумать — не поспоришь, но сама реализация на самом деле — не сложная, это единственное, что меня зацепило. Ну в том плане, что да, Профет быстрее, примерно такой же по результатам, все дела. И это тоже является своим плюсом.

                  Единственное, что меня зацепило, это фраза "гораздо больших затрат". Ну не гораздо совсем) Плюс те же визуализации — 1 кликом вставляются, ну это такое.

                  В целом — статья отличная, кроме вот этого вот нюанса ;)
                  • +3
                    Тут наверное имелось ввиду коммулятивные затраты: написание кода, подумать, трансформации, подбор параметров, построение графиков. А в Профете всего лишь 10 строк и задача по сути решена :)
                    • 0

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

          • +1
            MAPE — коварная метрика. Поскольку ошибка нормируется на таргет, при стремлении таргета к 0, ошибка может возрасти в сотни раз. Поэтому сравнительный график в начале статьи вызывает вопросы

            Допустим

            y = 60, y_hat =50: MAPE = (60 — 50)/60 = 1/6; MAE = 10

            y = 2, y_hat = 3: MAPE = (3-2)/2 = 1/2; MAE = 1

            y = 0.1, y_hat = 1: MAPE = (1-0.1)/0.1 = 9; MAE =0.9

            • 0
              А у меня
              m.plot(forecast)
              m.plot_components(forecast)
              напрочь отказываются что-либо рисовать
              Как увидеть графики? ;0((
              • 0
                впрочем, глянул на код fbPropheta разобрался
                надо было
                import matplotlib.pyplot as plt
                ...
                plt.savefig('plot.png', dpi=96)
                

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

              Самое читаемое