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

2022/12/22に公開約7,800字

統計的モデリングと情報量基準AIC

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

統計的モデリングとモデル評価


統計的モデリングとは

  • 統計的モデル①:情報抽出の「道具」である。

    • モデルを作成するために必要な要素(入力)
      • 理論
      • 経験的知識(ドメイン知識)
      • データ
    • 作成したモデルから得られる要素(出力)
      • 情報抽出 / 知識発見
      • 予測シミュレーション
      • 制御 / 意思決定
  • 統計的モデル②:「確率分布」で表現する。

    • 基本的な形

      連続型確率分布 離散型確率分布
      正規分布 コーシー分布 ピアソン分布 二項分布
      指数分布 二重指数分布 \chi{^2}分布 多項分布
      一様分布 ポアソン分布
    • プラスアルファの様々な情報を取り入れる必要がある。

      • 他の変数の情報
      • 過去の変動の情報
      • 時間経過による変化
  • ①、②共に、いかに有効活用するかで得られるモデルは異なる。特に、それぞれの課題設定にあった出力となっているかは重要である。

モデル評価

  • モデリング時に様々な情報から色々なモデルを構築することになる。
  • そのためにもどのモデルが適しているか評価し選択することが重要。
  • その尺度として用いられるのが情報量基準である。

Kullback-Leibler情報量


情報量基準とは

  1. ”モデルの良さ”を予測能力で評価する。
  2. 予測は点推定ではなく”予測分布”で行う。(信頼区間とか?)
  3. 分布の近さを”K-L情報量”で評価する。

K-L情報量によるモデルの評価

g(y):真の分布\\ f(y):モデルの分布
  • Kullback-Leibler情報量:真の分布とモデルの分布の乖離具合を測る尺度

    \begin{align}I(g;f) &= E_Y\log{[{\frac{g(Y)}{f(Y)}}]}\\ &=\int{\log{[\frac{g(x)}{f(x)}}}]dG(x)\\ &= \begin{cases}\int{\log{[\frac{g(x)}{f(x)}]}g(x)dx}~~連続分布モデル\\ ~~ \\ \sum_{i \in J}{g_i \log{[\frac{g_i}{f_i}]}}~~離散分布モデル\end{cases} \end{align}

  • K-L情報量の性質

    \begin{align}&(i)~I(g;f) ≥ 0\\ &(ii)~I(g;f)=0 \iff g(x)=f(x)\end{align}

  • K-L情報量は距離ではないことに注意。そのため、「対称性がある」という距離の性質は成り立たないということになる。d(g,f)≠d(f,g)


KL情報量の推定と最尤法


K-L情報量の推定

  • そもそもK-L情報量に用いられるg(y)というのは未知であるため、公式通りに計算というのはできない。
  • そのため、I(g(y);f(y))をデータから推定し、求めることになる。
  • その推定をするにあたり、K-L情報量の便利な性質がある。それが\logを用いていることである。そのため、gfの期待値を分離することができる。
    I(g;f) = E_Y\log{\frac{g}{f}} = E_Y\log{g}-E_Y\log{f}
  • この等式の中で注目すべき点がE_Y\log{f(Y)}の部分である。この項は、平均対数尤度と呼ばれる。
  • また、E_Y\log{g}についても便利な性質がある。それは、道ではあるもののfに関係することなく一定であるという点。
  • そういった点から、モデリング時には、K-L情報量を求める代わりに、平均対数尤度を使えるらしい。数式で関係を表すと下記の通り。
    E_Y\log{f(Y)}:大 \iff I(g;f):小
  • ちなみに、K-L情報量は0に近くづくほど良いです。

