Python机器学习-线性回归

酥酥 发布于 2022-04-01 956 次阅读


1 概述
1.1 线性回归大家族
回归是一种应用广泛的预测建模技术,这种技术的核心在于预测的结果是连续型变量。决策树,随机森林,支持向量机的分类器等分类算法的预测标签是分类变量,多以{0,1}来表示,而无监督学习算法比如PCA,KMeans并不求解标签,注意加以区别。回归算法源于统计学理论,它可能是机器学习算法中产生最早的算法之一,其在现实中的应用非常广泛,包括使用其他经济指标预测股票市场指数,根据喷射流的特征预测区域内的降水量,根据公司的广告花费预测总销售额,或者根据有机物质中残留的碳-14的量来估计化石的年龄等等,只要一切基于特征预测连续型变量的需求,我们都使用回归技术。


既然线性回归是源于统计分析,是结合机器学习与统计学的重要算法。通常来说,我们认为统计学注重先验,而机器学习看重结果,因此机器学习中不会提前为线性回归排除共线性等可能会影响模型的因素,反而会先建立模型以查看效果。模型确立之后,如果效果不好,我们就根据统计学的指导来排除可能影响模型的因素。我们的课程会从机器学习的角度来为大家讲解回归类算法,如果希望理解统计学角度的小伙伴们,各种统计学教材都可以满足你的需求。


回归需求在现实中非常多,所以我们自然也有各种各样的回归类算法。最著名的就是我们的线性回归和逻辑回归,从他们衍生出了岭回归,Lasso,弹性网,除此之外,还有众多分类算法改进后的回归,比如回归树,随机森林的回归,支持向量回归,贝叶斯回归等等。除此之外,我们还有各种鲁棒的回归:比如RANSAC,Theil-Sen估计,胡贝尔回归等等。考虑到回归问题在现实中的泛用性,回归家族可以说是非常繁荣昌盛,家大业大了。


回归类算法的数学相对简单,相信在经历了逻辑回归,主成分分析与奇异值分解,支持向量机这三个章节之后,大家不会感觉到线性回归中的数学有多么困难。通常,理解线性回归可以有两种角度:矩阵的角度和代数的角度。几乎所有机器学习的教材都是从代数的角度来理解线性回归的,类似于我们在逻辑回归和支持向量机中做的那样,将求解参数的问题转化为一个带条件的最优化问题,然后使用三维图像让大家理解求极值的过程。如果大家掌握了逻辑回归和支持向量机,这个过程可以说是相对简单的,因此我们在本章节中就不进行赘述了。相对的,在我们的课程中一直都缺乏比较系统地使用矩阵来解读算法的角度,因此在本堂课中,我将全程使用矩阵方式(线性代数的方式)为大家展现回归大家族的面貌。


学完这堂课之后,大家需要对线性模型有个相对全面的了解,尤其是需要掌握线性模型究竟存在什么样的优点和问题,并且如何解决这些问题。

1.2 sklearn中的线性回归
sklearn中的线性模型模块是linear_model,我们曾经在学习逻辑回归的时候提到过这个模块。linear_model包含了多种多样的类和函数,其中逻辑回归相关的类和函数在这里就不给大家列举了。今天的课中我将会为大家来讲解:普通线性回归,多项式回归,岭回归,LASSO,以及弹性网。

2.3 linear_model.LinearRegression
class sklearn.linear_model.LinearRegression (fit_intercept=True, normalize=False, copy_X=True, n_jobs=None)

线性回归的类可能是我们目前为止学到的最简单的类,仅有四个参数就可以完成一个完整的算法。并且看得出,这些参数中并没有一个是必填的,更没有对我们的模型有不可替代作用的参数。这说明,线性回归的性能,往往取决于数据本身,而并非是我们的调参能力,线性回归也因此对数据有着很高的要求。幸运的是,现实中大部分连续型变量之间,都存在着或多或少的线性联系。所以线性回归虽然简单,却很强大。
顺便一提,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() #会需要下载,大家可以提前运行试试看
X = pd.DataFrame(housevalue.data) #放入DataFrame中便于查看
y = housevalue.target
X.shape
y.shape
X.head()
housevalue.feature_names
X.columns = housevalue.feature_names
"""
MedInc:该街区住户的收入中位数
HouseAge:该街区房屋使用年代的中位数
AveRooms:该街区平均的房间数目
AveBedrms:该街区平均的卧室数目
Population:街区人口
AveOccup:平均入住率
Latitude:街区的纬度
Longitude:街区的经度
"""
菜菜的sklearn课堂直播间: https://live.bilibili.com/12582510
sklearn专题第九期:回归大家族
				
			

3. 分训练集和测试集

				
					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])
Xtrain.shape
#如果希望进行数据标准化,还记得应该怎么做吗?
#先用训练集训练标准化的类,然后用训练好的类分别转化训练集和测试集
				
			

4. 建模

				
					reg = LR().fit(Xtrain, Ytrain)
yhat = reg.predict(Xtest)
yhat
				
			

5. 探索建好的模型

				
					reg.coef_
[*zip(Xtrain.columns,reg.coef_)]
"""
MedInc:该街区住户的收入中位数
HouseAge:该街区房屋使用年代的中位数
AveRooms:该街区平均的房间数目
AveBedrms:该街区平均的卧室数目
Population:街区人口
AveOccup:平均入住率
Latitude:街区的纬度
Longitude:街区的经度
"""
reg.intercept_
				
			

建模的过程在sklearn当中其实非常简单,但模型的效果如何呢?接下来我们来看看多元线性回归的模型评估指标。

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参数来设置使用均方误差。

				
					from sklearn.metrics import mean_squared_error as MSE
