修正済み赤池情報量規準(AICc)の導出

Hurvich and Tsai (1989)を流し読みして得た理解まとめ

準備

まず，データは次のTrue Modelから生成されたと仮定する．

$\boldsymbol{y} = \boldsymbol{\mu} + \boldsymbol{\epsilon}$
$\text{where}$
$~~~~ \boldsymbol{y} = (y_1, \cdots, y_n)^T,$
$~~~~ \boldsymbol{\mu} = ({\mu}_1, \cdots, {\mu}_n)^T,$
$~~~~ \boldsymbol{\epsilon} = ({\epsilon}_1, \cdots, {\epsilon}_n)^T, {\epsilon}_i \overset{\text{i.i.d.}}{\sim} N(0, {{\sigma}_0}^2).$

Hurvich and Tsai (1989)ではTrue Modelの意味でOperating Modelという用語が使われているが，どういうニュアンスを含んでいるのかよくわからない

また，近似に用いるモデル族として次のものを考える．

$\boldsymbol{y} = \boldsymbol{h}(\boldsymbol{\theta}) + \boldsymbol{u}$
$\text{where}$
$~~~~ \boldsymbol{h}(\boldsymbol{\theta}) = (h_1(\boldsymbol{\theta}), \cdots, h_n(\boldsymbol{\theta}))^T.$
$~~~~ \boldsymbol{u} = (u_1, \cdots, u_n)^T, u_i \overset{\text{i.i.d.}}{\sim} N(0, {\sigma}^2)$

$\boldsymbol{\theta}$ は $m$ 次元ベクタ， $\boldsymbol{h}$ は連続で2階微分可能であると仮定する．

近似に用いるモデルとTrue Modelの具体例として，両者が線形モデルであるとすると， $\boldsymbol{h}(\boldsymbol{\theta})=X\boldsymbol{\theta}$ ， $\boldsymbol{\mu}=X_0\boldsymbol{\theta}_0$ となる．ここで $X$ と $X_0$ は各々に $n \times m$ と $n \times m_0$ のフルランクの行列， $\boldsymbol{\theta}_0$ は $m_0$ 次元ベクタである．

近似に用いるモデルとTrue Modelの乖離を測る有用な指標として，カルバック・ライブラ情報量がある．

$\Delta(\boldsymbol{\theta}, {\sigma}^2) = E_F\{ -2 \text{log} (g_{\boldsymbol{\theta},{\sigma}^2}(\boldsymbol{y})) \}$

ここで $F$ は期待値がTrue Modelに関するものであることを示し， $g_{\boldsymbol{\theta},{\sigma}^2}(\boldsymbol{y})$ は近似に用いるモデルに基づいた尤度関数である．ここでのカルバック・ライブラ情報量は定数になる項を無視し， $2$ をかけている．

カルバック・ライブラ情報量が小さいほど，対象となったモデルとTrue Modelの乖離が少ない．

正規分布についての尤度関数を次に示す．

$\displaystyle \begin{align} g_{\boldsymbol{\theta},{\sigma}^2}(\boldsymbol{y}) &= \prod_{i=1}^n \left\{ (2 \pi {\sigma}^2)^{-\frac{1}{2}} \text{exp} \left\{ -\frac{(y_i-h_i(\boldsymbol{\theta}))^2}{2 {\sigma}^2 } \right\} \right\} \\ &= (2 \pi {\sigma}^2)^{-\frac{n}{2}} \text{exp} \left\{ -\frac{1}{2{\sigma}^2} \sum_{i=1}^n (y_i - h_i (\boldsymbol{\theta}))^2 \right\} \\ &= (2 \pi {\sigma}^2)^{-\frac{n}{2}} \text{exp} \left\{ -\frac{(\boldsymbol{y} - \boldsymbol{h}(\boldsymbol{\theta}))^T(\boldsymbol{y} - \boldsymbol{h}(\boldsymbol{\theta}))}{2 {\sigma}^2} \right\} \end{align}$

