粒子フィルタ

出典: フリー百科事典『ウィキペディア(Wikipedia)』
移動: 案内検索

粒子フィルタ、または逐次モンテカルロ法 (Sequential Monte Carlo: SMC)とは、シミュレーションに基づく複雑なモデルの推定法である。

この手法はふつうベイズモデルを推定するのに用いられ、バッチ処理であるマルコフ連鎖モンテカルロ法 (MCMC) の逐次 (オンライン) 版である。またこの手法は重点サンプリング法にも似たところがある。 うまく設計すると、粒子フィルタはMCMCよりも高速である。拡張カルマンフィルタ無香カルマンフィルタ (Unscented カルマンフィルタ) に対して、十分なサンプル点があればベイズ最適推定に近付くためにより精度が高くなることから、これらの代わりに用いられることがある。手法を組み合わせ、カルマンフィルタを粒子フィルタの提案分布として使うこともできる。

目的[編集]

粒子フィルタの目的は、隠れたパラメータ列 x_k (k=0,1,2,3, \ldots) を観測値 y_k (k=0,1,2,3, \ldots) のみから推定することである。ベイズ推定値 x_k は一つ過去の確率分布 p(x_k|y_0,y_1,\ldots,y_k) から得られるのに対し、マルコフ連鎖法では過去の全てを含む確率分布 p(x_0,x_1,\ldots,x_k|y_0,y_1,\ldots,y_k) より求められる.

モデル[編集]

粒子フィルタではx_kと観測値 y_k は以下のように表せるとする。

ただし初期確率分布 p(x_0) を持つものとする。

  • 観測データ列 y_0, y_1, \ldotsx_0, x_1, \ldots が既知であるという条件の下で独立である。
換言すると、それぞれの y_kx_k のみに従属する。
y_k|x_k \sim p_{y|x}(y|x_k)

このようなものの一例として次のものがある。

x_k = f(x_{k-1}) + v_k \,
y_k = h(x_k) + w_k \,

ただし v_kw_k も異なるkについて互いに独立であり、かつkによらず同じ確率密度関数にしたがうものとする。また、f(\cdot)h(\cdot) は既知の関数である。 この2つの方程式は状態方程式とみなすことができ、カルマンフィルタの状態方程式と類似している。もし関数 f(\cdot)h(\cdot) が線形かつ v_k および w_kガウス分布 だとすると、カルマンフィルタによってベイズ推定と厳密に一致する解が得られる。そうでない場合はカルマンフィルタは1次近似に過ぎない。粒子フィルタも近似には変わりないが、十分な数の粒子があれば高い精度が得られる。

モンテカルロ近似[編集]

粒子法は他のサンプルリング法と同様に、フィルタリング確率分布p(x_k|y_0,\dots,y_k)を近似する点列を生成する。したがって P 個のサンプル点があれば、フィルタリング確率分布による期待値は

\int f(x_k)p(x_k|y_0,\dots,y_k)dx_k\approx\frac1P\sum_{L=1}^Pf(x_k^{(L)})

によって近似される。そして通常のモンテカルロ法と同様に f(\cdot)を適切に調整することで, 近似の程度に応じた次数までの モーメント を得ることができる。


Sampling Importance Resampling (SIR)[編集]

Sampling importance resampling (SIR) アルゴリズムは、(Gordon et al. 1993)による本来の粒子フィルタであるが、よく使われる粒子フィルタアルゴリズムである。ここではフィルタリング確率分布 p(x_k|y_0,\ldots,y_k)P個の重みつき粒子によって近似する。

\{(w^{(L)}_k,x^{(L)}_k)~:~L=1,\ldots,P\}.

重み w^{(L)}_k は以後の \sum_{L=1}^P w^{(L)}_k = 1 である粒子の相対事後分布(密度)の近似となっており、 \sum_{L=1}^P w^{(L)}_k = 1を満たす。

SIRは重点サンプリング の逐次版 (つまり、再帰的) と言える。 重点サンプリングにあるように、関数 f(\cdot) の推定値は重みつき平均

    \int f(x_k) p(x_k|y_0,\dots,y_k) dx_k \approx \sum_{L=1}^P w^{(L)} f(x_k^{(L)})

で近似できる。

粒子の個数が有限である場合、このアルゴリズムの性能は提案分布\pi(x_k|x_{0:k-1},y_{0:k})の選択に依存する。 最適な提案分布は目的分布

    \pi(x_k|x_{0:k-1},y_{0:k}) = p(x_k|x_{k-1},y_{k})

である。

しかしながら事前遷移確率

    \pi(x_k|x_{0:k-1},y_{0:k}) = p(x_k|x_{k-1})

を提案分布として用いることが多い。 なぜならそこからは容易に粒子をサンプルすることができるし、その後に重みを計算することもできるからである。

提案分布として事前遷移確率を用いるSIRフィルタは一般にブートストラップフィルタ・コンデンセーションアルゴリズムとして知られている。


唯一つを除いた全ての重みがゼロに近付くこと ― アルゴリズムの縮退 ― を防ぐために再サンプルが行なわれる。 このアルゴリズムの性能は再サンプリング方式の選びかたにもかかっている。 北川 (1996) によって提案された階層化再サンプル法は分散の意味で最適である。


SIR法の1ステップは以下の様になる。

1) L=1,\ldots,P について, 提案分布から粒子をサンプルする

