(Qiitaに投稿していた当時)Qiitaで重回帰と検索をかけても意外と数式での説明がなかったので今回は数式で攻めたいと思います。
例題として、The Boston Housing Dataset を使います。
crim
zn
indus
chas
nox
rm
age
dis
rad
tax
ptratio
black
lstat
medv
1
0.00632
18
2.31
0
0.538
6.575
65.2
4.0900
1
296
15.3
396.90
4.98
24.0
2
0.02731
0
7.07
0
0.469
6.421
78.9
4.9671
2
242
17.8
396.90
9.14
21.6
3
0.02729
0
7.07
0
0.469
7.185
61.1
4.9671
2
242
17.8
392.83
4.03
34.7
4
0.03237
0
2.18
0
0.458
6.998
45.8
6.0622
3
222
18.7
394.63
2.94
33.4
5
0.06905
0
2.18
0
0.458
7.147
54.2
6.0622
3
222
18.7
396.90
5.33
36.2
上のデータセットの列は、それぞれ次のようなことを意味しています。
column name
列名の意味するところ
crim
犯罪率
zn
宅地の割合
indus
非商用地の割合
chas
チャールズ川流域かどうか
nox
窒素酸化物濃度
rm
平均部屋数
age
築年数
dis
ビジネス地区への距離
rad
高速道路へのアクセス指数
tax
固定資産税
ptratio
学生と教師の割合
black
黒人の割合
lstat
低所得者の割合
mdev
住宅価格の中央値
このデータを用いて、次のようなことを考えます。
住宅価格の中央値を、目的変数として、その他変数(e.g., 犯罪率)を説明変数として予測モデルを立てたい。つまり、
住宅価格中央値 = β 0 + β 1 × 犯罪率 + β 2 × 住宅の割合 + ⋯ + β 13 × 低所得者層の割合 住宅価格中央値 = \beta_0 + \beta_1 \times 犯罪率 + \beta_2 \times 住宅の割合 + \cdots + \beta_{13} \times 低所得者層の割合 住 宅 価 格 中 央 値 = β 0 + β 1 × 犯 罪 率 + β 2 × 住 宅 の 割 合 + ⋯ + β 1 3 × 低 所 得 者 層 の 割 合 と表現できるような
β 0 , β 1 , ⋯ , β 13 \beta_0, \beta_1, \cdots, \beta_{13} β 0 , β 1 , ⋯ , β 1 3
を求めたい。
もう少し一般化したものを考えてみます。以下の表の目的変数 y y y p p p x 1 ⋯ x p x_1 \cdots x_p x 1 ⋯ x p
index
y y y x 1 x_1 x 1 x 2 x_2 x 2 ⋯ \cdots ⋯ x p x_p x p
1 1 1 y 1 y_1 y 1 x 11 x_{11} x 1 1 x 12 x_{12} x 1 2 ⋯ \cdots ⋯ x 1 p x_{1p} x 1 p
2 2 2 y 2 y_2 y 2 x 21 x_{21} x 2 1 x 22 x_{22} x 2 2 ⋯ \cdots ⋯ x 2 p x_{2p} x 2 p
⋮ \vdots ⋮ ⋮ \vdots ⋮ ⋮ \vdots ⋮ ⋮ \vdots ⋮ ⋮ \vdots ⋮ ⋮ \vdots ⋮
i i i y i y_i y i x i 1 x_{i1} x i 1 x i 2 x_{i2} x i 2 ⋯ \cdots ⋯ x i p x_{ip} x i p
⋮ \vdots ⋮ ⋮ \vdots ⋮ ⋮ \vdots ⋮ ⋮ \vdots ⋮ ⋮ \vdots ⋮ ⋮ \vdots ⋮
n n n y n y_n y n x n 1 x_{n1} x n 1 x n 2 x_{n2} x n 2 ⋯ \cdots ⋯ x n p x_{np} x n p
重回帰モデル(multiple regression model) は以下のように表されます。
y = β 0 + β 1 x 1 + ⋯ + β p x p + ϵ y = \beta_0 + \beta_1x_1 + \cdots + \beta_px_p + \epsilon y = β 0 + β 1 x 1 + ⋯ + β p x p + ϵ このとき、上の表においてi i i
y i = β 0 + β 1 x i 1 + ⋯ + β p x i p + ϵ i y_i = \beta_0 + \beta_1x_{i1} + \cdots + \beta_px_{ip} + \epsilon_i y i = β 0 + β 1 x i 1 + ⋯ + β p x i p + ϵ i と表せる。
ベクトルを用いて表すと、
[ y 1 ⋮ y n ] = [ 1 x 11 ⋯ x 1 p ⋮ ⋮ ⋱ ⋮ 1 x n 1 ⋯ x n p ] [ β 0 ⋮ β p ] + [ ϵ 1 ⋮ ϵ n ] \begin{bmatrix}
y_1 \\
\vdots \\
y_n
\end{bmatrix} =
\begin{bmatrix}
1 & x_{11} & \cdots & x_{1p} \\
\vdots & \vdots & \ddots & \vdots \\
1 & x_{n1} & \cdots & x_{np}
\end{bmatrix}
\begin{bmatrix}
\beta_0 \\
\vdots \\
\beta_p
\end{bmatrix} +
\begin{bmatrix}
\epsilon_1 \\
\vdots \\
\epsilon_n
\end{bmatrix} ⎣ ⎢ ⎢ ⎡ y 1 ⋮ y n ⎦ ⎥ ⎥ ⎤ = ⎣ ⎢ ⎢ ⎡ 1 ⋮ 1 x 1 1 ⋮ x n 1 ⋯ ⋱ ⋯ x 1 p ⋮ x n p ⎦ ⎥ ⎥ ⎤ ⎣ ⎢ ⎢ ⎡ β 0 ⋮ β p ⎦ ⎥ ⎥ ⎤ + ⎣ ⎢ ⎢ ⎡ ϵ 1 ⋮ ϵ n ⎦ ⎥ ⎥ ⎤ となり、
y = [ y 1 ⋮ y n ] \boldsymbol{y}
=
\begin{bmatrix}
y_1 \\
\vdots \\
y_n
\end{bmatrix} y = ⎣ ⎢ ⎢ ⎡ y 1 ⋮ y n ⎦ ⎥ ⎥ ⎤ X = [ 1 x 11 ⋯ x 1 p ⋮ ⋮ ⋱ ⋮ 1 x n 1 ⋯ x n p ] X=
\begin{bmatrix}
1 & x_{11} & \cdots & x_{1p} \\
\vdots & \vdots & \ddots & \vdots \\
1 & x_{n1} & \cdots & x_{np}
\end{bmatrix} X = ⎣ ⎢ ⎢ ⎡ 1 ⋮ 1 x 1 1 ⋮ x n 1 ⋯ ⋱ ⋯ x 1 p ⋮ x n p ⎦ ⎥ ⎥ ⎤ β = [ β 0 ⋮ β p ] \boldsymbol{\beta}=
\begin{bmatrix}
\beta_0 \\
\vdots \\
\beta_p
\end{bmatrix} β = ⎣ ⎢ ⎢ ⎡ β 0 ⋮ β p ⎦ ⎥ ⎥ ⎤ ϵ = [ ϵ 1 ⋮ ϵ n ] \boldsymbol{\epsilon}=
\begin{bmatrix}
\epsilon_1 \\
\vdots \\
\epsilon_n
\end{bmatrix} ϵ = ⎣ ⎢ ⎢ ⎡ ϵ 1 ⋮ ϵ n ⎦ ⎥ ⎥ ⎤ とおくと、
y = X β + ϵ \boldsymbol{y} = X\boldsymbol{\beta} + \boldsymbol{\epsilon} y = X β + ϵ と表すことができます。
確認ですが、
y \boldsymbol{y} y X X X 計画行列(design matrix) β \boldsymbol{\beta} β これを求めたい。 )ϵ \boldsymbol{\epsilon} ϵ
となります。
未知パラメータであるβ \boldsymbol{\beta} β
S ( β ) = ∑ i = 1 n ϵ i 2 = ϵ T ϵ = ( y − X β ) T ( y − X β ) S(\boldsymbol{\beta}) = \sum_{i=1}^{n}{\epsilon_i^2} = \boldsymbol{\epsilon}^{\mathsf{T}}\boldsymbol{\epsilon}
= (\boldsymbol{y}-X\boldsymbol{\beta})^{\mathsf{T}}(\boldsymbol{y}-X\boldsymbol{\beta}) S ( β ) = i = 1 ∑ n ϵ i 2 = ϵ T ϵ = ( y − X β ) T ( y − X β ) を最小にするβ \boldsymbol{\beta} β
S ( β ) = ( y − X β ) T ( y − X β ) = ( y T − ( X β ) T ) ( y − X β ) = ( y T − β T X T ) ( y − X β ) = y T y − β T X T y − y T X β + β T X T X β = y T y − 2 y T X β + β T X T X β \begin{aligned}
S(\boldsymbol{\beta}) & = (\boldsymbol{y}-X\boldsymbol{\beta})^{\mathsf{T}}(\boldsymbol{y}-X\boldsymbol{\beta}) \\
& = (\boldsymbol{y}^{\mathsf{T}}-(X\boldsymbol{\beta})^{\mathsf{T}})(\boldsymbol{y}-X\boldsymbol{\beta}) \\
& = (\boldsymbol{y}^{\mathsf{T}}-\boldsymbol{\beta}^{\mathsf{T}}X^{\mathsf{T}})(\boldsymbol{y}-X\boldsymbol{\beta}) \\
& = \boldsymbol{y}^{\mathsf{T}}\boldsymbol{y} - \boldsymbol{\beta}^{\mathsf{T}}X^{\mathsf{T}}\boldsymbol{y}-\boldsymbol{y}^{\mathsf{T}}X\boldsymbol{\beta} + \boldsymbol{\beta}^{\mathsf{T}}X^{\mathsf{T}}X\boldsymbol{\beta} \\
& = \boldsymbol{y}^{\mathsf{T}}\boldsymbol{y} -2\boldsymbol{y}^{\mathsf{T}}X\boldsymbol{\beta} + \boldsymbol{\beta}^{\mathsf{T}}X^{\mathsf{T}}X\boldsymbol{\beta}
\end{aligned} S ( β ) = ( y − X β ) T ( y − X β ) = ( y T − ( X β ) T ) ( y − X β ) = ( y T − β T X T ) ( y − X β ) = y T y − β T X T y − y T X β + β T X T X β = y T y − 2 y T X β + β T X T X β これをβ \boldsymbol{\beta} β
S ( β ) ∂ β = − 2 X T y + 2 X T X β \frac{S(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}= -2X^{\mathsf{T}}\boldsymbol{y} + 2X^{\mathsf{T}}X\boldsymbol{\beta} ∂ β S ( β ) = − 2 X T y + 2 X T X β S ( β ) ∂ β = 0 \frac{S(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = 0 ∂ β S ( β ) = 0
− 2 X T y + 2 X T X β = 0 X T X β = X T y \begin{aligned}
-2X^{\mathsf{T}}\boldsymbol{y} + 2X^{\mathsf{T}}X\boldsymbol{\beta} & = 0 \\
X^{\mathsf{T}}X\boldsymbol{\beta} = X^{\mathsf{T}}\boldsymbol{y}
\end{aligned} − 2 X T y + 2 X T X β X T X β = X T y = 0 もし ( X T X ) − 1 (X^{\mathsf{T}}X)^{-1} ( X T X ) − 1 が存在するなら、
β = ( X T X ) − 1 X T y \boldsymbol{\beta} = (X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}\boldsymbol{y} β = ( X T X ) − 1 X T y と求めることができる。
参考：データ解析 Rによる多変量解析入門 (7) 最尤法、モデル選択
重回帰式
y = X β + ϵ \boldsymbol{y} = X\boldsymbol{\beta} + \boldsymbol{\epsilon} y = X β + ϵ の誤差ベクトル ϵ \boldsymbol{\epsilon} ϵ
ϵ ∼ N n ( 0 , σ 2 I n ) \boldsymbol{\epsilon} \sim N_n(\boldsymbol{0}, \sigma^2 I_n) ϵ ∼ N n ( 0 , σ 2 I n ) このとき、
y ∼ N n ( X β , σ 2 I n ) \boldsymbol{y} \sim N_n(X\boldsymbol{\beta}, \sigma^2 I_n) y ∼ N n ( X β , σ 2 I n ) となり、確率密度関数は、
f ( y ∣ X , β , σ 2 ) = 1 ( 2 π σ 2 I n ) n / 2 exp [ − 1 2 σ 2 ( y − X β ) T ( y − X β ) ] 。 f(\boldsymbol{y} | X, \boldsymbol{\beta}, \sigma^2) = \frac{1}{(2\pi\sigma^2 I_n)^{n/2}}\exp \biggl[ - \frac{1}{2\sigma^2}(\boldsymbol{y}-X\boldsymbol{\beta})^{\mathsf{T}}(\boldsymbol{y}-X\boldsymbol{\beta}) \biggr]。 f ( y ∣ X , β , σ 2 ) = ( 2 π σ 2 I n ) n / 2 1 exp [ − 2 σ 2 1 ( y − X β ) T ( y − X β ) ] 。 尤度は、
L ( β , σ 2 ∣ X , y ) = f ( y ∣ X , β , σ 2 ) L(\boldsymbol{\beta}, \sigma^2 | X, \boldsymbol{y})=f(\boldsymbol{y} | X, \boldsymbol{\beta}, \sigma^2) L ( β , σ 2 ∣ X , y ) = f ( y ∣ X , β , σ 2 ) であるから、対数尤度は、
ℓ ( β , σ 2 ∣ X , y ) = log L ( β , σ 2 ∣ X , y ) = − n 2 log ( 2 π σ 2 ) − 1 2 σ 2 ( y − X β ) T ( y − X β ) \begin{aligned}
\ell (\boldsymbol{\beta}, \sigma^2| X, \boldsymbol{y}) & = \log L(\boldsymbol{\beta}, \sigma^2 | X, \boldsymbol{y})\\
& = - \frac{n}{2} \log (2\pi\sigma^2) - \frac{1}{2\sigma^2}(\boldsymbol{y}-X\boldsymbol{\beta})^{\mathsf{T}}(\boldsymbol{y}-X\boldsymbol{\beta})
\end{aligned} ℓ ( β , σ 2 ∣ X , y ) = log L ( β , σ 2 ∣ X , y ) = − 2 n log ( 2 π σ 2 ) − 2 σ 2 1 ( y − X β ) T ( y − X β ) となり、最尤推定量は、
β ^ = ( X T X ) − 1 X T y \hat{\boldsymbol{\beta}} = (X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}\boldsymbol{y} β ^ = ( X T X ) − 1 X T y σ ^ 2 = 1 n ( y − X β ^ ) T ( y − X β ^ ) \hat{\sigma}^2 = \frac{1}{n}(\boldsymbol{y}-X\hat{\boldsymbol{\beta}})^{\mathsf{T}}(\boldsymbol{y}-X\hat{\boldsymbol{\beta}}) σ ^ 2 = n 1 ( y − X β ^ ) T ( y − X β ^ ) 最小二乗法と同様の結果が得られた。
以下の実装では、説明変数同士を同じ尺度で評価するためにX X X 。
基準化(normalization)はいうまでもなく、平均値(μ \mu μ σ \sigma σ
z i = x i − μ σ z_i = \frac{x_i - \mu}{\sigma} z i = σ x i − μ import pandas as pd
import numpy as np
from numpy import linalg as la
df = pd. read_csv( "Boston.csv" , index_col= 0 )
y = df. iloc[ : , 13 ] . values
df = ( df - df. mean( ) ) / df. std( )
X = df. iloc[ : , : 13 ] . values
X = np. column_stack( ( np. ones( len ( X) ) , X) )
beta = np. dot( np. dot( la. inv( np. dot( X. T, X) ) , X. T) , y)
print ( beta) [ 2.25328063e+01 -9.29064566e-01 1.08263896e+00 1.41039433e-01
6.82414381e-01 -2.05875361e+00 2.67687661e+00 1.94853355e-02
-3.10711605e+00 2.66485220e+00 -2.07883689e+00 -2.06264585e+00
8.50108862e-01 -3.74733185e+00]import pandas as pd
import numpy as np
import sklearn. linear_model as lm
df = pd. read_csv( "Boston.csv" , index_col= 0 )
y = df. iloc[ : , 13 ] . values
df = ( df - df. mean( ) ) / df. std( )
X = df. iloc[ : , : 13 ] . values
model = lm. LinearRegression( )
model. fit( X, y)
print ( model. intercept_)
print ( model. coef_) 22.5328063241
[-0.92906457 1.08263896 0.14103943 0.68241438 -2.05875361 2.67687661
0.01948534 -3.10711605 2.6648522 -2.07883689 -2.06264585 0.85010886
-3.74733185]library( MASS)
X <- as.matrix( Boston[ , - 14 ] )
y <- as.vector( Boston[ , 14 ] )
X <- scale( X)
X <- cbind( rep( 1 , nrow( X) ) , X)
beta <- solve( t( X) %*% X) %*% t( X) %*% y
print( beta) [,1]
22.53280632
crim -0.92906457
zn 1.08263896
indus 0.14103943
chas 0.68241438
nox -2.05875361
rm 2.67687661
age 0.01948534
dis -3.10711605
rad 2.66485220
tax -2.07883689
ptratio -2.06264585
black 0.85010886
lstat -3.74733185library( MASS)
X <- as.matrix( Boston[ , - 14 ] )
y <- as.vector( Boston[ , 14 ] )
X <- scale( X)
model <- lm( y ~ X)
print( model) Call:
lm(formula = y ~ X)
Coefficients:
(Intercept) Xcrim Xzn Xindus Xchas
22.53281 -0.92906 1.08264 0.14104 0.68241
Xnox Xrm Xage Xdis Xrad
-2.05875 2.67688 0.01949 -3.10712 2.66485
Xtax Xptratio Xblack Xlstat
-2.07884 -2.06265 0.85011 -3.74733よって求めたかった予測モデルは、
住宅価格の中央値 = 22.53 − 0.93 × 犯罪率 + 1.08 × 宅地の割合 + ⋯ − 3.73 × 低所得者の割合 \begin{aligned}
住宅価格の中央値 = 22.53 - 0.93 \times 犯罪率 + 1.08 \times 宅地の割合 + \cdots -3.73 \times 低所得者の割合
\end{aligned} 住 宅 価 格 の 中 央 値 = 2 2 . 5 3 − 0 . 9 3 × 犯 罪 率 + 1 . 0 8 × 宅 地 の 割 合 + ⋯ − 3 . 7 3 × 低 所 得 者 の 割 合 のようになります。
PythonおよびRの実行結果から、ボストンの住宅価格には「低所得者層の割合」が一番影響を与えていることがわかる。当たり前といえば当たり前だ。
先ほどのモデルについて、β \boldsymbol{\beta} β ( X T X ) − 1 (X^{\mathsf{T}}X)^{-1} ( X T X ) − 1
β = ( X T X ) − 1 X T y \boldsymbol{\beta} = (X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}\boldsymbol{y} β = ( X T X ) − 1 X T y ( X T X ) − 1 (X^{\mathsf{T}}X)^{-1} ( X T X ) − 1 X X X X X X
説明変数間の相関が強いとき
n < p n < p n < p
他にもマルチコなどが考えられますが、私の理解が追いついていないので割愛します。(ごめんなさい)
仮に説明変数に全く同じものを追加して、先ほどのコードを実行してみます。
library( MASS)
X <- as.matrix( Boston[ , - 14 ] )
y <- as.vector( Boston[ , 14 ] )
X <- scale( X)
X <- cbind( X, X[ , 13 ] )
X <- cbind( rep( 1 , nrow( X) ) , X)
beta <- solve( t( X) %*% X) %*% t( X) %*% y solve.default(t(X) %*% X) でエラー:
Lapack routine dgesv: システムは正確に特異です: U[15,15] = 0となり、相関が1となる変数が2つ潜んでいるため、X X X ( X T X ) − 1 (X^{\mathsf{T}}X)^{-1} ( X T X ) − 1
> cor(X[, 13], X[, 14])
[1] 1では先程のケースではどうすればよいか。次の対処法が考えられる。
相関が強い説明変数同士を入れないようにする
擬似逆行列を用いいる。
( X T X ) (X^{\mathsf{T}}X) ( X T X )
太字部分が正則化で、今回はこれに注目してみます。これはどういうことかというと、( X T X ) − 1 (X^{\mathsf{T}}X)^{-1} ( X T X ) − 1 ( X T X + λ I p + 1 ) − 1 (X^{\mathsf{T}}X+ \lambda\boldsymbol{I_{p+1}})^{-1} ( X T X + λ I p + 1 ) − 1 I I I
β λ = ( X T X + λ I p + 1 ) − 1 X T y ⋯ ( ∗ ) \boldsymbol{\beta_{\lambda}} = (X^{\mathsf{T}}X+ \lambda\boldsymbol{I_{p+1}})^{-1}X^{\mathsf{T}}\boldsymbol{y} \cdots (*) β λ = ( X T X + λ I p + 1 ) − 1 X T y ⋯ ( ∗ ) (どうしてこれで逆行列が求まるようになるかは割愛。)
このように正則でない行列を正則にすることによって推定量を得る方法を正則化(regularization)というらしい。
λ ( > 0 ) \lambda(>0) λ ( > 0 ) ( ∗ ) (*) ( ∗ )
( ∗ ) (*) ( ∗ )
S λ ( β ) = ( y − X β ) T ( y − X β ) + λ β T β = ∣ ∣ y − X β ∣ ∣ 2 + λ ∣ ∣ β ∣ ∣ 2 \begin{aligned}
\boldsymbol{S_{\lambda}}(\boldsymbol{\beta})
& = (\boldsymbol{y}-X\boldsymbol{\beta})^{\mathsf{T}}(\boldsymbol{y}-X\boldsymbol{\beta}) + \lambda\boldsymbol{\beta^{\mathsf{T}}\beta} \\
& = ||\boldsymbol{y}-X\boldsymbol{\beta}||^2 + \lambda||\boldsymbol{\beta}||^2
\end{aligned} S λ ( β ) = ( y − X β ) T ( y − X β ) + λ β T β = ∣ ∣ y − X β ∣ ∣ 2 + λ ∣ ∣ β ∣ ∣ 2 を誤差関数としてβ \boldsymbol{\beta} β L2正則化 という。λ ∣ ∣ β ∣ ∣ 2 \lambda||\boldsymbol{\beta}||^2 λ ∣ ∣ β ∣ ∣ 2
S λ ( β ) = ∣ ∣ y − X β ∣ ∣ 2 + λ ∣ ∣ β ∣ ∣ 2 = ( y − X β ) T ( y − X β ) + λ β T β = ( y T − ( X β ) T ) ( y − X β ) + λ β T β = ( y T − β T X T ) ( y − X β ) + λ β T β = y T y − β T X T y − y T X β + β T X T X β + λ β T β = y T y − 2 β T X T y + β T X T X β + λ β T β \begin{aligned}
\boldsymbol{S_{\lambda}}(\boldsymbol{\beta})
& = ||\boldsymbol{y}-X\boldsymbol{\beta}||^2 + \lambda||\boldsymbol{\beta}||^2 \\
& = (\boldsymbol{y}-X\boldsymbol{\beta})^{\mathsf{T}}(\boldsymbol{y}-X\boldsymbol{\beta}) + \lambda\boldsymbol{\beta}^{\mathsf{T}}\boldsymbol{\beta} \\
& = (\boldsymbol{y}^{\mathsf{T}}-(X\boldsymbol{\beta})^{\mathsf{T}})(\boldsymbol{y}-X\boldsymbol{\beta})+ \lambda\boldsymbol{\beta}^{\mathsf{T}}\boldsymbol{\beta} \\
& = (\boldsymbol{y}^{\mathsf{T}}-\boldsymbol{\beta}^{\mathsf{T}}X^{\mathsf{T}})(\boldsymbol{y}-X\boldsymbol{\beta})+ \lambda\boldsymbol{\beta}^{\mathsf{T}}\boldsymbol{\beta} \\
& = \boldsymbol{y}^{\mathsf{T}}\boldsymbol{y} - \boldsymbol{\beta}^{\mathsf{T}}X^{\mathsf{T}}\boldsymbol{y} - \boldsymbol{y}^{\mathsf{T}}X\boldsymbol{\beta} + \boldsymbol{\beta}^{\mathsf{T}}X^{\mathsf{T}}X\boldsymbol{\beta} + \lambda\boldsymbol{\beta}^{\mathsf{T}}\boldsymbol{\beta} \\
& = \boldsymbol{y}^{\mathsf{T}}\boldsymbol{y} - 2 \boldsymbol{\beta}^{\mathsf{T}}X^{\mathsf{T}}\boldsymbol{y} + \boldsymbol{\beta}^{\mathsf{T}}X^{\mathsf{T}}X\boldsymbol{\beta}+ \lambda\boldsymbol{\beta}^{\mathsf{T}}\boldsymbol{\beta} \\
\end{aligned} S λ ( β ) = ∣ ∣ y − X β ∣ ∣ 2 + λ ∣ ∣ β ∣ ∣ 2 = ( y − X β ) T ( y − X β ) + λ β T β = ( y T − ( X β ) T ) ( y − X β ) + λ β T β = ( y T − β T X T ) ( y − X β ) + λ β T β = y T y − β T X T y − y T X β + β T X T X β + λ β T β = y T y − 2 β T X T y + β T X T X β + λ β T β これをβ \boldsymbol{\beta} β
S λ ( β ) ∂ β = − 2 X T y + 2 X T X β + 2 λ β \frac{\boldsymbol{S_{\lambda}}(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}= -2X^{\mathsf{T}}\boldsymbol{y} + 2X^{\mathsf{T}}X\boldsymbol{\beta} + 2\lambda\boldsymbol{\beta} ∂ β S λ ( β ) = − 2 X T y + 2 X T X β + 2 λ β S λ ( β ) ∂ β = 0 \frac{\boldsymbol{S_{\lambda}}(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = 0 ∂ β S λ ( β ) = 0
− 2 X T y + 2 X T X β + 2 λ β = 0 ( X T X + λ I p + 1 ) β = X T y β = ( X T X + λ I p + 1 ) − 1 X T y \begin{aligned}
-2 X^{\mathsf{T}}\boldsymbol{y} + 2X^{\mathsf{T}}X\boldsymbol{\beta} + 2\lambda\boldsymbol{\beta} = 0 \\
(X^{\mathsf{T}}X+\lambda\boldsymbol{I_{p+1}})\boldsymbol{\beta} = X^{\mathsf{T}} \boldsymbol{y} \\
\boldsymbol{\beta} = (X^{\mathsf{T}}X + \lambda \boldsymbol{I_{p+1}})^{-1}X^{\mathsf{T}}\boldsymbol{y}
\end{aligned} − 2 X T y + 2 X T X β + 2 λ β = 0 ( X T X + λ I p + 1 ) β = X T y β = ( X T X + λ I p + 1 ) − 1 X T y 最後に、相関が1となるようなデータが入っているときでも、L2正則化を行うことにより推定量が求まるのかRで実装してみます。
library( MASS)
X <- as.matrix( Boston[ , - 14 ] )
y <- as.vector( Boston[ , 14 ] )
X <- scale( X)
X <- cbind( X, X[ , 13 ] )
X <- cbind( rep( 1 , nrow( X) ) , X)
lambda <- 0.001
beta <- solve( t( X) %*% X + lambda* diag( ncol( X) ) ) %*% t( X) %*% y
print( beta) [,1]
22.53276179
crim -0.92905477
zn 1.08262398
indus 0.14101604
chas 0.68241750
nox -2.05872447
rm 2.67688066
age 0.01948267
dis -3.10708854
rad 2.66477986
tax -2.07876854
ptratio -2.06263709
black 0.85010730
lstat -1.87366482
-1.87366482今度はちゃんと求まりました！
今回は重回帰モデルとL2正則化について扱いました。次回はL1正則化、スパース推定について扱いたいと思います。
2016/01/09追記：この記事の続編は「Lassoの理論と実装 -スパースな解の推定アルゴリズム- 」になります。
