机器学习实战_线性回归&逻辑回归（一）

一些基础概念就随意带过了啦~~
线性回归预测模型：$\theta_0$是偏置项bias

$$\hat{y} = \theta _{0} + \theta _{1}x _{1}+\theta _{2}x _{2}+\dots+\theta _{n}x _{n}$$

上述公式可以写成更为简洁的向量形式:

$$\hat{y} = h _{\theta} (\mathbf{x})= \theta^T \cdot \mathbf{x}$$

• $\theta^T$表示向量$\theta$ 的转置（行向量变为了列向量）
• $\mathbf{x}$ 为每个样例中特征值的向量形式，包括 $x _{1}$到 $x_{n}$ ，而且$x_0$恒为1
• $\theta^T \cdot \mathbf{x}$ 表示$\theta^T$ 和 $\mathbf{x}$ 的点积
• $h_{\theta}$ 表示参数为$\theta$ 的假设函数

损失函数：以后的使用中为了公式的简洁，使用$MSE(\theta)$ 来代替$MSE (\mathbf{X},h _{\theta}) $

$$MSE (\mathbf{X},h _{\theta}) = \frac{1}{m} \sum\limits_{i=1}^m{\left(\theta^T \cdot \mathbf{x}^{(i)}-y^{(i)}\right)}^2$$

如何求解这个损失函数呢？
使均方误差最小化来进行模型求解的方法称为“最小二乘法”。在线性回归中，最小二乘法就是试图找到一条直线，使让本到直线上的欧氏距离之和最小。

推倒过程：

所以最小化损失可以写成（正态方程求解）：
$$\hat{\theta} = ({\mathbf{X}}^T\cdot\mathbf{X})^{-1}\cdot{\mathbf{X}}^T\cdot\mathbf{y}$$

这里的$X^TX$要满足满秩矩阵，然而，现实大多数任务不会满足这个条件。此时可解除多个θ，选择哪一个由学习算法的偏好决定，常见的做法是引入正则化。

让我们生成一些近似线性的数据来测试一下这个方程。

import numpy as np 
X = 2 * np.random.rand(100, 1)    # 产生100个[0,2]的数，列表示
y = 4 + 3 * X + np.random.randn(100, 1)

现在让我们使用上面的矩阵方程来计算$\hat{\theta}$ ，我们将使用Numpy的线性代数模块(np.linalg)中的inv()函数来计算矩阵的逆，以及dot()方法来计算矩阵的乘法。

# np.c_表示按列操作拼接，np.r_表示按行操作拼接
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)

我们生产数据的函数实际上是y = 4 + 3x_0 + 高斯噪声。让我们看一下最后的计算结果。

>>> theta_best
array([[4.21509616],[2.77011339]])    # 理论上为4,3；由于存在噪声，参数不可能达到到原始函数的值

利用$\hat{\theta}$ 进行预测

>>> X_new = np.array([[0],[2]])    #两个测试样本分别x=0，x=2
>>> X_new_b = np.c_[np.ones((2, 1)), X_new]     # 一行为一个样本x1=[1,0],x2=[1,2]
>>> y_predict = X_new_b.dot(theta.best)
>>> y_predict
array([[4.21509616],[9.75532293]])

使用下面的Scikit-Learn代码可以达到相同的效果：

# lin_reg.intercept_获得截距bias, lin_reg.coef_获得模型参数
>>> form sklearn.linear_model import LinearRegression
>>> lin_reg = LinearRegression()
>>> lin_reg.fit(X,y)
>>> lin_reg.intercept_, lin_reg.coef_
(array([4.21509616]),array([2.77011339]))
>>> lin_reg.predict(X_new)
array([[4.21509616],[9.75532293]])

这样的矩阵复杂度较高。同时，一旦你得到了线性回归模型（通过解正态方程或者其他的算法），进行预测是非常快的。求矩阵的逆是很费时间的，尤其是特征数量很多的情况下就更糟糕了。

