粒子フィルタ

粒子フィルタ（りゅうしフィルタ、英: particle filter）や逐次モンテカルロ法 (ちくじモンテカルロほう、英: sequential Monte Carlo; SMC)とは、シミュレーションに基づく複雑なモデルの推定法である。1993年1月に北川源四郎がモンテカルロフィルタの名称で、1993年4月にN.J. Gordonらがブートストラップフィルタの名称で同時期に同じものを発表した。 この手法はふつうベイズモデルを推定するのに用いられ、バッチ処理であるマルコフ連鎖モンテカルロ法 (MCMC) の逐次 (オンライン) 版である。またこの手法は重点サンプリング法にも似たところがある。 うまく設計すると、粒子フィルタはMCMCよりも高速である。拡張カルマンフィルタや無香カルマンフィルタ (Unscented カルマンフィルタ) に比べて、サンプル点が十分多くなるとベイズ最適推定に近付くことからより高い精度の解が得られるので、これらの代わりに用いられることがある。また手法を組み合わせて、カルマンフィルタを粒子フィルタの提案分布として使うこともできる。 粒子フィルタの目的は、観測不可能な状態列 x k {\displaystyle x_{k}} ( k = 0 , 1 , 2 , 3 , … ) {\displaystyle (k=0,1,2,3,\ldots )} の確率分布を観測値 y k {\displaystyle y_{k}} ( k = 0 , 1 , 2 , 3 , … ) {\displaystyle (k=0,1,2,3,\ldots )} から推定することである。ベイズ推定値 x k {\displaystyle x_{k}} は一つ過去の確率分布 p ( x k | y 0 , y 1 , … , y k ) {\displaystyle p(x_{k}|y_{0},y_{1},\ldots ,y_{k})} から得られるのに対し、マルコフ連鎖法では過去の全てを含む確率分布 p ( x 0 , x 1 , … , x k | y 0 , y 1 , … , y k ) {\displaystyle p(x_{0},x_{1},\ldots ,x_{k}|y_{0},y_{1},\ldots ,y_{k})} より求められる. 粒子フィルタでは観測不可能な状態 x k {\displaystyle x_{k}} と観測値 y k {\displaystyle y_{k}} は以下のように表せるとする。 観測不可能な状態列 x 0 , x 1 , … {\displaystyle x_{0},x_{1},\ldots } は以下を満たす1階マルコフ過程。つまり x k {\displaystyle x_{k}} は x k − 1 {\displaystyle x_{k-1}} で決まる。ただし初期値 x 0 {\displaystyle x_{0}} の確率分布 p ( x 0 ) {\displaystyle p(x_{0})} を持つものとする。なお、これは2つ前 x k − 2 {\displaystyle x_{k-2}} の状態を使えないという意味ではなく、必要なら2つ前の状態も x k − 1 {\displaystyle x_{k-1}} に含めて使うという意味である。 x k | x k − 1 ∼ p x k | x k − 1 ( x | x k − 1 ) {\displaystyle x_{k}|x_{k-1}\sim p_{x_{k}|x_{k-1}}(x|x_{k-1})} 観測データ列 y 0 , y 1 , … {\displaystyle y_{0},y_{1},\ldots } は x 0 , x 1 , … {\displaystyle x_{0},x_{1},\ldots } が既知であるという条件の下で独立である。換言すると y k {\displaystyle y_{k}} は x k {\displaystyle x_{k}} のみに従属する。 y k | x k ∼ p y | x ( y | x k ) {\displaystyle y_{k}|x_{k}\sim p_{y|x}(y|x_{k})} その上で、下記に従う状態方程式を状態空間モデルと呼ぶ。 x k = f k ( x k − 1 , v k ) {\displaystyle x_{k}=f_{k}(x_{k-1},v_{k})\,} y k = h k ( x k , w k ) {\displaystyle y_{k}=h_{k}(x_{k},w_{k})\,} ただし v k {\displaystyle v_{k}} も w k {\displaystyle w_{k}} も異なる k {\displaystyle k} について互いに独立であり、ある確率分布に従うノイズで、 v k {\displaystyle v_{k}} をシステムノイズ、 w k {\displaystyle w_{k}} を観測ノイズと呼ぶ。また、 f k ( ⋅ ) {\displaystyle f_{k}(\cdot )} と h k ( ⋅ ) {\displaystyle h_{k}(\cdot )} は既知の関数である。 カルマンフィルターの状態方程式は状態空間モデルの一種であり、 x k {\displaystyle x_{k}} と y k {\displaystyle y_{k}} が実数の列ベクトル、関数 f k ( ⋅ ) {\displaystyle f_{k}(\cdot )} と h k ( ⋅ ) {\displaystyle h_{k}(\cdot )} が線形、 v k {\displaystyle v_{k}} と w k {\displaystyle w_{k}} が多変量正規分布に従うとすると、下記形式で表現でき、 x k = F k x k − 1 + G k v k y k = H k x k + 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}}} x k {\displaystyle x_{k}} の確率分布は多変量正規分布になり、カルマンフィルタによってベイズ推定と厳密に一致する解が得られる。線形でない場合などは、カルマンフィルタは1次近似に過ぎない。粒子フィルタもモンテカルロ法による近似には変わりないが、十分な数の粒子があれば高い精度が得られる。 粒子法は他のサンプルリング法と同様に、フィルタリング確率分布 p ( x k | y 0 , … , y k ) {\displaystyle p(x_{k}|y_{0},\dots ,y_{k})} を近似する点列を生成する。したがって P {\displaystyle P} 個のサンプル点があれば、フィルタリング確率分布による期待値は ∫ f ( x k ) p ( x k | y 0 , … , y k ) d x k ≈ 1 P ∑ L = 1 P f ( x k ( L ) ) {\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)})} によって近似される。そして通常のモンテカルロ法と同様に f ( ⋅ ) {\displaystyle f(\cdot )} を適切に調整することで, 近似の程度に応じた次数までの モーメント を得ることができる。 Sampling importance resampling (SIR) アルゴリズムは、モンテカルロフィルタ（北川源四郎 1993）やブートストラップフィルタ（Gordon et al. 1993）による本来の粒子フィルタであるが、よく使われる粒子フィルタアルゴリズムである。ここではフィルタリング確率分布 p ( x k | y 0 , … , y k ) {\displaystyle p(x_{k}|y_{0},\ldots ,y_{k})} を P {\displaystyle P} 個の重みつき粒子によって近似する。 { ( w k ( L ) , x k ( L ) )   :   L = 1 , … , P } {\displaystyle \{(w_{k}^{(L)},x_{k}^{(L)})~:~L=1,\ldots ,P\}} . 重み w k ( L ) {\displaystyle w_{k}^{(L)}} は以後の ∑ L = 1 P w k ( L ) = 1 {\displaystyle \sum _{L=1}^{P}w_{k}^{(L)}=1} である粒子の相対事後分布(密度)の近似となっており、 ∑ L = 1 P w k ( L ) = 1 {\displaystyle \sum _{L=1}^{P}w_{k}^{(L)}=1} を満たす。 SIRは重点サンプリング の逐次版 (つまり、再帰的) と言える。 重点サンプリングにあるように、関数 f ( ⋅ ) {\displaystyle f(\cdot )} の推定値は重みつき平均 ∫ f ( x k ) p ( x k | y 0 , … , y k ) d x k ≈ ∑ L = 1 P w ( L ) f ( x k ( L ) ) {\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)})} で近似できる。 粒子の個数が有限である場合、このアルゴリズムの性能は提案分布 π ( x k | x 0 : k − 1 , y 0 : k ) {\displaystyle \pi (x_{k}|x_{0:k-1},y_{0:k})} の選択に依存する。 最適な提案分布は目的分布 π ( 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},y_{k})} である。 しかしながら事前遷移確率 π ( x k | x 0 : k − 1 , y 0 : k ) = p ( x k | x k − 1 ) {\displaystyle \pi (x_{k}|x_{0:k-1},y_{0:k})=p(x_{k}|x_{k-1})} を提案分布として用いることが多い。 なぜならそこからは容易に粒子をサンプルすることができるし、その後に重みを計算することもできるからである。 提案分布として事前遷移確率を用いるSIRフィルタは一般にブートストラップフィルタ・コンデンセーションアルゴリズムとして知られている。 唯一つを除いた全ての重みがゼロに近付くこと ― アルゴリズムの縮退 ― を防ぐために再サンプルが行なわれる。 このアルゴリズムの性能は再サンプリング方式の選びかたにもかかっている。 北川源四郎 (1993) によって提案された層化抽出法は分散の意味で最適である。 SIR法の1ステップは以下の様になる。 1) L = 1 , … , P {\displaystyle L=1,\ldots ,P} について, 提案分布から粒子をサンプルする x k ( L ) ∼ π ( x k | x 0 : k − 1 ( L ) , y 0 : k ) {\displaystyle x_{k}^{(L)}\sim \pi (x_{k}|x_{0:k-1}^{(L)},y_{0:k})} 2) L = 1 , … , P {\displaystyle L=1,\ldots ,P} について、重みを更新し、正規化定数を得る。 w