Particle Filter を学ぶ：理論編
本記事は 2 部構成となっております．
本記事は個人学習の整理として作成しています
TL;DR
- Particle Filterは，事後分布を大量の「粒子（仮説）」の集まりで表現する逐次ベイズ推定
- アルゴリズムは 3 ステップの繰り返し:予測（粒子を動かす）→ 重み付け（観測との整合性で重みを更新）→ リサンプリング（重みに比例して選び直す）
- リサンプリングを行わないと重みが少数の粒子に集中し，有効な粒子数が急速に減少する（縮退）．やり過ぎると多様性が失われる（サンプル貧化）．ESS（有効サンプルサイズ） を目安にリサンプリングを制御する
- 強みは非線形・非ガウスモデルにも柔軟に対応できること．弱みは計算コストが大きいこと，特に高次元では粒子数を増やさないと近似精度が落ちやすいこと
モチベーション
先日， Kaggle で ROGII - Wellbore Geology Prediction という穴掘りコンペが終わりました．
概要としては，目で見えない地中の掘削の状態をセンサ等から推定するというものでした．
そのコンペで Particle Filter という手法がよく使われていたので，その手法についてまとめてみました．
余談ですが，コンペの成績はかなりギリギリで銅メダル圏内に入れました．<3<3
0. イントロ
「観測にはノイズ・誤差が乗るが，状態をリアルタイムに推定したい」
こういう問題に対する定番の答えは Kalman Filter が有名です．
しかし，Kalman Filter は線形ガウス状態空間モデルという仮定の上に成り立っています．（非線形モデルに対応する EKF や UKF などもありますが，いずれも状態分布を主にガウス分布として近似します．）
Particle Filter は，非線形・非ガウス分布のモデルにも対応できる手法です．
この記事では，Particle Filterを
- 何を計算したいのか（状態空間モデルとベイズフィルタ）
- どうやって計算するのか（予測・重み付け・リサンプリング）
- 何にハマるのか（重みの縮退，サンプル貧化）
の順に解説します．
1. 問題設定：状態空間モデル
Particle Filter が扱うのは状態空間モデルです．
- 状態 $x_t$：時刻 $t$ におけるシステムの状態を表す，推定したい量（ロボットの位置，物体の速度，etc.）．多くの場合，直接観測できません
- 観測 $y_t$：手に入るがノイズを含む量（センサ値，etc.）
状態空間モデルの中心となるのは，状態遷移モデルと観測モデルの 2 つです．
状態遷移モデル
「前の状態 $x_{t-1}$ から次の状態 $x_t$ がどうなるか」「時間とともに状態がどう変化するか」を表す確率分布です．
$$
\begin{aligned}
x_t &\sim p(x_t \mid x_{t-1})
\end{aligned}
$$
この予測にはランダムなノイズが含まれます．例えばロボットの位置なら，モータの誤差や地面の摩擦などで，前の位置から正確に次の位置に行けるわけではありません．そのため，確率分布として扱います．
この分布を状態遷移モデルと呼びます．
確率モデルを表現する典型は $x_t = f(x_{t-1}) + v_t$（$v_t$ はノイズ）という形です．
観測モデル
「状態 $x_t$ から観測 $y_t$ がどうなるか」「今の状態から，センサや観測結果がどれくらい出やすいか」を表す確率分布です．
$$
\begin{aligned}
y_t &\sim p(y_t \mid x_t)
\end{aligned}
$$
例えば，センサの誤差やノイズで，状態が同じでも観測値は毎回少しずつ変わります．そのため，確率分布として扱います．
確率モデルを表現する典型は $y_t = h(x_t) + e_t$（$e_t$ はノイズ）という形です．
フィルタリング分布（事後分布）
「今の時刻 $t$ の状態 $x_t$ を，時刻 $t$ までの観測 $y_{1:t}$ を使って推定した確率分布」をフィルタリング分布と呼びます．
やりたいことは，観測列 $y_{1:t} = (y_1, \dots, y_t)$ が手に入るたびに，フィルタリング分布
$$
p(x_t \mid y_{1:t})
$$
を更新し続けることです．
「いま状態はどこにありそうか」を，点ではなく分布ごと追いかけるのがポイントです．
Kalman Filter では，状態遷移・観測モデルを線形とし，それぞれのノイズをガウス分布と仮定します．その結果，フィルタリング分布もガウス分布となるため，平均と共分散だけで表現できます．一方で，Particle Filter は「粒子の集まり」で表現します．
2. ベイズフィルタ：予測と更新の繰り返し
ベイズフィルタリングは，次の「予測」と「更新」の2ステップを繰り返すことで行われます．
予測：「状態遷移モデルで 1 時刻進める」操作
$$
p(x_t \mid y_{1:t-1}) = \int p(x_t \mid x_{t-1}) p(x_{t-1} \mid y_{1:t-1}) dx_{t-1}
$$
予測では，前時刻の状態を状態遷移モデルによって 1 時刻進めます．前時刻の状態 $x_{t-1}$ は 1 つに確定しているわけではなく，$p(x_{t-1} \mid y_{1:t-1})$ という分布で表されるので，積分で「すべての可能な前時刻の状態」を考慮します．
つまり，「前の状態がどれくらいありそうか」×「そこから現在の状態へ遷移する確率」を，前の状態のすべての可能性について合計しています．（周辺化）
更新：「観測でベイズ更新する」操作
$$
p(x_t \mid y_{1:t}) \propto p(y_t \mid x_t) p(x_t \mid y_{1:t-1})
$$
予測ステップでは，時刻 $t$ の状態がどこにありそうかという予測分布 $p(x_t \mid y_{1:t-1})$ を求めました．ここで新しい観測 $y_t$ が得られた場合に，「もし状態が $x_t$ なら，実際に得られた観測 $y_t$ がどれくらい起こりやすいか」を尤度 $p(y_t \mid x_t)$ で評価します．
そして，予測時点での $x_t$ のありそう度と，観測の尤度を掛け合わせることで，観測を考慮したフィルタリング分布 $p(x_t \mid y_{1:t})$ を求めます．
これはベイズの定理そのものです．
図にするとこうなります．「予測で分布を 1 時刻先へ進め，観測との整合性によって分布を更新する」を繰り返すのがベイズフィルタです．
これで話が終われば楽なのですが，一般の状態空間モデルでは，この予測・更新を厳密に計算することは難しい場合があります．
厳密に逐次計算できる代表的なケースとして，次の2つがあります．
| 条件 | 解ける理由 | アルゴリズム |
|---|---|---|
| 状態遷移・観測モデルが線形 & ノイズがガウス | 予測分布・フィルタリング分布が常にガウスのまま | Kalman Filter |
| 状態が離散で有限個 | 積分が有限和になる | HMM の forward algorithm |
一般の非線形・非ガウスな状態空間モデルでは，これらの積分やベイズ更新をそのまま計算できないことが多く，近似が必要になります．Particle Filter はその代表的な方法の1つです．
3. Particle Filter：分布を「粒子」で持つ
Particle Filter のアイデアはシンプルです．
複雑な確率分布そのものを数式で表現する代わりに，「この状態かもしれない」という多数の候補（粒子）と，それぞれの確からしさ（重み）で分布を近似します．
例えば，ある状態について
x=10 あたりに粒子がたくさんある
x=20 あたりには少ししかない
さらに x=10 付近の粒子の重みが大きい
のであれば，「状態は 10 付近にいる可能性が高い」という分布を表現できます．
この粒子による近似は，次のように書けます．
$$
p(x_t \mid y_{1:t}) \approx \sum_{i=1}^{N} w_t^{(i)} \delta\left(x_t - x_t^{(i)}\right),
\qquad \sum_i w_t^{(i)} = 1
$$
$\delta\left(x_t - x_t^{(i)}\right)$ は直感的には「$x_t^{(i)}$ の位置に粒子が 1 つある」ぐらいのイメージで，その位置に確率質量を置くことを表す記号です．式の細部よりも，分布を重み付きの粒子の集合として表していると理解すれば十分です．
$x_t^{(i)}$ が $i$ 番目の粒子（状態の仮説），$w_t^{(i)}$ がその重みです．粒子の配置と重みを合わせて，「どの状態がどれくらいありそうか」を表します．
こうすると，先ほどの複雑な予測・更新の計算を，粒子ごとの単純な操作で近似できるようになります．最初に粒子を用意した後，1 サイクルは次の 3 ステップです．
以下では，粒子を状態遷移モデルそのものから生成する，最も基本的な Particle Filter（bootstrap Particle Filter）を考えます．
⓪ 初期化 — 最初の粒子を用意する（初回のみ）
時刻 $t=0$ では，初期状態の分布 $p(x_0)$ から $N$ 個の粒子をサンプリングします．
$$
x_0^{(i)} \sim p(x_0), \quad i=1,\dots,N
$$
最初は全ての粒子を同じ重みとして，
$$
w_0^{(i)} = \frac{1}{N}
$$
とします．
① 予測 — 各粒子を状態遷移モデルで動かす
各粒子から，状態遷移モデルに従って次時刻の状態をサンプリングします．
$$
x_t^{(i)} \sim p(x_t \mid x_{t-1}^{(i)})
$$
実装上は「$f$ で動かしてノイズを足す」だけです．
解析的な積分の代わりに，粒子を実際に動かすことで予測分布を近似します．
② 重み付け — 観測との整合性を重みにする
観測 $y_t$ が来たら，各粒子について観測との整合性（尤度）を評価し，重みを更新します．
$$
w_t^{(i)} \propto w_{t-1}^{(i)} p(y_t \mid x_t^{(i)})
$$
観測とよく合う粒子は重く，合わない粒子は軽くなります．最後に総和が 1 になるよう正規化します．
直前にリサンプリングしている場合は，全粒子の $w_{t-1}^{(i)} = 1/N$ なので，実質的には観測尤度 $p(y_t \mid x_t^{(i)})$ だけで重みが決まります．
点推定が欲しければ，重み付き平均を使うことが多いです．ただし，事後分布が多峰な場合には，重み付き平均が山と山の間の確率の低い位置を指すこともあるため，粒子分布そのものを見ることも重要です．
③ リサンプリング — 重みに比例して選び直す
このまま①②だけを繰り返すと，やがて一部の粒子だけに重みが集中する「縮退」が起こります．そこで，必要に応じて重みに比例して粒子を選び直し，重みを $1/N$ にリセットします．重い粒子は複数回コピーされ，軽い粒子は消えます．
4. 重みの縮退とリサンプリング
縮退（degeneracy）：少数の粒子に重みが集中する
リサンプリングを行わずに重み更新を繰り返すと，各時刻の観測尤度が重みに掛け合わされていくため，粒子間の小さな尤度差が時間とともに蓄積します．その結果，少数の粒子に重みが集中し，多くの粒子が推定にほとんど寄与しなくなります．これが重みの縮退です．Particle Filter では，この問題を抑えるためにリサンプリングを行います．
縮退の度合いを測る定番の指標が 有効サンプルサイズ（ESS） です．重みを総和 1 に正規化したとき，ESS は次の式で定義されます．
$$
\mathrm{ESS} = \frac{1}{\sum_{i=1}^{N} \left(w_t^{(i)}\right)^2}
$$
重みが均等なら $\mathrm{ESS}=N$，1 粒子に集中すると $\mathrm{ESS}=1$ になります．「実質何粒子で推定しているか」という量です．
実際に非線形モデルで Particle Filter を回して，リサンプリングの有無で ESS と最大重みを比べたのが次の図です．リサンプリングなしでは数ステップで ESS がほぼ 1（最大重みがほぼ 1）まで低下するのが分かります．
そこで，ESS が閾値（よく使うのは $N/2$）を下回ったらリサンプリングする，が定石になります．毎ステップ無条件にリサンプリングする実装も多いですが，リサンプリング自体もモンテカルロ誤差を足す操作なので，必要なときだけやる方が筋が良いです．
系統リサンプリング（systematic resampling）
リサンプリングの実装にはいくつか流儀があります．素朴なのは重みを確率とする選び直し（多項リサンプリング）ですが，分散が大きめです．実務の定番は系統リサンプリングで，くじを $N$ 本独立に引く代わりに等間隔の $N$ 本を 1 回だけずらして引きます．
- 区間 $[0, 1/N)$ から乱数 $u_1$ を 1 つだけ引く
- $u_j = u_1 + (j-1)/N$ の $N$ 点で累積重みを引き当てる
$O(N)$ で実装でき，実用上は多項リサンプリングよりもばらつきが小さくなることが多いため，よく使われます．
サンプル貧化（sample impoverishment）：やり過ぎの副作用
リサンプリングは万能ではありません．重い粒子をコピーするだけなので，繰り返すと粒子の多様性が失われ，同じ値のコピーだらけになります（サンプル貧化）．特に状態遷移ノイズが小さいモデルで顕著です．
対策としてよく使われるのは:
- ESS 閾値でリサンプリング頻度を抑える
- リサンプリング後の粒子に小さなノイズ（ジッタ）を足して散らす
- 必要に応じて，粒子の多様性を保つため状態遷移ノイズを調整する
いずれも「現在の事後分布を正確に表すこと」と「粒子の多様性を維持すること」のトレードオフです．
5. まとめ
今回は，Particle Filter の基本的な考え方とアルゴリズムを解説しました．
- 事後分布を大量の「粒子（仮説）」の集まりで表現する逐次ベイズ推定
- アルゴリズムは 3 ステップの繰り返し:予測（粒子を動かす）→ 重み付け（観測との整合性で重みを更新）→ リサンプリング（重みに比例して選び直す）
- リサンプリングを行わないと重みが少数の粒子に集中し，有効な粒子数が急速に減少する（縮退）．やり過ぎると多様性が失われる（サンプル貧化）．ESS（有効サンプルサイズ） を目安にリサンプリングを制御する
- 強みは非線形・非ガウスモデルにも柔軟に対応できること．弱みは計算コストが大きいこと，特に高次元では粒子数を増やさないと近似精度が落ちやすいこと
また Python での実装例や，穴掘りコンペのまとめ記事などは別途 Qiita に投稿しようと思っています．
お楽しみに〜．