Zenn
⏳

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

2022/12/26に公開
#
時系列解析
tech

モデルの推定・選択

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

AICによるモデル選択例

（１）ヒストグラムのBin数の選択

  • ヒストグラムを作成する際に、適切なビン数をAICを用いて判断することが可能。
  • 考え方としては、ヒストグラムを多項分布として捉えるということ。データ数は一定であるものの、ビン数を変更することで、分布が変化するということになる。その際の分布を数式で表し、尤度を求めて、AICを算出するというイメージ。（数式については、動画を参照。）

（２）分布の形状の選択

  • 場面の例として、あるデータ点が20個あった時に、それがどんな分布に従っているかを考えるということがあったとする。
  • そういった時には、ピアソン分布族という 形状パラメータb を変更することで分布の形状が変わる分布を用いて、当てはまりのいい分布を見つける。
  • その際の当てはまりの指標として、AICを用いるというロジック。具体的には、bを定めた時の最尤推定値を求め、その時のAICの値を捉えるということ。

（３）Box-Cox変換の選択

  • 元のデータyに対して、Box-Cox変換を実施し、正規分布に従うようなモデリングをした時にそれを対象にしたAICで比較することは不可能。
  • 理由として、\lambdaの値が0どうかで変換の仕組み（数式）が変わってしまうため。
  • これから言えることは、AICを用いてモデル選択を行うには、土俵を合わせる必要があるということ。
  • 今回でいえば、Box-Cox適応後のモデルに対して、逆変換し、元の形に戻すということで、AICを用いた比較を行うことができる。
  • 逆変換する際に注意すべきことは、Box-Cox変換の性質上、逆変換時には、ヤコビアン
    の項を設ける必要があるということ。
    AIC'_z = AIC_z - 2\log{\biggm\vert{\frac{dh_\lambda}{dy}\biggm\vert}}

（４）分布の同一性

  • 複数のデータセットがあったときに、同じ分布に従うとしたときに、どれくらい似ているのかを測る。
  • その時の指標としては、分布は同じなので、平均と分散に着目するということである。
  • 例として、あるデータセットAとあるデータセットBがあり、それらのモデルをf(A|\mu_A,\sigma^2_A)とf(B|\mu_B,\sigma^2_B)とする。（分布は正規分布）
  • この時、計4つのパターンでモデリングを実施して、AICを算出してみる。
    1. 制約なし
    2. \sigma^2_A = \sigma^2_B
    3. \mu_A = \mu_B
    4. \sigma^2_A = \sigma^2_B,~~\mu_A = \mu_B
  • すると、どの仮定が今回のデータセットでは、優れたモデルであるかを判断できる。

重回帰モデル

回帰モデルにおけるAICの登場場面

  • そもそも回帰モデルを作成する時には、それぞれフェーズがあり、その時々で評価指標が異なる。
順序 フェーズ名 評価指標 該当モデル
1 モデル選択 AIC 単回帰モデル、非線形回帰モデル
2 変数選択 AIC 重回帰モデル、自己回帰モデル、非線形回帰モデル
3 次数選択 AIC 単回帰モデル、自己回帰モデル、非線形回帰モデル
4 パラメータ推定 最尤法、最小二乗法 すべてのモデルが該当

重回帰モデルとは

  • 目的変数を1つ定め、説明変数を任意の次数分用意して、残差を付け加えたもの。数式の方がわかりやすいの数式で表す。
    y_n = a_1x_{n1} + \cdots + a_mx_{nm} + \varepsilon_n,~~\varepsilon \backsim N(0, \sigma^2)
名称 数式
目的変数 y_n
説明変数 x_{n1},\cdots,x_{nm}
次数 m
回帰係数 a_i
残差 \varepsilon_n
残差分散 \sigma^2

  • この中で、推定の対象となるパラメータ\thetaは”回帰係数”と”残差分散”である。

    \theta = (a_1, \cdots, a_m,\sigma^2)^T

  • これを用いて、重回帰モデルを確率分布で表現すると、以下のようになる。

    p(y_n|\theta,x_{n1},\cdots,x_{nm})\backsim N(\sum_{j=1}^m{a_jx_{nj},\sigma^2})

  • これを対象に尤度→対数尤度を求めていく。

    L(\theta)=\Pi_{n=1}^N{p(y_n|\theta,x_{n1},\cdots,x_{nm})}~~尤度

    l(\theta) = \log{L(\theta)}=\sum_{n=1}^N{\log{p(y_n|\theta,x_{n1},\cdots,x_{nm})}}~~対数尤度

  • この対数尤度を用いて、最尤推定していくわけだが、今回（線形正規型回帰モデル）の場合は最小二乗法を解くことがそれにあたる。

  • モデルは項それぞれベクトル、行列で表すことができる。

    y=\begin{bmatrix}y_1\\ y_2\\ \vdots\\ y_N \end{bmatrix}~~ 目的変数ベクトル

    Z=\begin{bmatrix}x_{11}&\cdots&x_{1m}\\ x_{21}&\cdots&x_{2m}\\ \vdots&\ddots&\vdots\\ x_{N1}&\cdots&x_{Nm} \end{bmatrix}~~ 説明変数行列

    \varepsilon=\begin{bmatrix}\varepsilon_1\\ \varepsilon_2\\ \vdots\\ \varepsilon_N\end{bmatrix}~~残差ベクトル

    a=\begin{bmatrix}a_1\\ \vdots\\ a_m\end{bmatrix}~~回帰係数ベクトル

