状態空間モデルの推定方法の分類

新規作成：2016年2月28日
最終更新：2016年2月28日

カルマンフィルタに拡張カルマンフィルタ、粒子フィルタ、そんでもってギブスサンプラーにハミルトニアンモンテカルロ法。
呪文みたいですが、これらはすべて、状態空間モデルの推定に使われるアルゴリズムです。

細かいアルゴリズムの話はここではしません。そこは専門の書籍をご参照いただければと思います。
ただ、呪文のままでは、少々やりづらいですね。
ここではおおざっぱな分類について、ごく簡単に解説します。

その他の状態空間モデル関連のページはこちら。

 


広告

 

目次

1．分類1：どのようなデータに適用可能か
2．分類2：フィルタ系とMCMC系で分ける
3．フィルタ系とMCMC系の比較
　逐次処理とバッチ処理
　パラメタ推定をやるかやらないか
　　状態推定とパラメタ推定の違い
　　フィルタ系でのパラメタ推定の流れ
　　MCMC系でのパラメタ推定の流れ
4．フィルタ系とMCMC系の使い分け

 

1．分類1：どのようなデータに適用可能か

最初に、適用可能なデータによる分類の仕方を紹介します。
よく見られる分類なんですが、誤解を生みやすいかなと思うので、気を付けてください。

■分類1 単純なモデルしか推定できない手法
・カルマンフィルタ
・拡張カルマンフィルタ

■分類2　非正規、非線形でも適用できる手法
・粒子フィルタ
・ギブスサンプラー
・ハミルトニアンモンテカルロ法

カルマンフィルタは単純で、粒子フィルタやハミルトニアンモンテカルロ法は複雑なことができる、という分類です。

計算を単純にするためには、仮定を置きます。仮定を置けば、「仮定と異なる状況」は無視して計算を進めることができるので、計算が単純になります。
よく使われる仮定は、データが正規分布に従うという仮定と、線形であるという仮定です。

正規、線形である状態空間モデルを推定するにはカルマンフィルタが使われます。
正規分布だけれども、非線形であるデータに対して適用したい場合は拡張カルマンフィルタを使います。
非正規、非線形の状態空間モデルを推定する場合には、粒子フィルタやギブスサンプラー、ハミルトニアンモンテカルロ法を使います。
計算は難しくなる（計算量が多くなる）んですが、様々なデータに対して適用が可能です。

これが分類方法その1。適用可能なデータによる分類です。

ただし、粒子フィルタとギブスサンプラーはかなり違っているので、これを同じだとまとめてしまうのは良くないかなと思います。
あくまでも、計算方法は完全に無視して、「適用可能なデータ」のみを見た分類だとご理解ください。

 

2．分類2：フィルタ系とMCMC系で分ける

状態空間モデルの推定方法には「○○フィルタ」と呼ばれるものと、そうでないものがあります。これで分けるのが2つ目の分類です。
こっちのほうが、少し正確な分類の仕方であるように思います。

■フィルタ系
・カルマンフィルタ
・拡張カルマンフィルタ
・粒子フィルタ

■MCMC系
・ギブスサンプラー
・ハミルトニアンモンテカルロ法

名前を見れば判別できるので簡単ですね。
ちなみに「フィルタ系」や「MCMC系」とは私の造語なので気を付けてください。
分け方だけを参考にしていただければと思います。

次からは、○○フィルタとMCMC系の違いについてみていきます。

 

3．フィルタ系とMCMC系の比較

○○フィルタと呼ばれるものとMCMC系の比較の表を作ってみました。
ここでは、分類1で見た「どのような仮定を置くことが必要か」という観点は無視していることに注意してください。

また、スペースの都合上、専門用語を使っています。知らない用語があった場合は、別途文面で説明しているので、そこを見てからまた表に戻られると良いかと思います。

	フィルタ系	MCMC系
1	逐次処理	バッチ処理
1-2	一回の計算量は少なめ	毎回、多くの計算を行わなければならない
1-3	フィルタリングをした後でスムージングをする	フィルタリングなしで、最初からスムージングしている
2	状態推定のみを行う	状態推定だけでなく、モデルのパラメタも同時に推定する
2-2	別途、最尤法によるパラメタの推定が必要	パラメタ推定だけを別に行う必要がない

