「大人の教養・知識・気付き」を伸ばすブログ

一流の大人(ビジネスマン、政治家、リーダー…)として知っておきたい、教養・社会動向を意外なところから取り上げ学ぶことで“気付く力”を伸ばすブログです。データ分析・語学に力点を置いています。 →現在、コンサルタントの雛になるべく、少しずつ勉強中です(※2024年1月21日改訂)。

MENU

統計的機械学習の数理100問(13/20)

 いい加減時代の潮流に乗ろうということで機械学習を学びたいと思う。またRはともかくとしてPythonは未経験であるため、丁度良い書籍として

を用いることにする。

4. 情報量基準

  • モデル選択において適合性と簡潔性のバランスが取れた統計モデルの指標である情報量基準を扱う。

4.1 情報量基準

 情報量基準は観測データから統計モデルの妥当性を評価する場合の指標として定義される。情報量基準は統計モデルによるデータの説明力(適合力)と統計モデルの簡潔度合い(簡潔性)*1の両方を評価する。
 線形回帰を想定する*2と、n組の観測値(\boldsymbol{X}_i,Y_i)\in\mathbb{R}^{p}\times\mathbb{R}を用いてp個の変数候補から説明変数として実際に採用する変数を選ぶ問題として定式化される。すなわちN=\{1,2,\cdots,n\}についてS\in2^{N}を選ぶ問題として定式化される。このSを用いた場合の残差変動RSS(S)が適合性を、変数の個数k(S)=|S|が簡潔性を表す。このとき、\mathrm{AIC},\mathrm{BIC}はそれぞれ


\begin{aligned}
\mathrm{AIC}&=n\log\hat{\sigma}_k^2+2k,\\
\mathrm{BIC}&=n\log\hat{\sigma}_k^2+k\log n
\end{aligned}

により定義される。

4.2 有効推定量とFisher情報量行列

 \mathrm{AIC}を導出すべく最小二乗法による線形回帰の推定方法が不偏推定量の中でも分散を最小にする有効推定量であることを導く。
 観測データ(Y_1,\boldsymbol{X}_1),\cdots,(Y_n,\boldsymbol{X}_n),\boldsymbol{X}_i\in\mathbb{R}^{p+1},Y_i\in\mathbb{R}が定数\boldsymbol{\beta}\in\mathbb{R}^{p+1}および確率変数\varepsilon_1,\cdots,\varepsilon_n\sim N(0,\sigma^2)を用いて


\begin{aligned}
Y_i={}^{t}\boldsymbol{X}_i\boldsymbol{\beta}+\varepsilon_i
\end{aligned}

で表されるものとする。すなわち


\begin{aligned}
f(y|\boldsymbol{x},\boldsymbol{\beta})=\displaystyle{\frac{1}{\sqrt{(2\pi)^{\frac{p}{2}}}}\exp\left\{-\frac{1}{2\sigma^2}(y-{}^{t}\boldsymbol{x}\boldsymbol{\beta})^2\right\}}
\end{aligned}

である。
 ここで最尤法、すなわち尤度


\begin{aligned}
L(\boldsymbol{\beta})=\displaystyle{\prod_{i=1}^{n}f(y_i|\boldsymbol{x}_i,\boldsymbol{\beta})}
\end{aligned}

を最大にするような\boldsymbol{\beta}を考える。計算を簡単にすべく、また対数y=\log xxに関して単調増加であることを踏まえ、対数尤度


\begin{aligned}
l=\log L(\boldsymbol{\beta})=\displaystyle{\log\left(\prod_{i=1}^{n}f(y_i|\boldsymbol{x}_i,\boldsymbol{\beta})\right)}
\end{aligned}

を最大化することを考えることにする。すると


\begin{aligned}
l(\beta)=-\displaystyle{\frac{n}{2}\log(2\pi\sigma^2)}-\displaystyle{\frac{1}{2\sigma^2}}(Y-{}^{t}\boldsymbol{X}\boldsymbol{\beta})^2
\end{aligned}

であり、\sigma^2を固定すれば、対数尤度の最大化と(Y-\boldsymbol{X}\boldsymbol{\beta})^2の最小化は同値となる。また対数尤度l\sigma^2偏微分すれば


\begin{aligned}
\displaystyle{\frac{\partial l}{\partial\sigma^2}}=-\displaystyle{\frac{n}{2\sigma^2}}+\displaystyle{\frac{1}{2(\sigma^2)^2}}(Y-{}^{t}\boldsymbol{X}\boldsymbol{\beta})^2
\end{aligned}

であるから、\hat{\boldsymbol{\beta}}=({}^{t}\boldsymbol{X}\boldsymbol{X})^{-1}{}^{t}\boldsymbol{X}Yを代入することで\sigma^2最尤推定\hat{\sigma}^2は上記の偏微分0とおいて


