Prophet入門

Meta（Facebook）が公開した時系列予測のための OSS ライブラリです。インストール方法は Installation in Python を参照してください。基本的には pip install prophet を実行すればインストールできます。

1
2
3
4
5

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from prophet import Prophet

実験に使用するデータ #

1 年間（2020/1/1〜2020/12/31）のデータを用意します。曜日ごとに cos(weekday) の値が変わる 7 日周期、偶数月だけ 0.5 下がる段差、log(i + 1) / 3 でゆるやかに伸びるトレンド、rand() / 10 のノイズの 4 つを足した系列です。正解の成分が分かっているので、Prophet の分解結果と見比べられます。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

date = pd.date_range("2020-01-01", periods=365, freq="D")

# 予測対象
y = [
    np.cos(di.weekday())
    + di.month % 2 / 2
    + np.log(i + 1) / 3.0
    + np.random.rand() / 10
    for i, di in enumerate(date)
]

# トレンド成分
x = [18627 + i - 364 for i in range(365)]
trend_y = [np.log(i + 1) / 3.0 for i, di in enumerate(date)]
weekly_y = [np.cos(di.weekday()) for i, di in enumerate(date)]
seasonal_y = [di.month % 2 / 2 for i, di in enumerate(date)]
noise_y = [np.random.rand() / 10 for i in range(365)]

df = pd.DataFrame({"ds": date, "y": y})
df.index = date

# 実験に使用するデータ
plt.title("サンプルデータ")
sns.lineplot(data=df)
plt.show()

何が起きているか（ステップで見る） #

コードを読む前に、この手法が図の上で何をしているかを追ってみましょう。「再生」を押すか、スライダーを動かすと 1 手順ずつ進みます。

見るポイント：水準 → 傾き → 季節 と要素を足すほど、未来への伸ばし方が豊かになる。

時系列データを構成する要素 #

「時系列データ」という言葉にはいろいろな種類のデータが含まれますが、ここで扱うのは、タイムスタンプと数値の 2 列だけで、タイムスタンプに欠損がなく 1 日ごとの等間隔（evenly spaced）に並んでいるデータです。Prophet は欠損日があっても学習できますが、成分の図を見比べるにはまず等間隔の例で感覚をつかむのが分かりやすいので、この形に限定します。次のコードは y を作った 4 つの成分を上から順に 4 段で描き、5 段目に合成後の y を置きます。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

plt.figure(figsize=(14, 6))
plt.title("yを成分ごとに分解したもの")
plt.subplot(511)
plt.plot(x, trend_y, "-.", color="red", label="トレンド", alpha=0.9)
plt.subplot(512)
plt.plot(x, weekly_y, "-.", color="green", label="周期的な変動(週ごと)", alpha=0.9)
plt.subplot(513)
plt.plot(x, seasonal_y, "-.", color="orange", label="周期的な変動(月ごと)", alpha=0.9)
plt.subplot(514)
plt.plot(x, noise_y, "-.", color="k", label="ノイズ成分")
plt.subplot(515)
sns.lineplot(data=df)
plt.show()

Prophetで2021年の1~3月を予測する #

2020/1/1〜2020/12/31 のデータを使って、さらに先の 3 か月（2021 年 1〜3 月）を予測します。データは 1 年しかないので yearly_seasonality=False、週ごとに周期性が見られるので weekly_seasonality=True を指定し、daily_seasonality=True は元のコードに合わせてそのまま付けています。学習から予測までを train_and_forecast_pf という関数にまとめておきます。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

def train_and_forecast_pf(
    data,
    periods=90,
    yearly_seasonality=False,
    weekly_seasonality=True,
    daily_seasonality=True,
):
    """prophetを訓練して予測をする

    Args:
        data (pandas.DataFrame): 時系列データ
        periods (int, optional): 予測する期間の長さ. Defaults to 90.
        yearly_seasonality (bool, optional): 年単位の季節成分があるかどうか. Defaults to False.
        weekly_seasonality (bool, optional): 週単位の季節成分があるかどうか. Defaults to True.
        daily_seasonality (bool, optional): 日単位の季節成分があるかどうか. Defaults to True.

    Returns:
        _type_: 予測モデル, 予測結果
    """
    assert "ds" in data.columns and "y" in data.columns, "入力するデータにはds, yの列が必要です。"
    # モデルを訓練
    m = Prophet(
        yearly_seasonality=yearly_seasonality,
        weekly_seasonality=weekly_seasonality,
        daily_seasonality=daily_seasonality,
    )
    m.fit(df)

    # 将来を予測
    future = m.make_future_dataframe(periods=periods)
    forecast = m.predict(future)
    return m, forecast

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

