機械と学習する

統計解析、機械学習について学習したことをまとめていきます

トップ > 統計モデリング > 系列データの統計モデリング　〜パーティクルフィルタによる状態推定〜

系列データの統計モデリング　〜パーティクルフィルタによる状態推定〜

統計モデリング 系列データ パーティクルフィルタ

【概要】

  • パーティクルフィルタを利用して系列データのモデリングと状態推定をやってみます

【目次】

はじめに

統計モデリングを行うにあたって、独立同分布(i.i.d; independent and identically distributed)を仮定することが多いと思います。

しかし、i.i.dを仮定できないケースはよくあります。時系列などの系列データが代表的です。また、i.i.dを仮定することが多い対象でも、計測対象や計測装置の経時的な変化（劣化）など本来的には系列を考慮しないといけない場面はあると思います。

系列データの解析について、主に「予測にいかす統計モデリングの基本」を参考に確認しましたので、数回に分けてまとめていこうと思います。なお、状態空間モデル(State Space Model)を使ったモデリングが対象です。私個人的には状態空間モデルでのモデリング経験はあるのですが、平滑化やパラメータ推論あたりなどの理解が曖昧でしたので、このあたりの整理をメインに考えています。

今回は第1回として、状態空間モデルをパーティクルフィルタで推論してみます。

間違いや勘違いなど、なにかありましたら指摘いただけるとすごく助かります。

No. リンク
第1回 系列データの統計モデリング 〜パーティクルフィルタによる状態推定〜 - 機械と学習する
第2回 系列データの統計モデリング 〜ハイパーパラメータの推定〜 - 機械と学習する
第3回 系列データの統計モデリング 〜自己組織型状態空間モデルを利用した状態推定モデルのハイパーパラメータ推定〜 - 機械と学習する
第4回 系列データの統計モデリング 〜カルマンフィルタ の導出と実装〜 - 機械と学習する

【トップに戻る】

状態空間モデル

系列データ

以下のデータ\left( Y = \{y_1, \cdots, y_T \} \right)が与えられ、このデータを基になんらかの統計モデルを設計したいとします。

f:id:hippy-hikky:20201019095703p:plain
時系列データの例

大きく分けて二通りのアプローチが考えられると思います。

  1. 回帰モデルの設計
  2. 時系列データのモデリング

1のアプローチでは、横軸（ここでは時間t）を特徴として、データ全体になんらかの関数を当てはめます。例えば上記データは3次関数のフィッティングができそうな雰囲気を感じます（山が二つあるように見えるので）。このアプローチの場合、Yの前後関係に依存関係を仮定せず、それぞれのデータは独立であると仮定します。この場合の同時確率は次のようになります。


\begin{align}
 p(Y) = \prod^T_{t=1} p(y_t)
\end{align}

2のアプローチでは、Yの時間的な依存関係（順番の影響）を仮定し、時間tの値自体には興味を持ちません。この場合の同時確率は次のようになります。


\begin{align}
 p(Y) = \prod^T_{t=1} p(y_t | y_{1:t})
\end{align}

図の横軸が気温などなんらかの量の場合は1のアプローチが効きそうです。しかし時間の場合、「時間」の値はデータのインデックスなども含むので、この値を特徴として使うのは違うはずです。今回は、アプローチ2の時系列データとして統計モデルを設計することを考えます。

系列データのモデリング

アプローチ2を採用することに決めましたが、モデリングを考えていくにあたって、アプローチ1との比較を考えながら進めてみます。

グラフィカルモデルで回帰モデル（アプローチ1）を書いてみると以下のようになります。

f:id:hippy-hikky:20201019103758p:plain
回帰モデルのグラフィカルモデル

ここで、\tilde{x}=\mathrm{f}(x, \theta)は通常は明示しないですが、パラメータθと特徴量xから決定的に算出される値です。回帰モデルの場合は、\tilde{x}に観測ノイズが加わってy_nが観測されるとします。

一方、時系列データなど時間的な依存関係があるデータに対して、次の構造のモデルを考えます。

f:id:hippy-hikky:20201019094613p:plain
状態空間モデル

このようなモデルを「状態空間モデル(State Space Model)」と呼びます。x_tを観測できない内部状態として、x_tに観測ノイズが付加されてy_tが観測されるとします。先の回帰モデルと対比すると、内部状態が独立であるか、時間依存性があるかが違うものと理解できそうです。

上記のグラフィカルモデルから、状態空間モデルを定式化すると次のようになります。


\begin{align}
  x_t &= p(x_t | x_{t-1}) \\
  y_t &= p(y_t | x_t) 
\end{align}