平均対数尤度E\log{f}の推定

  • 平均対数尤度を求めることで、モデリングを評価し、選択につながるということであったが、平均対数尤度においても未知の分布を含んだ数式になっている。
    E_Y{\log{f(Y)}}=\int{f(y)dG(y)}
  • 上記の式におけるG(y)の部分が未知であり、推定の対象になる。(dG(y)=g(y)dy
    G(y) \rightarrow \hat{G}_n(y) = \frac{1}{n}\sum_{i=1}^{n}{I(y,X_i)}
  • この式を対数尤度として求めると、、、(nは元が平均の式なので、それを戻すためにつけている。)
l = n\int\log{f(y)}d\hat{G}_n(y) = \sum_{i=1}^{n}\log{f(X_i)}
  • 上記式を改めて平均化し、大数の法則を当てはめると、平均対数尤度に収束していく。
    \frac{1}{n}\sum_{i=1}^{n}\log{f(X_i)} \rightarrow E_Y{\log{f(Y)}}

最尤法

  • 今までのf(y)にパラメータ\thetaを与える。
    f(y~|~\theta), ~~ \theta \equiv (\theta_1, \cdots, \theta_k)'
  • この式を対象に最尤法を適応する。最尤法は近似的にK-L情報量を最小化することになる。
    対数尤度:l(\theta) = \sum_{i=1}^n{\log{f(X_i~|~\theta)}} \equiv \log{f(X~|~\theta)}

    最尤法:\max_\theta{l(\theta)} \rightarrow \hat\theta = \hat\theta(X)
    • 主な用語
      用語 数式
      最尤推定量 \hat\theta
      最大対数尤度 l(\hat\theta)
      最尤モデル f(y~\lvert ~\hat\theta)

最尤推定値の求め方

  1. 尤度方程式を解く
  2. 数値的最適化を用いる --> 時系列解析においては、データの性質上、主にこちら

最尤推定量の性質

  1. 尤度方程式\frac{\partial l(\theta)}{\partial\theta}=0\theta_0に収束する解をもつ
  2. \hat\theta_nn \rightarrow+\infinの時、\theta_0に確率収束
  3. \sqrt{n}(\hat\theta_n - \theta_0) \rightarrow N(0, I(\theta_0)^{-1})(Iはフィッシャー情報量)

バイアス補正とAICの導出


複数のモデルの比較

  • 出てくる用語
名称 数式
モデル M_1, \cdots,M_k
パラメータ \theta_1,\cdots,\theta_k
最大対数尤度 l_1(\hat\theta_1), \cdots, l_k(\hat\theta_k)
Fisher情報行列 I(・)
Hessianの期待値 J(・)
  • k個のモデルの中でどのモデルが有用かを最大対数尤度を用いて比較すれば良いのか?

    • 残念ながらそういうわけにはいかない。なぜなら、低次数モデルと高次数モデルでは、高次数モデルの方が良い結果を出してしまうため。
    • 通常の機械学習同様、高次数モデルが予測精度が高いのかというと、そういうわけではない。
  • 理由・原因・対策

項目 内容
理由 最大対数尤度l(\hat\theta)が平均対数尤度E\log{f(x~\lvert~\hat\theta)}の推定値としてバイアスを持ち、しかもバイアス量がモデルによって異なるため。
原因 パラメータ推定と平均対数尤度の推定に同じデータを2回用いたため。
対策 バイアスを評価し補正する。
  • どのようにバイアスを補正するのか?それは、バイアス補正をするための項b(G)を対数尤度加える(数式的にはマイナス)。

    • b(G)は、(対数尤度の平均値-平均対数尤度)の平均値というものなので、平均から戻すためにバイアス補正項にnがかけられる。
      \log{f(X~\lvert~\hat\theta(X))}-nb(G)
  • b(G)については、下記のように分解が可能である。詳細については、元の動画を参照して欲しいが、Dというのが対数尤度と平均対数尤度とのバイアスによるズレとしている。また、Dは3つに分解できる。

    b(G) = E_X[D] = E_X[D_1 + D_2 + D_3]

それぞれの点
点A 対数尤度グラフ上の\hat\theta
点B 対数尤度グラフ上の\theta_0
点C 平均対数尤度グラフ上の\theta_0
点D 平均対数尤度グラフ上の\hat\theta
D1・D2・D3のイメージ
D_1 点A-点B
D_2 点B-点C
D_3 点C-点D
  • D1・D2・D3のそれぞれの期待値を求めて、それを代入すると、下記のような式になる。

    b(G)=\frac{1}{n}tr[I(G)J(G)^{-1}]

  • この値を情報量基準ICの一般形の式に代入すると、TICが得られる。

    IC = -2\log{f(x~\lvert~\hat\theta)}+2nb(G)

    TIC = -2\log{f(x~\lvert~\hat\theta)}+2tr[I(G)J(G)^{-1}]

  • ただし、ここから使うのはAICとなる。詳しくは次の項目で。

    • kは自由パラメータ数(おそらくAIC固有の概念)であり、モデルの次元数を意味する。
      AIC = -2\log{f(x~\lvert~\hat\theta)}+2k

AICとTICの関係


nb(G)=tr[I(G)J(G)^{-1}]kになる理由

  • 複数作成したモデル(モデル族)内に真の分布を含む場合、下記の関係が成り立つらしい。

    I(\theta)=J(\theta)

  • 上記の性質が満たされるのであれば、下記の式が得られる。

    • この時、I_kは単位行列である。そのトレース(対角成分の和)をとるので、スカラーとなる。
      \begin{align}nb(G) &= tr[I(G)J(G)^{-1}]\\ &= tr[I_k]\\ &= k\end{align}
  • ちなみに、TICの方が、AICに比べて、多くを含んだ形で表現できているような気がする。(I(\theta)=J(\theta)という制約が必要ないため)

    • ただ、これは誤解で、I(\theta),J(\theta)ともに、推定をして求めることになるため、誤差が生じる。
    • また、高次モーメントを含んでいるため、分散が大きくなってしまう。

AICによるモデルの選択例


  • どういったふうに選択するのかは、動画を見て欲しい。基準としては、AICの値がいかに小さいか。その時の次数はいくつかというところである。

モデルの予測誤差分散

  • 予測誤差については、以下の等式が成り立つ。

    • バイアス:モデルの悪さ
    • 分散:モデルの不安定さ
      予測誤差=バイアス+分散
  • つまり、AICが最小のモデルというのは、、、

    • バイアスと分散を適度に小さくしたモデル
    • それは、期待予測誤差を最小にするモデルを選択したことになる。

その他の情報量基準


  • 情報量基準というのは、K-L情報量のバイアスを補正して導くもの。

AIC_c:情報量基準の有限修正

  • AICはデータ量が少ない時やデータ数に比べてパラメータ数が多い時、バイアスの補正がうまくいかない性質がある。
  • それへの対策として、有限修正を施す。具体的にはバイアス補正項を置き換える。そうすることで、データ少ない場合でもある程度の精度を出せるようにする。
    nb_c(G) = \frac{n(p+1)}{n-p-2}

    AIC_c = -2\log{f(X~\lvert~\hat\theta(X))}+2\frac{n(p+1)}{n-p-2}

GIC:一般情報量基準

  • 特徴として、統計的汎関数として定義できる任意の推定量に適用可能。つまり、最尤推定量に限らず、他の推定量を用いることができる。
    • 例)ロバスト推定量、正則化法、ベイズ推定の一部
  • また、高次モーメントを含むようなものでもうまく補正してくれる。
  • ただし、計算が複雑で面倒。
    GIC = -2\log{f(X~\lvert~\hat\theta(X))}+2nb(\hat{G}_{GIC})

