• 粒子フィルタは，以下に述べる非線形・非ガウス状態空間モデルを推定対象として想定している．
• まず，時点「$t$」 の状態「$x_t$」が状態「$x_{t-1}$」に依存して確率的に決まる場合を考える．
• さらにその状態「$x_t$」から確率的に観測値「$y_t$」が生成されているとする．
• なお，状態は未知(直接観測することが不可能) であり，観測値は既知 (直接観測が可能) であると仮定する．
• この状況を表すと図 1 となる．

図1 状態空間モデル

• 状態「$x_t$」の時間変化を関数「$f$」，状態「$x_t$」から観測値「$y_t$」への変換を関数「$h$」で表すと，ここまでの設定は以下の 2 式に集約される．

$\Large \displaystyle x_t=f(x_{t-1},\xi_s,\upsilon_t)\ \ \ \ \ \ \ \ \ $(2.1)「システム方程式」

$\Large \displaystyle y_t=h(x_t,\xi_m,\epsilon_t)\ \ \ \ \ \ \ \ \ \ \ \ $(2,2)「観測方程式」

• 式 (2.1) をシステム方程式，式 (2.2) を観測方程式と呼び，この 2 式をまとめて「非線形・非ガウス状態空間モデル」 もしくは「一般化状態空間モデル」という (以下，状態空間モデル)．
• なお，「$\upsilon_t$」はシステムノイズ，「$\epsilon_t$」は観測ノイズであり，「$\xi_s$」と「$\xi_m$」は、 それぞれシステム方程式と観測方程式のパラメータである．
• また，「$\theta'=[\xi_{s}', \xi_{m}']$」と表す．
• なお式 (2.1)，(2.2)の詳細は付録 A1 に記述する．
• 式 (2.1) は，状態 xt が条件付き確率分布で

$\Large \displaystyle x_t \sim \mathbb{P}(x_t|x_{t-1})\ \ \ \ \ \ \ \ \ \ \ \ \ $(2-3)

• と表され，式 (2.2) は，観測値「$y_t$」が条件付き確率分布で

$\Large \displaystyle y_t \sim \mathbb{P}(y_t|x_{t})\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $(2-4)

• と表されることを示している．
• なお，ここで ∼ は「条件付き確率分布 P(xt|·) から確率変数 xt がサンプリングされること」を示す．
• なお，本稿では確率分布を P(·)，確率密度関数を p(·) と表す．
• ただし，以下では誤解が生じない限り確率分布と確率密度関数を単に分布と表現する．
• 本稿で我々が考える問題は，時点 1 から時点 t までの観測値「$y_t$」を所与とした時，式 (2.1)，(2.2) の状態空間モデルを用いて， 「未知の状態「$x_t$」を推定する」ことである．
• このような問題を状態推定という．
• 以下，時点「$1$」から時点「$t$」までの観測値「$y_t$」の組を「$y_{1:t} = \{y_1, y_2, . . . , y_t\}$」と表すことにする．
• 状態推定では，状態「$x_t$」と観測値「$y_{1:t} = \{y_1, y_2, . . . , y_t\}$」を用いて， 予測分布 (predictive distribution)，フィルタ分布 (filtering distribution)，平滑化分布 (smoothing distribution)を それぞれ以下のように定義する．

　　　予測分布：$\Large \displaystyle x_{t|t-1} = \mathbb{P}(x_t|y_{1:(t-1)})\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $(2-5-1)

　フィルタ分布：$\Large \displaystyle x_{t|t} = \mathbb{P}(x_t|y_{1:t})\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $(2-5-2)

　　平滑化分布：$\Large \displaystyle x_{t|T} = \mathbb{P}(x_t|y_{1:T})\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $(2-5-3)

