机器学习入门——多项式回归

引言

我们之前介绍的线性回归法有个很大的局限性，它假设数据存在线性关系。但在实际中，具有线性关系的数据是不常见的。更常见的是非线性关系的数据。

本文介绍的多项式回归就是线性回归的一种改进，可以处理非线性关系的数据。

---

• 上一篇：机器学习入门——详解主成分分析
• 下一篇：机器学习入门——图解逻辑回归

多项式回归

在线性回归中我们通过直线$y=ax+b$来拟合上图中的数据，这里假设样本只有一个特征$x$。

但是有些数据，比如上图所示的这种数据，显然用这种二次曲线可以更好的拟合。如果所有的样本也只有一个特征，那么二次曲线方程可以写成：$y=a{x}^2+bx+c$。

我们可以从另外一个角度来理解这个式子，如果把$x^2$理解成一个特征，$x$理解成另外一个特征。这样这式子依然是线性回归的式子，但从$x$的角度来看是一个非线性的方程。

这样的方式就叫多项式回归，相当于我们为样本多添了特征，这些特征是原来样本的多项式项。

下面我们编程实现一下多项式回归。

import numpy as np
import matplotlib.pyplot as plt

x = np.random.uniform(-3,3,size=100)#生成-3到3的100个随机数
X = x.reshape(-1,1) #变成二维数组

y = 0.5 * x**2 + x  + 2 + np.random.normal(0,1,size=100) #定义y=0.5x^2+x+2的函数，加上噪音

plt.scatter(x,y)
plt.show()

首先生成二次项的数据并画图，从这些数据分布来看，显然是非线性的关系。但是我们先用之前介绍的线性回归来拟合一下。

from sklearn.linear_model import LinearRegression
lin_reg = LinearRegression()
lin_reg.fit(X,y)

y_predict = lin_reg.predict(X)
plt.scatter(x,y)
plt.plot(x,y_predict,color='r')
plt.show()

显然这用一根直线来拟合这种有弧度的曲线分布样本的效果是不太好的。

下面我们来看下如何用多项式回归来拟合。

就像上面所说的，我们对$X$中的每个数据进行平方，作为一个新的特征。

 X_new = np.hstack ([ X , X ** 2 ] )
 X_new 

可以看到在这个新的数据集中，每个样本有两个特征。现在用这个新的数据集来进行线性回归训练。

lin_reg2 = LinearRegression()
lin_reg2.fit(X_new,y)
y_predict2 = lin_reg2.predict(X_new)

plt.scatter(x,y)
plt.plot(x,y_predict2,color='r')
plt.show()

哇，结果好像坏掉了。其实并不是的，因为x是乱序的。我们想要生成一条平滑曲线，需要对x进行排序。

lin_reg2 = LinearRegression()
lin_reg2.fit(X_new,y)
y_predict2 = lin_reg2.predict(X_new)

plt.scatter(x,y)
plt.plot(np.sort(x),y_predict2[np.argsort(x)],color='r')
plt.show()

最终绘制的就是如上图所示的曲线。
这就是我们将原来的数据集X添加了一个特征(对所有样本的特征进行平方)，当我们添加了这个特征后，从x的维度来看，就形成了一条曲线。

显然这条曲线对我们的数据集的拟合程度是更好的。

回顾下我们的样本X_new = np.hstack([X,X**2])，第一列是$X$，第二列是$X^2$，也就是得到的$X$前的系数大概是$1$，$X^2$前的系数大概是$0.5$，这和我们实际方程式是一样的。

因为添加了噪音，所以无法完全一样。

截距是1.87，和方程里面的2也很接近。

这就是多项式回归，从这里可以看出，多项式回归完全是使用了线性回归的思路，关键在于为原来的样本通过多项式组合添加了新的特征。

这样我们就可以解决非线性问题。

sklearn中的多项式回归

import numpy as np
import matplotlib.pyplot as plt

x = np.random.uniform(-3,3,size=100)#生成-3到3的100个随机数
X = x.reshape(-1,1) #变成二维数组

y = 0.5 * x**2 + x  + 2 + np.random.normal(0,1,size=100) #定义y=0.5x^2+x+2的函数，加上噪音

采用和上小节一样的数据。我们在上小节使用线性回归主要改造了$X$，添加了特征。

sklearn中便是采用这种思路，也是为数据添加多项式特征。

