医療職からデータサイエンティストへ

統計学、機械学習に関する記事をまとめています。

facebookの時系列予測パッケージ{prophet}を使って、ブログアクセス数を予測する。

prophetはfacebookが無料で提供している時系列予測パッケージです。RでもPythonでも使うことができます。本家様サイトによると

Prophet is a procedure for forecasting time series data based on an additive model where non-linear trends are fit with yearly, weekly, and daily seasonality, plus holiday effects. It works best with time series that have strong seasonal effects and several seasons of historical data. Prophet is robust to missing data and shifts in the trend, and typically handles outliers well.

引用:公式サイト

簡単に意訳すると

prophetは、年、週、季節、イベント、その他祝日効果とともに非線形トレンドを当てはめた加法モデルに基づいて、予測を行う。これは強い周期性をもつデータにも対応できる。さらに、欠損やトレンドの変化、そして異常値にも対応が可能である。

というものです。今回はこの預言者(Prophet)を使って、本ブログのアクセス数を予測したいと思います。prophetでは内部でstanを使っているので、場合によってはstanの設定が必要かもしれません。設定の詳細については、以下の公式ドキュメントをご参照ください。

公式ドキュメント:

Quick Start | Prophet

prophetホワイトペーパー:

Forecasting at scale [PeerJ Preprints]

Prophetの理論

Prophetは既存の時系列解析手法とは、理論的な面で少し異なっています。(以下追記参照)

一般的な時系列解析では、自己相関を考慮して、時系列データが確率過程に沿って得られていることを前提に、その確率過程をモデル化することを目的とします。状態空間モデル然り、ARIMAモデル然りです。

www.medi-08-data-06.work

www.medi-08-data-06.work

しかし、prophetが行うのは、確率過程のモデリングではなく、単に点と点をスムーズにつなぐだけです。つまり、普通の線形回帰モデルと何ら変わらない仕組みなのです。そのため、以下の記事で紹介しているような現象の解釈は行えません。

www.medi-08-data-06.work

この時間の流れを無視する代償として、prophetは周期の変動やトレンドの変化などを柔軟にモデルに組み込むことができ、その預言者の名の通り、予測を得意とするのです。

具体的には、prophetは以下のようにモデリングされます。

y(t) = g(t) + s(t) + h(t) + \epsilon (t)

見るからに線形回帰モデルみたいですね。ここでg(t)はトレンド、s(t)は周期(フーリエ変換による)、h(t)は突発的なイベント、\epsilon(t)は誤差です。詳細な理論については、論文や他のサイト様をご参照ください。m(__)m

[R] fb Prophet の解剖で学ぶベイズ時系列モデリング - ill-identified diary

Prophet統計モデル概要 | 西陣に住むデータ分析屋のブログ


2019/07/23 追記

上記のprophetと既存の時系列モデルとの関係について、以下のようにご意見頂きました。

prophetは状態空間モデルの概念に内包されるモデルというのが正しいです。既存のモデルでよく使われる自己回帰項を考慮せず、単純な線形と周期に分解することで、解釈がしやすくモデリングしやすい時系列解析モデルを作れるというのがprophet特徴です。ご指摘ありがとうございますm(__)m


今回はRでの実践法を中心に書いていきますが、pythonでも同じことができます。

ブログアクセス数を予測してみる。

さて、それでは本題のブログアクセス数を予測してみます。2018/12から記事執筆現在2019/07までのデータになります。

今回使うデータはこちらからダウンロードできます。

medi_data/prophet_data at master · kojiro0208/medi_data · GitHub

まずはデータインポートして、ざっくり傾向を見てみます。

library(tidyverse)
library(lubridate)
library(prophet)

dat <- read_csv("Analytics.csv",skip = 6,n_max = 211) %>% 
  rename("ds"="日の指標",
         "y" = "ユーザー") 

dat %>% 
  ggplot(aes(x=ds,y=y))+
  geom_line()

時系列解析

