「Pythonではじめるベイズ機械学習入門」を読む　～第4章潜在変数モデル④(状態空間モデル)～

はじめに

Pythonによるベイズモデルの実装をきちんと学ぼうと思い，森賀新・木田悠歩・須山敦志著「Pythonではじめるベイズ機械学習入門」を読むことにした。

本記事は，第4章「潜在変数モデル」のうち状態空間モデルに関する読書メモである。

本書の紹介ページ

www.kspub.co.jp

4.3 状態空間モデル

状態空間モデル(State Space Model)は，時系列データのモデリングに用いられる潜在変数モデルである。
前の時点 $t$ から次の時点 $t+1$ の状態を得るための方程式を状態方程式(State Equation)と呼ぶ。
また時点 $t$ における状態から観測値を得るための方程式を観測方程式(Observation Equation)と呼ぶ。

状態空間モデルを用いた時系列分析の利点としては，表現力・説明力の高さが挙げられる。すなわち，トレンド，周期変動などの複数成分の状態成分を組合わせたモデルを作ることが可能で，観測データに影響を与えるさまざまな要因を分解して考えることができる。

4.3.1 モデル概要 : ローカルレベルモデル

本節では，状態空間モデルのうち最も単純なモデルともいえるローカルレベルモデル(Local Level Model)について説明している。

ローカルレベルモデルにおける，状態方程式，観測方程式，状態の初期値は以下の式で表される。

$\begin{align} \mu_t &= \mu_{t-1} + \eta_t, \quad \eta_t \sim \mathcal{N}(0, \sigma_{\eta}^2) \\ \\ y_t &= \mu_t + \varepsilon_t, \quad \varepsilon_t \sim \mathcal{N}(0, \sigma_{\varepsilon}^2) \\ \\ \mu_0 &\sim \mathcal{N}(m_0, \sigma_0 ^2) \\ \end{align}$

ここで $\mu_t$ は水準成分(レベル成分)と呼ばれるものである。
状態方程式は，「水準成分は，ランダムな値を足し合わせている」，すなわちランダムウォークを表している。そのままでは予測の役には立たないが，より実用的なモデルのベースになるモデルである。

4.3.2 実装：ローカルレベルモデル

サンプルコードを動かしながら，挙動を確認した。
ｖ
github.com

データの準備

データは，ハワイ島のマウナロア観測所における1966年から2019年の月次の二酸化炭素の測定データを利用する。なおこのデータは，サンプルコードにベタ打ちされている。

このうち，最後の10年分をテストデータに用いる。

事後分布の変分推論

実装には，TFP(TensorFlow Probability)のstsモジュールを利用する。これには以下のような特徴がある。

状態空間モデルに特化している。
変分推論法を用いた高速な推論が可能である。
初期値やシステムノイズ，観測誤差の事前分布を指定しなかった場合，観測値に基づき自動的に設定される。

事後分布を求めるには変分推論法を用いる。ELBOの負値のグラフは下図のようになる。

予測分布

将来予測の結果は，下図のようになる。

青線が実際の観測値，オレンジの破線が予測分布の平均，オレンジの帯が平均±標準偏差の区間を表す。
予測分布の平均は，一定の値になっている。これは，システムノイズ・観測誤差である $\eta, \varepsilon$ の平均が0であるためである。
このようにローカルレベルモデルでは，増加傾向や周期変動などを考慮していないので，観測データの特徴を捉えられていない。

4.3.3 モデル概要 : 構造時系列モデル

本節では，ローカルレベルモデルよりも実用的なモデルとして，構造時系列モデル(Structural Time Series Model)を紹介している。
構造時系列モデルには，以下のような成分が含まれる。

トレンド成分(Trend Component) : 増加傾向・減少傾向といったデータの大まかな変動傾向を表す。
季節成分(Seasonal Component) : 周期的な変動を表す。
不規則成分(Irregular Component) : 他の成分では表現できない変動を表す。

4.3.4 モデル概要 : トレンド成分

本節では，ローカル線形トレンドモデル(Local Linear Trend Model)を紹介している。

