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

状態空間モデルにおいて、観測値 y1, y2, .., yt が与えられたとき、状態 xi (i = 1, 2, ..., t-1, t, t+1, ...) を推定することができる。観測値 y1, y2, .., yt のもとで推定された、状態 xt の分布をフィルタ分布(フィルタリング分布)、xr (r < t) の分布を平滑化分布、xs (s > t) の分布を予測分布と呼ぶ。

フィルタ分布

いま、観測値 y1, y2, .., yt が与えられ、また状態方程式および観測方程式のパラメーターも既知とする。このとき、フィルタ分布は次のように求めることができる。

\[ p(\mathbf{x}_{t} | y_{1:t}) = p(\mathbf{x}_{t} | y_{t}, y_{t-1}) \]

ここで、条件付き確率の関係を用いて、

\[ p(\mathbf{x}_{t} | y_{1:t}) = \frac{p(\mathbf{x}_{t}, y_{t} | y_{1:t-1})}{p(y_{t}|y_{t-1})} = \frac{p(y_{t} | \mathbf{x}_{t}, y_{1:t-1}) p(\mathbf{x}_{t} | y_{1:t-1})}{p(y_{t}|y_{t-1})} \]

状態空間モデルにおいて、観測値 yt は、状態 xt のみに依存するので、

\[ p(\mathbf{x}_{t} | y_{1:t}) = \frac{p(y_{t} | \mathbf{x}_{t}) p(\mathbf{x}_{t} | y_{1:t-1})}{p(y_{t}|y_{t-1})} \]

このとき、p(xt | y1:t-1) を 1 期先予測分布、p(yt | yt-1) を 1 期先予測尤度という。

1 期先予測尤度と状態 xt の関係を示すために、状態 xt に関して周辺化を行う。(xt に関して周辺化を行うことで、xt が取りうるすべての値において、その値となる確率を積分しているので、xt に関して周辺化すると確率変数 xt に関しては 1 となる。よって次の式の 1 行目が成り立つ。)

\[ \begin{eqnarray} p(y_{t}|y_{t-1}) &=& \int p(y_{t}, \mathbf{x}_{t} | y_{1:t-1}) d\mathbf{x}_{t} \\ &=& \int p(y_{t} | \mathbf{x}_{t}, y_{1:t-1}) p(\mathbf{x}_{t} | y_{1:t-1}) d\mathbf{x}_{t} \\ &=& \int p(y_{t} | \mathbf{x}_{t}) p(\mathbf{x}_{t} | y_{1:t-1}) d\mathbf{x}_{t} \end{eqnarray} \]

以上の式変形により、観測値 y1, y2, .., yt が与えられたとき、フィルタ分布は次の式で表せるようになる。

\[ p(\mathbf{x}_{t} | y_{1:t}) = \frac{p(y_{t} | \mathbf{x}_{t}) p(\mathbf{x}_{t} | y_{1:t-1})}{\int p(y_{t} | \mathbf{x}_{t}) p(\mathbf{x}_{t} | y_{1:t-1}) d\mathbf{x}_{t} } \]

平滑化分布

観測値 y1, y2, .., yt が与えられ、また状態方程式および観測方程式のパラメーターも既知とする。このとき、状態 xr (r < t) における平滑化分布 p(xr | y1:t) を求めるには、時点 t までのフィルタ分布が求まり、時点 t から時間を遡って求めることになる。そのため、平滑化分布 p(xr | y1:t) と状態 xr+1 の関係を明らかにして、漸化式を作っていく必要がある。ここで、まず平滑化分布 p(xr | y1:t) に対して xr+1 で周辺化を行なってから、式変更して行なっていく。

\[ \begin{eqnarray} p(\mathbf{x}_{r} | y_{1:t}) &=& \int p(\mathbf{x}_{r}, \mathbf{x}_{r+1} | y_{1:t}) d\mathbf{x}_{r+1} \\ &=& \int p(\mathbf{x}_{r} | \mathbf{x}_{r+1}, y_{1:t}) p(\mathbf{x}_{r+1} | y_{1:t}) d\mathbf{x}_{r+1} \\ &=& \int p(\mathbf{x}_{r} | \mathbf{x}_{r+1}, y_{1:r}) p(\mathbf{x}_{r+1} | y_{1:t}) d\mathbf{x}_{r+1} \\ &=& \int \frac{p(\mathbf{x}_{r}, \mathbf{x}_{r+1} | y_{1:r})}{ p(\mathbf{x}_{r+1} | y_{1:r}) }p(\mathbf{x}_{r+1} | y_{1:t}) d\mathbf{x}_{r+1} \\ &=& \int \frac{p(\mathbf{x}_{r+1} | \mathbf{x}_{r}, y_{1:r}) p(\mathbf{x}_{r} | y_{1:r})}{p(\mathbf{x}_{r+1} | y_{1:r})}p(\mathbf{x}_{r+1} | y_{1:t}) d\mathbf{x}_{r+1} \\ &=& \int \frac{p(\mathbf{x}_{r+1} | \mathbf{x}_{r}) p(\mathbf{x}_{r} | y_{1:r})}{p(\mathbf{x}_{r+1} | y_{1:r})}p(\mathbf{x}_{r+1} | y_{1:t}) d\mathbf{x}_{r+1} \\ &=& p(\mathbf{x}_{r} | y_{1:r}) \int \frac{p(\mathbf{x}_{r+1} | \mathbf{x}_{r}) }{p(\mathbf{x}_{r+1} | y_{1:r})}p(\mathbf{x}_{r+1} | y_{1:t}) d\mathbf{x}_{r+1} \end{eqnarray} \]

予測分布

観測値 y1, y2, .., yt が与えられ、また状態方程式および観測方程式のパラメーターも既知とする。このとき、状態 xs (s > t) における予測分布 p(xs | y1:t) を求めるには、時点 t までのフィルタ分布が求まり、時点 t+1, t+2, ..., s-1, s のように逐次的に状態の予測分布を求めていく必要がある。そのため、予測分布 p(xs | y1:t) と状態 xs-1 の関係を明らかにして漸化式を作って必要がある。ここで、予測分布 p(xs | y1:t) に対して xs-1 で周辺化を行う。

\[ \begin{eqnarray} p(\mathbf{x}_{s} | y_{1:t}) &=& \int p(\mathbf{x}_{s}, \mathbf{x}_{s-1} | y_{1:t}) d\mathbf{x}_{s-1} \\ &=& \int p(\mathbf{x}_{s} | \mathbf{x}_{s-1}, y_{1:t}) p(\mathbf{x}_{s-1} | y_{1:t}) d\mathbf{x}_{s-1} \\ &=& \int p(\mathbf{x}_{s} | \mathbf{x}_{s-1}) p(\mathbf{x}_{s-1} | y_{1:t}) d\mathbf{x}_{s-1} \end{eqnarray} \]

References

  • 深谷 肇一. 状態空間モデルによる時系列解析とその生態学への応用. 日本生態学会誌. 2016, 66:375-389. DOI: 10.18960/seitai.66.2_375
  • 萩原 淳一郎, 瓜生 真也, 牧山 幸史, 石田 基広. 基礎からわかる時系列分析―Rで実践するカルマンフィルタ・MCMC・粒子フィルタ(第 6 章). 技術評論社. 2018.