ありがたいことに右肩上がりに伸びております(T_T)注意点として、日付のカラム名はds、値のカラム名はyとしておいてください。

ここで、コレログラムを書いて自己相関を見てみましょう。

dat$y%>%
  diff() %>% 
  acf()

時系列解析

何となく7日間の周期がありそうです。本来であれば、これをtsオブジェクトに直して、stl()で周期成分を....とやるのですが、このpophetはそんなことは必要ありません!

まずは全てデフォルト設定で実行してみます。

model <- prophet(dat)
future <- make_future_dataframe(model,7*4*2)
> future %>% head()
          ds
1 2018-12-18
2 2018-12-19
3 2018-12-20
4 2018-12-21
5 2018-12-22
6 2018-12-23
> future %>% tail()
            ds
262 2019-09-05
263 2019-09-06
264 2019-09-07
265 2019-09-08
266 2019-09-09
267 2019-09-10

pred <- predict(model,future)
plot(model,pred)

時系列解析

上から見ていきます。まずはprophed()に先ほどのデータを渡して、モデルを作成します。続いてmake_future_dataframe()にモデルと予測したい日数を入れます。今回は2ヶ月先(56日)までを予測の範囲とします。すると、日数の入ったデータセットが出来上がります。あとは、predict()にモデルと予測データセットを渡せば、完了です。plot()に渡すと一発で可視化してくれます。また、トレンドや周期はprophet_plot_components()を使うと

prophet_plot_components(model,pred)

#prophet_plot_components(model,pred)[[1]]などとすると一つづつ取り出せる。

時系列解析

トレンド、周期が可視化できます。これを見るとトレンドは今のところ右肩あがりで、周期は月曜日に一番アクセス数が多く、土曜日が一番少ないことが分かります。確かに、休日はアクセス数が少ない印象でした...

さて、もう一度プロットを見てみるとあまり当てはまりがよくないように見えます。ここで、もう一つ設定を加えてみます。

model <- prophet(dat,seasonality.mode = "multiplicative")
pred <- predict(model,future)
plot(model,pred)
prophet_plot_components(model,pred)

時系列解析時系列解析

先ほどより当てあまりが良いよくなりました。何を加えたかというと、prophet()の引数に’seasonality.mode = "multiplicative" ’を加え、デフォルト設定が”additive”(加法的)であったのを、周期が"multiplicative"、つまり乗法的に増えるということをモデルに組み込みました。これは、周期成分が以下のような感じに増えそうな場合に用いることが出来ます。

時系列解析

ブログのアクセス数の場合、増えれば増えるほど、分散は大きくなるので、乗法モデルが使えそうです。

交差検証(クロスバリデーション)による時系列モデルの評価

モデルを作ったら、予測精度の評価を行います。prophetは交差検証法によるモデル評価も簡単に行えます。時系列データの予測精度は、通常と少し変わっていて、1日先、2日先、3日先...とどこまで先の日数までなら精度よく予測できるかを調べます。こうすることで、1週間ぐらいまでなら大丈夫そうだけど、10日先以上の予測をするときは気をつけてねなどということができます。

また、交差検証する際は、トレーニングするデータをずらしながら調べるrolling originという手法を使います。prophetでは、予測を開始する時期(initial)、予測する期間(horizon)、データをずらす期間(period)を指定することができます。例えばinitial120日、horizon30日、period15日とするとこんな感じです。

時系列解析

1個目のモデルは、データから120日目から121~150日後までの予測値を評価します。トレーニング期間はhorizonの3倍なっているらしいので、30~120日までがモデル作成のトレーニング期間となります。そして、2個目のモデルでは、15日づつずらし、135~160日後を予測値の評価期間として、45~135日目までがトレーニング期間となります。これをデータの終わりまで続けます。

ここで見てわかるように、トレーニング期間は学習期間の3倍必要なので、Initialが短すぎると、最初の方のモデルは精度が悪く、periodで進むごとに使えるデータが増え、精度が上がっていくという複雑なことが起こるので、initialは十分に余裕を持って決めましょう。また、peridoは指定しないとデフォルトでhorizonの半分になります。