具体的に得られた尤度関数からカルバック・ライブラ情報量は次のようになる．

$\displaystyle \begin{align} \Delta (\boldsymbol{\theta}, {\sigma}^2) &= -2 E_F \left[ \text{log} \left\{ (2 \pi {\sigma}^2)^{-\frac{n}{2}} \text{exp} \left\{ -\frac{(\boldsymbol{y} - \boldsymbol{h}(\boldsymbol{\theta}))^T(\boldsymbol{y} - \boldsymbol{h}(\boldsymbol{\theta}))}{2 {\sigma}^2} \right\} \right\} \right] \\ &= n \text{log} (2 \pi {\sigma}^2) + \frac{1}{{\sigma}^2} E_F \left[ (\boldsymbol{\mu} + \boldsymbol{\epsilon} - \boldsymbol{h} (\boldsymbol{\theta}))^T(\boldsymbol{\mu} + \boldsymbol{\epsilon} - \boldsymbol{h} (\boldsymbol{\theta})) \right] \\ &= n \text{log} (2 \pi {\sigma}^2) + \frac{n {{\sigma}_0}^2}{{\sigma}^2} + \frac{(\boldsymbol{\mu} - \boldsymbol{h} (\boldsymbol{\theta}))^T(\boldsymbol{\mu} - \boldsymbol{h} (\boldsymbol{\theta}))}{{\sigma}^2} \end{align}$

補足蛇足
$\boldsymbol{\alpha} = \boldsymbol{\mu} - \boldsymbol{h}(\boldsymbol{\theta})$
$\displaystyle (\boldsymbol{\epsilon} - \boldsymbol{\alpha})^T(\boldsymbol{\epsilon} - \boldsymbol{\alpha}) = \sum_{i=1}^n ({\epsilon}_i + {\alpha}_i)^2 = \sum_{i=1}^n ({\epsilon}^2 +2 {\epsilon}_i \alpha + {{\alpha}_i}^2) = \sum_{i=1}^n {{\epsilon}_i}^2 + 2 \sum_{i=1}^n {\epsilon}_i {\alpha}_i + \sum_{i=1}^n {{\alpha}_i}^2$
$\displaystyle E\left[\sum_{i=1}^n {{\epsilon}_i}^2 \right] = n {{\sigma}_0}^2, ~~ E\left[2 \sum_{i=1}^n {\epsilon}_i {\alpha}_i \right] = 0, ~~ E\left[ \sum_{i=1}^n {{\alpha}_i}^2 \right] =(\boldsymbol{\mu} - \boldsymbol{h} (\boldsymbol{\theta}))^T(\boldsymbol{\mu} - \boldsymbol{h} (\boldsymbol{\theta}))$

妥当な規準として近似に用いるモデルを評価するため，カルバック・ライブラ情報量の期待値 $E_F \{ \Delta (\boldsymbol{\theta}, {\sigma}^2) \}$ を取り，パラメータとして得られたデータに基づいた最尤推定量 $(\hat{\boldsymbol{\theta}}, \hat{\sigma}^2)$ を使用する．

$\displaystyle \hat{\boldsymbol{\theta}} = \underset{\boldsymbol{\theta}}{\text{argmin}} (\boldsymbol{\mu} - \boldsymbol{h} (\boldsymbol{\theta}))^T(\boldsymbol{\mu} - \boldsymbol{h} (\boldsymbol{\theta}))$
$\displaystyle \hat{\sigma}^2 = \frac{1}{n}(\boldsymbol{\mu} - \boldsymbol{h} (\hat{\boldsymbol{\theta}}))^T(\boldsymbol{\mu} - \boldsymbol{h} (\hat{\boldsymbol{\theta}}))$

最尤推定量を適用してカルバック・ライブラ情報量は次のようになる．