EIC:Bootstrap 情報量基準

  • バイアス補正量を解析的にではなく、ブートストラップによって数値的に求める。

  • 特長としては以下3つ。

    1. 解析的近似が不要(今までのような計算)
    2. 計算実装も比較的容易
    3. 最尤推定量以外の様々な推定量やモデルに適用可能
  • 弱点としては、ブートストラップ法のため、データ生成・推定を繰り返す。すると、計算量が多くなってしまう。

ABIC:ベイズ型情報量基準

  • 良いモデルを求める方法には、モデル選択だけでなく、「パラメータに制約を課す」という方法が存在する。

  • そこで登場するのがベイズモデル。

  • f(x~\lvert~\theta)における\thetaに事前分布\pi{(\theta~|~\lambda)}を設けるというもの。\lambdaはハイパーパラメータ。これらを用いて、xの周辺分布を数式化する。

    p(x~|~\lambda)=\int{f(x~|~\theta)\pi(\theta~|~\lambda)d\theta}

  • この式をパラメータを\lambdaとするモデルとして、情報量基準を作成する。(p\lambdaの次元数)

    \begin{align} ABIC &= -2\log[\max_\lambda{p(x|\lambda)}]+2q\\ &= -2\max_{\lambda}\log{\int{f(x~|~\theta)\pi(\theta~|~\lambda)d\theta}} +2q \end{align}

  • 注意事項としては、\lambdaをパラメータに据えているので、本来の\thetaの良さを直接的に測っていることにはならない。

BIC

  • 周辺尤度を用いて、それを近似的に最大化するようなモデルを求めるというもの。これまでの情報量基準とは少しロジックが異なるらしい。
  • 周辺尤度については、事前分布\pi(\theta_m)を用いて、以下のように表すことができる。
    p(x) = \int{f(x|\theta_m)\pi(\theta_m)d\theta_m}
  • 以上を用いて、BICを表現すると、、、(\thickapprox:近似値)
    \begin{align}BICm = -2\log{p(x)} &=-2\log{\int{f(x|\theta_m)\pi(\theta_m)d\theta_m}}\\ &\thickapprox -2\log{f(x|\hat\theta)}+m\log{N} \end{align}

交差検証法(Cross Validation)

  1. 全データを推定用データと評価用データに分ける
  2. 推定用データでモデルを推定
  3. 評価用データでモデルを評価(予測2乗誤差など)
  4. 1の分割の仕方を変えて、すべての場合について、2,3を繰り返し、評価量の平均を求める。
  • 分割の仕方
    • Leave-one-out:1個のデータだけ評価に用いる
    • k分割法:全体のデータをk分割し、そのうちの1つを評価に用いる

AICに関する批判

  • よくある批判として「AICには真の次数への一致性がない」ということが挙げられる。
    • 一致性:サンプル数を増やすことで、本来のものに近似する
  • ただ、それは、AICにとって最も重要な問題ではない。以下がその理由。
    1. モデリングの目的は、「良い」モデルを求めることで、「真の」モデルを求めることではない。
    2. 次数の一致性は良いモデルを求めるための「必要条件」でも「十分条件」でもない。
    3. 「真」の次数は一般に存在しない。存在したとしても、真の次数が推定されたモデルが予測にも良いとは限らない。
    4. 「真」の次数より高くてもパラメータが一致性を持てば、モデルは一致する。

Discussion

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