from sklearn.preprocessing import PolynomialFeatures
poly = PolynomialFeatures(degree=2) #增加X^2特征
poly.fit(X)
X2 = poly.transform(X)# 对X进行转换
X2.shape

但是为什么输出了3个特征，我们来看下前10行样子。

多了一列1，其实1可以看成是$x^0=1$，而第二列就是原来的样本特征，第三列是$X^2$。

来验证一下第二列和第三列。

我们获取到了多项式特征的数据集后，也是直接调用线性回归的类即可：

from sklearn.linear_model import LinearRegression
# 还是调用线性回归的方法
lin_reg2 = LinearRegression()
lin_reg2.fit(X2,y)
y_predict2 = lin_reg2.predict(X2)

plt.scatter(x,y)
plt.plot(np.sort(x),y_predict2[np.argsort(x)],color='r')
plt.show()

这样就完成了多项式回归，这里也可以看一下系数。

上面我们生成的$X$只有一个特征，如果有两个特征呢

X = np.arange(1,11).reshape(-1,2)
poly = PolynomialFeatures(degree=2)
poly.fit(X)
X2 = poly.transform(X)
X2.shape #(5,6)

为什么现在生成了6个系数呢？

上图是我们生成的X。

对比分析的话，X2的第一列是全1，第二列是X的第一列，X2的第三列是X的第二列，X2的第4列是X第一列的平方，第6列是X第二列的平方。那第5列呢，其实是X两列各元素乘积的结果。

poly = PolynomialFeatures(degree=3)
poly.fit(X)
X3 = poly.transform(X)
X3.shape

如果传入的degree=3三的话，就是最多生成3次幂。从上图看到，生成了10列特征。
是怎么生成的呢，用下图来解释吧。

Pipeline

在这小节介绍下sklearn中常用的一个工具——Pipeline(管道)，它可以将很多模型算法串起来，形成一个类似管道的东西，你只要喂输入数据，经过管道里面的运算(特征提取、归一化、分类等)最终得到输出。

在多项式回归中，如果特征的次数相差太多的话，数据规模很容易相差巨大，因此需要标准化。标准化这个过程就可以添加到Pipeline中。除了标准化外，在上小节的代码中，我们每次将数据传入线性回归类之前，都要实例化一个PolynomialFeatures对象对数据进行多项式处理，这一过程也可以封装到管道中。

下面构造一个Pipeline

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

x = np.random.uniform(-3,3,size=100)
X = x.reshape(-1,1)
y = 0.5 * x**2 + x + 2 +np.random.normal(0,1,100)

# 接受一个元组列表(名称,类)
poly_reg = Pipeline([
    #多项式处理
    ('poly',PolynomialFeatures(degree=2)), 
    #标准化
    ('std_scaler',StandardScaler()), 
    #传入线性回归
    ('lin_reg',LinearRegression()) 
])

训练和预测很简单，只要执行

poly_reg.fit(X,y)
y_predict =poly_reg.predict(X)

即可，可以看到，相当于我们实现了一个多项式回归的类，可以避免重复代码。

plt.scatter(x,y)
plt.plot(np.sort(x),y_predict[np.argsort(x)],color='r')
plt.show()

这里通过画图来验证了一下，可以看到得到的结果是正确的的。

多项式回归中degree越高在训练数据上可以拟合的越好，那是不是每次设置一个很高的degree就万事大吉了呢，看完下一节你就知道答案了。

过拟合与欠拟合

我们之前生成的数据是这样的，然后我们先用线性回归来拟合，得到的结果是下面这样的：

我们说这条直线并不能很好的拟合这些数据点，然后我们就学习到了多项式回归。我们用2次曲线来拟合：

我们说这条曲线拟合的好一点。但是这只是看起来的结果，有没有一个指标能衡量这种拟合程度呢。

我们可以用均方误差来衡量，下面我们来看下代码：

from sklearn.metrics import mean_squared_error

y_predict = lin_reg.predict(X)
mean_squared_error(y,y_predict) #使用线性回归的均方误差

再看下用多项式拟合的均方误差：

y_predict =poly_reg.predict(X)
mean_squared_error(y,y_predict) #使多项式(2次幂)回归的均方误差

可以看到误差显然小于线性回归的误差。

这里用的是上小节中degree=2的多项式，如果我们将这个数值调高结果会如何呢？