• ただし，1 ≤ t ≤ T である．
• 予測分布は「$y_{1:(t-1)}$」(時点$1$から時点$t-1$までの観測値) を所与として状態「$x_t$」の分布を求めること，
• フィルタ分布は「$y_{1:t}$」(時点$1$から時点$t$までの観測値) を所与として状態「$x_t$」の分布を求めること，
• 平滑化分布は「$y_{1:T}$」(時点$1$から時点$T$までの観測値) を所与として状態「$x_t$」の分布を求めることである．

• 粒子フィルタで特徴的な発想は，推定すべき分布 (予測分布・フィルタ分布・平滑化分布) を「状態ベクトル「$x_i^t$」と その点に対応する『尤度』「$w_i^t$」 の集合で近似する」ところにある．
• 本来であれば，推定すべき分布の解析解が得られれば最も良いはずであるが，粒子フィルタが想定する非線形・非ガウス状態空間モデルでは、 そのような解析解が得られることはほとんどない．
• そのため，モンテカルロ法で状態ベクトル「$x_t^i$」を発生させ，さらにその尤度「$w_t^i$」を適切に評価して経験分布を構成する．

• たとえば，フィルタ分布「$p(x_t|y_{1:t})$」を経験分布で近似することを考えよう．
• 粒子フィルタでは以下のような記法が採用されることが多い．

$\Large \displaystyle p(x_t|y_{1:t})\simeq \sum_{i=1}^{M}\frac{w_t^i}{\sum_{i=1}^{M}w_t^i}\delta(x_t-x_{t|t-1}^i)$

• 今、時刻「$t$」のm次元ベクトルの観測時系列を「$y_t(t=1,...,T)$」とする
• この観測値に関連した観測値を表現するのに必要な変量を一つのベクトルに納めたものを「$x_t$」と書く事にする。
• この「$x_t$」は、状態ベクトルと呼ばれ、観測する事はできない量（統計では潜在変数という）である。
• 状態空間モデルと通常呼ばれる「線形・ガウス型の状態空間モデル」は、これらの確率変数間の関係を以下のような線形モデルと正規ノイズの二つの式で表したものである。

$\Large \displaystyle x_t=F_t x_{t-1}+G_t v_t\ \ \ \ \ $(1)「システムモデル」

$\Large \displaystyle y_t=H_t x_{t}+w_t\ \ \ \ \ \ \ \ \ \ \ \ $(2)「観測モデル」

• ただし、「$x_i$」は、k次元の状態ベクトルである。

• 大量の粒子をランダムに生成する。
• 粒子には位置や向きといった推定したい状態変数を好きなだけ持たせて構わない。
• それぞれの粒子は重みを持ち、系の実際の状態と粒子の状態が合致する可能性が重みによって (確率的に) 表される。
• 最初の重みは全て同じとする。

• 粒子の次の状態を予測する。
• あなたが予測した実際の系の振る舞いに基づいて、全ての粒子を移動させる。

• 観測値に基づいて粒子の重みを更新する。
• 観測値に近い粒子には観測値から遠い粒子よりも大きい重みが割り当てられる。

• 確率が非常に低い粒子を破棄して、確率の高い粒子のコピーで置き換える。

• もし必要なら、粒子の重み付き平均と重み付き共分散行列を計算して状態の推定値を得る。

• 以下の、初期値時刻「$t=t_0$」と「$N$」個の潜在変数を定めます。

$\Large \displaystyle \{ x_{t_0}^{(n)}\}_{n\in [1,...,N]}$

• 粒子フィルタではこの潜在変数の集合の元を粒子と呼び、粒子を更新していくことで潜在変数を推定していきます。
• 粒子の初期値は適当に1つずつ決めてもいいし，乱数で決めていってもいいし，通常は全て同じ値で初期化することさえ問題ありません．
• 今回は全て同じ値で初期化することにします
• (今回はハイパーパラメータ「$\theta_{sys}$」も「$\theta_{obs}$」もこの段階で何らかの値に定めます．色んな値に変えて遊びましょう．)

• 各粒子の1単位時間先の潜在変数をサンプリングします．

$\Large \displaystyle \forall n\ x_i^{(n)}\sim p(x_t^{(n)}|x_{t-1}^{(n)})$