x_tの式は「システムモデル」、y_tのモデルは観測モデルと呼ばれます。なお、内部状態x_tには、過去の状態を含めることができるので、1時刻前の状態に”のみ”依存するというモデルでは無いことに注意です。

この状態空間モデルの同時確率を考えます。


\begin{align}
   p(y_{1:T}, x_{0:T}) &= p(y_T | y_{1:T-1}, x_{0:T}) p(y_{1:T-1}, x_{0:T}) \\
      &= p(y_T | y_{1:T-1}, x_{0:T}) p(x_T | y_{1:T-1}, x_{0:T-1}) p(y_{1:T-1}, x_{0:T-1}) \\
      &= \prod^T_{t=1} p(y_t | y_{1:t-1}, x_{0:t}) p(x_t | y_{1:t-1}, x_{0:t-1}) \\
      &= \prod^T_{t=1} p(y_t | x_{t}) p(x_t | x_{t-1})
\end{align}

1行目から3行目までは条件付き確立と同時確率の基本的jな性質から展開します。最後の式は上記のグラフィカルモデルに基づいて依存関係がある変数だけを残した結果です。

【トップに戻る】

状態空間モデルの状態推定

これらのモデルから、予測したいものは内部状態x_tであることがわかります。ここではx_tの推論について考えていきます（ p(x_{1:T} | y_{1:T})）。

