時系列解析の基礎(32/XX) - 「大人の教養・知識・気付き」を伸ばすブログ

　以下の書籍

基礎からわかる時系列分析―Rで実践するカルマンフィルタ・MCMC・粒子フィルタ― Data Science Library

作者:萩原淳一郎,瓜生真也,牧山幸史
技術評論社

Amazon

を中心に時系列解析を勉強していきます。

前回

power-of-awareness.com

前回
16.　一般状態空間モデルの逐次解法
- 16.1　粒子フィルタ
  - 16.1.1　粒子フィルタリング
  - 16.1.2　粒子フィルタリングの導出
次回

16.　一般状態空間モデルの逐次解法

16.1　粒子フィルタ

　順次データが得られる場合、都度都度 $\mathrm{MCMC}$ で状態を推定するのは計算効率が悪い。そのためカルマンフィルタを発展させた様々な改良手法が提案されている。ここでは粒子フィルタ(逐次モンテカルロ法、 $\mathrm{SIR}$ 法とも呼ばれる。)を導入する。
　粒子フィルタでは事後分布の近似を直接求めた後で、必要に応じて各種統計量を算出する。

16.1.1　粒子フィルタリング

　以降では、重点サンプリングの発展版として粒子フィルタリングを導入する。
　重点サンプリングは、求めたい分布に関する当家利用を算出する際に、サンプリングが容易な代理分布(提案分布という。)から一旦標本を抽出し、求めたい分布の尤度で補正を施す。たとえば確率ベクトル $\boldsymbol{x}$ の確率密度を $p(\boldsymbol{x})$ とするとき、ある関数 $f(\boldsymbol{x})$ の期待値は

$\begin{aligned} E\left[f(\boldsymbol{x})\right]=\displaystyle{\int f(\boldsymbol{x})p(\boldsymbol{x})}d\boldsymbol{x} \end{aligned}$

で与えられ、 $p(\boldsymbol{x})$ 空のサンプリングが容易ならば、

$\begin{aligned} E\left[f(\boldsymbol{x})\right]\approx\displaystyle{\frac{1}{n}\sum_{i=1}^{n}f\left(\boldsymbol{x}^{(n)}\right)} \end{aligned}$

が得られる。しかし $p(\boldsymbol{x})$ からのサンプリングが容易だとは限らない。そこで代わりに提案分布 $q(\cdots)$ を用いることにすれば、

$\begin{aligned} E\left[f(\boldsymbol{x})\right]&=\displaystyle{\int f(\boldsymbol{x})p(\boldsymbol{x})}d\boldsymbol{x}\\ &=\displaystyle{\int f(\boldsymbol{x})\frac{p(\boldsymbol{x})}{q(\boldsymbol{x})}q(\boldsymbol{x})}d\boldsymbol{x} \end{aligned}$

において、 $w(\boldsymbol{x})=\displaystyle{\frac{p(\boldsymbol{x})}{q(\boldsymbol{x})}}$ とおけば、

$\begin{aligned} E\left[f(\boldsymbol{x})\right]&=\displaystyle{\int f(\boldsymbol{x})w(\boldsymbol{x})q(\boldsymbol{x})}d\boldsymbol{x}\\ &\approx\displaystyle{\frac{1}{n}\sum_{i=1}^{n}w\left(\boldsymbol{x}^{(n)}\right)f\left(\boldsymbol{x}^{(n)}\right)} \end{aligned}$

と近似できる。そこで提案分布からの標本を用いれば $E\left[f(\boldsymbol{x})\right]$ を近似することができることが分かった。ここで $w\left(\boldsymbol{x}\right)$ は重点関数と呼ばれ、代理である提案分布の影響を打ち消して求めたい分布の尤度で補正を施す意味を持つ。
　もし求めたい分布と大きく異なる提案分布から標本を発生させると、重点関数により尤度が低く補正されることで、近似計算において無駄が多くなる。したがって標本サイズが同じならば、提案分布 $q(\cdot)$ は求めたい分布 $p(\cdot)$ に似ている程、近似精度が向上する。
　この重点サンプリングの仕組みをフィルタリングに応用し、求めたいフィルタリング分布の離散近似における重みを逐次更新するようにまとめると、粒子フィルタリングのアルゴリズムが得られる。

16.1.2　粒子フィルタリングの導出

