今回は、
Lasso回帰は多くのパラメータを０にするので説明するのに役立つ
と巷では言われるLasso回帰を実装していきます！
重回帰、リッジ回帰について理解をしてから読むとより分かりやすいと思います。
もし、まだ読まれていなければ先にこちらの記事をオススメします。
個人的には精度や、なぜその変数を０にしたのかを、万人が納得できるように説明するのが難しいので好きではありません。笑
ですが、それも実装をして自分で使って見なければ分からないことだと思うので、やっていきます！
Contents
ラッソ(Lasso)回帰とは？
重回帰のパラメータを求める際に、正則化項を入れることで、重回帰よりもロバストな結果を返すモデルです。（一応？？）
これだけを見るとリッジ回帰と全く同じなのですが、正則化項の定義が異なります。
リッジ回帰ではL2ノルム、ラッソ回帰ではL1ノルムを用います。
L1ノルムを用いることで、様々な効果をもたらします。（良い意味でも悪い意味でも。）
L1ノルムを使うことで、微分をする際に、場合分けをしなければならないという点です。
数学がとても複雑になります。そのおかげで実装も面倒です。
場合分けをしなければならないのなら、そのついでに、パラメータ数を減らしてしまえ！ということで、複数の$w$の値を0にして回帰をすることができます。
そのため、解釈がしやすくなると言われています。
個人的にはLassoのパラメータから万人が納得できる解釈をすることはできないと思っています。
理由は、似たような働きをする特徴量$x_1, x_2$に対して、ラッソ回帰の重み$w_1$, $w_2$があったときに、
$w_1=0$,$w_2=2.5$のように決定された場合、それは、$w_1$と比較して $w_2$が意味のある数値だからと言い切ることはできないからです。
なぜなら、$w_1$, $w_2$の値については、訓練データのサンプリング方法や、後ほど出てくる繰り返し回数によって影響を受けるからです。
ただ、パラメータの値を、参考にすることはできると思います。
それを元に、もっと解釈がしやすいモデルで特徴量を削減して入力をするなどすれば分かりやすい結果を得ることはできると思います。
使いどころの問題です。
文字の説明
- $y$ : 目的変数 ( 従属変数 )
- $x$: 説明変数 ( 独立変数 )
- $ϕ_j$ : 基底関数 ( $x_j$, $x_j^2$, $ln{(x_j)}$ など )
- $ω_j$ : 回帰係数
- $ε$ : 撹乱項、誤差項
- $λ_1$：正則化項 ← New!!
どんな考え方をするの？
基本的には重回帰と全く同じです。$y$を求めるための$w$を求めます。
$y =ω_0+ω_1 ϕ_1 (x)+ω_2 ϕ_2 (x)+⋯+ω_{K-1} ϕ_{K-1} (x)+ε$
$=w^T ϕ+ε$
ただし、
$w= (ω_0,ω_1,⋯,ω_{K-1})^T, ϕ=(ϕ_0 (x),ϕ_1 (x) , ⋯,ϕ_{K-1} (x))^T$
すなわち、$y$を求めるために必要なものは、
- 重みとなる$ω$
- 独立関数$ϕ$
この二つです。
$ϕ$については、$x$の何乗までを入力にしてあげるかを関数の中で決めれるようにします。
なので、ここでは解説を飛ばします。
どうやって$ω$を求めるの？
良い感じの$f(x)$について考えるために誤差$E$を調べます。
ここで初めて、リッジ回帰との違いが現れます！
$E = \frac{1}{2}||w^T ϕ – y || + λ_1|w|$
のように、$λ_1|w|$のようにL1ノルムと呼ばれる、正則化項を加えます。
これによって、$w$に対して制限を加えます。
ただ、L1ノルムでは、L2ノルムのようにパラメータの大きさを制限するものではないので、パラメータが発散することがあり得ます。
それ以外はおおむね一緒です。
やはり、$E$が最小になってくれると嬉しいので、$E$が最小となる$ω$を計算します。
最初値を求めたいので、$E$を$ω$で微分した結果が0になればよいですが、絶対値の微分になるので、場合分けをしなければなりません。
その点について以下の記事が分かりやすかったです。
更新の部分については上記の記事を参考にしました。
ラッソ（lasso）回帰の実装
class MyLassoRegressor:
def __init__(self, alpha=1.0, max_iter=1000):
"""
alpha : float, default=1.0
Regularization strength
max_iter : int, default=1000
Iterate count to update the weight
"""
self.max_iter = int(max_iter) # New!!
self.alpha = float(alpha)
self.weight = None # New!!
self.dim = None
self.coef_ = None
self.intercept_ = None
self.n_feature = None
self.__name__ = "Lasso"
def _thresholding(self, numerator, denominator, n_sample): # New!!
# wk is positive
if numerator > self.alpha * n_sample:
return (numerator - self.alpha * n_sample) / denominator
# wk is negative
elif numerator < -self.alpha * n_sample:
return (numerator + self.alpha * n_sample) / denominator
# others
else:
return 0
def _update(self, n_sample, n_weight, X, y): # New!!
for k in range(n_weight):
if k == 0:
# update w0
self.weight[0] = (y - np.dot(X[:, 1:], self.weight
)).sum() / n_sample
else:
# update index
_k = [i for i in range(n_weight) if i not in [0, k]]
# # a part of renewal numerator
numerator = np.dot((y - np.dot(X[:, _k], self.weight[_k]) - self.weight[0]), X[:, k]).sum()
# renewal denominator
denominator = (X[:, k] ** 2).sum()
# update weights
self.weight[k] = self._thresholding(numerator, denominator, n_sample)
def fit(self, X, y, dim=1):
"""
Calculate weight of Lasso Regression
Parameters
X : array_like (n_sample, n_feature)
The matrix of explanatory variables
y : array_like (n_sample, 1)
Objective variable
dim : int, default 1
Dimension of multiple regression
"""
n_sample, n_feature = X.shape
# X (n_sample, dim, n_feature)
X = np.array([[x**i for i in range(dim + 1)] for x in X])
# delete duplicated ones array
X = X.reshape(n_sample, -1)[:, n_feature - 1:]
n_weight = X.shape
# initialize weight as 0
self.weight = np.zeros(n_weight)
# update
for _ in range(self.max_iter):
self._update(n_sample, n_weight, X, y)
# print(self.weight)
self.coef_ = self.weight
self.intercept_ = self.weight[0]
self.dim = dim
self.n_feature = n_feature
def predict(self, X):
"""
Calculate Lasso Regression
Parameters
X : array_like (n_sample, n_feature)
The matrix of explanatory variables in test dataset
Returns
Predicted objective variable
"""
n_sample, n_feature = X.shape
# X (n_sample, dim, n_feature)
X = np.array([[x**i for i in range(self.dim + 1)] for x in X])
# delete ones array
X = X.reshape(n_sample, -1)[:, n_feature:]
return (np.dot(X, self.coef_) + self.intercept_)
def mse(self, X, y):
"""
Calculate mean square error
Parameters
X : array_like (n_sample, n_feature)
The matrix of explanatory variables in test dataset
y : Ground truth of objective variable
Returns
Mean square error
"""
return np.mean((self.predict(X) - y)**2)
def title(self):
"""
Make title for plot
Return
string of equation
"""
Xn = [f"x_{k}" for k in range(1, self.n_feature + 1)]
title = "$f(x) = $" + f"{float(self.intercept_): .2f}"
for d in range(1, self.dim + 1):
for k in range(self.n_feature):
if self.coef_[(d - 1) * self.n_feature + k] > 0:
sign = "+"
elif self.coef_[(d - 1) * self.n_feature + k] < 0:
sign = "-"
else:
continue
if d == 1:
title += f"${sign}{abs(float(self.coef_[(d-1) * self.n_feature + k])):.2f}{Xn[k]}$"
else:
title += f"${sign}{abs(float(self.coef_[(d-1) * self.n_feature + k])):.2f}{Xn[k]}^{d}$"
return title
_thresholdindの部分で、各weightに対して場合分けを行い、更新をするか0にするかを決定します。
この計算によってweightが0になった理由を、定性的(数学をよく知らない偉そうな人が直感的に理解できる感じ)に説明できますか？？笑
とても難しいと思います。
これが僕がLasso回帰をあまり推していない理由のひとつです。
また、この更新では、weightが更新によって符号をまたがないことを暗黙に仮定しています。
つまり、更新をすればするほど、weightの絶対値は大きくなります。(ただし足したときの和は一定)
それがどういった影響を及ぼすのかも考えてみると面白いかも知れません。
それでも、精度が良ければ使えば良いと思うんですよ。
はい。では、精度を見ていきます。
ラッソ (lasso) 回帰を適応
import argparse
import pandas as pd
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from sklearn.model_selection import KFold
def main():
parser = argparse.ArgumentParser()
parser.add_argument("-d", "--dim", type=int, default=1, help="回帰の次元数．デフォルトは1")
parser.add_argument("-r", "--regressor", type=str, default="m", help="回帰モデルの名前．デフォルトは'm'")
parser.add_argument("-f", "--fname", type=str, default="data1.csv", help="読み込むファイルの名前．デフォルトは'data1.csv'")
parser.add_argument("-a", "--alpha", type=float, default=1.0, help="正則化項の係数．デフォルトは1.0")
args = parser.parse_args()
# define regression model
dim = args.dim
alpha = args.alpha
regressor = args.regressor
if regressor == "m":
model = MyMultipleRegressor()
elif regressor == "r":
model = MyRidgeRegressor(alpha=alpha)
elif regressor == "l":
model = MyLassoRegressor(alpha=alpha, max_iter=2)
# load data
fname = args.fname
data = pd.read_csv('../' + fname)
# split data
n_sample, n_cols = data.shape
feature_cols = list(data.columns[:-1])
target_col = list(data.columns[-1:])
kf = KFold(n_splits=2, shuffle=True)
for train_idx, test_idx in kf.split(data):
break
# train data
train_X = data.loc[train_idx, feature_cols].values.reshape(-1, n_cols - 1)
train_y = data.loc[train_idx, target_col].values.reshape(-1, 1)
model.fit(train_X, train_y, dim=dim)
train_mse = model.mse(train_X, train_y)
# test data
test_X = data.loc[test_idx, feature_cols].values.reshape(-1, n_cols - 1)
test_y = data.loc[test_idx, target_col].values.reshape(-1, 1)
test_mse = model.mse(test_X, test_y)
title = model.title()
if n_cols == 2:
# plot original data
plt.figure(figsize=(10, 8))
plt.xlabel(feature_cols[0], size=16)
plt.ylabel(target_col[0], size=16)
plt.title(f"{model.__name__} Regression\n" + title + f"\nTrain MSE : {train_mse:.4f}, Test MSE : {test_mse:.4f}", size=20)
plt.scatter(train_X, train_y, label="Train Observed")
plt.scatter(test_X, test_y, label="Test Observed")
# plot predicted data
X1 = np.linspace(data[feature_cols].min(), data[feature_cols].max(), 500)
y1_hat = model.predict(X1)
plt.plot(X1, y1_hat, c="r", label="Predicted")
plt.legend(fontsize=16)
plt.savefig(fname[:-4] + f"{model.__name__}_{dim}D.png", bbox_inches="tight")
plt.show()
elif n_cols == 3:
# plot original data
fig = plt.figure(figsize=(10, 8))
ax = Axes3D(fig)
ax.set_xlabel(feature_cols[0], size=14)
ax.set_ylabel(feature_cols, size=14)
ax.set_zlabel(target_col[0], size=14)
plt.title(f"{model.__name__} Regression\n" + title + f"\nTrain MSE : {train_mse:.4f}, Test MSE : {test_mse:.4f}", size=20)
ax.plot(data.loc[train_idx, feature_cols[0]].values, data.loc[train_idx, feature_cols].values,
data.loc[train_idx, target_col].values.flatten(),
marker="o", linestyle='None', label="Train Observed")
ax.plot(data.loc[test_idx, feature_cols[0]].values, data.loc[test_idx, feature_cols].values,
data.loc[test_idx, target_col].values.flatten(),
marker="o", linestyle='None', label="Test Observed")
# plot predicted data
X1 = np.linspace(data[feature_cols[0]].min(), data[feature_cols[0]].max(), 30)
X2 = np.linspace(data[feature_cols].min(), data[feature_cols].max(), 30)
X1, X2 = np.meshgrid(X1, X2)
X3 = np.concatenate([X1.reshape(-1, 1), X2.reshape(-1, 1)], axis=1)
y3_hat = model.predict(X3).reshape(X1.shape)
ax.plot_wireframe(X1, X2, y3_hat, color="red", label="Predictied")
plt.legend(loc='center left', fontsize=16)
plt.savefig(fname[:-4] + f"{model.__name__}_{dim}D.png", bbox_inches="tight")
plt.show()
if __name__ == "__main__":
main()
データを見て比較
data2（正則化項による差を確認）
左上から順に正則化項は
$λ_1 = 3175, 3225, 3275$です。
$f(x)$の独立変数は、どの式でも$(1, x, x^2, x^3)$を入力しました。
これを見ると確かに、一部の重みが0となることでパラメータが削減されたことが分かります！
そして、見てください。この正則化項の微妙な調整を。
この数字を出すために何度調整したことでしょうか。$λ_1 = 1e-4~ 1e6$で調べました。
この3225付近という微妙な数値がフィッティングの努力を物語っていると思います。
それなのに、この回帰曲線です。
なんか、むなしくなりませんか！？
実装が間違っていることを疑いました。
そこで、https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Lasso.html
ドキュメントにあるものでも同じパラメータを設定して実行をしました。
なんということでしょう。
全く同じ数値がでました。
この問題が偶然難しかったのかなとおもったので、他のデータでも試しました。
data1
data3
他のデータでも似た結果になりました。
さすがに偶然じゃないですね。
これが、僕があまりLasso回帰を推していない理由です。
最後に
Lasso回帰はいかがでしたでしょうか？
いろいろガッカリしたかもしれません。（恣意的に笑）
それでも、パラメータがかなり多いときに削減をしてくれるのは嬉しいことだと思います。
削減をしてもそれなりの精度が出せるように調整できれば使用する価値はあると思います。
やはり、使用目的にあったモデルの選択が重要という結論は変わらないです。
以上で全3回となる、回帰のフルスクラッチ実装終了です。
自分で実装をすることでかなり理解が深まりました。
これからも、学習をしたことをブログにまとめていきたいと思います。
こんな感じのブログを書いていきます。
もし興味を持っていただけたら、フォローお願いします。
タイピングもままならない完全にプログラミング初心者から
たった二ヶ月で
応用も簡単にできる…！！
という状態になるまで、一気に成長させてくれたオススメのプログラミングスクールをご紹介します！
