机器学习（3）

机器学习（3）

本章主要介绍线性回归模型的原理

一，线性回归

线性模型其实就是对输入的特征值加权求和，再加上一个偏置项的常数

其中，是预测值，n是特征的数量，是第i个特征值，是第j个模型参数

这也可以更为简洁的向量化形式表达

其中，是的转置向量（为行向量，而不再是列向量），X是实例的特征向量。

我们如何来训练模型呢？训练模型就是设置模型参数直到模型最适应训练集的过程。要到达这个目的，我们首先需要知道怎么衡量模型对训练数据的拟合程度是好还是差。回归模型最常见的性能指标是均方根误差（RMSE）,因此在训练线性回归模型时，需要找到最小的RMSE的的值。在实践中将均方误差（MSE）最小化比最小化RMSE更为简单，二者效果相同。线性回归模型的MSE成本函数：

可以直接改写成：

1.1，标准方程

推导：https://blog.csdn.net/joob000/article/details/81295144

为了得到使成本函数最小的值，有一个闭式解方法–也就是一个直接得出结果的数学方程，即标准方程：

其中，是使成本函数最小的值，y是目标值向量。

我们生成一些线性数据来预测这个公式：

import numpy as np
import matplotlib.pyplot as plt
X = 2 * np.random.rand(100,1)#随机生成一些x点
y = 4 + 3 * X + np.random.randn(100,1)#y = 4 + 3 * X + 高斯噪声
X_b = np.c_[np.ones((100, 1)), X]  # 添加 x0 = 1 到每一个实例中
theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)#按照标准方程计算theta值
print(theta_best)
#得到结果3.87和3.17，然而我们期望的是4和3，不过已经很接近了。接下来我们开始预测
X_new = np.array([[0], [2]])#生成x的值为0和2
#np.r_是按列连接两个矩阵，就是把两矩阵上下相加，要求列数相等。
#np.c_是按行连接两个矩阵，就是把两矩阵左右相加，要求行数相等
X_new_b = np.c_[np.ones((2, 1)), X_new]  # 增加 x0 = 1 到每一个实例
y_predict = X_new_b.dot(theta_best)#dot函数返回两个数组的点积

plt.plot(X_new, y_predict, "r-")
plt.plot(X, y, "b.")
plt.axis([0, 2, 0, 15])
plt.show()

结果如图（因数据随机，所以每次运行都略有差别）：

Scikit-Learn的等效代码为：

import numpy as np
X = 2 * np.random.rand(100,1)
y = 4 + 3 * X + np.random.randn(100,1)
X_new = np.array([[0], [2]])
from sklearn.linear_model import LinearRegression
lin_reg = LinearRegression()
lin_reg.fit(X, y)
print(lin_reg.intercept_, lin_reg.coef_)
print(lin_reg.predict(X_new))

结果为：

1.2，计算复杂度

标准方程求逆的矩阵是一个nxn的矩阵（n是特征数量），对这种矩阵求逆的计算复杂度通常为之间。换句话说，如果将特征数量翻倍，那么计算时间将乘以大约倍 到倍之间。

特征数量比较大时标准方程的计算将极其缓慢。

好的一面是，相对于训练集中的实例数量来说，方程是线性的，所以能够有效的处理大量的数据集，只要内存足够。同样，线性回归模型一经训练，预测就非常快速，因为计算复杂度相对于想要预测的实例数量和特征数量来说都是线性的。

二，梯度下降

梯度下降法的基本思想可以类比为一个下山的过程。假设这样一个场景：一个人被困在山上，需要从山上下来(i.e. 找到山的最低点，也就是山谷)。但此时山上的浓雾很大，导致可视度很低。因此，下山的路径就无法确定，他必须利用自己周围的信息去找到下山的路径。这个时候，他就可以利用梯度下降算法来帮助自己下山。具体来说就是，以他当前的所处的位置为基准，寻找这个位置最陡峭的地方，然后朝着山的高度下降的地方走，同理，如果我们的目标是上山，也就是爬到山顶，那么此时应该是朝着最陡峭的方向往上走。然后每走一段距离，都反复采用同一个方法，最后就能成功的抵达山谷。

