移動平均と時系列データの簡単な線形回帰モデル

いまいち時系列データを重点的に扱ったことがなく、いろんな機械学習の教材でもサラッとpandasのrolling関数くらいにとどめて終わっていることが多い印象で、Kaggleのチュートリアルで時系列データと機械学習を扱った教材があるため、それをアウトプットしていきます

なんせKaggleは英語に他ならないため、（私含む、）英語わかんないのよ！という方の一助となれば幸いです。
（ちなみに以下のURLに準拠して進めています。データは自作。）
Linear Regression With Time Series

移動平均の簡単な理論おさらい

たしかKaggleチュートリアルに後程出てきたはずですが、一応先に簡単な移動平均の理論だけ触れておきます
統計検定にも出たことありますしね...

移動平均の意味

区間（一年とか一週間）ごとにその規則性などが見出せそうなときに効果的です。
季節性のなにか要因があるのか？とか一年区間で見れば常に単調増加しているとみなせるよ！みたいなことがわかったりします。

移動平均の計算

区間が偶数か奇数かで分かれますが、中央値とかの計算とほぼ同じと考えて良いと思っています

奇数

区間2k+1に対して中央のn+1からそれぞれ前後にnだけ広げたものの平均を取ります

k=1のとき
$$ s = \frac{x_{n-1} + x_n + x_{n+1}}{3}$$

偶数

偶数はt~t+2k-1, t+1~t+2k区間で区切り、それぞれの平均をさらに平均を取得するイメージです

k=2
$$ s = \frac{\frac{x_t + x_{t+1} + x_{t+2} + x_{t+3}}{4} + \frac{x_{t+1} + x_{t+2} + x_{t+3} + x_{t+4}}{4}}{2}$$

pythonなどではdf.rolling(window=~, periods=~)みたいな感じです（今回は使わない）

線形回帰で時系列データを予測する

では、本題。

（このチュートリアルにはCalenderFourierみたいなフーリエ変換を使う高度なアプローチもあり、技術的にも数学的にも理解したいのですが、まだまだ先になりそうな予感）

時系列データの最初の一歩として、じつは非常にシンプルなアプローチでいろいろわかることがあります。
機械学習は8割が前処理と特徴量エンジニアリングみたいなことを言われていますが、こういう一見簡単だけどスルーしがちなことも着実にしておくといいのかもしれませんね。

データの自作

一応そのままKaggleデータセットを使うのもいいですが、簡単なデータなので自作しました
ので、準備として以下。

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import seaborn as sns

from scipy import stats


rng = np.random.default_rng()

x = np.linspace(100, 400, 300)
# 単調増加に乱数を加えたもの
y = np.round(x + rng.integers(-100, 100, size=300), 0)
# 時系列データの作成
time_series = pd.date_range(start='2019-4-1', periods=300) 

data = pd.DataFrame({'Sales': y}, index=time_series, dtype='int')
plt.plot('Time', 'Sales', data=data)

意外と時系列データのカラム作るの知らなかったので勉強なった。。

ダミー変数でアプローチ

機械学習では時系列のカラムを変数として扱うことも多々ありますが、それをダミー変数で連番にしたものを使って予測するアプローチがあります。

まずはコードみましょう

data['Time'] = np.arange(len(data))

_, ax = plt.subplots(1, 1, figsize=(15, 5))
ax.plot('Time', 'Sales', data=data)
# ci: 回帰を行うときの信頼区間 (%)。ブーストラップ（無作為復元抽出）に使う
sns.regplot(x='Time', y='Sales', data=data, ax=ax, ci=None, scatter_kws=dict(color='0.25'))

seabornのregplotは2次元のデータプロットに対して線形回帰モデルを当て込む関数です
非常にシンプルですが、このグラフからこの連番を割り当てた特徴量'Time'により、時間依存性のmodelを作成することができました。

大まかに単調増加していることがこれだけで推測できますね。