MSE(yhat,Ytest)
y.max()
y.min()
cross_val_score(reg,X,y,cv=10,scoring="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去掉负号的数字。除了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来衡量。

				
					#调用R2
from sklearn.metrics import r2_score
r2_score(yhat,Ytest)
r2 = reg.score(Xtest,Ytest)
r2
				
			

我们现在踩到了线性回归的大坑二号:相同的评估指标不同的结果。
为什么结果会不一致呢?这就是回归和分类算法的不同带来的坑。
在我们的分类模型的评价指标当中,我们进行的是一种 if a == b的对比,这种判断和if b == a其实完全是一种概念,所以我们在进行模型评估的时候,从未踩到我们现在在的这个坑里。然而看R2的计算公式,R2明显和分类模型的指标中的accuracy或者precision不一样,R2涉及到的计算中对预测值和真实值有极大的区别,必须是预测值在分子,真实值在分母,所以我们在调用metrcis模块中的模型评估指标的时候,必须要检查清楚,指标的参数中,究竟是要求我们先输入真实值还是先输入预测值。

				
					#使用shift tab键来检查究竟哪个值先进行输入
r2_score(Ytest,yhat)
#或者你也可以指定参数,就不必在意顺序了
r2_score(y_true = Ytest,y_pred = yhat)
cross_val_score(reg,X,y,cv=10,scoring="r2").mean()
				
			

我们观察到,我们在加利福尼亚房屋价值数据集上的MSE其实不是一个很大的数(0.5),但我们的 不高,这证明我们的模型比较好地拟合了一部分数据的数值,却没有能正确拟合数据的分布。让我们与绘图来看看,究竟是不是这样一回事。我们可以绘制一张图上的两条曲线,一条曲线是我们的真实标签Ytest,另一条曲线是我们的预测结果yhat,两条曲线的交叠越多,我们的模型拟合就越好。

				
					import matplotlib.pyplot as plt
sorted(Ytest)
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()
				
			

可见,虽然我们的大部分数据被拟合得比较好,但是图像的开头和结尾处却又着较大的拟合误差。如果我们在图像右侧分布着更多的数据,我们的模型就会越来越偏离我们真正的标签。这种结果类似于我们前面提到的,虽然在有限的数据集上将数值预测正确了,但却没有正确拟合数据的分布,如果有更多的数据进入我们的模型,那数据标签被预测错误的可能性是非常大的。

思考
在sklearn中,一个与 非常相似的指标叫做可解释性方差分数(explained_variance_score,EVS),它也是衡量 1 – 没有捕获到的信息占总信息的比例,但它和 略有不同。虽然在实践应用中EVS应用不多,但感兴趣的小伙伴可以探索一下这个回归类模型衡量指标。现在,来看一组有趣的情况:

				
					import numpy as np
rng = np.random.RandomState(42)
X = rng.randn(100, 80)
y = rng.randn(100)
cross_val_score(LR(), X, y, cv=5, scoring='r2'
				
			

多重共线性与相关性
多重共线性如果存在,则线性回归就无法使用最小二乘法来进行求解,或者求解就会出现偏差。幸运的是,不能存在多重共线性,不代表不能存在相关性——机器学习不要求特征之间必须独立,必须不相关,只要不是高度相关或者精确相关就好。

多重共线性 Multicollinearity 与 相关性 Correlation
多重共线性是一种统计现象,是指线性模型中的特征(解释变量)之间由于存在精确相关关系或高度相关关系,多重共线性的存在会使模型无法建立,或者估计失真。多重共线性使用指标方差膨胀因子(variance inflationfactor,VIF)来进行衡量(from statsmodels.stats.outliers_influence import variance_inflation_factor),通常当我们提到“共线性”,都特指多重共线性。相关性是衡量两个或多个变量一起波动的程度的指标,它可以是正的,负的或者0。当我们说变量之间具有相关性,通常是指线性相关性,线性相关一般由皮尔逊相关系数进行衡量,非线性相关可以使用斯皮尔曼相关系数或者互信息法进行衡量。

在现实中特征之间完全独立的情况其实非常少,因为大部分数据统计手段或者收集者并不考虑统计学或者机器学习建模时的需求,现实数据多多少少都会存在一些相关性,极端情况下,甚至还可能出现收集的特征数量比样本数量多的情况。通常来说,这些相关性在机器学习中通常无伤大雅(在统计学中他们可能是比较严重的问题),即便有一些偏差,只要最小二乘法能够求解,我们都有可能会无视掉它。毕竟,想要消除特征的相关性,无论使用怎样的手段,都无法避免进行特征选择,这意味着可用的信息变得更加少,对于机器学习来说,很有可能尽量排除相关性后,模型的整体效果会受到巨大的打击。这种情况下,我们选择不处理相关性,只要结果好,一切万事大吉。
然而多重共线性就不是这样一回事了,它的存在会造成模型极大地偏移,无法模拟数据的全貌,因此这是必须解决的问题。为了保留线性模型计算快速,理解容易的优点,我们并不希望更换成非线性模型,这促使统计学家和机器学习研究者们钻研出了多种能够处理多重共线性的方法,其中有三种比较常见的:

这三种手段中,第一种相对耗时耗力,需要较多的人工操作,并且会需要混合各种统计学中的知识和检验来进行使用。在机器学习中,能够使用一种模型解决的问题,我们尽量不用多个模型来解决,如果能够追求结果,我们会尽量避免进行一系列检验。况且,统计学中的检验往往以“让特征独立”为目标,与机器学习中的”稍微有点相关性也无妨“不太一致。

第二种手段在现实中应用较多,不过由于理论复杂,效果也不是非常高效,因此向前逐步回归不是机器学习的首选。在本周的课程中,我们的核心会是使用第三种方法:改进线性回归来处理多重共线性。为此,一系列算法,岭回归,Lasso,弹性网就被研究出来了。接下来,我们就来看看这些改善多重共线性问题的算法。
4.2 岭回归
4.2.1 岭回归解决多重共线性问题
在线性模型之中,除了线性回归之外,最知名的就是岭回归与Lasso了。这两个算法非常神秘,他们的原理和应用都不像其他算法那样高调,学习资料也很少。这可能是因为这两个算法不是为了提升模型表现,而是为了修复漏洞而设计的(实际上,我们使用岭回归或者Lasso,模型的效果往往会下降一些,因为我们删除了一小部分信息),因此在结果为上的机器学习领域颇有些被冷落的意味。这一节我们就来了解一下岭回归。

4.2.2 linear_model.Ridge
在sklearn中,岭回归由线性模型库中的Ridge类来调用:
class sklearn.linear_model.Ridge (alpha=1.0, fit_intercept=True, normalize=False, copy_X=True, max_iter=None,
tol=0.001, solver=’auto’, random_state=None)
和线性回归相比,岭回归的参数稍微多了那么一点点,但是真正核心的参数就是我们的正则项的系数 ,其他的参数是当我们希望使用最小二乘法之外的求解方法求解岭回归的时候才需要的,通常我们完全不会去触碰这些参数。所以大家只需要了解 的用法就可以了。
之前我们在加利佛尼亚房屋价值数据集上使用线性回归,得出的结果大概是训练集上的拟合程度是60%,测试集上的拟合程度也是60%左右,那这个很低的拟合程度是不是由多重共线性造成的呢?在统计学中,我们会通过VIF或者各种检验来判断数据是否存在共线性,然而在机器学习中,我们可以使用模型来判断——如果一个数据集在岭回归中使用各种正则化参数取值下模型表现没有明显上升(比如出现持平或者下降),则说明数据没有多重共线性,顶多是特征之间有一些相关性。反之,如果一个数据集在岭回归的各种正则化参数取值下表现出明显的上升趋势,则说明数据存在多重共线性。
接下来,我们就在加利佛尼亚房屋价值数据集上来验证一下这个说法:

				
					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])
