15
Help us understand the problem. What are the problem?

More than 1 year has passed since last update.

posted at

Dartsを使ってみた

  • 製造業出身のデータサイエンティストがお送りする記事
  • 今回は時系列モデルを統一的なAPIで使用できるライブラリー(Darts)があったので使ってみました。

はじめに

過去に時系列解析に関する記事を整理しておりますので、興味ある方はそちらも参照して頂けますと幸いです。

Dartsとは

時系列モデルを分析する際は、様々なライブラリのAPIを理解して、それに合わせてデータ整形し、モデルを評価する必要がありました。
しかし、スイスの企業Unit8が2020年6月末に公開したDartsは上記のような課題を解決するライブラリです。簡単に言うと時系列データに関する様々なモデルをscikit-learnのように統一的に扱うことができます。

現在は下記モデルが対応しているようです。

  • 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
  • Transformer
  • N-BEATS

Dartsを使ってみた

ライブラリーのインストールは下記で簡単にできます。

pip install u8darts

まずは「AirPassengers」のデータを使って分析を行っていきます。

# 必要なライブラリーのインポート
import darts
from darts import TimeSeries
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

from darts.metrics import mape, mase
from darts.utils.statistics import check_seasonality, plot_acf, plot_residuals_analysis

import warnings
warnings.filterwarnings("ignore")
import logging
logging.disable(logging.CRITICAL)

# https://stat.ethz.ch/R-manual/R-devel/library/datasets/html/AirPassengers.html
df = pd.read_csv('AirPassengers.csv')
# データの中身を確認
df.head()

スクリーンショット 2021-01-14 13.37.12.png

次にDartsのTimeSeries型に変換します。

# DartsのTimeSeriesに変換
series = TimeSeries.from_dataframe(df, time_col='Month', value_cols='#Passengers')

ExponentialSmoothing(指数平滑法)

初めに、HyndmanとAthanasopoulosによる指数平滑法を試してみます。今回は各手法の詳細は省きます。

# Jan-59以前と以後に分割
train, val = series.split_after(pd.Timestamp('Jan-59'))
from darts.models import ExponentialSmoothing

# モデル生成
model_es = ExponentialSmoothing()
# 学習
model_es.fit(train) 
# 予測 (predictには予測数を入れることに注意)
historical_fcast_es = model_es.predict(len(val))
# 可視化
series.plot(label='actual', lw=1)
historical_fcast_es.plot(label='forecast', lw=1)
plt.legend()
plt.xlabel('Year')

スクリーンショット 2021-01-14 13.43.04.png

NaiveSeasonal

次にNaiveSeasonalと言う手法を試してみます。パラメータKが重要らしいので、適切に設定しないと上手く予測できないようです。

# 自己相関を求める
plot_acf(train, m = 12, alpha = .05)

# NaiveSeasonal
from darts.models import NaiveSeasonal
# モデル生成
seasonal_model = NaiveSeasonal(K=12)
# 学習
seasonal_model.fit(train)
# 予測
seasonal_forecast = seasonal_model.predict(len(val))
# 可視化
series.plot(label='actual', lw=1)
seasonal_forecast.plot(label='forecast', lw=1)
plt.legend()
plt.xlabel('Year')

スクリーンショット 2021-01-14 13.47.02.png
スクリーンショット 2021-01-14 13.47.08.png

12ヶ月周期を考慮したモデルを採用しております。

NaiveDrift

次は、NaiveDriftを使ってみます。

# NaiveDrift
from darts.models import NaiveDrift
# モデル生成
drift_model = NaiveDrift()
# 学習
drift_model.fit(train) 
# 予測
drift_forecast = drift_model.predict(len(val))
# 可視化
series.plot(label='actual', lw=1)
drift_forecast.plot(label='forecast', lw=1)
plt.legend()
plt.xlabel('Year')

スクリーンショット 2021-01-14 13.50.51.png

適切に予測できていないので、NaiveSeasonalとNaiveDriftを組み合わせたモデルを使ってみます。

combined_forecast = drift_forecast + seasonal_forecast - train.last_value()

series.plot()
combined_forecast.plot(label='combined')
drift_forecast.plot(label='drift')
plt.legend()
plt.xlabel('Year')

スクリーンショット 2021-01-14 13.52.12.png

Prophet

次はFacebookが開発した時系列予測のオープンソースソフトウェア(OSS)ライブラリ「Prophet」です。詳細は割愛します。

# Prophet
from darts.models import Prophet
# モデル生成
model_P = Prophet()
# 学習
model_P.fit(train) 
# 予測
prediction = model_P.predict(len(val))
# 可視化
series.plot(label='actual', lw=1)
prediction.plot(label='forecast', lw=1)
plt.legend()
plt.xlabel('Year')

スクリーンショット 2021-01-14 13.53.51.png

Theta(シータ法)

次はHyndman & Billahの時系列予測手法「シータ法」です。

# Search for the best theta parameter, by trying 50 different values
# Theta
from darts.models import Theta
thetas = 2 - np.linspace(-10, 10, 50)

best_mape = float('inf')
best_theta = 0

