粒子フィルタ

粒子フィルタ（りゅうしフィルタ、）や逐次モンテカルロ法 (ちくじモンテカルロほう、)とは、シミュレーションに基づく複雑なモデルの推定法である。1993年1月に北川源四郎がモンテカルロフィルタの名称で、1993年4月にN.J. Gordonらがブートストラップフィルタの名称で同時期に同じものを発表した。
この手法はふつうベイズモデルを推定するのに用いられ、バッチ処理であるマルコフ連鎖モンテカルロ法 (MCMC) の逐次 (オンライン) 版である。またこの手法は重点サンプリング法にも似たところがある。
うまく設計すると、粒子フィルタはMCMCよりも高速である。拡張カルマンフィルタや無香カルマンフィルタ (Unscented カルマンフィルタ) に比べて、サンプル点が十分多くなるとベイズ最適推定に近付くことからより高い精度の解が得られるので、これらの代わりに用いられることがある。また手法を組み合わせて、カルマンフィルタを粒子フィルタの提案分布として使うこともできる。
目的.
粒子フィルタの目的は、観測不可能な状態列 formula_1 formula_2 の確率分布を観測値 formula_3 formula_2 から推定することである。ベイズ推定値 formula_1 は一つ過去の確率分布 formula_6 から得られるのに対し、マルコフ連鎖法では過去の全てを含む確率分布 formula_7 より求められる.
状態空間モデル.
粒子フィルタでは観測不可能な状態 formula_1 と観測値 formula_3 は以下のように表せるとする。
観測不可能な状態列 formula_10 は以下を満たす1階マルコフ過程。つまり formula_1 は formula_12 で決まる。ただし初期値 formula_13 の確率分布 formula_14 を持つものとする。なお、これは2つ前 formula_15 の状態を使えないという意味ではなく、必要なら2つ前の状態も formula_12 に含めて使うという意味である。
観測データ列 formula_18 は formula_10 が既知であるという条件の下で独立である。換言すると formula_3 は formula_1 のみに従属する。
その上で、下記に従う状態方程式を状態空間モデルと呼ぶ。
ただし formula_25 も formula_26 も異なる formula_27 について互いに独立であり、ある確率分布に従うノイズで、formula_25 をシステムノイズ、formula_26 を観測ノイズと呼ぶ。また、formula_30 と formula_31 は既知の関数である。
カルマンフィルターの状態方程式は状態空間モデルの一種であり、formula_1 と formula_3 が実数の列ベクトル、関数 formula_30 と formula_31 が線形、formula_25 と formula_26 が多変量正規分布に従うとすると、下記形式で表現でき、
formula_1 の確率分布は多変量正規分布になり、カルマンフィルタによってベイズ推定と厳密に一致する解が得られる。線形でない場合などは、カルマンフィルタは1次近似に過ぎない。粒子フィルタもモンテカルロ法による近似には変わりないが、十分な数の粒子があれば高い精度が得られる。
モンテカルロ近似.
粒子法は他のサンプルリング法と同様に、フィルタリング確率分布formula_40を近似する点列を生成する。したがって formula_41 個のサンプル点があれば、フィルタリング確率分布による期待値は
によって近似される。そして通常のモンテカルロ法と同様に formula_43を適切に調整することで, 近似の程度に応じた次数までの モーメント を得ることができる。 
Sampling Importance Resampling (SIR).
"Sampling importance resampling (SIR)" アルゴリズムは、モンテカルロフィルタ（北川源四郎 1993）やブートストラップフィルタ（Gordon et al. 1993）による本来の粒子フィルタであるが、よく使われる粒子フィルタアルゴリズムである。ここではフィルタリング確率分布 formula_44 をformula_41個の重みつき粒子によって近似する。
"重み" formula_47 は以後の formula_48 である粒子の相対事後分布(密度)の近似となっており、
formula_48を満たす。
SIRは重点サンプリング の逐次版 (つまり、再帰的) と言える。
重点サンプリングにあるように、関数 formula_43 の推定値は重みつき平均
で近似できる。
粒子の個数が有限である場合、このアルゴリズムの性能は提案分布formula_52の選択に依存する。
最適な提案分布は目的分布
である。
しかしながら事前遷移確率
を提案分布として用いることが多い。
なぜならそこからは容易に粒子をサンプルすることができるし、その後に重みを計算することもできるからである。
提案分布として事前遷移確率を用いるSIRフィルタは一般にブートストラップフィルタ・コンデンセーションアルゴリズムとして知られている。
唯一つを除いた全ての重みがゼロに近付くこと ― アルゴリズムの縮退 ― を防ぐために再サンプルが行なわれる。
このアルゴリズムの性能は再サンプリング方式の選びかたにもかかっている。
北川源四郎 (1993) によって提案された層化抽出法は分散の意味で最適である。
SIR法の1ステップは以下の様になる。
このときformula_59
ならば更新式は
のように簡単になる。
"Sequential Importance Resampling" 等の用語は SIR フィルターの意味で用いられることがある。
逐次的 Importance Sampling (SIS).
再サンプルが無い点を除いて SIR と同様である。
直接法.
直接法は他の粒子フィルタに比べて簡単で、合成と棄却を利用する。
formula_27 番目の一つのサンプル formula_69 を formula_70 から生成するのに以下の手順を踏む。
目的はformula_27における P 個の粒子を、formula_92 だけから生成することである。 
これにはマルコフ方程式が formula_12 だけから formula_1 を生成できるように記述できなければならない。
このアルゴリズムは formula_27 における粒子を生成するためにformula_92におけるformula_41個の粒子の合成を利用しており、formula_27においてformula_41個の粒子が生成されるまでステップ 2--6 以降を繰り返す。
これは formula_69 を2次元配列とみなすと、より視覚的に理解できる。
一つの次元は formula_27 であり、もう一つの次元は粒子番号formula_73である。
例えば formula_103 は formula_27 におけるformula_73番目の粒子であり、formula_106 と書ける
(上記のアルゴリズムの記述に用いた通り)。
ステップ 3 は formula_92 においてランダムに選ばれた粒子formula_108 から"潜在的な" formula_1 を生成し、ステップ 6で棄却 / 受領の判定が行われる。言い替えれば、formula_1 の値はそれまでに生成された formula_12 を用いて生成される。