状態空間モデルの考え方と適用例

前書き

年をとっても思考力は高めていきたいAIエンジニアです。
衰えないように意識して生きていきたいと思います。

今回は時系列分析で利用される状態空間モデルについて適用例も含めてまとめていきます。

状態空間モデル

とある探偵漫画の主人公は自分の得た情報を元に推理を行い、真実を明らかにしています。

一般人が漫画のような推理をすることは難しいですが、"自分が得た情報から何かしらを予測する"ということはよく行われると思います。

例えば

場面	得た情報	予測
対人関係	挨拶の時に相手の表情がいつもより明るかった	良いことでもあったのかな
旅先	道にゴミが多くて壁に落書きが多い	この都市、治安悪そうだな
野球	変化球が2球続けてボール球になった	このピッチャー、変化球のコントロールがあまりよくないな

といった感じです。

人間の場合、自分の得た情報に加えて自分の知識や過去の経験を元に予測していると思います。

---

これと似たようなことを数式を使って表現したものとして状態空間モデルというものが存在します。
※状態空間モデルはあくまで数式を使って表現することなので直接予測する手法ではないのですが、モデル化することで様々な予測手法と組み合わせることが可能になります。
（予測するための前準備みたいなものと考えてください。）

状態空間モデルでは、自分の得た情報を観測値、予測する対象のことを状態と呼びます。

人間が
・自分の得た情報
・知識や経験
を元に予測していくのと同じように、
これら2つのものと予測対象(状態) の関係を数式を使ってモデル化していきます。

---

状態空間モデルでは2つの数式を使ってモデル化します。

1つ目の式は
・観測値(得た情報) と状態(予測対象)
を関連づける式です。

人間が予測のために自分の得た情報を活用する部分を表現します。
この式を観測方程式と呼びます。

簡単な文字式で表現すると
観測値 = 状態 + (観測時の)ノイズ
となります。

---

ここでサラッとノイズが追加されていますが、これは観測値にはノイズが伴うという点を前提にしているからです。
人間で例えると、普段は視力1.0の人も体調が悪い日は0.9とかに落ちてしまうと思います。
そうなると全く同じものを見たとしてもわずかに異なって見えると思います。
このわずかに異なる部分をノイズで表現します。

状態空間モデルではこの観測時に伴うノイズのことを観測ノイズと呼びます。

---

2つ目の式が
・過去の状態と現在の状態
を関連づける式です。

人間が現状の予測をするために過去の情報を活用する部分を表現します。
この式を状態方程式と呼びます。

同じく簡単な文字式で表現すると
現在の状態 = 前の状態 + (状態変化に伴う)ノイズ
となります。

---

ここでもノイズが登場します。
これは状態が遷移するときに何らか不確実なものが伴うという点を前提にしています。
100%確定的に起こる現象というものが存在すればこのノイズも不要になるかもしれませんが、現実的にあり得ないという前提でノイズ項を追加しています。

状態空間モデルではこの状態遷移に伴うノイズのことをプロセスノイズと呼びます。

---

時間の流れによって変化していく状態(予測対象) に対して、
・観測値(得た情報) と状態(予測対象) の関係
・状態の遷移
の2つに分けて考えることで複雑なシステムも表現しやすくなります。

また、それぞれの関係にノイズが含まれることを前提とすることで、ノイズの影響を考慮した本質的な表現が可能となります。

---

状態空間モデルのイメージ

状態空間モデル_要点まとめ

・状態空間モデル
観測値(得た情報) と状態(予測対象) の関係を数式を使ってモデル化(表現) したもの。
観測値と状態の関係を表す観測方程式と状態の遷移を表す状態方程式の2つの式を使って表現する。

・観測方程式
観測値 = 状態 + 観測ノイズ

・状態方程式
現在の状態 = 前の状態 + プロセスノイズ

状態空間モデルの適用例

ここまでの説明は状態空間モデルがどんなものなのかという直感的な説明にとどめてきたので、ここからは具体的な数式を出しつつ書き進めていきます。

線形回帰モデル

まずは、そもそも状態が時間変化しないことを前提としたモデルに適用してみましょう。
（時間変化しないということは時間変化することの特殊な場合なので状態空間モデルでも表現できます。）

わかりやすくするため単回帰モデルで説明していきます。
単回帰モデルは体重と身長の関係などを

$$体重=係数×身長+切片 体重 = 係数 \times 身長 + 切片$$

という形で直線でエイヤーと近似する手法でした。

私が中学生の時に習った変数を使い

$$y=ax+b y = ax + b$$

で進めていきます。

---

今回、予測したい対象(状態) は回帰係数$aa$と切片$bb$になります。

時刻$kk$での状態(回帰係数$aa$と切片$bb$)をそれぞれ$aka_{k}$, $bkb_{k}$とします。
今考えている問題では状態は時間変化しないので、
状態方程式は