shift関数で特徴量を作成

次に前日と当日の関係性を見てみる、という手法があります
（時間のラグをみることからLag featuresと書かれていました。）

こちらも非常にシンプルで、前日のデータを反映したカラムをつくり、線形モデルを当て込んでみるという流れです。

data['Lag_1'] = data['Sales'].shift(1)

fig, ax = plt.subplots(figsize=(10, 10))
ax = sns.regplot(x='Lag_1', y='Sales', data=data, ci=None, scatter_kws=dict(color='0.25'))
ax.set_aspect('equal')

実際にこんなにうまくいくことも少ないですが、前日のデータと当日のデータには正の相関関係がみられるため、　この前日の販売数も重要な特徴量と言える可能性があります。

Kaggle曰く、このLag_1はserial dependence（意味合いは上記のようなラグとの関連性和訳わからない...　連続的依存性？とか？　）のmodel構築に役立ちます.

今回の例で言えば、大体前日の売り上げが高ければ今日もうりあげは大体良さげ！ってのが予測できるようになります。

ちょっと実践。

先ほどは描画した程度にとどめましたが、実際にはskleanとか使うので、それをサラッと扱って終わります
だいたい同じです。

データ準備

rng = np.random.default_rng()

x = np.linspace(100, 400, 300)
y = np.round(x + rng.integers(-90, 110, size=300), 0)
time_series = pd.date_range(start='2019-4-1', periods=300) 

data = pd.DataFrame({'Sales': y}, index=time_series, dtype='int')
data

Timeダミー変数を用意し、LinearReg

# 描画のための準備（なくてもいいです）
plot_params = dict(
    color="0.75",
    style=".-",
    markeredgecolor="0.25",
    markerfacecolor="0.25",
    legend=False,
)

from sklearn.linear_model import LinearRegression


data['Time'] = np.arange(len(data))

X = data.loc[:, ['Time']]
y = data['Sales']

model = LinearRegression()
model.fit(X, y)

y_pred = model.predict(X)
y_pred = pd.Series(y_pred, index=X.index)

plt.figure(figsize=(10, 5))
ax = y.plot(**plot_params)
ax = y_pred.plot(ax=ax,  linewidth=3)
ax.set_title('Time plot of Sales data!')
plt.show()

急にunpacking（**）とか使っちゃってますが、「ax = y.plot()」でも全然動きます

Lag features

LinearRegでは欠損値のあるデータはエラーが出ちゃうので、何かしら補完（0とかbackfillingという欠損値のあるデータの次の日のデータを埋め込むなど）がありますが、今回はせいぜい1ssのみなので、dropnaで削除してしまいます。

data['Lag_1'] = data['Sales'].shift(1)

# lineaRegでは欠損値ありのデータ使えないため、注意
X = data.loc[:, ['Lag_1']]
X.dropna(inplace=True)

y = data['Sales']
y, X = y.align(X, join='inner') # xが欠損しているため、yもそれに合わせる

model = LinearRegression()
model.fit(X, y)

y_pred = model.predict(X)
y_pred = pd.Series(y_pred, index=X.index)

_, ax = plt.subplots(1, 1, figsize=(8, 8))
ax.plot(X['Lag_1'], y, '.')
ax.plot(X['Lag_1'], y_pred)
ax.set_aspect('equal')
ax.set_title('Lag plot of Sales data!')
plt.show()

終わり

時系列データは割と扱うことも少なく、そして調べた感じ踏み込んだものが日本語のサイトにはないので、みんな嫌いなのかな〜と思ったりしていますw
このKaggleチュートリアルは適宜作れればい〜な〜なんて思っていますが、割と高度な感じもするので理解するのに時間がかかりそう。。

今回参考にさせていただいた記事や資料

Linear Regression With Time Series
Seaborn で散布図・回帰モデルを可視化する
pandasで時系列データの作成
[Python]pandas.DataFrame.alignについての解説