def PolynomialRegression(degree):
    # 接受一个元组列表(名称,类)
    return Pipeline([
    #多项式处理
    ('poly',PolynomialFeatures(degree=degree)), #设置为我们传入的参数
    #标准化
    ('std_scaler',StandardScaler()), 
    #传入线性回归
    ('lin_reg',LinearRegression()) 
])

为例避免重复代码，将生成多项式回归管道的代码也放到一个函数中，接收degree参数。然后我们尝试下将degree设成10，结果会如何

poly10_reg = PolynomialRegression(10)
poly10_reg.fit(X,y)

y10_predict =poly10_reg.predict(X)
mean_squared_error(y,y10_predict)

可以看到，结果确实好了一点，我们此时来看下它的曲线。

plt.scatter(x,y)
plt.plot(np.sort(x),y10_predict[np.argsort(x)],color='r')
plt.show()

可以看到，这条线有波浪起伏。

上面看到传入了10结果会好，那直接传入100呢

poly100_reg = PolynomialRegression(100)

poly100_reg.fit(X,y)
y100_predict =poly100_reg.predict(X)
mean_squared_error(y,y100_predict)

结果好了一倍，那图像是怎么的

可以看到图像的上下起伏更加明显，好像特地去拟合了一些点，造成了这种图像。

事实上，这还不是真实的曲线，因为这只是存在的数据的预测值的连接结果。

下面尝试画出真实曲线。

X_plot = np.linspace(-3,3,100).reshape(100,1) #从-3到3之间均匀取值
y_plot = poly100_reg.predict(X_plot)
plt.scatter(x,y)
plt.plot(X_plot[:,0],y_plot,color='r')
plt.axis([-3,3,-1,10])
plt.show()

这就是将degree=100之后多项式回归拟合的结果。为什么degree越高，拟合的越好呢，因为我们有这么多样本点，我们总能找到一根曲线能尽可能将所有的样本点都拟合起来。

看到这个曲线第一个感觉应该是被吓到了吧，毕竟为了拟合边缘的样本点，而形成的垂线一样的曲线也太夸张了。

这种曲线真的能反应样本的走势吗

随便拿一根近似垂直的线来说(其实这么多多根垂线都是连接在一起的，不过这个坐标轴上看不到而已)，按照这根垂直的线的说法，应该也会有很多样本存在这根垂线之上，我们就标记出一些黑点，问题是实际的数据真的会像黑点这样分布吗，我们知道数据是从二次方程中生成的，显然这些黑点是不存在的。

这种情况就是过拟合，也就是针对训练数据拟合的过好，而用在测试数据上反而效果很差。

而我们之前直接用线性回归拟合的直线，情况和上面的正好相反，它甚至不能很好的反应训练数据的趋势。此时我们说它是欠拟合。就是拟合的不够好，我们需要重新训练模型。

为什么要训练数据集与测试数据集

在上面我们提到了训练数据和测试数据，那本小节就来看下为什么要训练数据集与测试数据集。

上面我们说这条曲线是过拟合的。虽然对已知样本的预测能力很高，但是如果碰到新的样本，

假设新样本的$x$是2.8左右，按照预测，得到的位置如上图紫色圆圈所示。显然这个紫色的样本点和我们训练用的样本点不在一个趋势上。

直观上感觉这个预测结果是错误的。
也就是在过拟合的情况下，虽然这条曲线能拟合原来的样本点很好，但是无法预测新的样本点。
这种情况下，我们称我们得到的这个模型的泛化能力(预测新样本的能力)是非常低的。

而我们训练模型的目的就是预测新样本，那此时要怎么评估呢。

很简单，使用训练集和测试集分离的数据集。

我们使用训练集来训练模型，用模型没见过的测试数据来评估模型。如果模型在测试数据上也表现很好，我们就说这个模型泛化能力是很强的。

下面我们来看下我们之前的模型对于测试数据根据degree不同，误差是如何变化的。

from sklearn.model_selection import train_test_split
X_train,X_test,y_train,y_test = train_test_split(X,y,random_state=777)

先把数据集分为训练集和测试集。

lin_reg = LinearRegression()
# 用训练数据训练
lin_reg.fit(X_train,y_train)
# 用测试数据测试
y_predict = lin_reg.predict(X_test)
mean_squared_error(y_test,y_predict)