　フィルタリング分布 $p(\boldsymbol{x}_{1:t}|y_{1:t})$ の離散分布を得るべく、まず同時事後分布 $p(\boldsymbol{x}_{0:t}|y_{1:t})$ の離散近似を求め、それを周辺化することとする。同時事後分布 $p(\boldsymbol{x}_{0:t}|y_{1:t})$ の離散近似は、提案分布を $q(\boldsymbol{x}_{0:t}|y_{1:t})$ とした重点サンプリングで得ることとする。このとき同時事後分布と提案分布の比率 $\omega_t$ は、

$\begin{aligned} \omega_t&=\displaystyle{\frac{p(\boldsymbol{x}_{0:t}|y_{1:t})}{q(\boldsymbol{x}_{0:t}|y_{1:t})}}\\ &=\displaystyle{\frac{p(\boldsymbol{x}_{0:t}|y_t,y_{1:t-1})}{q(\boldsymbol{x}_{0:t}|y_{1:t})}}\\ &=\displaystyle{\frac{\displaystyle{\frac{p(\boldsymbol{x}_{0:t},y_t|y_{1:t-1})}{p(y_t|y_{1:t-1})}}}{q(\boldsymbol{x}_{0:t}|y_{1:t})}}\\ &\propto\displaystyle{\frac{p(\boldsymbol{x}_{0:t},y_t|y_{1:t-1})}{q(\boldsymbol{x}_{0:t}|y_{1:t})}}\\ &=\displaystyle{\frac{p(\boldsymbol{x}_t,\boldsymbol{x}_{0:t-1},y_t|y_{1:t-1})}{q(\boldsymbol{x}_{0:t}|y_{1:t})}}\\ &=\displaystyle{\frac{p(\boldsymbol{x}_t,y_t|\boldsymbol{x}_{0:t-1},y_{1:t-1})p(\boldsymbol{x}_{0:t-1}|y_{1:t-1})}{q(\boldsymbol{x}_{0:t}|y_{1:t})}} \end{aligned}$

を満たす。分母の提案分布に関して $\mathrm{Markov}$ 性を仮定して、 $q(\boldsymbol{x}_{0:t}|y_{1:t})=$ $q(\boldsymbol{x}_t|\boldsymbol{x}_{0:t-1},y_{1:t})$ $q(\boldsymbol{x}_{0:t-1}|y_{1:t-1})$ と分解できると仮定すれば、

$\begin{aligned} \omega_t&=\displaystyle{\frac{p(\boldsymbol{x}_t,y_t|\boldsymbol{x}_{0:t-1},y_{1:t-1})p(\boldsymbol{x}_{0:t-1}|y_{1:t-1})}{q(\boldsymbol{x}_{0:t}|y_{1:t})}}\\ &=\displaystyle{\frac{p(\boldsymbol{x}_t,y_t|y_{1:t-1})p(\boldsymbol{x}_{0:t-1}|\boldsymbol{x}_{0:t-1},y_{1:t-1})}{q(\boldsymbol{x}_t|\boldsymbol{x}_{0:t-1},y_{1:t})q(\boldsymbol{x}_{0:t-1}|y_{1:t})}}\\ &=\displaystyle{\frac{p(\boldsymbol{x}_t,y_t|\boldsymbol{x}_{0:t-1},y_{1:t-1})}{q(\boldsymbol{x}_t|\boldsymbol{x}_{0:t-1},y_{1:t})}\frac{p(\boldsymbol{x}_{0:t-1}|y_{1:t-1})}{q(\boldsymbol{x}_{0:t-1}|y_{1:t-1})}}\\ &=\displaystyle{\frac{p(\boldsymbol{x}_t,y_t|\boldsymbol{x}_{0:t-1},y_{1:t-1})}{q(\boldsymbol{x}_t|\boldsymbol{x}_{0:t-1},y_{1:t})}\omega_{t-1}}\\ &=\displaystyle{\frac{p(y_t|\boldsymbol{x}_t,\boldsymbol{x}_{0:t-1},y_{1:t-1})p(\boldsymbol{x}_t|y_{1:t-1})}{q(\boldsymbol{x}_t|\boldsymbol{x}_{0:t-1},y_{1:t})}\omega_{t-1}} \end{aligned}$

を得る。
　条件付き独立性

$\begin{aligned} p(\boldsymbol{x}_t|\boldsymbol{x}_{0:t-1},y_{1:t-1})&=p(\boldsymbol{x}_t|\boldsymbol{x}_{t-1}),\\ p(y_t|\boldsymbol{x}_{0:t},y_{1:t-1})&=p(y_t|\boldsymbol{x}_t) \end{aligned}$

より、

