萩原淳一郎，基礎からわかる時系列分析 ―Rで実践するカルマンフィルタ・MCMC・粒子フィルターをざっと読んで状態空間モデルを勉強した．

線形・ガウス型の状態空間モデルの場合はカルマンフィルタで解析的に解けるが，一般状態空間モデルの場合は，フィルタリング分布・予測分布・平滑化分布の積分計算を必ずしも解析的に解けない．そこで，MCMCを用いた一括解法や，粒子フィルタによる逐次解法が提案されている．本書のレポジトリは，hagijyun/tsbook - GitHub．

カルマンフィルタ，MCMC，粒子フィルタの関係を理解することができた．次はPRMLの下巻を読み返して，MCMCの理論的背景を押さえておきたい．

確率・統計に関する基礎

  • 周辺化消去は，同時分布から外したい変数の影響を，あらかじめ平均的に考慮することに相当する．なるほど，腑に落ちた．
  • 定常過程と非定常過程．
    • 強定常：確率過程の確率分布が時点\(t\)によらず一定．
    • 弱定常：確率過程の平均・分散・自己共分散・自己相関係数が時点\(t\)によらず一定．

時系列分析ひとめぐり

  • 五数要約：最小値・25%値・中央値・75%値・最大値．Rだとsummary()で，Python（pandas）だとDataFrame.describe()で確認できる．
  • 自己相関係数．\(\frac{\mathrm{Cov}[X_{t}, X_{t-k} ]}{\sqrt{\mathrm{Var}[X_t]}\sqrt{\mathrm{Var}[X_{t-k}]}}\)．
  • 時系列データの成分分解：レベル＋傾き＋周期．
  • ホルト・ウィンタース法．指数加重型移動平均（EWMA：Exponentially Weighted Moving Average）の一種．Rでは，HoltWinters()関数で実施できる．
  • 残差を確認して，モデルの妥当性を分析する．残差の自己相関により，周期性を上手く捉えられたか確認できる．
  • 状態空間モデルの解法の分類：
    • 線形・ガウス型状態空間モデル：
      • 一括型解法：ウィナーフィルタ．
      • 逐次型解法：カルマンフィルタ．
    • 一般状態空間モデル：
      • 一括型解法：MCMCを利用した解法．
      • 逐次型解法：粒子フィルタ．

状態空間モデル

  • 方程式による表現：
    • 状態方程式1：\(\mathbf{x}_t = f ( \mathbf{x}_{t-1}, \mathbf{w}_{t-1})\)
    • 観測方程式：\(y_t = h ( \mathbf{x}_{t}, v_{t})\)
  • 確率分布による表現：
    • 状態方程式に相当：\(p(\mathbf{x}_t \mid \mathbf{x}_{0:t-1}, y_{0:t-1}) = p(\mathbf{x}_t \mid \mathbf{x}_{t-1})\)
    • 観測方程式に相当：\(p(y_t \mid \mathbf{x}_{0:t}, y_{0:t-1}) = p(y_t\mid \mathbf{x}_{t})\)
  • 状態空間モデルの同時分布：\(p(\mathbf{x}_{0:T}, y_{0:T}) = p(\mathbf{x}_0) \prod_{t=1}^{T} p(y_t\mid \mathbf{x}_{t}) p(\mathbf{x}_t \mid \mathbf{x}_{t-1})\)
  • ARIMAモデルと状態空間モデルは密接な関係がある．