#使用岭回归来进行建模
reg = Ridge(alpha=1).fit(Xtrain,Ytrain)
reg.score(Xtest,Ytest)
#交叉验证下,与线性回归相比,岭回归的结果如何变化?
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)
				
			

可以看出,加利佛尼亚数据集上,岭回归的结果轻微上升,随后骤降。可以说,加利佛尼亚房屋价值数据集带有很轻微的一部分共线性,这种共线性被正则化参数 消除后,模型的效果提升了一点点,但是对于整个模型而言是杯水车薪。在过了控制多重共线性的点后,模型的效果飞速下降,显然是正则化的程度太重,挤占了参数 本来的估计空间。从这个结果可以看出,加利佛尼亚数据集的核心问题不在于多重共线性,岭回归不能够提升模型表现。另外,在正则化参数逐渐增大的过程中,我们可以观察一下模型的方差如何变化:

				
					#模型方差如何变化?
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()
				
			

可以发现,模型的方差上升快速,不过方差的值本身很小,其变化不超过 上升部分的1/3,因此只要噪声的状况维持恒定,模型的泛化误差可能还是一定程度上降低了的。虽然岭回归和Lasso不是设计来提升模型表现,而是专注于解决多重共线性问题的,但当 在一定范围内变动的时候,消除多重共线性也许能够一定程度上提高模型的泛化能力。
但是泛化能力毕竟没有直接衡量的指标,因此我们往往只能够通过观察模型的准确性指标和方差来大致评判模型的泛化能力是否提高。来看看多重共线性更为明显一些的情况:

				
					from sklearn.datasets import load_boston
from sklearn.model_selection import cross_val_score
X = load_boston().data
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()
#查看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 选取最佳的正则化参数取值
既然要选择a的范围,我们就不可避免地要进行最优参数的选择。在各种机器学习教材中,总是教导使用岭迹图来判断正则项参数的最佳取值。传统的岭迹图长这样,形似一个开口的喇叭图(根据横坐标的正负,喇叭有可能朝右或者朝左):

这一个以正则化参数为横坐标,线性模型求解的系数w为纵坐标的图像,其中每一条彩色的线都是一个系数w 。其目标是建立正则化参数与系数w之间的直接关系,以此来观察正则化参数的变化如何影响了系数 的拟合。岭迹图认为,线条交叉越多,则说明特征之间的多重共线性越高。我们应该选择系数较为平稳的喇叭口所对应的 a取值作为最佳的正则化参数的取值。绘制岭迹图的方法非常简单,代码如下:

				
					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()
				
			

难得一提的机器学习发展历史:过时的岭迹图
其实,岭迹图会有这样的问题不难理解。岭回归最初始由Hoerl和Kennard在1970提出来用来改进多重共线性问题的模型,在这片1970年的论文中,两位作者提出了岭迹图并且向广大学者推荐这种方法,然而遭到了许多人的批评和反抗。大家接受了岭回归,却鲜少接受岭迹图,这使得岭回归被发明了50年之后,市面上关于岭迹图的教材依然只有1970年的论文中写的那几句话。
1974年,Stone M发表论文,表示应当在统计学和机器学习中使用交叉验证。1980年代,机器学习技术迎来第一次全面爆发(1979年ID3决策树被发明出来,1980年之后CART树,adaboost,带软间隔的支持向量,梯度提升树逐渐诞生),从那之后,除了统计学家们,几乎没有人再使用岭迹图了。在2000年以后,岭迹图只是教学中会被略微提到的一个知识点(还会被强调是过时的技术),在现实中,真正应用来选择正则化系数的技术是交叉验证,并且选择的标准非常明确——我们选择让交叉验证下的均方误差最小的正则化系数 a。

				
					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 = ["住户收入中位数","房屋使用年代中位数","平均房间数目"
     ,"平均卧室数目","街区人口","平均入住率","街区的纬度","街区的经度"]
Ridge_ = RidgeCV(alphas=np.arange(1,1001,100)
        #,scoring="neg_mean_squared_error"
        ,store_cv_values=True
        #,cv=5
       ).fit(X, y)
       #无关交叉验证的岭回归结果
Ridge_.score(X,y)
#调用所有交叉验证的结果
Ridge_.cv_values_.shape
#进行平均后可以查看每个正则化系数取值下的交叉验证结果
Ridge_.cv_values_.mean(axis=0)
#查看被选择出来的最
				
			
				
					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])
#线性回归进行拟合
reg = LinearRegression().fit(Xtrain,Ytrain)
(reg.coef_*100).tolist()
#岭回归进行拟合
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既然不能够从根本解决多重共线性引起的最小二乘法无法使用的问题,那我们为什么要坚持最小二乘法呢?明明有其他更快更好的求解方法,比如坐标下降就很好呀。下面两篇论文解释了了scikit-learn坐标下降求解器中使用的迭代方式,以及用于收敛控制的对偶间隙计算方式,感兴趣的大家可以进行阅读。