実際にやってみましょう。先ほどの加法モデル、乗法モデルを比較してみます。

#加法モデル
model_add <- prophet(dat)
df_cv_add <- cross_validation(model_add,initial = 120,horizon = 56,period = 10, units = 'days')
plot_cross_validation_metric(df_cv_add ,metric = "mape")+
  geom_hline(yintercept = 0.5,color="red")

#乗法モデル
model_multi <- prophet(dat,seasonality.mode = "multiplicative")
df_cv_multi <- cross_validation(model_multi,initial = 120,horizon = 56,period = 10, units = 'days')
plot_cross_validation_metric(df_cv_multi ,metric = "mape")+
  geom_hline(yintercept = 0.5,color="red")

時系列解析

時系列解析

モデルの作成は先ほどと同様です。次のcross_validationにはモデルとinitial、horizon、periodを指定します。今回はhorizonを2ヶ月先の56日としたので、本来はinitialが160日ぐらい欲しいところですが、データセット211日分と少ないため、今回は仕方なく120日としています。そして、periodを10日に設定したので、4つのモデルが出来ます。そして、次の plot_cross_validation_metricで結果を可視化できます。評価指標にはMSE(平均二乗誤差)・RMSE(MSEの平方根をとったもの)・MAE(平均絶対誤差)・MAPE(平均絶対パーセント誤差)などがあり、今回はよく使われるMAPEを使いました。ちなみにprophetの可視化では全てにおいてggplotライクに記述できます。

結果を見てみると、先の日数になるにつれて少しだけ予測精度が悪くなっていますが、どちらもおよそ50%の範囲の誤差があるようです。結果を取り出すには以下のようにします。

performance_metrics(df_cv_add )
> df_p_add  %>% head()
  horizon       mse     rmse      mae      mape  coverage
1  6 days  568.6869 23.84716 20.13109 0.3465458 0.8863636
2  7 days  574.5968 23.97075 20.04307 0.3400142 0.9090909
3  8 days  594.6151 24.38473 20.07582 0.3265833 0.9090909
4  9 days  743.7558 27.27189 22.25053 0.3359939 0.8636364
5 10 days 1148.1578 33.88448 27.48582 0.3757320 0.7727273
6 11 days 1155.8110 33.99722 27.83652 0.3838891 0.7727273

時系列データの交差検証は奥が深いので、もし間違ったことを書いていれば、ご指摘ください。

触れなかった機能

上記の記事では触れていませんでしたが、少し重要な機能があります。それは、prophetの変化点検知です。上のトレンドグラフを見るとトレンドが変化しているところがあるかと思います。実はこのトレンド変化点の検知がprophetを使う意義の一つなのです。

prophetは内部ではstanを使って、ベイズ推定を行なっています。そして、変化点の数を適当に決めて、ラプラス分布という正規分布の先の尖ったバージョンを事前分布として、スパース推定のようなことを行なっています。詳細については以下が分かりやすいです。

[R] fb Prophet の解剖で学ぶベイズ時系列モデリング - ill-identified diary

そして、prophetでは、変化点の数やパラメーターを指定することでトレンド変化の調整をすることができます。何がともあれやってみましょう。

model <- prophet(dat,changepoint.prior.scale = 0.1,n.changepoints = 25)
pred <- predict(model,future)

plot(model,pred)+
  add_changepoints_to_plot(model)

#おまけ
tibble(points=model$changepoints,
       change = model$params$delta[1,]) %>%
  ggplot()+
  geom_bar(aes(x=points,y=change),stat = "identity")

時系列解析時系列解析

初めのprophet()にchangepoint.prior.scale = 0.1を指定して、トレンドの柔軟性を調整しています。デフォルトは0.05になっていて、これを大きくすることで、トレンドが柔軟になります。そして、n.changepointsで変化が予測される数を指定できます。こちらのデフォルトは25です。changepoint.prior.scale の方を調整してモデルを評価するのが一般的みたいです。そして、plot()add_changepoints_to_plot()を加えることで、変換点を可視化できます。便利ですねー。おまけとして、それぞれの変化点でどんなトレンド変化があったのかも可視化しました。