• このステップは各粒子ごとに独立なため，並列に実行することが可能です．
• 各粒子で同じ値のものがあったとしても，サンプリングが確率的に揺らぐようなモデルである限りは1期先の予測はバラバラな値に予測されます．
• こういう事情から粒子が同じ値をもつことを許しても何の問題もないわけです．

• 所与の「$y_t$」と先に得られた潜在変数の予測値「$x_t^{(n)}$」を観測モデルに代入することで予測確率を計算します．

$\Large \displaystyle \mu_t^{(n)}=p(y_t|x_t^{(n)})$

• 各粒子の予測確率を計算したらこれらを正規化します．

$\Large \displaystyle \tilde{\mu}_t^{(n)}=\frac{\mu_t^{(n)}}{\sum_{n=1}^N\mu_t^{(n)}}$

• 所与の「$y_t$」と先に得られた潜在変数の予測値「$x_t^{(n)}$」を観測モデルに代入することで予測確率を計算します．

$\Large \displaystyle \forall n\ c_t^{(n)}\sim Categorical(c_t^{(n)}|)$

$\Large \displaystyle x_t^{(n)}=c_t^{(n)^T} \begin{pmatrix} x_t^{(1)} \\ \vdots \\ x_t^{(N)} \end{pmatrix} $

• 「$c_t^{(n)} \in \{ 0,1\}^N$」は、ある成分の値が1でそれ以外が0の値を持つ「$N$」次元のベクトルです．
• 「$c_t^{(n)} $」は、ステップ2で計算した重み「$\tilde{\mu}_t$」をパラメータとするカテゴリカル分布から生成します．
• このステップによって実現確率の低い粒子は実現確率のより高い粒子に取って代わられる可能性が出てきます．
• 極限的には最適な粒子が大多数を占めるようになっていく感じです．
• 「$t\rightarrow t+1$」として、ステップ1に戻ります．

$\large \displaystyle \large 1:Algorithm Particle_filter(X_{t-1},u_t, z_t):\\ \large 2:\ \ \ \bar{X_t}=X_t = 0\\ \large 3:\ \ \ for\ M=1\ to\ M\ do\\ \large 4:\ \ \ \ \ \ sample\ x_t^{[m]}\sim p(x_t|u_t, x_{t-1}^{[m]})\\ \large 5:\ \ \ \ \ \ w_t^{[m]}=p(z_t|x_{t}^{[m]})\\ \large 6:\ \ \ \ \ \ \bar{X_t}=\bar{X_t}+\lt x_t^{[m]},w_t^{[m]}\gt\\ \large 7:\ \ \ endfor\\ \large 8:\ \ \ for\ M=1\ to\ M\ do\\ \large 9:\ \ \ \ \ \ draw\ i\ with \ probalility\propto w_t^{[i]}\\ \large 10:\ \ \ \ add\ x_t^{[i]}\ to\ X_t\\ \large 11:\ endfor\\ \large 12:\ return\ X_t $

• 図3の右最上段のパネルに示すように、一般に条件分布「$p(x_t|y_{1:j})$」は、非ガウス分布なので、その分布は何らかの形での近似が必要である。
• 代表的な方法は、大別して以下の二つに分けられる。
• 解析的表現の補助利用のあるもの
• 数値的表現のみ

• Byte0: 73
• Byte1: 40
• Byte2: 64
• bit1: 1=人の動きを検知、0=人の動きを不検知
• Byte3: 00
• Byte4: 42
• 動体を検出してからの時間（秒）
• Byte5: 01
• bit6: 明るい
• bit7: 暗い