$\displaystyle \Delta (\hat{\boldsymbol{\theta}}, \hat{\sigma}^2) = n \text{log} (2 \pi \hat{\sigma}^2) + \frac{n {{\sigma}_0}^2}{\hat{\sigma}^2} + \frac{(\boldsymbol{\mu} - \boldsymbol{h} (\hat{\boldsymbol{\theta}}))^T(\boldsymbol{\mu} - \boldsymbol{h} (\hat{\boldsymbol{\theta}}))}{\hat{\sigma}^2}$

Hurvich and Tsai (1989)では，この段階で $n\text{log}(2 \pi)$ が定数として無視されるが理由は不明． $n\text{log}(2 \pi)$ を無視すると尤度との関係がわかり難くなるので良いことはない気がするが…

近似に用いるモデル族が複数あるとき， $E_F \{ \Delta (\hat{\boldsymbol{\theta}}, \hat{\sigma}^2) \}$ を最小化するものが，ある意味でTrue Modelに最も近いものとして選択され得る．

実際にはTrue Modelは未知であるから，それに基づく $E_F \{ \Delta (\hat{\boldsymbol{\theta}}, \hat{\sigma}^2) \}$ も同様に未知である．しかし，何らかの追加の仮定が置かれることで，この量は推定され得る．推定され得る1例として，次に示す赤池情報量規準は近似的に $E_F \{ \Delta (\hat{\boldsymbol{\theta}}, \hat{\sigma}^2) \}$ の不偏推定量を提供する．

$AIC = n(\text{log} ( 2 \pi \hat{\sigma}^2) + 1) + 2 (m + 1)$

ここで $m$ は近似に用いるモデルのパラメータの次数である．また， $n(\text{log} ( 2 \pi \hat{\sigma}^2)+1)$ は最大対数尤度 $l \times -2$ と等しい．

AICcの導出

ここでの目標は次の $E_F \{ \Delta (\hat{\boldsymbol{\theta}}, \hat{\sigma}^2) \}$ の推定量を得ることであり，それが最終的にAICcを導く．

$\displaystyle \begin{align} E_F \{ \Delta (\hat{\boldsymbol{\theta}}, \hat{\sigma}^2) \} &= E_F \left\{ n \text{log} (2 \pi \hat{\sigma}^2) + \frac{n {{\sigma}_0}^2}{\hat{\sigma}^2} + \frac{(\boldsymbol{\mu} - \boldsymbol{h} (\hat{\boldsymbol{\theta}}))^T(\boldsymbol{\mu} - \boldsymbol{h} (\hat{\boldsymbol{\theta}}))}{\hat{\sigma}^2} \right\} \\ &= E_F \left\{ n \text{log} (2 \pi \hat{\sigma}^2) \right\} + E_F \left\{ \frac{n {{\sigma}_0}^2 }{\hat{\sigma}^2} \right\} + E_F \left\{ \frac{(\boldsymbol{\mu} - \boldsymbol{h}(\hat{\boldsymbol{\theta}}))^T (\boldsymbol{\mu} - \boldsymbol{h}(\hat{\boldsymbol{\theta}}))}{\hat{\sigma}^2} \right\} \end{align}$

まず， $(SSE / {{\sigma}_0}^2) \sim {\chi}^2(n-m)$ を定理として受け入れ， $\hat{\sigma}^2$ の定義より $SSE = n \hat{\sigma}^2$ であることから， $n \hat{\sigma}^2 / {{\sigma}_0}^2$ は ${\chi}^2(n-m)$ に従い， ${{\sigma}_0}^2/(n \hat{\sigma}^2)$ は $\text{Inv-}{\chi}^2(n-m)$ に従う．ここで証明なしに $(SSE / {{\sigma}_0}^2) \sim {\chi}^2(n-m)$ を定理として受け入れたが，これは調べればWikiにも出てくる；自由度あたりが浅学過ぎて証明はわからない．書籍として少なくともPuntanen et al. (2013)には載ってる．
このことから $E_F \{ \Delta (\hat{\boldsymbol{\theta}}, \hat{\sigma}^2) \}$ の右辺第2項は次のように求まる．