使用坐标下降法求解Lasso
“Regularization Path For Generalized linear Models by Coordinate Descent”, Friedman, Hastie &Tibshirani, J Stat Softw, 2010 ( Paper).
“An Interior-Point Method for Large-Scale L1-Regularized Least Squares,” S. J. Kim, K. Koh, M. Lustig, S.Boyd and D. Gorinevsky, in IEEE Journal of Selected Topics in Signal Processing, 2007 (Paper)

				
					#岭回归进行拟合
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()
				
			

这样就不会报任何错误了。

				
					#加大正则项系数,观察模型的系数发生了什么变化
Ridge_ = Ridge(alpha=10**4).fit(Xtrain,Ytrain)
(Ridge_.coef_*100).tolist()
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()
#将系数进行绘图
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()
				
			

可见,比起岭回归,Lasso所带的L1正则项对于系数的惩罚要重得多,并且它会将系数压缩至0,因此可以被用来做特征选择。也因此,我们往往让Lasso的正则化系数 在很小的空间中变动,以此来寻找最佳的正则化系数。


4.3.3 选取最佳的正则化参数取值
class sklearn.linear_model.LassoCV (eps=0.001, n_alphas=100, alphas=None, fit_intercept=True,normalize=False, precompute=’auto’, max_iter=1000, tol=0.0001, copy_X=True, cv=’warn’, verbose=False,n_jobs=None, positive=False, random_state=None, selection=’cyclic’)

				
					from sklearn.linear_model import LassoCV
#自己建立Lasso进行alpha选择的范围
alpharange = np.logspace(-10, -2, 200,base=10)
#其实是形成10为底的指数函数
#10**(-10)到10**(-2)次方
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,跨行求均值
#而在这里,我们返回的是,每一个alpha取值下,每一折交叉验证的结果
#因此我们要求每个alpha下的交叉验证均值,就是axis=1,跨列求均值
#最佳正则化系数下获得的模型的系数结果
lasso_.coef_
lasso_.score(Xtest,Ytest)
#与线性回归相比如何?
reg = LinearRegression().fit(Xtrain,Ytrain)
reg.score(Xtest,Ytest)
#使用lassoCV自带的正则化路径长度和路径中的alpha个数来自动建立alpha选择的范围
ls_ = LassoCV(eps=0.00001
      ,n_alphas=300
      ,cv=5
       ).fit(Xtrain, Ytrain)
ls_.alpha_
ls_.alphas_ #查看所有自动生成的alpha取值
ls_.alphas_.shape
ls_.score(Xtest,Ytest)
ls_.coef_
				
			

到这里,岭回归和Lasso的核心作用就为大家讲解完毕了。时间缘故无法为大家将坐标下降法展开来解释,Lasso作为线性回归家族中在改良上走得最远的算法,还有许多领域等待我们去探讨。比如说,在现实中,我们不仅可以适用交叉验证来选择最佳正则化系数,我们也可以使用BIC( 贝叶斯信息准则)或者AIC(Akaike informationcriterion,艾凯克信息准则)来做模型选择。同时,我们可以不使用坐标下降法,还可以使用最小角度回归来对lasso进行计算。当然了,这些方法下做的模型选择和模型计算,其实在模型效果上表现和普通的Lasso没有太大的区别,不过他们都在各个方面对原有的Lasso做了一些相应的改进(比如说提升了本来就已经很快的计算速度,增加了模型选择的维度,因为均方误差作为损失函数只考虑了偏差,不考虑方差的存在)。除了解决多重共线性这个核心问题之外,线性模型还有更重要的事情要做:提升模型表现。这才是机器学习最核心的需求,而Lasso和岭回归不是为此而设计的。下一节,让我们来认识一下为了提升模型表现而做出的改进:多项式回归。


5 非线性问题:多项式回归
5.1 重塑我们心中的“线性”概念
在机器学习和统计学中,甚至在我们之前的课程中,我们无数次提到”线性“这个名词。首先我们本周的算法就叫做”线性回归“,而在支持向量机中,我们也曾经提到最初的支持向量机只能够分割线性可分的数据,然后引入了”核函数“来帮助我们分类那些非线性可分的数据。我们也曾经说起过,比如说决策树,支持向量机是”非线性“模型。所有的这些概念,让我们对”线性“这个词非常熟悉,却又非常陌生——因为我们并不知道它的真实含义。在这一小节,我将来为大家重塑线性的概念,并且为大家解决线性回归模型改进的核心之一:帮助线性回归解决非线性问题。

5.1.2 数据的线性与非线性
从线性关系这个概念出发,我们有了一种说法叫做“线性数据”。通常来说,一组数据由多个特征和标签组成。当这些特征分别与标签存在线性关系的时候,我们就说这一组数据是线性数据。当特征矩阵中任意一个特征与标签之间的关系需要使用三角函数,指数函数等函数来定义,则我们就说这种数据叫做“非线性数据”。对于线性和非线性数据,最简单的判别方法就是利用模型来帮助我们——如果是做分类则使用逻辑回归,如果做回归则使用线性回归,如果效果好那数据是线性的,效果不好则数据不是线性的。当然,也可以降维后进行绘图,绘制出的图像分布接近一条直线,则数据就是线性的。

可以看得出,这些数据都不能由一条直线来进行拟合,他们也没有均匀分布在某一条线的周围,那我们怎么判断,这些数据是线性数据还是非线性数据呢?在这里就要注意了,当我们在回归中绘制图像时,绘制的是特征与标签的关系图,横坐标是特征,纵坐标是标签,我们的标签是连续型的,所以我们可以通过是否能够使用一条直线来拟合图像判断数据究竟属于线性还是非线性。然而在分类中,我们绘制的是数据分布图,横坐标是其中一个特征,纵坐标是另一个特征,标签则是数据点的颜色。因此在分类数据中,我们使用“是否线性可分”(linearly separable)这个概念来划分分类数据集。当分类数据的分布上可以使用一条直线来将两类数据分开时,我们就说数据是线性可分的。反之,数据不是线性可分的。


总结一下,对于回归问题,数据若能分布为一条直线,则是线性的,否则是非线性。对于分类问题,数据分布若能使
用一条直线来划分类别,则是线性可分的,否则数据则是线性不可分的。