• こんな問題を考えよう。移動している複数の物体があって、それらを追跡したいとする。
• 追跡する物体は複数のジェット戦闘機や複数のミサイルかもしれないし、球場でクリケットをしている人々かもしれない。
• なんでもいい。これまでに学んだフィルタにこの問題を扱えるものはあるだろうか？
• 残念ながら、どのフィルタも理想的ではない。この問題の特徴を考えよう:
• 多峰性 (マルチモーダル): ゼロ個、一つ、あるいは複数の物体を同時に追跡したい。
• 遮蔽 (オクルージョン): 物体が他の物体を隠し、複数のオブジェクトから観測値が一つしか得られない場合がある。
• 非線形な観測値: レーダーは物体への距離を観測するので、(x,y,z) の座標に変換するとき非線形な平方根が必要になる。
• ガウス分布に従わないノイズ: 背景が大きく変化すると、コンピュータービジョンのアルゴリズムが背景を物体と誤認する場合がある。
• 連続: 物体の位置と速度 (つまり状態) は時間の経過とともに滑らかに変化する。
• 多変量: 位置、速度、旋回速度といった複数の属性を追跡したい。
• 未知のプロセスモデル: 系のプロセスモデルを私たちは知らない。

• 離散ベイズフィルタ:
• 離散ベイズフィルタは最も多くの制約に対応できる。
• 多峰性であり、非線形な観測値を処理でき、非線形な振る舞いに対応できるよう拡張することもできる。
• しかし、このフィルタは離散的かつ単変量である。
• カルマンフィルタ:
• カルマンフィルタはガウス分布に従うノイズを持った単峰性の線形系に対して最適な推定値を生成する。
• この問題ではどの仮定も成り立たない。
• 無香料カルマンフィルタ:
• UKF は非線形、連続、多変量の問題に対応できる。しかし多峰性ではなく、遮蔽を処理することもできない。
• また UKF はガウス分布と少しだけ異なる分布に従うノイズなら処理できるが、ノイズがガウス分布から大きく離れている場合や問題が非常に非線形な場合は上手く処理できない。
• 拡張カルマンフィルタ:
• EKF は基本的に UKF と同じ強みと弱みを持つ。
• ただし EKF はガウス分布に従わないノイズと非線形性からの影響をさらに受けやすい。

• 無香料カルマンフィルタの章で、ガウス分布を非線形な系で変換したときの効果を説明するために次のような図を示した:

非線形な変換を適用したガウス分布

• 左のグラフは次のガウス分布に独立に従う 3,000 個の点を示す:

$\Large \displaystyle \mu= \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \Sigma= \begin{bmatrix} 32 & 15 \\ 15 & 40 \end{bmatrix} $

• 右のグラフはそれらの点を次の関数に通したときに得られる点を示す:

$\Large \displaystyle x=x+y$

$\Large \displaystyle y=0.1x^2+y^2$

• ランダムに生成された有限個の点 (サンプル) を使って結果を計算する手法をモンテカルロ法 (Monte Carlo method, MC) と呼ぶ。
• アイデアは簡単だ: 状態を表現するサンプル (点) を十分多く生成し、モデル化している系にそれらの点を通し、変換された点から結果を計算する。
• 簡単に言えばこれが粒子フィルタの処理である。本書を通して使ってきたベイズフィルタのアルゴリズムが数千個の粒子に適用される。
• 一つ一つの粒子は系の状態の可能性を表す。数千個の粒子から状態の推定値を取り出すときは重み付き統計量が使われる。

• 大量の粒子をランダムに生成する。
• 粒子には位置や向きといった推定したい状態変数を好きなだけ持たせて構わない。
• それぞれの粒子は重みを持ち、系の実際の状態と粒子の状態が合致する可能性が重みによって (確率的に) 表される。
• 最初の重みは全て同じとする。

• 粒子の次の状態を予測する。
• あなたが予測した実際の系の振る舞いに基づいて、全ての粒子を移動させる。

• 観測値に基づいて粒子の重みを更新する。
• 観測値に近い粒子には観測値から遠い粒子よりも大きい重みが割り当てられる。

• 確率が非常に低い粒子を破棄して、確率の高い粒子のコピーで置き換える。

• もし必要なら、粒子の重み付き平均と重み付き共分散行列を計算して状態の推定値を得る。