\begin{aligned}
\hat{\sigma}^2=\displaystyle{\frac{1}{n}}(Y-{}^{t}\boldsymbol{X}\hat{\boldsymbol{\beta}})^2
\end{aligned}

を得る。
 次に\hat{\boldsymbol{\beta}}が有効推定量であることを示す。一般に対数尤度l\boldsymbol{\beta}の各成分で微分した値から成り立つベクトル\nabla lの分散共分散行列をnで割ったJ\mathrm{Fisher}情報量行列という。


\begin{aligned}
\nabla l=\displaystyle{\frac{\displaystyle{\nabla\prod_{i=1}^{n}f(y_i|\boldsymbol{x}_i,\boldsymbol{\beta})}}{\displaystyle{\prod_{i=1}^{n}f(y_i|\boldsymbol{x}_i,\boldsymbol{\beta})}}}
\end{aligned}

が成立する。
 もし\nablaによる微分Yによる積分が可換ならば、


\begin{aligned}
\displaystyle{\int \displaystyle{\prod_{i=1}^{n}f(y_i|\boldsymbol{x}_i,\boldsymbol{\beta})}}dy=1
\end{aligned}

の両辺を\boldsymbol{\beta}偏微分することで


\begin{aligned}
\displaystyle{\int\nabla\displaystyle{\prod_{i=1}^{n}f(y_i|\boldsymbol{x}_i,\boldsymbol{\beta})}}dy=0
\end{aligned}

を得る。
 また


\begin{aligned}
E\left[\nabla l\right]&=\displaystyle{\int\frac{\nabla\displaystyle{\prod_{i=1}^{n}f(y_i|\boldsymbol{x}_i,\boldsymbol{\beta})}}{\displaystyle{\prod_{i=1}^{n}f(y_i|\boldsymbol{x}_i,\boldsymbol{\beta})}} }dy\\
&=\displaystyle{\int\nabla\displaystyle{\prod_{i=1}^{n}f(y_i|\boldsymbol{x}_i,\boldsymbol{\beta})}}dy=0
\end{aligned}

および


\begin{aligned}
0=\nabla\left[E[\nabla l]\right]&=\nabla\displaystyle{\int(\nabla l)\displaystyle{\prod_{i=1}^{n}f(y_i|\boldsymbol{x}_i,\boldsymbol{\beta})}}dy\\
&=\displaystyle{\int(\nabla^2 l)\prod_{i=1}^{n}f(y_i|\boldsymbol{x}_i,\boldsymbol{\beta})}dy+\displaystyle{\int(\nabla l)\left\{\nabla \displaystyle{\prod_{i=1}^{n}f(y_i|\boldsymbol{x}_i,\boldsymbol{\beta})}\right\}}dy\\
&=E\left[\nabla^2 l\right]+E[(\nabla l)^2]
\end{aligned}

が成り立つ。これは


\begin{aligned}
J=\displaystyle{\frac{1}{n}V[\nabla l]}=\displaystyle{\frac{1}{n}E[(\nabla l)^2]}=-\displaystyle{\frac{1}{n}E[\nabla ^2l]}
\end{aligned}

を意味する。


Cramér-Raoの不等式 任意の不偏推定量\tilde{\boldsymbol{\beta}}に関する分散共分散行列\mathbb{V}[\tilde{\boldsymbol{\beta}}]\in\mathbb{R}^{(p+1)\times(p+1)}について以下が成り立つ:

\begin{aligned}
\mathbb{V}[\tilde{\boldsymbol{\beta}}]\geq(n J)^{-1}
\end{aligned}

が成り立つ。