可以看到，使用线性回归得到的误差是2.10。
那使用多项式回归的结果是怎样呢

poly2_reg = PolynomialRegression(degree=2)
poly2_reg.fit(X_train,y_train)
y2_predict = poly2_reg.predict(X_test)
mean_squared_error(y_test,y2_predict)

先是2项式回归，结果是0.80，说明在这个数据集中，使用2阶多项式泛化能力比线性模型的要好。

那我们使用10阶多项式呢。

poly10_reg = PolynomialRegression(degree=10)
poly10_reg.fit(X_train,y_train)
y10_predict = poly10_reg.predict(X_test)
mean_squared_error(y_test,y10_predict)

可以看到，在测试集上的误差比2阶还要高，说明泛化能力变差了。

这就是我们使用测试集来测试模型的泛化能力的一个方式。
可以看到degree越大，模型越复杂，下面我们看一下模型复杂度与在训练集和测试集上准确率的关系。

大概示意图如上。随着模型复杂度逐渐增加，训练集上的准确率越来越好，而测试集上的准确率先是上增，然后是下降。

上面我们说了过拟合，还有一种叫欠拟合。所谓欠拟合，就是拟合的不够好，甚至在训练数据上也表现不好。

比如对于这种数据，我们使用一根直线去拟合，显然是无法很好的描述这些数据的趋势。

学习曲线

我们上小节看到了这样的模型复杂度的曲线，我们学习的目标其实是要找到测试集上最好的点所对应的模型参数。

但是上面这个图像只是理论上的示例，对于不同的数据、不同的模型图形都是不同的。

对于过拟合与欠拟合，还有另外一种曲线可以可视化的方式来观察到，就是我们本小节的重点——学习曲线。

学习曲线描述的是随着训练样本的逐渐增多，算法训练出的模型的表现能力。

下面来看一个实例。

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(777)

x = np.random.uniform(-3.0,3.0,size=100)
X = x.reshape(-1,1)
y = 0.5 * x**2 + x + 2 + np.random.normal(0,1,size=100)

plt.scatter(x,y)
plt.show()

首先生成数据如上。

现在来看下学习曲线是怎么绘制的。

from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error

X_train,X_test,y_train,y_test = train_test_split(X,y,random_state=777)

train_score = []
test_score = []

# 每次多一点训练数据
for i in range(1,X_train.shape[0]+1):
    lin_reg = LinearRegression()
    # 每次取前i个样本
    lin_reg.fit(X_train[:i],y_train[:i])
    # 在训练集上的预测结果
    y_train_predict = lin_reg.predict(X_train[:i])
    train_score.append(mean_squared_error(y_train[:i],y_train_predict))
    # 在测试集上的预测结果 ，这里要在整个测试集上进行判断
    y_test_predict = lin_reg.predict(X_test)
    test_score.append(mean_squared_error(y_test,y_test_predict))


# 绘制学习曲线
# 横坐标表示样本数逐渐增加
# 纵坐标表示对应的误差
plt.plot([i for i in range(1,X_train.shape[0]+1)],np.sqrt(train_score),label='train')
plt.plot([i for i in range(1,X_train.shape[0]+1)],np.sqrt(test_score),label='test')
plt.legend()
plt.show()

  

上面就是我们使用线性回归得到的学习曲线，可以看到，在训练集上的误差是逐渐升高后趋于平缓的，而在测试集上的误差先是升高，然后又降低，最终趋于平缓。

为了方便绘制不同算法的学习曲线，我们抽出上面的代码形成一个函数：

def plot_learning_curve(algo,X_train,X_test,y_train,y_test,title=''):
    train_score = []
    test_score = []

    # 每次多一点训练数据
    for i in range(1,len(X_train)+1):
        # 每次取前i个样本
        algo.fit(X_train[:i],y_train[:i])
        # 在训练集上的预测结果
        y_train_predict = algo.predict(X_train[:i])
        train_score.append(mean_squared_error(y_train[:i],y_train_predict))
        # 在测试集上的预测结果 ，这里要在整个测试集上进行判断
        y_test_predict = algo.predict(X_test)
        test_score.append(mean_squared_error(y_test,y_test_predict))


    # 绘制学习曲线
    # 横坐标表示样本数逐渐增加
    # 纵坐标表示对应的误差
    plt.plot([i for i in range(1,len(X_train)+1)],np.sqrt(train_score),label='train')
    plt.plot([i for i in range(1,len(X_train)+1)],np.sqrt(test_score),label='test')
    plt.legend()
    plt.title(title)
    plt.axis([0,len(X_train) + 1,0,4])
    plt.show()

