【入門】データ同化・状態空間モデルをわかりやすく説明してみた

本記事では、データ同化の基本的な内容を紹介します。

自分も勉強中のため必ずしも正しい情報が含まれているわけではないことを留意してください…

むしろ、間違えている部分がありましたら本記事のコメント欄にご指摘よろしくお願いします！

また、わかりやすさを求めるために本記事を更新していきます。

データ同化とは
データ同化の定式化（一般状態空間モデル）
逐次ベイズフィルタ
参考資料
- 参考文献
まとめ : データ同化の基礎

データ同化とは

データ同化とは、『数値シミュレーションの考え方』と『データを利用する統計学の考え方』を融合する考え方のことです。

このように説明されることが多いですが、定義が結構揺れている気がします…

とりあえずは、『数値シミュレーション+データ解析 = データ同化』だと思っておけば良いと思います。

具体的な定式化方法はこれから説明します！

データ同化の定式化（一般状態空間モデル）

基本的にデータ同化という学問で扱うのは、以下のような時系列の確率モデル（一般状態空間モデル）です.

定義 : 一般状態空間モデル

一般状態空間モデルは、システムモデルと呼ばれる状態ベクトル$\mathbf{x}$の確率分布と観測モデルと呼ばれる観測ベクトル$\mathbf{y}$の確率分布で構成されます。

システムモデル

$$\mathbf{x}_{t} \sim p(\mathbf{x}_{t} \mid \mathbf{x}_{t-1})$$

観測モデル

$$\mathbf{y}_{t} \sim p(\mathbf{y}_{t} \mid \mathbf{x}_{t})$$

この後、数値シミュレーションとの関係を説明する際に、状態ベクトルや観測ベクトルのイメージが理解できると思うので、『状態ベクトル？観測ベクトル？』という方も安心してください。

少なくとも定義からわかることは以下の二つですね。

状態$\mathbf{x}_{t-1}$から次の状態$\mathbf{x}_{t}$が確率$p(\mathbf{x}_{t} \mid \mathbf{x}_{t-1})$で生成される
観測ベクトル$\mathbf{y}_{t}$は状態$\mathbf{x}_{t}$から確率$p(\mathbf{y}_{t}\mid \mathbf{x}_{t})$で生成される

つまり、一般状態空間モデルでは、以下のマルコフ性が成り立つことを仮定します！

\begin{align} p(\mathbf{x}_{t} \mid \mathbf{x}_{1:t-1}, \mathbf{y}_{1:t-1}) &= p(\mathbf{x}_{t} \mid \mathbf{x}_{t-1}) \\ p(\mathbf{y}_{t} \mid \mathbf{x}_{1:t}, \mathbf{y}_{1:t-1}) &= p(\mathbf{y}_{t} \mid \mathbf{x}_{t}) \end{align}

数値シミュレーションとシステムモデルの関係

仮にある状態$\mathbf{x}(t)$が以下の微分方程式に従うとします。

＊状態ベクトルは考える対象に依存します。例えば$\mathbf{x} = (\text{温度}, \text{湿気}, \text{高さ}, \ldots)$とかをイメージしてください。

$$\frac{d\mathbf{x}(t)}{dt} = \mathbf{g}(\mathbf{x}(t), t),~~\mathbf{x}(0) = \mathbf{x}_{0}$$

この方程式を数値的に解く場合以下のように差分化して方程式を解くことが多いです。

\begin{align} \mathbf{x}_{t} &= \mathbf{x}_{t-1} + \mathbf{g}(\mathbf{x}_{t-1}, t) \Delta t \\ &\equiv \mathbf{f}_{t}(\mathbf{x}_{t-1}) \end{align}

このような定式化をシミュレーションモデルと呼ぶことにします。

しかし、現実的には離散化による誤差や、そもそも微分方程式自体が不完全な可能性もあるため、シミュレーションモデルが厳密とは限りません…

それ以外にも初期値・境界条件等の不確実性も一般には含まれていたり…

そのような非厳密性・不確実性を考慮し、以下のようにシミュレーションモデルにノイズを加えたものを定義します。

$$\mathbf{x}_{t} = \mathbf{f}_{t}(\mathbf{x}_{t-1}) + \mathbf{v}_{t}$$

ここで、$\mathbf{v}_{t}$はある確率密度$q(\mathbf{v}_{t})$に従い、システムノイズと呼ばれます（より一般には、システムノイズも含む$\mathbf{f}_{t}(\mathbf{x}_{t-1}, \mathbf{v}_{t})$となりますが今回は簡単のために上記のように表せることを仮定します）