那线性回归在非线性数据上的表现如何呢?我们来建立一个明显是非线性的数据集,并观察线性回归和决策树的而回归在拟合非线性数据集时的表现:
1. 导入所需要的库

				
					import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.tree import DecisionTreeRegressor
				
			

2. 创建需要拟合的数据集

				
					rnd = np.random.RandomState(42) #设置随机数种子
X = rnd.uniform(-3, 3, size=100) #random.uniform,从输入的任意两个整数中取出size个随机数
#生成y的思路:先使用NumPy中的函数生成一个sin函数图像,然后再人为添加噪音
y = np.sin(X) + rnd.normal(size=len(X)) / 3 #random.normal,生成size个服从正态分布的随机数
#使用散点图观察建立的数据集是什么样子
plt.scatter(X, y,marker='o',c='k',s=20)
plt.show()
#为后续建模做准备:sklearn只接受二维以上数组作为特征矩阵的输入
X.shape
X = X.reshape(-1, 1)
				
			

3. 使用原始数据进行建模

				
					#使用原始数据进行建模
LinearR = LinearRegression().fit(X, y)
TreeR = DecisionTreeRegressor(random_state=0).fit(X, y)
#放置画布
fig, ax1 = plt.subplots(1)
#创建测试数据:一系列分布在横坐标上的点
line = np.linspace(-3, 3, 1000, endpoint=False).reshape(-1, 1)
#将测试数据带入predict接口,获得模型的拟合效果并进行绘制
ax1.plot(line, LinearR.predict(line), linewidth=2, color='green',
    label="linear regression")
ax1.plot(line, TreeR.predict(line), linewidth=2, color='red',
    label="decision tree")
#将原数据上的拟合绘制在图像上
ax1.plot(X[:, 0], y, 'o', c='k')
#其他图形选项
ax1.legend(loc="best")
ax1.set_ylabel("Regression output")
ax1.set_xlabel("Input feature")
ax1.set_title("Result before discretization")
plt.tight_layout()
plt.show()
#从这个图像来看,可以得出什么结果?
				
			

线性模型拟合非线性数据
但是相反的,线性模型若用来拟合非线性数据或者对非线性可分的数据进行分类,那通常都会表现糟糕。通常如果我们已经发现数据属于非线性数据,或者数据非线性可分的数据,则我们不会选择使用线性模型来进行建模。改善线性模型在非线性数据上的效果的方法之一时进行分箱,并且从下图来看分箱的效果不是一般的好,甚至高过一些非线性模型。在下一节中我们会详细来讲解分箱的效果,但很容易注意到,在没有其他算法或者预处理帮忙的情况下,线性模型在非线性数据上的表现时很糟糕的。

模型在线性和非线性数据集上的表现为我们选择模型提供了一个思路:当我们获取数据时,我们往往希望使用线性模型来对数据进行最初的拟合(线性回归用于回归,逻辑回归用于分类),如果线性模型表现良好,则说明数据本身很可能是线性的或者线性可分的,如果线性模型表现糟糕,那毫无疑问我们会投入决策树,随机森林这些模型的怀抱,就不必浪费时间在线性模型上了。
不过这并不代表着我们就完全不能使用线性模型来处理非线性数据了。在现实中,线性模型有着不可替代的优势:计算速度异常快速,所以也还是存在着我们无论如何也希望使用线性回归的情况。因此,我们有多种手段来处理线性回归无法拟合非线性问题的问题,接下来我们就来看一看。


5.2 使用分箱处理非线性问题
让线性回归在非线性数据上表现提升的核心方法之一是对数据进行分箱,也就是离散化。与线性回归相比,我们常用的一种回归是决策树的回归。我们之前拟合过一条带有噪音的正弦曲线以展示多元线性回归与决策树的效用差异,我们来分析一下这张图,然后再使用采取措施帮助我们的线性回归。
1. 导入所需要的库

				
					import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.tree import DecisionTreeRegressor
				
			

2. 创建需要拟合的数据集

				
					rnd = np.random.RandomState(42) #设置随机数种子
X = rnd.uniform(-3, 3, size=100) #random.uniform,从输入的任意两个整数中取出size个随机数
#生成y的思路:先使用NumPy中的函数生成一个sin函数图像,然后再人为添加噪音
y = np.sin(X) + rnd.normal(size=len(X)) / 3 #random.normal,生成size个服从正态分布的随机数
#使用散点图观察建立的数据集是什么样子
plt.scatter(X, y,marker='o',c='k',s=20)
plt.show()
#为后续建模做准备:sklearn只接受二维以上数组作为特征矩阵的输入
X.shape
X = X.reshape(-1, 1)
				
			

3. 使用原始数据进行建模

				
					#使用原始数据进行建模
LinearR = LinearRegression().fit(X, y)
TreeR = DecisionTreeRegressor(random_state=0).fit(X, y)
#放置画布
fig, ax1 = plt.subplots(1)
#创建测试数据:一系列分布在横坐标上的点
line = np.linspace(-3, 3, 1000, endpoint=False).reshape(-1, 1)
#将测试数据带入predict接口,获得模型的拟合效果并进行绘制
ax1.plot(line, LinearR.predict(line), linewidth=2, color='green',
    label="linear regression")
ax1.plot(line, TreeR.predict(line), linewidth=2, color='red',
    label="decision tree")
#将原数据上的拟合绘制在图像上
ax1.plot(X[:, 0], y, 'o', c='k')
#其他图形选项
ax1.legend(loc="best")
ax1.set_ylabel("Regression output")
ax1.set_xlabel("Input feature")
ax1.set_title("Result before discretization")
plt.tight_layout()
plt.show()
#从这个图像来看,可以得出什么结果?
				
			

从图像上可以看出,线性回归无法拟合出这条带噪音的正弦曲线的真实面貌,只能够模拟出大概的趋势,而决策树却通过建立复杂的模型将几乎每个点都拟合出来了。此时此刻,决策树正处于过拟合的状态,对数据的学习过于细致,而线性回归处于拟合不足的状态,这是由于模型本身只能够在线性关系间进行拟合的性质决定的。为了让线性回归在类似的数据上变得更加强大,我们可以使用分箱,也就是离散化连续型变量的方法来处理原始数据,以此来提升线性回归的表现。来看看我们如何实现:
4. 分箱及分箱的相关问题

				
					from sklearn.preprocessing import KBinsDiscretizer