然后我们用刚才的线性回归算法测试一下：

下面我们看看使用多项式回归，学习曲线是怎样的。

from sklearn.preprocessing import PolynomialFeatures
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline

def PolynomialRegression(degree):
    # 接受一个元组列表(名称,类)
    return Pipeline([
    #多项式处理
    ('poly',PolynomialFeatures(degree=degree)), #设置为我们传入的参数
    #标准化
    ('std_scaler',StandardScaler()), 
    #传入线性回归
    ('lin_reg',LinearRegression()) 
])

for i in range(2,21): #传入degree从2到20
	plot_learning_curve(PolynomialRegression(i),X_train,X_test,y_train,y_test,title='degree =' + str(i))

这就是2阶多项式回归得到的学习曲线，我们可以看到，从趋势上和线性回归是一致的。

不过有一个重大区别是这两根曲线稳定的点（误差）是在1.0左右。而线性回归稳定的点是在1.5左右。

说明2阶多项式回归拟合的结果比线性回归要好。

当degree=20的时候，可以看到虽然在训练集上拟合的很好，稳点的点比2阶还要低，但是在测试集上的跳跃性非常大，从最后一个横坐标对应的误差值来看，反而有一个向上的趋势。说明是过拟合了。

比较欠拟合（线性回归）和最佳拟合（2阶多项式回归）的学习曲线可以看到，欠拟合的曲线在训练集和测试集上的表现都要弱于最佳拟合的情况。而且这两根曲线相差较大。

而对于过拟合的曲线，训练集上的误差是最小的，但是测试集上的误差非常大。

验证数据集与交叉验证

验证数据集

上面我们把数据集分为训练数据和测试数据也还存在一个问题，就是我们的模型可能会对测试数据集过拟合。

因为我们是根据测试数据来评估模型的好坏的，一旦模型不好，我们就会调整参数，直到我们的模型在测试集上表现也不错。那么这种方式就很有可能使得我们的模型对测试集过拟合。

解决这个问题的方法也很简单，再多分出一个数据集作为验证数据集。

现在是这样的流程：我们再训练数据上进行学习，然后在验证集上继续验证；如果有问题就调整参数，直到验证集的表现很好。然后我们最终通过测试数据来进行测试。

在这个过程中，测试数据是没有参与模型的创建的，只有训练数据和验证数据参与了模型的训练。相当于我们把测试数据当成真实的线上数据。

现在验证数据集成了调整超参数使用的数据集；
测试数据集作为衡量最终模型性能的数据集。

交叉验证

有一点要注意的是，在只有一份验证数据集的情况下，可能验证数据集中存在极端的数据，使得我们训练的模型有可能过拟合这些极端数据。

为了解决这个问题，就需要用到交叉验证(Cross Validation)。

交叉验证，就是将 训练数据分词k份,这里假设k=3，然后使用其中2份作为训练数据，剩下的那份作为验证数据。
这样每一种训练数据集/验证数据集的搭配就会产生一个模型，每个模型都会在验证集上产生一个性能指标，我们平均这个性能指标就得到了当前模型的性能指标。

现在我们相当于有了k份验证数据集，并且有了一个求平均的过程，可以减少某个验证集中有极端数据的影响。

下面我们看使用交叉验证方式的实例，首先看一下只拆分为训练集和测试集的情况。

import numpy as np
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.neighbors import KNeighborsClassifier

# 手写数字
digits = datasets.load_digits()
X = digits.data
y = digits.target

X_train,X_test,y_train,y_test = train_test_split(X,y,test_size=0.4,random_state=777)

best_score,best_p,best_k = 0,0,0
# 找到knn中最好的k
for k in range(2,11):
    for p in range(1,6):
        knn_clf = KNeighborsClassifier(weights='distance',n_neighbors=k,p=p)
        knn_clf.fit(X_train,y_train)
        score = knn_clf.score(X_test,y_test)
        if score > best_score:
            best_score,best_p,best_k = score,p,k

print('best k = ', best_k)
print('best p = ', best_p)
print('best best_score = ', best_score)