梯度下降(推导和最小二乘法一样，得到偏导结果即可)

当我们使用梯度下降的时候，应该确保所有的特征有着相近的尺度范围（例如：使用Scikit_Learn的 StandardScaler类），否则它将需要很长的时间成能够收敛。

批量梯度下降：使用梯度下降的过程中，你需要计算每一个$\theta_j$ 下代价函数的梯度
代价函数的偏导数:
$$\frac{\partial }{\partial \theta_j}MSE(\theta)=\frac{2}{m} \sum\limits_{i=1}^m{\left(\theta^T \cdot \mathbf{x}^{(i)}-y^{(i)}\right)}{x_j}^{(i)}$$

为了避免单独计算每一个梯度，你也可以使用下面的公式来一起计算它们。梯度向量记为$\nabla_{\theta}MSE(\theta)$ ，其包含了代价函数所有的偏导数(每个模型参数只出现一次)。

代价函数的梯度向量：

在这个方程中每一步计算时都包含了整个训练集$\mathbf{X}$ ，这也是为什么这个算法称为批量梯度下降：每一次训练过程都使用所有的的训练数据。因此，在大数据集上，其会变得相当的慢（但是我们接下来将会介绍更快的梯度下降算法）。然而，梯度下降的运算规模和特征的数量成正比。训练一个数千数量特征的线性回归模型使用梯度下降要比使用正态方程快的多.

最后更新梯度步长：($\eta$ 为学习速率)
$$\theta^{(next\ step)}=\theta - \eta\nabla_{\theta}MSE(\theta)$$

让我们看一下这个算法的应用

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 * gradiens    # 最后输出的是所有系数矩阵

让我们看一下最后的结果：

>>> theta
array([[4.21509616],[2.77011339]])

不同学习率的迭代结果：

为了找到一个好的学习率，你可以使用网格搜索（详见第二章）。当然，你一般会限制迭代的次数，以便网格搜索可以消除模型需要很长时间才能收敛这一个问题。一个简单的解决方法是：设置一个非常大的迭代次数，但是当梯度向量变得非常小的时候，结束迭代。非常小指的是：梯度向量小于一个值$\varepsilon$ （称为容差）。这时候可以认为梯度下降几乎已经达到了最小值。

随机梯度下降：能跳出局部最优，但不一定是全局最优

却保所有的特征有着相近的尺度范围；却保所有的特征有着相近的尺度范围；却保所有的特征有着相近的尺度范围；（再强调一遍）

批量梯度下降的最要问题是计算每一步的梯度时都需要使用整个训练集，这导致在规模较大的数据集上，其会变得非常的慢。与其完全相反的随机梯度下降，在每一步的梯度计算上只随机选取训练集中的一个样例。

虽然随机性可以很好的跳过局部最优值，但同时它却不能达到最小值。解决这个难题的一个办法是逐渐降低学习率。 开始时，走的每一步较大（这有助于快速前进同时跳过局部最小值），然后变得越来越小，从而使算法到达全局最小值。 这个过程被称为模拟退火，因为它类似于熔融金属慢慢冷却的冶金学退火过程。 决定每次迭代的学习率的函数称为learning schedule。 如果学习速度过快的降低，你可能会陷入局部最小值，甚至在到达最小值的半路就停止了。 如果学习速度降低得太慢，你可能在最小值的附近长时间摆动，同时如果过早停止训练，最终只会出现次优解

下面的代码使用一个简单的learning schedule来实现随机梯度下降：

n_epochs = 50 
t0, t1 = 5, 50  #learning_schedule的超参数

def learning_schedule(t):
    return t0 / (t + t1)
    
# np.random.randn返回2行1列符合标准正态分布的数；np.random.rand返回[0,1）的随机数；randint返回范围内的整数
theta = np.random.randn(2,1)