我们同时可以假设这座山最陡峭的地方是无法通过肉眼立马观察出来的，而是需要一个复杂的工具来测量，同时，这个人此时正好拥有测量出最陡峭方向的能力。所以，此人每走一段距离，都需要一段时间来测量所在位置最陡峭的方向，这是比较耗时的。那么为了在太阳下山之前到达山底，就要尽可能的减少测量方向的次数。这是一个两难的选择，如果测量的频繁，可以保证下山的方向是绝对正确的，但又非常耗时，如果测量的过少，又有偏离轨道的风险。所以需要找到一个合适的测量方向的频率，来确保下山的方向不错误，同时又不至于耗时太多！

首先，我们有一个可微分的函数。这个函数就代表着一座山。我们的目标就是找到这个函数的最小值，也就是山底。根据之前的场景假设，最快的下山的方式就是找到当前位置最陡峭的方向，然后沿着此方向向下走，对应到函数中，就是找到给定点的梯度 ，然后朝着梯度相反的方向，就能让函数值下降的最快！因为梯度的方向就是函数之变化最快的方向
所以，我们重复利用这个方法，反复求取梯度，最后就能到达局部的最小值，这就类似于我们下山的过程。而求取梯度就确定了最陡峭的方向，也就是场景中测量方向的手段。那么为什么梯度的方向就是最陡峭的方向呢？

2.1，微分

看待微分的意义，可以有不同的角度，最常用的两种是：

• 函数图像中，某点的切线的斜率

• 函数的变化率
  几个微分的例子：

例子都是单变量的微分，当一个函数有多个变量的时候，就有了多变量的微分，即分别对每个变量进行求微分

2.2，梯度

梯度实际上就是多变量微分的一般化。
下面这个例子：

我们可以看到，梯度就是分别对每个变量进行微分，然后用逗号分割开，梯度是用<>包括起来，说明梯度其实一个向量。

梯度是微积分中一个很重要的概念，之前提到过梯度的意义

• 在单变量的函数中，梯度其实就是函数的微分，代表着函数在某个给定点的切线的斜率
• 在多变量函数中，梯度是一个向量，向量有方向，梯度的方向就指出了函数在给定点的上升最快的方向

这也就说明了为什么我们需要千方百计的求取梯度！我们需要到达山底，就需要在每一步观测到此时最陡峭的地方，梯度就恰巧告诉了我们这个方向。梯度的方向是函数在给定点上升最快的方向，那么梯度的反方向就是函数在给定点下降最快的方向，这正是我们所需要的。所以我们只要沿着梯度的方向一直走，就能走到局部的最低点！

2.3，梯度下降算法的数学解释

上面我们花了大量的篇幅介绍梯度下降算法的基本思想和场景假设，以及梯度的概念和思想。下面我们就开始从数学上解释梯度下降算法的计算过程和思想！

此公式的意义是：J是关于Θ的一个函数，我们当前所处的位置为Θ0点，要从这个点走到J的最小值点，也就是山底。首先我们先确定前进的方向，也就是梯度的反向，然后走一段距离的步长，也就是α，走完这个段步长，就到达了Θ1这个点！

下面就这个公式的几个常见的疑问：

• α是什么含义？
  α在梯度下降算法中被称作为学习率或者步长，意味着我们可以通过α来控制每一步走的距离，以保证不要步子跨的太大扯着蛋，哈哈，其实就是不要走太快，错过了最低点。同时也要保证不要走的太慢，导致太阳下山了，还没有走到山下。所以α的选择在梯度下降法中往往是很重要的！α不能太大也不能太小，太小的话，可能导致迟迟走不到最低点，太大的话，会导致错过最低点！

• 为什么要梯度要乘以一个负号？
  梯度前加一个负号，就意味着朝着梯度相反的方向前进！我们在前文提到，梯度的方向实际就是函数在此点上升最快的方向！而我们需要朝着下降最快的方向走，自然就是负的梯度的方向，所以此处需要加上负号