#将数据分箱
enc = KBinsDiscretizer(n_bins=10 #分几类?
           ,encode="onehot") #ordinal
X_binned = enc.fit_transform(X)
#encode模式"onehot":使用做哑变量方式做离散化
#之后返回一个稀疏矩阵(m,n_bins),每一列是一个分好的类别
#对每一个样本而言,它包含的分类(箱子)中它表示为1,其余分类中它表示为0
X.shape
X_binned
#使用pandas打开稀疏矩阵
import pandas as pd
pd.DataFrame(X_binned.toarray()).head()
#我们将使用分箱后的数据来训练模型,在sklearn中,测试集和训练集的结构必须保持一致,否则报错
LinearR_ = LinearRegression().fit(X_binned, y)
LinearR_.predict(line) #line作为测试集
line.shape #测试
X_binned.shape #训练
#因此我们需要创建分箱后的测试集:按照已经建好的分箱模型将line分箱
line_binned = enc.transform(line)
line_binned.shape #分箱后的数据是无法进行绘图的
line_binned
LinearR_.predict(line_binned).shape
				
			

5. 使用分箱数据进行建模和绘图

				
					#准备数据
enc = KBinsDiscretizer(n_bins=10,encode="onehot")
X_binned = enc.fit_transform(X)
line_binned = enc.transform(line)
#将两张图像绘制在一起,布置画布
fig, (ax1, ax2) = plt.subplots(ncols=2
               , sharey=True #让两张图共享y轴上的刻度
               , figsize=(10, 4))
#在图1中布置在原始数据上建模的结果
ax1.plot(line, LinearR.predict(line), linewidth=2, color='green',
    label="linear regression")
ax1.plot(line, TreeR.predict(line), linewidth=2, color='red',
    label="decision tree")
ax1.plot(X[:, 0], y, 'o', c='k')
ax1.legend(loc="best")
ax1.set_ylabel("Regression output")
ax1.set_xlabel("Input feature")
ax1.set_title("Result before discretization")
#使用分箱数据进行建模
LinearR_ = LinearRegression().fit(X_binned, y)
TreeR_ = DecisionTreeRegressor(random_state=0).fit(X_binned, y)
#进行预测,在图2中布置在分箱数据上进行预测的结果
ax2.plot(line #横坐标
    , LinearR_.predict(line_binned) #分箱后的特征矩阵的结果
    , linewidth=2
    , color='green'
    , linestyle='-'
    , label='linear regression')
ax2.plot(line, TreeR_.predict(line_binned), linewidth=2, color='red',
    linestyle=':', label='decision tree')
#绘制和箱宽一致的竖线
ax2.vlines(enc.bin_edges_[0] #x轴
     , *plt.gca().get_ylim() #y轴的上限和下限
     , linewidth=1
     , alpha=.2)
#将原始数据分布放置在图像上
ax2.plot(X[:, 0], y, 'o', c='k')
#其他绘图设定
ax2.legend(loc="best")
ax2.set_xlabel("Input feature")
ax2.set_title("Result after discretization")
plt.tight_layout()
plt.show()
				
			

从图像上可以看出,离散化后线性回归和决策树上的预测结果完全相同了——线性回归比较成功地拟合了数据的分布,而决策树的过拟合效应也减轻了。由于特征矩阵被分箱,因此特征矩阵在每个区域内获得的值是恒定的,因此所有模型对同一个箱中所有的样本都会获得相同的预测值。与分箱前的结果相比,线性回归明显变得更加灵活,而决策树的过拟合问题也得到了改善。但注意,一般来说我们是不使用分箱来改善决策树的过拟合问题的,因为树模型带有丰富而有效的剪枝功能来防止过拟合。在这个例子中,我们设置的分箱箱数为10,不难想到这个箱数的设定肯定会影响模型最后的预测结果,我们来看看不同的箱数会如何影响回归的结果:


6. 箱子数如何影响模型的结果

				
					enc = KBinsDiscretizer(n_bins=5,encode="onehot")
X_binned = enc.fit_transform(X)
line_binned = enc.transform(line)
fig, ax2 = plt.subplots(1,figsize=(5,4))
LinearR_ = LinearRegression().fit(X_binned, y)
print(LinearR_.score(line_binned,np.sin(line)))
TreeR_ = DecisionTreeRegressor(random_state=0).fit(X_binned, y)
ax2.plot(line #横坐标
    , LinearR_.predict(line_binned) #分箱后的特征矩阵的结果
    , linewidth=2
    , color='green'
    , linestyle='-'
    , label='linear regression')
ax2.plot(line, TreeR_.predict(line_binned), linewidth=2, color='red',
    linestyle=':', label='decision tree')
ax2.vlines(enc.bin_edges_[0], *plt.gca().get_ylim(), linewidth=1, alpha=.2)
ax2.plot(X[:, 0], y, 'o', c='k')
ax2.legend(loc="best")
ax2.set_xlabel("Input feature")
ax2.set_title("Result after discretization")
plt.tight_layout()
plt.show()
				
			

7. 如何选取最优的箱数

				
					#怎样选取最优的箱子?
from sklearn.model_selection import cross_val_score as CVS
import numpy as np
pred,score,var = [], [], []
binsrange = [2,5,10,15,20,30]
for i in binsrange:
  #实例化分箱类
  enc = KBinsDiscretizer(n_bins=i,encode="onehot")
  #转换数据
  X_binned = enc.fit_transform(X)
  line_binned = enc.transform(line)
  #建立模型
  LinearR_ = LinearRegression()
  #全数据集上的交叉验证
  cvresult = CVS(LinearR_,X_binned,y,cv=5)
  score.append(cvresult.mean())
  var.append(cvresult.var())
  #测试数据集上的打分结果
  pred.append(LinearR_.fit(X_binned,y).score(line_binned,np.sin(line)))