for epoch in range(n_epochs):
    for i in range(m):
        random_index = np.random.randint(m)    # 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 * gradiens

>>> theta
array([[4.21076011],[2.748560791]])

由于每个实例的选择是随机的，有的实例可能在每一代中都被选到，这样其他的实例也可能一直不被选到。如果你想保证每一代迭代过程，算法可以遍历所有实例，一种方法是将训练集打乱重排，然后选择一个实例，之后再继续打乱重排，以此类推一直进行下去。但是这样收敛速度会非常的慢。

通过使用Scikit-Learn完成线性回归的随机梯度下降，你需要使用SGDRegressor类，这个类默认优化的是均方差代价函数。下面的代码迭代了50代，其学习率eta为0.1，使用默认的learning schedule(与前面的不一样)，同时也没有添加任何正则项（penalty = None）：

##因为这个函数需要的y是一个行向量，所以压扁;另外，numpy.flatten() 与 numpy.ravel()将多维数组降位一维，前者会进行拷贝处理
from sklearn.linear_model import SGDRegressor
sgd_reg = SGDRregressor(n_iter=50, penalty=None, eta=0.1)
sgd_reg.fit(X,y.ravel()) 

结果很接近

>>> sgd_reg.intercept_, sgd_reg.coef_
(array([4.18380366]),array([2.74205299]))

小批量梯度下降：小批量梯度下降中，它则使用一个随机的小型实例集，小批量梯度下降在参数空间上的表现比随机梯度下降要好的多，尤其在有大量的小型实例集时，主要利用了矩阵运算的硬件优化

下图显示了训练期间三种梯度下降算法在参数空间中所采用的路径。 他们都接近最小值，但批量梯度的路径最后停在了最小值，而随机梯度和小批量梯度最后都在最小值附近摆动。 但是，不要忘记，批次梯度需要花费大量时间来完成每一步，但是，如果你使用了一个较好的learning schedule，随机梯度和小批量梯度也可以得到最小值。

多项式回归

如果你的数据实际上比简单的直线更复杂呢？ 令人惊讶的是，你依然可以使用线性模型来拟合非线性数据。 一个简单的方法是对每个特征进行加权后作为新的特征，然后训练一个线性模型在这个扩展的特征集。 这种方法称为多项式回归。

首先，我们根据一个简单的二次方程（并加上一些噪声，如图）来生成一些非线性数据

m = 100
X = 6 * np.random.rand(m, 1) - 3
y = 0.5 * X**2 + X + 2 + np.random.randn(m, 1)

于是，我们使用Scikit-Learning的PolynomialFeatures类进行训练数据集的转换，让训练集中每个特征的平方（2次多项式）作为新特征（在这种情况下，仅存在一个特征）：

>>> from sklearn.preprocessing import PolynomialFeatures
>>> poly_features = PolynomialFeatures(degree=2,include_bias=False)
>>> X_poly = poly_features.fit_transform(X)    # 转换特征，包含原始特征和二次项特征
>>> X[0]
array([-0.75275929])
>>> X_poly[0]
array([-0.75275929, 0.56664654])

现在包含原始特X并加上了这个特征的平方X^2。现在你可以在这个扩展训练集上使用LinearRegression模型进行拟合

>>> lin_reg = LinearRegression()
>>> lin_reg.fit(X_poly, y)
>>> lin_reg.intercept_, lin_reg.coef_
(array([ 1.78134581]), array([[ 0.93366893, 0.56456263]]))

还是不错的，模型预测函数$\hat{y}=0.56x_1^2+0.93x_1+1.78$ ，事实上原始函数为$y=0.5x_1^2+1.0x_1+2.0$ 再加上一些高斯噪声。

请注意，当存在多个特征时，多项式回归能够找出特征之间的关系（这是普通线性回归模型无法做到的）。 这是因为LinearRegression会自动添加当前阶数下特征的所有组合。例如，如果有两个特征a,b，使用3阶（degree=3）的LinearRegression时，不仅仅只有$a^2,a^3,b^2$,同时也会有它们的其他组合项$ab,a^2b,ab^2$。