p(x_{1:T} | y_{1:T})を考えて行くにあたって、事後周辺分布 p(x_{t} | y_{1:t'})を考えます。ここはあまり詳しく理解できてないのですが、x_{1:T}をそのまま扱うよりも計算が楽になるであろうということは想像できます。

事後周辺分布 p(x_{t} | y_{1:t'})は、t'の取り方で3種類の事後分布を考えることができます。


\begin{align}
   p(x_t | y_{1:t-1}) \cdots (1) \\
   p(x_t | y_{1:t}) \cdots (2) \\
   p(x_t | y_{1:T}) \cdots (3)
\end{align}

(1)は予測分布、(2)はフィルタ分布、(3)は平滑化分布と呼ばれています。平滑化分布については、次回で扱うとして、今回は予測分布とフィルタ分布を考えていきます。

予測分布とフィルタ分布を、確率の基本性質と、グラフィカルモデルを利用して展開すると、以下の図のようになります（手書きでごめんなさい）。この辺りは（全体的にですが）、参考文献[1]に丁寧に書かれているので是非そちらを参照してください。

f:id:hippy-hikky:20201019120943j:plain:h200

展開結果を確認すると、予測分布とフィルタ分布は互いに逐次更新していけば良いことがわかります。

f:id:hippy-hikky:20201019121014p:plain

【トップに戻る】

パーティクルフィルタ

前章で事後分布の導出をみていきました。しかし、実際にこのモデルを使って推論を行うにあたって、確率分布の表現と積分の計算をどのように実現するのかという問題があります。 ここではこの問題を、パーティクルフィルタ(Particle Filtering)を利用して近似推論していきます。パーティクルフィルタは、粒子フィルタ、モンテカルロフィルタなどとも呼ばれます。 パーティクルフィルタについては、参考文献[1]にも詳しく書かれていますが、[3]では疑似コードなども書かれてますし説明もとてもわかりやすいのでお勧めです。

パーティクルフィルタとは

パーティクルフィルタは、確率分布の表現を確率分布に従って得られたサンプル（パーティクル）の集合として表現します。サンプルによる確率分布表現のイメージとしては、ヒストグラムを思い浮かべてもらえば良いのかなと思います。

パーティクルの集合を使った分布の表現は次のようになります。

  • 予測分布（p(x_t | y_{1:t-1})）

\begin{align}
   p(x_t | y_{1:t-1}) \simeq \frac{1}{N} \sum^{N}_{i=1} \delta(x_t - x^{(i)}_{t|t-1})
\end{align}
  • フィルタ分布（p(x_t | y_{1:t})）

\begin{align}
   p(x_t | y_{1:t}) \simeq \frac{1}{N} \sum^{N}_{i=1} \delta(x_t - x^{(i)}_{t|t})
\end{align}

ここで、Nはパーティクルの数（ハイパーパラメータ）です。 x^{(i)}_{t|t-1}はi番目の予測パーティクルを示しており、x^{(i)}_{t|t}はフィルタパーティクルを表しています。\delta(\cdot)はデルタ関数を表しています。

次に、どのようにして予測パーティクル、フィルタパーティクルを得るかということですが、状態空間モデルとして定義したシステムモデルと観測モデルを利用します。

予測パーティクルは、1時刻前のフィルタ粒子（x^{i}_{t-1 | t-1}）をシステムモデル（p(x_t | x_{t-1})）に従って動かします。こうすることで予測パーティクルの集合X_{t|t-1}を得ることができます。（この部分の証明は参考文献[1]に丁寧に書かれているので参考にしてください）

次に、フィルタ分布p(x_t | y_{1:t})については、確率の基本性質を利用して展開します。


\begin{align}
   p(x_t | y_{1:t}) \propto p(y_t | x_t)  p(x_t | y_{1:t-1})
\end{align}

ここで、 p(x_t | y_{1:t-1})について、時刻tの状態は予測パーティクルの集合X_{t|t-1}で表現されています。各パーティクルでフィルタ分布を考えるとつぎのようになります。


\begin{align}
   p(x_t=x^{(i)}_{t|t-1} | y_{1:t}) \propto p(y_t | x^{(i)}_{t|t-1})  p(x^{(i)}_{t|t-1} | y_{1:t-1})
\end{align}

ここで、右辺のp(x^{(i)}_{t|t-1} | y_{1:t-1})について、パーティクルはそれぞれ重みみたいな物は持っておらず全て平等です（ヒストグラムをとるときに、データの重みは同じですよね。同じビン幅の中に何個データがあったかをカウントするわけなので。というイメージを持ってます。）。そのため、ここは\frac{1}{N}となります。 一方、p(y_t | x^{(i)}_{t|t-1})は、パーティクルの尤度となり、観測モデルに基づいて計算できます。

すると、


\begin{align}
   p(x_t=x^{(i)}_{t|t-1} | y_{1:t}) &\propto p(y_t | x^{(i)}_{t|t-1})  p(x^{(i)}_{t|t-1} | y_{1:t-1}) \\
     &=\frac{p(y_t | x^{(i)}_{t|t-1})}{\sum p(y_t | x^{(j)}_{t|t-1})} = w^{i}_t
\end{align}

となり、[tex:w^{i}t]がi番目のパーティクルの重みとなります。 パーティクルはそれぞれ平等であるため、[tex:w^{i}t]に従って復元抽出することで、フィルタ分布を得ることができます。

パーティクルフィルタのアルゴリズム

以上をアルゴリズムとして表現すると次のようになります。

  • X_{0|0} 初期分布を適当に設定
  • for t in (1,...,T) データ毎に繰り返し
    • for i in (1,...,N) パーティクル毎に繰り返し
      • x^{i}_{t|t-1} \sim p(x^{i}_{t|t-1}) | x^{i}_{t-1|t-1}, \theta), システムモデルに基づいたサンプル
      • w^{i}_{t} \propto p(y_t | x^{i}_{t|t-1}), 重み付け
    • X_{t|t} \leftarrow resampling(X_{t|t-1}, w_t), 重みに基づいた復元抽出

アルゴリズムとしてはとてもシンプルです。

【トップに戻る】

実装

実際に単純な時系列データの状態推定をしてみます。使うデータは冒頭に示した1次元の時系列データです。

システムモデルは、1時刻前の状態にガウスノイズが付加された物とします。観測ノイズはガウスノイズとします。（こういうモデルを1階差分モデルと呼ぶという理解で良いんですよね？）


\begin{align}
  x_t &\sim p(x_t | x_{t-1}) = \mathcal{N}(x_t | x_{t-1}, \alpha^2\sigma^2) \\
  y_t &\sim p(y_t | x_t)  = \mathcal{N}(y_t | x_{t}, \sigma^2)
\end{align}

\alpha, \sigmaはパラメータで、ここでは既知とします（これらのパラメータの推論については次回以降で扱います）。

このモデルに従った実装をgistにて公開していますので、必要な方は参考にしてもらえたらと思います。 ここでは、結果だけを掲載します。

f:id:hippy-hikky:20201020223207p:plain
パーティクルフィルタによる推論結果．青は観測値，オレンジは真の状態，緑の点がパーティクルの集合を示し，緑の太線はパーティクルの中心を表す．

パーティクルの集合がオレンジの線をだいたい覆っていることがわかりますね。

内部で利用しているParticleFilterクラスは次の通りです。

【トップに戻る】

おわりに

ということで、パーティクルフィルタを実装して時系列データの推論をやってみました。

パーティクルフィルタは証明もそんなに複雑じゃなく、実装もシンプルで使いやすいです。こんなにシンプルなのに、実際に使っても強力な印象ですごく面白いアルゴリズムだと思います。（最近使われているのかは、サーベイ不足で知らないですけど）

次は、既知と仮定したパラメータの推論に焦点を当てて行きたいと思います。

【トップに戻る】

参考資料

【トップに戻る】