ガウス過程で時系列予測する

概要

ガウス過程を活用した時系列モデル（特に予測）ってあまり聞かないなと思い，今回試してみました． 結論から言うと，特徴量を工夫すれば使えないこともないですが，あえてガウス過程を使う理由はなさそうです．

私のガウス過程に対する浅学により，知識のある方からすると「そんなの当たり前だろ」と言われてしまいそうな結論になってしまいましたが，ご笑覧いただければ幸いです．

実行コードはこちらです．

github.com

ガウス過程とは

まだまだ筆者も勉強中なので詳しい説明ができなくて恐縮なのですが， 共分散行列などが解法にででくるところから推測するに，同じような特徴量を持つときは同じような目的変数になるだろうみたいな予測をするのだと認識しています．

以下のページが分かりやすかったです．

datachemeng.com

使用データ

prophetの公式ドキュメントで使用されている以下のデータを使用しました．

github.com

wgetが入っている方は以下をjupyterとかで実行すればダウンロードできます．

wget https://raw.githubusercontent.com/facebook/prophet/main/examples/example_wp_log_peyton_manning.csv

グラフに可視化すると以下のようになります． 年ごとの動きが似ており，年間の季節性があると考えることができます． なお，ガウス過程が平均０の変動を仮定しているみたいだったので，ダウンロードしたデータを標準化したものを扱うこととします（標準化しなくても良い方法あればご教示ください）．

df = pd.read_csv('./example_wp_log_peyton_manning.csv', parse_dates=['ds'])

# 2015年以降をテストデータ，それ以外を学習データとする
df['is_train'] = df['ds'] < pd.to_datetime('2015-01-01')

# ガウス過程は目的変数の平均０を仮定しているようなので，標準化する
df['y'] = (df['y'] - df['y'].mean()) / df['y'].std()

plt.figure(figsize=(12, 8))
plt.scatter(df['ds'], df['y'], color='#999999', label='y_true')
plt.xlabel('ds')
plt.grid()
plt.legend(loc='upper left')
plt.show()

使用するデータ

三角関数で季節成分を表現

年次の季節性が確認できたので，ここでもprophetの真似をして 以下のフーリエ級数のような式をN＝10まで特徴量として生成します．

ここで，は各三角関数の重みです． 回帰モデルの際は必要ですが，今回はガウス過程なので必要ありません． またpは１周期の期間です．今回は年次の周期性を表現したいので，とします．

コードで表現すると以下のようになります．

def add_fourier(df: pd.DataFrame, period: float, order: int) -> pd.DataFrame:
    df[f'fourier_O{order}_sin'] = np.sin(2 * np.pi * order * np.arange(len(df)) / period)
    df[f'fourier_O{order}_cos'] = np.cos(2 * np.pi * order * np.arange(len(df)) / period)
    return df

for i in range(10):
    df = add_fourier(df, period=365.25, order=i+1)

詳しく知りたい方はprophetの論文を読んでみてください．

peerj.com

学習・予測を実行してみる

今回，ガウス過程の計算にはGPyというライブラリを使用します．

github.com

カーネルを適当に設定して，とりあえず動くか見てみます． 必要なライブラリをimportします．

import GPy
from GPy import kern

カーネルの設計や，学習予測は次のように行うようです．

kernel = kern.RBF(input_dim=20) + kern.Bias(input_dim=20)

x_columns = [c for c in df.columns.tolist() if 'fourier' in c]

m = GPy.models.GPRegression(
    X=df.loc[df['is_train']==True, x_columns].values,
    Y=df.loc[df['is_train']==True, 'y'].values[:, None],
    normalizer=None
)

yhat, yhat_var = m.predict(df.loc[:, x_columns].values)
yhat_std = yhat_var ** .5

予測結果を可視化してみましょう．

plt.figure(figsize=(12, 8))
plt.scatter(df['ds'], df['y'], color='#999999', label='y_true', s=10)
plt.plot(df['ds'], yhat, label='y_pred', color='#377eb8')
plt.fill_between(
    df['ds'],
    yhat.reshape(-1) - yhat_std.reshape(-1),
    yhat.reshape(-1) + yhat_std.reshape(-1),
    color='#377eb8',
    alpha=.2,
    label='1 sigma'
)
plt.vlines(pd.to_datetime('2014-12-31'), -4, 5.5, label='End of training', color='k')
plt.xlabel('ds')
plt.grid(linestyle='--')
plt.legend(loc='upper left')
plt.show()

予測結果

予測期間（End of trainingより右側）において，少しスケールが正解値y_trueと予測値y_predで異なってしまっていますが，概ねの形はうまく予測できているのが確認できます．

終わりに

特徴量に無理やり三角関数を導入したあたりで，「それガウス過程でやる意味ある？」という疑問が出てきました． LSTMやProphetのようなモデルで十分そうなことをわざわざ難しくしすぎてしまった感が出てしまいました．

とはいえ，一応ガウス過程でも時系列予測ができそうなことは確認できました．

今後の個人的な課題は，ガウス過程やカーネル法を導入することで何が嬉しいのか勉強することです． 時系列でも使えそうだなということ自体，私の勘違いの部分が大きかったので，この分野は全体的に知識が浅いなと痛感しました．

参考(本文に登場しなかったもの)

nykergoto.hatenablog.jp

gpy.readthedocs.io