備忘録的な

プログラミングや機械学習に関する備忘録

時系列予測ツールProphetの使い方備忘録

Facebookが提供している時系列予測ツールのOSSProphetの使い方を備忘録的に残しておきます*1
動作確認したProphetのバージョンは1.0です.

基本的な使い方

library(prophet)
# 必須カラムはds(日付)とy(目的変数)
# ds列のフォーマットはYYYY-MM-DDかYYYY-MM-DD HH:MM:SS
df <- read.csv('../examples/example_wp_log_peyton_manning.csv') 
# モデルのfit
m <- prophet(df)
# 予測結果格納用データフレームの作成
# 学習期間も含むので注意
future <- make_future_dataframe(m, periods = 365)
# 予測
forecast <- predict(m, future)
# 結果の描画
plot(m, forecast)
# 各成分の描画
prophet_plot_components(m, forecast)

飽和するトレンドのモデリング

Prophetはロジスティック曲線を使って飽和するトレンドをモデリングできます.
トレンドの上限はcapを,下限はfloorをdfの新たな列として加えます.これらの値は全行で同じ値である必要はありません.

df$cap <- 8.5
df$floor <- 1.5
m <- prophet(df, growth = 'logistic')
future <- make_future_dataframe(m, periods = 1826)
future$cap <- 8.5
future$floor <- 1.5
fcst <- predict(m, future)

トレンドの変化点のモデリング

Prophetはトレンドの変化点の候補を多めに(初期値は25)見積もり,L1正則化のような形で変化点の絞り込みを行います.変化点の候補はデフォルトでは時系列データの先頭から80%の時点から抽出されます(予測が最後の方の変化に引きずられないようにするため).
変化点候補の数はn.changepointsパラメータで,変化点抽出レンジはchangepoint.rangeパラメタ(0~1)で調整できます.
抽出された変化点は以下のようにして可視化できます.

plot(m, forecast) + add_changepoints_to_plot(m)

トレンド変化への追従度合いはchangepoint.prior.scaleパラメータで変更できます.初期値は0.05で大きくすると追従しやすくなります.このパラメータは後述する交差検証でチューニングすることができます.

変化点は以下のように明示的にしていすることもできます.

m <- prophet(df, changepoints = c('2014-01-01'))

季節成分のモデリング

ここでいう季節とは四季のことではなく,年次や週次の周期性のことです.
Prophetは年次,週次,日次の周期性モデリングをデフォルトでサポートしており,以下のように書くことができます.

m <- prophet(df, weekly.seasonality = TRUE, yearly.seasonality = TRUE, daily.seasonality=FALSE)

季節成分はどのくらい高周波成分に追従するかというハイパーパラメータを持っています.これを調整する場合は,bool値ではなく数値を指定します.yearly.seasonalityの初期値は10,weekly.seasonalityの初期値は3です(daily.seasonalityの初期値は分かりませんでした).
季節成分は以下のように可視化することができます.

prophet:::plot_yearly(m)

月次の周期性のように新たな季節成分を取り入れたいときは以下のようにします.fourier.orderは先の高周波成分への追従性に対するハイパーパラメータです.

m <- prophet(weekly.seasonality=FALSE, yearly.seasonality=FALSE)
m <- add_seasonality(m, name='monthly', period=30.5, fourier.order=5)
m <- fit.prophet(m, df)

条件付き季節性

Prophetのチュートリアルで使われる,Wikipediaペイトン・マニング(アメフトの選手)の検索数において,週次成分がオンシーズンとオフシーズンで異なるという仮説を立てた場合,以下の用の考慮することができます.

m <- prophet(weekly.seasonality=FALSE)
m <- add_seasonality(m, name='weekly_on_season', period=7, fourier.order=3, condition.name='on_season')
m <- add_seasonality(m, name='weekly_off_season', period=7, fourier.order=3, condition.name='off_season')
m <- fit.prophet(m, df)

デフォルトのweekly.seasonalityをオフにし,新たに独自の週次成分を2つ追加します.
こららの成分はdfに新たに追加されたbool値の変数,on_seaon,off_seaonがTRUEのときのみ有効になります.

祝日やイベントのモデリング

以下のようなcsvファイルを用意します.
必須な列はヘッダ名がholidayの列とdsの列です.
holiday列にイベント名を,ds列にそのイベントの発生日を入れます.
このデータには,予測期間の情報も含まれている必要があります.

lower_window,upper_window列は,そのイベントの前後への影響日数を記載します.
例えばクリスマスというholidayのイベント(2021-12-25)に対し,クリスマスイブも考慮したい場合は,
lower_windowを-1とします.

この表の例では,superbowlsの日付はplayoffの日付と重複しています.
このような場合,playoffの効果にsuperbowlsの効果が加算されます.

holiday ds lower_window upper_window
playoff 2008-01-13 0 1
playoff 2009-01-03 0 1
playoff 2010-01-16 0 1
playoff 2010-01-24 0 1
playoff 2010-02-27 0 1
playoff 2011-01-08 0 1
playoff 2013-01-12 0 1
playoff 2014-01-12 0 1
playoff 2014-01-19 0 1
playoff 2014-02-02 0 1
playoff 2015-01-11 0 1
playoff 2016-01-17 0 1
playoff 2016-01-24 0 1
playoff 2016-02-07 0 1
superbowls 2010-02-07 0 1
superbowls 2014-02-02 0 1
superbowls 2016-02-07 0 1

このデータを以下のように読み込みます.

df <- read.csv("example_wp_log_peyton_manning.csv")
holidays <- read.csv("holidays.csv")
m <- prophet(df, holidays = holidays)

