#はじめに
本記事は'The elements of statistical learning'の輪読ゼミ用資料です
#参考文献
An Introduction to Statistical Learning (P203~)
#発表箇所
-3.2.4 Multiple Outputs (P.56-)
3.3 Subset Selection (P.57-)
-3.3.1 Best-Subset Selection (P.57-)
-3.3.2 Forward- and Backward-Stepwise Selection (P.58-)
-3.3.3 Forward-Stagewise Regression (P.60-)
-3.3.4 Prostate Cancer Data Example(Continued) (P.61-)
##Multiple Outputs
これまでは多変数のinputに対して一つのレスポンスが返ってくるケースのみを扱ってきた。
この節では多変数のoutputが返ってくるモデルについて扱う。
\begin{align}
Y_k &= \beta_{0k} + \sum_{j=1}^{p}X_j\beta_{jk} +\varepsilon_k \quad(k=1,2,\cdots,K)\tag{3.34}\\
&=f_k(X) + \varepsilon_k \tag{3.35}
\end{align}
この時観測されたデータの数を$N$とすると
$$\rm{Y=XB+E}\tag{3.36}$$
と表記できる。
###RSS(B)
この時の$RSS(B)$について考えたい。
\begin{align}
RSS(B) &= \sum_{k=1}^{K}\sum_{i=1}^{N}(y_{ik}-f_k(x_i))^2 \tag{3.37} \\
&=\rm tr[(Y-XB)^T(Y-XB)]\tag{3.38}
\end{align}
※別解（間違っているような気がするので飛ばしてください。）
\begin{align}
RSS(B) &= \sum_{k=1}^{K}\sum_{i=1}^{N}(y_{ik}-f_k(x_i))^2 \tag{3.37} \\
&=\sum_{k=1}^{K}\sum_{i=1}^{N}\varepsilon_{ik}^2 \\
\end{align}
ここで両辺の期待値を取る。
\begin{align}
E[RSS(B)] &= E\left[\sum_{k=1}^{K}\sum_{i=1}^{N}\varepsilon_{k}^2\right] \\
&=\sum_{k=1}^{K}\sum_{i=1}^{N}E[\varepsilon_{ik}^2] \\
\end{align}
ここで
Cov(\varepsilon_i,\varepsilon_j)=\left\{
\begin{array}\\
0 & (i\not=j) \\
\sigma_i^2&(i=j)
\end{array}
\right.
より
\begin{align}
RSS(B) &= E\left[\rm{tr[E^TE]}\right] \\
&=E\left[\rm{tr}[(Y-XB)^T(Y-XB)]\right] \\
&=\rm tr\left[(Y-XB)^T(Y-XB)\right]\\
\end{align}
###推定値
回帰係数$B$と目的変数$Y$について
B=(\beta_1,\beta_2,\cdots,\beta_K) \\
Y=(y_1,y_2,\cdots,y_K)
とすると
$$\hat{\beta_i} = (X^TX)^{-1}X^Ty_i$$
となるため
$$\hat{B} = (X^TX)^{-1}X^TY \tag{3.39}$$
###相関がある場合のRSS
今までは暗黙の内に異なる要素での誤差項の相関はないものとして考えられていた。
ここでは誤差項同士に相関があった場合のケースを考えてみる事にする。
最終的な形は(3.40)式になる。
$$RSS(B;\Sigma)=\sum_{i=1}^N(y_i-f(x_i))^T\Sigma^{-1}(y_i-f(x_i)) \tag{3.40}$$
まずは、$N=1$の時のMulti Outputsに関して考える。
ここで
$$Y = BX + \varepsilon$$
それぞれ$Y\in\mathbb{R}^{K\times 1}$,$X\in\mathbb{R}^{(p+1)\times 1}$,$B\in\mathbb{R}^{K\times (p+1)}$,$\varepsilon\in \mathbb{R}^{K\times 1}$である。
また$\varepsilon$は確率変数であり標準正規分布に従う
$$\varepsilon \sim N(0,I_K)$$
\begin{align}
RSS(B)&=\sum_{j=1}^{K}\left(y_j-f_j(x)\right)^2 \\
&=\sum_{j=1}^{K}\varepsilon_j^2 \\
&=\varepsilon \varepsilon^T \tag{*}
\end{align}
ここで求めたい$\varepsilon'\sim N(0,\Sigma)$について考える。
$$\varepsilon' = \varepsilon A+\mu$$
この時
$$E(\varepsilon')=\mu$$
\begin{align}
Var(\varepsilon')&=E\left[(\varepsilon A -\mu-E[\varepsilon A -\mu])^T(\varepsilon A -\mu-E[\varepsilon A -\mu])\right] \\
&=E[A^T \varepsilon^T \varepsilon A] \\
&=A^TE[\varepsilon^T \varepsilon] A \\
&=A^TA \\
&=\Sigma
\end{align}
ここで$\mu=0$となる事がわかり
$$\varepsilon = \varepsilon'A^{-1}$$
であり$(*)$に代入すれば
\begin{align}
RSS(B)&=\varepsilon' A^{-1}(\varepsilon' A^{-1})^T \\
&=\varepsilon' A^{-1}(A^{-1})^T\varepsilon'^T \\
&=\varepsilon' (A^TA)^{-1} \varepsilon'^T \\
&=\varepsilon' \Sigma^{-1} \varepsilon'^T
\end{align}
ここで$\varepsilon' = (y-f(x))^T$であり、この$RSS$に関して$N$個の観測値を足し合わせればよいので
$$RSS(B;\Sigma)=\sum_{i=1}^N(y_i-f(x_i))^T\Sigma^{-1}(y_i-f(x_i))$$
となり(3.44)を得る事が出来た。
###Exercise 3.11
相関ありの場合でも得られるBの推定値は相関がない時と同じく$\hat{B} = (X^TX)^{-1}X^TY$である事を導く。
ここで$$f(x_i) = Bx_i$$
とする。つまりここでは$B\in \mathbb{R}^{K\times (p+1)}$と置いている（この本ではこの$B$を転置したものを$B$としている事が多い）
\begin{align}
\frac{\partial}{\partial B} RSS(B;\Sigma)&=\frac{\partial}{\partial B}\sum_{i=1}^N(y_i-Bx_i)^T\Sigma^{-1}(y_i-Bx_i)
\\
&= \sum_{i=1}^N \frac{\partial}{\partial B} \left( y_i^T\Sigma^{-1} y_i - y_i^T \Sigma^{-1}Bx_i - x_i^TB^T\Sigma^{-1}y_i + x_i^TB^T\Sigma^{-1}Bx_i \right) \\
&= \sum_{i=1}^{N}\left(0-(y_i^T \Sigma^{-1})^Tx_i^T-(x_i(\Sigma^{-1}y_i)^T)^T + \frac{\partial}{\partial B}tr\left[B^T\Sigma^{-1} B\right]\frac{\partial}{\partial U}x_i^TUx_i\right) \\
&=\sum_{i=1}^{N}\left(-2\Sigma^{-1}y_ix_i^T+\left(B^T\left(\Sigma^{-1}+(\Sigma^{-1})^T\right)\right)^Tx_ix_i^T\right) \\
&=\sum_{i=1}^{N}\left( -2\Sigma^{-1}y_ix_i^T + 2\Sigma^{-1}Bx_ix_i^T \right) \\
&=\sum_{i=1}^{N}\left( 2\Sigma^{-1}(Bx_i-Y_i)x_i^T \right)
\end{align}
上の式を解く過程の中で以下の変形を用いた
\frac{\partial x^TBy}{\partial B} = xy^T \\
\frac{\partial f\left(U=g(A)\right)}{\partial A} = \frac{\partial tr(U)}{\partial A}\frac{\partial x^TUx}{\partial U}
(自信なし)\\
\frac{\partial tr(B^TAB)}{\partial B} = B(A^T+A) \\
\frac{\partial tr(B^TAB)}{\partial B^T} = \left( \frac{\partial tr(B^TAB)}{\partial B} \right)^T \\
x\in\mathbb{R}^{n\times 1},A\in\mathbb{R}^{n\times n},B\in\mathbb{R}^{n\times m}
$\frac{\partial RSS}{\partial B}=0$となる$B$が求めたい$\hat{B}$なので
\begin{align}
& \sum_{i=1}^{N}\left( 2\Sigma^{-1}(Bx_i-Y_i)x_i^T \right) =0 \\
&\Leftrightarrow \sum_{i=1}^{N}\left((Bx_i-Y_i)x_i^T \right) =0 \\
&\Leftrightarrow \left(BX-Y\right)X^T = 0 \\
&\Leftrightarrow BXX^T = YX^T \\
&\Leftrightarrow B=YX^T(XX^T)^{-1}\quad(*) \\
\end{align}
(3.39)式で書かれている$B,X,Y$は$(*)$で示された$B,X,Y$を転置した形なので(3.39)式に合わせて書き直せば
\hat{B}=\left( Y^TX(X^TX)^{-1} \right) ^T \Leftrightarrow \hat{B}=(X^TX)^{-1}X^TY
よって題意は示された。
#Subset Selection（部分集合選択）
最小二乗法によって回帰係数を決定しようとした時に低バイアス高バリアンス（過学習）の状況に陥る事がある。
例としてサイン波の回帰問題を考える。
$$y=x^0\beta_0+x^1\beta_1+\cdots+x^M\beta_M$$
import numpy as np
import matplotlib.pyplot as plt
# y = b0*x^0+....bM*x^M を、引数xの配列数分求める
def y(w, x, M):
X = np.empty((M + 1, x.size))
for i in range(M + 1):
X[i, :] = x ** i
return np.dot(w.T, X)
# ランダムシードを固定
np.random.seed(0)
# 多項式の最大べき乗数(x^0+...+x^M)
M = 3
# 訓練データ数
N = 11
# 訓練データの列ベクトル
x = np.linspace(0, 1, N).reshape(N, 1)
# 訓練データtの列ベクトル
t = np.sin(2*np.pi*x.T) + np.random.normal(0, 0.2, N)
t = t.reshape(N, 1)
# 行列Xを作成
X = np.empty((N, M+1))
for i in range(M+1):
X[:, i] = x.reshape(1, N) ** i
# 係数wの列ベクトルを求める
w = np.linalg.solve(np.dot(X.T, X), np.dot(X.T, t))
# 求めた係数wを元に、新たな入力x2に対する予測値yを求める
x2 = np.linspace(0, 2, 100)
y1 = y(w, x2, M)
# 結果の表示
plt.xlim(0, 1.3)
plt.ylim(-1.5, 1.5)
plt.scatter(x, t)
plt.plot(x2, y1.T)
plt.title(f'M={M}')
plt.show()
よって高次元の回帰問題には過学習の恐れがある為、Outputに対して影響力が高い変数のみを抽出し、それらの変数のみで回帰問題を解いた方がバリアンスを抑える可能性が高く、汎用性能が高くなる事が予想される。
以下ではSubsetの効率的な選び方として知られてるいくつかの手法を紹介する。
##Best-subset Selection
もっとも基本的で初歩的な手法としてBest-subset Seletionが知られる。
###Algorithm
$(i)$すべての$k \in [0,1,2,\cdots,p] $がとりえるすべての変数の組み合わせについて$RSS(\beta)$を求める。
###メリット
- 取り得るすべての$k$で必ず最適な$RSS$を求める事が出来る。
###デメリット
- 組み合わせの試行回数が$2^p$回あり、変数の数が膨大であるとコンピュータの計算処理能力を考えると現実的な手法ではなくなる。
- 適切な計算処理を行えば$p$が30から40までが利用可能であるとされている。$(2^{40}=1.10\times10^{12})$
###ex)プロ野球選手のホームラン数(Best-subset Selection)
> library(ISLR)
> fix(Hitters)
> names(Hitters)
[1] "AtBat" "Hits" "HmRun" "Runs" "RBI" "Walks"
[7] "Years" "CAtBat" "CHits" "CHmRun" "CRuns" "CRBI"
[13] "CWalks" "League" "Division" "PutOuts" "Assists" "Errors"
[19] "Salary" "NewLeague"
> library(leaps)
> regfit.full=regsubsets(Salary~.,Hitters)
> summary(regfit.full)
1 subsets of each size up to 8
Selection Algorithm: exhaustive
AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun CRuns
1 ( 1 ) " " " " " " " " " " " " " " " " " " " " " "
2 ( 1 ) " " "*" " " " " " " " " " " " " " " " " " "
3 ( 1 ) " " "*" " " " " " " " " " " " " " " " " " "
4 ( 1 ) " " "*" " " " " " " " " " " " " " " " " " "
5 ( 1 ) "*" "*" " " " " " " " " " " " " " " " " " "
6 ( 1 ) "*" "*" " " " " " " "*" " " " " " " " " " "
7 ( 1 ) " " "*" " " " " " " "*" " " "*" "*" "*" " "
8 ( 1 ) "*" "*" " " " " " " "*" " " " " " " "*" "*"
CRBI CWalks LeagueN DivisionW PutOuts Assists Errors NewLeagueN
1 ( 1 ) "*" " " " " " " " " " " " " " "
2 ( 1 ) "*" " " " " " " " " " " " " " "
3 ( 1 ) "*" " " " " " " "*" " " " " " "
4 ( 1 ) "*" " " " " "*" "*" " " " " " "
5 ( 1 ) "*" " " " " "*" "*" " " " " " "
6 ( 1 ) "*" " " " " "*" "*" " " " " " "
7 ( 1 ) " " " " " " "*" "*" " " " " " "
8 ( 1 ) " " "*" " " "*" "*" " " " " " "
6回まではCRBI(生涯出場打席数)が変数として選ばれていたものの7回以降はCRBIは選ばれない
##Forward-Stepwise Selection
###Algorithm
$(i)$全ての変数の係数を0として切片の係数のみで回帰する
$(ii)RSS$が最も低くなるような変数を一つ選び加える
$(iii)$全ての変数が加わるまで$(ii)$を繰り返す
###メリット
- 試行回数が$1+\frac{p(p+1)}{2}$なのでBest-subset Secectionに比べると圧倒的に計算が楽であり次元数が膨大になっても計算可能である。
- バリアンスが低くなりやすい
###デメリット
- 変数を一つ加える際には今まで選んだ変数を変える事はない為、各$k$に対して最適な$RSS$になっているとは限らない(sub-optimal)
###ex)プロ野球選手のホームラン数(Best-subset Selection)
> regfit.fwd=regsubsets(Salary~.,data=Hitters , method="forward")
> summary(regfit.fwd)
1 subsets of each size up to 8
Selection Algorithm: forward
AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun CRuns
1 ( 1 ) " " " " " " " " " " " " " " " " " " " " " "
2 ( 1 ) " " "*" " " " " " " " " " " " " " " " " " "
3 ( 1 ) " " "*" " " " " " " " " " " " " " " " " " "
4 ( 1 ) " " "*" " " " " " " " " " " " " " " " " " "
5 ( 1 ) "*" "*" " " " " " " " " " " " " " " " " " "
6 ( 1 ) "*" "*" " " " " " " "*" " " " " " " " " " "
7 ( 1 ) "*" "*" " " " " " " "*" " " " " " " " " " "
8 ( 1 ) "*" "*" " " " " " " "*" " " " " " " " " "*"
CRBI CWalks LeagueN DivisionW PutOuts Assists Errors NewLeagueN
1 ( 1 ) "*" " " " " " " " " " " " " " "
2 ( 1 ) "*" " " " " " " " " " " " " " "
3 ( 1 ) "*" " " " " " " "*" " " " " " "
4 ( 1 ) "*" " " " " "*" "*" " " " " " "
5 ( 1 ) "*" " " " " "*" "*" " " " " " "
6 ( 1 ) "*" " " " " "*" "*" " " " " " "
7 ( 1 ) "*" "*" " " "*" "*" " " " " " "
8 ( 1 ) "*" "*" " " "*" "*" " " " " " "
##Exercise 3.9（未完成）
$X=QR$と分解すると(詳細はP.55)
\begin{align}
\hat{\beta} &= (X^TX)^{-1}X^Ty \\
&= (R^TQ^TQR)^{-1}R^TQ^Ty \\
&= R^{-1}(R^T)^{-1}R^TQ^Ty \\
&= R^{-1}Q^Ty\\
\\
\hat{y} &= X\hat{\beta} \\
&= QRR^{-1}Q^Ty \\
&=QQ^Ty
\end{align}
より$$R\hat{\beta}=Q^Ty$$
という関係が導く事ができる。ここで$R$は上三角行列なので後退法により簡単に(逆行列の計算をなしに)$\hat{\beta}$を求める事が可能になる。
ここで$X_1$を$N\times q$行列とし、$X_1$で選ばれた変数に加えてもう一つ新たに変数を加えた行列を$X_2$($\in \mathbf{R}^{N\times (q+1)}$)とする。この時$X_2$は
X_2=
\left[
\begin{matrix}
X_1&x_q
\end{matrix}
\right]\quad(x_q\in \mathbf{R}^{N \times 1})
であり
\begin{align}
RSS(\beta) &= \left(y-X_2(X_2^TX_2)^{-1}X_2^Ty\right)^T\left(y-X_2(X_2^TX_2)^{-1}X_2^Ty\right)
\end{align}
\begin{align}
X_2(X_2^TX_2)^{-1}X_2^Ty &=
\begin{matrix}
(X_1 & x_q)
\end{matrix}
\left[
\left(
\begin{matrix}
X_1^T \\
x_q^T
\end{matrix}
\right)
\left(
\begin{matrix}
X_1 & x_q
\end{matrix}
\right)
\right]^{-1}
\left(
\begin{matrix}
X_1^T \\
x_q^T
\end{matrix}
\right)
y
\\
&=
\begin{matrix}
(X_1 & x_q)
\end{matrix}
\left[
\begin{matrix}
X_1^TX_1 & X_1^Tx_q \\
x_q^TX_1 & x_q^Tx_q
\end{matrix}
\right]^{-1}
\left(
\begin{matrix}
X_1^T \\
x_q^T
\end{matrix}
\right)
y \\
&=\cdots
\end{align}
\begin{align}
\left(
\begin{matrix}
A & B \\
C & D
\end{matrix}
\right)^{-1}
=
\left(
\begin{matrix}
A-BD^{-1}C)^{-1} & -A^{-1}BS^{-1} \\
-S^{-1} CA^{-1} & S^{-1}
\end{matrix}
\right)
\\
(S=D-CA^{-1}B)
\end{align}
後者のやり方を行えば$QR$分解を実行する回数はp回になるので（上記の方法では$1+\frac{p(p+1)}{2}$回)の為計算上楽になるような気はするがコンピュータ上でどうなるのかは不明
##Backward-Stepwise Selection
###Algorithm
$(i)$全ての変数を考慮したモデル（full model)で回帰する
$(ii)RSS$に最も寄与しない変数を一つ削る
$(iii)$全ての変数がなくなるまで$(ii)$を繰り返す
###メリット
- 試行回数が$1+\frac{p(p+1)}{2}$なのでBest-subset Secectionに比べると圧倒的に計算が楽であり次元数が膨大になっても計算可能である。
###デメリット
- 変数を一つ加える際には今まで選んだ変数を変える事はない為、各$k$に対して最適な$RSS$になっているとは限らない
- N>pのみでしか扱えない（なぜ？）
##Exercise 3.10
##Forward Stagewise Selection
###Algorithm
$(i)\bar y$に等しい切片から始める（それ以外の係数はゼロ）
$(ii)$残差と最も相関がある係数を選出する
$(iii)$相関がなくなるまで$(ii)$を繰り返す
###メリット
- 高次元のSubset Selectionで高い効果を出す。
###デメリット
- ステップ数大きく、他の手法に比べて時間がかかる
##まとめ
- 基本的な手法選びとしては
- $p$の数が少ない時：Best-subset Selection
- $p$の数が多い時：Forward stepwise Selection
- $p$がとても多い時：Forward stagewise Selection
- Backward Selectionは使い道が分からない。(Nが多い時？）
- 今日述べたすべてのsubset selection に言える事であるが、扱う変数の数を増やせば増やす程$RSS$は低くなってしまう為、これを指標に適切な変数$k$を特定する事はできない。
- バリアンスが低いモデル（過学習していないモデル）を探すには、データをテストデータと訓練データに分け、テストデータに対してcross-validationなどを用いた誤差測定を行う事で見つけ出す事ができる（詳細は7章）