2.4，梯度下降算法的实例

我们已经基本了解了梯度下降算法的计算过程，那么我们就来看几个梯度下降算法的小实例，首先从单变量的函数开始

（1）单变量函数的梯度下降

我们假设有一个单变量的函数

函数的微分

初始化，起点为

学习率为

根据梯度下降的计算公式

我们开始进行梯度下降的迭代计算过程：

如图，经过四次的运算，也就是走了四步，基本就抵达了函数的最低点，也就是山底

（2）多变量函数的梯度下降

我们假设有一个目标函数

现在要通过梯度下降法计算这个函数的最小值。我们通过观察就能发现最小值其实就是 (0，0)点。但是接下来，我们会从梯度下降算法开始一步步计算到这个最小值！
我们假设初始的起点为：

初始的学习率为：

函数的梯度为：

进行多次迭代：

我们发现，已经基本靠近函数的最小值点

2.5 梯度下降算法的实现

下面我们将用python实现一个简单的梯度下降算法。场景是一个简单的线性回归的例子：假设现在我们有一系列的点，如下图所示

我们将用梯度下降法来拟合出这条直线！

首先，我们需要定义一个代价函数，在此我们选用均方误差代价函数

此公示中

• m是数据集中点的个数
• ½是一个常量，这样是为了在求梯度的时候，二次方乘下来就和这里的½抵消了，自然就没有多余的常数系数，方便后续的计算，同时对结果不会有影响
• y 是数据集中每个点的真实y坐标的值
• h 是我们的预测函数，根据每一个输入x，根据Θ 计算得到预测的y值，即

我们可以根据代价函数看到，代价函数中的变量有两个，所以是一个多变量的梯度下降问题，求解出代价函数的梯度，也就是分别对两个变量进行微分

明确了代价函数和梯度，以及预测的函数形式。我们就可以开始编写代码了。但在这之前，需要说明一点，就是为了方便代码的编写，我们会将所有的公式都转换为矩阵的形式，python中计算矩阵是非常方便的，同时代码也会变得非常的简洁。

为了转换为矩阵的计算，我们观察到预测函数的形式

我们有两个变量，为了对这个公式进行矩阵化，我们可以给每一个点x增加一维，这一维的值固定为1，这一维将会乘到Θ0上。这样就方便我们统一矩阵化的计算

然后我们将代价函数和梯度转化为矩阵向量相乘的形式

2.6 coding time

首先，我们来看看这个算法的快速实现

import numpy as np
X = 2 * np.random.rand(100,1)
y = 4 + 3 * X + np.random.randn(100,1)
X_b = np.c_[np.ones((100, 1)), X]
eta = 0.1#学习率
n_iterations = 1000#迭代次数
m = 100#样本数
theta = np.random.randn(2,1)

for iteration in range(n_iterations):
    gradients = 2/m * X_b.T.dot(X_b.dot(theta) - y)
    theta = theta - eta * gradients
print(theta)

计算结果：

2.7 随机梯度下降

批量梯度下降的主要问题是它要用整个训练集来计算每一步的梯度，所以训练集很大时，算法会特别慢。

由于批梯度下降每跟新一个参数的时候，要用到所有的样本数，所以训练速度会随着样本数量的增加而变得非常缓慢。随机梯度下降正是为了解决这个办法而提出的。它是利用每个样本的损失函数对θ求偏导得到对应的梯度，来更新θ：

随机梯度下降是通过每个样本来迭代更新一次，对比上面的批量梯度下降，迭代一次需要用到所有训练样本（往往如今真实问题训练数据都是非常巨大），一次迭代不可能最优，如果迭代10次的话就需要遍历训练样本10次。但是，SGD伴随的一个问题是噪音较BGD要多，使得SGD并不是每次迭代都向着整体最优化方向。

我们可以从图中看出SGD迭代的次数较多，在解空间的搜索过程看起来很盲目。但是大体上是往着最优值方向移动。