(\because まず

\begin{aligned}
\beta_i=\displaystyle{\int\tilde{\beta}_i\prod_{i=1}^{n}f(y_i|\boldsymbol{x}_i,\boldsymbol{\beta})}dy_i
\end{aligned}

の両辺を\beta_j偏微分すると


\begin{aligned}
\displaystyle{\int\tilde{\beta}_i\frac{\partial}{\partial\beta_j}\prod_{i=1}^{n}f(y_i|\boldsymbol{x}_i,\boldsymbol{\beta})}dy_i=\begin{cases}
1,&i=j,\\
0,&i\neq j
\end{cases}
\end{aligned}

である。これを分散共分散行列の形式で書けば


\begin{aligned}
\displaystyle{\int\tilde{\beta}{}^{t}\left\{\displaystyle{\frac{\nabla \displaystyle{\prod_{i=1}^{n}f(y_i|\boldsymbol{x}_i,\boldsymbol{\beta})}}{\displaystyle{\prod_{i=1}^{n}f(y_i|\boldsymbol{x}_i,\boldsymbol{\beta})}}\displaystyle{\prod_{i=1}^{n}f(y_i|\boldsymbol{x}_i,\boldsymbol{\beta})}}\right\}}=E\left[\tilde{\boldsymbol{\beta}}{}^{t}(\nabla l)\right]=I
\end{aligned}

となる。
またE[\nabla l]=0より、これは


\begin{aligned}
E[(\tilde{\boldsymbol{\beta}}-\boldsymbol{\beta}){}^{t}(\nabla l)]=I
\end{aligned}

を意味する。そこで\tilde{\boldsymbol{\beta}}-\boldsymbol{\beta}\nabla lをつなぎ合わせたサイズ2(p+1)の確率ベクトルの分散共分散行列は


\begin{aligned}
\begin{bmatrix}
V[\tilde{\boldsymbol{\beta}}]&I\\
I&n J
\end{bmatrix}
\end{aligned}

である。ここでV[\tilde{\boldsymbol{\beta}}],Jはともに分散共分散行列であるから非負である。したがって


\begin{aligned}
\begin{bmatrix}
V[\tilde{\boldsymbol{\beta}}]-(nJ)^{-1}&O\\
O&n J
\end{bmatrix}=\begin{bmatrix}
I&-(nJ)^{-1}\\
O&I
\end{bmatrix}\begin{bmatrix}
V[\tilde{\boldsymbol{\beta}}]&I\\
I&nJ
\end{bmatrix}=\begin{bmatrix}
I&O\\
\ -(nJ)^{-1}&I
\end{bmatrix}
\end{aligned}

である。右辺が非負定値であるから、左辺も非負定値である。任意の\boldsymbol{x}\in\mathbb{R}^nについて{}^{t}\boldsymbol{X}A\boldsymbol{X}\geq0ならば任意のB\in\mathbb{R}^{n\times m}\boldsymbol{Y}\in\mathbb{R}^{m}について


\begin{aligned}
{}^{t}\boldsymbol{Y}{}^{t}BAB\boldsymbol{Y}\geq0
\end{aligned}

が成立する。これは\mathbb{V}[\tilde{\boldsymbol{\beta}}]-(nJ)^{-1}が非負定値であることを意味する。実際、{}^{\forall}\boldsymbol{X},{}^{\forall}\boldsymbol{Y}\in\mathbb{R}^{p+1}について


\begin{aligned}
{}^{t}\boldsymbol{X}\left\{\mathbb{V}\left[\tilde{\boldsymbol{\beta}}\right]-(nJ)^{-1}\right\}\boldsymbol{X}+{}^{t}\boldsymbol{Y}J\boldsymbol{Y}\geq0
\end{aligned}

ならば、\boldsymbol{Y}=\boldsymbol{0}においても不等号が成り立つ。 \blacksquare)

4.3 Kullback-Leibler情報量

 \mathbb{R}上の確率密度関数f,gに対して


\begin{aligned}
D(f,g)=\displaystyle{\int_{-\infty}^{\infty}f(x)\log\frac{f(x)}{g(x)}}dx
\end{aligned}

\mathrm{Kullback}-\mathrm{Leibler}情報量といい、任意の\mathbb{R}上の\mathrm{Borel}集合について


\begin{aligned}
\displaystyle{\int_S f(x)dx}\gt0\Rightarrow\displaystyle{\int_S g(x)dx}\gt0,
\end{aligned}

すなわち絶対連続ならばDは定義できる。


線形回帰の\mathrm{Kullback}-\mathrm{Leibler}情報量 標本(Z_i,\boldsymbol{X}_i)に対して、母数\boldsymbol{\beta}\in\mathbb{R}^{p+1}に関する対数尤度

\begin{aligned}
\displaystyle{-\sum_{i=1}^{n}\log f(z_i|\boldsymbol{X}_i,\boldsymbol{\beta})}
\end{aligned}

は任意の\boldsymbol{\beta}^{\prime}\in\mathbb{R}^{p+1}について


\begin{aligned}
&\displaystyle{\frac{n}{2}\log(2\pi\sigma^2)}+\displaystyle{\frac{1}{2\sigma^2}(Z-\boldsymbol{X}\boldsymbol{\beta})^2}\\
&-\displaystyle{\frac{1}{\sigma^2}{}^{t}(\boldsymbol{\beta}-\boldsymbol{\beta}^{\prime})}{}^{t}\boldsymbol{X}(Z-\boldsymbol{X}\boldsymbol{\beta}^{\prime})\\
&+\displaystyle{\frac{1}{2\sigma^2}{}^{t}(\boldsymbol{\beta}-\boldsymbol{\beta}^{\prime})}{}^{t}\boldsymbol{X}\boldsymbol{X}(\boldsymbol{\beta}-\boldsymbol{\beta}^{\prime})
\end{aligned}

