状態空間モデル

今回は状態空間モデルというものを紹介します。

状態空間モデルは統計学や信号処理、制御工学ではお馴染みのモデルで、観測されたデータからモデルの内部の状態を推定することを可能とするものです。代表的なものに、制御や信号処理の分野ではカルマンフィルタがあり、機械学習の分野では隠れマルコフモデルがあります。基本的に扱う対象は時系列データです。

カルマンフィルタは観測が不可能な状態を、観測可能なデータから推測することで、制御対象を適切に制御するための情報を獲得することに用いられています。また隠れマルコフモデルは、観測された音の信号から、それがどんな言葉であったのかを状態量として推定することに用いられています。

状態空間モデル
- 状態空間モデルの例：線形ガウス状態空間モデル
- 線形ガウス状態空間モデルと線形回帰
状態空間モデルでの推定

時刻 $t$ における状態を $\bf x_t$ 、観測データを $\bf y_t$ とすると状態空間モデルは以下の2つのモデルによって表現されます。

${\bf x_t}〜p({\bf x_t|x_{t-1},θ_1})$

${\bf y_t}〜p({\bf y_t|x_t,θ_2})$

$x〜p$ は $x$ が確率分布 $p$ から生起することを表しています。

1つ目の式をシステムモデル、2つ目の式を観測モデルと呼びます。これらを合わせて状態空間モデルと呼びます。考え方は簡単で、状態 $\bf x_t$ は1時刻前の状態 $\bf x_{t-1}$ に依存しながら生起し、観測 $\bf y_t$ は、同時刻の状態 $\bf x_t$ に依存して得られるということです。 $\bf θ$ はいずれも確率分布のパラメータです。初期状態 $t=0$ では ${\bf x_0}〜p({\bf x_0})$ と生起しますがこれは普通省略して書かれています。以下の図が状態空間モデルをグラフィカルに表現したものになります。

$\bf x$ が観測できないような状態量（例えば大型プラントの中のセンサーを取り付けられないような場所の物理量）の場合に、他の物理量をセンサーで観測し、その観測値 $\bf y$ から $\bf x$ を逆算したいような場合に状態空間モデルが使われます。

f:id:s0sem0y:20161024234515p:plain

（一方で現代制御では、状態量の観測ができるような場合に、それに基づいて変化する出力を制御するために状態フィードバックの理論などが提案されました。）

状態空間モデルの例：線形ガウス状態空間モデル

もう一度、状態空間モデルの式に立ち返りましょう。

${\bf x_t}〜p({\bf x_t|x_{t-1},θ_1})$

${\bf y_t}〜p({\bf y_t|x_t,θ_2})$

このとき、システムモデルと観測モデルについて、状態は1時刻前の状態から線形変換を受けるが、そこにはガウスノイズが生じていると考えます。また観測の際にも、現在の状態から観測値が線形変換で得られるのだがその際にもガウスノイズが混入するとしましょう。その際には以下の式で表現できます。

${\bf x_t}={\bf F_t x_{t-1}+v_t}$

${\bf y_t}={\bf H_t x_t + w_t}$

$\bf F_t,H_t$ はそれぞれ時刻 $t$ における線形変換（行列）です。これらの行列は刻一刻と変化するかもしれませんし、実際には時刻によらず常に一定かもしれません（現代制御をやっている方は馴染みが深いでしょう）。

$\bf v_t,w_t$ はそれぞれガウス分布に従うノイズであり、本来の変換から僅かに値がずれる原因だと考えます。ノイズのみが確率的に変化しており、これがガウス分布に従うのならば、結局状態量も観測もガウス分布に従い、確率のモデルとしては平均 $\bf μ$ 、分散共分散行列 $\bf Σ$ の正規分布を $N({\bf μ,Σ})$ として

${\bf x_t}〜N({\bf F_t x_{t-1},Σ_v})$

${\bf y_t}〜N({\bf H_t x_t,Σ_w})$

と表現されることになります。 $\bf Σ_v,Σ_w$ はそれぞれ $\bf v,w$ の共分散行列です。すなわち、ノイズの分散によって、本来変換されるべき値（平均）からずれた値が生起するということです。

線形ガウス状態空間モデルと線形回帰

線形ガウス状態空間モデルの式に立ち返りましょう。

${\bf x_t}={\bf F_t x_{t-1}+v_t}$

${\bf y_t}={\bf H_t x_t + w_t}$

この観測モデル（2式目）は、時刻 $t$ を止めてしまえば、 $\bf y_t$ と $\bf x_t$ の線形回帰モデルになっています。すなわち、時刻が変わる毎に、その線形回帰モデルが刻一刻と変化していく、動的な線形回帰モデルを表現していると解釈できます。

