「統計学実践ワークブック」 演習問題etc Ch.30 「モデル選択」

当記事は「統計学実践ワークブック(学術図書出版社)」の読解サポートにあたってChapter.30の「モデル選択」に関して演習問題を中心に解説を行います。AIC(Akaike information criterion)やBIC(Bayesian information criterion)はよく出てくるトピックなので、抑えておくと良いと思います。

日本統計学会公式認定 統計検定準1級対応 統計学実践ワークブック

3,080円(04/20 02:24時点)

Amazon

本章のまとめ

本章のまとめ

AIC

BIC

演習問題解説

問30.1

$$
\large
\begin{align}
L &= \max_{\mathbf{\beta},\sigma^2} (2 \pi)^{-\frac{n}{2}} |\sigma^2 I_{n}|^{-\frac{1}{2}} \exp \left( -\frac{1}{2}(\mathbf{y}-X\mathbf{\beta})^{\mathrm{T}}(\sigma^2 I_{n})^{-1}(\mathbf{y}-X\mathbf{\beta}) \right) \\
&= \max_{\mathbf{\beta},\sigma^2} (2 \pi)^{-\frac{n}{2}} (\sigma^2)^{-\frac{n}{2}} \exp \left( -\frac{1}{2 \sigma^2}(\mathbf{y}-X\mathbf{\beta})^{\mathrm{T}}(\mathbf{y}-X\mathbf{\beta}) \right) \quad (1)
\end{align}
$$

上記のように表した$(1)$式に対して、$\mathbf{\beta},\sigma^2$の最尤推定量である$\displaystyle \hat{\mathbf{\beta}}=(X^{\mathrm{T}}X)^{-1}X^{\mathrm{T}}\mathbf{y},\hat{\sigma}^2=\frac{S_e}{n}$を代入した際に$L$は最大値をとる。よって$(1)$式は下記のように変形できる。
$$
\large
\begin{align}
L &= \max_{\mathbf{\beta},\sigma^2} (2 \pi)^{-\frac{n}{2}} (\sigma^2)^{-\frac{n}{2}} \exp \left( -\frac{1}{2 \sigma^2}(\mathbf{y}-X\mathbf{\beta})^{\mathrm{T}}(\mathbf{y}-X\mathbf{\beta}) \right) \quad (1) \\
&= (2 \pi)^{-\frac{n}{2}} \left( \frac{S_e}{n} \right)^{-\frac{n}{2}} \exp \left( -\frac{n}{2 S_e}(\mathbf{y}-X(X^{\mathrm{T}}X)^{-1}X^{\mathrm{T}}\mathbf{y})^{\mathrm{T}}(\mathbf{y}-X(X^{\mathrm{T}}X)^{-1}X^{\mathrm{T}}\mathbf{y}) \right) \\
&= \left( \frac{2 \pi S_e}{n} \right)^{-\frac{n}{2}} \exp \left( -\frac{n\cancel{S_e}}{2\cancel{S_e}}\right) = \left( \frac{2 \pi S_e}{n} \right)^{-\frac{n}{2}} \exp \left( -\frac{n}{2}\right) \quad (30.2)
\end{align}
$$

上記の変形では$S_e=(\mathbf{y}-X(X^{\mathrm{T}}X)^{-1}X^{\mathrm{T}}\mathbf{y})^{\mathrm{T}}(\mathbf{y}-X(X^{\mathrm{T}}X)^{-1}X^{\mathrm{T}}\mathbf{y})$を用いたので、以下この式が成立することを確認する。
$$
\large
\begin{align}
& (\mathbf{y}-X(X^{\mathrm{T}}X)^{-1}X^{\mathrm{T}}\mathbf{y})^{\mathrm{T}}(\mathbf{y}-X(X^{\mathrm{T}}X)^{-1}X^{\mathrm{T}}\mathbf{y}) \\
&= \mathbf{y}^{\mathrm{T}}\mathbf{y} – 2\mathbf{y}^{\mathrm{T}}X(X^{\mathrm{T}}X)^{-1}X^{\mathrm{T}}\mathbf{y} + \mathbf{y}^{\mathrm{T}}X\cancel{(X^{\mathrm{T}}X)^{-1}}\cancel{X^{\mathrm{T}}X}(X^{\mathrm{T}}X)^{-1}X^{\mathrm{T}}\mathbf{y} \\
&= \mathbf{y}^{\mathrm{T}}\mathbf{y} – 2\mathbf{y}^{\mathrm{T}}X(X^{\mathrm{T}}X)^{-1}X^{\mathrm{T}}\mathbf{y} + \mathbf{y}^{\mathrm{T}}X(X^{\mathrm{T}}X)^{-1}X^{\mathrm{T}}\mathbf{y} \\
&= \mathbf{y}^{\mathrm{T}}\mathbf{y} – \mathbf{y}^{\mathrm{T}}X(X^{\mathrm{T}}X)^{-1}X^{\mathrm{T}}\mathbf{y} \\
&= \mathbf{y}^{\mathrm{T}} (I_{n} – X(X^{\mathrm{T}}X)^{-1}X^{\mathrm{T}}) \mathbf{y} \\
&= S_e
\end{align}
$$

問30.2

下記を実行することで、説明変数の組み合わせごとのAIC、BICの値が得られる。

import numpy as np

n = 20.
x_combination = ["nothing", "x1", "x2", "x3", "x1-x2", "x1-x3", "x2-x3", "x1-x2-x3"]
k = np.array([0., 1., 1., 1., 2., 2., 2., 3.])
res = np.array([[53180., 10.881], [41499., 10.633], [30338., 10.320], [42420., 10.607], [29727., 10.3], [26252., 10.176], [26923., 10.201], [23893., 10.081]])

AIC = n*(res[:,1]+np.log(2*np.pi/n)+1)+2*(k+2)
BIC = n*(res[:,1]+np.log(2*np.pi/n)+1)+(k+2)*np.log(n)

for i in range(res.shape[0]):
    print("Variable: {}, AIC: {:.3f}, BIC: {:.3f}".format(x_combination[i],AIC[i],BIC[i]))

・実行結果

Variable: nothing, AIC: 218.463, BIC: 220.454
Variable: x1, AIC: 215.503, BIC: 218.490
Variable: x2, AIC: 209.243, BIC: 212.230
Variable: x3, AIC: 214.983, BIC: 217.970
Variable: x1-x2, AIC: 210.843, BIC: 214.826
Variable: x1-x3, AIC: 208.363, BIC: 212.346
Variable: x2-x3, AIC: 208.863, BIC: 212.846
Variable: x1-x2-x3, AIC: 208.463, BIC: 213.442

上記の計算結果より、AICによる変数選択で選ばれる説明変数の組み合わせは「$x_1,x_3$」、BICによる変数選択で選ばれる説明変数の組み合わせは「$x_2$」であることがそれぞれわかる。

参考

・準$1$級関連まとめ
https://www.hello-statisticians.com/toukeikentei-semi1