Published on

重回帰モデルの理論と実装 -なぜ正則化が必要か-

Authors

導入

(Qiitaに投稿していた当時)Qiitaで重回帰と検索をかけても意外と数式での説明がなかったので今回は数式で攻めたいと思います。

例題として、The Boston Housing Datasetを使います。

crimzninduschasnoxrmagedisradtaxptratioblacklstatmedv
10.00632182.3100.5386.57565.24.0900129615.3396.904.9824.0
20.0273107.0700.4696.42178.94.9671224217.8396.909.1421.6
30.0272907.0700.4697.18561.14.9671224217.8392.834.0334.7
40.0323702.1800.4586.99845.86.0622322218.7394.632.9433.4
50.0690502.1800.4587.14754.26.0622322218.7396.905.3336.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,,β13\beta_0, \beta_1, \cdots, \beta_{13}

を求めたい。

モデル式の確認

もう少し一般化したものを考えてみます。以下の表の目的変数 yypp 個の説明変数 x1xpx_1 \cdots x_p で予測するモデルを作成したい。

indexyyx1x_1x2x_2\cdotsxpx_p
11y1y_1x11x_{11}x12x_{12}\cdotsx1px_{1p}
22y2y_2x21x_{21}x22x_{22}\cdotsx2px_{2p}
\vdots\vdots\vdots\vdots\vdots\vdots
iiyiy_ixi1x_{i1}xi2x_{i2}\cdotsxipx_{ip}
\vdots\vdots\vdots\vdots\vdots\vdots
nnyny_nxn1x_{n1}xn2x_{n2}\cdotsxnpx_{np}

重回帰モデル(multiple regression model) は以下のように表されます。

y=β0+β1x1++βpxp+ϵy = \beta_0 + \beta_1x_1 + \cdots + \beta_px_p + \epsilon

このとき、上の表においてii番目データは

yi=β0+β1xi1++βpxip+ϵiy_i = \beta_0 + \beta_1x_{i1} + \cdots + \beta_px_{ip} + \epsilon_i

と表せる。 ベクトルを用いて表すと、

[y1yn]=[1x11x1p1xn1xnp][β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=[y1yn]\boldsymbol{y} = \begin{bmatrix} y_1 \\ \vdots \\ y_n \end{bmatrix} X=[1x11x1p1xn1xnp]X= \begin{bmatrix} 1 & x_{11} & \cdots & x_{1p} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n1} & \cdots & x_{np} \end{bmatrix} β=[β0βp]\boldsymbol{\beta}= \begin{bmatrix} \beta_0 \\ \vdots \\ \beta_p \end{bmatrix} ϵ=[ϵ1ϵn]\boldsymbol{\epsilon}= \begin{bmatrix} \epsilon_1 \\ \vdots \\ \epsilon_n \end{bmatrix}

とおくと、

y=Xβ+ϵ\boldsymbol{y} = X\boldsymbol{\beta} + \boldsymbol{\epsilon}

と表すことができます。 確認ですが、

  • y\boldsymbol{y}:n次元の観測ベクトル
  • XX:n×(p+1)次元の計画行列(design matrix)
  • β\boldsymbol{\beta}:(p+1)次元の回帰係数ベクトル (これを求めたい。)
  • ϵ\boldsymbol{\epsilon}:n次元の誤差ベクトル

となります。

最小二乗法による推定

未知パラメータであるβ\boldsymbol{\beta}を最小二乗法により推定していきます。 二乗誤差である

S(β)=i=1nϵi2=ϵTϵ=(yXβ)T(yXβ)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})

を最小にするβ\boldsymbol{\beta}を求める。式を展開して、

S(β) =(yXβ)T(yXβ)=(yT(Xβ)T)(yXβ)=(yTβTXT)(yXβ)=yTyβTXTyyTXβ+βTXTXβ=yTy2yTXβ+βTXTXβ\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}

これをβ\boldsymbol{\beta}で微分して、

S(β)β=2XTy+2XTXβ\frac{S(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}= -2X^{\mathsf{T}}\boldsymbol{y} + 2X^{\mathsf{T}}X\boldsymbol{\beta}

S(β)β=0\frac{S(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = 0とおいて、

2XTy+2XTXβ=0XTXβ=XTy\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}

もし(XTX)1(X^{\mathsf{T}}X)^{-1}が存在するなら、

β=(XTX)1XTy\boldsymbol{\beta} = (X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}\boldsymbol{y}

と求めることができる。

最尤法による推定

参考:データ解析 Rによる多変量解析入門 (7) 最尤法、モデル選択

重回帰式

y=Xβ+ϵ\boldsymbol{y} = X\boldsymbol{\beta} + \boldsymbol{\epsilon}

の誤差ベクトル ϵ\boldsymbol{\epsilon} が多変量正規分布に従うものと仮定する。

ϵNn(0,σ2In)\boldsymbol{\epsilon} \sim N_n(\boldsymbol{0}, \sigma^2 I_n)

このとき、

yNn(Xβ,σ2In)\boldsymbol{y} \sim N_n(X\boldsymbol{\beta}, \sigma^2 I_n)

となり、確率密度関数は、

f(yX,β,σ2)=1(2πσ2In)n/2exp[12σ2(yXβ)T(yXβ)]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]。

尤度は、

L(β,σ2X,y)=f(yX,β,σ2)L(\boldsymbol{\beta}, \sigma^2 | X, \boldsymbol{y})=f(\boldsymbol{y} | X, \boldsymbol{\beta}, \sigma^2)

であるから、対数尤度は、

(β,σ2X,y)=logL(β,σ2X,y)=n2log(2πσ2)12σ2(yXβ)T(yXβ)\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}