for theta in thetas:
    model = Theta(theta)
    model.fit(train)
    pred_theta = model.predict(len(val))
    res = mape(val, pred_theta)

    if res < best_mape:
        best_mape = res
        best_theta = theta

# Theta
# モデル生成
best_theta_model = Theta(best_theta)
# 学習
best_theta_model.fit(train) 
# 予測
pred_best_theta = best_theta_model.predict(len(val))
# 可視化
series.plot(label='actual', lw=1)
pred_best_theta.plot(label='forecast', lw=1)
plt.legend()
plt.xlabel('Year')

スクリーンショット 2021-01-14 13.55.39.png

Transformer Model

2017年12月頃にGoogleチームにより発表され世間の注目を集めた、時系列データ処理に革新をもたらしたモデルです。
データの与え方を少し工夫する必要がありますが、これも簡単に使えます。

#Transformer Model
from darts.models import TransformerModel

# Normalize the time series 
scaler = Scaler()
train_scaled = scaler.fit_transform(train)
val_scaled = scaler.transform(val)
series_scaled = scaler.transform(series)

N_EPOCHS = 1000
NR_EPOCHS_VALIDATION_PERIOD = 10

# Number of previous time stamps taken into account.
SEQ_LENGTH = 10
# number of output time-steps to predict
OUTPUT_LEN = 1
# Number of layers in encoder/decoder
NUM_LAYERS = 6

my_model = TransformerModel(
    batch_size = 32,
    input_length = SEQ_LENGTH,
    input_size = 1,
    output_length = OUTPUT_LEN,
    output_size = 1,
    n_epochs = N_EPOCHS,
    model_name = 'air_transformer',
    log_tensorboard=True,
    nr_epochs_val_period = NR_EPOCHS_VALIDATION_PERIOD,
    model = None,
    d_model = 64,
    nhead = 32,
    num_encoder_layers = 6,
    num_decoder_layers = 6,
    dim_feedforward = 2048,
    dropout = 0.1,
    activation = "relu",
    custom_encoder = None,
    custom_decoder = None,
    random_state=42,
)

my_model.fit(training_series=train_scaled, val_training_series=val_scaled, verbose=True)

def eval_model(model, n, series, val_series):
    pred_series = model.predict(n=n)
    plt.figure(figsize=(8,5))
    series.plot(label='actual')
    pred_series.plot(label='forecast')
    plt.title('MAPE: {:.2f}%'.format(mape(pred_series, val_series)))
    plt.legend();

eval_model(my_model, 26, series_scaled, val_scaled)

スクリーンショット 2021-01-14 13.59.11.png

LSTM

LSTM(long short-term memory:長短期記憶)とは、ニューラルネットワークの中間層の構造の一つで、自身の出力を、再帰的に入力するような構造を持ったものです。

# LSTM
# Number of previous time stamps taken into account.
SEQ_LENGTH = 12
# Number of features in last hidden state
HIDDEN_SIZE = 25
# number of output time-steps to predict
OUTPUT_LEN = 1
# Number of stacked rnn layers.
NUM_LAYERS = 1

LSTM_model = RNNModel(
    model='LSTM',
    output_length=OUTPUT_LEN,
    hidden_size=HIDDEN_SIZE,
    n_rnn_layers=NUM_LAYERS,
    input_length=SEQ_LENGTH,
    dropout=0.4,
    batch_size=16,
    n_epochs=400,
    optimizer_kwargs={'lr': 1e-3}, 
    model_name='Air_RNN',
    log_tensorboard=True,
    random_state=42
)

LSTM_model.fit(train_scaled, val_training_series=val_scaled, verbose=True)

eval_model(LSTM_model, 26, series_scaled, val_scaled)

スクリーンショット 2021-01-14 14.06.34.png

bestLSTM_model = RNNModel.load_from_checkpoint(model_name='Air_RNN', best=True)
eval_model(bestLSTM_model, 26, series_scaled, val_scaled)

スクリーンショット 2021-01-14 14.08.20.png

GRU

GRUは,LSTMの代替となる手法です。GRUのメリットは、LSTMと比較してパラメーター数が少ないので計算時間が抑えられるといった点があります。性能としてはLSTMと同程度らしいです。

# GRU
gru_model = RNNModel(
    model='GRU',
    output_length=OUTPUT_LEN*4,
    input_length=SEQ_LENGTH,
    hidden_size=HIDDEN_SIZE,
    n_rnn_layers=NUM_LAYERS,
    batch_size=64,
    n_epochs=1500,
    dropout=0.2,
    model_name='Air_GRU_out12',
    log_tensorboard=True,
    random_state=42
)

gru_model.fit(train_scaled, val_training_series=val_scaled, verbose=True)

eval_model(gru_model, 26, series_scaled, val_scaled)

スクリーンショット 2021-01-14 14.12.21.png

さいごに

最後まで読んで頂き、ありがとうございました。
こんな便利なライブラリーがあるとは知りませんでした。PyCaretといい、世の中便利なライブラリーが多いですね。きちっと理論を理解し、自分でも実装できるようになる必要があるとは思いますが、製造現場で普及させるには簡単に使えることも重要かなと思っております。

訂正要望がありましたら、ご連絡頂けますと幸いです。

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Sign upLogin
15
Help us understand the problem. What are the problem?