x^{(L)}_k \sim \pi(x_k|x^{(L)}_{0:k-1},y_{0:k})
2) L=1,\ldots,P について、重みを更新し、正規化定数を得る。

        \hat{w}^{(L)}_k = w^{(L)}_{k-1} \frac{p(y_k|x^{(L)}_k) p(x^{(L)}_k|x^{(L)}_{k-1})} {\pi(x_k^{(L)}|x^{(L)}_{0:k-1},y_{0:k})}

このとき\pi(x_k^{(L)}|x^{(L)}_{0:k-1},y_{0:k}) = p(x^{(L)}_k|x^{(L)}_{k-1}) ならば更新式は


            \hat{w}^{(L)}_k = w^{(L)}_{k-1} p(y_k|x^{(L)}_k),

のように簡単になる。

3) L=1,\ldots,Pについて、正規化された重みを計算する。

w^{(L)}_k = \frac{\hat{w}^{(L)}_k}{\sum_{J=1}^P \hat{w}^{(J)}_k}
4) 有効粒子数の推定量を計算する。

\hat{N}_\mathit{eff} = \frac{1}{\sum_{L=1}^P\left(w^{(L)}_k\right)^2}


5) もし有効粒子数が与えられた閾値より少なかったら \hat{N}_\mathit{eff} < N_{thr} 再サンプルを実行する。
a) 現在の粒子の集合から、重みに比例した確率で P 個の粒子を描く。現在の粒子の集合を新しい粒子の集合で置き換える。
b) L=1,\ldots,P について、 w^{(L)}_k = 1/P とする。


Sequential Importance Resampling 等の用語は SIR フィルターの意味で用いられることがある。


逐次的 Importance Sampling (SIS)[編集]

  • 再サンプルが無い点を除いて SIR と同様である。


直接法[編集]

直接法は他の粒子フィルタに比べて簡単で、合成と棄却を利用する。 k 番目の一つのサンプル xp_{x_k|y_{1:k}}(x|y_{1:k}) から生成するのに以下の手順を踏む。

1) n = 1とする。
2) \{1,..., P\}上の一様分布からLを生成する
3) テスト粒子 \hat{x} を確率分布 p_{x_k|x_{k-1}}(x|x_{k-1|k-1}^{(L)}) から生成する
4) \hat{y} の確率を \hat{x} を用いて p_{y|x}(y_k|\hat{x}) から生成する。ただし、y_k は観測値である
5) 別の数u[0, m_k]上の一様分布からを生成する。ただし、m_k = \sup_x p_{y|x}(y_k|x)
6) u\hat{y} を比較
6a) u が大きければ step 2 以降を繰り返す
6b) u が小さかったら \hat{x}x{k|k}^{(p)} として保存して n を1増やす
7) n > P ならば終了

目的はkにおける P 個の粒子を、k-1 だけから生成することである。 これにはマルコフ方程式が x_{k-1} だけから x_k を生成できるように記述できなければならない。 このアルゴリズムは k における粒子を生成するためにk-1におけるP個の粒子の合成を利用しており、kにおいてP個の粒子が生成されるまでステップ 2--6 以降を繰り返す。

これは x を2次元配列とみなすと、より視覚的に理解できる。 一つの次元は k であり、もう一つの次元は粒子番号Lである。 例えば x(k,L)k におけるL番目の粒子であり、x_k^{(L)} と書ける (上記のアルゴリズムの記述に用いた通り)。 ステップ 3 は k-1 においてランダムに選ばれた粒子x_{k-1}^{(L)} から潜在的な x_k を生成し、ステップ 6で棄却 / 受領の判定が行われる。言い替えれば、x_k の値はそれまでに生成された x_{k-1} を用いて生成される。

その他の粒子フィルタ[編集]

関連項目[編集]

参考文献[編集]

  • モンテカルロフィルタおよび平滑化について, by 北川源四郎. 統計数理, Vol.44, No.1, pp. 31--48, 1996
  • Sequential Monte Carlo Methods in Practice, by A Doucet, N de Freitas and N Gordon. Springer, 2001. ISBN 978-0387951461
  • On Sequential Monte Carlo Sampling Methods for Bayesian Filtering, by A Doucet, C Andrieu and S. Godsill, Statistics and Computing, vol. 10, no. 3, pp. 197-208, 2000 CiteSeer link
  • Tutorial on Particle Filters for On-line Nonlinear/Non-Gaussian Bayesian Tracking (2001); S. Arulampalam, S. Maskell, N. Gordon and T. Clapp; CiteSeer link
  • Particle Methods for Change Detection, System Identification, and Control, by Andrieu C., Doucet A., Singh S.S., and Tadic V.B., Proceedings of the IEEE, Vol 92, No 3, March 2004. SMC Homepage link
  • A Tutorial on Bayesian Estimation and Tracking Techniques Applicable to Nonlinear and Non-Gaussian Processes, A.J. Haug, The MITRE Corporation; Mitre link
  • Distributed Implementations of Particle Filters, Anwer Bashi, Vesselin Jilkov, Xiao Rong Li, Huimin Chen, Available from the University of New Orleans
  • A survey of convergence results on particle filtering methods for practitioners, Crisan, D. and Doucet, A., IEEE Transactions on Signal Processing 50(2002)736-746 doi:10.1109/78.984773
  • Beyond the Kalman Filter: Particle Filters for Tracking Applications, by B Ristic, S Arulampalam, N Gordon. Artech House, 2004. ISBN 1-58053-631-x

外部リンク[編集]