粒子フィルタ

粒子フィルタ（Particle Filter）とは

粒子フィルタ、または逐次モンテカルロ法（Sequential Monte Carlo; SMC）は、シミュレーションに基づいた複雑なモデルの推定手法です。この手法は、1993年に北川源四郎が「モンテカルロフィルタ」として、また同時期にN.J. Gordonらが「ブートストラップフィルタ」として発表しました。

概要

粒子フィルタは、主にベイズモデルを推定するために使用され、バッチ処理であるマルコフ連鎖モンテカルロ法（MCMC）の逐次（オンライン）版と見なすことができます。重点サンプリング法とも関連性があります。適切に設計された粒子フィルタは、MCMCよりも高速に処理できる場合があります。また、拡張カルマンフィルタや無香カルマンフィルタと比較して、十分なサンプル点があればベイズ最適推定に近づき、より高精度な解を得られることがあります。そのため、これらのフィルタの代替として使用されることがあります。

目的

粒子フィルタの主な目的は、観測できない状態列 \( x_k \) の確率分布を、観測値 \( y_k \) から推定することです。ベイズ推定では、過去の観測値に基づいて現在の状態 \( x_k \) の確率分布 \( p(x_k | y_0, y_1, ..., y_k) \) を求めますが、マルコフ連鎖法では過去全ての状態と観測値を含む確率分布 \( p(x_0, x_1, ..., x_k | y_0, y_1, ..., y_k) \) を計算します。

状態空間モデル

粒子フィルタでは、観測不可能な状態 \( x_k \) と観測値 \( y_k \) が以下の関係を持つと仮定します。

観測不可能な状態列 \( x_0, x_1, ... \) は、一次マルコフ過程に従います。つまり、\( x_k \) は一つ前の状態 \( x_{k-1} \) によって決定されます。初期状態 \( x_0 \) は確率分布 \( p(x_0) \) に従います。必要であれば、\( x_{k-1} \) に過去の状態を含めることも可能です。

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

観測データ列 \( y_0, y_1, ... \) は、状態 \( x_0, x_1, ... \) が既知であるという条件下で独立です。したがって、\( 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 \) は互いに独立なノイズであり、それぞれシステムノイズと観測ノイズと呼ばれます。また、\( f_k(·) \) と \( h_k(·) \) は既知の関数です。

カルマンフィルタの状態方程式は、状態空間モデルの一種であり、\( x_k \) と \( y_k \) が実数の列ベクトル、\( f_k(·) \) と \( h_k(·) \) が線形関数、\( v_k \) と \( w_k \) が多変量正規分布に従う場合、以下の形式で表現できます。

\( x_k = F_k x_{k-1} + G_k v_k \)
\( y_k = H_k x_k + w_k \)

この場合、\( x_k \) の確率分布は多変量正規分布となり、カルマンフィルタによって厳密なベイズ推定解が得られます。しかし、線形でない場合はカルマンフィルタは一次近似に過ぎず、粒子フィルタはモンテカルロ法による近似であるものの、十分な数の粒子があれば高精度な解が得られます。

モンテカルロ近似

粒子フィルタは、他のサンプリング法と同様に、フィルタリング確率分布 \( p(x_k | y_0, ..., y_k) \) を近似する点列を生成します。\( P \) 個のサンプル点がある場合、フィルタリング確率分布による期待値は以下の式で近似できます。

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

関数 \( f(·) \) を適切に調整することで、近似の程度に応じた次数までのモーメントが得られます。

Sampling Importance Resampling (SIR)

SIRアルゴリズムは、モンテカルロフィルタやブートストラップフィルタとして知られる基本的な粒子フィルタアルゴリズムです。フィルタリング確率分布 \( p(x_k | y_0, ..., y_k) \) を \( P \) 個の重み付き粒子 \( \{ (w_k^{(L)}, x_k^{(L)}) : L = 1, ..., P \} \) で近似します。重み \( w_k^{(L)} \) は粒子の相対事後分布の近似であり、\( \sum_{L=1}^{P} w_k^{(L)} = 1 \) を満たします。

SIRは重点サンプリングの逐次版と見なすことができ、関数 \( f(·) \) の推定値は以下の重み付き平均で近似できます。

\( \int f(x_k) p(x_k | y_0, ..., 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フィルタは、一般にブートストラップフィルタまたはコンデンセーションアルゴリズムと呼ばれます。アルゴリズムの縮退を防ぐために、再サンプリングが行われます。再サンプリング方式の選択もアルゴリズムの性能に影響し、層化抽出法が分散の意味で最適とされています。

SIR法のステップ

1. 提案分布から粒子をサンプリング：\( x_k^{(L)} \sim \pi(x_k | x_{0:k-1}^{(L)}, y_{0:k}) \) （\( L = 1, ..., P \)）
2. 重みを更新：\( \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. 重みを正規化：\( w_k^{(L)} = \frac{\hat{w}_k^{(L)}}{\sum_{J=1}^{P} \hat{w}_k^{(J)}} \) （\( L = 1, ..., P \)）
4. 有効粒子数の推定：\( \hat{N}_{eff} = \frac{1}{\sum_{L=1}^{P} (w_k^{(L)})^2} \)
5. 有効粒子数が閾値より小さければ再サンプリングを実行
a. 重みに比例した確率で \( P \) 個の粒子を描画
b. 重みを均一化：\( w_k^{(L)} = 1/P \) （\( L = 1, ..., P \)）

その他の粒子フィルタ

• - 逐次的 Importance Sampling (SIS)
• - 直接法
• - Auxiliary particle filter
• - Gaussian particle filter
• - Unscented particle filter
• - Monte Carlo particle filter
• - Gauss-Hermite particle filter
• - Cost Reference particle filter

参考文献

• - モンテカルロフィルタおよび平滑化について, 北川源四郎, 1996
• - Sequential Monte Carlo Methods in Practice, Doucet, de Freitas, Gordon, 2001
• - On Sequential Monte Carlo Sampling Methods for Bayesian Filtering, Doucet, Andrieu, Godsill, 2000
• - Tutorial on Particle Filters for On-line Nonlinear/Non-Gaussian Bayesian Tracking, Arulampalam, Maskell, Gordon, Clapp, 2001
• - Particle Methods for Change Detection, System Identification, and Control, Andrieu C., Doucet A., Singh S.S., and Tadic V.B., 2004
• - A Tutorial on Bayesian Estimation and Tracking Techniques Applicable to Nonlinear and Non-Gaussian Processes, A.J. Haug
• - Distributed Implementations of Particle Filters, Anwer Bashi, Vesselin Jilkov, Xiao Rong Li, Huimin Chen
• - A survey of convergence results on particle filtering methods for practitioners, Crisan, D. and Doucet, A., 2002
• - Beyond the Kalman Filter: Particle Filters for Tracking Applications, Ristic, Arulampalam, Gordon, 2004

関連項目

• - アンサンブルカルマンフィルタ
• - カルマンフィルタ
• - 反復ベイズ推定
• - 北川源四郎

外部リンク

• - 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

もう一度検索