【時系列】状態空間モデルをわかりやすく解説|カルマンフィルタの仕組み

こんにちは、青の統計学です。

今回は、状態空間モデルについて解説いたします。

MMMと並び広告効果の予測に使われたりと実務での応用も可能な時系列モデルですが、チューニングや実装の難易度が高いという点もあります。

状態の概念から詳しく解説していきます。

状態空間モデル|状態とは

まずは、状態の概念から理解しましょう。

下の図のように、時刻nまでのデータが手元にあると考えます。

\(n\)が現在なので、現在と過去が張る空間が左側になります。

そして、我々が解きたい問いは「将来の予測をする」なので、右側が将来と現在\(n\)が張る空間です。

状態とは、時刻nまでの情報のうち、そのモデルが将来の予測に必要な情報を集約したもの、です。

下は1次のARモデルなので、1期前のデータしか必要ないので状態\(x_n\)については、以下のようになります。

$$x_n=y_n$$

状態空間モデルの構成

では、実際に状態空間モデルの構成をみてみましょう。

大きく分けて、システムモデル(状態方程式と呼んだりします)と観測モデルの2つに分かれています。

まず、システムモデルです。

$$x_n=F_nx_{n-1}+G_nv_n$$

\(x_n\)：m次元の状態ベクトルです。

\(v_n\)：k次元のホワイトノイズ。

\(F_n\)：状態遷移行列

システムモデルは、直接観測できない状態\(x_n\)がどのように時間変化するのかを表現するモデルです。

たとえば、マーケティングチャネルの収益分配最適化を目的とした場合(MMMが有名ですが)、説明変数にはSEOやリスティング、TV CM経由のコンバージョンが入ると思います。

こうした説明変数に加えて、状態ベクトルには、\(状態_t\)が加わるよというイメージです。

$$

\begin{equation}
x_t=\left(\begin{array}{c}
S E O_t \
\text { Listing }_t \
\text { TVCM }_t \
\text { CampaignFLG }_t \
\text { 状態 }_t
\end{array}\right)
\end{equation}

$$

まとめると、説明変数+時間とともに動く切片で構成されています。

次は、観測モデルです。

$$y_n=H_nx_{n}+w_n$$

\(y_n\)：l次元の時系列

\(w_n\)：l次元のホワイトノイズ。

多くの場合、観測モデルが回帰モデル、そして状態\(x_n\)が回帰係数、そしてシステムモデルが回帰係数がどう変化するのかを表現したモデルとして扱われます。

ARモデルを状態空間表現にしてみる

では、時系列表現の中でも基礎的なARモデルを状態空間表現に直してみましょう。

ARモデルについて復習したい方は、以下のコンテンツをご覧ください。

【時系列】ARモデルをわかりやすく解説|Yule-Walker法や最尤法も

$$y_n=H_nx_{n}+w_n$$

$$x_n=F_nx_{n-1}+G_nv_n$$

$$x_n=\begin{equation}\begin{bmatrix}y_n \\ y_{n-1}\\ \vdots \\ y_{n-m+1} \end{bmatrix}\end{equation}$$

$$G=\begin{equation}\begin{bmatrix}1 \\ 0\\ \vdots \\ 0 \end{bmatrix}\end{equation}$$

$$H=\begin{equation}\begin{bmatrix}1 & 0 & \cdots & 0 \end{bmatrix} \end{equation}$$

$$Q=σ^2,R=0$$

$$F=\begin{equation}\begin{bmatrix}a_1 & a_2 & \cdots & a_m \\ 1 & & & \\ & \ddots & & \\ & & 1 & \\ \end{bmatrix} \end{equation}$$

ARモデルの場合は、時間によってシステムモデルが変わらないということがわかりました。

システムノイズ\(v_n\)と観測ノイズ\(w_n\)の相関が1になっている場合には、イノベーション表現というモデル表現があるのですが、やや込み入っており取り上げません。