が成り立つ。

(\because u\in\mathbb{R},\boldsymbol{x}\in\mathbb{R}^{p+1},\boldsymbol{\beta},\boldsymbol{\beta}^{\prime}\in\mathbb{R}^{p+1}に対して

\begin{aligned}
\log f(u|\boldsymbol{x},\boldsymbol{\beta})&=-\displaystyle{\frac{1}{2}\log(2\pi\sigma^2)}-\displaystyle{\frac{1}{2\sigma^2}(u-\boldsymbol{x}\boldsymbol{\beta})^2},\\
(u-\boldsymbol{x}\boldsymbol{\beta})^2=&\{(u-\boldsymbol{x}\boldsymbol{\beta})-\boldsymbol{x}(\boldsymbol{\beta}-\boldsymbol{\beta}^{\prime})\}^2\\
=&(U-\boldsymbol{x}\boldsymbol{\beta})^2-2{}^{t}(\boldsymbol{\beta}-\boldsymbol{\beta}^{\prime}){}^{t}\boldsymbol{X}(u-\boldsymbol{x}\boldsymbol{\beta})\\
&+{}^{t}(\boldsymbol{\beta}-\boldsymbol{\beta}^{\prime}){}^{t}\boldsymbol{x}(\boldsymbol{\beta}-\boldsymbol{\beta}^{\prime})
\end{aligned}

である。したがって


\begin{aligned}
\log f(u|\boldsymbol{x},\boldsymbol{\beta})=&-\displaystyle{\frac{1}{2}}\log(2\pi\sigma^2)-\displaystyle{\frac{1}{2\sigma^2}}(u-\boldsymbol{x}\boldsymbol{\beta})^2\\
&+\displaystyle{\frac{1}{\sigma^2}}{}^{t}(\boldsymbol{\beta}-\boldsymbol{\beta}^{\prime}){}^{t}\boldsymbol{x}(u-{}^{t}\boldsymbol{x}\boldsymbol{\beta})\\
&-\displaystyle{\frac{1}{2\sigma^2}{}^{t}(\boldsymbol{\beta}-\boldsymbol{\beta}^{\prime}){}^{t}\boldsymbol{x}\boldsymbol{x}(\boldsymbol{\beta}-\boldsymbol{\beta}^{\prime})}
\end{aligned}

が成立し、i=1,2,\cdots,nに関して総和を取ることで


\begin{aligned}
\displaystyle{-\sum_{i=1}^{n}\log f(z_i|\boldsymbol{x}_i,\boldsymbol{\beta})}=&\displaystyle{\frac{n}{2}\log(2\pi\sigma^2)}+\displaystyle{\frac{1}{2\sigma^2}(z-\boldsymbol{x}\boldsymbol{\beta})}\\
&-\displaystyle{\frac{1}{\sigma^2}{}^{t}(\boldsymbol{\beta}-\boldsymbol{\beta}^{\prime}){}^{t}\boldsymbol{x}(u-{}^{t}\boldsymbol{x}\boldsymbol{\beta})}\\
&+\displaystyle{\frac{1}{2\sigma^2}{}^{t}(\boldsymbol{\beta}-\boldsymbol{\beta}^{\prime}){}^{t}\boldsymbol{x}\boldsymbol{x}(\boldsymbol{\beta}-\boldsymbol{\beta}^{\prime})}
\end{aligned}

が成り立つ。 \blacksquare)
 Z_1,\cdots,Z_nがそれぞれf(z_1|\boldsymbol{x}_1),\cdots,f(z_n|\boldsymbol{x}_n)から発生したと仮定するとき、Z_i=z_i,i=1,2,\cdots,nに関する対数尤度の平均は、Z-\boldsymbol{X}\boldsymbol{\beta}の平均が0(Z-\boldsymbol{X}\boldsymbol{\beta})^2の平均がn\sigma^2であることに注意すれば


\begin{aligned}
\ -E_Z\left[\displaystyle{\sum_{i=1}^{n}\log f(z_i|\boldsymbol{x}_i,\boldsymbol{\beta}^{\prime})}\right]=\displaystyle{\frac{n}{2}\log(2\pi\sigma^2e)}+\displaystyle{\frac{1}{2\sigma^2}\left(\boldsymbol{X}(\boldsymbol{\beta}-\boldsymbol{\beta}^{\prime})\right)^2}
\end{aligned}

を得る。この値は