$\displaystyle E_F \left\{ \frac{n {{\sigma}_0}^2}{\hat{\sigma}^2} \right\} = n^2 E_F \left\{ \frac{{{\sigma}_0}^2}{n \hat{\sigma}^2} \right\} = \frac{n^2}{n-m-2}$

以下，近似に用いるモデル族にTrue Modelが含まれているという強い仮定を置く．この仮定の下でTrue Modelの平均応答 $\mu$ は $\mu = \boldsymbol{h} ( \boldsymbol{\theta}^{*} )$ となる．ここで $\boldsymbol{\theta}^{*}$ は $m \times 1$ 次元の未知のベクタである．
次に， $\boldsymbol{h} (\hat{\boldsymbol{\theta}})$ の $\boldsymbol{\theta} = \boldsymbol{\theta}^{*}$ におけるテイラ展開による1次近似を考える．

$\boldsymbol{h} (\hat{\boldsymbol{\theta}}) \approx \boldsymbol{h} (\boldsymbol{\theta}^{*}) + V (\hat{\boldsymbol{\theta}} - \boldsymbol{\theta}^{*})$

ここで $V$ は $\hat{\boldsymbol{\theta}} = \boldsymbol{\theta}^{*}$ における $\partial \boldsymbol{h} / \partial \boldsymbol{\theta}$ である．追加の仮定として線形モデルを考えると $V$ は計画行列になり， $\hat{\boldsymbol{\theta}} - \boldsymbol{\theta}^{*}$ は多変量正規分布 $N_m(\boldsymbol{0}, {{\sigma}_0}^2 (V^T V)^{-1})$ に従う(Madisen and Thyregod (2010), p.49, Theorem 3.3)．Hurvich and Tsai (1989)では $\boldsymbol{h}$ を線形モデルに限定した書き方はしていないが，非線形モデルの場合でも同様の正規性が導入できるのかは浅学過ぎてよくわからない．

以上から $E_F \{ \Delta (\hat{\boldsymbol{\theta}}, \hat{\sigma}^2) \}$ の右辺第3項の期待値内は近似的に $F$ 分布に帰着できる．

$\displaystyle \begin{align} \frac{(\boldsymbol{\mu} - \boldsymbol{h}(\hat{\boldsymbol{\theta}}))^T (\boldsymbol{\mu} - \boldsymbol{h}(\hat{\boldsymbol{\theta}}))}{\hat{\sigma}^2} &= \frac{(\boldsymbol{h}(\boldsymbol{\theta}^{*}) - \boldsymbol{h}(\hat{\boldsymbol{\theta}}))^T (\boldsymbol{h}(\boldsymbol{\theta}^{*}) - \boldsymbol{h}(\hat{\boldsymbol{\theta}}))}{\hat{\sigma}^2} \\ &\approx \frac{ ( \hat{\boldsymbol{\theta}} - \boldsymbol{\theta}^{*})^T V^T V ( \hat{\boldsymbol{\theta}} - \boldsymbol{\theta}^{*} ) }{ \hat{\sigma}^2 } \\ &= \frac{nm}{n-m} \frac{ ( ( \hat{\boldsymbol{\theta}} - \boldsymbol{\theta}^{*})^T V^T V ( \hat{\boldsymbol{\theta}} - \boldsymbol{\theta}^{*} ) / {{\sigma}_0}^2 ) (1/m) }{ (n \hat{\sigma}^2 / {{\sigma}_0}^2) (1/(n-m)) } \end{align}$

$\displaystyle \Rightarrow \frac{n-m}{nm} \frac{(\boldsymbol{\mu} - \boldsymbol{h}(\hat{\boldsymbol{\theta}}))^T (\boldsymbol{\mu} - \boldsymbol{h}(\hat{\boldsymbol{\theta}}))}{\hat{\sigma}^2} \approx \frac{ ( ( \hat{\boldsymbol{\theta}} - \boldsymbol{\theta}^{*})^T V^T V ( \hat{\boldsymbol{\theta}} - \boldsymbol{\theta}^{*} ) / {{\sigma}_0}^2 ) (1/m) }{ (n \hat{\sigma}^2 / {{\sigma}_0}^2) (1/(n-m)) } \sim F(m, n-m)$