状態空間モデルにおける状態の推定

  • 状態の周辺事後分布：\(p(\mathbf{x}_{t'} \mid y_{1:t}) = \int p(\mathbf{x}_{0}, \mathbf{x}_{1}, \dots \mathbf{x}_{t'}, \dots \mid y_{1:t}) d\mathbf{x}_{t'以外}\)
    • \(t' < t\)：平滑化分布
    • \(t' = t\)：フィルタリング分布
    • \(t' > t\)：予測分布
  • フィルタリング分布：\(p(\mathbf{x}_{t} \mid y_{1:t}) = p(\mathbf{x}_{t} \mid y_{1:t-1}) \frac{p(y_{t} \mid \mathbf{x}_{t})}{p(y_{t} \mid y_{1:t-1})}\)
    • モデルの状態方程式\(p(y_{t} \mid \mathbf{x}_{t})\)2と，以下の一期先予測分布\(p(\mathbf{x}_{t} \mid y_{1:t-1})\)，一期先予測尤度\(p(y_{t} \mid y_{1:t-1})\)からなる．直感的には，一期先予測分布を観測方程式で補正して，一期先予測尤度で正規化したもの．
    • 一期先予測分布： \(p(\mathbf{x}_{t} \mid y_{1:t-1}) = \int p(\mathbf{x}_{t} \mid \mathbf{x}_{t-1}) p(\mathbf{x}_{t} \mid y_{1:t-1}) d\mathbf{x}_{t-1}\)
    • 一期先予測尤度3：\(p(y_{t} \mid y_{1:t-1}) = \int p(y_t \mid \mathbf{x}_{t} ) p(\mathbf{x}_{t} \mid y_{1:t-1}) d\mathbf{x}_{t-1}\)
  • 予測分布：
    • 定式化：\(p(\mathbf{x}_{t+k} \mid y_{1:t}) = \int p(\mathbf{x}_{t+k} \mid \mathbf{x}_{t+k-1}) p(\mathbf{x}_{t+k-1} \mid y_{1:t}) d\mathbf{x}_{t+k-1}\)
    • 一期先予測分布をもとに，状態方程式に基いて時間順方向への遷移を\(k-1\)解繰り返したもの．
  • 平滑化分布：
    • 定式化：\(p(\mathbf{x}_{t} \mid y_{1:T}) = p(\mathbf{x}_{t} \mid y_{1:t}) \int \frac{p(\mathbf{x}_{t+1} \mid \mathbf{x}_{t})}{p(\mathbf{x}_{t+1} \mid y_{1:t})} p(\mathbf{x}_{t+1} \mid y_{1:T}) d \mathbf{x}_{t+1}\)
    • \(p(\mathbf{x}_{t} \mid y_{1:t})\)はフィルタリング分布，\(p(\mathbf{x}_{t+1} \mid \mathbf{x}_{t})\)は状態方程式，\(p(\mathbf{x}_{t+1} \mid y_{1:t})\)は予測分布，\(p(\mathbf{x}_{t+1} \mid y_{1:T})\)は1時点先の平滑化分布．
    • 平滑化分布は，フィルタリング分布を再帰的に補正したもの．
  • パラメータ集合を\(\mathbf{\theta}\)とすると，対数尤度は\(l(\mathbf{\theta}) = \sum_{t=1}^{T} \mathrm{log} p( y_t \mid y_{1:t-1} ; \mathbf{\theta})\)．
    • 一期先予測尤度が含まれているため，特定のデータに基づきながら，過去の情報のみから見洗の値を当てる際の予測能力を踏まえた指標になっている．
    • この式には状態が登場しないため，周辺尤度あるいは対数周辺尤度と呼ばれることがある．
  • 状態空間モデルにおけるパラメータの扱い．
    • \(\mathbf{\theta}\)が確率変数でない：\(\mathrm{argmax}_{\mathbf{\theta}} l(\mathbf{\theta})\)を準ニュートン法やEMアルゴリズムで求める．
    • \(\mathbf{\theta}\)が確率変数：\(\mathrm{argmax}_{\mathbf{\theta}} p(\mathbf{\theta} \mid y_{1:t}) \propto \mathrm{argmax}_{\mathbf{\theta}} \{ p(y_{1:t} \mid \mathbf{\theta} ) + p(\mathbf{\theta} ) \}\)をMAP推定で求める．
  • パラメータを確率変数として扱う場合は，たとえもともと線形ガウス型の状態空間モデルであっても，一般状態空間モデルになる．つまり，ウィナーフィルタやカルマンフィルタでは直接解を求められない．

線形・ガウス型状態空間モデルの一括解法

  • ウィナーフィルタは定常な時系列を前提としている．
  • ウィナーフィルタはスペクトル分解を使って定式化できる．詳細は割愛する．

線形・ガウス型状態空間モデルの逐次解法

  • カルマンフィルタは，推定すべきデータの真値と点推定値の間の平均二乗誤差を最小にする逐次推定法である．
    • 状態方程式：\(\mathbf{x}_t = \mathbf{G}_t \mathbf{x}_{t-1} + \mathbf{w}_t\)，\(\mathbf{w}_t \sim \mathcal{N} (\mathbf{0}, \mathbf{W}_t)\)
    • 観測方程式：\(y_t = \mathbf{F}_t \mathbf{x}_t + v_t\)，\(v_t \sim \mathcal{N} (0, V_t)\)
    • 初期状態：\(\mathbf{x}_0 \sim \mathcal{N} (\mathbf{m}_0, \mathbf{C}_0)\)
    • パラメータ：\(\mathbf{\theta} = \{\mathbf{G}_t, \mathbf{F}_t, \mathbf{W}_t, V_t, \mathbf{m}_0, \mathbf{C}_0 \}\)
  • 上記の前提でフィルタリング分布・予測分布・平滑化分布を計算すると，逐次更新式を得る．詳細は割愛する．

一般状態空間モデルの一括解法

  • フィルタリング分布，予測分布，および平滑化分布には多次元の積分計算が含まれる．線形・ガウス型状態空間モデルの場合は，幸運にも 解析的に解けたが，一般状態空間モデルの場合は必ずしも解析的に解けないため，工夫が必要：
    • グリッドフィルタ：数値積分を活用するアプローチ．
    • 重点サンプリング：ラプラス近似を利用して，一般状態空間モデルにおける分布を線形・ガウス型状態空間モデルにおける分布で近似するアプローチ．
    • MCMC（マルコフ連鎖モンテカルロ法；Markov chain Monte Carlo）：求めたい分布からのサンプルを直接生成するアプローチ．
  • ギブス法によるMCMC：
    1. \(\mathbf{\theta^{(0)}}\)を初期化する．
    2. FOR \(i = \{1, \dots, I \}\):
      1. \(p ( \mathbf{x}_{0:T} \mid \mathbf{\theta} \! = \! \mathbf{\theta}^{(i-1)}, y_{1:T})\)から\(\mathbf{x}_{0:T}^{(i)}\)をサンプリングする．
      2. \(p ( \mathbf{\theta} \mid \mathbf{x}_{0:T} \! = \! \mathbf{x}_{0:T}^{(i)}, y_{1:T})\)から\(\mathbf{\theta}^{(i)}\)をサンプリングする．
  • 1-1の試行において，線形・ガウス型状態空間モデルに従う場合は，FFBS（前向きフィルタ後ろ向きサンプリング；Forward Filtering Backward Sampling）で効率的にサンプリングできる．
  • FFBSアルゴリズムは，時間順方向にカルマンフィルタをかけたあとで，時間逆方向にカルマン平滑化に基づいたサンプリングを行う．シミュレーション平滑化と呼ばれることもある．

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

  • カルマンフィルタの拡張：
    • 拡張カルマンフィルタ
    • アンセンテッドカルマンフィルタ
    • アンサンブルカルマンフィルタ
    • 粒子フィルタ：本書で紹介する．
  • 重点サンプリング
    • \(p(x)\)からサンプリング可能な場合：\(\mathrm{E}[ f(\mathbf{x})] \approx \frac{1}{N}\sum_{n=1}^{N} f (\mathbf{x}^{(n)})\)
    • \(p(x)\)の代わりに\(q(x)\)からサンプリングする場合：\(\mathrm{E}[ f(\mathbf{x})] \approx \frac{1}{N}\sum_{n=1}^{N} w(\mathbf{x}^{(n)}) f (\mathbf{x}^{(n)})\)
    • ただし，\(w(\mathbf{x}^{(n)}) = \frac{p(\mathbf{x}^{(n)})}{q(\mathbf{x}^{(n)})}\)を満たし，重点関数と呼ばれる．
  • 粒子フィルタリング：
    1. 時点\(t-1\)のフィルタリング分布：\(\{ \mathbf{x}_{t-1}^{(n)}, w_{t-1}^{(n)} \}_{n=1}^N\)
    2. 時点\(t\)において，FOR \(n = \{1, \dots, N \}\)：
      1. 実現値を更新：\(q \left( \mathbf{x}_{t} \mid \mathbf{x}_{0:t-1}^{(n)}, y_{1:t} \right)\)から\(\mathbf{x}_{t}^{(n)}\)をサンプリング．
      2. 重みを更新：\(w_t^{(n)} \leftarrow \frac{p\left(\mathbf{x}_t^{(n)} \mid \mathbf{x}_{t-1}^{(n)} \right) p \left( y_t \mid y_{1:t-1}, \mathbf{x}_t^{(n)} \right)}{q\left( \mathbf{x}_t^{(n)} \mid \mathbf{x}_{0:t-1}^{(n)}, y_{1:t} \right)}\)．
    3. 重みの規格化：\(w_t^{(n)} \leftarrow \frac{w_t^{(n)}}{\sum_{n=1}^N w_{t}^{(n)}}\)
    4. \(w_t^{(n)}\)大きさに応じて，\(x_t^{(n)}\)を復元抽出．
    5. 時点\(t\)のフィルタリング分布：\(　\{ \mathbf{x}_t^{(n)}, w_{t}^{(n)} \}_{n=1}^{N}\)．
  • フィルタリング分布の一般式と，次のように対応している．
    • 一期先予測分布：\(w_{t-1}^{(n)} p \left(\mathbf{x}_t^{(n)} \mid \mathbf{x}_{t-1}^{(n)} \right)\)
    • 一期先予測尤度：重みの規格化．
    • フィルタリング分布：一期先予測分布\(w_{t-1}^{(n)} p \left(\mathbf{x}_t^{(n)} \mid \mathbf{x}_{t-1}^{(n)} \right)\)を\(p\left( y_t \mid y_{1:t-1}, \mathbf{x}_{t}^{(n)} \right)\)で補正4．
  • 提案分布\(q \left(\mathbf{x}_{t}^{(n)} \mid \mathbf{x}_{t-1}^{(n)} \right)\)として，状態方程式\(p \left(\mathbf{x}_{t}^{(n)} \mid \mathbf{x}_{t-1}^{(n)} \right)\)を用いると，重みの更新式は非常に簡単になる．この条件の粒子フィルタをブートストラップフィルタ，あるいはモンテカルロフィルタと呼ぶ．
  • 粒子数が有限の場合，リサンプリングにより近似性能の劣化を防ぐことができる．
  • 粒子予測：重みを更新せず，実現値\(x_{t+k}^{(n)}\)を状態方程式で逐次的に更新することで，未来の点を予測できる．これは，粒子フィルタリングにおいて観測値\(y_{t+k}\)が存在しない（欠測）場合と等価．
  • 粒子平滑化：
    • 北川アルゴリズム：過去の粒子に対して，現時点の規範でリサンプリングする．原理的には任意の過去の時点まで平滑化可能だが，リサンプリングによる粒子の退化が問題となる．SIR平滑化，固定ラグ平滑化，廉価版平滑化とも呼ばれる．
    • 前向きフィルタ後ろ向きシミュレーション（FFBS：Forward Filtering Backward Simulation）アルゴリズム：フィルタリング分布の\(w_t^{(n)}\)を再帰的に補正する．
  • 推定精度向上のためのテクニック：
    • 補助粒子フィルタ：現在の観測値の影響を考慮したリサンプリング処理が，実現値や重み更新より先に行われる．
    • リュウ・ウエストフィルタ：補助粒子フィルタの適用例．粒子フィルタを用いて，状態とパラメータをあわせて推定する．カーネル平滑化を用いて，パラメータの推定値の分散を増大させることなく実現値をリフレッシュする．
    • ラオ・ブラックウェル化：一般化線形モデルにおいて，部分的に線形・ガウス型状態空間モデルが成り立っているとき，カルマンフィルタと粒子フィルタを組み合わせる手法．例えば，線形・ガウス型状態空間モデルのパラメータが未知の場合などがこれに該当する．

一般状態空間モデルにおける応用的な分析例

  • 本章では，構造変化について扱う．
    • 変化点が既知：カルマンフィルタ
    • 変化点が未知：粒子フィルタ，MCMC
  • 未知の変化点を実時間で検出するために，粒子フィルタを用いる．状態雑音の分散に対する時変の倍率が特定の閾値を超えた際に，構造変化が発生したと考える．

  1. システム方程式とも呼ばれる． ↩

  2. 書籍中，\(p(y_t \mid \mathbf{x}_t)\)を尤度と呼ぶことがあり，一期先予測尤度と混同してハマった． ↩

  3. なぜ尤度と呼ぶんだろう．まだ腑に落ちてない． ↩

  4. 一般的な状態空間モデルの定式化では一期先予測分布を観測方程式で補正していた． ↩