#绘制图像
plt.figure(figsize=(6,5))
plt.plot(binsrange,pred,c="orange",label="test")
plt.plot(binsrange,score,c="k",label="full data")
plt.plot(binsrange,score+np.array(var)*0.5,c="red",linestyle="--",label = "var")
plt.plot(binsrange,score-np.array(var)*0.5,c="red",linestyle="--")
plt.legend()
plt.show()
				
			

在工业中,大量离散化变量与线性模型连用的实例很多,在深度学习出现之前,这种模式甚至一度统治一些工业中的机器学习应用场景,可见效果优秀,应用广泛。对于现在的很多工业场景而言,大量离散化特征的情况可能已经不是那么多了,不过大家依然需要对“分箱能够解决线性模型无法处理非线性数据的问题”有所了解。


5.3 多项式回归PolynomialFeatures
5.3.1 多项式对数据做了什么
除了分箱之外,另一种更普遍的用于解决”线性回归只能处理线性数据“问题的手段,就是使用多项式回归对线性回归进行改进。这样的手法是机器学习研究者们从支持向量机中获得的:支持向量机通过升维可以将非线性可分数据转化为线性可分,然后使用核函数在低维空间中进行计算,这是一种“高维呈现,低维解释”的思维。那我们为什么不能让线性回归使用类似于升维的转换,将数据由非线性转换为线性,从而为线性回归赋予处理非线性数据的能力呢?当然可以。

接下来,我们就来看看线性模型中的升维工具:多项式变化。这是一种通过增加自变量上的次数,而将数据映射到高维空间的方法,只要我们设定一个自变量上的次数(大于1),就可以相应地获得数据投影在高次方的空间中的结果。这种方法可以非常容易地通过sklearn中的类PolynomialFeatures来实现。我们先来简单看看这个类是如何使用的。
class sklearn.preprocessing.PolynomialFeatures (degree=2, interaction_only=False, include_bias=True)

				
					from sklearn.preprocessing import PolynomialFeatures
import numpy as np
#如果原始数据是一维的
X = np.arange(1,4).reshape(-1,1)
X
#二次多项式,参数degree控制多项式的次方
poly = PolynomialFeatures(degree=2)
#接口transform直接调用
X_ = poly.fit_transform(X)
X_
X_.shape
#三次多项式
PolynomialFeatures(degree=3).fit_transform(X)
				
			
				
					#三次多项式,不带与截距项相乘的x0
PolynomialFeatures(degree=3,include_bias=False).fit_transform(X)
#为什么我们会希望不生成与截距相乘的x0呢?
#对于多项式回归来说,我们已经为线性回归准备好了x0,但是线性回归并不知道
xxx = PolynomialFeatures(degree=3).fit_transform(X)
xxx.shape
rnd = np.random.RandomState(42) #设置随机数种子
y = rnd.randn(3)
y
#生成了多少个系数?
LinearRegression().fit(xxx,y).coef_
#查看截距
LinearRegression().fit(xxx,y).intercept_
#发现问题了吗?线性回归并没有把多项式生成的x0当作是截距项
#所以我们可以选择:关闭多项式回归中的include_bias
#也可以选择:关闭线性回归中的fit_intercept
#生成了多少个系数?
LinearRegression(fit_intercpet=False).fit(xxx,y).coef_
#查看截距
LinearRegression(fit_intercpet=False).fit(xxx,y).intercept_
				
			

不过,这只是一维状况的表达,大多数时候我们的原始特征矩阵不可能会是一维的,至少也是二维以上,很多时候还可能存在上千个特征或者维度。现在我们来看看原始特征矩阵是二维的状况:

				
					X = np.arange(6).reshape(3, 2)
X
#尝试二次多项式
PolynomialFeatures(degree=2).fit_transform(X)
				
			
				
					#尝试三次多项式
PolynomialFeatures(degree=3).fit_transform(X)
				
			
				
					PolynomialFeatures(degree=2).fit_transform(X)
PolynomialFeatures(degree=2,interaction_only=True).fit_transform(X)
#对比之下,当interaction_only为True的时候,只生成交互项
				
			

从之前的许多次尝试中我们可以看出,随着多项式的次数逐渐变高,特征矩阵会被转化得越来越复杂。不仅是次数,当特征矩阵中的维度数(特征数)增加的时候,多项式同样会变得更加复杂:

				
					#更高维度的原始特征矩阵
X = np.arange(9).reshape(3, 3)
X
PolynomialFeatures(degree=2).fit_transform(X)
PolynomialFeatures(degree=3).fit_transform(X)
X_ = PolynomialFeatures(degree=20).fit_transform(X)
X_.shape
				
			

如此,多项式变化对于数据会有怎样的影响就一目了然了:随着原特征矩阵的维度上升,随着我们规定的最高次数的上升,数据会变得越来越复杂,维度越来越多,并且这种维度的增加并不能用太简单的数学公式表达出来。因此,多项式回归没有固定的模型表达式,多项式回归的模型最终长什么样子是由数据和最高次数决定的,因此我们无法断言说某个数学表达式”就是多项式回归的数学表达”,因此要求解多项式回归不是一件容易的事儿,感兴趣的大家可以自己去尝试看看用最小二乘法求解多项式回归。接下来,我们就来看看多项式回归的根本作用:处理非线性问题。


5.3.2 多项式回归处理非线性问题
之前我们说过,是希望通过这种将数据投影到高维的方式来帮助我们解决非线性问题。那我们现在就来看一看多项式转化对模型造成了什么样的影响:

				
					from sklearn.preprocessing import PolynomialFeatures as PF
from sklearn.linear_model import LinearRegression
import numpy as np
rnd = np.random.RandomState(42) #设置随机数种子
X = rnd.uniform(-3, 3, size=100)
y = np.sin(X) + rnd.normal(size=len(X)) / 3
#将X升维,准备好放入sklearn中
X = X.reshape(-1,1)
#创建测试数据,均匀分布在训练集X的取值范围内的一千个点
line = np.linspace(-3, 3, 1000, endpoint=False).reshape(-1, 1)
#原始特征矩阵的拟合结果
LinearR = LinearRegression().fit(X, y)
#对训练数据的拟合
LinearR.score(X,y)
#对测试数据的拟合
LinearR.score(line,np.sin(line))
#多项式拟合,设定高次项
d=5
#进行高此项转换
poly = PF(degree=d)
X_ = poly.fit_transform(X)
line_ = PF(degree=d).fit_transform(line)
#训练数据的拟合
LinearR_ = LinearRegression().fit(X_, y)
LinearR_.score(X_,y)
#测试数据的拟合
LinearR_.score(line_,np.sin(line))
				
			