ローカルトレンドモデルでは，水準成分 $\mu_t$ に加えて，傾き成分 $\delta_t$ が現れる。
水準成分 $\mu_t$ は，前の時点の水準成分 $\mu_{t-1}$ に $\delta_t-1$ が加わっている。すなわち $\delta$ は $\mu$ の傾きを表している。

式で表すと以下のようになる。

$\begin{align} \mu_t &= \mu_{t-1} + \delta_{t-1}+ \eta_{\mu,t}, \quad \eta_{\mu,t} \sim \mathcal{N}(0, \sigma_{\eta_\mu}^2) \\ \\ \delta_t &= \delta_{t-1} + \eta_{\delta, t}, \quad \eta_{\delta, t} \sim \mathcal{N}(0, \sigma_{\eta_\mu}^2) \\ \\ y_t &= \mu_t + \varepsilon_t, \quad \varepsilon_t \sim \mathcal{N}(0, \sigma_{\varepsilon}^2) \\ \end{align}$

4.3.5 モデル概要 : 季節成分

季節成分は，周期的な変動を表現する状態成分である。
季節成分 $\gamma$ は， $s$ を周期としたときに，1周期分の合計は0になるという制約を設ける。このような制約を設けないと，季節成分が一位に定まらないためである。

$\begin{align} \sum_{i=1}^{s} \gamma_i = 0 \\ \end{align}$

この式より，時刻 $t$ における季節成分 $\gamma_t$ について，以下が成り立つことが分かる。

$\begin{align} \gamma_t &+ \gamma_{t-1} + \cdots \gamma_{t-s+1} = 0 \\ \\ \Rightarrow \gamma_t &= - \sum_{t-s+1}^{t-1} \gamma_i \\ \end{align}$

さらに季節成分には，季節成分の揺らぎをあらわすノイズ $\eta_{\gamma, t}$ を加える。

$\begin{align} \gamma_t &= - \sum_{t-s+1}^{t-1} \gamma_i + \eta_{\gamma, t} \\ \end{align}$

以上より，トレンド成分と季節成分を組み込んだモデル式は以下のようになる。

$\begin{align} \mu_t &= \mu_{t-1} + \delta_{t-1} + \eta_{\mu,t}, \quad \eta_{\mu,t} \sim \mathcal{N}(0, \sigma_{\eta_\mu}^2) \\ \\ \delta_t &= \delta_{t-1} + \eta_{\delta, t}, \quad \eta_{\delta, t} \sim \mathcal{N}(0, \sigma_{\eta_\mu}^2) \\ \\ \gamma_t &= - \sum_{t-s+1}^{t-1} \gamma_i + \eta_{\gamma, t}, \quad \eta_{\gamma, t} \sim \mathcal{N}(0, \sigma_{\eta_{\gamma}}^2 ) \\ y_t &= \mu_t + \gamma_t + \varepsilon_t, \quad \varepsilon_t \sim \mathcal{N}(0, \sigma_{\varepsilon}^2) \\ \end{align}$

4.3.6 実装：構造時系列モデル

構造時系列モデルについて，変分推論法を用いてパラメータの推定を行なう。
ELBOの負値の推移は下図のようになる。

また，ローカル線形トレンド＋季節成分モデルによる予測結果は以下の通りである。

ローカルレベルモデルに比べて，時系列データの特徴を捉えられている。

更に得られたトレンド成分と季節成分は下図のようになる。

まとめと感想

潜在変数モデルのうち，状態空間モデルの実装方法について学んだ。

今回の事例では，状態方程式・観測方程式ともに線形の式であり，ノイズも正規分布だったので，解析的に解けるパターンになると考えられる。
本書は事後分布をMCMCや変分推論で解く方法を紹介しているので変分推論で解いていたが，状態空間モデルにおいてこれらの近似推論が力を発揮するのは，式が非線形の場合であったり，ノイズが正規分布に従わないパターンであると考えられる。
今後はこのような問題設定における問題を解いてみたい。

本記事を最後まで読んでくださり，どうもありがとうございました。