非線形状態空間モデルの粒子フィルタ

非線形状態空間モデル †

標準形

xt+1=f(xt,wt):状態推移式
yt=h(xt,et)　：観測式
wtはシステムノイズ、etは観測ノイズで相互に独立とする。

観測結果{y0,y1,...yt}をYt={yt}で表わし、観測後の条件付き密度関数は次のように表わされる。(ベイズの定理より)

p(xt|Yt)=p(yt|xt)p(xt|Yt-1)/p(yt|Yt-1)
p(yt|Yt-1)=∫p(yt|xt)P(xt|Yt-1)dxt

また、状態xt+1のYtによる事後確率密度関数は、次のように表わされる。

p(xt+1|Yt)=∫p(xt+1|xt)P(xt|Yt)dxt

初期分布は次式で表わす。

P(x0|Y-1)=P(x0)

線形のカルマンフィルタでは、wt,etをガウス分布と仮定することで、期待値ベクトルと共分散行列を用いて、上記の逐次式を解析的な解(逐次形)表現することができた。詳細：カルマンフィルタを参照ください。

最近のシミュレーションによる方法として粒子フィルタが注目されている。より良い推定値を少ない計算量で求めるために、状態ベクトルxtを２つ分けて、次のように表わす。

xt=(xt(k),Xt(n))
xt(k)は、条件付きの線形ダイナミクスをもつ状態変数である。
Xt(n)は、非線形状態ベクトルである。

前者は、線形部分であり、有限次元の最適フィルタと推定できる。後者は粒子フィルタを持ちいて推定する。

このような方法はRao-Blackwellizationとも言われ、ベイズフィルタの研究（2001)で行われている。以下では、線形部分構造をもつ非線形状態空間モデルの一般的な方法を紹介する。

技術的には、以下の応用分野が考えられる。

標準的な粒子フィルタをまず紹介する。

粒子フィルタは、逐次的に事後密度関数を推定する問題を解いて、近似解を与える。以下では、大文字Xt={xt}で、t期までの状態の集合を表わす。Yt={yt}も同様である。

最も関心をもつのは、フィルタリング事後密度　P(xt|Yt)であるが、この密度は、粒子と呼ばれるたくさんのサンプルの集合{xt(i)}i=1,Nで近似される。（そこで粒子フィルタと呼ぶ）。

P*(xt|Yt)=Σqt*(i)δ(xt-xt(i))
qt*(i)は、正規化された重要度の重み

このqt(i)は、次式で修正(update)される。

qt(i)=P(yt|xt(i))qt-1(i)

この式の意味は、大きい尤度をもつサンプルにより大きい重要度が与えられるということである。

完全な粒子フィルタのアルゴリズムは、次式で与えられる。詳細な導出は、参考(2),(3)を見ること。

2.重要度の推定と正規化を行う

qt(i)=P(yt|　xt|t-1(i))　　:重要度の推定
qt*(i)=qt(i)/[Σ qt(i)]　：正規化

次のモデルを考える

xt+1(n)=ft(xt(n)) +wt(n) :非線形部分
Xt+1(k)=Atxt(K) + Wt(k)
yt=h(xt(n)) +et
Wt(k)は平均0、共分散Qｔ(k)、Wt(n)は平均0、共分散Qｔ(n)のガウス型システムノイズ
etは、平均0、共分散Rtのガウス型の観測ノイズ

ベイズの定理より

P(xt(k),Xt(n)|Yt)=P(xt(K)|Xt(n),Yt)・P(Xt(n)|Yt)

右辺の前の項は、線形部分構造の条件付き確率であるので、カルマンフィルターによって計算できる。後の項は、非線形であるので、粒子フィルタで近似計算ができる。

それぞれの部分にわけて、計算できるので便利である。