また、線形トレンドだけでなく、ロジスティック関数的なトレンドを指定することもできます。これは、ある一定までいくと頭打ちになることが予測されるデータに用います。おそらくそんなことはありませんが、200アクセスぐらいで頭打ちになると仮定してみましょう。

dat$cap <- 200
model <- prophet(dat,growth = "logistic")
future$cap <- 200
pred <- predict(model,future)
plot(model,pred)+
  add_changepoints_to_plot(model)

時系列解析

prophet()にgrowth="logistic"を指定します。注意点としては、モデルを作成するデータ、予測をするデータそれぞれにcapという変数を作り、そこに予測される頭打ち数を加えます。しかし、このモデルでは頭打ち数に依存するため、明らかに上限が決まっている場合にのみ有効ですね。

イベント効果を使って予測精度をあげてみる。

さて、今までは曜日の周期とトレンドのみでモデルを作りました。prophetの良いところはドメイン知識に基づき、外部変数やその他のイベント効果を”簡単に”モデルに加えることができます。アクセス数がお休みに影響されるのであれば、祝日の効果は入れておいたほうが良さそうです。また、感覚的にはブログの投稿日から2日間ぐらいはアクセス数が伸びている感覚です。この2つをモデルに加えてみます。あらかじめデータセットを作っておきました。

article_dat <- read_csv("article_dat.csv")

> article_dat %>% head()
holiday ds upper_window lower_window
article_post 2018-12-18 2 0
article_post 2018-12-19 2 0
article_post 2018-12-22 2 0
article_post 2018-12-23 2 0
article_post 2018-12-23 2 0
article_post 2018-12-24 2 0
syukujitsu_dat <- read_csv("syukujitsu_dat.csv") 
syukujitsu_dat %>% head() 
holiday ds upper_window lower_window
syukujitsu 2018-12-24 0 0
syukujitsu 2019-01-01 0 0
syukujitsu 2019-01-14 0 0
syukujitsu 2019-02-11 0 0
syukujitsu 2019-03-21 0 0
syukujitsu 2019-04-29 0 0

モデルに加える際には、変数名はこれと同じにしてください。holidayには、任意のイベント名、dsにはそのイベントの日付、upper_window、 lower_windowにはそのイベントの効果が前後何日に影響を及ぼすかと指定できます。今回のブログ投稿日は投稿から2日ぐらいはアクセス数伸びる印象なので、2としました。祝日はおそらくその日しか影響しないので、前後とも0にしてあります。

注意点として、このデータセットは予測したい未来の日付まで含める必要があるので、あらかじめそこでイベント発生が予測されるものしか入れることができませんん。祝日は未来まで決まっていますが、ブログ投稿日は私次第(笑)なので、データのない未来の日付には毎週日曜日に投稿することを前提しています。

さて、モデルに加えてみます。

#データフレーム を一つにする。
holiday_dat <- bind_rows(syukujitsu_dat,article_dat)

model <- prophet(dat,holidays=holiday_dat,seasonality.mode = "multiplicative")
pred <- predict(model,future)
plot(model,pred)+
  add_changepoints_to_plot(model)
prophet_plot_components(model,pred)

時系列解析時系列解析

簡単ですね!prophet()のholidayにそれぞれのイベントデータフレームを結合したものを渡すだけです。結果を見てみると、突発的な変動にうまく対応しているようにも見えます。次に曜日効果をみてみると、土曜と日曜のマイナス効果がほぼ同じになりました。これは日曜日にブログ投稿をすることが多かったので、その効果がイベントに吸収され、本来のお休み効果になったことが考えられますね。