可参考knn详解和统计学习方法中的相关理论。

以knn算法为例，实现手写数字识别，上面打印出了最好的参数。其中p是明可夫斯基距离的参数。

接下来我们看一下交叉验证是如何实现的。

from sklearn.model_selection import cross_val_score

knn_clf = KNeighborsClassifier()
cross_val_score(knn_clf,X_train,y_train,cv=3)

可以看到，它打印出了每个验证集的准确率。

下面我们就使用交叉验证来进行调参。

digits = datasets.load_digits()
X = digits.data
y = digits.target

X_train,X_test,y_train,y_test = train_test_split(X,y,test_size=0.4,random_state=777)

best_score,best_p,best_k = 0,0,0
# 找到knn中最好的k
for k in range(2,11):
    for p in range(1,6):
        knn_clf = KNeighborsClassifier(weights='distance',n_neighbors=k,p=p)
        knn_clf.fit(X_train,y_train)
        scores = cross_val_score(knn_clf,X_train,y_train)
        score = np.mean(scores) #求均值
        if score > best_score:
            best_score,best_p,best_k = score,p,k

print('best k = ', best_k)
print('best p = ', best_p)
print('best best_score = ', best_score)

可以看到，得到的结果和train_test_split是不一样的，一般来说，我们可以更加相信交叉验证的结果。因为在
train_test_split可能会过拟合测试集中的极端数据。

我们交叉验证的过程只是为了找到最好的参数，要计算我们模型的准确率，还是要在train_test_split中的测试集去考量：

best_knn_clf = KNeighborsClassifier(weights='distance',n_neighbors=2,p=2)
best_knn_clf.fit(X_train,y_train)
best_knn_clf.score(X_test,y_test)

这样就得到了用我们最佳参数的knn算法的准确率。

总结上面的过程，我们在训练数据中使用交叉验证得到了最好的参数，然后在没有看过的测试数据上进行最终的测试。

其实这个过程，在网格搜索中已经帮我们进行了。下面我们就来回顾一下网格搜索。

from sklearn.model_selection import GridSearchCV
knn_clf = KNeighborsClassifier()

param_grid = [
    {
        "weights" : ['distance'],
        'n_neighbors' : [i for i in range(2,11)],
        'p' : [i for i in range(6)]
    }
]

# GridSearchCV 后面的CV就是交叉验证的意思
grid_search = GridSearchCV(knn_clf,param_grid,verbose=1)
grid_search.fit(X_train,y_train)

执行网格搜索后，后面打印出的信息意味着用了5折交叉验证，k有9种取值，p有5种取值，一种就是有45种组合。需要搜索5 x 45 =225次。

网格搜索得到的最好参数和我们上面得到的一致。

偏差与方差平衡

我们先来看一下什么是偏差(bias)和方差(variance)。以打靶为例。

我们看最左下角的那个靶子，我们的目标是中心的红点，我们所有的点完全偏离了中心的位置，这种情况就叫高偏差。

我们再看右上角的那个靶子，虽然射出的子弹有些命中了红点(低偏差)，但是太过分散了，这种情况就叫高方差。

总结一下，如果子弹很集中，说明是低方差的；如果子弹很少击中靶心，说明是高偏差的。因此，总共有四种情况。

类比我们机器学习的过程，我们要训练的机器学习模型，都是为了解决一个问题。问题本身我们可以看成是靶子。

我们根据数据来拟合这个模型，我们拟合出来的模型，就是打出去的子弹。
一般来说，我们模型的误差来自三方面：偏差+方差+不可避免的误差。

不可避免的误差就是客观存在的误差，比如我们采集的数据本身就可能是有噪音的。
因此，我们重点要关注的是和我们模型相关的偏差和方差。

我们说一个模型有偏差，很可能的原因是我们对问题本身的假设不正确。
比如对于非线性的数据使用线性模型，欠拟合就是这样的一个例子。

而方差的表现就在于，数据的一点扰动多会极大地影响我们的模型。主要原因是模型过于复杂，比如上面介绍的高阶多项式回归。
换句话说，我们的模型学习到了很多的噪音。过拟合就是一个例子。

对于机器学习算法来说，有些算法天生就是高方差的算法，比如KNN这种非参数学习算法，因为不对数据进行任何假设，只能根据现有训练数据进行预测，所以极大地依赖于训练数据的准确性。

