[mathjax]
タイトルの通り。Losso回帰と違って損失関数を偏微分するだけで出来そうなのでやってみる。
Ridge回帰は線形回帰の1種だけれども、損失関数として最小二乗法をそのまま使わず、
(L_2)ノルムの制約を付けたものを使う((L_2)正則化)。
データとモデル
教師データ(boldsymbol{y})、訓練データ(boldsymbol{x})があるとする。
(または目的変数(boldsymbol{y})、説明変数(boldsymbol{x})があるとする。)
例えば(p)次の属性データが(n)個あり、それらと結果の対応が分かっている状況。
begin{eqnarray}
boldsymbol{y} &=& begin{pmatrix}
y_1 \\
y_2 \\
vdots \\
y_p
end{pmatrix} ,
boldsymbol{x} &=& begin{pmatrix}
x_{11} & x_{21} & cdots & x_{n1} \\
x_{12} & x_{22} & cdots & x_{n2} \\
vdots & vdots & ddots & vdots \\
x_{1p} & x_{2p} & cdots & x_{np}
end{pmatrix}
end{eqnarray}
モデルは以下。特徴ベクトル(boldsymbol{w})は訓練データの重み。
特徴空間において損失を最小化する特徴ベクトルを求める問題。
begin{eqnarray}
boldsymbol{y} &=& boldsymbol{w} boldsymbol{x} + k \\
boldsymbol{w} &=& begin{pmatrix}
w_1 & w_2& cdots &w_p
end{pmatrix}
end{eqnarray}
損失関数
普通の2乗損失に正則化項((L_2)ノルムを定数倍した値)を付けたものを損失関数として利用する。
正則化項の係数はハイパーパラメータとして調整する値。逆数なのはsklearnに従う。
begin{eqnarray}
L(boldsymbol{w}) = |boldsymbol{y} - boldsymbol{w} boldsymbol{x}|^2 +C |boldsymbol{w}|^2
end{eqnarray}
特徴ベクトルは以下。(mathjaxでargminが出せない...)
begin{eqnarray}
newcommand{argmin}[1]{underset{#1}{operatorname{arg},operatorname{min}};}
boldsymbol{w} = argmin w L(boldsymbol{w}) = argmin w |boldsymbol{y} - boldsymbol{w} boldsymbol{x}|^2 + C |boldsymbol{w}|^2
end{eqnarray}
特徴ベクトルを求める
勾配=0と置けば上の式の解を得られる。
損失関数が微分可能だからできる技。
begin{eqnarray}
frac{partial L(boldsymbol{w})}{partial boldsymbol{w}} &=& 2 boldsymbol{w}^T (boldsymbol{y} - boldsymbol{w} boldsymbol{x}) + C boldsymbol{w} \\
&=& 0
end{eqnarray}
変形する。
begin{eqnarray}
2 boldsymbol{x}^T (boldsymbol{x}boldsymbol{w}-boldsymbol{y}) + C boldsymbol{w} &=& 0 \\
boldsymbol{x}^T (boldsymbol{x}boldsymbol{w}-boldsymbol{y}) + C boldsymbol{w} &=& 0 \\
boldsymbol{x}^T boldsymbol{x} boldsymbol{w} -boldsymbol{x}^T boldsymbol{y} + Cboldsymbol{w} &=& 0 \\
(boldsymbol{x}^T boldsymbol{x} +C E) boldsymbol{w} &=& boldsymbol{x}^T boldsymbol{y} \\
boldsymbol{w} &=& (boldsymbol{x}^T boldsymbol{x} + C E)^{-1} boldsymbol{x}^T boldsymbol{y}
end{eqnarray}
テストデータを作る
練習用にsklearnのbostonデータを使ってみる。
ボストンの住宅価格が目的変数、属性データが説明変数として入ってる。
import pandas as pd
import numpy as np
from pandas import Series,DataFrame
import matplotlib.pyplot as plt
from sklearn.datasets import load_boston
boston = load_boston()
boston_df = DataFrame(boston.data)
boston_df.columns = boston.feature_names
print(boston_df.head())
boston_df[\"PRICE\"] = DataFrame(boston.target)
# CRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX PTRATIO B LSTAT PRICE
# 0 0.00632 18.0 2.31 0.0 0.538 6.575 65.2 4.0900 1.0 296.0 15.3 396.90 4.98 24.0
# 1 0.02731 0.0 7.07 0.0 0.469 6.421 78.9 4.9671 2.0 242.0 17.8 396.90 9.14 21.6
# 2 0.02729 0.0 7.07 0.0 0.469 7.185 61.1 4.9671 2.0 242.0 17.8 392.83 4.03 34.7
# 3 0.03237 0.0 2.18 0.0 0.458 6.998 45.8 6.0622 3.0 222.0 18.7 394.63 2.94 33.4
# 4 0.06905 0.0 2.18 0.0 0.458 7.147 54.2 6.0622 3.0 222.0 18.7 396.90 5.33 36.2
散布図行列を表示してみる。
PRICEと関係がありそうなZN,RM,AGE,DIS,LSTATの5個を使ってみる。
pg = sns.pairplot(boston_df)
plt.show()
pg.savefig(\'boston_fig.png\')
特徴ベクトルを自力で計算する
これを自力で計算してみる。(C=0.01)、(C=0)、(C=100)としてみた。
begin{eqnarray}
boldsymbol{w} &=& (boldsymbol{x}^T boldsymbol{x} + C E)^{-1} boldsymbol{x}^T boldsymbol{y}
end{eqnarray}
X_df = boston_df.drop(columns=[\'CRIM\',\'INDUS\',\'CHAS\',\'NOX\',\'RAD\',\'TAX\',\'PTRATIO\',\'B\',\'PRICE\'])
X = X_df.values
y = boston.target.T
C1 = 0.01
C2 = 0
C3 = 100
e = np.identity(5)
w1 = np.dot( np.linalg.inv(np.dot(X.T , X) + C1 * e), np.dot(X.T,y))
w2 = np.dot( np.linalg.inv(np.dot(X.T , X) + C2 * e), np.dot(X.T,y))
w3 = np.dot( np.linalg.inv(np.dot(X.T , X) + C3 * e), np.dot(X.T,y))
print(w1)
# [ 0.05338557 5.40396159 -0.01209002 -0.83723303 -0.63725397]
print(w2)
# [ 0.05338539 5.40403743 -0.01209427 -0.83728837 -0.63725093]
print(w3)
# [ 0.05612977 4.76664789 0.02374402 -0.38576708 -0.66137596]
(C=0)のとき、つまり最小二乗法のとき。
sklearnを使う
sklearnのridge回帰モデルを使うと以下みたいになる。
from sklearn.linear_model import Ridge
from sklearn.model_selection import train_test_split
Xf_train,Xf_test,yf_train,yf_test = train_test_split(X,y,random_state=0)
ridge = Ridge().fit(Xf_train,yf_train)
print(f\"accuracy for training data:{ridge.score(Xf_train,yf_train):.2}\")
print(f\"accuracy for test data:{ridge.score(Xf_test,yf_test):.2f}\")
# accuracy for training data:0.68
# accuracy for test data:0.58
print(ridge.coef_)
# [ 0.06350701 4.3073956 -0.02283312 -1.06820241 -0.73188192]
出てきた特徴ベクトルを並べてみる
自力で計算したものとsklearnに計算してもらったものを並べてみる。
似てるのか似ていないのかよくわからない .. けど、RMの寄与度が高いというのは似ている。
# 自力で計算 (C=100)
# [ 0.05612977 4.76664789 0.02374402 -0.38576708 -0.66137596]
# sklearnで計算
# [ 0.06350701 4.3073956 -0.02283312 -1.06820241 -0.73188192]
自力で計算したモデルの正答率を求めてみないとなんとも...
そして、正規化項の係数の大小がどう影響するのか、あまり良くわからなかった..。
(L_2)ノルムの制約を付けると、パラメタの大小が滑らかになると言いたかったのだけども。
あと、訓練データに対して68%、テストデータに対して58%という感じで、
大して成績が良くない...。
