PythonでProphetを使いたい(3:デフォルトの周期性)

時系列=トレンド + 周期性 + その他のうちの周期性、さらにその中のデフォルトの周期性の確認とチューニングをします。

▼こちらで作ったデータ、モデルをさらにチューニングしていきます。 esu-ko.hatenablog.com

デフォルトの周期性の確認

デフォルトでは年ごとの影響、曜日の影響、時間の影響を設定できます。元のモデルは曜日影響だけ載せたので、これの中身を見ていきます。

曜日周期はweeklyの中に格納されています。

pred['weekly'][0:6].plot()
df['season'][0:6].plot()

比較的、形としてはとらえられているようです。 続けて、推定された曜日ごとの数値と、実測値に占める割合を見てみます。

week = pred[['ds','weekly','weekly_lower','weekly_upper']]
week.assign(dow=(
    week['ds'].dt.dayofweek.apply(lambda x : str(x)+"_")
    + week['ds'].dt.day_name()
    )
).groupby('dow').mean().sort_index(ascending=False).plot.barh(y='weekly')

(pred['weekly']/ts).plot()

月曜日は+の影響目、土日にはマイナスの影響をうけています。

また、noiseは置いておいて、真ん中ごろ(外部影響で上がったあと）から徐々に曜日影響の占める割合が減ってきており、トレンドで増加した影響が強くなってきているのがわかります。

周期性にもトレンドの影響を乗せる

この影響を乗せるために、seasonality_mode='multiplicative'を足します。 こうすることで、観測値 = トレンド + トレンド * 周期性 + その他のモデルになります

mdl_df['cap'] = 6
mdl = Prophet(
    yearly_seasonality=False,
    weekly_seasonality=True,
    daily_seasonality=False,
    growth = 'logistic',
    changepoint_prior_scale=0.5,
    #seasonalityを増加させる
    seasonality_mode='multiplicative'
    
)
mdl.fit(
    mdl_df
)
future = mdl.make_future_dataframe(periods=0,freq='D')
future['cap'] = 6
pred  = mdl.predict(future)

可視化

ts.plot()
pred.yhat.plot()

(pred['weekly'] * pred['trend']).plot()

曜日影響にトレンドの影響をのせた値を推定しました。

PythonでProphetを使いたい(2:トレンドと転換)

時系列はトレンド + 周期性 + その他に分けて分析します。

今回はProphetのモデルのトレンドに仮定を追加したり、変化点を取り出すことで上記の式のトレンド部分の予測の調整/分析します。

▼こちらの記事で、Prophetのデフォルト予測まで実施しました。 esu-ko.hatenablog.com

Prophetのトレンドはデフォルトでは線形という仮定を置いており、予測の形は下記のようになっています。

しかし、この切りあがり停滞する形の外部要因にあたりがついており、このような線形の上昇がない、という前提をいれたい場合があります。

予測の調査

このような前提をいれたい場合、トレンドの上限にcapをつけます。
付け方は、データにcapというカラムを追加するのと、mdlにgrowth = 'logistic'を設定します。

今回はデータをもとに6くらいで停滞する、という判断がされた、という前提にします。

mdl_df['cap'] = 6
#mdl_df['floor'] = 0
mdl = Prophet(
    yearly_seasonality=False,
    weekly_seasonality=True,
    daily_seasonality=False,
    growth = 'logistic' #非線形のトレンドを仮定
)
mdl.fit(
    mdl_df
)
future = mdl.make_future_dataframe(periods=365,freq='D')
future['cap'] = 6
#future['floor'] = 0
pred  = mdl.predict(future)

※加工する場合はfloorという前提をつけます。

可視化すると下記のようになり、上昇後は横ばいになっていることがうかがえます。

pred.yhat.plot()
ts.plot()

トレンド転換の仕組み

Prophetのトレンド転換は、転換となる点をばらまいて、その転換の強さを推定しています。

import numpy as np

#転換点
cp = mdl.changepoints

#delta:転換の強さ
deltas = mdl.params['delta'].mean(0)

#changepoint_prior_scale:フィットさせる強さ
mdl = Prophet(
    yearly_seasonality=False,
    weekly_seasonality=True,
    daily_seasonality=False,
    growth = 'logistic',
    changepoint_prior_scale=0.5 #これ
)

時系列データを見るときに転換点はいつであった可能性が高いのか、を定量的に把握したいことがあり、それを抽出できると非常に便利です。

# 転換点と転換可能性のdf
cp_df = pd.DataFrame(
    {
        "ds":cp,
        "deltas":deltas
        
    }
)
cp_df.set_index('ds').plot.barh()

これをみると6月ごろに大きな転換が、前半にちょいちょいトレンドが動いていたことがわかります。

トレンド動きと、転換点の強さを実測データに重ねてみると、議論をするときに便利です。

#上位1/4のみ可視化する
#75パーセンタイルを閾値
th = cp_df['deltas'].quantile(0.75)

#抽出して予測データと結合
f = lambda x : x if abs(x) >= th else np.nan
w_cp_df = pred[['ds','trend','yhat']].merge(cp_df.assign(deltas = cp_df.deltas.apply(f)),how='left')

# 以下可視化
ts.plot()
plt.plot(w_cp_df['trend'])

#転換位置にvlinesで線を引く
y_min_v = 0 #線の下限
y_max_v = 10 #線の上限
plt.vlines(w_cp_df[w_cp_df['deltas'].notna()].index,ymin=y_min_v,ymax=y_max_v,linestyle="dotted",color='green')

トレンド転換の位置の日付がわかったので、これをベースに施策/外部変化がなかったといった定性調査ができそうです。

ProphetをPythonで使いたい(1:基本編)

Prophetは時系列モデルを簡単に扱える手法です。 Facebookから発表され、Pythonからも使用することができます。

モデリングでいじれる部分がかなり多いので、何本かの記事に分けて試していきます。

使い方

sklearn-likeな呼び出しになっています。
predict前の準備のみ注意が必要です。

from fbprophet import Prophet

mdl = Prophet()
#データをfit
mdl.fit(mdl_df)

#predict
future = mdl.make_future_dataframe(periods=0,freq='D')
pred  = mdl.predict(future)

使ってみる

テストデータをつくる

テストデータをつくって試してみます。 今回は

• 祝日
• 曜日周期
• 外部効果でベースがあがる
• よくわからないトレンド
• 別の変量効果

といろいろな要素を盛り込んでみます。

import math
import random
import pandas as pd
import jpholiday

N  = 365

#祝日効果
dt = pd.Series(pd.date_range(start='2019-01-01',periods=N))
f = lambda x : random.random() * 5 if jpholiday.is_holiday_name(x) else 0
hol = dt.apply(f)

#外部効果で突然上がる効果
trend1 = [ 2/(1 + math.exp(-1*(i - N/2) * 0.1)) for i in  range(N)]

#変な周期性効果
trend2 = [ math.cos(i * 0.25) for i in range(N)]

#曜日効果
season_base = [ random.random() for i in range(7)]
season = [season_base[i%7] for i in range(N)]

#別の変数
x = pd.Series([random.random() * 0.01 for i  in range(N)]).cumsum()
y = x.apply(lambda x : x * 1.25 + random.gauss(mu=0.5,sigma=0.25))

#ノイズ
noise = [random.random() for i in range(N)]

df = pd.DataFrame(
    {
        "trend1":trend1,
        "trend2":trend2,
        "season":season,
        "noise":noise,
        "hol":hol,
        "y":y
    }
)


#上記の要素を足し合わせる

ts = df.apply(lambda x : sum(x),axis=1)
ts.plot()

データの読み込みとモデリング

ほとんどいじらずにフィッティングしてみます。

つっこむデータはds、yというカラム名で入れる必要があります。

mdl_df = pd.DataFrame(
    {
        "ds":pd.date_range(start='2019-01-01',periods=365),
        'y':ts
        
    }
    
)

フィッティングしていきます。 一年分のデータかつ,時間情報はないので、weekly以外の季節性をオフにしています。(勝手にオフにはしてくれるのですが、わかりやすくするため)

from fbprophet import Prophet

mdl = Prophet(
    yearly_seasonality=False,
    weekly_seasonality=True,
    daily_seasonality=False
)
mdl.fit(
    mdl_df
)
future = mdl.make_future_dataframe(periods=0,freq='D')
pred  = mdl.predict(future)

pred.head()
pred.yhat.plot()
ts.plot()

なんとなくこれでも、緩やかな形はとらえてくれてはいます。実務上、とてもざっくりみるだけならこれでもよいこともあるかもしれません。

ただし、外部効果としていれてロジスティックの形などはほとんどとれていません。時系列の分析としては、まだまだ仮定を入れていく必要がありそうです。

予測してみる

periodsの期間を長くとるだけです。

future = mdl.make_future_dataframe(periods=365,freq='D')
pred  = mdl.predict(future)
pred.yhat.plot()
ts.plot()