\begin{aligned}
\displaystyle{-\sum_{i=1}^{n}\int_{-\infty}^{\infty}f(z|\boldsymbol{X}_i,\boldsymbol{\beta})\log f(z|\boldsymbol{X}_i,\boldsymbol{\beta}^{\prime})}dz
\end{aligned}

と書け、これは\boldsymbol{\beta}^{\prime}に依らない定数であるから、\mathrm{Kullback}-\mathrm{Leibler}情報量の和


\begin{aligned}
E_Z\left[\displaystyle{\sum_{i=1}^{n}\log\frac{f(z_i|\boldsymbol{x}_i,\boldsymbol{\beta})}{f(z_i|\boldsymbol{x}_i,\boldsymbol{\beta}^{\prime})}}\right]=
\displaystyle{\sum_{i=1}^{n}\int_{-\infty}^{\infty}f(z|\boldsymbol{x}_i,\boldsymbol{\beta})\log\frac{f(z|\boldsymbol{x}_i,\boldsymbol{\beta})}{f(z|\boldsymbol{x}_i,\boldsymbol{\beta}^{\prime})}}dz=\displaystyle{\frac{1}{2\sigma^2}\|\boldsymbol{X}(\boldsymbol{\beta}^{\prime}-\boldsymbol{\beta})\|^2}
\end{aligned}

を最小するように\boldsymbol{\beta}^{\prime}を選べばよい。

4.4 AICの導出

 ここでは赤池情報量基準\mathrm{AIC}を導出する。
 真の母数(ベクトル)を\boldsymbol{\beta}とおき、その不偏推定量の中で\mathrm{Kullback}-\mathrm{Leibler}情報量の和


\begin{aligned}
E_{Z}\left[\displaystyle{\sum_{i=1}^{n}\log\left(\frac{f(z|x_i,\beta)}{f(z|x_i,\hat{\beta})}\right)}\right]=
\displaystyle{\sum_{i=1}^{n}\int_{-\infty}^{\infty}f(z|x_i,\beta)\log\left(\frac{f(z|x_i,\beta)}{f(z|x_i,\hat{\beta})}\right)dz}
\end{aligned}

を最小にするような母数の推定量\hat{\beta}を選ぶ。
 確率変数U,V\in\mathbb{R}^nに対して\mathrm{Schwarz}の不等式から


\begin{aligned}
\left(E[|{}^{t}UV|]\right)^2\leq E\left\|U\right\|^2E\left\|V\right\|^2
\end{aligned}

が成り立つ。\mathrm{Cramér}-\mathrm{Rao}の不等式から任意の\betaの推定量\tilde{\beta}について対数尤度をlとおけば、実際、t\in\mathbb{R}に関する2次方程式


\begin{aligned}
E\left[(tU+V)^2\right]=t^2E[\|U\|^2]+2t E[{}^{t}UV]+E[\|V\|^2]=0
\end{aligned}

において解は存在しても高々1つであるから、判別式を考えれば、U=\boldsymbol{X}({}^{t}\boldsymbol{X}\boldsymbol{X})^{-1}\nabla l, V=\boldsymbol{X}(\tilde{\boldsymbol{\beta}}-\boldsymbol{\beta})とおけば


\begin{aligned}
\{E[{}^{t}(\tilde{\boldsymbol{\beta}}-\boldsymbol{\beta})\nabla l]\}^2\leq 
E\left[\|\boldsymbol{X}({}^{t}\boldsymbol{X}\boldsymbol{X})^{-1}\nabla l\|^2\right]
E\left[(\boldsymbol{X}(\tilde{\boldsymbol{\beta}}-\boldsymbol{\beta}))^2\right]
\end{aligned}

が得られる。
 E[(\tilde{\boldsymbol{\beta}}-\boldsymbol{\beta}){}^{t}(\nabla l)]=Iにおいて


\begin{aligned}
\mathrm{tr}\left\{E\left[(\tilde{\boldsymbol{\beta}}-\boldsymbol{\beta}){}^{t}\left(\nabla l\right)\right]\right\}=
\mathrm{tr}\left\{E\left[{}^{t}\left(\nabla l\right)(\tilde{\boldsymbol{\beta}}-\boldsymbol{\beta})\right]\right\}=
E[{}^{t}(\tilde{\boldsymbol{\beta}}-\boldsymbol{\beta})(\nabla l)]
\end{aligned}

であり、またこの右辺のトレースを取ればp+1であるから


\begin{aligned}
E[{}^{t}(\tilde{\boldsymbol{\beta}}-\boldsymbol{\beta})(\nabla l)]=p+1
\end{aligned}

が成り立つ。さらに