如果我们将这个过程可视化:

				
					import matplotlib.pyplot as plt
d=5
#和上面展示一致的建模流程
LinearR = LinearRegression().fit(X, y)
X_ = PF(degree=d).fit_transform(X)
LinearR_ = LinearRegression().fit(X_, y)
line = np.linspace(-3, 3, 1000, endpoint=False).reshape(-1, 1)
line_ = PF(degree=d).fit_transform(line)
#放置画布
fig, ax1 = plt.subplots(1)
#将测试数据带入predict接口,获得模型的拟合效果并进行绘制
ax1.plot(line, LinearR.predict(line), linewidth=2, color='green'
    ,label="linear regression")
ax1.plot(line, LinearR_.predict(line_), linewidth=2, color='red'
    ,label="Polynomial regression")
#将原数据上的拟合绘制在图像上
ax1.plot(X[:, 0], y, 'o', c='k')
#其他图形选项
ax1.legend(loc="best")
ax1.set_ylabel("Regression output")
ax1.set_xlabel("Input feature")
ax1.set_title("Linear Regression ordinary vs poly")
plt.tight_layout()
plt.show()
#来一起鼓掌,感叹多项式回归的神奇
#随后可以试试看较低和较高的次方会发生什么变化
#d=2
#d=20
				
			

从这里大家可以看出,多项式回归能够较好地拟合非线性数据,还不容易发生过拟合,可以说是保留了线性回归作为线性模型所带的“不容易过拟合”和“计算快速”的性质,同时又实现了优秀地拟合非线性数据。到了这里,相信大家对于多项式回归的效果已经不再怀疑了。多项式回归非常迷人也非常神奇,因此一直以来都有各种各样围绕着多项式回归进行的讨论。在这里,为大家梳理几个常见问题和讨论,供大家参考。


5.3.3 多项式回归的可解释性
线性回归是一个具有高解释性的模型,它能够对每个特征拟合出参数 以帮助我们理解每个特征对于标签的作用。当我们进行了多项式转换后,尽管我们还是形成形如线性回归的方程,但随着数据维度和多项式次数的上升,方程也变得异常复杂,我们可能无法一眼看出增维后的特征是由之前的什么特征组成的(之前我们都是肉眼看肉眼判断)。不过,多项式回归的可解释性依然是存在的,我们可以使用接口get_feature_names来调用生成的新特征矩阵的各个特征上的名称,以便帮助我们解释模型。来看下面的例子:

				
					import numpy as np
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
X = np.arange(9).reshape(3, 3)
X
poly = PolynomialFeatures(degree=5).fit(X)
#重要接口get_feature_names
poly.get_feature_names()
				
			

使用加利佛尼亚房价数据集给大家作为例子,当我们有标签名称的时候,可以直接在接口get_feature_names()中输入标签名称来查看新特征究竟是由原特征矩阵中的什么特征组成的:

				
					from sklearn.datasets import fetch_california_housing as fch
import pandas as pd
housevalue = fch()
X = pd.DataFrame(housevalue.data)
y = housevalue.target
housevalue.feature_names
X.columns = ["住户收入中位数","房屋使用年代中位数","平均房间数目"
     ,"平均卧室数目","街区人口","平均入住率","街区的纬度","街区的经度"]
poly = PolynomialFeatures(degree=2).fit(X,y)
poly.get_feature_names(X.columns)
X_ = poly.transform(X)
#在这之后,我们依然可以直接建立模型,然后使用线性回归的coef_属性来查看什么特征对标签的影响最大
reg = LinearRegression().fit(X_,y)
coef = reg.coef_
[*zip(poly.get_feature_names(X.columns),reg.coef_)]
#放到dataframe中进行排序
coeff = pd.DataFrame([poly.get_feature_names(X.columns),reg.coef_.tolist()]).T
coeff.columns = ["feature","coef"]
coeff.sort_values(by="coef")
				
			

可以发现,不仅数据的可解释性还存在,我们还可以通过这样的手段做特征工程——特征创造。多项式帮助我们进行了一系列特征之间相乘的组合,若能够找出组合起来后对标签贡献巨大的特征,那我们就是创造了新的有效特征,对于任何学科而言发现新特征都是非常有价值的。在加利佛尼亚房屋价值数据集上来再次确认多项式回归提升模型表现的能力:

				
					#顺便可以查看一下多项式变化之后,模型的拟合效果如何了
poly = PolynomialFeatures(degree=4).fit(X,y)
X_ = poly.transform(X)
reg = LinearRegression().fit(X,y)
reg.score(X,y)
from time import time
time0 = time()
reg_ = LinearRegression().fit(X_,y)
print("R2:{}".format(reg_.score(X_,y)))
print("time:{}".format(time()-time0))
#假设使用其他模型?
from sklearn.ensemble import RandomForestRegressor as RFR
time0 = time()
print("R2:{}".format(RFR(n_estimators=100).fit(X,y).score(X,y)))
print("time:{}".format(time()-time0))
				
			

到这里,多项式回归就全部讲解完毕了。多项式回归主要是通过对自变量上的次方进行调整,来为线性回归赋予更多的学习能力,它的核心表现在提升模型在现有数据集上的表现。


6 结语
本章之中,大家学习了多元线性回归,岭回归,Lasso和多项式回归总计四个算法,他们都是围绕着原始的线性回归进行的拓展和改进。其中岭回归和Lasso是为了解决多元线性回归中使用最小二乘法的各种限制,主要用途是消除多重共线性带来的影响并且做特征选择,而多项式回归解决了线性回归无法拟合非线性数据的明显缺点,核心作用是提升模型的表现。除此之外,本章还定义了多重共线性和各种线性相关的概念,并为大家补充了一些线性代数知识。回归算法属于原理简单,但操作困难的机器学习算法,在实践和理论上都还有很长的路可以走,希望大家继续探索,让线性回归大家族中的算法真正称为大家的武器。