$\begin{aligned} \omega_t&=\displaystyle{\frac{p(y_t|\boldsymbol{x}_t,y_{1:t-1})p(\boldsymbol{x}_t|\boldsymbol{x}_{t-1})}{q(\boldsymbol{x}_t|\boldsymbol{x}_{0:t-1},y_{1:t})}\omega_{t-1}}\\ &=\displaystyle{\frac{p(y_t|\boldsymbol{x}_t,y_{1:t-1})p(\boldsymbol{x}_t|y_{1:t-1})}{q(\boldsymbol{x}_t|\boldsymbol{x}_{0:t-1},y_{1:t})}\omega_{t-1}} \end{aligned}$

が得られる。この式に $\omega_t$ の規格化を加えることでアルゴリズムの原型が完成する。この仕組みに基づき、提案分布からの標本に対する重みの計算を時間順に繰り返す。これにより、同時事後分布の離散近似

$\begin{aligned} \hat{p}(\boldsymbol{x}_{0:t}|y_{1:t})=\displaystyle{\sum_{i=1}^{n}\omega_{t}^{(i)}\delta_{\boldsymbol{x}_{0:t}^{(i)}}} \end{aligned}$

が得られる。ここで $\delta_{\cdot}$ は $\cdot$ における単位確率質量を表す。
　この離散近似に関する時点 $t$ での周辺化は、時点 $t$ における重みと実現値がフィルタリング分布の離散近似

$\begin{aligned} \hat{p}(\boldsymbol{x}_t|y_{1:t})=\displaystyle{\sum_{i=1}^{n}\omega_{t}^{(i)}\delta_{\boldsymbol{x}_t^{i}}} \end{aligned}$

で与えられる。そして実用性を増すべくリサンプリングを加えたアルゴリズムを用いる。

粒子フィルタリングのアルゴリズム　 $t=1,\cdots,T$ について

時点 $\tau=t$ におけるフィルタリング分布として、実現値と重みをそれぞれ $\left\{\boldsymbol{x}_{t-1}^{(i)}, \omega_{t-1}^{(i)}\right\}^{n}_{i=1}$ とする。
時点 $\tau=\tau+1$ とし、 $i=1$ とする。
提案分布 $q\left(\boldsymbol{x}_t|\boldsymbol{x}_{0:t-1}^{(i)},y_{1:t}\right)$ から標本 $\boldsymbol{x}_{t}^{(i)}$ を得る。
重み $\omega_{t-1}$ について、
$\begin{aligned}\omega_t^{(i)}\leftarrow\omega_{t-1}^{(i)}\displaystyle{\frac{p\left(\boldsymbol{x}_{t}^{(i)}|\boldsymbol{x}_{t-1}^{(i)}\right)p\left(y_t|y_{1:t-1}\boldsymbol{x}_{t}^{(i)}\right)}{q\left(\boldsymbol{x}_{t}^{(i)}|\boldsymbol{x}_{0:t-1}^{(i)},y_{1:t}\right)}}\end{aligned}$
とする。
$i\lt n$ ならば $i=i+1$ として $2.$ に戻る。そうでなければ $6.$ に移る。
$i=1,\cdots,n$ について $\omega_{t}^{(i)}\leftarrow\displaystyle{\frac{\omega_{t}^{(i)}}{\displaystyle{\sum_{i=1}^{n}\omega_{t}^{(i)}}}}$ と規格化する。
集合 $\{1,\cdots,n\}$ から $\omega_{t}^{(i)},$ $i=1,\cdots,n$ に比例した確率で復元抽出を行ない、リサンプリング用のインデックス $\boldsymbol{k}=$ $\{k_1,\cdots,k_i,\cdots,k_n\}$ を構築する。
$\boldsymbol{x}_t^{(\boldsymbol{k})}=$ $\left\{\boldsymbol{x}_t^{(k_1)}\right.,$ $\cdots,$ $\boldsymbol{x}_t^{(k_i)},$ $\cdots,$ $\left.\boldsymbol{x}_t^{(k_n)}\right\}$ に対してインデックスを振り直し実現値 $\boldsymbol{x}_t^{(i)}$ とする。
重み $\omega_t^{(i)}=\displaystyle{\frac{1}{n}},i=1,\cdots,n$ と置き直す。
$\tau\lt T$ ならば $1.$ に戻る。そうでなければ終了する。

　時点 $0$ におけるフィルタリング分布は事前分布に相当し、通常 $\boldsymbol{x}_0^{(i)}$ は事前分布から抽出された実現値とし、 $\omega_0^{(i)}=\displaystyle{\frac{1}{n}}$ とする。
　　このアルゴリズムは、フィルタリング分布の定式化と整合している。

