「大人の教養・知識・気付き」を伸ばすブログ

一流の大人(ビジネスマン、政治家、リーダー…)として知っておきたい、教養・社会動向を意外なところから取り上げ学ぶことで“気付く力”を伸ばすブログです。データ分析・語学に力点を置いています。 →現在、コンサルタントの雛になるべく、少しずつ勉強中です(※2024年1月21日改訂)。

MENU

時系列解析の基礎(32/XX)

 以下の書籍

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

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について

  1. 時点\tau=tにおけるフィルタリング分布として、実現値と重みをそれぞれ\left\{\boldsymbol{x}_{t-1}^{(i)}, \omega_{t-1}^{(i)}\right\}^{n}_{i=1}とする。
  2. 時点\tau=\tau+1とし、i=1とする。
  3. 提案分布q\left(\boldsymbol{x}_t|\boldsymbol{x}_{0:t-1}^{(i)},y_{1:t}\right)から標本\boldsymbol{x}_{t}^{(i)}を得る。
  4. 重み\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}

    とする。
  5. i\lt nならばi=i+1として2.に戻る。そうでなければ6.に移る。
  6. i=1,\cdots,nについて\omega_{t}^{(i)}\leftarrow\displaystyle{\frac{\omega_{t}^{(i)}}{\displaystyle{\sum_{i=1}^{n}\omega_{t}^{(i)}}}}と規格化する。
  7. 集合\{1,\cdots,n\}から\omega_{t}^{(i)},i=1,\cdots,nに比例した確率で復元抽出を行ない、リサンプリング用のインデックス\boldsymbol{k}=\{k_1,\cdots,k_i,\cdots,k_n\}を構築する。
  8. \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)}とする。
  9. 重み\omega_t^{(i)}=\displaystyle{\frac{1}{n}},i=1,\cdots,nと置き直す。
  10. \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の情報が考慮されないため、有限の粒子では推定性能が劣化する可能性がある。
 リサンプリングは有限の粒子を効率的に再利用するための処置であり、粒子数が無限の場合には原理的に不要ではあるものの、粒子数が有限の場合には現実的に必要となる。リサンプリングすることでフィルタリング分布において重みの大きい粒子を優先的に積み増すことになるため、次の時点で重みが大きくなると見込まれる領域に多くの粒子が配置されることになり、分布の近似精度が劣化せずに済む。ただしリサンプリングは粒子の多様性を減らすため、過度に行うと推定精度を悪化させる。
 最後にブートストラップ・フィルタにおいて、リサンプリングを毎時点必ず行うとすれば、粒子の重みはすべて同じであり、フィルタリングにおいて一期先予測を行なってから観測値に基づく尤度で補正を施す仕組みがより直接的に表現できる。

プライバシーポリシー お問い合わせ