\begin{aligned}
E\left\|\boldsymbol{X}({}^{t}\boldsymbol{X}\boldsymbol{X})^{-1}\nabla l \right\|^2&=E\left[\mathrm{tr}\left({}^{t}(\nabla l)({}^{t}\boldsymbol{X}\boldsymbol{X})^{-1}{}^{t}\boldsymbol{X}\boldsymbol{X}({}^{t}\boldsymbol{X}\boldsymbol{X})^{-1}\nabla l\right)\right]\\
&=\mathrm{tr}\left( ({}^{t}\boldsymbol{X}\boldsymbol{X})^{-1}E\left[\nabla l\right]{}^{t}\boldsymbol{X}\boldsymbol{X}\right)\\
&=\mathrm{tr}\left( ({}^{t}\boldsymbol{X}\boldsymbol{X})^{-1}\sigma^{-2}{}^{t}\boldsymbol{X}\boldsymbol{X}\right)\\
&=\mathrm{tr}\left(\sigma^{-2}I\right)\\
&=\displaystyle{\frac{p+1}{\sigma^2}}
\end{aligned}

である。以上から


\begin{aligned}
E\left[\|\boldsymbol{X}(\tilde{\boldsymbol{\beta}}-\boldsymbol{\beta})\|^2\right]&=E\left[\mathrm{tr}{}^{t}(\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}){}^{t}\boldsymbol{X}\boldsymbol{X}(\hat{\boldsymbol{\beta}}-\boldsymbol{\beta})\right]\\
&=\mathrm{tr}\left(\mathbb{V}\left[\hat{\boldsymbol{\beta}}\right]{}^{t}\boldsymbol{X}\boldsymbol{X}\right)\\
&=\mathrm{tr}(\sigma^2 I)\\
&=(p+1)\sigma^2
\end{aligned}

が得られ、等号が成立する。
 \mathrm{AIC}は、第2項をその最小値で置き換えた


\begin{aligned}
\displaystyle{\frac{n}{2}\log(2\pi\sigma^2)}+\displaystyle{\frac{1}{2}}p+1
\end{aligned}

を最小にすることを目的とする。特に変数選択の問題では、p個すべての変数を用いるのではなく、0\leq k\leq p個の変数を選択する。すなわち


\begin{aligned}
n\log \sigma_k^2+k
\end{aligned}

を最小にするようなkを選択する。
 このとき\sigma_k^2=\displaystyle{\min_{k(S)=s}\sigma^2(S)}は未知である。そこで変数の部分集合S\subset\{1,\cdots,p\}に対して\hat{\sigma}^2(S)で代用する方法が考えられる。しかし\log\hat{\sigma}^2(S)\log\sigma^2(S)と比較して平均的に小さくなる。


標本誤差の平均値の評価 k(S)を集合Sに含まれる元の数だとして

\begin{aligned}
E\left[\log\hat{\sigma}^2(S)\right]=\log\sigma^2(S)-\displaystyle{\frac{k(S)+2}{n}}+O\left(\displaystyle{\frac{1}{n^2}}\right)
\end{aligned}

が成り立つ。