大きく2つの観点があります。
計算の仕方（逐次処理かバッチ処理か）という観点と、私たちがやらなければいけない作業（パラメタ推定を別途やるかやらないか）がどう変わるかという観点です。

逐次処理とバッチ処理

最も大きい違いは、逐次処理かバッチ処理かの違いでしょう。
フィルタ系は逐次処理、MCMC系はバッチ処理です。

バッチ処理とは「データをまとめていっぺんに使って、まとめていっぺんに計算すること」くらいの意味合いだととらえてください。
例えば、データが増えるたびに1を足すという処理があったとします。
データが来るたびに1を足すのが逐次処理。
データが5つ来たならば、「1を足す処理を最後にまとめて5回やる」のがバッチ処理です。

当然ですが、逐次処理のほうが、一回一回の計算量は少なくなります。
一方、バッチ処理の場合は、「1を足す処理を最後にまとめて5回やる」という処理を、時と場合によっては「最後に5を足す」と言い換えることもできます。このようなやり替えができるのがバッチ処理の良いところです。

この計算方式の違いにより、フィルタリング、スムージングという処理に大きな違いが出てきます。
まずは、フィルタ系の計算方法を見ていきます。1期前のデータを使って、次の期間のデータを予測する状態空間モデルなのだと思っておいてください。例えばローカルレベルモデルなどです。
1．1期前の状態推定値を使って、次の状態を予測する
2．予測された状態推定値を、当期のデータを使って修正する→これがフィルタリング

例を挙げます。2000年の状態が10だとわかっていたとしましょう。湖の中にいる魚の数が10尾だと推定されている、とイメージしてください。
　　　　　2000年　　2001年　　2002年　　2003年　　2004年
状態　　　 10
データ

 
予測をすると、2001年の状態推定値も10だと予測されました
　　　　　2000年　　2001年　　2002年　　2003年　　2004年
状態　　　 10　　 →　10
データ

 
しかし、2001年に潜水艦を使って魚の数を調べてみると、5尾しかいないようでした。
　　　　　2000年　　2001年　　2002年　　2003年　　2004年
状態　　　 10　　 →　10
データ　　　　　　　　5

 
というわけで、2001年の状態推定値を補正します。ちょっと減って、7になりました。
　　　　　2000年　　2001年　　2002年　　2003年　　2004年
状態　　　 10　　 →　7
データ　　　　　　　　5

これがフィルタリングです。
潜水艦で調べても、数え漏らしがあるでしょう。なので、観測結果と状態の推定値は一致していません。いい感じの塩梅で補正をしてくれるのが、フィルタリングの良いところです。

 
次、2002年の状態を推定するときには、2001年の状態がわかっていればよいです。2000年のデータや状態は全く不要です。
　　　　　2001年　　2002年　　2003年　　2004年
状態　　　 7
データ

あとはこの繰り返し。
データが来るたびに、「1期前の状態推定値から、1期後の状態推定値を予測する」という処理と「当期のデータを使って、当期の状態推定値を修正する」という処理を逐次的に行っています。
これがフィルタ系でやっていることです（厳密には「当期のデータを使って、当期の状態推定値を修正する」の部分だけがフィルタリングです。「予測→フィルタリング」を繰り返し行うのがフィルタ系ということです）。

○○フィルタと呼ばれるアルゴリズムはこれをやっています。なので、逐次的に計算が行われます。

しかし、このやり方だと、状態の推定に使われているのが「1期前の状態推定値」と「当期のデータ」だけです。
データを長い期間とっていたとしても、使われるのはこの2つだけ。
これってちょっともったいなくないですか？

というわけで、フィルタリングをした後に、もう一度スムージングをします。
下記のようにデータが得られていたとしましょう。
　　　　　2000年　　2001年　　2002年　　2003年　　2004年
状態　　　 10　　 　　？　　　　　　　　　　　　　　
データ　　　　　　　　5　　　　　12　　　　　15　　　　13

 
このとき、2002年以降は、データが大きな数値をとっています。
ということは2001年の状態も、ちょっと大きめにとっておくのがよさそうな気がしますね。
　　　　　2000年　　2001年　　2002年　　2003年　　2004年
状態　　　 10　　 　　9　　　　　　　　　　　　　　
データ　　　　　　　　5　　　　　12　　　　　15　　　　13