また、トレンドと周期に加えて、holidayというイベント効果が追加されました。しかし、これは祝日効果とブログ投稿効果を足しわせた効果なので、わかりやくするためにそれぞれの効果を分けて書いてみます。

#ブログ投稿と祝日の効果を分けて表示
pred %>% 
  select(ds,article_post,syukujitsu) %>% 
  gather(key = "holiday",value = "y",-ds) %>% 
  ggplot()+
  geom_line(aes(x=ds,y=y,color=holiday))

時系列解析

予想通り、ブログ投稿は変動をプラス側に、祝日は変動をマイナス側に動かす要因となるようです。このモデルで予測精度を評価した結果がこちらです。

時系列解析

あら、イマイチのようですね....特に20日を越えた辺りから精度が落ちています。

外部変数を使ってみる。

いや、まだまだ諦めません。今度はブログの記事数をモデルに入れてみます。この考え方を応用すれば、他の時系列データを入れることができます。ブログの記事数のデータを作っておきました。

dat_with_articles <- read_csv("dat_with_articles.csv")

> dat_with_articles %>% head() 
ds y articles
2018-12-18 0 1
2018-12-19 0 2
2018-12-20 0 2
2018-12-21 0 2
2018-12-22 0 3
2018-12-23 0 5
future_with_articles <- read_csv("future_with_articles.csv")
> future_with_articles %>% head() 
ds articles
2018-12-18 1
2018-12-19 2
2018-12-20 2
2018-12-21 2
2018-12-22 3
2018-12-23 5

外部変数は、元のデータセット、そしてmake_future_dataframe()で作ったデータセットに新たな変数として加えます。今回も未来の記事数は週に1度更新することを前提とし、articlesとして変数に加えました。あとはこれを、モデルに加えるだけです。

model <- prophet(holidays=syukujitsu_dat,seasonality.mode = "multiplicative") %>% 
  add_regressor(name="articles") %>% 
  fit.prophet(dat_with_articles)
pred <- predict(model,future_with_articles)

plot(model,pred)+
  add_changepoints_to_plot(model)
prophet_plot_components(model,pred)

時系列解析時系列解析

holidayには祝日効果のみ加えました。そして、ポイントは、add_regressor()fit.prophet()です。流れとしては、prophet()でモデルの型を決める、add_reggressorで加えたい変数名を指定する、fit.prophet()でデータセットを渡す、です。グラフには、トレンド、周期、イベント効果に加えて、外部変変数効果も加わりました。少し見にくいですが...

さて、このモデルの予測精度を評価した結果がこちらです。

時系列解析

おお!さっきより改善しました。このモデルは使えそうです。トレンドと周期だけのモデルの精度が良いのは、ブログを始めた初期はアクセスがほとんどなく、イベントに影響されなかったためだと思われます。今回の結果からも記事数や祝日は影響を与えることが予想されるので、これを最終モデルとしましょう!

予測によれば、2ヶ月後には250アクセス近く行くようです...頑張ろう

まとめ

prophetは統計の専門知識よりも、ドメイン知識を使ってイベントや外部変数を柔軟にモデルが作れる印象です。データサイエンススキルが一般化する流れの中では、このprophetは大いに活躍するのではないでしょうか?

今回は紹介できませんでしが、周期変動を自分で指定し時刻周期を扱ったり、サンプリングをmcmcで行ったり、とても柔軟に対応できます。また、2ヶ月後に予測結果の確認と共に紹介するかもしれません(~~)

追記:2020/1/8

やってみました。

www.medi-08-data-06.work

※本記事は筆者が個人的に学んだことをまとめた記事なります。所属する組織の意見・見解とは無関係です。また、数学の記法や詳細な理論、用語等で誤りがあった際はご指摘頂けると幸いです。

参考

時系列解析ライブラリProphet 公式ドキュメント翻訳1(概要&特徴編) - Qiita

[R] fb Prophet の解剖で学ぶベイズ時系列モデリング - ill-identified diary

Facebookによる時系列予測AI Prophet紹介 | 西陣に住むデータ分析屋のブログ