このシステムノイズを加えることで上述の非厳密性・不確実性をカバーするという作戦です！！

シミュレーションモデルを背後に仮定した場合、システムモデルは確率密度関数の変換式より以下のように表せます。

\begin{align} p(\mathbf{x}_{t} \mid \mathbf{x}_{t-1}) &= p(\mathbf{v}_{t} = \mathbf{f}_{t}^{-1}(\mathbf{x}_{t}, \mathbf{x}_{t-1}) \left\|\frac{\partial \mathbf{v}_{t}}{\partial \mathbf{x}_{t}} \right\| \\ &= q(\mathbf{x}_{t}-f_{t}(\mathbf{x}_{t-1})) \end{align}

一般化状態モデルは、ノイズありのシミュレーションモデルを含むより広いクラスだということもわかりますね！

数値シミュレーションと観測モデルの関係

システムモデルにより得られた状態$\mathbf{x}_{1}, \ldots, \mathbf{x}_{T}$がそのまま観測データとなるケースもありますが、$\mathbf{x}$は一般に超高次元なので観測されるのは$\mathbf{x}$の一部または何らかに変換された量になることが一般的です。

どのレベルのことを知りたいかに依存しますが、皆さんがホットミルクを作る場合は超高次元の水分子等の動きを調べるのではなく、温度や総量を観測して美味しいホットミルクを作りますよね…そんな感じです（イメージが悪いかもしれない）

そのような観測プロセスを以下のように定式化します。

$$\mathbf{y}_{t} = \mathbf{h}_{t}(\mathbf{x}_{t})$$

注意しなくてはならないのが、システムモデルによって決定される$\mathbf{x}$は確率変数なので、$\mathbf{h}_{t}(\mathbf{x})$も確率変数となることです。

後に詳しく説明しますが、確率的な形式で表しておくと、確率分布$p(\mathbf{h}(\mathbf{x}_{t})$と実際の観測値を比較し、状態の確率分布$p(\mathbf{x}_{t})$を修正することができ、嬉しいというわけです（確率的な定式化をしなくても同じようなことができるらしいですが不勉強でよくわかりません）

また同様に、モデルが再現できない部分の影響や測定誤差を考慮し、観測プロセスを以下のように変形します（システムモデル同様、より一般には$h_{t}(\mathbf{x}_{t}, \mathbf{w}_{t})$となりますが簡単のため以下の形式で話を進めます）

$$\mathbf{y}_{t} = \mathbf{h}_{t}(\mathbf{x}_{t}) + \mathbf{w}_{t}$$

ここで、$\mathbf{w}_{t}$は確率分布$r(\mathbf{w}_{t})$に従う確率変数で観測ノイズと呼ばれます。

このような確率的な観測プロセスを仮定した場合、観測モデルは以下のように表せます。

\begin{align}p(\mathbf{y} \mid \mathbf{x}_{t}) &= p(\mathbf{w}_{t} = \mathbf{h}_{t}^{-1}(\mathbf{y}_{t}, \mathbf{x}_{t}) \left\|\frac{\partial \mathbf{w}_{t}}{\partial \mathbf{y}_{t}} \right\| \\ &= r(\mathbf{y}_{t}-\mathbf{h}_{t}(\mathbf{x}_{t}) \end{align}

一般化状態空間モデルはノイズを考慮した観測プロセスを含むことがわかります！

システムモデル・観測モデルによる分類

一般化状態空間モデルは、かなり一般化されたモデルであることがわかりました。

そのため、そのクラスの中でもシステムモデルや観測モデルの形によって呼び方が異なります。

その一部を下記に簡単にまとめておきます。

線形・ガウス状態空間モデル（or 状態空間モデル）: 行列$F_{t}$, $G_{t}$, $H_{t}$に対して

\begin{align} \mathbf{x}_{t} &= F_{t} \mathbf{x}_{t-1} + G_{t} \mathbf{v}_{t},~~~\mathbf{v}_{t} \sim \mathcal{N}(\mathbf{0}, Q_{t}) \\ \mathbf{y}_{t} &= H_{t}\mathbf{x}_{t} + \mathbf{w}_{t},~~~\mathbf{w}_{t} \sim \mathcal{N}(\mathbf{0}, R_{t}) \end{align}

線形・非ガウス状態空間モデル（上述の例）

\begin{align} \mathbf{x}_{t} = \mathbf{f}_{t}(\mathbf{x}_{t-1}) + \mathbf{v}_{t},~~~\mathbf{v}_{t} \sim q \\ \mathbf{y}_{t} = \mathbf{h}_{t}(\mathbf{x}_{t})+ \mathbf{w}_{t},~~~\mathbf{w}_{t} \sim r \end{align}

非線形・非ガウス状態空間モデル

\begin{align} \mathbf{x}_{t} = \mathbf{f}_{t}(\mathbf{x}_{t-1}, \mathbf{v}_{t}),~~~\mathbf{v}_{t} \sim q \\ \mathbf{y}_{t} = \mathbf{h}_{t}(\mathbf{x}_{t}, \mathbf{w}_{t}),~~~\mathbf{w}_{t} \sim r \end{align}

他にも分類はありますが、今回は三つを覚えておけば十分です。

逐次ベイズフィルタ

『数値モデルとシミュレーションの関係』で少し述べましたが、データ同化の目的は、観測データ集合$\mathbf{y}_{1:t}$が与えられた時に状態ベクトル$\mathbf{x}_{1:t}$を修正し、推定制度を上げることでした。

つまり、目的は観測データ集合$\mathbf{y}_{1:T}$が与えられたときに、その時刻までの状態ベクトル$\mathbf{x}_{1:t}$の確率分布を推定することです。

このような推定を行うために便利な条件付き分布が三つ存在します。

（一期先）予測分布: $p(\mathbf{x}_{t} \mid \mathbf{y}_{1:t-1})$

フィルタ分布: $p(\mathbf{x}_{t} \mid \mathbf{y}_{1:t})$

（固定区間）平滑化分布: $p(\mathbf{x}_{t} \mid \mathbf{y}_{1:T})$

この三つ分布間には便利な漸化式が存在し、この三つの分布を知ることで任意の条件付き確率が厳密に計算できます。

この漸化式を詳しく説明していきます！

＊観測データを使用し、状態ベクトルを修正することで将来のデータ$\mathbf{x}_{t+1}, \mathbf{y}_{t+1}, \ldots) $の推定制度を上げることができます。

（一期先）予測分布の漸化式

（一期先）予測分布の意味は、$t-1$までのデータに基づいた状態ベクトル$\mathbf{x}_{t}$の推定を表しています。

そして一期先予測分布は一般化状態空間モデルを仮定した場合、$t-1$のフィルタ分布$p(\mathbf{x}_{t-1} \mid \mathbf{y}_{1:t-1}$を使って以下のように表せます。

$$p(\mathbf{x}_{t} \mid \mathbf{y}_{1:t-1}) = \int p(\mathbf{x}_{t} \mid \mathbf{x}_{t-1})p(\mathbf{x}_{t-1} \mid \mathbf{y}_{t-1})$$

一般化状態空間モデルのマルコフ性を利用すると導出できます。

（一期先）と括弧付きで限定していますが、一期先ではなく$t-1$までのデータに基づいて$\mathbf{x}_{t+1}, \mathbf{x}_{t+2}, \ldots$を長期先を予測するようなケース場合もあり、同様にフィルタ分布を使って計算できます！

フィルタ分布の漸化式

フィルタ分布の意味は、$t$までのデータに基づいた状態ベクトル$\mathbf{x}_{t}$の推定を表しています。

使用できる観測データ数が予測分布よりも多いため推定精度はフィルタ分布の方が良いです。

そしてフィルタ分布も一般化状態空間モデルを仮定した場合、（一期先）予測分布$p(\mathbf{x}_{t} \mid \mathbf{y}_{1:t-1})$を用いて以下のように表すことができます。

$$p(\mathbf{x}_{t} \mid \mathbf{y}_{1:t}) = \frac{p(\mathbf{y}_{t} \mid \mathbf{x}_{t}p(\mathbf{x}_{t} \mid \mathbf{y}_{1:t-1})}{\int p(\mathbf{y}_{t} \mid \mathbf{x}_{t})p(\mathbf{x}_{t} \mid \mathbf{y}_{1:t}}$$

分子は、（一期先）予測分布$\mathbf{x}_{t} \mid \mathbf{y}_{1:t-1})$と尤度関数$p(\mathbf{y}_{t} \mid \mathbf{x}_{t})$積で表されています。

これは、予測分布からデータによく適合する状態により大きい確率重みをかけたものがフィルタ分布と解釈可能です。

$p(\mathbf{y}_{t} \mid \mathbf{x}_{t})$を尤度関数と呼びましたが、他の尤度関数と区別するため、一般に一時点尤度関数と呼ばれます。

予測分布とフィルタ分布に関しては、紹介した漸化式を用いれば終時間$T$時の予測分布$p(\mathbf{x}_{T} \mid \mathbf{y}_{1:T-1})$とフィルタ分布$p(\mathbf{x}_{T} \mid \mathbf{y}_{1:T})$が求まります。

（固定区間）平滑化分布の漸化式

（固定区間）平滑化分布の意味は、$T$までのデータに基づいた状態ベクトル$\mathbf{x}_{t}$の推定を表しています。

使用できる観測データ数が最も多いため推定精度は予測分布・フィルタ分布より良いです。

平滑化分布も一般化状態空間モデルを仮定した場合、以下の漸化式が成り立ちます。

$$p(\mathbf{x}_{t} \mid \mathbf{y}_{1:T}) = p(\mathbf{x}_{t} \mid \mathbf{y}_{1:t}) \int \frac{p(\mathbf{x}_{t+1} \mid \mathbf{x}_{t})p(\mathbf{x}_{t+1}\mid \mathbf{y}_{1:T})}{p(\mathbf{x}_{t+1} \mid \mathbf{y}_{1:t})}d\mathbf{x}_{t+1}$$

この漸化式から、平滑化分布$p(\mathbf{x}_{t} \mid \mathbf{y}_{1:T})$を求めるためには、$p(\mathbf{x}_{T} \mid \mathbf{y}_{1:T-1})$, $p(\mathbf{x}_{T} \mid \mathbf{y}_{1:T})$を求め、$p(\mathbf{x}_{T-1} \mid \mathbf{y}_{1:T}), p(\mathbf{x}_{T-2} \mid \mathbf{y}_{1:T}) \ldots,$のように時間を遡り求めることになります。

今回も（固定区間）と限定句がついていることから、固定区間でない平滑化分布があることも予想できますね。

実際、$L$個先のデータを使用して$p(\mathbf{x} \mid \mathbf{y}_{1:t+L})$を利用して推定を行う『固定ラク平滑化分布』という方法もあります。

ここまでで、予測分布・フィルタ分布・平滑化分布を求める漸化式を説明しました。

このような漸化式を使えば区間$1 \sim T$の任意の条件付き確率を厳密に求めることができます。

このような方法は逐次ベイズフィルタと呼ばれています。

逐次ベイズフィルタの問題点

逐次ベイズフィルタの問題点は予測分布・フィルタ分布・平滑化分布に現れる状態ベクトル$\mathbf{x}$の積分です。

一般に状態ベクトルは超高次元になるため、特定のケースを除いて解析的に評価するのは難しく、数値計算に頼ることになります。

しかし、数値計算に頼ったとしても超高次元の数値積分はなかなか難しくさまざまな問題が生じます。

そのため、各システムモデル・観測モデルに応じて使用する手法が様々で各モデルに対して最適な手法を使用する必要があります！

具体的に逐次ベイズフィルタを実行する方法をいくつか紹介します

カルマンフィルタ : 線形・ガウス状態空間モデルに対して解析的に予測分布・フィルタ分布・平滑化分布を評価する
粒子フィルタ（パーティクルフィルタ）: 非線形・非ガウス状態空間モデルに対して数値的に予測分布・フィルタ分布・平滑化分布を評価する
アンサンブルカルマンフィルタ : システムモデルが非線形・観測モデルが線形の状態空間モデルに対して解析的手法と数値的手法を組み合わせて予測分布・フィルタ分布・平滑化分布を評価する

パラメータの最適化

システムノイズ・観測ノイズの従う確率分布のパラメータ、またはシミュレーションモデル内のパラメータも観測データから推定することができます。

具体的には、得られた観測データ$\mathbf{y}_{1:T}$に対する対数尤度を最大化するようにあるパラメータ$\theta$を決定すれば良いです。

\begin{align} l(\theta) &\equiv \log p(\mathbf{y}_{1:T} \mid \theta) \\ &= \sum_{t=1}^{T}\log p(\mathbf{y}_{t} \mid \mathbf{y}_{1:t-1}, \theta) \end{align}

この形式見ればわかるように$p(\mathbf{y}_{t} \mid \mathbf{y}_{1:t-1}, \theta)$は以下のようにフィルタ分布の分母の計算で現れるため追加で計算は必要ありません！

\begin{align}&\int p(\mathbf{y}_{t} \mid \mathbf{x}_{t})p(\mathbf{x}_{t} \mid \mathbf{y}_{1:t-1})d\mathbf{x}_{t} \\ = &\int p(\mathbf{y}_{t}, \mathbf{x}_{t} \mid \mathbf{y}_{1:t-1}) \\ = &p(\mathbf{y}_{t} \mid \mathbf{y}_{1:t-1}) \end{align}

しかし、対数尤度を計算することはできますが、最大化するための最適化は新たに直面する問題です。

他にも拡大された状態ベクトルを再定義することでパラメータを推定する方法もありますが今回は省略します。