はじめに
今回は状態空間モデルと呼ばれる非常に広いクラスのモデルを扱います。 状態空間モデルは、ARモデルを代表とする時系列モデルや空間的に隣接している局所構造を表すモデルを記述し、統一的に扱える非常に便利な表現を提供します。
状態空間モデルでのとある設定が、わざわざ固有名詞を持って（別の表現で）モデリングされているケースも多々あることを考えると、状態空間モデルを勉強すれば固有名詞を与えられるほどの重要なモデルを包含した体系を学ぶことができるというわけで、これは抑えておくと非常に良いということです。
まずはじめに状態空間モデルの一般的な表現を学び、それに対して具体的な設定を与えるとどのようなモデルになるのかを見ていくという形式にします。また、状態空間モデルに対する推論手法として非常に重要なカルマンフィルタとパーティクルフィルタにも触れてみます。
状態空間モデル
状態空間モデル
状態空間モデルとは「システム方程式」と「観測方程式」の2本の式から構成されます。
状態 $x _ t$ に関して
$$ x _ {t+1} = f( x _ t) + { \eta} _ t $$
と表される式をシステム方程式と呼びます。これは時刻のインデックス $t$ が1つ進むと状態がどのように変換されるかを $f$ で表現し、更にその時間発展にはノイズ $ \eta _ t$ が混入するということを表しています。 今、 $t$ は時刻のインデックスであると述べましたが、実際のところこれは地理的な情報を表すインデックスであっても構いません。そうだとすると、上記のシステム方程式の意味は、ある地点 $t$ とその1つ隣の場所 $t + 1$ との関係を表していることになります。
そして、上記のようなシステム方程式に従っている状態 $ x _ t$ を我々は直接観測することは無いという問題設定を考えます。我々が観測するのは $ x _ t$ を何らかの方法 $h$ で観測して定まる値にノイズ $ \epsilon _ t$ が混入していると考え
$$ y _ t = h(x _ {t}) + {\bf \epsilon} _ t $$
という観測方程式を定義することにします。このように2本立てになったモデルを状態空間モデルと呼びます。
AR(1) モデル
今、具体例を見るために、$f(x _ t) = a x _ t + c$ と、$h(x _ t) = x _ t$ という関数を状態空間モデルに与えてみることにします。
$$ \begin{align} x _ {t+1} &= a x _ t + c + \eta _ t \\ y _ t &= x _ t + \epsilon _ t \end{align} $$
という二本立てのモデルになります。さて、このモデルは状態 $x$ は時間発展で $a$ 倍され、更に $c$ が加算され、更にノイズ $\eta _ t $ が毎時刻混入する模様です。また、観測値 $y$ は状態 $x$ をそのまま観測するのですが、ここにも観測ノイズ $\epsilon _ t $ が入っているようです。
今、観測方程式に関して実は一切観測時にノイズが混入しないということにしてしまえば、$y _ t = x _ t$ と書くことができます。するともはや、常に状態と観測は同じであって二本立てにする意味はなくなるので、ここから $x _ t$ を排除して、観測データの時系列 $y _ t$ だけを考えることにすると、
$$ y _ {t + 1} = a y _ t + c + \eta _ t $$
という式になります。$\eta _ t$ がホワイトノイズだとすると、これはいわゆるAR(1)モデルにほかなりません。
AR(p) モデル
さて、AR(1)モデルだけを見てる限り、状態空間モデルは単に表現を複雑にしているだけの厄介な存在にしか見えません。しかし、状態空間表現の良い点は、その表記の統一性にあります。
いま、$z _ t = (x _ t, x _ {t-1}) ^ T$ というベクトルを新たな状態として取ってみましょう。すると、次の時刻の状態は $z _ {t+1} = (x _ {t+1}, x _ t) ^ T$ ということになります。また、状態がベクトルになったのでノイズもベクトル $\hat \eta _ t$ として準備してやることにします。そしてこの $z _ t$ と $z _ {t+1}$ の関係を表すシステム方程式は形式的に
$$ z _ {t+1} = f(z _ t) + \hat \eta _ t $$
と書き表されます。状態をどのように定義したかは隠蔽されていますが、状態空間モデルにおけるシステム方程式は常にこの形式になります。さて、このシステム方程式の時間発展に関して $f (z _ t) = A z _ t$ という行列による変換であることにすれば
$$ z _ {t+1} = Az _ t + \hat \eta _ t $$
と表せます。時間発展が線形変換であることを、勝手に定義した $z _ t$ という状態に関して与えてあげただけです。いま、線形変換行列 $A$ を具体的に
$$ A = \begin{pmatrix} a & b \\ 1 & 0 \\ \end{pmatrix} $$
という成分を持つということにしてしまいましょう。また、ノイズベクトルは $\hat \eta _ t = (\eta _ t, 0) ^ T$ であるとして、勝手に定義した $z _ t$ という状態の成分は $(x _ t, x _ {t-1}) ^ T$ であったのでしたから、成分表示でシステム方程式を書いてやると
$$ \begin{pmatrix} x _ {t+1} \\ x _ t \\ \end{pmatrix} \simeq \begin{pmatrix} a & b \\ 1 & 0 \\ \end{pmatrix} \left( \begin{array}{c} x _ t\\ x _ {t-1 } \end{array} \right) + \left( \begin{array}{c} \eta _ t\\ 0 \end{array} \right) $$
と書くことができます。ここで更にシステム方程式を細分化してやると
$$ \begin{align} x _ {t+1} &= a x _ t + b x _ {t-1} + \eta _ t \\ x _ t &= x _ t \end{align} $$
という連立方程式がシステム方程式の正体という事がわかります。見ての通り2個めの式は単なる自明な式であり、意味を持っているのは1つ目の式だけです。状態を $z _ t = (x _ {t}, x _ {t-1} ) ^ T$ と、とりあえずベクトルにしておきながら、実際に表現したかったのは1つ目の式だけで、残りは帳尻合わせに $A$ の成分を決めておきましたということになります。
今観測方程式 $y _ t = h(z _ t) + \epsilon _ t$ の方は、$h(z _ t) = (1, 0) \cdot z _ t$ と $\epsilon _ t = 0$ と与えてあげると、観測は単に $z _ t$ の第一成分だけを見るという意味になり
$$ y _ t = x _ t $$
というのが観測方程式の正体になります。この観測方程式を、システム方程式の成分表示の方に戻してやると
$$ \begin{align} y _ {t+1} & = a y _ t + b y_ {t-1} + \eta _ t \\ y _ t &= y _ t \end{align} $$
と表され、二個目の式は全く同様に自明な式であり、1個目の式が AR(2) モデルを表現していることとなります。
同様に $z _ t = (x _ t, x _ {t-1}, x _ {t-2}, ..., x _ {t - p + 1}) ^ T$ という状態を考えてあげて、帳尻合わせの行列 $A$ を見繕ってやると簡単にAR(p) モデルを仕立て上げることができます。
状態空間モデルに対する推論
さて、ここまで見ても、やはり状態空間モデルは書き方を単に複雑にしている感が否めないのではないでしょうか。 しかし、くどいかもしれませんが、状態空間モデルの良いところは色々なモデルが統一的に扱えるという点であり、中身がどのようなものであろうとも常に状態 $x _ t$ と観測 $y _ t$ に関する2つの式
$$ \begin{align} x _ {t+1} &= f(x _ t) + \eta _ t \\ y _ t &= h(x _ t) + \epsilon _ t \end{align} $$
で表現してやれるのでした。したがって、具体的なモデルとしてどのようなモデルを与えるのかを一旦忘れて、上記の状態空間モデルに関する解法や性質を調べ上げることができれば、その結果を具体的なモデルにいつでも流用ができるというわけです。
当然それはすでに調べられており、状態空間モデルで、観測 $y _ t$ だけしかデータとして得られない場合に、状態 $x _ t$ を推論するアルゴリズムが世の中には提供されています。その1つがよく知られているパーティクルフィルタというわけです。
ここで重要なのはARモデルの例では、状態 $x _ t$ と 観測 $y _ t$ は同一のものである（状態を直接観測できている）という設定のモデルであり、イマイチ状態空間モデルにおいて観測から状態を推論できるという恩恵が感じにくいですが、仮に $y _ t$ は $x _ t$ そのものでなく、間に変換 $h$ が入っており、かつその際にはノイズが混入するという場合（すなわち状態空間モデルの一般的な形式）であっても、パーティクルフィルタは観測 $y _ t$ の系列から $x _ t$ を推論することができるということです。
無論、ノイズの性質はモデルに仮定する必要があり、例えばガウス分布に従っているとか、コーシ分布に従っているとか、そういう情報は与えなければなりません。しかし、そのような仮定の下で、妥当な推論をしてくれるというのはなかなかに強力です。パーティクルフィルタはモンテカルロサンプリングの一種で、それなりに計算量を要しますが、今ではSLAMなど多くのアプリケーションで利用される代表的な手法となっています。
一方で、すでに見たように、システム方程式の時間発展 $f$ が線形変換 $f: x \rightarrow Ax$ であり、同様に観測方程式の $h: x \rightarrow Bx$ という線形変換であるとすると
$$ \begin{align} x _ {t+1} &= Ax _ t + \eta _ t \\ y _ t &= Bx _ t + \epsilon _ t \end{align} $$
という形式で書くことができるようになります。もしも、そんな都合よく線形変換であると言えない場合でも、非線形変換 $f(x)$ をテイラー展開して一次近似を各時刻で実施できるとすれば、$A$ は $f$ のヤコビ行列で置き換えることで一応各時刻局所的には線形システム方程式として表すことができ、観測方程式も同様の処理を実施することで、状態空間モデルをひとまずは上記の形式に落ち着かせることができます。
状態空間モデルが線形変換で構成されており、かつノイズがガウス分布に従っていると仮定できる場合には、観測 $y _ t$ から $x _ t$ を推論する代数計算によって推論するアルゴリズムが利用でき、モンテカルロ法に頼ることが無いため非常に軽い処理で実現できます。その正体がよく知られるカルマンフィルタであり、高速性を要求される組み込み制御の分野で活躍しています（一応補足ですが、状態空間モデルに対してヤコビ行列による一次近似を用いた場合には、特別に拡張カルマンフィルタと呼ばれることが多いです。これはもともと、非線形システムの状態推論をする上でカルマンフィルタの中で一次近似をしているのだという主張から、線形モデルにしか使えないカルマンフィルタを拡張したという意味で使われますが、実際のところ、モデル側を線形近似しているのだというストーリーで考えるとアルゴリズムの中身に変更はありません。）。
まとめ
状態空間モデルは既存のモデルを多く包含した一般的なモデル表現形式を提供します。また、状態空間モデルの形式に持っていくことで、パーティクルフィルタやカルマンフィルタなどの優れた既存の推論アルゴリズムを活用することが可能になります。
多くの場合でパーティクルフィルタが活用されていますが、状態空間モデルのノイズがガウス分布であると考えられるケースではカルマンフィルタ及び拡張カルマンフィルタが計算量の観点から積極的に利用されます。また、無香料カルマンフィルタは二次近似相当の精度を有しています（これは適当な説明で、拡張カルマンフィルタが状態空間モデルの非線形性を近似しているのに対し、無香料カルマンフィルタはモデルを素通りさせておいたあと期待値を二次まで近似している。モデル自体の二次近似はきついが統計量ならなんとか…という考え）。またアンサンブルカルマンフィルタはノイズ項はガウス分布である前提が置かれますが、非線形性をそのまま扱うことができ、一般にパーティクルフィルタより少ない計算量（パーティクルの数を減らせる）で推論が可能です。
