時系列解析初心者が始める時系列解析(8)

2023/01/16に公開約7,300字

局所定常ARモデル

参考:https://elf-c.he.u-tokyo.ac.jp/courses/386


概要


  • いよいよ、非定常時系列モデルに移っていく。今までは定常時系列モデルであったが、世の中のデータが何でもかんでも定常性を持ったデータなことはないので、非定常性を持ったデータにも対応できるモデルを構築する必要がある。

  • そこで今回はその第1歩として、区分的定常モデルというのを学んでいく。次回は状態空間モデルとなる。

  • 区分的定常モデルは、今までに学んだ手法を応用して、非定常時系列モデルを構築していくというようなイメージ。


レベルシフトの検出


レベルシフトとは?

  • ある時系列データがあった時に、それがどっかのタイミングで構造が変化したと捉えて、その点の前と後で分けて、モデリングしていくということ。(式が変わるということではなく、パラメータを増やすということ。)

  • 例えば、シートベルトの義務化の法案が施行される前と後で、事故率の時系列データに変化点があったのかという問題があったとする。その時に、定常時系列モデルとしてモデリングした時の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

\begin{align}&l(k, N_j, m_j, (a_{ji}, i=1,\dots,m_j), \hat\sigma_j^2;j=1,\dots,k)\\ &= -\frac{1}{2}\sum_{j=1}^k\{N_j\log 2\pi \hat{\sigma}_j^2 + N_j\}\\ &= -\frac{N-m}{2}(\log2\pi+1) - \frac{1}{2}\sum_{j=1}^k{N_j \log{\hat\sigma_j^2}} \end{align}
  • 比較的簡単な式になったと思われる。この式の要素を用いて、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つの手法で実装していく。

    1. 任意個の区間への自動分割・・・大まかに分割
    2. 変化点の精密推定・・・変化点っぽい所にあたりをつけて、1カ所を精密に推定

任意の区間への自動分割


スイッチするか、プールするか

  • 定常性を仮定する区間の長さの最小値をLとして、ARの最大次数Mを決める。

  • この時、長さLの新しいデータが得られる度に、下記2つのモデルで算出したAICの値を比較し、小さい方をCurrent~modelとして選択する。

    • Switched~model:ある長さLのデータでモデリングした際のAICをAIC_1(Old model)とし、新たに得られた長さLのデータでモデリングした際のAICをAIC_2(Current model)として、全体のモデルのAICをAIC = AIC_1 + AIC_2とする。これは、前者と後者のデータの構造が異なるという仮定がある。

    • Pooled~model:前者と後者に構造の違いが見受けられないのであれば、一色たんにしてモデリングしても問題ないはずなので、前者、後者データまとめてモデリングする。実際に計算するとすれば、前者、後者のデータ(Householder変換済)をまとめて、Householder変換して、それを用いて推定を行う。

  • この時、先ほどの節で説明した数式を用いて計算するわけだが、その中に含まれるN_jは次のように、N_j = c_jL(D_j区間の長さ)と表される。(c_jは整数)

  • この計算手法では、Lを基準にサンプリングして、SwitchedPooledかで、どちらが当てはまりがいいかをその都度考え、それをデータ全体まで施す。(動画の具体例を見てみるとわかりやすい。)


変化時点の精密な推定


前半のモデルと後半のモデルの足し合わせ

  • 考え方はシンプルで、変化点を含むであろう空間を捉えて、その中の時点を軸に、前半モデルと後半モデルを構築するというもの。

  • 最初からドンピシャの変化点を捉えられるわけはないので、1時点ずつずらしながら、前半と後半のモデリングを実行して、それぞれAICを算出し、比較して、最小値となったところを変化点と捉える。AIC_j = AIC_j^1 + AIC_j^2

  • 計算式は、範囲に注意は必要なものの、通常のARモデルと同様。

  • ここで課題になるのが、計算スピードを速くする(=計算量を抑える)ことを備えること。

  • 結局はHouseholder法を用いるのだが、少し工夫する。それは、データの逐次追加である。単純にARモデルを1から構築するよりもベースに対して、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の中の1つと考える。
      p(j|y) = \frac{p(y|j)p(j)}{\sum_{j=1}^L{p(y|j)p(j)}}
  • これを計算すると、波形にどんな山(衝撃)が発生したかを可視化できる。すると、今節の冒頭で述べたような、誤った観測を見破ることができる。


統計的制御


統計的制御とは?

  • この話は今までと大きく変わります。

  • 統計的制御と検索すると”統計的プロセス制御”というのが出てくるが、これは別物。

  • 統計的制御の定義は不明だが、利用例として挙げられるのは、船舶の自動操舵がある。波の波形の変化を適切に捉えることが出来れば、どういった方向に進むかを制御できる。

  • こういった現象には理論があるはずで、理論モデルを作成したいとこだが、「複雑」であったり、「外乱が強い」とそれを適切に理論に落とし込むことが難しい。ただ、統計的には、その内容を踏まえてモデリングするので、制御が可能となる。

数学的な説明について

  • 数学的には、次章の状態空間モデルが根幹にあるため、ここで説明することが叶わない。(現状、何を言っているかの整理がつかない)

  • 動画では、数式が出ているので、興味ある方は動画をご参照ください。


船舶のオートパイロット・NADCON


  • 上記を実際に実装した事例を紹介されている。事例なので、動画をご参照ください。

  • 実装した内容はシステム化されており、実際の船に搭載されているらしい。

Discussion

ログインするとコメントできます