粒子フィルタ

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

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

目的

粒子フィルタの目的は、観測不可能な状態列 ${\displaystyle x_{k}}$ ${\displaystyle (k=0,1,2,3,\ldots )}$ の確率分布を観測値 ${\displaystyle y_{k}}$ ${\displaystyle (k=0,1,2,3,\ldots )}$ から推定することである。ベイズ推定値 ${\displaystyle x_{k}}$ は一つ過去の確率分布 ${\displaystyle p(x_{k}|y_{0},y_{1},\ldots ,y_{k})}$ から得られるのに対し、マルコフ連鎖法では過去の全てを含む確率分布 ${\displaystyle p(x_{0},x_{1},\ldots ,x_{k}|y_{0},y_{1},\ldots ,y_{k})}$ より求められる.

状態空間モデル

粒子フィルタでは観測不可能な状態 ${\displaystyle x_{k}}$ と観測値 ${\displaystyle y_{k}}$ は以下のように表せるとする。

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

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

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

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

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

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

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

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

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

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

${\displaystyle x_{k}}$ の確率分布は多変量正規分布になり、カルマンフィルタによってベイズ推定と厳密に一致する解が得られる。線形でない場合などは、カルマンフィルタは1次近似に過ぎない。粒子フィルタもモンテカルロ法による近似には変わりないが、十分な数の粒子があれば高い精度が得られる。

モンテカルロ近似

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

${\displaystyle \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)})}$

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

Sampling Importance Resampling (SIR)

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

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

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

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

${\displaystyle \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)})}$

で近似できる。

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

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

である。

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

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

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

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

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

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

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

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

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

${\displaystyle {\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})}}}$

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

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

のように簡単になる。

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

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

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

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

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

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

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

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

逐次的 Importance Sampling (SIS)

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

直接法

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

1) ${\displaystyle n=1}$とする。

2) ${\displaystyle \{1,...,P\}}$上の一様分布から${\displaystyle L}$を生成する

3) テスト粒子 ${\displaystyle {\hat {x}}}$ を確率分布 ${\displaystyle p_{x_{k}|x_{k-1}}(x|x_{k-1|k-1}^{(L)})}$ から生成する

4) ${\displaystyle {\hat {y}}}$ の確率を ${\displaystyle {\hat {x}}}$ を用いて ${\displaystyle p_{y|x}(y_{k}|{\hat {x}})}$ から生成する。ただし、${\displaystyle y_{k}}$ は観測値である

5) 別の数${\displaystyle u}$を${\displaystyle [0,m_{k}]}$上の一様分布からを生成する。ただし、${\displaystyle m_{k}=\sup _{x}p_{y|x}(y_{k}|x)}$

6) ${\displaystyle u}$ と ${\displaystyle {\hat {y}}}$ を比較

6a) ${\displaystyle u}$ が大きければ step 2 以降を繰り返す

6b) ${\displaystyle u}$ が小さかったら ${\displaystyle {\hat {x}}}$ を ${\displaystyle x{k|k}^{(p)}}$ として保存して ${\displaystyle n}$ を1増やす

7) ${\displaystyle n>P}$ ならば終了

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

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