(\because m\geq1,U\sim\chi^2_m,V_1,\cdots,V_m\sim N(0,1),i.i.d.として、i=1,2,\cdots,mについて

\begin{aligned}
E\left[e^{t V_i^2}\right]&=\displaystyle{\int_{-\infty}^{\infty}e^{t v_i^2}\frac{1}{\sqrt{2\pi}}e^{-\frac{v_i^2}{2}}}dv_i=\displaystyle{\int_{-\infty}^{\infty}\frac{1}{\sqrt{2\pi}}\exp\left\{-\displaystyle{\frac{(1-2t)v_i^2}{2}}\right\} }dv_i=(1-2t)^{-\frac{1}{2}},\\
E\left[e^{t U}\right]&=\displaystyle{\int_{-\infty}^{\infty}\cdots\int_{-\infty}^{\infty}e^{t(v_1^2+\cdots+v_m^2)}\frac{1}{\sqrt{2\pi}}e^{-\frac{v_1^2+\cdots+v_m^2}{2}} }dv_1\cdots dv_m=(1-2t)^{-\frac{m}{2}}
\end{aligned}

が成り立つ。またE\left[e^{tU}\right]=1+t\cdot E[U]+\displaystyle{\frac{t^2}{2}E[U^2]}+\cdotsであるから、各n=1,2,\cdotsにおいて


\begin{aligned}
E\left[U^n\right]=\left.\displaystyle{\frac{d E[e^{tU}]}{dt^n}}\right|_{t=0}=m(m+2)\cdot\cdots\cdot(m+2n-2)
\end{aligned}

である。さらに\mathrm{Taylor}の定理から


\begin{aligned}
E\left[\log\displaystyle{\frac{U}{m}}\right]=E\left[\displaystyle{\frac{U}{m}}-1\right]-\displaystyle{\frac{1}{2}}E\left[\left(\displaystyle{\frac{U}{m}}-1\right)^2\right]+\cdots
\end{aligned}

が成立する。ここでE[U]=m,E[U^2]=m(m+2)より


\begin{aligned}
E\left[\log\displaystyle{\frac{U}{m}}\right]&=E\left[\displaystyle{\frac{U}{m}}-1\right]-\displaystyle{\frac{1}{2}}E\left[\left(\displaystyle{\frac{U}{m}}-1\right)^2\right]+\cdots\\
&=-\displaystyle{\frac{m^2}{2}}\left(E\left[U^2\right]-2m E\left[U\right]+m^2\right)+\cdots\\
&=-\displaystyle{\frac{m^2}{2}}\left(m(m+2)-2m^2+m^2\right)+\cdots\\
&=-\displaystyle{\frac{1}{m}}
\end{aligned}

が成り立つ。
 この式の右辺の第n\geq 3項について二項定理から


\begin{aligned}
E\left[(U-m)^n\right]&=\displaystyle{\sum_{i=1}^{n}{}_{n}C_{j}E[U^{j}]}(-m)^{n-j},\\
&=\displaystyle{\sum_{j=0}^{n}(-1)^{n-j}{}_{n}C_{j}m^{n-j}m(m+2)\cdot\cdots\cdot(m+2j-2)}
\end{aligned}

が得られる。
 まずm^{n-j}m(m+2)\cdot\cdots\cdot(m+2j-2)mに関する多項式と見れば、最高次の項の係数は1n-1次の係数は2\{1+2+\cdots+(j-1)\}=j(j-1)であるから、E\left[(U-m)^n\right]n次およびn-1次の項の係数は


\begin{aligned}
\displaystyle{\sum_{i=0}^{n}{}_{n}C_{j}(-1)^{n-j}}&=\displaystyle{\sum_{i=0}^{n}{}_{n}C_{j}(-1)^{n-j}1^{j}}=0,\\
\displaystyle{\sum_{i=0}^{n}{}_{n}C_{j}(-1)^{n-j}j(j-1)}&=\displaystyle{\sum_{j=2}^{n}\frac{n!}{(n-j)!(j-2)!}}(-1)^{n-j-2}\\
&=n(n-1)\displaystyle{\sum_{i=0}^{n-2}{}_{n-2}C_{i}(-1)^{n-2-i}1^{i}}=0
\end{aligned}

を得る。したがって


\begin{aligned}
E\left[(U-m)^n\right]=O\left(\displaystyle{\frac{1}{m^2}}\right)
\end{aligned}

が成り立つ。
 また\displaystyle{\frac{RSS(S)}{\sigma^2(S)}}=\displaystyle{\frac{n\hat{\sigma}^2(S)}{\sigma^2(S)}}\sim\chi_{n-k(S)-1}^2および\mathrm{Taylor}の定理の適用結果から、m=n-k(S)-1を適用して


\begin{aligned}
\displaystyle{\log\frac{n}{n-k(S)-1}}&=\displaystyle{\frac{k(S)+1}{n-k(S)-1}}+O\left(\left(\displaystyle{\frac{1}{n-k(S)-1}}\right)^2\right)\\
&=\displaystyle{\frac{k(S)+1}{n}}+O\left(\displaystyle{\frac{1}{n^2}}\right)
\end{aligned}
であり、また

\begin{aligned}
E\left[\log\left(\displaystyle{\frac{\displaystyle{\frac{\hat{\sigma^2(S)}}{n-k(S)-1}}}{\displaystyle{\frac{\sigma^2}{n}}}}\right)\right]&=\displaystyle{\frac{1}{n-k(S)-1}}+O\left(\displaystyle{\frac{1}{n^2}}\right)\\
&=-\displaystyle{\frac{1}{n}}+O\left(\displaystyle{\frac{1}{n^2}}\right)
\end{aligned}

であるから、


\begin{aligned}
E\left[\log\left(\displaystyle{\frac{\hat{\sigma}^2(S)}{\sigma^2}}\right)\right]=-\displaystyle{\frac{1}{n}}-\displaystyle{\frac{k(S)+1}{n}}+O\left(\displaystyle{\frac{1}{n^2}}\right)=\displaystyle{-\frac{k(S)+2}{n}}+O\left(\displaystyle{\frac{1}{n^2}}\right)
\end{aligned}

が得られる。 \blacksquare)

 以上から、O\left(\displaystyle{\frac{1}{n^2}}\right)の項を除いて


