状態空間モデルとは、時系列データにおいて、見えない状態を推定し、その推定した状態から観測値を求めるモデルのことです。状態を介して観測値を予測することにより自由度の高いモデルとなっています。状態の表現力の高さは凄まじく、トレンドや季節成分、更には自己回帰モデルを表すこともできるのです。そんな願ったりかなったりなモデルですが、自由度が高いことの代償として、学習の難易度も高くなっているかと思います。
自由度の高さに加えて、定常・非定常の両方を扱えるという利点もあります。自己回帰モデルでは定常的な時系列データしか扱えなかったため、かなりの進歩です。
この記事では、状態空間モデルってどんな感じなのかを数式で見て、その後にトレンド推定モデルをPythonで実装します。本当は、季節性と自己回帰を含んだ状態の状態空間モデルまでを書きたかったのですが、だいぶ長くなったので次回に回します。なので、予測精度は微妙なトレンド推定モデルを作るところまでです。今回のファイルはこちらに置いてあります。
状態空間モデルとは
状態空間モデルでは、観測値と状態の式をそれぞれ作ります。そのことで、より柔軟な時系列モデリングが可能になっているのです。世の天然と呼ばれる方の発言を分析しても発言のパターンはわからないかもしれませんが、発言の背後にある考え方は似通っていると思われます。発言の背後の状態がわかることで、天然の発言を予測できるようになりますね、きっと。それでは、観測値と状態の２つの式を作ります。
\begin{eqnarray}
\mathbf {y}_t & = & \mathbf {H}(\mathbf{x}_t) + \mathbf {w}_t \\
\mathbf {x}_t & = & \mathbf {F}(\mathbf{x}_{t-1}) + \mathbf {G}(\mathbf{v}_t)
\end{eqnarray}
状態空間モデルは、この２つの式で表されます。実にシンプルです。では、これらの式で使っている変数の定義を書いていきます。
- \( \mathbf { y_t } \)：時刻 \( t \) の観測値。ベクトルの場合があるが、大抵はスカラー
- \( \mathbf { x_t } \)：時刻 \( t \) の状態ベクトル
- \( \mathbf { H } \)：状態を観測値に変換する関数。線形、非線形どちらの場合もあり、線形の場合は行列で表し、観測行列と呼ばれる
- \( \mathbf { F } \)：状態を次の状態に変換する関数。線形、非線形どちらの場合もあり、線形の場合は行列で表し、状態推移行列と呼ばれる
- \( \mathbf { w_t } \)：時刻 \( t \) の観測ノイズ
- \( \mathbf { v_t } \)：時刻 \( t \) の状態ノイズ
- \( \mathbf { G } \)：時刻 \( t \) の状態ノイズを変換する線形 or 非線形関数
上記の状態変数が未知のパラメータであり、この値を特定することが状態空間モデルのミッションとなります。状態と観測値以外のパラメータは、予め設定します。いや、それどう設定するんだよってなると思いますが、これはもう慣れですね。
状態空間モデルの式が定義できたので、式の中身に関してじっくり見てみましょう。
状態空間モデルの式からわかること
状態空間モデルの数式の意味を見ていきます。まずは観測値を推定する式に注目しましょう。
\begin{eqnarray}
\mathbf{y}_t = \mathbf{H}(\mathbf{x}_t) + \mathbf{w}_t
\end{eqnarray}
上の式を見ると、観測値は同時刻の状態の値から求められていることがわかります。その横についている \( \mathbf{w}_t \) はノイズであり、この式の大部分を担っているのは \( \mathbf{x}_t \) です。また、過去の観測値の値が使われておらず、これが自己回帰モデルと違います。つまり、状態の値から観測値が求まるということだけを覚えておけばOKなわけです。
次に状態を推定する式を見ていきます。
\begin{eqnarray}
\mathbf{x}_t = \mathbf{F}( \mathbf{x}_{t-1}) + \mathbf{G}(\mathbf{v}_t)
\end{eqnarray}
上の式を見ると、状態の値は直前の状態の値から推定されていることがわかります。それに、一定値のノイズが加わっています。状態 \( \mathbf{x}_t \) にくっついている状態遷移行列（線形の場合）で時系列データの様々な形式を扱えるようにしており、これにより、トレンド、季節周期、自己回帰などを表すことができるのです。この遷移行列は予め設定するもので、状態ベクトル \( \mathbf{x}_t \) を求めることこそがミッションであり、そこに全集中しましょう。
以上の話から、状態空間モデルのロジックを肌身に感じることができましたね。武者震いがしてくる方もいるでしょう。これまでの話を図でまとめると、以下のようになります。
観測値は、同時刻の状態から推定され、状態は直前の過去の状態から推定されます。実はたったこれだけの話なのであり、考え方は非常にシンプルですね。
状態空間モデルの実装
状態空間モデルがどのようなものかなんとなくわかってきたので、Pythonを使って実践により状態空間モデルの世界を体験していきます。状態の変換関数には、線形、非線形がありますが、今回は線形のモデルを構築し、かつ、状態と誤差項が正規分布に従うことを仮定する、線形ガウスモデルに挑戦します。
また、今回扱う状態空間のモデルは、自由度が高いpykalmanです。pipで簡単にインストールできます。それでは、実際にやっていきましょう。まずは時系列データを取得します。
時系列データの取得
状態空間モデルを構築するのに使う時系列データを用意します。今回扱うデータは、おなじみの月別の飛行機の乗客数のデータです。
必要なモジュールのインポート
from pathlib import Path import numpy as np import pandas as pd import requests import io import copy from tqdm import tqdm import matplotlib.pyplot as plt import seaborn as sns import scipy from scipy.optimize import minimize from pykalman import KalmanFilter
データの取得
飛行機の乗客数のデータを指定したURLから取得します。
url = "https://www.analyticsvidhya.com/wp-content/uploads/2016/02/AirPassengers.csv" stream = requests.get(url).content passengers = pd.read_csv(io.StringIO(stream.decode('utf-8')))
月別の飛行機の乗客数のデータが取得できましたね。グラフ化して時系列の様子を見てみます。
passengers.plot(
x='Month',
y='#Passengers',
figsize=(8, 5),
xlabel=''
)
グラフを見ると、今回扱うデータは段々増加しています。つまり、トレンドがあるデータであり、非定常なデータだと思われます。なので、自己回帰モデルでは、変換しない限り扱えません。
データの加工
Pykalmanで扱える時系列データにするために加工しましょう。その際に、モデル構築用と検証用データに分けておきます。
y = pd.Series(
passengers['#Passengers'].values,
index=pd.to_datetime(passengers['Month'], infer_datetime_format=True)
)
y = y.astype(float)
n_train = 120
train_data, test_data = y.values[:n_train], y.values[n_train:]
以上で準備は終了です。次からは、状態空間モデルの構築に入ります。
状態空間モデルの構築 〜トレンド推定モデル〜
状態空間モデルは、状態の推移行列を拡張していくことで複雑なモデルを構築できます。ですが、いきなり複雑なモデルから始めると、理解せずにただプログラムを実行するマシーンとかしてしまうので、簡単なモデルから試します。それは、トレンド推定モデルです。
トレンドは、状態を直近の状態で表すモデルです。式で表すと、以下の形になります。
\begin{eqnarray}
\mathbf { y_t } & = & \mathbf { tre }_t + \mathbf {w}_t , \quad \mathbf{w}_t \sim \mathbf{ N }(0, \sigma_{w}^2 ) \\
\mathbf {tre}_t & = & \sum_{ i=1 }^k { c_i^{ (k) } \mathbf {tre}_{ t-i } } + \mathbf{v}_{t}, \quad \mathbf{v}_{t} \sim \mathbf{N}(0, \sigma_{tre}^2 )
\end{eqnarray}
上記の式で、次数の \(k\) と係数 \( c \) を指定すれば、トレンドモデルの形は決定です。今回はよく使われる次数2の以下の式で実装します。
\[
\mathbf{tre}_t = 2\mathbf{tre}_{t-1} - \mathbf{tre}_{t-2} + \mathbf{v}_t
\]
トレンド推定モデルの式が決まったので、最初に定義した状態空間モデルの各行列は、以下のようになります。
\begin{eqnarray}
\mathbf{F} = \left[
\begin{matrix}
2 & -1 \\
1 & 0
\end{matrix} \right], \quad
\mathbf{x}_t & = & \left[
\begin{matrix}
\mathbf{tre}_{t - 1} \\
\mathbf{tre}_{t - 2}
\end{matrix} \right], \quad
\mathbf{G} = \left[
\begin{matrix}
1 \\
0
\end{matrix} \right] \\
\mathbf{H}^T & = & \left[
\begin{matrix}
1 \\
0
\end{matrix} \right]
\end{eqnarray}
トレンド推定モデルのパラメターの設定
上の定義に合わせてトレンド推定モデルのパラメータを愚直に作っていきましょう。まぁ、行列とベクトルを作っていくだけですが。
n_dim_obs = 1 # 観測値の次元数 n_dim_trend = 2 # トレンドの次元数（状態の次元数） n_dim_state = n_dim_trend F = np.array([ [2, -1], [1, 0] ], dtype=float) G = np.array([ [1], [0] ], dtype=float) H = np.array([ [1, 0] ], dtype=float) Q = np.eye(1) * 10 Q = G.dot(Q).dot(G.T) state_mean = np.zeros(n_dim_state) # 状態の平均値ベクトルの初期値 state_cov = np.ones((n_dim_state, n_dim_state)) # 状態の分散共分散行列の初期値
あとは、実際の \( tre_t \) の値を求めるだけです。これはPykalmanさんに任せれば秒で求まります。実際に、どのようなやり方で求めているかは、紙面上書ききれない（便利な言葉で逃げました）ので、ぜひブログの最後に挙げた参考文献を御覧ください。
トレンド推定モデルの構築
パラメータの設定ができたので、状態 \( tre_t \) を推定して、状態空間モデルを完成させます。トレンド推定モデルは状態が線形回帰に従うと仮定しており、状態の平均と共分散を求めれば状態の値の分布が決まります。ゴールは目の前です。ただ、そんなに気負う必要はなく、実はサクッとできてしまいます。
kf = KalmanFilter(
n_dim_obs=n_dim_obs,
n_dim_state=n_dim_state,
initial_state_mean=state_mean,
initial_state_covariance=state_cov,
transition_matrices=F,
transition_covariance=Q,
observation_matrices=H,
observation_covariance=1.0,
)
以上です！
トレンド推定モデルができました。ついに、念願の状態の平均と共分散が求まったのです。果報は寝て待てといいますが、流行る気持ちを抑えながら、トレンド推定モデルの推定結果を見ていきましょう。
トレンド推定モデルの推定結果を確認 〜平滑化〜
トレンド推定モデルの推定結果を見ていきます。まず確認するのは、与えられたデータの期間内の状態の推定ができているかです。状態空間モデル界隈の中では、これを平滑化と呼びます。平滑化では、与えられた期間のデータを使って、その期間内の状態を推定します。平滑化の良いところは、与えらた期間内に欠損値があっても期間内のどの時刻の状態でも推定ができるので、欠損値の時系列の値を補えることです。では、平滑化してみます。
state_means, state_covs = kf.smooth(train_data)
上記のように、１行で平滑化した状態の平均値と共分散を推定しています。
print('状態の平均値 : \n{} \n\n 状態の共分散 : \n{}'.format( state_means[:5], state_covs[:5] ))
上記の結果を見ると、状態の平均値は時間 \( t \) ごとに \( x_t, x_{t-1} \) が求められていることがわかります。また、状態の共分散は時間 \( t \) ごとに \( t, t-1 \) の共分散行列を出しているので、３次元になっています。
これで終わりかと思って気を緩めてはいけません。先程はあくまでも状態の推定値であり、観測値ではないです。この結果を使って、観測値に変換します。正規表現の線形変換なので、観測値も正規表現になります。なので、観測値は、観測行列に推定した状態の平均値をかけることで求めます。
ovsevation_means_predicted = state_means.dot(H.T)
これで、平滑化した観測地の値が求まりました。この中身を見ますが、せっかくなので、学習データと重ねて表示します。
plt.figure(figsize=(8, 6)) plt.plot(train_data, label="observation") plt.plot(ovsevation_means_predicted, '--', label="forecast") plt.tight_layout() plt.savefig(result_dir_path.joinpath('passenger_trend_smooth.png'))
見事に学習データと平滑化した観測値が重なっていることがわかります。めでたし、めでたしですが、話はこれで終わりません。時系列モデルの真骨頂である、予測を行ってみます。
トレンド推定モデルの結果確認 〜フィルタリングと１期先予測〜
与えらたデータの期間内の観測値を推定することができたので、予測を行います。状態空間モデルで予測で第一に行うことは、与えらた時系列データの最新の時間である \( t \) の状態の推定です（これをフィルタリングと呼ぶ）。次に、推定した時間 \( t \) 状態の平均ベクトル、共分散行列を使って、１期先の状態を予測します。これを繰り返すことで、どんどん予測していけるのです。
train_data, test_data = y[:n_train], y[n_train:]
state_means, state_covs = kf.smooth(train_data)
ovsevation_means_predicted = np.dot(state_means, kf.observation_matrices.T)
current_state = state_means[-1]
current_cov = state_covs[-1]
pred_means = np.array([])
for i in range(len(test_data)):
current_state, current_cov = kf.filter_update(
current_state, current_cov, observation=None
)
pred_mean = kf.observation_matrices.dot(current_state)
pred_means = np.r_[pred_means, pred_mean]
変数 pred_means に予測した結果を格納しています。上のコードを見ればわかるように、はじめに時間 \( t \) の状態の平均ベクトル・共分散行列を求め（current_state, current_cov）、その値を使って１期先の予測を行っています。この結果をグラフ化しましょう。
plt.figure(figsize=(8, 6)) plt.plot(y.values, label="observation") plt.plot( np.hstack([ ovsevation_means_predicted.flatten(), np.array(pred_means).flatten() ]), '--', label="forecast" ) plt.tight_layout() plt.savefig(result_dir_path.joinpath('passenger_trend_predict.png'))
予測結果を見てもらえば分かる通り、最初の予測の値まではある程度正しそうですが、真の値からどんどん離れています。ちょっとした値の振れ幅や急激な変化に全く対応できていません。まぁ、モデルが簡単なのでこの結果は妥当ですね。
状態空間モデルでは状態の共分散行列も推定しているので、そこから信頼区間を求めることができます。先程のコードに信頼区間を求める記述を書き加えるだけなので、一気に書いてしまいましょう。
train_data, test_data = y[:n_train], y[n_train:]
state_means, state_covs = kf.smooth(train_data)
ovsevation_means_predicted = np.dot(state_means, kf.observation_matrices.T)
ovsevation_covs_predicted = kf.observation_matrices.dot(
np.abs(state_covs)
).transpose(1, 0, 2).dot(kf.observation_matrices.T)
lowers, uppers = scipy.stats.norm.interval(
0.95,
ovsevation_means_predicted.flatten(),
scale=np.sqrt(ovsevation_covs_predicted.flatten())
)
current_state = state_means[-1]
current_cov = state_covs[-1]
pred_means = np.array([])
tmp_lowers = []
tmp_uppers = []
for i in range(len(test_data)):
current_state, current_cov = kf.filter_update(
current_state, current_cov, observation=None
)
pred_mean = current_state.dot(kf.observation_matrices.T)
pred_cov = kf.observation_matrices.dot(
np.abs(current_cov)
).dot(kf.observation_matrices.T)
pred_means = np.r_[pred_means, pred_mean]
lower, upper = scipy.stats.norm.interval(
0.95, pred_mean, scale=np.sqrt(pred_cov)
)
tmp_lowers.append(lower)
tmp_uppers.append(upper)
lowers = np.hstack([lowers, np.array(tmp_lowers).flatten()])
uppers = np.hstack([uppers, np.array(tmp_uppers).flatten()])
plt.figure(figsize=(8, 6))
plt.plot(y.values, label="observation")
plt.plot(
np.hstack([
ovsevation_means_predicted.flatten(),
pred_means.flatten()
]),
'--', label="forecast"
)
plt.fill_between(range(len(y)), uppers, lowers, alpha=0.5, label="credible interval")
plt.legend()
plt.tight_layout()
plt.savefig(result_dir_path.joinpath('passenger_trend_predict.png'))
予測結果がだめだめなので、信頼区間がとんでもないことになっていますね。この結果を改善するために、状態を複雑化していきます。追加するのは、季節周期と自己回帰です。これをやっていきたいのですがすでに長くなっているので、次の記事で書きます。
参考文献
今回の記事では、状態空間モデルの式の説明と実装に重きを置いたため、フィルタリング、１期先予測、平滑化をどうやっているかは書けませんでした。詳しく知りたい方は下の本をぜひご覧ください。
