粒子フィルタ
粒子フィルタ(りゅうしフィルタ、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{align}
x_k &= F_k x_{k-1} + G_k v_k \\
y_k &= H_k x_k + w_k
\end{align}
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\frac1P\sum_{L=1}^Pf(x_k^{(L)})
によって近似される。そして通常のモンテカルロ法と同様に f(\cdot)を適切に調整することで, 近似の程度に応じた次数までの モーメント を得ることができる。
Sampling importance resampling (SIR) アルゴリズムは、モンテカルロフィルタ(北川源四郎 1993)やブートストラップフィルタ(Gordon et al. 1993)による本来の粒子フィルタであるが、よく使われる粒子フィルタアルゴリズムである。ここではフィルタリング確率分布 p(x_k|y_0,\ldots,y_k) をP個の重みつき粒子によって近似する。
\{(w^{(L)}_k,x^{(L)}_k)~:~L=1,\ldots,P\}.
重み w^{(L)}_k は以後の \sum_{L=1}^P w^{(L)}_k = 1 である粒子の相対事後分布(密度)の近似となっており、
\sum_{L=1}^P w^{(L)}_k = 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^{(L)}_k \sim \pi(x_k|x^{(L)}_{0:k-1},y_{0:k})
2) L=1,\ldots,P について、重みを更新し、正規化定数を得る。
\hat{w}^{(L)}_k = w^{(L)}_{k-1} \frac{p(y_k|x^{(L)}_k) p(x^{(L)}_k|x^{(L)}_{k-1})} {\pi(x_k^{(L)}|x^{(L)}_{0:k-1},y_{0:k})}
このとき\pi(x_k^{(L)}|x^{(L)}_{0:k-1},y_{0:k}) = p(x^{(L)}_k|x^{(L)}_{k-1})
\hat{w}^{(L)}_k = w^{(L)}_{k-1} p(y_k|x^{(L)}_k),
のように簡単になる。
3) L=1,\ldots,Pについて、正規化された重みを計算する。
w^{(L)}_k = \frac{\hat{w}^{(L)}_k}{\sum_{J=1}^P \hat{w}^{(J)}_k}
4) 有効粒子数の推定量を計算する。
\hat{N}_\mathit{eff} = \frac{1}{\sum_{L=1}^P\left(w^{(L)}_k\right)^2}
5) もし有効粒子数が与えられた閾値より少なかったら \hat{N}_\mathit{eff} < N_{thr} 再サンプルを実行する。
a) 現在の粒子の集合から、重みに比例した確率で P 個の粒子を描く。現在の粒子の集合を新しい粒子の集合で置き換える。
b) L=1,\ldots,P について、 w^{(L)}_k = 1/P とする。
Sequential Importance Resampling 等の用語は SIR フィルターの意味で用いられることがある。
再サンプルが無い点を除いて 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} を用いて生成される。