フィルタリングだと2001年の状態推定値は7でしたが、ちょっと大きめに9となりました。
このように、データをまとめて使って、状態を推定するのがスムージングです。
フィルタリングした後にわざわざスムージングをするのは、「データすべてを使って状態を推定したいから」です。

で、MCMC系は、バッチ処理なのでした。
ということは、最初から、データすべてを使って状態を推定しています。
ということなので、MCMC系のアルゴリズムは、最初からスムージングをしているようなものなのだということです。

※ちょっとざっくり書いてしまったので、ご指摘があれば、受け付けます。

パラメタ推定をやるかやらないか

フィルタ系は、状態推定とは別に、パラメタの推定も必要です。
まずは、状態推定とパラメタ推定の違いから説明します。

状態推定とパラメタ推定の違い

状態とは、湖の中にいる魚の数や、景気の動向を無視した企業体力などを指します。
ところで、状態ってどうやって推定しましょうか。
先のフィルタリングの説明では、データを使って「いい感じの塩梅で」補正するという大変適当な書き方をしていました。
「いい感じの塩梅」を決めるのがパラメタです。良いパラメタを使うといい感じの塩梅になります。ダメなパラメタを使うと、ダメな感じで補正されます。

もうちょっとまともに説明します。
状態空間モデルの場合、観測誤差の分散の大きさが一つのパラメタになります。
観測誤差が大きければ、観測結果、先の例では潜水艦で調べた魚の数などが信用できなくなります。乗組員の目が悪かったのでしょう。
もし、観測誤差が大きければ「データを使っても、あまり補正をかけないようにする」ことが最善となります。

例えば、こんなデータがあったとします。
　　　　　2000年　　2001年
状態　　　 10　　 →　10
データ　　　　　　　　5

観測誤差が大きければ、たとえデータが「５」であったとしても、2001年の状態推定値は「10」のままになります。観測誤差が大きいんだから、データはあまり信用しないんですね。
逆に、観測誤差が小さければ、極端な話、観測誤差がないのだとしたら、2001年の状態推定値は、2001年のデータと同じ「5」となります。

というわけで、状態を推定するのに必要となるパラメタがあって、そのパラメタを決めなければならないという話でした。

フィルタ系でのパラメタ推定の流れ

ローカルレベルモデルなどでフィルタリングをするのに必要となるものは、以下の通りです。
・1期前の状態推定値
・当期のデータ
・状態を推定するのにつかわれるパラメタ

パラメタは、引数として事前に与えられています。
パラメタを与えなければ、フィルタリングできません。

パラメタを推定する流れは以下の通りです。
1．テキトーなパラメタを使ってモデルを作る
2．テキトーなパラメタを使ったモデルの結果の「良さ」を測る
3．良さが大きくになるようにパラメタを少しずつ修正していく
4．良さが最大のパラメタを採用

ちなみに「良さ」の基準には尤度、あるいはその対数をとった対数尤度が使用されます。
尤度が最も大きくなるパラメタを採用することを最尤推定法と呼びます。

dlmの使いかたという記事において、状態空間モデルの推定には以下の4つのステップを踏むと書きました。
カルマンフィルタの例ですが、ほかも大体一緒です。
Step1　状態空間モデルの「型」を決める
Step2　その「型」にいれるパラメタを推定する　←今ここ
Step3　Step2で推定されたパラメタを「型」に入れてカルマンフィルタを回す
Step4　カルマンフィルタの結果を使ってスムージングする

フィルタ系では、事前にパラメタを推定してから、フィルタリングを行います。

ちなみに、パラメタを推定するためには、過去、現在、未来、すべてのデータが使われます。
当期のデータだけを使うフィルタリングとは大きく異なる点に注意してください。

MCMC系でのパラメタ推定の流れ