ゆえに

$\displaystyle \begin{align} E_F \left\{ \frac{ (\boldsymbol{\mu} - \boldsymbol{h}( \hat{\boldsymbol{\theta}} ) )^T (\boldsymbol{\mu} - \boldsymbol{h}( \hat{\boldsymbol{\theta}} ) ) }{\hat{\sigma}^2} \right\} &= \frac{nm}{n-m} E_F \left\{ \frac{n-m}{nm} \frac{ (\boldsymbol{\mu} - \boldsymbol{h}( \hat{\boldsymbol{\theta}} ) )^T (\boldsymbol{\mu} - \boldsymbol{h}( \hat{\boldsymbol{\theta}} ) ) }{\hat{\sigma}^2} \right\} \\ &{ \approx ~} \frac{nm}{n-m} E_F \left\{ \frac{ ( ( \hat{\boldsymbol{\theta}} - \boldsymbol{\theta}^{*})^T V^T V ( \hat{\boldsymbol{\theta}} - \boldsymbol{\theta}^{*} ) / {{\sigma}_0}^2 ) (1/m) }{ (n \hat{\sigma}^2 / {{\sigma}_0}^2) (1/(n-m)) } \right\} \\ &= \frac{nm}{n-m} \frac{n-m}{n-m-2} = \frac{nm}{n - m - 2} \end{align}$

最終的に全体としてまとめると $E_F \{ \Delta (\hat{\boldsymbol{\theta}}, \hat{\sigma}^2) \}$ は次のようになる．

$\displaystyle \begin{align} E_F \{ \Delta (\hat{\boldsymbol{\theta}}, \hat{\sigma}^2) \} ~ &{\approx} ~ E_F \{ n \text{log} (2 \pi \hat{\sigma}^2) \} + \frac{n^2}{n-m-2} + \frac{nm}{n - m - 2} \\ &= E_F \{ n ( \text{log} (2 \pi \hat{\sigma}^2) + 1) \} + \frac{2 n (m + 1) }{n - m - 2} \end{align}$

そしてここからAICcは次のように定義される．

$\displaystyle \text{AICc} = n ( \text{log} (2 \pi \hat{\sigma}^2) + 1) + \frac{2 n (m + 1) }{n - m - 2} = \text{AIC} + \frac{2 (m + 1) (m + 2) }{n - m - 2}$

AICcの期待値は $E_F \{ \Delta (\hat{\boldsymbol{\theta}}, \hat{\sigma}^2) \}$ を近似する．

$\displaystyle E_F \{ \text{AICc} \} = E_F \{ n ( \text{log} (2 \pi \hat{\sigma}^2) + 1) \} + \frac{2 n (m + 1) }{n - m - 2} \approx E_F \{ \Delta (\hat{\boldsymbol{\theta}}, \hat{\sigma}^2) \}$

上のAICcの罰則項は少し見慣れないかもしれないが，正規線形モデル等で推定される分散もはじめから考慮に入れたパラメータ数として $k=m+1$ とすれば，比較的よく見かけるAICcになる．

$\displaystyle \text{AICc} = n (\text{log} (2 \pi \hat{\sigma}^2) + 1) + \frac{2 n k}{n - k - 1} = \text{AIC} + \frac{2 k (k + 1) }{n - k - 1}$

前のBlogでAICcの利用法に関して，情報量規準で意図される罰則項という観点から述べたが，AICcの導出過程から罰則項の分母は $\text{Inv-}{\chi}^2$ 分布と $F$ 分布の自由度に由来するものであるため，明確に $n>m+1$ とならない場合の使用は悪用である．

今回も思いの外長くなった．Hurvich and Tsai (1989)を読んでざっくりは理解できた気がしていたが，これを書いてて理解の穴が多分にあることに気がついた．さらなる研鑽が必要である(小並感)．なにか改善点・間違いに気がついたら修正する．