時系列解析初心者が始める時系列解析(8)
局所定常ARモデル
参考：https://elf-c.he.u-tokyo.ac.jp/courses/386
概要
-
いよいよ、非定常時系列モデルに移っていく。今までは定常時系列モデルであったが、世の中のデータが何でもかんでも定常性を持ったデータなことはないので、非定常性を持ったデータにも対応できるモデルを構築する必要がある。
-
そこで今回はその第１歩として、区分的定常モデルというのを学んでいく。次回は状態空間モデルとなる。
-
区分的定常モデルは、今までに学んだ手法を応用して、非定常時系列モデルを構築していくというようなイメージ。
レベルシフトの検出
レベルシフトとは？
-
ある時系列データがあった時に、それがどっかのタイミングで構造が変化したと捉えて、その点の前と後で分けて、モデリングしていくということ。（式が変わるということではなく、パラメータを増やすということ。）
-
例えば、シートベルトの義務化の法案が施行される前と後で、事故率の時系列データに変化点があったのかという問題があったとする。その時に、定常時系列モデルとしてモデリングした時のAICと、レベルシフトを施してモデリングした時のAICを比較して、後者の方がAICが優れている（値が小さい）のであれば、そのデータには変化点が存在したということになる。今回の例で言えば、法案の施行には効果があったと捉えることができるのである。
数学的に説明してみる
-
ここからはもう少し数学的な説明に落とし込んでみる。まず、レベルシフトを仮定した時の時系列は、
y_n \sim N(\mu_n,\sigma^2) \mu_n -
この場合のモデリングとしては、下記のようになる。これは、正規分布の確率密度関数の公式に平均、分散を代入した結果である。
p(y_n|\mu_n,\sigma^2) = (2\pi\sigma^2)^{-\frac{1}{2}} \exp\bigg\{-\frac{(y_n - \mu_n)^2}{2\sigma^2}\bigg\} -
定常時系列モデルであれば、この式を用いて尤度を算出するが、今回はレベルシフトということなので、もう1つ考えることがある。それは、変化点として時刻
k+1 k k+1 k
\mu_n = \begin{cases} \theta_1 ~~~~~ n≤k \\ \theta_2 ~~~~~ n>k \end{cases} -
以上の条件を用いて、尤度の式を構築すると下記のようになる。
\Pi k
L(\theta_1,\theta_2,\sigma^2_k) = \Pi_{n=1}^{k}{p(y_n|\theta_1,\sigma^2_k)} \Pi_{n=k+1}^{N}{p(y_n|\theta_2,\sigma^2_k)} -
この式の対数をとり、対数尤度の式を構築する。
l(\theta_1,\theta_2,\sigma^2_k) = -\frac{N}{2}\log{(2\pi\sigma_k^2) - \frac{1}{2\sigma^2_k}}\bigg\{ \sum_{n=1}^k(y_n - \theta_1)^2 + \sum_{n=k+1}^N{(y_n-\theta_2)^2} \bigg\} -
上記式をパラメータそれぞれで偏微分し、
=0
\begin{align} \hat\theta_1 &= \frac{1}{k}\sum_{n=1}^k{y_n} \\ \hat\theta_2 &= \frac{1}{N-k}\sum_{n=k+1}^N{y_n} \\ \hat\sigma_k^2 &= \frac{1}{N}\bigg\{ \sum_{n=1}^k(y_n-\hat\theta_1)^2 + \sum_{n=k+1}^N(y_n -\hat\theta_2)^2 \bigg\} \end{align} -
式(1)(2)(3)を代入することで、最大対数尤度を推定するための式が構築できる。
l(\hat\theta_1,\hat\theta_2,\hat\sigma_k^2) = -\frac{N}{2}\log(2\pi\hat\sigma_k^2) -\frac{N}{2} -
この式におけるパラメータの数は、3つとなるため、AICの計算式は下記のようになる。
AIC_k = N\log(2\pi\hat\sigma^2_k) + N + 2 \times 3 -
ここまでできたら、様々なモデルに対してのAICを算出してみる。その時のAICの動きが変化点の部分で最小値を取るようなことがあれば、変化点を軸にレベルシフトを行なった方がAICは低くなるような結果が得られる。
局所定常ARモデル
局所定常って？
- 文字通り、局所（一部）的には定常性が捉えられるデータのこと。地震波のように、初期微動と主要動を合わせてみると、平均も分散も異なるが、それを初期微動と主要動とそれぞれ分けてみることで、それぞれに定常性が見られるということ。
局所定常ARモデルを数学的に
-
先ほどの地震波のような初期微動や主要動を区間とすると、区間
D_j D_1 D_2 -
区間
D_j -
D_j AR(m_j) m_j ≤ M
y_n = \sum_{i=1}^{m_j}a_{ij}y_{n-i} + v_{nj} \sim N(0,\sigma_j^2), ~~~~~ n \in D_j -
このような考え方でモデリングをすることで、以下のようなメリットが得られる。
- 平均値の変化に対応できる。
- 分散の変化に対応できる。
- スペクトル・共分散構造（波形）の変化に対応できる。
- 任意の時点でのデータ構造の変化に対応できる。
- AICによる自動推定が可能になる。（？）
-
ARモデルということで、結局は尤度を求めるようにしていく必要があるわけだが、その際はどのようにデータを扱っていくのかということになる。
-
まず、時間そのものの分割である。全体時間を
N k k D
N = N_1 + N_2 + \cdots + N_k -
次に、一つの区間に焦点を当てた時の区間
D
\begin{align} D_i &= [n_{i0}, n_{i1}]\\ n_{i0} &= \sum_{j=1}^{i-1}{N_j} + 1\\ n_{i1} &= \sum_{j=1}^i{N_j} \end{align} -
この時、実際に推定する必要があるパラメータは、通常のARモデルに比べて増える格好となる。
\theta = (k, N_j, m_j, (a_{ji}, i=1,\dots,m), \sigma_j^2;j=1,\dots,k)^T -
尤度の式は以下の通り。
\begin{align} L(\theta) &= p(y_1,\dots,y_N|\theta)\\ &= \Pi_{j=1}^k\Pi_{n=n_{j0}}^{n_{j1}}p(y_n|y_1,\dots,y_{n-1},\theta)\\ &\cong \Pi_{j=1}^k(2\pi\sigma_j^2)^{-\frac{N_j}{2}}\exp\bigg\{-\frac{1}{2\sigma_j^2} \sum_{n=n_{j0}}^{n_{j1}}(y_n - \sum_{i=1}^{m_j}{a_{ji}y_{n-i}})^2 \bigg\} \end{align} -
さらに、これの対数を取ると、、、
l(k, N_j, m_j, (a_{ji}, i=1,\dots,m_j), \sigma_j^2;j=1,\dots,k)\\ = -\frac{1}{2}\sum_{j=1}^k\bigg\{N_j\log{2\pi\sigma_j^2} + \frac{1}{\sigma_j^2} \sum_{n=n_{j0}}^{n_{j1}}\bigg(y_n - \sum_{i=1}^{m_j}{a_{ji}y_{n-i}}\bigg)^2 \bigg\} -
ここで、
j
\hat\sigma_j^2= \frac{1}{N_j} \sum_{n=n_{j0}}^{n_{j1}}\bigg(y_n - \sum_{i=1}^{m_j}{a_{ji}y_{n-i}}\bigg)^2
- 比較的簡単な式になったと思われる。この式の要素を用いて、AICの算出式を構築する。
AIC = (N-m)(\log2\pi+1) + \sum_{j=1}^k{N_j\log\hat\sigma_j^2} + 2\sum_{j=1}^k(m_j + 1)
現実的な話
-
これまでに記述した式で求めて、AICまで算出するというのは、非現実的。
-
理由として、全ての区分数(
k N_1,\cdots,N_k -
というわけで、主に2つの手法で実装していく。
- 任意個の区間への自動分割・・・大まかに分割
- 変化点の精密推定・・・変化点っぽい所にあたりをつけて、１カ所を精密に推定
任意の区間への自動分割
スイッチするか、プールするか
-
定常性を仮定する区間の長さの最小値を
L M -
この時、長さ
L Current~model -
Switched~model L AIC_1 L AIC_2 AIC = AIC_1 + AIC_2 -
Pooled~model
-
-
この時、先ほどの節で説明した数式を用いて計算するわけだが、その中に含まれる
N_j N_j = c_jL D_j c_j -
この計算手法では、
L Switched Pooled
変化時点の精密な推定
前半のモデルと後半のモデルの足し合わせ
-
考え方はシンプルで、変化点を含むであろう空間を捉えて、その中の時点を軸に、前半モデルと後半モデルを構築するというもの。
-
最初からドンピシャの変化点を捉えられるわけはないので、１時点ずつずらしながら、前半と後半のモデリングを実行して、それぞれAICを算出し、比較して、最小値となったところを変化点と捉える。
AIC_j = AIC_j^1 + AIC_j^2 -
計算式は、範囲に注意は必要なものの、通常のARモデルと同様。
-
ここで課題になるのが、計算スピードを速くする(=計算量を抑える)ことを備えること。
-
結局はHouseholder法を用いるのだが、少し工夫する。それは、データの逐次追加である。単純にARモデルを１から構築するよりもベースに対して、1つデータを追加するというふうにした方が抑えられるという理論。
-
詳細については動画をご参照ください。
地震データへの適用例
- 主に、事例のお話なので、動画をご参照ください。
- 内容としては、北海道の浦河というところの観測データを東、南、上下方向で波形を見て、局所定常ARモデルを当てはめて、どのタイミングで地震が発生しているのかを観測しています。
- これは、緊急地震速報に応用されており、東日本大震災では、新幹線の緊急停止を、地震がくる1分前には実施することができ、脱線を防ぐことができたとのこと。
変化点の事後確率
-
前節の地震観測の応用事例には、いくつか欠点があった。
- 震源推定の誤差：少数の地震波に基づいて推定しいるため、一部の到着時刻推定の誤差（外れ値）が大きな影響を及ぼしてしまう。
- 緊急地震速報の誤動作：引き続いて2件の小さな地震が発生した時、それを同一の地震と見做してしまう。（波形が繋がっているように見えるため。）
-
こういった課題に対処するための方法として、変化点の事後確率を推定の要素として、盛り込むというものがある。
-
では、どうやって盛り込んでいくのか。まず考えることとしては、AICの性質である。AICは、元を辿ると尤度であり、それをバイアス補正して、対数を取ったものを考えられる。
AIC_j = -2 \times (バイアス補正した対数尤度) -
この考え方を用いて、
j y
p(y|j) \equiv \exp\bigg\{-\frac{1}{2}AIC_j\bigg\} -
また、この時の
p(j) -
p(j) = \frac{1}{L} L
p(j|y) = \frac{p(y|j)p(j)}{\sum_{j=1}^L{p(y|j)p(j)}}
-
-
これを計算すると、波形にどんな山（衝撃）が発生したかを可視化できる。すると、今節の冒頭で述べたような、誤った観測を見破ることができる。
統計的制御
統計的制御とは？
-
この話は今までと大きく変わります。
-
統計的制御と検索すると”統計的プロセス制御”というのが出てくるが、これは別物。
-
統計的制御の定義は不明だが、利用例として挙げられるのは、船舶の自動操舵がある。波の波形の変化を適切に捉えることが出来れば、どういった方向に進むかを制御できる。
-
こういった現象には理論があるはずで、理論モデルを作成したいとこだが、「複雑」であったり、「外乱が強い」とそれを適切に理論に落とし込むことが難しい。ただ、統計的には、その内容を踏まえてモデリングするので、制御が可能となる。
数学的な説明について
-
数学的には、次章の状態空間モデルが根幹にあるため、ここで説明することが叶わない。（現状、何を言っているかの整理がつかない）
-
動画では、数式が出ているので、興味ある方は動画をご参照ください。
船舶のオートパイロット・NADCON
-
上記を実際に実装した事例を紹介されている。事例なので、動画をご参照ください。
-
実装した内容はシステム化されており、実際の船に搭載されているらしい。
Discussion