時系列予測で使えるpythonライブラリ一覧
本記事では、時系列予測に利用できるpythonのライブラリの使い方について説明をします。
パッとライブラリを使うことを目指すため具体的なアルゴリズムの説明は省きます。
※説明が間違えている場合があればご指摘いただけると助かります。
目次
利用データ
ライブラリ
Prophet
PyFlux
Pyro
Pytorch
Lightgbm
補足:Darts
まとめ
ソースコード
このブログで記載されているソースコードはGitHubに上げておいたのでもしよろしければ参考にしてください。
github.com
利用データ
今回用いるデータはkaggleのM5 Forecasting - Accuracyと呼ばれるコンペティションで利用されたデータを用います。
作成したランダムなデータよりも実データのほうが予測をしている感があるからです。
予測に使うデータはwalmartの売上データです。
下図はそのデータをplotしたものです。
売上が0になっている日がありますが、この日付を見ると12/25になっています。
アメリカではクリスマスは休日なのでwalmartもお休みのようです。
こちらのデータの最後28個(28日間)のデータを予測します。
下の図は最後の56個のデータをplotしたものです。
曜日による周期があるように見えます。
では、こちらのデータを様々なライブラリで予測してみようと思います。
ライブラリ一覧
Prophet
Prophetはfacebookが開発している時系列解析ライブラリです。
トレンドの変化を確認したり、モデル作成時に周期やイベント効果を追加することができます。
公式ドキュメント
facebook.github.io
参考ブログ
時系列解析ライブラリProphet 公式ドキュメント翻訳１（概要＆特徴編） - Qiita
データの準備
Prophetでは、データをpandasのDataFrameで渡す必要があります。
そして、このDataFrameには2つのカラムが必要です。
2つのカラムはそれぞれ"y","ds"と名前をつける必要があります。
y : 予測したい値(数値)
ds : YYYY-MM-DD フォーマットになっている日付データ。(datestamp,strどちらでも可)
具体的には以下の様になっている必要があります。
| y | ds |
|---|---|
| 123 | 2018-1-1 |
| 111 | 2018-1-2 |
| 144 | 2018-1-3 |
| 156 | 2018-1-4 |
モデルの定義/学習
model = Prophet() # 月と週の周期を追加 model.add_seasonality(name='monthly', period=28, fourier_order=5) model.add_seasonality(name='weekly', period=7, fourier_order=5) # 学習 model.ft(data)
予測
予測をする際は、どれくらい先まで予測をするのかを決めることができます。
future_data = model.make_future_dataframe(periods=28, freq='D') forecast_data = model.predict(future_data) #予測
予測結果
予測結果を確認する方法はいくつかありますが、今回はforecast_dataに格納されているデータを用いて確認をします。
また、forecast_dataにはトレンドや周期情報なども格納されています。
plt.figure(figsize=(20,5)) plt.plot(df['y'][-56:],label='true') plt.plot(forecast_data['yhat'][-84:],label='pred') plt.legend() plt.show()
精度としてはすごく良いものではないかもしれないですが形は捉えることができています。
簡単に時系列予測をするときは便利なライブラリです。
以下のように簡単にトレンドや周期も確認できます。
model.plot_components(forecast_data)
PyFlux
※現在、PyFluxはメンテナンスを停止しており、Python3.6以降には対応していないようです。
利用の際はご注意ください。(他のライブリラを利用することを推奨します。)
PyFluxは時系列予測のためのライブラリで、確率分布を用いたアプローチをしているようです。
ドキュメント
pyflux.readthedocs.io
参考記事
時系列データの予測ライブラリ--PyFlux-- - Qiita
データの準備
PyFluxでは時系列データをそのまま扱うことができます。
なので、事前に必要な準備は特にありません。
モデルの定義
pyfluxではARIMAやARIMAXなどのモデルが使えますが、今回はARIMAモデルを使いたいと思います。
以下のように非常にシンプルに書けます。
arima_model = pf.ARIMA(data=data, ar=8, ma=8, target='sunspot.year', family=pf.Normal())
学習
学習も非常にシンプル。
summaryの確認もできます。
x = model.fit("MLE")
x.summary()
Normal ARIMA(8,0,8) ======================================================= ================================================== Dependent Variable: Series Method: MLE Start Date: 8 Log Likelihood: -18498.7318 End Date: 1940 AIC: 37033.4637 Number of observations: 1933 BIC: 37133.6666 ========================================================================================================== Latent Variable Estimate Std Error z P>|z| 95% C.I. ======================================== ========== ========== ======== ======== ========================= Constant 34480.769 3412.7463 10.1035 0.0 (27791.7863 | 41169.7517) AR(1) -0.4931 0.0673 -7.3239 0.0 (-0.625 | -0.3611) AR(2) -0.1625 0.0295 -5.5125 0.0 (-0.2203 | -0.1048) AR(3) -0.1608 0.0289 -5.5687 0.0 (-0.2174 | -0.1042) AR(4) -0.0973 0.0269 -3.6126 0.0003 (-0.1501 | -0.0445) AR(5) -0.2247 0.0335 -6.7046 0.0 (-0.2904 | -0.159) AR(6) -0.0474 0.0274 -1.7302 0.0836 (-0.1011 | 0.0063) AR(7) 0.7559 0.0363 20.8272 0.0 (0.6847 | 0.827) AR(8) 0.4058 0.0662 6.1256 0.0 (0.276 | 0.5357) MA(1) 0.9191 0.0639 14.3768 0.0 (0.7938 | 1.0444) MA(2) 0.6648 0.0479 13.8649 0.0 (0.5708 | 0.7587) MA(3) 0.6405 0.045 14.2294 0.0 (0.5523 | 0.7287) MA(4) 0.5536 0.0441 12.564 0.0 (0.4672 | 0.6399) MA(5) 0.6257 0.0455 13.762 0.0 (0.5366 | 0.7149) MA(6) 0.504 0.0364 13.8621 0.0 (0.4328 | 0.5753) MA(7) -0.1486 0.0541 -2.744 0.0061 (-0.2547 | -0.0425) MA(8) -0.1788 0.0408 -4.3821 0.0 (-0.2588 | -0.0988) Normal Scale 3467.0383 ==========================================================================================================
予測
予測結果も簡易的に確認できます。
(あまりうまく予測はできていないようですが?…)
model.plot_predict(h=28, past_values=56, figsize=(15,5))
Pyro
PyroはPyTorchをベースにした確率プログラミングのフレームワークです。
なのでpyroを利用するにはまずPyTorchをインストールする必要があります。
公式ドキュメント
pyro.ai
参考ブログ
Pyro on PyTorch の時系列モデリングが超進化していた【HMM】 - HELLO CYBERNETICS
データの準備
pyroでは予測データと共変量(特徴量みたいなイメージ)の準備をします。
今回は、日にちをone hotした値をday_of_monthに格納して入れます。
T0 = 0 # start T2 = data.size(-2) + 28 # end 28日後まで予測 time = torch.arange(T0, float(T2), device="cpu") / 365 # 共変量(特徴量みたいなイメージ) covariates = torch.cat([ time.unsqueeze(-1), day_of_month, ], dim=-1)
モデルの定義
pyro.contrib.forecastにあるForecastingModelクラスを継承して実装します。
ここで必要になるのは、model(self, zero_data, covariates) というメソッドです。
model メソッドは戻り値が不要で、内部で self.predict(noise_dist, prediction) が呼ばれていればOKです。
class Model(ForecastingModel): def model(self, zero_data, covariates): assert zero_data.size(-1) == 1 # univariate duration = zero_data.size(-2) time, feature = covariates[..., 0], covariates[..., 1:] bias = pyro.sample("bias", dist.Normal(0, 10)) # construct a linear trend; we know that the sales are increasing # through years, so a positive-support prior should be used here trend_coef = pyro.sample("trend", dist.LogNormal(-2, 1)) trend = trend_coef * time # set prior of weights of the remaining covariates weight = pyro.sample("weight", dist.Normal(0, 1).expand([feature.size(-1)]).to_event(1)) regressor = (weight * feature).sum(-1) # encode the additive monthly seasonality with pyro.plate("month",28,dim=-1): seasonal = pyro.sample("seasonal_1", dist.Normal(0, 7)) seasonal = periodic_repeat(seasonal, duration, dim=-1) # make prediction prediction = bias + trend + seasonal + regressor prediction = prediction.unsqueeze(-1) # Now, we will use heavy tail noise because the data has some outliers # (such as Christmas day) dof = pyro.sample("dof", dist.Uniform(1, 10)) noise_scale = pyro.sample("noise_scale", dist.LogNormal(-2, 1)) noise_dist = dist.StudentT(dof.unsqueeze(-1), 0, noise_scale.unsqueeze(-1)) self.predict(noise_dist, prediction)
学習
学習では、Forecasterに上記で作成したモデル,データ,共変量,ハイパーパラメータを渡すことで可能です。
戻り値は学習が終了済の予測推論ができるようになっているインスタンスになっています。
seedの設定もできるようです。
pyro.set_rng_seed(127) forecaster_options = { "learning_rate": 0.1, "learning_rate_decay": 0.1, "clip_norm": 10, "num_steps": 1001, "log_every": 100, } forecaster = Forecaster(Model(), data, covariates[:-28], **forecaster_options)
予測
forecasterにデータと共変量(予測先も含めた)とサンプルサイズを渡すことで予測が可能です。
サンプルサイズだけ予測結果を出す(quantileを取ることで予測区間を確認することができる)ため今回は平均値を予測値としました。
samples = forecaster(data, covariates, num_samples=2000).exp().squeeze(-1).cpu() pred = samples.mean(0) quantiles = [0.005, 0.025, 0.165, 0.25, 0.5, 0.75, 0.835, 0.975, 0.995] q_pred = np.quantile(samples,quantiles,axis=0)
※赤いところが0.5〜99.5%の予測区間
Pytorch(LSTM)
Pytorchは皆さんご存知ニューラルネットの作成に利用されるライブラリです。(もちろんTensorflow,kerasでもOK)
ニューラルネットにも時系列に特化したネットワークが存在します。RNNやTransformerなどもありますが、今回はLSTMのみの実装とします。
公式ドキュメント
pytorch.org
データの準備
データは少し変わった(?)用意の仕方をします。
以下のような時系列データがあった場合、
赤い線を説明変数、緑の点を目的変数として扱います。
今回は28日後までを予測させるため28日後を目的変数としました。
もちろん説明変数には時系列データ以外にも曜日等のデータを特徴量として加えることができます。
def sliding_windows(data, seq_length , train): if train: x = [] y = [] for i in range(len(data)-seq_length*2): _x = data[i:(i+seq_length)] _y = data[i+(seq_length*2)] x.append(_x) y.append(_y) return np.array(x),np.array(y)[:,0].reshape(-1,1) else: x = [] for i in range(len(data)-seq_length+1): _x = data[i:(i+seq_length)] x.append(_x) return np.array(x) seq_length = 28 x, y = sliding_windows(data, seq_length, True) trainX = Variable(torch.Tensor(np.array(x[0:-28]))) trainY = Variable(torch.Tensor(np.array(y[0:-28]))) valX = Variable(torch.Tensor(np.array(x[-28:]))) valY = Variable(torch.Tensor(np.array(y[-28:])))
モデルの定義
私自身、ニューラルネットはあまり詳しくないですが、モデルの定義方法や学習方法は一般的なニューラルネットを学習するときと変わらないです。
class LSTM(nn.Module): def __init__(self, num_classes, input_size, hidden_size, num_layers): super(LSTM, self).__init__() self.device = torch.device('cuda' if torch.cuda.is_available() else 'cpu') self.num_classes = num_classes self.num_layers = num_layers self.input_size = input_size self.hidden_size = hidden_size self.batch_size = 1 self.LSTM = nn.LSTM(input_size=input_size, hidden_size=hidden_size, num_layers=num_layers,batch_first=True,dropout = 0.25) self.fc = nn.Linear(hidden_size, num_classes) self.dropout = nn.Dropout(p=0.2) def forward(self, x): h_1 = Variable(torch.zeros(self.num_layers, x.size(0), self.hidden_size).to(self.device)) c_1 = Variable(torch.zeros(self.num_layers, x.size(0), self.hidden_size).to(self.device)) _, (hn, cn) = self.LSTM(x, (h_1, c_1)) y = hn.view(-1, self.hidden_size) final_state = hn.view(self.num_layers, x.size(0), self.hidden_size)[-1] out = self.fc(final_state) return out device = torch.device('cuda' if torch.cuda.is_available() else 'cpu') num_epochs = 200 learning_rate = 1e-3 input_size = x.shape[-1] hidden_size = 512 num_layers = 3 num_classes = 1 lstm = LSTM(num_classes, input_size, hidden_size, num_layers) lstm.to(device) criterion = torch.nn.MSELoss().to(device) optimizer = torch.optim.Adam(lstm.parameters(), lr=learning_rate,weight_decay=1e-5) scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau(optimizer, patience=100, factor =0.5 ,min_lr=1e-7, eps=1e-08)
学習
for epoch in tqdm(range(num_epochs)): lstm.train() outputs = lstm(trainX.to(device)) optimizer.zero_grad() # obtain the loss function loss = criterion(outputs, trainY.to(device)) loss.backward() optimizer.step() #Evaluate on test lstm.eval() valid = lstm(valX.to(device)) vall_loss = criterion(valid, valY.to(device)) scheduler.step(vall_loss) if epoch % 50 == 0: print("Epoch: %d, loss: %1.5f valid loss: %1.5f " %(epoch, loss.cpu().item(),vall_loss.cpu().item()))
予測
x = sliding_windows(data, seq_length, False)
testX = Variable(torch.Tensor(np.array(x)))
train_predict = lstm(testX.to(device))
ニューラルネットでも簡単(?)に時系列予測ができます。
Lightgbm(GBDT)
Lightgbmが時系列予測ライブラリ？と思った方もいるかもしれないですが、特徴量をうまく入れることで予測は可能です。
しかし、上記のライブラリのように過去のデータから時系列的に予測するということは行っていません。あくまで、過去のデータを特徴量として扱うことで予測をします。
ドキュメント
lightgbm.readthedocs.io
データの準備
以下が時系列系でよく利用する特徴量です。(今回は28日後までを予測するためlagは28を最小にしています。)
これ以外にも曜日,日付,平日,祝日,大型連休などの特徴量が効くことが多いです。
時系列予測では時系列データだけで予測しようとしがちですが使える特徴量はできるだけ使うのが無難だと思います。
for lag in [28,35,42,364]: df[f'lag_{lag}'] = df["target"].shift(lag) for lag in [28,35,42,364]: for win in [7,14]: df[f'win{win}_lag{lag}'] = df['target'].transform(lambda x: x.shift(lag).rolling(win).mean())
学習
時系列データ以外と同様の方法で学習させます。
trn_data = lgb.Dataset(train_df.drop('target',axis=1), label=train_df['target']) val_data = lgb.Dataset(valid_df.drop('target',axis=1), label=valid_df['target']) lgb_params = { 'boosting_type': 'gbdt', 'metric': 'rmse', 'objective': 'regression', 'seed': 127, 'early_stopping_rounds' : 500, 'learning_rate': 0.03 } lgb_model = lgb.train(lgb_params, trn_data, num_boost_round = 10000, valid_sets = [trn_data, val_data], verbose_eval=500)
予測
val_pred = lgb_model.predict(valid_df.drop('target',axis=1)) test_pred = lgb_model.predict(test_df.drop('target',axis=1))
こんなシンプルな特徴量だけでもそれっぽい予測ができました。さすがGBDT。
補足：Darts
上記以外のライブラリにDartsと呼ばれるライブラリがあるようです。
以下のモデルに対応しており、非常に便利そうです。
- Exponential smoothing
- ARIMA & auto-ARIMA
- Facebook Prophet
- Theta method
- FFT (Fast Fourier Transform)
- Recurrent neural networks (vanilla RNNs, GRU, and LSTM variants)
- Temporal convolutional network
詳しくは以下の記事を参考にしてみてください。
まとめ
時系列予測に使えるpythonのライブラリを紹介しました。
簡単に予測をすることができるライブラリが多くあって非常に便利でした。
特にProphetとPyFluxは非常に簡単に予測をすることができ、最初にパッと使うライブラリとして向いていると思います。
高度な予測をするには、Pyro,Pytorch,Lightgbm等を利用するのが良いかと思います。(過学習に気を付けて😨!!)