而有一些算法天生是高偏差的，比如线性回归这种参数学习算法，因为这种算法对数据具有极强的假设，比如线性回归就假设我们的训练数据是线性的。

在大多数算法中，我们可以通过调整参数来调整偏差和方差。比如knn中的k，或者我们再线性回归中使用多项式回归。

偏差和方差通常是互相矛盾的，降低偏差会提高方差；而降低方差会提高偏差。

一般情况我们很难找到低偏差同时低方差的算法，所以我们要在方差和偏差之间找到一个平衡点。

针对算法来说，机器学习的主要挑战来自于方差。换句话说，我们很容易让我们的模型非常复杂，从而过拟合。

那么解决高方差(过拟合)的手段有哪些呢

• 降低模型复杂度
• 降维：减少数据维度
• 增加样本数量
• 使用验证集
• 模型正则化

下面我们就来探讨下模型正则化。

模型正则化

模型正则化(Regularization)目的是限制参数的大小。

这个图是我们多项式回归中过拟合的一个例子。从上图可以看出，在这个模型中，参数前面的系数是非常大的。

from sklearn.preprocessing import PolynomialFeatures
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.linear_model import LinearRegression



lin_reg = LinearRegression()

def PolynomialRegression(degree):
    # 接受一个元组列表(名称,类)
    return Pipeline([
    #多项式处理
    ('poly',PolynomialFeatures(degree=degree)), #设置为我们传入的参数
    #标准化
    ('std_scaler',StandardScaler()), 
    #传入线性回归
    ('lin_reg',lin_reg)
])

poly100_reg = PolynomialRegression(100)

poly100_reg.fit(X,y)
y100_predict =poly100_reg.predict(X)

X_plot = np.linspace(-3,3,100).reshape(100,1) #从-3到3之间均匀取值
y_plot = poly100_reg.predict(X_plot)
plt.scatter(x,y)
plt.plot(X_plot[:,0],y_plot,color='r')
plt.axis([-3,3,-1,10])
plt.show()

我们设置degree=100，然后看我们的系数是多少。

可以看到参数的系数大都是非常大的，使得我们的线条非常陡峭。

那我们看如何通过模型正则化来解决这个问题。

我们线性回归问题的损失函数是上面的样子，目标是是均方误差尽可能的小，为了使得我们的参数$\theta$不能过大，我们加入模型正则化。

为了使得加入正则化后的式子尽可能小，除了要使均方误差小之外，还要使${\theta }^2$尽可能小。这样就能防止$\theta$过大。

上式中的$\alpha$来控制正则化的强度，如果$\alpha =0$相当于没有正则化。
这种正则化的方式还叫做$L_2$范数，加上$L_2$范数的线性回归通常有一个名称——岭回归。

下面我们来看一下，使用岭回归的方式来限制模型参数大小后的结果是怎样的。

import numpy as np 
import matplotlib.pyplot as plt

np.random.seed(42)

x = np.random.uniform(-3.0,3.0,size=100)
X = x.reshape(-1,1)
y = 0.5 * x + 3 + np.random.normal(0,1,size=100)

plt.scatter(x,y)
plt.show()

首先生成训练数据，然后用我们上面学到的多项式回归来进行预测。

from sklearn.preprocessing import PolynomialFeatures
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error

np.random.seed(666)
X_train,X_test,y_train,y_test = train_test_split(X,y)


def PolynomialRegression(degree):
    # 接受一个元组列表(名称,类)
    return Pipeline([
	    #多项式处理
	    ('poly',PolynomialFeatures(degree=degree)), #设置为我们传入的参数
	    #标准化
	    ('std_scaler',StandardScaler()), 
	    #传入线性回归
	    ('lin_reg',LinearRegression())
])

poly_reg = PolynomialRegression(20)

poly_reg.fit(X_train,y_train)
y_predict =poly_reg.predict(X_test)
mean_squared_error(y_test,y_predict)

从输出可以看出，均方误差很大。

X_plot = np.linspace(-3,3,100).reshape(100,1) #从-3到3之间均匀取值
y_plot = poly_reg.predict(X_plot)
plt.scatter(x,y)
plt.plot(X_plot[:,0],y_plot,color='r')
plt.axis([-3,3,-1,10])
plt.show()