季節成分,イベントの影響度調整

季節成分やイベントの影響がややデータに対し過剰適応気味だと判断した場合,seasonality_prior_scale,holidays_prior_scaleを小さくすることで対応できます.初期値は10です.

m <- prophet(df, holidays = holidays, holidays.prior.scale = 0.05)

イベントごとに影響度を調整したい場合,データフレームholidaysにprior_scaleの列を追加します.
年次や週次の成分ごとに影響度を調整したい場合,add_seasonalityのハイパーパラメータ,prior.scaleで調整します.

m <- add_seasonality(m, name='weekly', period=7, fourier.order=3, prior.scale=0.1)

回帰項の追加

いわゆる外生変数は以下のように追加できる(元のdfに新たな列nlf_sundayが追加されているものとする).

m.add_regressor('nfl_sunday')
m.fit(df)

外生変数が0/1の場合はイベントとして扱うこともできる.0/1以外の変数も加えることができる(アイスクリームの売り上げ予測に気温を使うなど)が,その値の未来の値が手に入るかをよく考える必要がある.
外生変数の係数はprophet::regressor_coefficientsで確認することができる.

季節成分(やイベント成分)の乗算

デフォルトでは季節成分やイベント成分はトレンド成分に加算されますが,例えばトレンドの上昇に伴い季節成分の増減も大きくなっていくような場合,以下のようにすると季節成分の(デフォルトではイベント成分や外生成分も同様に)影響を乗算することができます.

m <- prophet(df, seasonality.mode = 'multiplicative')

ある成分は乗算,別の成分は加算したい場合は以下のようにします(が本当にそのようなモデリングが妥当かはよく考える必要があります).

m <- prophet(seasonality.mode = 'multiplicative')
m <- add_seasonality(m, 'quarterly', period = 91.25, fourier.order = 8, mode = 'additive')
m <- add_regressor(m, 'regressor', mode = 'additive')

不確定区間

トレンドの不確定区間*2(デフォルトは0.8)は以下のように変更することができます.

m <- prophet(df, interval.width = 0.95)

季節成分の不確定区間もみたいばあい,以下のようにMAP推定ではなくMCMCサンプリングを使う必要があります.

m <- prophet(df, mcmc.samples = 300)

生の事後予測サンプルはpredictive_samples(m, future)で得ることができます.ちなみに,WindowsPythonのPyStanは非常に遅いので,MCMCサンプリングをしたい場合は,Rを使うか,LinuxPythonを使うのがよいそうです.

外れ値の扱い

外れ値があると予測値の不確定区間が非常に大きくなってしまいます.Prophetは欠損値を扱うことができるので,外れ値の影響を軽減したければ,外れ値をNAにすることで対応できます.

outliers <- (as.Date(df$ds) > as.Date('2010-01-01')
             & as.Date(df$ds) < as.Date('2011-01-01'))
df$y[outliers] = NA
m <- prophet(df)

時間軸や月次の予測

時間軸の予測をするさいに,もしデータが午前中のデータしか含まないような場合,futureデータフレームもどうようの形式に合わせてあげればよい.

future2 <- future %>% 
  filter(as.numeric(format(ds, "%H")) < 6)
fcst <- predict(m, future2)

月次のデータを扱う場合,dtをYYYY-MM-01のような形で記載すればよい.このとき,予測を月次で行うには以下のように記載する.

future <- make_future_dataframe(m, periods = 120, freq = 'month')

交差検証によるハイパーパラメータチューニング

ハイパーパラメータをチューニングするために交差検証をすることができます.

df.cv <- cross_validation(m, initial = 730, period = 180, horizon = 365, units = 'days')

initialが訓練期間,horizonが検証期間,periodがスライド幅です.ペイトン・マニングのサンプルデータ(2007/12/10~2016/1-20)に対して上記のコードを実行した場合,Cutoff date(検証期間の初日)は,
'10/2/15,'10/8/14,'11/2/10,'11/8/9,'12/2/5,'12/8/3,
'13/1/30,'13/7/29,'14/1/25,'14/7/24,'15/1/20
の11か所になります('10/2/15から180日遡るとデータ始点まで618日となり730日が確保できないため).

Rでは,unitsはweeksより短い必要があるようです.

以下のように独自のカットオフ日を設定することもできます.

cutoffs <- as.Date(c('2013-02-15', '2013-08-15', '2014-02-15'))
df.cv2 <- cross_validation(m, cutoffs = cutoffs, horizon = 365, units = 'days')

交差検証の結果の要約統計量はperformance_metricsで,プロットはplot_cross_validation_metricで確認することができます.

df.p <- performance_metrics(df.cv)
plot_cross_validation_metric(df.cv, metric = 'mape')

Pythonですと交差検証の並列実行ができるようですが,Rには用意されていないようです.

Prophetのドキュメントで推奨されているチューニングすべきハイパーパラメータは,changepoint,seasonality,holidaysのprior scaleと,seasonality_mode(additive or multiplicative)です.

年次,週次,日次の周期性を見るかどうかのオプションはデータの取得期間で決めればよいと記載されています(複数年のデータであっても年次の周期性が明確にみられなければFALSEにしてもよい気がしますが).

モデルのセーブとロード

saveRDS(m, file="model.RDS")  # Save model
m <- readRDS(file="model.RDS")  # Load model

*1:結局ほとんどドキュメントの和訳になってしまいました

*2:信頼区間(Confidence interval)ではなく不確定区間(Unsertainty interval)と使おうという提言があるみたいです.