y=Za+\varepsilon
  • 実際にAICを求めると、、、
    • k：パラメータ数。+1というのは、残差分散\sigma^2のこと。
      \begin{align}AIC_m &= -2l_m(\hat{\theta})+2k\\ l_m(\hat\theta)&=-\frac{N}{2}\log{2\pi\hat\sigma^2_m} - \frac{N}{2}\\ k&=m+1 \end{align}

      AIC_m = N\log{2\pi\hat\sigma^2_m}+N+2(m+1)

多項式回帰モデル

多項式回帰モデルとは？

  • 一般には、y_i = f(x_i) + \varepsilon_iと表すことができ、f(x),y_iについては、下記のように表すことができる。（ちなみに、\varepsilonは何かしらの分布に従うが、今回はN(0,\sigma^2)とする）

    f(x) = a_0 + a_1x + \cdots + a_mx^m

    y_i \backsim N(\mu_i, \sigma^2), \mu_i=a_0 + a_1x_i + \cdots + a_mx_i^m

  • 多項式回帰モデルについては、\muや\sigma^2、y_iが従う分布、mで表される次数によって、データ点それぞれの時のグラフの移動が変わる。

    • 実際のグラフは3次元で表され、横軸がx、縦軸がy、高さがf(y|x)となる。

  • 多項式回帰でのAICの計算は下記のようになる。

    • m+2はパラメータ数である。重回帰モデルの場合、切片を考慮していないため、m+1となっていたが、多項式回帰では切片があるため、パラメータ数が不一致。ちなみに、重回帰モデルで切片を考慮しない理由関しては、データを標準化することで、無視できるものになるかららしい。
      l_m(\hat\theta)=-\frac{N}{2}\log{2\pi\hat\sigma^2_m}-\frac{N}{2} ~~対数尤度

      AIC_m=N\log2\pi\hat\sigma^2_m + N + 2(m+2)

回帰モデルの変数選択と予測誤差

多変量データでは変数選択が必要になる

  • 多変量の場合、用いる説明変数も任意に変更できる。
  • 例として、目的変数を1月の「気温」、説明変数をその時の「緯度」「経度」「標高」とする。データ点としては、都道府県。
  • 考え方としては、それぞれの説明変数を組み合わせていき、その時のAICを算出し、比較する。今回で言えば、説明変数は計３つあり、1つ選ぶ時、2つ選ぶ時、すべて選ぶ時と場合分けをして、算出し、比較していく。
    1. 定数項だけではどうなるのか？を考え、AICを算出
    2. 定数項＋「緯度」とした時のAICを算出
    3. 定数項＋「経度」とした時のAICを算出
    4. continue...

予測誤差

  • モデル構築に用いたデータセット以外で、データを用意し、モデルで回帰予測をし、予測値を算出。
  • 実測値と予測値の残差を計算し、ばらつき(分散)を計算して、評価する。

最小二乗法、Householder法

最小二乗法の解法

  • 回帰モデルy = Za + \varepsilonにおいて、残差が最小となる値をもとめる。\|\varepsilon\|^2 = \|y-Za\|^2

  • この時の回帰係数を推定する従来の方法は、以下の式である。

    \hat{a} = (Z^TZ)^{-1}Z^Ty

  • ただ、近年は直交変換に基づく方法が流行り。

    • 直交変換：直交行列を用いた線形変換のこと。性質として、v,wというベクトルに対して、直交行列Aを用いた時に、(Av) \cdot (Aw)=v \cdot wが成り立つ。
      • 直交行列：直交行列とは、ある行列Aにおいて、A^{-1}=A^Tが成り立つ行列のことを指す。

  • いくつかある直交変換を用いた解法の中でも、時系列解析では、Householder法がいいらしい。

Householder法

  • U:任意の直交行列　を用いて、直交変換を行う。

    \begin{align}\|\varepsilon\|^2_N &= \|y-Za\|^2_N\\ &= \|U(y-Za)\|^2_N\\ &= \|Uy-UZa\|^2_N\\ \end{align}

  • そもそもHouseholder法とは？と思ったので、少し調べてみた。

    • Householder法には、Householder行列というものが必要。

      • Householder行列とは、下記のように表せる。

        H=E-\frac{2uu^T}{\|u\|^2}

        • Eは単位行列
        • u=x-y,xはある列ベクトルで、yはそれを鏡映変換したベクトル。
        • 鏡映変換：ベクトルの対称移動。何を基準に対称移動しているかというと、uというベクトルに直交するベクトルMを基準としている。

      • また、Householder行列には、下記の性質が成り立つ。

        H=H^T=H^{-1} \dashrightarrow 対称行列\bigcap直交行列

    • Householder行列を使った直交変換を行うのが、Householder変換であり、それを用いて、最小二乗法を求めるのが、Householder法。

    • どういったことができるのかは、こちらのURLが参考になるかと。

      https://cattech-lab.com/science-tools/lecture-mini-qr-decomposition/

  • そんなHouseholder法で式展開していくと、このようになる。（Sは変換後の値）