\begin{aligned}
E\left[\log\hat{\sigma}_k^2+\displaystyle{\frac{k+2}{n}}\right]\approx\log\sigma_k^2
\end{aligned}

が成り立つ。そこで\log\sigma_k^2\log\hat{\sigma}_k^2+\displaystyle{\frac{k+2}{n}}に置き換えることで


\begin{aligned}
n\log\hat{\sigma}_k^2+2k
\end{aligned}

を得、これを最小にするようなkを選ぶのが\mathrm{AIC}である。

4.5 実証

4.5.1 Rスクリプト

 さまざまなパッケージ内に\mathrm{AIC}を求める関数が存在する。今回は参考書P.81を基にしたものと、\mathrm{lm}(\cdot)に内蔵された\mathrm{AIC}(\cdot)とで計算してみる*3

# 残差変動を最小にするような変数の組み合わせを取得する

RSS_min <- function(mt_x,vc_y,mt_t){
  m <- ncol(mt_t)
  RSS_min <- Inf # 初期値
  
  for(j in 1:m){
    q <- mt_t[,j]
    
    # 回帰
    ls_reg <- lm(vc_y~mt_x[,q])
    RSS <- sum((ls_reg$fitted.values -vc_y)^2)/n
    
    # もし残差変動がより小さいならばその値を控えておく
    if(RSS < RSS_min){
      RSS_min <- RSS
      set_q <- q
      set_AIC <- AIC(ls_reg)
    }
  }

  return(list(value = RSS_min, set = set_q, AIC = set_AIC))
}

################################
### Bostonデータを用いた実証 ###
################################

library(MASS)

df <- Boston # テストデータ
mt_x <- as.matrix(df[,c(1,3,5:8,10:13)])
vc_y <- df[[14]]

p <- ncol(mt_x)
n <- length(vc_y)
AIC_min <- Inf

df_AIC_results <- data.frame() # 各結果(AICの値)を保存
ls_variables <- list() # 各結果(変数)

for(k in 1:p){
  mt_t <- combn(1:p, k)
  res <- RSS_min(mt_x, vc_y, mt_t)
  
  nm_AIC <- n * log(res$value / n) + 2 * k # AIC
  if(nm_AIC < AIC_min){
    AIC_min <- nm_AIC
    set_min <- res$set
  }
  
  # 各結果を保存しておく
  df_AIC_results <- rbind(df_AIC_results,c(k,nm_AIC,res$AIC)) 
  ls_variables[[k]] <- res$set
}

print(AIC_min) # 最小のときのAIC
print(set_min) # 最小のときの説明変数
4.5.2 Pythonスクリプト
###################################################
### データセットからAICを用いたモデル選択を行う ###
###################################################
import numpy as np
from sklearn.linear_model import LinearRegression
import itertools

res = LinearRegression()

### 残差変動が最も小さくなる説明変数の組み合わせを計算する
def RSS_min(X, y, T):
    S_min = np.inf # 初期値
    m = len(T)
    
    for j in range(m):
        q = T[j]
        res.fit(X[:,q], y) # 回帰
        y_hat = res.predict(X[:,q]) # 予測値
        S = np.linalg.norm(y_hat - y) ** 2 # 残差変動
        if S < S_min:
            S_min = S
            set_q = q
    return(S_min,set_q)

### 実証
# テストデータを取得
from sklearn.datasets import load_boston

boston = load_boston()
X = boston.data[:,[0,2,4,5,6,7,9,10,11,12]]
y = boston.target

n, p = X.shape # Xの行数と列数
AIC_min = np.inf # 初期値

for k in range(1,p+1,1):
    T = list(itertools.combinations(range(p), k)) # p個からk個を選ぶ組み合わせ
    S_min, set_q = RSS_min(X, y, T)
    AIC = n * np.log(S_min/n)+2 * k # AIC
    if AIC < AIC_min:
        AIC_min = AIC
        set_min = set_q
print(AIC_min, set_min)

*1:具体的には説明変数の数。

*2:これは説明上の都合でそうしたのであり、情報量基準は一般のモデルに対して適用できる。

*3:前述したように\mathrm{AIC}は真のモデルからの距離に近いもの(Kullback-Leibler情報量)として測っている。が、真のモデルの位置合いが不明なので相対値でしかない。そのため、計算パッケージによっては算出ロジックが相違する可能性がある。そのため絶対的な値を異なるパッケージの出力を比較しても意味が無い。同じパッケージ(計算ロジック)で計算したものを比較すること。実際、以下のスクリプトでは、自作した\mathrm{AIC}の計算結果とdf_AIC_resultsの3列目に保存した値は値が異なる。しかしこれらで選択したモデルは同一である。

プライバシーポリシー お問い合わせ