import numpy as np
import matplotlib.pyplot as plt
X = 2 * np.random.rand(100,1)
y = 4 + 3 * X + np.random.randn(100,1)
X_b = np.c_[np.ones((100, 1)), X]
X_new = np.array([[0], [2]])
X_new_b = np.c_[np.ones((2, 1)), X_new]
theta_path_sgd = []
m = len(X_b)
np.random.seed(42)
n_epochs = 50
t0, t1 = 5, 50  # learning schedule hyperparameters

def learning_schedule(t):
    return t0 / (t + t1)

theta = np.random.randn(2,1)  # random initialization

for epoch in range(n_epochs):
    for i in range(m):
        if epoch == 0 and i < 20:                    # not shown in the book
            y_predict = X_new_b.dot(theta)           # not shown
            style = "b-" if i > 0 else "r--"         # not shown
            plt.plot(X_new, y_predict, style)        # not shown
        random_index = np.random.randint(m)
        xi = X_b[random_index:random_index+1]
        yi = y[random_index:random_index+1]
        gradients = 2 * xi.T.dot(xi.dot(theta) - yi)
        eta = learning_schedule(epoch * m + i)
        theta = theta - eta * gradients
        theta_path_sgd.append(theta)                 # not shown
print(theta)
plt.plot(X, y, "b.")                                 # not shown
plt.xlabel("$x_1$", fontsize=18)                     # not shown
plt.ylabel("$y$", rotation=0, fontsize=18)           # not shown
plt.axis([0, 2, 0, 15])                              # not shown
                          # not shown
plt.show()

结果如下：

使用Scikit-Leaarn实现随机梯度下降

import numpy as np
import matplotlib.pyplot as plt
X = 2 * np.random.rand(100,1)
y = 4 + 3 * X + np.random.randn(100,1)
from sklearn.linear_model import SGDRegressor
sgd_reg = SGDRegressor(max_iter=50, tol=-np.infty, penalty=None, eta0=0.1, random_state=42)
sgd_reg.fit(X, y.ravel())
print(sgd_reg.intercept_, sgd_reg.coef_)

结果为：

2.8 小批量梯度下降

小批量梯度下降，是对批量梯度下降以及随机梯度下降的一个折中办法。其思想是：每次迭代使用 batch_size个样本来对参数进行更新。
优点：
  （1）通过矩阵运算，每次在一个batch上优化神经网络参数并不会比单个数据慢太多。
  （2）每次使用一个batch可以大大减小收敛所需要的迭代次数，同时可以使收敛到的结果更加接近梯度下降的效果。(比如上例中的30W，设置batch_size=100时，需要迭代3000次，远小于SGD的30W次)
  （3）可实现并行化。
缺点：
  （1）batch_size的不当选择可能会带来一些问题。
batcha_size的选择带来的影响：
  （1）在合理地范围内，增大batch_size的好处：
    a. 内存利用率提高了，大矩阵乘法的并行化效率提高。
    b. 跑完一次 epoch（全数据集）所需的迭代次数减少，对于相同数据量的处理速度进一步加快。
    c. 在一定范围内，一般来说 Batch_Size 越大，其确定的下降方向越准，引起训练震荡越小。
  （2）盲目增大batch_size的坏处：
    a. 内存利用率提高了，但是内存容量可能撑不住了。
    b. 跑完一次 epoch（全数据集）所需的迭代次数减少，要想达到相同的精度，其所花费的时间大大增加了，从而对参数的修正也就显得更加缓慢。
    c. Batch_Size 增大到一定程度，其确定的下降方向已经基本不再变化。

import numpy as np
import matplotlib.pyplot as plt
X = 2 * np.random.rand(100,1)
y = 4 + 3 * X + np.random.randn(100,1)
m = 100
X_b = np.c_[np.ones((100, 1)), X]
theta_path_mgd = []
n_iterations = 50
minibatch_size = 20
np.random.seed(42)
theta = np.random.randn(2,1)  # random initialization
t0, t1 = 200, 1000
def learning_schedule(t):
    return t0 / (t + t1)

