机器学习-Sklearn-12(回归类大家族-上——多元线性回归、岭回归、Lasso)(解决多重共线性)
1 概述
1.1 线性回归大家族
回归是一种应用广泛的预测建模技术,这种技术的核心在于预测的结果是连续型变量。其在现实中的应用
非常广泛,包括使用其他经济指标预测股票市场指数,根据喷射流的特征预测区域内的降水量,根据公司的广告花费预测总销售额,或者根据有机物质中残留的碳-14的量来估计化石的年龄等等,只要一切基于特征预测连续型变量的需求,我们都使用回归技术。
既然线性回归是源于统计分析,是结合机器学习与统计学的重要算法。通常来说,我们认为统计学注重先验,而机器学习看重结果,因此机器学习中不会提前为线性回归排除共线性等可能会影响模型的因素,反而会先建立模型以查看效果。模型确立之后,如果效果不好,我们就根据统计学的指导来排除可能影响模型的因素。我们的课程会从机器学习的角度来为大家讲解回归类算法,如果希望理解统计学角度的小伙伴们,各种统计学教材都可以满足你的需求。
回归需求在现实中非常多,所以我们自然也有各种各样的回归类算法。最著名的就是我们的线性回归和逻辑回归,从他们衍生出了岭回归,Lasso,弹性网,除此之外,还有众多分类算法改进后的回归,比如回归树,随机森林的回归,支持向量回归,贝叶斯回归等等。除此之外,我们还有各种鲁棒的回归:比如RANSAC,Theil-Sen估计,胡贝尔回归等等。考虑到回归问题在现实中的泛用性,回归家族可以说是非常繁荣昌盛,家大业大了。
将全程使用矩阵方式(线性代数的方式)为大家展现回归大家族的面貌。
1.2 sklearn中的线性回归
sklearn中的线性模型模块是linear_model,我们曾经在学习逻辑回归的时候提到过这个模块。linear_model包含了多种多样的类和函数,其中逻辑回归相关的类和函数在这里就不给大家列举了。今天的课中我将会为大家来讲解:普通线性回归,多项式回归,岭回归,LASSO,以及弹性网。
2 多元线性回归LinearRegression
2.1 多元线性回归的基本原理
记得在逻辑回归和SVM中,我们都是先定义了损失函数,然后通过最小化损失函数或损失函数的某种变化来将求解参数向量,以此将单纯的求解问题转化为一个最优化问题。在多元线性回归中,我们的损失函数如下定义:
在L2范式上开平方,就是我们的损失函数。这个式子,也正是sklearn当中,用在类Linear_model.LinerRegression背后使用的损失函数。我们往往称呼这个式子为SSE(Sum of Sqaured Error,误差平方和)或者RSS(Residual Sum of Squares 残差平方和)。在sklearn所有官方文档和网页上,我们都称之为RSS残差平方和,因此在我们的课件中我们也这样称呼。
2.2 最小二乘法求解多元线性回归的参数
现在问题转换成了求解让RSS最小化的参数向量w,这种通过最小化真实值和预测值之间的RSS来求解参数的方法叫做最小二乘法。求解极值的第一步往往是求解一阶导数并让一阶导数等于0,最小二乘法也不能免俗。因此,我们现在残差平方和RSS上对参数向量w求导。这里的过程涉及到少数矩阵求导的内容,需要查表来确定,感兴趣可以走维基百科去查看矩阵求导的详细公式的表格:
除了多元线性回归的推导之外,这里还需要提到一些在上面的推导过程中不曾被体现出来的问题。在统计学中,**使用最小二乘法来求解线性回归的方法是一种”无偏估计“的方法,这种无偏估计要求因变量,也就是标签的分布必须服从正态分布。**这是说,我们的y必须经由正太化处理(比如说取对数,或者使用在第三章《数据预处理与特征工程》中提到的类QuantileTransformer或者PowerTransformer)。在机器学习中,我们会先考虑模型的效果,如果模型效果不好,那我们可能考虑改变因变量的分布(即改变标签的分布)。
2.3 linear_model.LinearRegression
线性回归的类可能是我们目前为止学到的最简单的类,仅有四个参数就可以完成一个完整的算法。并且看得出,这些参数中并没有一个是必填的,更没有对我们的模型有不可替代作用的参数。这说明,线性回归的性能,往往取决于数据本身,而并非是我们的调参能力,线性回归也因此对数据有着很高的要求。幸运的是,现实中大部分连续型变量之间,都存在着或多或少的线性联系。所以线性回归虽然简单,却很强大。
顺便一提,sklearn中的线性回归可以处理多标签问题,只需要在fit的时候输入多维度标签就可以了。
来做一次回归试试看吧:
#1. 导入需要的模块和库
from sklearn.linear_model import LinearRegression as LR #线性回归
from sklearn.model_selection import train_test_split #划分训练集和测试集
from sklearn.model_selection import cross_val_score #交叉验证
from sklearn.datasets import fetch_california_housing as fch #加利福尼亚房屋价值数据集
import pandas as pd
#2. 导入数据,探索数据
housevalue = fch() #会需要下载,大家可以提前运行试试看
housevalue
housevalue.data #该数据集的特征矩阵
X = pd.DataFrame(housevalue.data) #放入DataFrame中便于查看
X
X.shape
y = housevalue.target
y
y.min()
y.max()
#因为是回归这些查看是必须的,因为是连续型的标签y
#此外y并不是房价本身,而是一个分布在0到5之间的连续型评分值
y.shape
X.head()
housevalue.feature_names #这里的feature_names其实就是特征的名字
"""
MedInc:该街区住户的收入中位数
HouseAge:该街区房屋使用年代的中位数
AveRooms:该街区平均的房间数目
AveBedrms:该街区平均的卧室数目
Population:街区人口
AveOccup:平均入住率
Latitude:街区的纬度
Longitude:街区的经度
"""
X.columns
X.columns = housevalue.feature_names
X.head()
#3. 分训练集和测试集
Xtrain, Xtest, Ytrain, Ytest = train_test_split(X,y,test_size=0.3,random_state=420)
Xtrain.shape[0]
range(Xtrain.shape[0])
Xtest.shape[0]
range(Xtest.shape[0])
#恢复索引
for i in [Xtrain, Xtest]:
i.index = range(i.shape[0])
Xtrain.head()
Xtest.head()
Xtrain.shape
#如果要对数据进行标准化,在划分完训练集合测试集后就可以做了。
#如果希望进行数据标准化,还记得应该怎么做吗?
#先用训练集训练(fit)标准化的类,然后用训练好的类分别转化(transform)训练集和测试集
#4. 建模
reg = LR(normalize=True).fit(Xtrain, Ytrain)
#通常情况下希望fit_intercept = Ture(默认为True),即计算截距
yhat = reg.predict(Xtest) #预测yhat
yhat
yhat.min()
yhat.max()
#发现预测的结果是在真实值范围之外的
#说明拟合效果不是特别好
#5. 探索建好的模型
reg.coef_ #w,系数向量
#查看了w,发现系数都很小,这样可以尝试将线性回归中的normalize参数设置为True.
#将normalize=True设置训练后,发现yhat.min()和yhat.max()没有变,此外系数w也没有变化。
#将特征名称与W系数相结合。zip()是一个匹配函数
[*zip(Xtrain.columns,reg.coef_)]
#w系数越大,说明对应的特征越重要。
#w系数越小,说明对应的特征越不重要。
"""
MedInc:该街区住户的收入中位数
HouseAge:该街区房屋使用年代的中位数
AveRooms:该街区平均的房间数目
AveBedrms:该街区平均的卧室数目
Population:街区人口
AveOccup:平均入住率
Latitude:街区的纬度
Longitude:街区的经度
"""
reg.intercept_ #截距,即矩阵中的w0,常数项。
3 回归类的模型评估指标
回归类算法的模型评估一直都是回归算法中的一个难点,但不像我们曾经讲过的无监督学习算法中的轮廓系数等等评估指标,回归类与分类型算法的模型评估其实是相似的法则——找真实标签和预测值的差异。只不过在分类型算法中,这个差异只有一种角度来评判,那就是是否预测到了正确的分类,而在我们的回归类算法中,我们有两种不同的角度来看待回归的效果:
第一,我们是否预测到了正确的数值。
第二,我们是否拟合到了足够的信息。
这两种角度,分别对应着不同的模型评估指标。
3.1 是否预测了正确的数值
回忆一下我们的RSS残差平方和,它的本质是我们的预测值与真实值之间的差异,也就是从第一种角度来评估我们回归的效力,所以RSS既是我们的损失函数,也是我们回归类模型的模型评估指标之一。但是,RSS有着致命的缺点:它是一个无界的和,可以无限地大。我们只知道,我们想要求解最小的RSS,从RSS的公式来看,它不能为负,所以RSS越接近0越好,但我们没有一个概念,究竟多小才算好,多接近0才算好?为了应对这种状况,sklearn中使用RSS的变体,均方误差MSE(mean squared error)来衡量我们的预测值和真实值的差异:
**均方误差,本质是在RSS的基础上除以了样本总量,得到了每个样本量上的平均误差。**有了平均误差,我们就可以将平均误差和我们的标签的取值范围在一起比较,以此获得一个较为可靠的评估依据。在sklearn当中,我们有两种方式调用这个评估指标,一种是使用sklearn专用的模型评估模块metrics里的类mean_squared_error,另一种是调用交叉验证的类cross_val_score并使用里面的scoring参数来设置使用均方误差。
#MSE均方误差有两种方法
#方法1,使用sklearn专用的模型评估模块metrics里的类mean_squared_error
from sklearn.metrics import mean_squared_error as MSE
MSE(yhat,Ytest)
Ytest.mean() #平均预测就错了20%,是一个比较大的值,模型预测效果不是很好。
y.max()
y.min()
#方法2,调用交叉验证的类cross_val_score并使用里面的scoring参数来设置使用均方误差。
cross_val_score(reg,X,y,cv=10,scoring="mean_squared_error")
#cv为10,进行10次交叉验证。
#这里会报错
#为什么报错了?来试试看!
#因为参数scoring对应的是“neg_mean_squared_error”,即负均方根误差
import sklearn
sorted(sklearn.metrics.SCORERS.keys())
cross_val_score(reg,X,y,cv=10,scoring="neg_mean_squared_error")
线性回归的大坑一号:均方误差为负。
#欢迎来的线性回归的大坑一号:均方误差为负。
#我们在决策树和随机森林中都提到过,虽然均方误差永远为正,但是sklearn中的参数scoring下,均方误差作为评判标准时,却是计算”负均方误差“(neg_mean_squared_error)。
#这是因为sklearn在计算模型评估指标的时候,会考虑指标本身的性质,均方误差本身是一种误差,所以被sklearn划分为模型的一种损失(loss)。
#在sklearn当中,所有的损失都使用负数表示,因此均方误差也被显示为负数了。真正的均方误差MSE的数值,其实就是neg_mean_squared_error去掉负号的数字。
cross_val_score(reg,X,y,cv=10,scoring="neg_mean_squared_error").mean()
#将10个负的均方根误差一起求得一个均值
#方法1的mse为0.5309012639324575
#方法2的mse为0.55095242969566
#二者较为接近。
除了MSE,我们还有与MSE类似的MAE(Mean absolute error,绝对均值误差):
其表达的概念与均方误差完全一致,不过在真实标签和预测值之间的差异外我们使用的是L1范式(绝对值)。**现实使用中,MSE和MAE选一个来使用就好了。**在sklearn当中,我们使用命令from sklearn.metrics import mean_absolute_error来调用MAE,同时,我们也可以使用交叉验证中的scoring = “neg_mean_absolute_error”,以此在交叉验证时调用MAE。
3.2 是否拟合了足够的信息
对于回归类算法而言,只探索数据预测是否准确是不足够的。除了数据本身的数值大小之外,我们还希望我们的模型能够捕捉到数据的”规律“,比如数据的分布规律,单调性等等,而是否捕获了这些信息并无法使用MSE来衡量。
来看这张图,其中红色线是我们的真实标签,而蓝色线是我们的拟合模型。这是一种比较极端,但的确可能发生的情况。这张图像上,前半部分的拟合非常成功,看上去我们的真实标签和我们的预测结果几乎重合,但后半部分的拟合却非常糟糕,模型向着与真实标签完全相反的方向去了。对于这样的一个拟合模型,如果我们使用MSE来对它进行判断,它的MSE会很小,因为大部分样本其实都被完美拟合了,少数样本的真实值和预测值的巨大差异在被均分到每个样本上之后,MSE就会很小。但这样的拟合结果必然不是一个好结果,因为一旦我的新样本是处于拟合曲线的后半段的,我的预测结果必然会有巨大的偏差,而这不是我们希望看到的。所以,我们希望找到新的指标,除了判断预测的数值是否正确之外,还能够判断我们的模型是否拟合了足够多的,数值之外的信息。
R^2可以使用三种方式来调用:
一种是直接从metrics中导入r2_score,输入预测值和真实值后打分。
第二种是直接从线性回归LinearRegression的接口score来进行调用。
第三种是在交叉验证中,输入"r2"来调用。
from sklearn.linear_model import LinearRegression as LR #线性回归
from sklearn.model_selection import train_test_split #划分训练集和测试集
from sklearn.model_selection import cross_val_score #交叉验证
from sklearn.datasets import fetch_california_housing as fch #加利福尼亚房屋价值数据集
import pandas as pd
#导入数据
housevalue = fch() #会需要下载,大家可以提前运行试试看
X = pd.DataFrame(housevalue.data) #放入DataFrame中便于查看
X.columns = housevalue.feature_names
y = housevalue.target
Xtrain, Xtest, Ytrain, Ytest = train_test_split(X,y,test_size=0.3,random_state=420)
#恢复索引
for i in [Xtrain, Xtest]:
i.index = range(i.shape[0])
reg = LR(normalize=True).fit(Xtrain, Ytrain)
#通常情况下希望fit_intercept = Ture(默认为True),即计算截距
yhat = reg.predict(Xtest) #预测yhat
yhat
#第一种:调用R2
from sklearn.metrics import r2_score
r2_score(yhat,Ytest)
#R^2越接近1越好
#第二种:
#线性回归LinearRegression的接口score来进行调用
#输入的参数是Xtest和Ytest
r2 = reg.score(Xtest,Ytest)
r2
#发现与metrics调出的R^2不一样
线性回归的大坑二号:相同的评估指标不同的结果。
#这是因为,在我们的分类模型的评价指标当中,我们进行的是一种 if a == b的对比,这种判断和if b == a其实完全是一种概念,所以我们在进行模型评估的时候,从未踩到我们现在在的这个坑里。
#然而看R2的计算公式,R2明显和分类模型的指标中的accuracy或者precision不一样,R2涉及到的计算中对预测值和真实值有极大的区别,必须是预测值在分子,真实值在分母。
#所以我们在调用metrcis模块中的模型评估指标的时候,必须要检查清楚,指标的参数中,究竟是要求我们先输入真实值还是先输入预测值,之前r2_score中的参数就输反了。
#使用shift tab键来检查究竟哪个值先进行输入
#查看后,先y_true,后y_pred,
r2_score(Ytest,yhat)
#或者你也可以指定参数,就不必在意顺序了
r2_score(y_true = Ytest,y_pred = yhat)
#第三种:
#交叉验证R^2的均值结果也不是很理想
cross_val_score(reg,X,y,cv=10,scoring="r2").mean()
#我们观察到,我们在加利福尼亚房屋价值数据集上的MSE其实不是一个很大的数(0.5),但我们的r2不高,这证明我们的模型比较好地拟合了一部分数据的数值,却没有能正确拟合数据的分布。
#让我们与绘图来看看,究竟是不是这样一回事。我们可以绘制一张图上的两条曲线,一条曲线是我们的真实标签Ytest,另一条曲线是我们的预测结果yhat,两条曲线的交叠越多,我们的模型拟合就越好。
import matplotlib.pyplot as plt
sorted(Ytest) #从小到大排序
#为什么要排序:test是胡乱分布在X轴上的,X轴是有顺序的,所以需要从小到大排序
plt.plot(range(len(Ytest)),sorted(Ytest),c="black",label= "Data")
plt.plot(range(len(yhat)),sorted(yhat),c="red",label = "Predict")
plt.legend()
plt.show()
#R^2是衡量回归类模型评估指标的第一评估指标。
#如果R^2不行,MSE均方根误差越小,该模型也是不行的。
#可见,虽然我们的大部分数据被拟合得比较好,但是图像的开头和结尾处却又着较大的拟合误差。
#如果我们在图像右侧分布着更多的数据,我们的模型就会越来越偏离我们真正的标签。
#这种结果类似于我们前面提到的,虽然在有限的数据集上将数值预测正确了,但却没有正确拟合数据的分布,如果有更多的数据进入我们的模型,那数据标签被预测错误的可能性是非常大的。
了线性回归的三号大坑:负的R^2
#现在查看一组有趣的情况
#使用能够numpy自己创建一组x和一组y
import numpy as np
rng = np.random.RandomState(42) #随机数,参数42表示随机数种子
X = rng.randn(100, 80)
y = rng.randn(100)
cross_val_score(LR(), X, y, cv=5, scoring='r2')
#当我们的r2显示为负的时候,这证明我们的模型对我们的数据的拟合非常糟糕,模型完全不能使用。
#所有,一个负的r2是合理的。
#当然了,现实应用中,如果你发现你的线性回归模型出现了负的r2,不代表你就要接受他了,首先检查你的建模过程和数据处理过程是否正确,也许你已经伤害了数据本身,也许你的建模过程是存在bug的。
#如果是集成模型的回归,检查你的弱评估器的数量是否不足,随机森林,提升树这些模型在只有两三棵树的时候很容易出现负的r2。
#如果你检查了所有的代码,也确定了你的预处理没有问题,但你的r2也还是负的,那这就证明,线性回归模型不适合你的数据,试试看其他的算法吧。
4 多重共线性:岭回归与Lasso
4.1 最熟悉的陌生人:多重共线性
在第二节中我们曾推导了多元线性回归使用最小二乘法的求解原理,我们对多元线性回归的损失函数求导,并得出求解系数 的式子和过程:
在最后一步中我们需要左乘(X^T)*X的逆矩阵,而逆矩阵存在的充分必要条件是特征矩阵不存在多重共线性。
多重共线性这个词对于许多人来说都不陌生,然而却很少有人能够透彻理解这个性质究竟是什么含义,会有怎样的影响。这一节我们会来深入讲解,什么是多重共线性,我们是如何一步步从逆矩阵必须存在推导到多重共线性不能存在的。本节需要大量的线性代数知识作为支撑,讲解会比较细致,如果依然感觉对其中的数学原理理解困难,建议以专门讲线性代数的数学教材为辅学习本节。
逆矩阵存在的充分必要条件:
首先,我们需要先理解逆矩阵存在与否的意义和影响。一个矩阵什么情况下才可以有逆矩阵呢?
来看逆矩阵的计算公式:
分字上A*是伴随矩阵,任何矩阵都可以有伴随矩阵,因此这一部分不影响逆矩阵的存在性。而分母上的行列式 |A| 就不同了,位于分母的变量不能为0,一旦为0则无法计算出逆矩阵。**因此逆矩阵存在的充分必要条件是:矩阵的行列式不能为0,对于线性回归而言,即是说 |(X^T)*X| 不能为0。**这是使用最小二乘法来求解线性回归的核心条件之一。
行列式不为0的充分必要条件:
那行列式要不为0,需要满足什么条件呢?在这里,我们来复习一下线性代数中的基本知识。
假设我们的特征矩阵X结构为(m,n),则(X^T)X就是结构为(n,m)的矩阵乘以结构为(m,n)的矩阵,从而得到结果为(n,n)的方阵。
因此以下所有的例子都将以方阵进行举例,方便大家理解。首先区别一下矩阵和行列式:
重要定义:矩阵和行列式
任何矩阵都可以有行列式。以一个33的行列式为例,我们来看看行列式是如何计算的:
这个式子乍一看非常混乱,其实并非如此,我们把行列式 |A| 按照下面的方式排列一下,就很容易看出这个式子实际上是怎么回事了:
三个特征的特征矩阵的行列式就有六个交互项,在现实中我们的特征矩阵不可能是如此低维度的数据,因此使用这样的方式计算行列式就变得异常困难。在线性代数中,我们可以通过行列式的计算将一个行列式整合成一个梯形的行列式:
行列式变换只能是选择行变换或列变换,不能即变换行同时也变换列。
我们可以对矩阵做初等行变换和列变换,包括交换行/列顺序,将一列/一行乘以一个常数后加减到另一列/一行上,来将矩阵化为梯形矩阵。
梯形的行列式表现为,所有的数字都被整合到对角线的上方或下方(通常是上方),虽然具体的数字发生了变化(比如由 变成了 ),但是行列式的大小在初等行变换/列变换的过程中是不变的。对于梯形行列式,行列式的计算要容易很多:
不难发现,由于梯形行列式下半部分为0,整个矩阵的行列式其实就是梯形行列式对角线上的元素相乘。并且此时此刻,只要对角线上的任意元素为0,整个行列式都会为0。那只要对角线上没有一个元素为0,行列式就不会为0了。在这里,我们来引入一个重要的概念:满秩矩阵。重要定义:满秩矩阵
满秩矩阵:A是一个n行n列的矩阵,若A转换为梯形矩阵后,没有任何全为0的行或者全为0的列,则称A为满秩矩阵。简单来说,只要对角线上没有一个元素为0,则这个矩阵中绝对不可能存在全为0的行或列。
举例来说,下面的矩阵就不是满秩矩阵,因为它的对角线上有一个0,因此它存在全为0的行。
即是说,矩阵满秩(即转换为梯形矩阵后对角线上没有0)是矩阵的行列式不为0的充分必要条件。矩阵满秩的充分必要条件:
一个矩阵要满秩,则转换为梯形矩阵后的对角线上没有0,那什么样的矩阵在转换为梯形矩阵后对角线上才没有0呢?
来看下面的三个矩阵:
我们可以对矩阵做初等行变换和列变换,包括交换行/列顺序,将一列/一行乘以一个常数后加减到另一列/一行上,来将矩阵化为梯形矩阵。对于上面的两个矩阵我们可以有如下变换:
如此就转换成了梯形矩阵。我们可以看到,矩阵A明显不是满秩的,它有全零行所以行列式会为0。而矩阵B和C没有全零行所以满秩。而矩阵A和矩阵B的区别在于,A中存在着完全具有线性关系的两行(1,1,2和2,2,4),而B和C中则没有这样的两行。而矩阵B虽然对角线上每个元素都不为0,但具有非常接近于0的元素0.02,而矩阵C的对角线上没有任何元素特别接近于0。矩阵A中第一行和第三行的关系,被称为**“精确相关关系”,即完全相关**,一行可使另一行为0。在这种精确相关关系下,矩阵A的行列式为0,则矩阵A的逆不可能存在。在我们的最小二乘法中,如果矩阵(X^T)*X中存在这种精确相关关系,则逆不存在,最小二乘法完全无法使用,线性回归会无法求出结果。
矩阵B中第一行和第三行的关系不太一样,他们之间非常接近于”精确相关关系“,但又不是完全相关,一行不能使另一行为0,这种关系被称为”高度相关关系“。在这种高度相关关系下,矩阵的行列式不为0,但是一个非常接近0数,矩阵A的逆存在,不过接近于无限大。在这种情况下,最小二乘法可以使用,不过得到的逆会很大,直接影响我们对参数向量w的求解:
这样求解出来的参数向量 会很大,因此会影响建模的结果,造成模型有偏差或者不可用。**精确相关关系和高度相关关系并称为"多重共线性"。**在多重共线性下,模型无法建立,或者模型不可用。相对的,矩阵C的行之间结果相互独立,梯形矩阵看起来非常正常,它的对角线上没有任何元素特别接近于0,因此其行列式也就不会接近0或者为0,因此矩阵C得出的参数向量w就不会有太大偏差,对于我们拟合而言是比较理想的。
从上面的所有过程我们可以看得出来,一个矩阵如果要满秩,则要求矩阵中每个向量之间不能存在多重共线性,这也构成了线性回归算法对于特征矩阵的要求。
多重共线性与相关性:
多重共线性如果存在,则线性回归就无法使用最小二乘法来进行求解,或者求解就会出现偏差。幸运的是,不能存在多重共线性,不代表不能存在相关性——机器学习不要求特征之间必须独立,必须不相关,只要不是高度相关或者精确相关就好。
多重共线性 Multicollinearity 与 相关性 Correlation:
多重共线性是一种统计现象,是指线性模型中的特征(解释变量)之间由于存在精确相关关系或高度相关关系。多重共线性的存在会使模型无法建立,或者估计失真。
多重共线性使用指标方差膨胀因子(variance inflation factor,VIF)来进行衡量(from statsmodels.stats.outliers_influence import variance_inflation_factor),VIF大于10就认为存在严重多重共线性,严格一些大于5也认为存在多重共线性,通常当我们提到“共线性”,都特指多重共线性。
相关性是衡量两个或多个变量一起波动的程度的指标,它可以是正的,负的或者0。当我们说变量之间具有相关性,通常是指线性相关性,线性相关一般由皮尔逊相关系数进行衡量,非线性相关可以使用斯皮尔曼相关系数或者互信息法进行衡量。
在现实中特征之间完全独立的情况其实非常少,因为大部分数据统计手段或者收集者并不考虑统计学或者机器学习建模时的需求,现实数据多多少少都会存在一些相关性,极端情况下,甚至还可能出现收集的特征数量比样本数量多的情况。**通常来说,这些相关性在机器学习中通常无伤大雅(在统计学中他们可能是比较严重的问题),即便有一些偏差,只要最小二乘法能够求解,我们都有可能会无视掉它。**毕竟,想要消除特征的相关性,无论使用怎样的手段,都无法避免进行特征选择,这意味着可用的信息变得更加少,对于机器学习来说,很有可能尽量排除相关性后,模型的整体效果会受到巨大的打击。这种情况下,我们选择不处理相关性,只要结果好,一切万事大吉。
然而多重共线性就不是这样一回事了,它的存在会造成模型极大地偏移,无法模拟数据的全貌,因此这是必须解决的问题。
为了保留线性模型计算快速,理解容易的优点,我们并不希望更换成非线性模型,这促使统计学家和机器学习研究者们钻研出了多种能够处理多重共线性的方法,其中有三种比较常见的:
这三种手段中,第一种相对耗时耗力,需要较多的人工操作,并且会需要混合各种统计学中的知识和检验来进行使用。在机器学习中,能够使用一种模型解决的问题,我们尽量不用多个模型来解决,如果能够追求结果,我们会尽量避免进行一系列检验。况且,统计学中的检验往往以“让特征独立”为目标,与机器学习中的”稍微有点相关性也无妨“不太一致。第二种手段在现实中应用较多,不过由于理论复杂,效果也不是非常高效,因此向前逐步回归不是机器学习的首选。
我们的核心会是使用第三种方法:改进线性回归来处理多重共线性。为此,一系列算法,岭回归,Lasso,弹性网就被研究出来了。接下来,我们就来看看这些改善多重共线性问题的算法。
4.2 岭回归
4.2.1 岭回归解决多重共线性问题
在线性模型之中,除了线性回归之外,最知名的就是岭回归与Lasso了。这两个算法非常神秘,他们的原理和应用都不像其他算法那样高调,学习资料也很少。**这可能是因为这两个算法不是为了提升模型表现,而是为了修复漏洞而设计的(实际上,我们使用岭回归或者Lasso,模型的效果往往会下降一些,因为我们删除了一小部分信息),因此在结果为上的机器学习领域颇有些被冷落的意味。**这一节我们就来了解一下岭回归。
岭回归,又称为吉洪诺夫正则化(Tikhonov regularization)。通常来说,大部分的机器学习教材会使用代数的形式来展现岭回归的原理,这个原理和逻辑回归及支持向量机非常相似,都是将求解 的过程转化为一个带条件的最优化问题,然后用最小二乘法求解。然而,岭回归可以做到的事其实可以用矩阵非常简单地表达出来。
岭回归在多元线性回归的损失函数上加上了正则项,表达为系数 的L2范式(即系数 的平方项)乘以正则化系数。如果你们看其他教材中的代数推导,正则化系数会写作 ,用以和Lasso区别,不过在sklearn中由于是两个不同的算法,因此正则项系数都使用 来代表。岭回归的损失函数的完整表达式写作:
这个操作看起来简单,其实带来了巨大的变化。在线性回归中我们通过在损失函数上对w求导来求解极值,在这里虽然加上了正则项,我们依然使用最小二乘法来求解。假设我们的特征矩阵结构为(m,n),系数w的结构是(1,n),则可以有:
最后得到的这个行列式还是一个梯形行列式,然而它的已经不存在全0行或者全0列了,除非:
重要:
4.2.2 linear_model.Ridge 岭回归
在sklearn中,岭回归由线性模型库中的Ridge类来调用:
和线性回归相比,岭回归的参数稍微多了那么一点点,但是真正核心的参数就是我们的正则项的系数alpha,其他的参数是当我们希望使用最小二乘法之外的求解方法求解岭回归的时候才需要的,通常我们完全不会去触碰这些参数。所以大家只需要了解alpha的用法就可以了。
之前我们在加利佛尼亚房屋价值数据集上使用线性回归,得出的结果大概是训练集上的拟合程度是60%,测试集上的拟合程度也是60%左右,那这个很低的拟合程度是不是由多重共线性造成的呢?在统计学中,我们会通过VIF或者各种检验来判断数据是否存在共线性,然而在机器学习中,我们可以使用模型来判断——如果一个数据集在岭回归中使用各种正则化参数取值下模型表现没有明显上升(比如出现持平或者下降),则说明数据没有多重共线性,顶多是特征之间有一些相关性。反之,如果一个数据集在岭回归的各种正则化参数取值下表现出明显的上升趋势,则说明数据存在多重共线性。
接下来,我们就在加利佛尼亚房屋价值数据集上来验证一下这个说法:
import numpy as np
import pandas as pd
from sklearn.model_selection import cross_val_score #交叉验证
from sklearn.linear_model import Ridge, LinearRegression, Lasso
from sklearn.model_selection import train_test_split as TTS
from sklearn.datasets import fetch_california_housing as fch
import matplotlib.pyplot as plt
housevalue = fch()
X = pd.DataFrame(housevalue.data)
y = housevalue.target
X.columns = ["住户收入中位数","房屋使用年代中位数","平均房间数目"
,"平均卧室数目","街区人口","平均入住率","街区的纬度","街区的经度"]
#这里最好不要使用中文,因为中文可能会报错
X.head()
Xtrain,Xtest,Ytrain,Ytest = TTS(X,y,test_size=0.3,random_state=420)
#数据集索引恢复
for i in [Xtrain,Xtest]:
i.index = range(i.shape[0])
Xtrain
#使用岭回归来进行建模
reg = Ridge(alpha=1).fit(Xtrain,Ytrain)
reg.score(Xtest,Ytest)
#加利福尼亚房屋价值数据集中应该不是共线性问题
#因为加上了正则项alpha,就是用来处理共线性的问题,分数应该有一定的提升,但是分数并没有明显变化。
#交叉验证下,与线性回归相比,岭回归的结果如何变化?
alpharange = np.arange(1,1001,100)
ridge, lr = [], []
for alpha in alpharange:
reg = Ridge(alpha=alpha)
linear = LinearRegression()
regs = cross_val_score(reg,X,y,cv=5,scoring = "r2").mean()
linears = cross_val_score(linear,X,y,cv=5,scoring = "r2").mean()
ridge.append(regs)
lr.append(linears)
plt.plot(alpharange,ridge,color="red",label="Ridge")
plt.plot(alpharange,lr,color="orange",label="LR")
plt.title("Mean")
plt.legend()
plt.show()
#细化一下学习曲线
alpharange = np.arange(1,201,10)
ridge, lr = [], []
for alpha in alpharange:
reg = Ridge(alpha=alpha)
linear = LinearRegression()
regs = cross_val_score(reg,X,y,cv=5,scoring = "r2").mean()
linears = cross_val_score(linear,X,y,cv=5,scoring = "r2").mean()
ridge.append(regs)
lr.append(linears)
plt.plot(alpharange,ridge,color="red",label="Ridge")
plt.plot(alpharange,lr,color="orange",label="LR")
plt.title("Mean")
plt.legend()
plt.show()
#通过观察图片,说明加利福尼亚州房屋价值数据存在轻微的共线性
#R^2只上升了还不到0.001
#R^2越小预测与真实偏差越大,即越小越差
#如果想要验证多重共线性,可以使用VIF,VIF可以判断哪些特征是存在多重共线性的。
#可以看出,加利佛尼亚数据集上,岭回归的结果轻微上升,随后骤降。可以说,加利佛尼亚房屋价值数据集带有很轻微的一部分共线性,这种共线性被正则化参数alpha消除后,模型的效果提升了一点点,但是对于整个模型而言是杯水车薪。
#在过了控制多重共线性的点后,模型的效果飞速下降,显然是正则化的程度太重,挤占了参数w本来的估计空间。从这个结果可以看出,加利佛尼亚数据集的核心问题不在于多重共线性,岭回归不能够提升模型表现。
#在正则化参数逐渐增大的过程中,我们可以观察一下模型的方差如何变化:
#模型方差如何变化?
alpharange = np.arange(1,1001,100)
ridge, lr = [], []
for alpha in alpharange:
reg = Ridge(alpha=alpha)
linear = LinearRegression()
varR = cross_val_score(reg,X,y,cv=5,scoring="r2").var()
varLR = cross_val_score(linear,X,y,cv=5,scoring="r2").var()
ridge.append(varR)
lr.append(varLR)
plt.plot(alpharange,ridge,color="red",label="Ridge")
plt.plot(alpharange,lr,color="orange",label="LR")
plt.title("Variance")
plt.legend()
plt.show()
#随着这个方差的急速上升,R^2的急速下降,
#证明随着正则化系数的逐渐增大,模型的泛化能力逐渐降低。
#数目没有alpha的时候,模型的泛化能力还比较好。
#即alpha比较小时,模型在测试集上的表现更好。
#我们调整出的泛化能力点是错误的,因为alpha越小,分数反而比我们调整出来的高。
可以发现,模型的方差上升快速,不过方差的值本身很小,其变化不超过 上升部分的1/3,因此只要噪声的状况维持恒定,模型的泛化误差可能还是一定程度上降低了的。虽然岭回归和Lasso不是设计来提升模型表现,而是专注于解决多重共线性问题的,但当 在一定范围内变动的时候,消除多重共线性也许能够一定程度上提高模型的泛化能力。
但是泛化能力毕竟没有直接衡量的指标,因此我们往往只能够通过观察模型的准确性指标和方差来大致评判模型的泛化能力是否提高。来看看多重共线性更为明显一些的情况:
#这里使用波士顿房价数据集
from sklearn.datasets import load_boston
from sklearn.model_selection import cross_val_score
X = load_boston().data
X.shape
y = load_boston().target
Xtrain,Xtest,Ytrain,Ytest = TTS(X,y,test_size=0.3,random_state=420)
#先查看方差的变化
alpharange = np.arange(1,1001,100)
ridge, lr = [], []
for alpha in alpharange:
reg = Ridge(alpha=alpha)
linear = LinearRegression()
varR = cross_val_score(reg,X,y,cv=5,scoring="r2").var()
varLR = cross_val_score(linear,X,y,cv=5,scoring="r2").var()
ridge.append(varR)
lr.append(varLR)
plt.plot(alpharange,ridge,color="red",label="Ridge")
plt.plot(alpharange,lr,color="orange",label="LR")
plt.title("Variance")
plt.legend()
plt.show()
#通过观察图片,岭回归虽则alpha不断的变大,方差逐渐减小。
#这个说明波士顿房价数据是存在多重共线性的。
#否则alpha为什么会对方差有这么大的影响。
#查看R2的变化
alpharange = np.arange(1,1001,100)
ridge, lr = [], []
for alpha in alpharange:
reg = Ridge(alpha=alpha)
linear = LinearRegression()
regs = cross_val_score(reg,X,y,cv=5,scoring = "r2").mean()
linears = cross_val_score(linear,X,y,cv=5,scoring = "r2").mean()
ridge.append(regs)
lr.append(linears)
plt.plot(alpharange,ridge,color="red",label="Ridge")
plt.plot(alpharange,lr,color="orange",label="LR")
plt.title("Mean")
plt.legend()
plt.show()
#细化学习曲线
alpharange = np.arange(100,300,10)
ridge, lr = [], []
for alpha in alpharange:
reg = Ridge(alpha=alpha)
#linear = LinearRegression()
regs = cross_val_score(reg,X,y,cv=5,scoring = "r2").mean()
#linears = cross_val_score(linear,X,y,cv=5,scoring = "r2").mean()
ridge.append(regs)
lr.append(linears)
plt.plot(alpharange,ridge,color="red",label="Ridge")
#plt.plot(alpharange,lr,color="orange",label="LR")
plt.title("Mean")
plt.legend()
plt.show()
可以发现,比起加利佛尼亚房屋价值数据集,波士顿房价数据集的方差降低明显,偏差也降低明显,可见使用岭回归还是起到了一定的作用,模型的泛化能力是有可能会上升的。
遗憾的是,没有人会希望自己获取的数据中存在多重共线性,因此发布到scikit-learn或者kaggle上的数据基本都经过一定的多重共线性的处理的,要找出绝对具有多重共线性的数据非常困难,也就无法给大家展示岭回归在实际数据中大显身手的模样。我们也许可以找出具有一些相关性的数据,但是大家如果去尝试就会发现,基本上如果我们使用岭回归或者Lasso,那模型的效果都是会降低的,很难升高,这恐怕也是岭回归和Lasso一定程度上被机器学习领域冷遇的原因。
4.2.3 为岭回归选取最佳的正则化参数取值
既然要选择alpha的范围,我们就不可避免地要进行最优参数的选择。在各种机器学习教材中,总是教导使用岭迹图来判断正则项参数的最佳取值。传统的岭迹图长这样,形似一个开口的喇叭图(根据横坐标的正负,喇叭有可能朝右或者朝左):
这一个以正则化参数为横坐标,线性模型求解的系数w为纵坐标的图像,其中每一条彩色的线都是一个系数w。其目标是建立正则化参数与系数w之间的直接关系,以此来观察正则化参数的变化如何影响了系数w的拟合。岭迹图认为,线条交叉越多,则说明特征之间的多重共线性越高。我们应该选择系数较为平稳的喇叭口所对应的alpha取值作为最佳的正则化参数的取值。绘制岭迹图的方法非常简单,代码如下:
import numpy as np
import matplotlib.pyplot as plt
from sklearn import linear_model
#创造10*10的希尔伯特矩阵
X = 1. / (np.arange(1, 11) + np.arange(0, 10)[:, np.newaxis])
y = np.ones(10)
#计算横坐标
n_alphas = 200
alphas = np.logspace(-10, -2, n_alphas)
#建模,获取每一个正则化取值下的系数组合
coefs = []
for a in alphas:
ridge = linear_model.Ridge(alpha=a, fit_intercept=False)
ridge.fit(X, y)
coefs.append(ridge.coef_)
#绘图展示结果
ax = plt.gca()
ax.plot(alphas, coefs)
ax.set_xscale('log')
ax.set_xlim(ax.get_xlim()[::-1]) #将横坐标逆转
plt.xlabel('正则化参数alpha')
plt.ylabel('系数w')
plt.title('岭回归下的岭迹图')
plt.axis('tight')
plt.show()
其中涉及到的希尔伯特矩阵长这样:
然而,我非常不建议大家使用岭迹图来作为寻找最佳参数的标准。
有这样的两个理由:
1、岭迹图的很多细节,很难以解释。比如为什么多重共线性存在会使得线与线之间有很多交点?当 很大了之后看上去所有的系数都很接近于0,难道不是那时候线之间的交点最多吗?
2、岭迹图的评判标准,非常模糊。哪里才是最佳的喇叭口?哪里才是所谓的系数开始变得”平稳“的时候?一千个读者一千个哈姆雷特的画像?未免也太不严谨了。所以我们应该使用交叉验证来选择最佳的正则化系数。在sklearn中,我们有带交叉验证的岭回归可以使用,我们来看一看:
可以看到,这个类于普通的岭回归类Ridge非常相似,不过在输入正则化系数alpha的时候我们可以传入元祖作为正则化系数的备选,非常类似于我们在画学习曲线前设定的for i in 的列表对象。来看RidgeCV的重要参数,属性和接口:
这个类的使用也非常容易,依然使用我们之前建立的加利佛尼亚房屋价值数据集:
import numpy as np
import pandas as pd
from sklearn.linear_model import RidgeCV, LinearRegression
from sklearn.model_selection import train_test_split as TTS
from sklearn.datasets import fetch_california_housing as fch #加利佛尼亚房屋价值数据集
import matplotlib.pyplot as plt
housevalue = fch()
X = pd.DataFrame(housevalue.data)
y = housevalue.target
X.columns = ["住户收入中位数","房屋使用年代中位数","平均房间数目"
,"平均卧室数目","街区人口","平均入住率","街区的纬度","街区的经度"]
X
#参数scoring和cv都是默认值:
#交叉验证是将全部数据导入的,所以就不用分训练集和测试集了
Ridge_ = RidgeCV(alphas=np.arange(1,1001,100)
#,scoring="neg_mean_squared_error" #scoring默认是R^2作为交叉验证的默认评估指标。
,store_cv_values=True
#,cv=5 #cv默认是none,即留一交叉验证,只有一个样本是做测试集,其他所以样本都是训练集。
).fit(X, y) #实例化和训练
#无关交叉验证的岭回归结果
Ridge_.score(X,y)
#调用所有交叉验证的结果
Ridge_.cv_values_ #二维数组
#查看结构,20640个样本,10表示10个不同alpha取值下交叉验证的结果。
#因为是留一交叉验证,所以是20640个样本,每一个样本所对应的10个交叉验证结果。
Ridge_.cv_values_.shape
#进行平均后可以查看每个正则化系数取值下的交叉验证结果
#这10个是不同alpha取值下的交叉验证结果
Ridge_.cv_values_.mean(axis=0)
#查看被选择出来的最佳正则化系数
Ridge_.alpha_
#为什么是101,因为下图是当最高R^2下的alpha值:
#为什么是R^2,因为R^2作为岭回归交叉验证的默认评估指标。
#交叉验证下,与线性回归相比,岭回归的结果如何变化?
alpharange = np.arange(1,1001,100)
ridge, lr = [], []
for alpha in alpharange:
reg = Ridge(alpha=alpha)
linear = LinearRegression()
regs = cross_val_score(reg,X,y,cv=5,scoring = "r2").mean()
linears = cross_val_score(linear,X,y,cv=5,scoring = "r2").mean()
ridge.append(regs)
lr.append(linears)
plt.plot(alpharange,ridge,color="red",label="Ridge")
plt.plot(alpharange,lr,color="orange",label="LR")
plt.title("Mean")
plt.legend()
plt.show()
#说明alpha=101时模型效果最好
#这里将scoring为负的均方根误差来试一下:
#交叉验证是将全部数据导入的,所以就不用分训练集和测试集了
Ridge_ = RidgeCV(alphas=np.arange(1,1001,100)
,scoring="neg_mean_squared_error" #scoring默认是R^2作为交叉验证的默认评估指标。
,store_cv_values=True
#,cv=5 #cv默认是none,即留一交叉验证,只有一个样本是做测试集,其他所以样本都是训练集。
).fit(X, y) #实例化和训练
#无关交叉验证的岭回归结果
Ridge_.score(X,y) #依然返回R^2
#调用所有交叉验证的结果
Ridge_.cv_values_.shape #二维数组
#进行平均后可以查看每个正则化系数取值下的交叉验证结果
#这10个是不同alpha取值下的交叉验证结果
Ridge_.cv_values_.mean(axis=0)
#这些结果都是很小的数,已经不是R^2了。这些结果是均方误差
#查看被选择出来的最佳正则化系数
Ridge_.alpha_
#发现依然选择了101
#这说明scoring为"neg_mean_squared_error"与scoring默认是R^2作为交叉验证的默认评估指标的结果是相同的,都选择了101
#这里将cv修改为5试一下:
#交叉验证是将全部数据导入的,所以就不用分训练集和测试集了
Ridge_ = RidgeCV(alphas=np.arange(1,1001,100)
#,scoring="neg_mean_squared_error" #scoring默认是R^2作为交叉验证的默认评估指标。
,store_cv_values=True
,cv=5 #cv默认是none,即留一交叉验证,只有一个样本是做测试集,其他所以样本都是训练集。
).fit(X, y) #实例化和训练
#发现报错了,因为cv参数输入数字(即cv不在是默认值)
#RidgeCV是不会为我们保留交叉验证的结果的。
#所以store_cv_values不能为Ture,在cv不是默认值的情况下。
#交叉验证是将全部数据导入的,所以就不用分训练集和测试集了
Ridge_ = RidgeCV(alphas=np.arange(1,1001,100)
#,scoring="neg_mean_squared_error" #scoring默认是R^2作为交叉验证的默认评估指标。
,store_cv_values=False
,cv=5 #cv默认是none,即留一交叉验证,只有一个样本是做测试集,其他所以样本都是训练集。
).fit(X, y) #实例化和训练
#无论上面参数如何调整,该值都不会变,因为与交叉验证无关。
#无关交叉验证的岭回归结果
Ridge_.score(X,y) #依然返回R^2
#调用所有交叉验证的结果
#报错是因为cv不是默认值,不会保留交叉验证的结果。
Ridge_.cv_values_
#这个也自然无法运行出来结果了,因为cv不是默认值,不会保留交叉验证的结果。
#进行平均后可以查看每个正则化系数取值下的交叉验证结果
#这10个是不同alpha取值下的交叉验证结果
Ridge_.cv_values_.mean(axis=0)
#这些结果都是很小的数,已经不是R^2了。这些结果是均方误差
#查看被选择出来的最佳正则化系数
Ridge_.alpha_
#发现alpha依然是101
这说明模型在选择最佳alpha时,都会根据评估指标最高的那个值来进行选择。不论参数如何调整,都会选择出来最佳正则化系数alpha(即参数调整不影响选择出最佳正则化系数Ridge_.alpha_),参数调整会运行Ridge_.cv_values_的结果,不过是不会展示趋势的。
4.3 Lasso
4.3.1 Lasso与多重共线性
除了岭回归之外,最常被人们提到还有模型Lasso。Lasso全称最小绝对收缩和选择算子(least absolute shrinkage and selection operator),由于这个名字过于复杂所以简称为Lasso。和岭回归一样,Lasso是被创造来作用于多重共线性问题的算法,不过Lasso使用的是系数w的L1范式(L1范式则是系数w的绝对值)乘以正则化系数alpha,所以Lasso的损失函数表达式为:
许多博客和机器学习教材会说,Lasso与岭回归非常相似,都是利用正则项来对原本的损失函数形成一个惩罚,以此来防止多重共线性。这种说法不是非常严谨,我们来看看Lasso的数学过程。当我们使用最小二乘法来求解Lasso中的参数w,我们依然对损失函数进行求导:
岭回归 vs Lasso:
岭回归可以解决特征间的精确相关关系导致的最小二乘法无法使用的问题,而Lasso不行。
重要:
通过增大alpha,我们可以为w的计算增加一个负项,从而限制参数估计中w的大小,而防止多重共线性引起的参数w被估计过大导致模型失准的问题。**Lasso不是从根本上解决多重共线性问题,而是限制多重共线性带来的影响。**何况,这还是在我们假设所有的系数都为正的情况下,假设系数w无法为正,则很有可能我们需要将我们的正则项参数alpha设定为负,因此alpha可以取负数,并且负数越大,对共线性的限制也越大。
所有这些让Lasso成为了一个神奇的算法,尽管它是为了限制多重共线性被创造出来的,然而世人其实并不使用它来抑制多重共线性,反而接受了它在其他方面的优势。我们在讲解逻辑回归时曾提到过,L1和L2正则化一个核心差异就是他们对系数w的影响:两个正则化都会压缩系数w的大小,对标签贡献更少的特征的系数会更小,也会更容易被压缩。不过,L2正则化只会将系数压缩到尽量接近0,但L1正则化主导稀疏性,因此会将系数压缩到0。这个性质,让Lasso成为了线性模型中的特征选择工具首选,接下来,我们就来看看如何使用Lasso来选择特征。
4.3.2 Lasso的核心作用:特征选择
sklearn中我们使用类Lasso来调用lasso回归,众多参数中我们需要比较在意的就是参数alpha,正则化系数。另外需要注意的就是参数positive。当这个参数为"True"的时候,是我们要求Lasso回归出的系数必须为正数,以此来保证我们的alpha一定以增大来控制正则化的程度。需要注意的是,在sklearn中我们的Lasso使用的损失函数是:
接下来,我们就来看看lasso如何做特征选择:
import numpy as np
import pandas as pd
from sklearn.linear_model import Ridge, LinearRegression, Lasso
from sklearn.model_selection import train_test_split as TTS
from sklearn.datasets import fetch_california_housing as fch
import matplotlib.pyplot as plt
housevalue = fch()
X = pd.DataFrame(housevalue.data)
y = housevalue.target
X.columns = ["住户收入中位数","房屋使用年代中位数","平均房间数目"
,"平均卧室数目","街区人口","平均入住率","街区的纬度","街区的经度"]
X.head()
Xtrain,Xtest,Ytrain,Ytest = TTS(X,y,test_size=0.3,random_state=420)
#划分了训练集和测试集后一定要恢复索引
for i in [Xtrain,Xtest]:
i.index = range(i.shape[0])
#对比线性回归、岭回归和Lasso
#线性回归进行拟合
reg = LinearRegression().fit(Xtrain,Ytrain)
(reg.coef_*100).tolist()
#reg.coef_的结果是系数w值,原本系数太小,所以乘以100
#岭回归进行拟合
Ridge_ = Ridge(alpha=0).fit(Xtrain,Ytrain)
(Ridge_.coef_*100).tolist()
#Lasso进行拟合
lasso_ = Lasso(alpha=0).fit(Xtrain,Ytrain)
(lasso_.coef_*100).tolist()
#可以对比这三个拟合的结果,发现几乎没有什么区别,只是最后几位数组不同。
#此外Lasso出现了三个红色警告
#这三条分别是这样的内容:
#1. 正则化系数为0,这样算法不可收敛!如果你想让正则化系数为0,请使用线性回归吧
#2. 没有正则项的坐标下降法可能会导致意外的结果,不鼓励这样做!
#3. 目标函数没有收敛,你也许想要增加迭代次数,使用一个非常小的alpha来拟合模型可能会造成精确度问题!
#看到这三条内容,大家可能会比较懵——怎么出现了坐标下降?这是由于sklearn中的Lasso类不是使用最小二乘法来进行求解,而是使用坐标下降。考虑一下,Lasso既然不能够从根本解决多重共线性引起的最小二乘法无法使用的问题,那我们为什么要坚持最小二乘法呢?明明有其他更快更好的求解方法,比如坐标下降就很好呀。Lasso不是从根本上解决多重共线性问题,而是限制多重共线性带来的影响。
#如果我们的确希望取到0,那我们可以使用一个比较很小的数,比如0.01,或者10*e^-3这样的值
#Lasson就不会报错了
#岭回归进行拟合
Ridge_ = Ridge(alpha=0.01).fit(Xtrain,Ytrain)
(Ridge_.coef_*100).tolist()
#Lasso进行拟合
lasso_ = Lasso(alpha=0.01).fit(Xtrain,Ytrain)
(lasso_.coef_*100).tolist()
#可以发现改变alpha,Lasso的值就会产生较大变化,而岭回归产生变化较小,Lasso对alpha特别敏感
#加大正则项系数,观察模型的系数发生了什么变化
Ridge_ = Ridge(alpha=10**4).fit(Xtrain,Ytrain)
(Ridge_.coef_*100).tolist()
#岭回归会将系数w压缩到接近于0,但始终不会像lasso那样直接压缩等于0.
lasso_ = Lasso(alpha=10**4).fit(Xtrain,Ytrain)
(lasso_.coef_*100).tolist()
#看来10**4对于Lasso来说是一个过于大的取值
lasso_ = Lasso(alpha=1).fit(Xtrain,Ytrain)
(lasso_.coef_*100).tolist()
#再次说明Lasso对alpha很敏感
#可以将特征系数w压缩到0
#将系数进行绘图
plt.plot(range(1,9),(reg.coef_*100).tolist(),color="red",label="LR")
plt.plot(range(1,9),(Ridge_.coef_*100).tolist(),color="orange",label="Ridge")
plt.plot(range(1,9),(lasso_.coef_*100).tolist(),color="k",label="Lasso")
plt.plot(range(1,9),[0]*8,color="grey",linestyle="--")
plt.xlabel('w') #横坐标是每一个特征所对应的系数
plt.legend()
plt.show()
#该图横坐标是八个系数,纵坐标是系数所对应的w的值。
#对于线性回归,w系数该是多少就是多少,w系数值该是0就是0,该是多少就是多少。,不存在压缩这一说法。
#对于岭回归,可以看到3和4已经压缩到了很小的位置,但始终不会到0。
#对于lasso,除了1号w以外,其他w都别压缩到0了。
可见,比起岭回归,Lasso所带的L1正则项对于系数的惩罚要重得多,并且它会将系数压缩至0,因此可以被用来做特征选择。
也因此,我们往往让Lasso的正则化系数 在很小的空间中变动,以此来寻找最佳的正则化系数。
使用Lasso进行特征选择,寻找某一个alpha,让w系数压缩为0,压缩为0的w系数所对应的特征,都可以舍弃掉,
只要找到一个适当的alpha,就可以去掉w系数为0所对应的特征,这样就看让模型轻量化。
4.3.3 选取最佳的正则化参数取值
使用交叉验证的Lasso类的参数看起来与岭回归略有不同,这是由于Lasso对于alpha的取值更加敏感的性质决定的。之前提到过,由于Lasso对正则化系数的变动过于敏感,因此我们往往让 在很小的空间中变动。这个小空间小到超乎人们的想象(不是0.01到0.02之间这样的空间,这个空间对lasso而言还是太大了),因此我们设定了一个重要概念“正则化路径”,用来设定正则化系数的变动:重要概念:正则化路径 regularization path
和岭回归的交叉验证类相似,除了进行交叉验证之外,LassoCV也会单独建立模型。它会先找出最佳的正则化参数,然后在这个参数下按照模型评估指标进行建模。需要注意的是,LassoCV的模型评估指标选用的是均方误差,而岭回归的模型评估指标是可以自己设定的,并且默认是R^2。
来看看将这些参数和属性付诸实践的代码:
#Lasso选取最佳的正则化参数取值
from sklearn.linear_model import LassoCV
from sklearn.model_selection import cross_val_score #交叉验证
import numpy as np
import pandas as pd
from sklearn.linear_model import Ridge, LinearRegression, Lasso
from sklearn.model_selection import train_test_split as TTS
from sklearn.datasets import fetch_california_housing as fch
import matplotlib.pyplot as plt
housevalue = fch()
X = pd.DataFrame(housevalue.data)
y = housevalue.target
X.columns = ["住户收入中位数","房屋使用年代中位数","平均房间数目"
,"平均卧室数目","街区人口","平均入住率","街区的纬度","街区的经度"]
X.head()
Xtrain,Xtest,Ytrain,Ytest = TTS(X,y,test_size=0.3,random_state=420)
#划分了训练集和测试集后一定要恢复索引
for i in [Xtrain,Xtest]:
i.index = range(i.shape[0])
#自己建立Lasso进行alpha选择的范围
#np.logspace()是用来建立指数函数的,一共四个参数:
#参数1:起始点1.00000000e-10;
#参数2:终止点1.00000000e-02;
#参数3:起始点到终止点之间想要选取多少个值;
#参数4:表示对数的底数。base=10表示对数以10为底。
#其实是形成10为底的指数函数
#10**(-10)到10**(-2)次方
alpharange = np.logspace(-10, -2, 200,base=10)
alpharange
alpharange.shape
Xtrain.head()
lasso_ = LassoCV(alphas=alpharange #自行输入的alpha的取值范围
,cv=5 #交叉验证的折数
).fit(Xtrain, Ytrain)
#查看被选择出来的最佳正则化系数
lasso_.alpha_
#调用所有交叉验证的结果
lasso_.mse_path_
lasso_.mse_path_.shape #返回每个alpha下的五折交叉验证结果
lasso_.mse_path_.mean(axis=1) #有注意到在岭回归中我们的轴向是axis=0吗?
#注意:(重要)
#在岭回归当中,我们是留一验证,因此我们的交叉验证结果返回的是,每一个样本在每个alpha下的交叉验证结果
#因此我们要求每个alpha下的交叉验证均值,就是axis=0,跨行求均值
#而在Lasso这里,我们返回的是,每一个alpha取值下,每一折交叉验证的结果
#因此我们要求每个alpha下的交叉验证均值,就是axis=1,跨列求均值
lasso_.mse_path_.mean(axis=1).shape
#最佳正则化系数下获得的模型的系数结果
lasso_.coef_
#将测试数据导入Lasso模型
lasso_.score(Xtest,Ytest)
#与线性回归相比如何?
reg = LinearRegression().fit(Xtrain,Ytrain)
reg.score(Xtest,Ytest)
#使用lassoCV自带的正则化路径长度和路径中的alpha个数来自动建立alpha选择的范围
#eps表示正则化路径长度,是alpha的最小值除以alpha的最大值的比例
#n_alphas表示300个alpha的取值
#cv表示模型的折数
ls_ = LassoCV(eps=0.00001
,n_alphas=300
,cv=5
).fit(Xtrain, Ytrain)
ls_.alpha_
ls_.alphas_ #查看所有自动生成的alpha取值
#应该有300个
ls_.alphas_.shape
ls_.score(Xtest,Ytest)
#没有好于线性回归,线性回归本身模型比较好
ls_.coef_
#依然没有任何特征的系数w被压缩到0
到这里,岭回归和Lasso的核心作用就为大家讲解完毕了。时间缘故无法为大家将坐标下降法展开来解释,Lasso作为线性回归家族中在改良上走得最远的算法,还有许多领域等待我们去探讨。比如说,在现实中,我们不仅可以适用交叉验证来选择最佳正则化系数,我们也可以使用BIC( 贝叶斯信息准则)或者AIC(Akaike information criterion,艾凯克信息准则)来做模型选择。同时,我们可以不使用坐标下降法,还可以使用最小角度回归来对lasso进行计算。
当然了,这些方法下做的模型选择和模型计算,其实在模型效果上表现和普通的Lasso没有太大的区别,不过他们都在各个方面对原有的Lasso做了一些相应的改进(比如说提升了本来就已经很快的计算速度,增加了模型选择的维度,因为均方误差作为损失函数只考虑了偏差,不考虑方差的存在)。
除了解决多重共线性这个核心问题之外,线性模型还有更重要的事情要做:提升模型表现。这才是机器学习最核心的需求,而Lasso和岭回归不是为此而设计的。下一节,让我们来认识一下为了提升模型表现而做出的改进:多项式回归。