绘制了我们的模型可以看到，存在很陡峭的地方。
为了后面方便绘制，把绘图代码抽出来了。

def plot_model(model):
    X_plot = np.linspace(-3,3,100).reshape(100,1) #从-3到3之间均匀取值
    y_plot = model.predict(X_plot)
    plt.scatter(x,y)
    plt.plot(X_plot[:,0],y_plot,color='r')
    plt.axis([-3,3,0,6])
    plt.show()

下面使用岭回归来看一下。

from sklearn.linear_model import Ridge

def RidgeRegression(degree,alpha=1):
    # 接受一个元组列表(名称,类)
    return Pipeline([
    #多项式处理
    ('poly',PolynomialFeatures(degree=degree)), #设置为我们传入的参数
    #标准化
    ('std_scaler',StandardScaler()), 
    #传入线性回归
    ('ridge_reg',Ridge(alpha=alpha))
])

ridge1_reg = RidgeRegression(20,alpha=0.0001)

ridge1_reg.fit(X_train,y_train)
y1_predict =ridge1_reg.predict(X_test)
mean_squared_error(y_test,y1_predict)

sklearn中的岭回归就是Ridge，哪怕我们设了一个很小的$\alpha =0.0001$，结果也比之前好了很多。

然后我们绘制一下这个岭回归的图像，可以看到和之前的线条相比，缓和了不少。
我们再尝试增加$\alpha$，看看会有什么样的结果。

将$\alpha$设成1之后，结果是好一点了，但是线条似乎优化的不够明显。我们知道我们的数据其实是通过线性函数生成的，只不过增加了一些噪音。 那么我们尝试更大的$\alpha =100$，看能否获得更好的结果。

从输出可以看到，误差反而上升了一点，说明我们的正则化可能过头了。但是得到的图像是更加平滑了，没有陡峭的地方。

如果我们将$\alpha$设成一个很大的数，那么为了是损失函数尽可能的小，只能让参数都为零的。得到的就是上面这样一根直线。

LASSO回归

LASSO回归实际上是加了$L_1$正则，就是使得参数的绝对值之和尽可能小，当然还有一个超参数$\alpha$。

下面我们编程来看下LASSO回归的样子，数据集还是和上小节一样。

from sklearn.linear_model import Lasso


def RidgeRegression(degree,alpha):
    # 接受一个元组列表(名称,类)
    return Pipeline([
    #多项式处理
    ('poly',PolynomialFeatures(degree=degree)), #设置为我们传入的参数
    #标准化
    ('std_scaler',StandardScaler()), 
    #传入线性回归
    ('lasso_reg',Lasso(alpha=alpha))
])

lasso1_reg = RidgeRegression(20,alpha=0.01)

lasso1_reg.fit(X_train,y_train)
y1_predict =lasso1_reg.predict(X_test)
mean_squared_error(y_test,y1_predict)

初始我们用0.01来进行测试，可以看到得到的误差是比较小的。

看图像的话是十分平缓的，下面我们尝试增加$\alpha$。

当$\alpha =0.1$时，误差小了一点，同时图像更加接近于一条直线。

当$\alpha =1$时，此时也正则化过头了。

我们来看$\alpha =0.1$时的这条曲线，几乎就是一条斜线了。而岭回归大都情况下还是曲线。

这是为什么呢，其实是由$L_1$正则化的特殊性决定的。

LASSO趋向于使得一部分$\theta$值变成0，所以还可以用作特征选择。

那些$\theta$值为0的特征，LASSO认为是没用的；LASSO选择的特征是$\theta$值非零的特征。

比较Ridge和LASSO

比较LASSO和Ridge的式子，我们可以想到一些其他的公式。

比如均方误差(MSE)和绝对值误差(MAE)，以及欧拉距离和曼哈顿距离。

在knn中我们学习了明可夫斯基距离，它的公式如下：

我们把它进行泛化，就可以提炼成下面的形式：

就得到了$L_p$范数，如同上面所说，LASSO其实就是$L_1$范数，而Ridge是$L_2$范数（有一点区别就是都不需要开根）。

下面我们介绍一个新的概念——弹性网(Elastic Net)，它其实就是结合了$L_1$正则和$L_2$正则。
它的公式如下：

引入了一个新的超参数$r$，来控制这两种正则的比例。

参考

1.Python3入门机器学习