\begin{align}\|Uy-UZa\|^2_N &= \| \begin{bmatrix}S_{1,m+1} \\ \vdots \\ S_{1,m+1} \\ S_{1,m+1} \\ 0 \end{bmatrix} - \begin{bmatrix}S_{11} & \cdots & S_{11}\\ & \ddots & \vdots \\ & & S_{11} \\ && \\ &0 \end{bmatrix} \begin{bmatrix}a_1 \\ \vdots \\ a_m \end{bmatrix} \|^2_N\\ &= \|\begin{bmatrix}S_{1,m+1} \\ \vdots \\ S_{m,m+1} \end{bmatrix} - \begin{bmatrix}S_{11} & \cdots & S_{1m}\\ & \ddots & \vdots \\ & & S_{mm}\end{bmatrix} \begin{bmatrix}a_1 \\ \vdots \\ a_m \end{bmatrix} \|^2_m + S^2_{m+1,m+1} \end{align}
  • （8）のS^2_{m+1,m+1}については、今回求めたい回帰係数aには関係なく、無視できる。
  • つまり、ノルム内の式の値を最小にする、0にするという方程式を解くことが目標となる。よって、ノルム内の式を=0とおくと以下のように表すことができる。
\begin{bmatrix}S_{11} & \cdots & S_{1m}\\ & \ddots & \vdots \\ & & S_{mm}\end{bmatrix} \begin{bmatrix}a_1 \\ \vdots \\ a_m \end{bmatrix}=\begin{bmatrix}S_{1,m+1} \\ \vdots \\ S_{m,m+1} \end{bmatrix}
  • この式を解くには、後退代入という、上三角化された係数行列と右辺のベクトルを用いて、一番下（今回であれば、m）から解いて、上（今回であれば、1）までを解くという方法を用いる。（詳細については、動画をご参照ください。）

Householder法のメリット

  • ここまで長々と式を書き連ねてきたが、従来の方法より、何が有用であるのかを紹介する。
    1. 精度が良い。
    2. 計算が効率的になる。
    3. 上三角行列がモデル推定に必要な全情報を持つ。
      • 特に、データ数が10,000(=N)×100(=m)あるようなデータでも、一度その計算をしてしまえば、100×100になり、その後はそれをベースに考えられるので、再計算時のボリュームを抑えることができる。
    4. 時系列モデリングに便利
    5. データの併合、モデルの併合等が容易

部分回帰モデル

部分回帰モデルとは

  • 通常の回帰モデル同様、目的変数y_nと説明変数x_{nj}(j=1,\cdots,m)がある。この時、説明変数m個をすべて用いるのではなく、k個を選んで、モデルを推定する。
  • 次数kには、_mC_k個の部分回帰モデルがあり、組み合わせが多いと計算が大変。

Householder法による推定

  • 通常の多項式回帰では、変数選択時、定数項から順に変数を加えていって、その時々のAICの変化を見ていくが、部分回帰の場合は、別途説明変数の優先順位をつけて、その順に追加し、AICを比べていく。
  • ただ、この方法では、結局組み合わせの数に依存して、計算量が多くなってしまう。

Hocking-Leslieアルゴリズム

  • k個以下の変数の部分集合を削除したモデルの残差分散\sigma^2_kが、\sigma^2_{k+1}以下ならば、k+1個以上の変数を一つでも削除したモデルの残差分散が\sigma^2_kより小さくなることはない。
  • 上記のような性質が成り立つと、任意の変数選択を行った時の残差分散より大きくなるとわかっているものを計算するのは、意味がないため、それらの組み合わせを省くことができる。
  • これにより、次元数が多くなっても、計算量を大きく抑えることができる。
筆者のコメント
  • 詳細なアルゴリズムについて、説明したかったが、参考になる記事がなく、断念。詳しい方がいたら、教えて欲しい。

AICによるモデル選択上の注意

  1. AICの絶対値には意味がない。ただし、AICの差には意味がある。実際に、データをスケール変換したとしても、差に変化は生じない。
  2. AICは小数点以下２桁程度の表示で十分。それ以降はほとんど同じと見做せる。
  3. AICが最小点以上でほぼ２ずつ増加する場合、モデル族が真のモデルを含むことになる。
  4. 次数の増加に伴ってAICが減少し続ける場合、モデル族が真のモデルを含まない。

Discussion