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

追加された行はこの色です。
削除された行はこの色です。
非線形状態空間モデルの粒子フィルタへ行く。
#freeze
*非線形状態空間モデル [#nfe35c91]
-標準形
 xt+1=f(xt,wt):状態推移式
 yt=h(xt,et)　：観測式
 wtはシステムノイズ、etは観測ノイズで相互に独立とする。
*事後確率密度関数の逐次表現 [#wd91f907]
観測結果{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をガウス分布と仮定することで、期待値ベクトルと共分散行列を用いて、上記の逐次式を解析的な解(逐次形)表現することができた。詳細：カルマンフィルタを参照ください。
*粒子フィルタ [#rb2832f9]
最近のシミュレーションによる方法として粒子フィルタが注目されている。より良い推定値を少ない計算量で求めるために、状態ベクトルxtを２つ分けて、次のように表わす。
 xt=(xt(k),Xt(n))
 xt(k)は、条件付きの線形ダイナミクスをもつ状態変数である。
 Xt(n)は、非線形状態ベクトルである。
前者は、線形部分であり、有限次元の最適フィルタと推定できる。後者は粒子フィルタを持ちいて推定する。
--このような方法はRao-Blackwellizationとも言われ、ベイズフィルタの研究（2001)で行われている。
以下では、線形部分構造をもつ非線形状態空間モデルの一般的な方法を紹介する。

技術的には、以下の応用分野が考えられる。
-位置追跡システム、ナビゲーション
-金融時系列分析
*粒子フィルタの紹介 [#haaf67a3]
標準的な粒子フィルタをまず紹介する。
-粒子フィルタは、逐次的に事後密度関数を推定する問題を解いて、近似解を与える。
以下では、大文字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)を見ること。
-1.初期化：粒子i=1~Nの事前確率を初期化、X0|-1(i)を与える。P(x0) の近似。
-2.重要度の推定と正規化を行う
 qt(i)=P(yt|　xt|t-1(i))　　:重要度の推定
 qt*(i)=qt(i)/[Σ qt(i)]　：正規化
-３．観測による修正(update)：
 Prob{xt|t(i)=xt|t-1(j)｝=qt(j)
-４.1期予測修正:
 xt+1|t(i)=P(xt+1|t　|　xt|t(i))
-５．t-->t+1として、ステップ２に戻る。
*非線形状態空間モデルの推定 [#z457c637]
次のモデルを考える
 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)
右辺の前の項は、線形部分構造の条件付き確率であるので、カルマンフィルターによって計算できる。後の項は、非線形であるので、粒子フィルタで近似計算ができる。

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

-補題１：線形部分
線形部分のフィルター値xt|t(k)と１期予測値xt+1|tの確率密度は、
 P(xt(k)|Xt,Yt)=N(x*t|t(k),Pt|t)
 P(xt+1(k)|xt+1,Yt)=N(x*t+1,pt+1|t)
のように、正規分布N(平均,共分散)の形で与えられ、そのアルゴリズムは、下記のカルマンフィルタで求められる。
--フィルタリング値
 x*t|t(k)=x*t|t-1(k)+Kt(yt-ht-Ctx*t|t-1)
 Pt|t-1=Pt|t-1 - KtCtPt|t-1
 St=CtPt|t-1Ct'+Rt
 Kt=Pt|t-1Ct'(St)^(-1)
--1期予測
 x*t+1|t(k)=Atx*t|t(k)
 Pt+1|t=AtPt|t(At')+Qt(k) 
-補題2
事後密度関数は
 P(xt(n)|Yt)={P(yt|Xt(n),Yt-1)p(xt(n)|Xt(n),Yt-1)/P(yt|Yt-1)}p(xt-1(n)|yt-1)
上式であたえられる。p(xt-1(n)|yt-1)は、前の計算から与えられるがこの右辺の残り２つの確率は次式で計算できる。
 p(yt|Xt(n),Yt-1)=N{ht-Ctx*t|t-1(k),CtPt|t-1Ct'+Rt)
 p(xt+1(n)|xt(n),Yt)=N(ft(n),Qt(n)
-以上の２つの補題から、全体のアルゴリズムを構成できる。

--詳細は、参考(4)の研究を見てください。

*参考 [#x5837960]
-(1)[[Arnaud Doucet, Neil J. Gordon, and Vikram Krishnamurthy;Particle Filters for State Estimation of Jump Markov Linear Systems,IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 49, NO. 3, MARCH 2001>http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.20.6159&rep=rep1&type=pdf]]
-(2)[[Arnaud Doucet；On Sequential Simulation-Based Methods for Bayesian Filtering,(1998)>http://stargate.uwaterloo.ca/~jzelek/teaching/SYDE777/tr310.pdf]]
-(3)[[On sequential Monte Carlo sampling methods for Bayesian filtering,Statistics and Computing (2000) 10, 197–208>http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.12.5875&rep=rep1&type=pdf]]
-（４）[[Thomas Sch¨on, Fredrik Gustafsson, Member, IEEE, and Per-Johan Nordlund;Marginalized Particle Filters for Mixed Linear/Nonlinear State-space Models>http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.4.2903&rep=rep1&type=pdf]]
非線形状態空間モデルの粒子フィルタ の変更点

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