となり、最尤推定量は、

β^=(XTX)1XTy\hat{\boldsymbol{\beta}} = (X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}\boldsymbol{y} σ^2=1n(yXβ^)T(yXβ^)\hat{\sigma}^2 = \frac{1}{n}(\boldsymbol{y}-X\hat{\boldsymbol{\beta}})^{\mathsf{T}}(\boldsymbol{y}-X\hat{\boldsymbol{\beta}})

最小二乗法と同様の結果が得られた。

実装

以下の実装では、説明変数同士を同じ尺度で評価するためにXXの各説明変数毎に基準化を行っています。 基準化(normalization)はいうまでもなく、平均値(μ\mu)を引いて標準偏差(σ\sigma)で割ることにより、平均値0分散1にすることです。

zi=xiμσz_i = \frac{x_i - \mu}{\sigma}

Python

行列演算による実装

  • コード
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]

scikit-learnによる実装

  • コード
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]

R

行列演算による実装

  • コード
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

lm関数による実装

  • コード
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.530.93×犯罪率+1.08×宅地の割合+3.73×低所得者の割合\begin{aligned} 住宅価格の中央値 = 22.53 - 0.93 \times 犯罪率 + 1.08 \times 宅地の割合 + \cdots -3.73 \times 低所得者の割合 \end{aligned}

のようになります。

PythonおよびRの実行結果から、ボストンの住宅価格には「低所得者層の割合」が一番影響を与えていることがわかる。当たり前といえば当たり前だ。

正則化

なぜ正則化するのか

最小二乗推定量を求められないケース

先ほどのモデルについて、β\boldsymbol{\beta}を求めるには、(XTX)1(X^{\mathsf{T}}X)^{-1}が存在する必要があります。

β=(XTX)1XTy\boldsymbol{\beta} = (X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}\boldsymbol{y}

(XTX)1(X^{\mathsf{T}}X)^{-1}が存在しないとはどのようなケースか考えてみると、それはXXの階数がフルランクでない時だとわかります。XXの階数がフルランクにならないのは次のような原因が考えられます。

  • 説明変数間の相関が強いとき
  • n<pn < pとなるとき

他にもマルチコなどが考えられますが、私の理解が追いついていないので割愛します。(ごめんなさい)

仮に説明変数に全く同じものを追加して、先ほどのコードを実行してみます。

  • Rコード
library(MASS)
X <- as.matrix(Boston[, -14])
y <- as.vector(Boston[,  14])
X <- scale(X)
X <- cbind(X, X[, 13]) # 13番目の説明変数を14番目の説明変数をとして追加
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つ潜んでいるため、XXの階数がフルランクでなくなり、(XTX)1(X^{\mathsf{T}}X)^{-1}を求めることできなくなります。

> cor(X[, 13], X[, 14])
[1] 1

最小二乗推定量を求められないケースへの対処法

では先程のケースではどうすればよいか。次の対処法が考えられる。

  • 相関が強い説明変数同士を入れないようにする
  • 擬似逆行列を用いいる。
  • (XTX)(X^{\mathsf{T}}X)の逆行列が求まるようにしてあげる

太字部分が正則化で、今回はこれに注目してみます。これはどういうことかというと、(XTX)1(X^{\mathsf{T}}X)^{-1}(XTX+λIp+1)1(X^{\mathsf{T}}X+ \lambda\boldsymbol{I_{p+1}})^{-1}のように何か足すことによって逆行列が求まるようにするということ。(I Iは単位行列。)つまり求める推定量は次のようになる。

βλ=(XTX+λIp+1)1XTy()\boldsymbol{\beta_{\lambda}} = (X^{\mathsf{T}}X+ \lambda\boldsymbol{I_{p+1}})^{-1}X^{\mathsf{T}}\boldsymbol{y} \cdots (*)

(どうしてこれで逆行列が求まるようになるかは割愛。) このように正則でない行列を正則にすることによって推定量を得る方法を正則化(regularization)というらしい。 λ(>0)\lambda(>0)は正則化パラメータと呼ばれる。 特に、()(*)はL2正則化による推定量である(後述)。

L2正則化

導出

()(*)を得るには、

Sλ(β)=(yXβ)T(yXβ)+λβTβ=yXβ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}

を誤差関数としてβ\boldsymbol{\beta}を推定してやればよい。これをL2正則化という。λβ2\lambda||\boldsymbol{\beta}||^2はペナルティ項と呼ばれ、これを設けることにより、過学習を抑え汎化能力を高めることができる。実際に導出してみる。

Sλ(β)=yXβ2+λβ2=(yXβ)T(yXβ)+λβTβ=(yT(Xβ)T)(yXβ)+λβTβ=(yTβTXT)(yXβ)+λβTβ=yTyβTXTyyTXβ+βTXTXβ+λβTβ=yTy2βTXTy+βTXTXβ+λβ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}

これをβ\boldsymbol{\beta}で微分して、

Sλ(β)β=2XTy+2XTXβ+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λ(β)β=0\frac{\boldsymbol{S_{\lambda}}(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = 0とおいて、

2XTy+2XTXβ+2λβ=0(XTX+λIp+1)β=XTyβ=(XTX+λIp+1)1XTy\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}

実装

最後に、相関が1となるようなデータが入っているときでも、L2正則化を行うことにより推定量が求まるのかRで実装してみます。

  • Rコード
library(MASS)
X <- as.matrix(Boston[, -14])
y <- as.vector(Boston[,  14])
X <- scale(X)
X <- cbind(X, X[, 13]) # 13番目の説明変数を14番目の説明変数をとして追加
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の理論と実装 -スパースな解の推定アルゴリズム-」になります。

参考