t = 0
for epoch in range(n_iterations):
    shuffled_indices = np.random.permutation(m)
    X_b_shuffled = X_b[shuffled_indices]
    y_shuffled = y[shuffled_indices]
    for i in range(0, m, minibatch_size):
        t += 1
        xi = X_b_shuffled[i:i+minibatch_size]
        yi = y_shuffled[i:i+minibatch_size]
        gradients = 2/minibatch_size * xi.T.dot(xi.dot(theta) - yi)
        eta = learning_schedule(t)
        theta = theta - eta * gradients
        theta_path_mgd.append(theta)
print(theta)

结果为：

下图显示了三种梯度下降算法的收敛过程：

三，多项式回归

如果数据比简单的直线更为复杂，我们还想用线性回归拟合非线性回归数据，一个简单的方法就是将每个特征的幂次方添加为一个新特征，然后在这个拓展过的特征集上训练线性模型。

import numpy as np
import numpy.random as rnd
import matplotlib.pyplot as plt
np.random.seed(42)
m = 100
X = 6 * np.random.rand(m, 1) - 3
y = 0.5 * X**2 + X + 2 + np.random.randn(m, 1)
plt.plot(X, y, "b.")
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.axis([-3, 3, 0, 10])
plt.show()

结果为：

显然，直线永远不可能拟合这个数据。所以我们使用Scikit-Learn的PolynomialFeatures类来对训练数据进行转换，将每个特征的平方（二次多项式）作为新特征加入训练集。

from sklearn.preprocessing import PolynomialFeatures
poly_features = PolynomialFeatures(degree=2, include_bias=False)
X_poly = poly_features.fit_transform(X)
print(X[0])
print(X_poly[0])

结果为：

X_poly现在包含原本的特征X和该特征的平方。现在对这个扩展后的训练集匹配一个LinearRegression模型。

import numpy as np
import numpy.random as rnd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

np.random.seed(42)
m = 100
X = 6 * np.random.rand(m, 1) - 3
y = 0.5 * X**2 + X + 2 + np.random.randn(m, 1)
from sklearn.preprocessing import PolynomialFeatures
poly_features = PolynomialFeatures(degree=2, include_bias=False)
X_poly = poly_features.fit_transform(X)
lin_reg = LinearRegression()
lin_reg.fit(X_poly, y)
print(lin_reg.intercept_, lin_reg.coef_)

X_new=np.linspace(-3, 3, 100).reshape(100, 1)
X_new_poly = poly_features.transform(X_new)
y_new = lin_reg.predict(X_new_poly)
plt.plot(X, y, "b.")
plt.plot(X_new, y_new, "r-", linewidth=2, label="Predictions")
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.legend(loc="upper left", fontsize=14)
plt.axis([-3, 3, 0, 10])
plt.show()

结果如图：

还不错，模型预估为：，而实际上原来的函数时
注意，当存在多个特征时，多项式回归能够发现特征与特征之间的关系。这是因为PolynomialFeatures会在给定多项式阶数下，添加所有特征组合。例如，有两个特征a和b,阶数degree=3,PolynomialFeatures不只会添加特征a^2、b2、a^3和b3,还会添加组合ab 、a^2*b以及a*b2.

四，学习曲线

高阶多项式回归对训练数据的拟合，很可能会比简单的线性回归要好。例如，下图使用了一个300阶多项式模型来处理训练数据，并将结果与一个纯线性模型和一个二次模型进行对比。

这个高阶多项式回归模型严重的过度拟合了训练数据，而线性模型则是拟合不足。这个案例中泛化结果最好的是二次模型。那么，应该如何确定模型的复杂程度？怎么样判断模型是欠拟合还是过拟合。

前面我们使用交叉验证来评估模型的泛化性能。如果模型在训练集上表现良好，但是交叉验证的泛化表现非常糟糕，那么模型就是过度拟合。如果在二者上表现都不佳，那就是拟合不足。

还有一种方法是观察学习曲线，这个曲线绘制的是模型在训练集和验证集上，关于“训练集大小”的性能函数。要生成这个曲线，只需要在不同大小的训练子集上多次训练模型即可。

import numpy as np
import numpy.random as rnd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

