(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 × 住宅の割合 + ⋯ + β 13 × 低所得者層の割合 と表現できるような
β 0 , β 1 , ⋯ , β 13 \beta_0, \beta_1, \cdots, \beta_{13} β 0 , β 1 , ⋯ , β 13
を求めたい。
もう少し一般化したものを考えてみます。以下の表の目的変数 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 11 x 12 x_{12} x 12 ⋯ \cdots ⋯ x 1 p x_{1p} x 1 p 2 2 2 y 2 y_2 y 2 x 21 x_{21} x 21 x 22 x_{22} x 22 ⋯ \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 11 ⋮ 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 11 ⋮ 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 :n次元の観測ベクトルX X X :n×(p+1)次元の計画行列(design matrix) β \boldsymbol{\beta} β :(p+1)次元の回帰係数ベクトル (これを求めたい。 )ϵ \boldsymbol{\epsilon} ϵ :n次元の誤差ベクトルとなります。
未知パラメータであるβ \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 σ )で割ることにより、平均値0分散1にすることです。
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.74733185
library( 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} 住宅価格の中央値 = 22.53 − 0.93 × 犯罪率 + 1.08 × 宅地の割合 + ⋯ − 3.73 × 低所得者の割合 のようになります。
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 ) は正則化パラメータと呼ばれる。 特に、( ∗ ) (*) ( ∗ ) はL2正則化による推定量である(後述)。
( ∗ ) (*) ( ∗ ) を得るには、
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の理論と実装 -スパースな解の推定アルゴリズム- 」になります。