統計的機械学習の数理100問（13/20） - 「大人の教養・知識・気付き」を伸ばすブログ

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

統計的機械学習の数理100問 with R (with R) (機械学習の数理100問シリーズ)

作者:讓, 鈴木
共立出版

Amazon

統計的機械学習の数理100問 with Python (機械学習の数理100問シリーズ)

作者:讓, 鈴木
共立出版

Amazon

を用いることにする。

前回
4.　情報量基準
次回

前回

power-of-awareness.com

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 x$ が $x$ に関して単調増加であることを踏まえ、対数尤度

$\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$ に関する多項式と見れば、最高次の項の係数は $1$ で $n-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)

次回

power-of-awareness.com

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

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

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

前回

4. 情報量基準

4.1 情報量基準

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

4.3 Kullback-Leibler情報量

4.4 AICの導出

4.5 実証

4.5.1 Rスクリプト

4.5.2 Pythonスクリプト

次回

4.　情報量基準

4.1　情報量基準

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

4.3　Kullback-Leibler情報量

4.4　AICの導出

4.5　実証

4.5.1　Rスクリプト

4.5.2　Python スクリプト