MCMC系では、状態推定とパラメタ推定を分けません。
MCMC系はバッチ処理です。みんなまとめていっぺんにやってしまいます。

というわけで、MCMC系では、パラメタの推定と状態の推定を分けることなくまとめてやってしまいます。
一発で済む代わりに、その一発に時間がかかるのがMCMCです。

 

4．フィルタ系とMCMC系の使い分け

まずはフィルタ系を使うと良いだろうという局面を紹介します。
・手早く状態を推定する必要があるとき

フィルタ系は、パラメタ推定と状態推定を完全に別に分けて計算することができます。
フィルタリングするだけですと、1期前の状態推定値と当期のデータ「だけ」で状態が推定できてしまいます。
必要とする情報も少ないのです。

状態空間モデルの応用例として真っ先に上がるのが、アポロ計画ですね。
あれは拡張カルマンフィルタを使っていたと思うんですが、早く計算を行うのが必須の要件でした。
具体的には、「ロケットの居場所」という状態を、各種のセンサー（観測値）から推定するという問題です。
センサーの数値を入力すると、即座にロケットの場所を計算し、正しい軌道にのせるように修正します。

この場合、観測誤差の大きさなどのパラメタは、過去の実験結果などをもとにして、事前に推定しておきます。
その「パラメタの推定結果」は更新することなく使います。これを新しい観測値にそのまま適用して「状態の推定値」のみを更新します。
逐次的に素早く計算ができます。

ほかにも、画像処理などで「カメラなどのセンサーから得られた観測値」をもとにして「障害物との距離という状態」を逐次的に推定したい場合などにも、フィルタ系のアルゴリズムはその効果を発揮するでしょう。
ですので、フィルタ系のアルゴリズムは制御系、工学系の教科書や論文でよく出てきます。カルマンフィルタの本が統計学の棚ではなく制御工学の棚に置いてあるのはこれが理由です。

一方のMCMC系に関しては、データからの知識発見というところでより効果を発揮するのではないかなと思います。
私はアザラシの個体数を状態空間モデルで推定したことがあります（論文はこちら）。
天候が悪くて観測できない、あるいは観測結果がぶれてしまうという場合でも、アザラシの個体数を評価することができるので、資源の保全などに役立てることができます。
これを手早く逐次的に計算する価値はあまりありませんね。

マーケティングなどで、お客様の購買データを分析して、お店においてどんな商品の配置にするかを決める、といった場合でも、逐次処理にするメリットはなさそうです（毎秒商品の陳列を自動的に決めることができる未来なお店だと別でしょうが。あるいはWebでの購買の最適化などにはフィルタ系のほうが良いかもしれません）。

こういった場合だと、パラメタ推定もまとめてやってくれるMCMC系を使ったほうが便利かなと思います。
ただし、カルマンフィルタだと、dlmパッケージという大変便利な計算ツールがあります。逐次処理でなくても、楽に計算したい場合はカルマンフィルターを使うと良いかもしれません。

最後に、こういう分類においてはいつも付きまとうのですが、「分け方」には当事者のバックグラウンドや個性などが大きく関与してきます。
異論、反論、ご意見ご感想は、メール（logics.of.blue★gmail.com）や記事のコメントなどでいただければ幸いです。

 

おまけ：計算手法別の参考書

おまけとして、計算手法別の参考書を上げておきます。

フィルタ系の参考書

カルマンフィルタの基礎   カルマンフィルタ、拡張カルマンフィルタ、UKF（これも非線形で使えるように拡張されたカルマンフィルタ）の解説が載っています。

予測にいかす統計モデリングの基本―ベイズ統計入門から応用まで (KS理工学専門書)   こちらは粒子フィルタの解説が載っています。

MCMC系の参考書

基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門   ベイズ推定について書かれた本です。ハミルトニアンモンテカルロ法の解説が載っています。 巻末の付録にStanというソフトの使い方も書かれています。

岩波データサイエンス Vol.1 特集「ベイズ推論とMCMCのフリーソフト」   こちらはソフトの使い方が詳しく載っています。Stanだけでなく、ギブスサンプラーができるBUGS言語の紹介もあります。 状態空間モデルの計算例も載っています。

 


広告