# 予測結果を確認する
periods = 90
m, forecast = train_and_forecast_pf(
    df,
    periods=periods,
    yearly_seasonality=False,
    weekly_seasonality=True,
    daily_seasonality=True,
)
fig = m.plot(forecast)
plt.title("prophetの予測結果")
plt.axvspan(18627, 18627 + periods, color="coral", alpha=0.4, label="予測期間")
plt.legend()
plt.show()

Initial log joint probability = -32.1541
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
      99       772.276   5.98161e-05       56.7832           1           1      135   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     131        772.59    0.00128893       157.592   1.465e-05       0.001      217  LS failed, Hessian reset 
     181       772.678   3.78737e-05       49.0389   6.852e-07       0.001      326  LS failed, Hessian reset 
     199       772.681   1.42622e-06       43.2231      0.6929      0.6929      350   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     230       772.681   6.80165e-06       56.0478   7.185e-08       0.001      432  LS failed, Hessian reset 
     245       772.681   4.06967e-08       48.5475      0.1802      0.8285      454   
Optimization terminated normally: 
  Convergence detected: relative gradient magnitude is below tolerance

実行結果の読み方 #

青い予測線は学習期間では黒い点の帯のまん中を通り、7 日周期の cos 波と対数トレンドの伸びを再現しています。偶数月に 0.5 下がる段差は Prophet の季節性では表せないので、黒い点が線から交互に上下にずれる月が見えます。coral 色の予測期間（2021 年 1〜3 月）では、週次の波を保ったまま最後のトレンドの傾きで直線的に伸び、薄い青の帯が先に行くほど広がります。帯の幅がノイズ（0〜0.1）より明らかに広いのは、トレンドの変化点が未来にも起こり得ると仮定して不確実性を積み上げているためです。ログの Optimization terminated normally は Stan の最適化が収束したことを示し、LS failed, Hessian reset は途中の探索のやり直しなので結果に問題はありません。

周期性の指定の影響 #

同じデータで yearly_seasonality=True にすると何が起きるかを比べます。1 年分しかないデータに年周期を学習させる、という誤った設定の例です。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

# 予測結果を確認する
periods = 90
m, forecast = train_and_forecast_pf(
    df,
    periods=periods,
    yearly_seasonality=True,
)
fig = m.plot(forecast)
plt.title("prophetの予測結果")
plt.axvspan(18627, 18627 + periods, color="coral", alpha=0.4, label="予測期間")
plt.legend()
plt.show()

Initial log joint probability = -32.1541
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
      99       1076.54   0.000445309       68.8033           1           1      133   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     199       1078.13   0.000151685       92.7241           1           1      256   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     246       1078.14   1.78997e-06       84.0649    1.52e-08       0.001      353  LS failed, Hessian reset 
     261       1078.14   3.82403e-08       101.692      0.2973           1      372   
Optimization terminated normally: 
  Convergence detected: relative gradient magnitude is below tolerance

ログの log prob が前の 772 から 1078 へ上がっているのは、学習データへの当てはまりが良くなったことを意味しますが、1 年分のデータで年周期を入れた当てはまりの良さは過学習の証拠です。年周期を使うなら最低 2 年、できれば 3 年以上のデータを用意し、それまでは yearly_seasonality=False にしておきます。

使いどころと落とし穴 #

Prophet は ds と y を渡すだけで週次・年次の季節性とトレンドを分解してくれるので、売上やアクセス数の最初のベースラインを 1 時間で出したい場面に向きます。ただし季節性はフーリエ級数で表すため、この例の「偶数月だけ下がる」ような周期でも滑らかでもない段差は捕まえられず、add_regressor で月の偶奇を外生変数として渡すか、休日効果 として日付を指定する必要があります。データが 1 年未満のときに yearly_seasonality を有効にすると、ここで見たように 1 年分の形を丸写しした過学習になります。予測区間（yhat_lower / yhat_upper）はトレンドの不確実性だけを既定で含み、季節性の不確実性は mcmc_samples を 0 より大きくしないと入らない点も、区間を意思決定に使うときには忘れないでください。