ARモデルを表現したように、状態空間モデルは線形の時系列モデルを統一的に取り扱うことができるメタモデルであることがわかります。

状態推定とカルマンフィルタ

さて、状態空間モデルではどのように状態推定をしていくのでしょうか？

今回紹介するカルマンフィルタは、ノイズの含まれる観測値からシステムの状態を推定することを目的としたアルゴリズムです。

観測変数\(Y\)に基づき、状態\(x_n\)を知りたい→つまり条件付き分布\(p(x_n|Y_j)\)を知りたいです。

ノイズはホワイトノイズ、状態と観測変数も正規分布に従うことから以下のような関係が成り立ちます。

$$p(x_n|Y_j)〜N(x_{n|j},V_{n|j})$$

上のように条件付き分布\(p(x_n|Y_j)\)は、条件付き期待値と条件付き分散(共分散行列)がわかれば良いです。

カルマンフィルタの手順

カルマンフィルタの手順は、大きく分けて「予測フェーズ」と「フィルタフェーズ」にわかれます。

予測フェーズは一期先の条件付き期待値と条件付き分散の予測を状態方程式をもとに行い、フィルタフェーズではデータを更新します。

データをもとに更新というのは、上で言うと\(y_{n+1|n}\)で\(n+1\)時点でデータの予測をしていましたが、実際のデータ\(y_{n+1}\)を与えると言うことです。

①初期分布を仮定する

②予測フェーズ(predict)

ここでは、一期先予測をすることにします。

条件付き期待値は以下のように計算します。

状態方程式通りに、一期先を予測しているだけですね。

$$x_{n|n-1}=F_nx_{n-1|n-1}$$

共分散行列は、以下のようになります。

この共分散行列は、時刻\(n\)における予測の不確かさを表しています。

$$V_{n|n-1}=F_nV_{n-1|n-1}F_n^T+G_nQ_nG_n^T$$

③フィルタフェーズ(update)

フィルタでは、カルマンゲイン\(K_n\)という値が使われ、予測と観測の情報を組み合わせて状態推定を更新するために使用されます。

以下のように計算されます。

$$K_n=V_{n|n-1}H_n^T(H_nV_{n|n-1}H_n^T+R_n)^{-1}$$

お忘れの方も多いかもしれませんが、\(H_n\)は観測行列です。

これは、状態変数から観測変数へのマッピングを行います。

そして、\(R_n\)は観測ノイズの共分散行列です。

これは、観測データに含まれるノイズの不確かさを表しています。

カルマンゲインは、観測に基づく状態推定値の更新にどれだけの重みを置くかを決定します。

観測の不確かさが大きい場合、ゲインは小さくなり、新しい観測による影響が小さくなります。

一方で、予測の不確かさが大きい場合、ゲインは大きくなり、新しい観測による影響が大きくなります。

では、このカルマンゲインを状態更新に使います。

$$x_{n|n}=x_{n|n-1}+K_n(y_n-H_nx_{n|n-1})$$

ちなみに、式変形して見ると、\(x_{n|n}\)は\(x_{n|n-1}\)と予測誤差の一次結合として表されていることがよくわかります。

$$x_{n|n}=K_ny_n+(I-KH_n)x_{n|n-1}$$

最後に共分散行列の更新です。

$$V_{n|n}=(I-K_nH_n)V_{n|n-1}$$

このカルマンフィルタを反復することで、システムモデルの状態を最適にしていきます。

一期先の予測とフィルタを行う部分の計算を行うだけで全データを使ったことになる(=効率的)と言うことです。

ここまで、状態空間モデルの中身を解説してきましたが、単なる重回帰モデルの拡張ではないことが分かりましたでしょうか。

状態方程式と観測方程式の両方を同時推定する時点で、モデルを規定→残差を最小化するアプローチとは異なります。

また、線形回帰だと各説明変数は独立というきつめの仮定を置いている時点で、時系列データを扱うのは難しいです。