粒子フィルタ

粒子フィルタ（りゅうしフィルタ、英: particle filter）や逐次モンテカルロ法 (ちくじモンテカルロほう、英: sequential Monte Carlo; SMC)とは、シミュレーションに基づく複雑なモデルの推定法である。1993年1月に北川源四郎がモンテカルロフィルタの名称で、1993年4月にN.J. Gordonらがブートストラップフィルタの名称で同時期に同じものを発表した。

この手法はふつうベイズモデルを推定するのに用いられ、バッチ処理であるマルコフ連鎖モンテカルロ法 (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}$ は以下のように表せるとする。

観測不可能な状態列 $x_{0},x_{1},\ldots$ は以下を満たす1階マルコフ過程。つまり $x_{k}$ は $x_{k-1}$ で決まる。ただし初期値 $x_{0}$ の確率分布 $p(x_{0})$ を持つものとする。なお、これは2つ前 $x_{k-2}$ の状態を使えないという意味ではなく、必要なら2つ前の状態も $x_{k-1}$ に含めて使うという意味である。

$x_{k}|x_{k-1}\sim p_{x_{k}|x_{k-1}}(x|x_{k-1})$

観測データ列 $y_{0},y_{1},\ldots$ は $x_{0},x_{1},\ldots$ が既知であるという条件の下で独立である。換言すると $y_{k}$ は $x_{k}$ のみに従属する。

$y_{k}|x_{k}\sim p_{y|x}(y|x_{k})$

その上で、下記に従う状態方程式を状態空間モデルと呼ぶ。

$x_{k}=f_{k}(x_{k-1},v_{k})\,$

$y_{k}=h_{k}(x_{k},w_{k})\,$

ただし $v_{k}$ も $w_{k}$ も異なる $k$ について互いに独立であり、ある確率分布に従うノイズで、$v_{k}$ をシステムノイズ、$w_{k}$ を観測ノイズと呼ぶ。また、$f_{k}(\cdot )$ と $h_{k}(\cdot )$ は既知の関数である。

カルマンフィルターの状態方程式は状態空間モデルの一種であり、$x_{k}$ と $y_{k}$ が実数の列ベクトル、関数 $f_{k}(\cdot )$ と $h_{k}(\cdot )$ が線形、$v_{k}$ と $w_{k}$ が多変量正規分布に従うとすると、下記形式で表現でき、

${\begin{aligned}x_{k}&=F_{k}x_{k-1}+G_{k}v_{k}\\y_{k}&=H_{k}x_{k}+w_{k}\end{aligned}}$

$x_{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 {\frac {1}{P}}\sum _{L=1}^{P}f(x_{k}^{(L)})$

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

Sampling Importance Resampling (SIR)

Sampling importance resampling (SIR) アルゴリズムは、モンテカルロフィルタ（北川源四郎 1993）やブートストラップフィルタ（Gordon et al. 1993）による本来の粒子フィルタであるが、よく使われる粒子フィルタアルゴリズムである。ここではフィルタリング確率分布 $p(x_{k}|y_{0},\ldots ,y_{k})$ を$P$個の重みつき粒子によって近似する。

$\{(w_{k}^{(L)},x_{k}^{(L)})~:~L=1,\ldots ,P\}$.

重み $w_{k}^{(L)}$ は以後の $\sum _{L=1}^{P}w_{k}^{(L)}=1$ である粒子の相対事後分布(密度)の近似となっており、 $\sum _{L=1}^{P}w_{k}^{(L)}=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フィルタは一般にブートストラップフィルタ・コンデンセーションアルゴリズムとして知られている。

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

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

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

$x_{k}^{(L)}\sim \pi (x_{k}|x_{0:k-1}^{(L)},y_{0:k})$

2) $L=1,\ldots ,P$ について、重みを更新し、正規化定数を得る。

${\hat {w}}_{k}^{(L)}=w_{k-1}^{(L)}{\frac {p(y_{k}|x_{k}^{(L)})p(x_{k}^{(L)}|x_{k-1}^{(L)})}{\pi (x_{k}^{(L)}|x_{0:k-1}^{(L)},y_{0:k})}}$

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

${\hat {w}}_{k}^{(L)}=w_{k-1}^{(L)}p(y_{k}|x_{k}^{(L)}),$

のように簡単になる。

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

$w_{k}^{(L)}={\frac {{\hat {w}}_{k}^{(L)}}{\sum _{J=1}^{P}{\hat {w}}_{k}^{(J)}}}$

4) 有効粒子数の推定量を計算する。

${\hat {N}}_{\mathit {eff}}={\frac {1}{\sum _{L=1}^{P}\left(w_{k}^{(L)}\right)^{2}}}$

5) もし有効粒子数が与えられた閾値より少なかったら ${\hat {N}}_{\mathit {eff}}<N_{thr}$ 再サンプルを実行する。

a) 現在の粒子の集合から、重みに比例した確率で $P$ 個の粒子を描く。現在の粒子の集合を新しい粒子の集合で置き換える。

b) $L=1,\ldots ,P$ について、 $w_{k}^{(L)}=1/P$ とする。

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

逐次的 Importance Sampling (SIS)

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

直接法

直接法は他の粒子フィルタに比べて簡単で、合成と棄却を利用する。 $k$ 番目の一つのサンプル $x$ を $p_{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}$ を用いて生成される。

その他の粒子フィルタ

• Auxiliary particle filter
• Gaussian particle filter
• Unscented particle filter
• Monte Carlo particle filter
• Gauss-Hermite particle filter
• Cost Reference particle filter

関連項目

• アンサンブルカルマンフィルタ
• カルマンフィルタ、ガウス分布の解析的な推定器
• 反復ベイズ推定
• 北川源四郎

参考文献

• モンテカルロフィルタおよび平滑化について, 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.

外部リンク

• Sequential Monte Carlo Methods (Particle Filtering) homepage on University of Cambridge
• Dieter Fox's MCL Animations
• Nando de Freitas' free software
• Rob Hess' free software