$${ak=ak−1bk=bk−1 \left\{ \begin{aligned} a_{k} &= a_{k-1} \\ b_{k} &= b_{k-1} \\ \end{aligned} \right.$$

観測値を$yky_{k}$, 観測ノイズを$hk∼N(0,σh2)h_{k} \sim N(0, \sigma^{2}_{h})$, 入力値を$xkx_{k}$とすると
観測方程式は

$$yk=akxk+bk+hk \begin{aligned} y_{k} &= a_{k}x_{k} + b_{k} + h_{k} \\ \end{aligned}$$

まとめると、

$${ak=ak−1bk=bk−1yk=akxk+bk+hk \left\{ \begin{aligned} a_{k} &= a_{k-1} \\ b_{k} &= b_{k-1} \\ y_{k} &= a_{k}x_{k} + b_{k} + h_{k} \end{aligned} \right.$$

と表現できます。

---

観測方程式と元の式の違いは、時刻を表す添え字$kk$がある点と観測ノイズ$hkh_{k}$がある点です。
時間変化しない想定なので$kk$は無視して大丈夫です。

添字の$kk$を無視すると

$$y=ax+b+h y = ax + b + h$$

観測ノイズ$hh$を切片$bb$とまとめて$bnewb_{new}$という文字で置き換えると

$$y=ax+bnew y = ax + b_{new}$$

となり元の式と同じ形になります。

---

時間変化を許容した線形回帰モデル

ここまでの説明では(時間変化を全く考慮しない) 一般的な線形回帰モデルを状態空間モデルに適用しました。
線形回帰モデルの前提として、パラメータが時間変化しないという点がありました。

ここから、ほんの少し拡張してあげることにしましょう。

1900年の人間の体重と身長の関係は2000年の人間の体重と身長の関係と全く同じでしょうか？
大方同じだと思いますが、社会が変化したこともあり少しは変化しているかもしれません。

そこで、微妙な時間変化をノイズで許容してみます。
時刻$kk$での係数のプロセスノイズを$ik∼N(0,σi2)i_{k} \sim N(0, \sigma_{i}^{2})$, 切片のプロセスノイズを$jk∼N(0,σj2)j_{k} \sim N(0, \sigma_{j}^{2})$とすると
状態方程式は

$${ak=ak−1+ikbk=bk−1+jk \left\{ \begin{aligned} a_{k} &= a_{k-1} + i_{k}\\ b_{k} &= b_{k-1} + j_{k}\\ \end{aligned} \right.$$

観測方程式は（今までと同じで）

$$yk=akxk+bk+hk \begin{aligned} y_{k} &= a_{k}x_{k} + b_{k} + h_{k} \\ \end{aligned}$$

まとめると、

$${ak=ak−1+ikbk=bk−1+jkyk=akxk+bk+hk \left\{ \begin{aligned} a_{k} &= a_{k-1} + i_{k}\\ b_{k} &= b_{k-1} + j_{k}\\ y_{k} &= a_{k}x_{k} + b_{k} + h_{k} \end{aligned} \right.$$

となります。

プロセスノイズを追加しただけですが、係数や切片が時間変化する場合も考慮して表現することができました。

観測誤差と時間変化を考慮した線形回帰モデル

さて、ここまで線形回帰の入力値の$xx$には触れていませんでしたが、この入力値もどこかで測定した値であることがほとんどでしょう。
例えば、身長を入力とする場合は身長の測定時に誤差があると考える方が自然だと思います。

この考えも上記までの式を少し拡張することで表現できます。

$xkx_{k}$を観測値として捉え、観測ノイズを$lk∼N(0,σl2)l_{k} \sim N(0, \sigma^{2}_{l})$として観測方程式を表現してみます。

$${xknew=xk+lkyk=akxknew+bk+hk \left\{ \begin{aligned} x^{new}_{k} &= x_{k} + l_{k}\\ y_{k} &= a_{k}x^{new}_{k} + b_{k} + h_{k} \end{aligned} \right.$$

状態方程式と観測方程式をまとめると、

$${ak=ak−1+ikbk=bk−1+jkxknew=xk+lkyk=akxknew+bk+hk \left\{ \begin{aligned} a_{k} &= a_{k-1} + i_{k}\\ b_{k} &= b_{k-1} + j_{k}\\ x^{new}_{k} &= x_{k} + l_{k}\\ y_{k} &= a_{k}x^{new}_{k} + b_{k} + h_{k} \end{aligned} \right.$$

となります。

このようにまとめることで、入力値$xx$の観測誤差、係数$aa$と切片$bb$の時間変化を考慮してモデル化することができます。

ローカルレベルモデル

ここからは状態が時間変化することを前提としたモデルを利用していきます。
温度管理された室内の温度のモデル化してみましょう。

状態が真の室内温度、観測値は温度計で測定した温度ということになります。
温度計にも測定の誤差はありますが、測定時の温度がそのまま観測値となります。
また室内は温度管理されているので測定のタイミングで温度の変化はないと考えることができます。

時刻$kk$での状態を$xkx_{k}$, プロセスノイズを$pk∼N(0,σp2)p_{k} \sim N(0, \sigma_{p}^{2})$とすると
状態方程式は

$$xk=xk−1+pk x_{k} = x_{k-1} + p_{k}$$

観測値を$yky_{k}$, 観測ノイズを$ok∼N(0,σo2)o_{k} \sim N(0, \sigma_{o}^{2})$とすると
観測方程式は

$$yk=xk+ok y_{k} = x_{k} + o_{k}$$

まとめると、

$${xk=xk−1+pkyk=xk+ok \left\{ \begin{aligned} x_{k} &= x_{k-1} + p_{k} \\ y_{k} &= x_{k} + o_{k} \\ \end{aligned} \right.$$

と表現できます。

ローカル線形トレンドモデル

ここからは少し難しくして、長期的に見た人口の推移をモデル化することを考えてみます。
人口は長期的に見ると右肩上がり(増加傾向) とか右肩下がり(減少傾向) とか傾向がありそうです。
一定の値からランダムに推移するだけのローカルレベルモデルで表現するのは無理がありそうですね(^^;

一定の基準となる値(基準値) に傾向(傾き、トレンド) を取り入れることで上手く表現できそうです。

時刻$kk$での
基準値を$xkx_{k}$, 基準値のプロセスノイズを$pk∼N(0,σp2)p_{k} \sim N(0, \sigma_{p}^{2})$,
傾きを$aka_{k}$, 傾きのプロセスノイズを$qk∼N(0,σq2)q_{k} \sim N(0, \sigma_{q}^{2})$とすると
状態方程式は

$${xk=xk−1+ak−1+pkak=ak+qk \left\{ \begin{aligned} x_{k} &= x_{k-1} + a_{k-1} + p_{k} \\ a_{k} &= a_{k} + q_{k} \\ \end{aligned} \right.$$

観測方程式は（ローカルレベルモデルと同じ変数を使用して）

$$yk=xk+ok y_{k} = x_{k} + o_{k}$$

まとめると、

$${xk=xk−1+ak−1+pkak=ak+qkyk=xk+ok \left\{ \begin{aligned} x_{k} &= x_{k-1} + a_{k-1} + p_{k} \\ a_{k} &= a_{k} + q_{k} \\ y_{k} &= x_{k} + o_{k} \\ \end{aligned} \right.$$

と表現できます。

ローカル季節性モデル

前節では長期的な予測を行うことを想定していましたが、短期的な予測をしたい場合もあると思います。
ここでは飲食店の日別の売上の推移を考えてみることにします。
飲食店の売上の場合、（平日よりも休日の方が売上が上がるなど）曜日ごとに売上の基準値が大きく変わることが想定されます。
このような場合には周期的に変わる季節成分を取り入れる必要がありそうです。

時刻$kk$での
基準値を$xkx_{k}$, 基準値のプロセスノイズを$pk∼N(0,σp2)p_{k} \sim N(0, \sigma_{p}^{2})$,
傾きを$aka_{k}$, 傾きのプロセスノイズを$qk∼N(0,σq2)q_{k} \sim N(0, \sigma_{q}^{2})$,
季節成分を$sks_{k}$, 季節の周期を$mm$, 季節成分のプロセスノイズを$rk∼N(0,σr2)r_{k} \sim N(0, \sigma_{r}^{2})$とすると

状態方程式は

$${xk=xk−1+ak−1+pkak=ak+qksk=−∑i=1m−1sk−i+rk \left\{ \begin{aligned} x_{k} &= x_{k-1} + a_{k-1} + p_{k} \\ a_{k} &= a_{k} + q_{k} \\ s_{k} &= -\sum_{i = 1}^{m - 1} s_{k - i} + r_{k} \\ \end{aligned} \right.$$

観測方程式は全体の基準値とトレンドと季節成分の合計に観測ノイズが乗り

$$yk=xk+sk+ok y_{k} = x_{k} + s_{k} + o_{k}$$

まとめると

$${xk=xk−1+ak−1+pkak=ak+qksk=−∑i=1m−1sk−i+rkyk=xk+sk+ok \left\{ \begin{aligned} x_{k} &= x_{k-1} + a_{k-1} + p_{k} \\ a_{k} &= a_{k} + q_{k} \\ s_{k} &= -\sum_{i = 1}^{m - 1} s_{k - i} + r_{k} \\ y_{k} &= x_{k} + s_{k} + o_{k} \\ \end{aligned} \right.$$

と表現できます。

後書き

久々に数学よりの部分をやりました。
状態空間モデルは柔軟性が高く、個人的に結構便利なモデルだと考えております。
今後予測手法についても書いていこうと思います。

参考

・状態空間モデル
http://www.pleiades-publishing.co.jp/pdf/pdf08.html

https://logics-of-blue.com/状態空間モデルとは/