np.random.seed(42)
m = 100
X = 6 * np.random.rand(m, 1) - 3
y = 0.5 * X**2 + X + 2 + np.random.randn(m, 1)
from sklearn.preprocessing import PolynomialFeatures
poly_features = PolynomialFeatures(degree=2, include_bias=False)
X_poly = poly_features.fit_transform(X)
lin_reg = LinearRegression()
lin_reg.fit(X_poly, y)
print(lin_reg.intercept_, lin_reg.coef_)

X_new=np.linspace(-3, 3, 100).reshape(100, 1)
X_new_poly = poly_features.transform(X_new)
y_new = lin_reg.predict(X_new_poly)


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

def plot_learning_curves(model, X, y):
    X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, random_state=10)
    train_errors, val_errors = [], []
    for m in range(1, len(X_train)):
        model.fit(X_train[:m], y_train[:m])
        y_train_predict = model.predict(X_train[:m])
        y_val_predict = model.predict(X_val)
        train_errors.append(mean_squared_error(y_train[:m], y_train_predict))
        val_errors.append(mean_squared_error(y_val, y_val_predict))

    plt.plot(np.sqrt(train_errors), "r-+", linewidth=2, label="train")
    plt.plot(np.sqrt(val_errors), "b-", linewidth=3, label="val")
    plt.legend(loc="upper right", fontsize=14)   # not shown in the book
    plt.xlabel("Training set size", fontsize=14) # not shown
    plt.ylabel("RMSE", fontsize=14)              # not shown

lin_reg = LinearRegression()
plot_learning_curves(lin_reg, X, y)
plt.axis([0, 80, 0, 3])                         # not shown in the book
plt.show()         

结果如下：

首先，我们来看训练数据上的性能，当模型只包含一两个实例时，模型可以完美拟合，这是为什么曲线是从0开始的，但是随着新的实例被添加进训练集中，模型不能完美拟合训练数据了，因为数据有噪声，并且根本就不是线性的。所以训练集的误差一路上升，直到抵达一个高地，从这一点开始，添加新实例到训练集中不再使平均误差上升或下降。然后我们再来看看验证集的性能表现。当训练集实例非常少时，模型不能很好的泛化，这是为什么验证集误差值一开始非常大，随着模型经历更多的训练数据，它开始学习，因此验证集误差慢慢下降。

这条学习曲线时典型的模型拟合不足。两条曲线均到达高地，非常接近，而且相当高。

下面我们看看在同样数据集上，一个10阶多项式模型的学习曲线：

from sklearn.pipeline import Pipeline
polynomial_regression = Pipeline([
        ("poly_features", PolynomialFeatures(degree=10, include_bias=False)),
        ("lin_reg", LinearRegression()),
    ])

plot_learning_curves(polynomial_regression, X, y)
plt.axis([0, 80, 0, 3])           # not shown
plt.show()

结果如图：

训练数据误差远低于线性回归模型，两条曲线之间有一定差距。这意味着该模型在训练数据上表现比验证集上要好很多，这是过度拟合的标志。

但是，如果你使用更大的训练集，那么这两条曲线会越来越近，所以改进过度拟合的方法之一就是提供更多的训练数据，直到验证误差接近训练误差。

五，正则线性模型

据集上，一个10阶多项式模型的学习曲线：

from sklearn.pipeline import Pipeline
polynomial_regression = Pipeline([
        ("poly_features", PolynomialFeatures(degree=10, include_bias=False)),
        ("lin_reg", LinearRegression()),
    ])

plot_learning_curves(polynomial_regression, X, y)
plt.axis([0, 80, 0, 3])           # not shown
plt.show()

结果如图：

[外链图片转存中…(img-zwIObxfj-1567400052040)]

训练数据误差远低于线性回归模型，两条曲线之间有一定差距。这意味着该模型在训练数据上表现比验证集上要好很多，这是过度拟合的标志。

但是，如果你使用更大的训练集，那么这两条曲线会越来越近，所以改进过度拟合的方法之一就是提供更多的训练数据，直到验证误差接近训练误差。