1期先予測分布： $\omega_{t-1}^{(i)}p\left(\boldsymbol{x}_t^{(i)}|\boldsymbol{x}_{t-1}^{(i)}\right)$
1期先予測尤度：重みの規格化により考慮する
フィルタリング分布： $p\left(y_t|y_{1:t-1},\boldsymbol{x}_t^{(i)}\right)$ により補正する

　このアルゴリズムではところどころ留意点がある。
　まず $3.$ において実現値を抽出すべく $q(\cdot)$ が用いられる。提案分布は求めたい分布の代理であるから、その影響はキャンセルアウトされている。提案分布は好きなものが利用でき得るものの、求めたい分布の密度が $0$ でないところで $0$ になってはいけない。また粒子数は有限であるから、基本的に求めたい分布に可能な限り類似した分布を用いるべきである。
　重み $\omega_t^{(i)}$ のばらつきを最小にする最適な提案分布では、 $q(\cdot)$ は状態空間モデルにおける条件付き独立の性質から $\boldsymbol{x}_{t-1}$ および $y_t$ 飲みに依存する。このため

$\begin{aligned} q\left(\boldsymbol{x}_t^{(i)}|\boldsymbol{x}_{0:t-1}^{i},y_{1:t}\right)=p\left(\boldsymbol{x}_t^{(i)}|\boldsymbol{x}_{t-1}^{(i)},y_t\right) \end{aligned}$

が成り立つ。
　実際にどの分布を提案分布にするかは問題にも依存するが、良く用いられるのは状態方程式である。その場合、

$\begin{aligned} \omega_t^{(i)}&\leftarrow\omega_{t-1}^{(i)}\displaystyle{\frac{p\left(\boldsymbol{x}_{t}^{(i)}|\boldsymbol{x}_{t-1}^{(i)}\right)p\left(y_t|y_{1:t-1}\boldsymbol{x}_{t}^{(i)}\right)}{q\left(\boldsymbol{x}_{t}^{(i)}|\boldsymbol{x}_{0:t-1}^{(i)},y_{1:t}\right)}}\\ &=\omega_{t-1}^{(i)}\displaystyle{\frac{p\left(\boldsymbol{x}_{t}^{(i)}|\boldsymbol{x}_{t-1}^{(i)}\right)p\left(y_t|y_{1:t-1}\boldsymbol{x}_{t}^{(i)}\right)}{p\left(\boldsymbol{x}_t^{(i)}|\boldsymbol{x}_{t-1}^{(i)},y_t\right)}}\\ &=\omega_{t-1}^{(i)}\displaystyle{\frac{p\left(\boldsymbol{x}_{t}^{(i)}|\boldsymbol{x}_{t-1}^{(i)}\right)p\left(y_t|y_{1:t-1}\boldsymbol{x}_{t}^{(i)}\right)}{p\left(\boldsymbol{x}_t^{(i)}|\boldsymbol{x}_{t-1}^{(i)}\right)}}\\ &=\omega_{t-1}^{(i)}p\left(y_t|y_{1:t-1}\boldsymbol{x}_{t}^{(i)}\right) \end{aligned}$

である。この形態の粒子フィルタはブートストラップ・フィルタ、モンテカルロ・フィルタとも呼ばれる。
提案分布を状態方程式に設定する場合、観測値 $y_t$ の情報が考慮されないため、有限の粒子では推定性能が劣化する可能性がある。
　リサンプリングは有限の粒子を効率的に再利用するための処置であり、粒子数が無限の場合には原理的に不要ではあるものの、粒子数が有限の場合には現実的に必要となる。リサンプリングすることでフィルタリング分布において重みの大きい粒子を優先的に積み増すことになるため、次の時点で重みが大きくなると見込まれる領域に多くの粒子が配置されることになり、分布の近似精度が劣化せずに済む。ただしリサンプリングは粒子の多様性を減らすため、過度に行うと推定精度を悪化させる。
　最後にブートストラップ・フィルタにおいて、リサンプリングを毎時点必ず行うとすれば、粒子の重みはすべて同じであり、フィルタリングにおいて一期先予測を行なってから観測値に基づく尤度で補正を施す仕組みがより直接的に表現できる。

次回

https://power-of-awareness.com/entry/2023/03/27/050000power-of-awareness.com

前回

16. 一般状態空間モデルの逐次解法

16.1 粒子フィルタ

16.1.1 粒子フィルタリング

16.1.2 粒子フィルタリングの導出

次回

16.　一般状態空間モデルの逐次解法

16.1　粒子フィルタ

16.1.1　粒子フィルタリング

16.1.2　粒子フィルタリングの導出