そして、 $x_t$ が1つ前の時刻の状態 $x_{t-1}$ に依存していることを手がかりに（ここも線形回帰になっている）、時間的に変化していくシステムを記述しようということです。

もしもこれを普通の線形回帰と見立てたければ、時間依存を無くしてしまえばいいわけです（例えば $\bf F_t=F_{t-1}$ でかつ、 $\bf x_t$ は時刻に対して独立。いつでもデータ $\bf x$ は個別に手に入って、そのデータ毎に $\bf y$ を線形回帰しているというわけです）。

線形回帰だけでも非常に威力の高いものであることはご存知かと思われますが、線形ガウス状態空間モデルだけでも、それを包含する非常に表現力の高いモデルになっていることがわかります。実際このような比較的単純（線形変換でノイズはガウス分布）なモデルでも応用価値は高く、線形システム理論に基づく制御工学では頻繁に現れますし、生態学などにも応用されています。

状態空間モデルでの推定

状態空間モデルでの推定を行う場合には、以下の2つが考えられます。

1.状態空間モデルで過去、現在、未来の状態量 $\bf x$ あるいは観測値 $\bf y$ を推定

2.あるシステムを状態空間モデルで記述する際のパラメータ $θ$ を推定

「1.」は状態空間モデルが記述できた際に、それを応用しようという場合で、「2.」はあるシステムをこれから記述してみようという場合です。

線形ガウス状態空間モデルでは、パラメータはノイズの分散にすぎませんでした。従って、ノイズの分散を適当に設定しておき、「1.」の用途としてすぐに使うことも考えられますし、システム自体にどんな大きさのノイズが生じているかを知りたければ「2.」を最尤推定などで行うことになります。

線形ガウス状態空間モデルで表現できなさそうなものに対してはどうでしょうか。「1.」では条件付き確率分布を計算することが（ガウス分布のような綺麗な式でなければ）困難であるために、実は簡単に利用することはできません。また最尤推定をしようにも、その計算自体が困難であるため、システムを記述しようという試み「2.」も困難を極めます。これに対する具体的な解決手法が様々提案されており、研究の主題でもあるというわけです。

カルマンフィルタ

状態量と観測値が連続値である場合に、線形ガウス状態空間モデルを用いて現在の観測値から現在の状態量を推定する手法として、カルマンフィルタがあります。線形ガウス状態空間モデルは、全ての確率分布がガウス分布に従うために、解析的に計算を行うことができるため、上記で述べたような困難は訪れません。状態量の推定もパラメータの推定（最尤推定）も厳密解が得られます。

問題はむしろ対象が本当に線形モデルなのかということであり、非線形なモデルを線形近似する手法が提案されています。

大規模なシステムの全ての状態量をセンサーで観測できるとは限りません。限られたセンサーから、観測できない他の状態量を推定することができれば、その推定値を使い適切にシステムの制御を行うことに役立てられます。通常このようにセンサーなどで観測できない情報を得る計算装置を「オブザーバ」と呼び、オブザーバの予測値を使い制御工学にて応用されています。

カルマンスムーザー

状態量と観測値が連続値である場合に、線形がガウス状態空間モデルを用いて現在の観測値からそれまで（過去）の状態量を推定する手法として、カルマンスムーザーがあります。

これも状態量の推定が厳密に行えますし、モデルのパラメータ推定はカルマンフィルタと同様です。（そもそもカルマンフィルタとカルマンスムーザーは予測したい状態量の時刻が異なっているだけ）

隠れマルコフモデル

状態量と観測値が離散値である場合に、状態空間モデルを用いて観測値から状態量を推定する手法に隠れマルコフモデルがあります。

隠れマルコフモデルはノイズがガウス分布であるという仮定はありませんが、状態が線形変換を受けているという仮定の下で、パラメータをEMアルゴリズムによって推定することができます。通常隠れマルコフモデルに用いるEMアルゴリズムは「Baum-Welchアルゴリズム」と呼ばれています。

観測されたデータ列から、状態量の遷移がどのように行われたのかもダイナミックプログラミングによって推定が可能で、通常隠れマルコフモデルに用いるダイナミックプログラミングをビタビアルゴリズムと呼びます。

例えば音声信号を観測データとして、その音声に対応する言葉を状態量としておけば、隠れマルコフモデルは音声認識装置として働くことになります。

粒子フィルタ

一般的な状態空間モデルに対して、観測値から未来の状態、現在の状態、過去の状態を予測する手法として粒子フィルタがあります。粒子フィルタはサンプリング手法の一種で、多数のサンプリングを行い、その標本分布から確率分布のベイズ更新で推定するため、計算量は大きくなります。現在は計算機の性能が向上しているため応用されるようになりました。