• LinearRegression(degree=d)把一个包含n个特征的数组转换为一个包含$\frac{(n+d)!}{d!n!}$ 特征的数组，n!表示n的阶乘，等于$1 * 2 * 3 \cdots * n$ 。小心大量特征的组合爆炸！*

学习曲线

如果你使用一个高阶的多项式回归，你可能发现它的拟合程度要比普通的线性回归要好的多。例如，使用一个300阶的多项式模型去拟合之前的数据集，并同简单线性回归、2阶的多项式回归进行比较

当然，这种高阶多项式回归模型在这个训练集上严重过拟合了，线性模型则欠拟合。在这个训练集上，二次模型有着较好的泛化能力。那是因为数据的产生是使用了二次模型，但是一般我们不知道这个数据生成函数是什么，那我们该如何决定我们模型的复杂度呢？你如何告诉我你的模型是过拟合还是欠拟合？

之前，我们可以使用交叉验证来估计一个模型的泛化能力。如果一个模型在训练集上表现良好，通过交叉验证指标却得出其泛化能力很差，那么你的模型就是过拟合了。如果在这两方面都表现不好，那么它就是欠拟合了。这种方法可以告诉我们，你的模型是太复杂了还是太简单了。

另一种方法是观察学习曲线：画出模型在训练集上的表现，同时画出以训练集规模为自变量的训练集函数。为了得到图像，需要在训练集的不同规模子集上进行多次训练。下面的代码定义了一个函数，用来画出给定训练集后的模型学习曲线：

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)
    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_predict, y_train[:m]))
        val_errors.append(mean_squared_error(y_val_predict, y_val))
plt.plot(np.sqrt(train_errors), "r-+", linewidth=2, label="train")    # 训练损失
plt.plot(np.sqrt(val_errors), "b-", linewidth=3, label="val")    # 验证损失

我们一起看一下简单线性回归模型的学习曲线

lin_reg = LinearRegression()
plot_learning_curves(lin_reg, X, y)

上面的曲线表现了一个典型的欠拟合模型，两条曲线都到达高原地带并趋于稳定，并且最后两条曲线非常接近，同时误差值非常大。

如果你的模型在训练集上是欠拟合的，添加更多的样例是没用的。你需要使用一个更复杂的模型或者找到更好的特征。

现在让我们看一个在相同数据上10阶多项式模型拟合的学习曲线

from sklearn.pipeline import Pipeline

polynomial_regression = Pipeline((
    ("poly_features", PolynomialFeatures(degree=10, include_bias=False)),
    ("sgd_reg", LinearRegression()),
))

plot_learning_curves(polynomial_regression, X, y)

• 在训练集上，误差要比线性回归模型低的多。
• 图中的两条曲线之间有间隔，这意味模型在训练集上的表现要比验证集上好的多，这也是模型过拟合的显著特点。当然，如果你使用了更大的训练数据，这两条曲线最后会非常的接近。

改善模型过拟合的一种方法是提供更多的训练数据，直到训练误差和验证误差相等

在统计和机器学习领域有个重要的理论：一个模型的泛化误差由三个不同误差的和决定：

• 偏差：泛化误差的这部分误差是由于错误的假设决定的。例如实际是一个二次模型，你却假设了一个线性模型。一个高偏差的模型最容易出现欠拟合。
• 方差：这部分误差是由于模型对训练数据的微小变化较为敏感，一个多自由度的模型更容易有高的方差（例如一个高阶多项式模型），因此会导致模型过拟合。
• 不可约误差：这部分误差是由于数据本身的噪声决定的。降低这部分误差的唯一方法就是进行数据清洗（例如：修复数据源，修复坏的传感器，识别和剔除异常值）。

下图依次为欠拟合，过拟合，较合适。