《Hands-On Machine Learning with Scikit-Learn & TensorFlow》读书笔记 第四章 训练模型（回归）

第4章 训练模型（回归）

首先我们将以一个简单的线性回归模型为例，讨论两种不同的训练方法来得到模型的最优解：

• 直接使用封闭方程进行求根运算，得到模型在当前训练集上的最优参数（即在训练集上使损失函数达到最小值的模型参数）

• 使用迭代优化方法：梯度下降（GD），在训练集上，它可以逐渐调整模型参数以获得最小的损失函数，最终，参数会收敛到和第一种方法相同的的值。同时，我们也会介绍一些梯度下降的变体形式：批量梯度下降（Batch GD）、小批量梯度下降（Mini-batch GD）、随机梯度下降（Stochastic GD），在第二部分的神经网络部分，我们会多次使用它们。

接下来，我们将研究一个更复杂的模型：多项式回归，它可以拟合非线性数据集，由于它比线性模型拥有更多的参数，于是它更容易出现模型的过拟合。因此，我们将介绍如何通过学习曲线去判断模型是否出现了过拟合，并介绍几种正则化方法以减少模型出现过拟合的风险。

最后，我们将介绍两个常用于分类的模型：Logistic回归和Softmax回归

---

线性回归

线性模型更一般化的描述指通过计算输入变量的加权和，并加上一个常数偏置项（截距项）来得到一个预测值。如公式 4-1：

公式 4-1：线性回归预测模型

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

$\hat{y}$表示预测结果
$n$ 表示特征的个数
$x_{i}$ 表示第$i$个特征的值
$\theta_{j}$ 表示第$j$个参数（包括偏置项 $\theta _{0}$ 和特征权重值$\theta _{1},\theta _{2},\dots,\theta _{n}$）

上述公式可以写成更为简洁的向量形式，如公式 4-2：
公式 4-2：线性回归预测模型（向量形式）

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

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

在训练集 $\mathbf{X}$上使用公式 4-3 来计算线性回归假设 $h_{\theta}$的均方差（MSE）。

公式 4-3：线性回归模型的 MSE 损失函数

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

为了找到最小化损失函数的 $\theta$ 值，可以采用公式解，换句话说，就是可以通过解正态方程直接得到最后的结果。

公式 4-4：正态方程

$\hat{\theta} = ({\mathbf{X}}^T\cdot\mathbf{X})^{-1}\cdot{\mathbf{X}}^T\cdot\mathbf{y}$

$\hat{\theta}$指最小化损失$\theta$ 的值
$\mathbf{y}$ 是一个向量，其包含了 $y^{(1)}$ 到 $y^{(m)}$的值

---

代码实现

创建数据，y = 4 + 3x + 随机噪声

import numpy as np

X = 2 * np.random.rand(100, 1) 
y = 4 + 3 * X + np.random.randn(100, 1)

1、根据公式 4-4 计算 theta 值

X_b = np.c_[np.ones((100, 1)), X] # add x0 = 1 to each instance 
theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)
>>>theta_best
array([[ 3.93182756],
       [ 3.00188366]])

与目标 theta =（4，3） 比较接近

预测当 x=0 和 x=2 时的预测值

X_new = np.array([[0], [2]])

X_new_b = np.c_[np.ones((2, 1)), X_new] # add x0 = 1 to each instance
y_predict = X_new_b.dot(theta_best)
>>>y_predict
array([[ 3.93182756],
       [ 9.93559488]])

绘制模型图像

%matplotlib inline
import matplotlib.pyplot as plt

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

计算复杂度

正态方程需要计算矩阵 ${\mathbf{X}}^T\cdot\mathbf{X}$的逆，它是一个 $n * n$ 的矩阵（n 是特征的个数）。这样一个矩阵求逆的运算复杂度大约在 $O(n^{2.4})$ 到 $O(n^3)$ 之间，具体值取决于计算方式。换句话说，如果你将你的特征个数翻倍的话，其计算时间大概会变为原来的 5.3$（2^{2.4}$）到 8$（2^3）$倍。当特征的个数较大的时候（例如：特征数量为 100000），正态方程求解将会非常慢。

有利的一面是，这个方程在训练集上对于每一个实例来说是线性的，其复杂度为 $O(m)$，因此只要有能放得下它的内存空间，它就可以对大规模数据进行训练。同时，一旦你得到了线性回归模型（通过解正态方程或者其他的算法），进行预测是非常快的。因为模型中计算复杂度对于要进行预测的实例数量和特征个数都是线性的。 换句话说，当实例个数变为原来的两倍多的时候（或特征个数变为原来的两倍多），预测时间也仅仅是原来的两倍多。

也可以使用 sklearn 完成上述过程

from sklearn.linear_model import LinearRegression

lin_reg = LinearRegression()
lin_reg.fit(X, y)
>>>lin_reg.intercept_, lin_reg.coef_
(array([ 3.93182756]), array([[ 3.00188366]]))

接下来，我们将介绍另一种方法去训练模型。这种方法适合在特征个数非常多，训练实例非常多，内存无法满足要求的时候使用。

2、梯度下降法求 theta 值
为方便计算运行时间，定义一个装饰器

from functools import wraps
import time

def func_timer(function):
    # 用装饰器实现函数计时
    def function_timer(*args, **kwargs):
        t0 = time.time()
        result = function(*args, **kwargs)
        t1 = time.time()
        print('[Function: {name} finished, spent time: {time:.8f}s]'.format(name = function.__name__,time = t1 - t0))
        return result
    return function_timer

3、Batch Gradient Descent

@func_timer
def batchGD(X, y, eta = 0.1, n_iterations = 1000):  # eta: learning rate
    [m,n] = X.shape
    X_b = np.c_[np.ones((m, 1)), X]  # add x0 = 1 to each instance 
    theta = np.random.randn((n+1),1)   # random initialization
    for iteration in range(n_iterations):
        gradients = 2/m * X_b.T.dot(X_b.dot(theta) - y) 
        theta = theta - eta * gradients
    return theta
>>>batchGD(X,y,0.1,10000)
[Function: batchGD finished, spent time: 0.06436276s]
array([[ 3.93182756],
       [ 3.00188366]])

Stochastic Gradient Descent
当损失函数很不规则时，随机梯度下降算法能够跳过局部最小值。因此，随机梯度下降在寻找全局最小值上比批量梯度下降表现要好。

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

t0, t1 = 5, 50
def learning_schedule(t):
    return t0 / (t + t1)

@func_timer
def stochGD(X, y, n_epochs = 500):
    [m,n] = X.shape
    X_b = np.c_[np.ones((m, 1)), X] 
    theta = np.random.randn((n+1),1)
    for epoch in range(n_epochs):# random initialization
        for i in range(m):
            random_index = np.random.randint(m) 
            xi = X_b[random_index:random_index+1] # shape = (1,2)
            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
    return theta
>>>stochGD(X, y)
[Function: stochGD finished, spent time: 0.48928499s]
array([[ 3.94101146],
       [ 3.0024468 ]])

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

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

from sklearn.linear_model import SGDRegressor 
sgd_reg = SGDRegressor(n_iter=50, penalty=None, eta0=0.1) 
sgd_reg.fit(X, y.ravel())
>>>sgd_reg.intercept_, sgd_reg.coef
(array([ 3.95507846]), array([ 3.02741407]))

4、Mini-batch Gradient Descent

在迭代的每一步，批量梯度使用整个训练集，随机梯度时候用仅仅一个实例，在小批量梯度下降中，它则使用一个随机的小型实例集。它比随机梯度的主要优点在于你可以通过矩阵运算的硬件优化得到一个较好的训练表现，尤其当你使用 GPU 进行运算的时候。

小批量梯度下降在参数空间上的表现比随机梯度下降要好的多，尤其在有大量的小型实例集时。作为结果，小批量梯度下降会比随机梯度更靠近最小值。但是，另一方面，它有可能陷在局部最小值中（在遇到局部最小值问题的情况下，和我们之前看到的线性回归不一样）。 下图显示了训练期间三种梯度下降算法在参数空间中所采用的路径。 他们都接近最小值，但批量梯度的路径最后停在了最小值，而随机梯度和小批量梯度最后都在最小值附近摆动。 但是，不要忘记，批次梯度需要花费大量时间来完成每一步，但是，如果你使用了一个较好的learning schedule，随机梯度和小批量梯度也可以得到最小值。

Table 4-1. Comparison of algorithms for Linear Regression

Algorithm	Large m	Out-of-core support	Large n	Hyperparams	Scaling required	Scikit-Learn
Normal Equation	Fast	No	Slow	0	No	LinearRegression
Batch GD	Slow	No	Fast	2	Yes	n/a
Stochastic GD	Fast	Yes	Fast	>=2	Yes	SGDRegressor
Mini-batch GD	Fast	Yes	Fast	>=2	Yes	n/a

---

多项式回归

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

让我们看一个例子。 首先，我们根据一个简单的二次方程，并加上一些噪声来生成一些非线性数据：

m1 = 100 
X1 = 6 * np.random.rand(m1, 1) - 3 
y1 = 0.5 * X1**2 + X1 + 2 + np.random.randn(m1, 1)

使用 Scikit-Learning 的PolynomialFeatures类进行训练数据集的转换，让训练集中每个特征的平方（2 次多项式）作为新特征。

from sklearn.preprocessing import PolynomialFeatures

poly_features = PolynomialFeatures(degree=2, include_bias=False)
X1_poly = poly_features.fit_transform(X1)
>>>X1[0]
array([-0.73346809])
>>>X1_poly[0]
array([-0.73346809,  0.53797544])

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

lin_reg = LinearRegression()

lin_reg.fit(X1_poly, y1)
lin_reg.intercept_, lin_reg.coef_
>>>(array([ 1.77482187]), array([[ 0.9749555 ,  0.53042132]]))
(array([ 1.77482187]), array([[ 0.9749555 ,  0.53042132]]))

绘制曲线

X2 = np.linspace(-3,3,600)
y2 = lin_reg.coef_[0][1] * X2**2 + lin_reg.coef_[0][0] * X2 + lin_reg.intercept_[0]
plt.plot(X1, y1, "b.") 
plt.plot(X2, y2, "r") 
plt.axis([-3, 3, 0, 10]) 
plt.show()

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

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

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

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

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

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, X1, y1)

当训练集只有一两个样本的时候，模型能够非常好的拟合它们，这也是为什么曲线是从零开始的原因。但是当加入了一些新的样本的时候，训练集上的拟合程度变得难以接受，出现这种情况有两个原因，一是因为数据中含有噪声，另一个是数据根本不是线性的。因此随着数据规模的增大，误差也会一直增大，直到达到高原地带并趋于稳定，在之后，继续加入新的样本，模型的平均误差不会变得更好或者更差。我们继续来看模型在验证集上的表现，当以非常少的样本去训练时，模型不能恰当的泛化，也就是为什么验证误差一开始是非常大的。当训练样本变多的到时候，模型学习的东西变多，验证误差开始缓慢的下降。但是一条直线不可能很好的拟合这些数据，因此最后误差会到达在一个高原地带并趋于稳定，最后和训练集的曲线非常接近。

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

在相同数据上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])
plt.show()

这幅图像和之前的有一点点像，但是其有两个非常重要的不同点：

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

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

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

---

线性模型的正则化

降低模型的过拟合的好方法是正则化这个模型（即限制它）：模型有越少的自由度，就越难以拟合数据。例如，正则化一个多项式模型，一个简单的方法就是减少多项式的阶数。

对于一个线性模型，正则化的典型实现就是约束模型中参数的权重。 接下来我们将介绍三种不同约束权重的方法：Ridge 回归，Lasso 回归和 Elastic Net。

岭（Ridge）回归

岭回归（也称为 Tikhonov 正则化）是线性回归的正则化版：在损失函数上直接加上一个正则项 $\alpha\sum_{i=1}^n\theta_i^2$。这使得学习算法不仅能够拟合数据，而且能够使模型的参数权重尽量的小。注意到这个正则项只有在训练过程中才会被加到损失函数。当得到完成训练的模型后，我们应该使用没有正则化的测量方法去评价模型的表现。

岭回归损失函数

$J(\theta)=MSE(\theta)+\alpha\frac{1}{2}\sum\limits_{i=1}^n\theta_i^2$

值得注意的是偏差 $\theta_0$ 是没有被正则化的（累加运算的开始是 i=1 而不是 i=0）。如定义 $\mathbf{w}$ 作为特征的权重向量（$\theta_1$到 $\theta_n$），那么正则项可以简写成 $\frac{1}{2}{({\parallel \mathbf{w}\parallel_2})}^2$，其中 $\parallel \cdot \parallel_2$ 表示权重向量的 $\ell_2$ 范数。对于梯度下降来说仅仅在均方差梯度向量加上一项 $\alpha\mathbf{w}$。

在使用岭回归前，对数据进行放缩（可以使用StandardScaler）是非常重要的，算法对于输入特征的数值尺度（scale）非常敏感。大多数的正则化模型都是这样的。

岭回归的封闭方程的解

$\hat{\theta} = ({\mathbf{X}}^T\cdot\mathbf{X}+\alpha\mathbf{A})^{-1}\cdot{\mathbf{X}}^T\cdot\mathbf{y}$

使用 Scikit-Learn 来进行封闭方程的求解:

>>> from sklearn.linear_model import Ridge
>>> ridge_reg = Ridge(alpha=1, solver="cholesky")
>>> ridge_reg.fit(X, y)
>>> ridge_reg.predict([[1.5]])
array([[ 1.55071465]]

使用随机梯度法进行求解：

>>> sgd_reg = SGDRegressor(penalty="l2")
>>> sgd_reg.fit(X, y.ravel())
>>> sgd_reg.predict([[1.5]])
array([[ 1.13500145]])

penalty参数指的是正则项的惩罚类型。指定“l2”表明你要在损失函数上添加一项：权重向量 $\ell_2$ 范数平方的一半，这就是简单的岭回归。

---

Lasso 回归

Lasso 回归（也称 Least Absolute Shrinkage，或者 Selection Operator Regression）是另一种正则化版的线性回归：就像岭回归那样，它也在损失函数上添加了一个正则化项，但是它使用权重向量的 $\ell_1$ 范数而不是权重向量 $\ell_2$ 范数平方的一半。

Lasso 回归的损失函数

$J(\theta)=MSE(\theta)+\alpha\sum\limits_{i=1}^n\left|\theta_i \right|$

Lasso 回归的一个重要特征是它倾向于完全消除最不重要的特征的权重（即将它们设置为零）。换句话说，Lasso回归自动的进行特征选择同时输出一个稀疏模型（即，具有很少的非零权重）。

你可以从下图中知道为什么会出现这种情况：在左上角图中，后背景的等高线（椭圆）表示了没有正则化的均方差损失函数（$\alpha=0$），白色的小圆圈表示在当前损失函数上批量梯度下降的路径。前背景的等高线（菱形）表示$\ell_1$惩罚，黄色的三角形表示了仅在这个惩罚下批量梯度下降的路径（$\alpha\rightarrow\infty$）。注意路径第一次是如何到达$\theta_1=0$，然后向下滚动直到它到达$\theta_2=0$。在右上角图中，等高线表示的是相同损失函数再加上一个 $\alpha=0.5$ 的 $\ell_1$ 惩罚。这幅图中，它的全局最小值在 $\theta_2=0$ 这根轴上。批量梯度下降首先到达 $\theta_2=0$，然后向下滚动直到达到全局最小值。 两个底部图显示了相同的情况，只是使用了 $\ell_2$ 惩罚。 规则化的最小值比非规范化的最小值更接近于 $\theta=0$，但权重不能完全消除。

下面是一个使用 Scikit-Learn 的Lasso类的小例子。也可以使用SGDRegressor(penalty=”l1”)来代替它。

>>> from sklearn.linear_model import Lasso
>>> lasso_reg = Lasso(alpha=0.1)
>>> lasso_reg.fit(X, y)
>>> lasso_reg.predict([[1.5]])
array([ 1.53788174]

---

弹性网络（ElasticNet）

弹性网络介于 Ridge 回归和 Lasso 回归之间。它的正则项是 Ridge 回归和 Lasso 回归正则项的简单混合，同时你可以控制它们的混合率 r，当 r=0 时，弹性网络就是 Ridge 回归，当 r=1 时，其就是 Lasso 回归。

弹性网络损失函数:

$J(\theta)=MSE(\theta)+r\alpha\sum\limits_{i=1}^n\left|\theta_i \right|+\frac{1-r}{2}\alpha\sum\limits_{i=1}^n\theta_i^2$

那么我们该如何选择线性回归，岭回归，Lasso 回归，弹性网络呢？一般来说有一点正则项的表现更好，因此通常你应该避免使用简单的线性回归。岭回归是一个很好的首选项，但是如果你的特征仅有少数是真正有用的，你应该选择 Lasso 和弹性网络。就像我们讨论的那样，它两能够将无用特征的权重降为零。一般来说，弹性网络的表现要比 Lasso 好，因为当特征数量比样本的数量大的时候，或者特征之间有很强的相关性时，Lasso 可能会表现的不规律。下面是一个使用 Scikit-Learn ElasticNet的简单样本：

>>> from sklearn.linear_model import ElasticNet
>>> elastic_net = ElasticNet(alpha=0.1, l1_ratio=0.5)
>>> elastic_net.fit(X, y)
>>> elastic_net.predict([[1.5]])
array([ 1.54333232])

---

Early stopping

对于迭代学习算法，有一种非常特殊的正则化方法，就像梯度下降在验证错误达到最小值时立即停止训练那样。我们称为早期停止法。随着训练的进行，算法一直学习，它在训练集上的预测误差（RMSE）自然而然的下降。然而一段时间后，验证误差停止下降，并开始上升。这意味着模型在训练集上开始出现过拟合。一旦验证错误达到最小值，便提早停止训练。这种简单有效的正则化方法被 Geoffrey Hinton 称为“完美的免费午餐”。

from sklearn.preprocessing import StandardScaler
np.random.seed(42)
m = 100
X = 6 * np.random.rand(m, 1) - 3
y = 2 + X + 0.5 * X**2 + np.random.randn(m, 1)

X_train, X_val, y_train, y_val = train_test_split(X[:50], y[:50].ravel(), test_size=0.5, random_state=10)

poly_scaler = Pipeline([
        ("poly_features", PolynomialFeatures(degree=90, include_bias=False)),
        ("std_scaler", StandardScaler()),
    ])

X_train_poly_scaled = poly_scaler.fit_transform(X_train)
X_val_poly_scaled = poly_scaler.transform(X_val)
# 数据预处理，增加高阶属性，并归一化

sgd_reg = SGDRegressor(max_iter=1,
                       penalty=None,
                       eta0=0.0005,
                       warm_start=True,
                       learning_rate="constant",
                       random_state=42)

n_epochs = 500
train_errors, val_errors = [], []
for epoch in range(n_epochs):
    sgd_reg.fit(X_train_poly_scaled, y_train)
    y_train_predict = sgd_reg.predict(X_train_poly_scaled)
    y_val_predict = sgd_reg.predict(X_val_poly_scaled)
    train_errors.append(mean_squared_error(y_train_predict, y_train))
    val_errors.append(mean_squared_error(y_val_predict, y_val))

best_epoch = np.argmin(val_errors)
best_val_rmse = np.sqrt(val_errors[best_epoch])

# 做注释
plt.annotate('Best model',
             xy=(best_epoch, best_val_rmse),
             xytext=(best_epoch, best_val_rmse + 1),
             ha="center",
             arrowprops=dict(facecolor='black', shrink=0.05),
             fontsize=16,
            )

best_val_rmse -= 0.03  # just to make the graph look better
plt.plot([0, n_epochs], [best_val_rmse, best_val_rmse], "k:", linewidth=2)
plt.plot(np.sqrt(val_errors), "b-", linewidth=3, label="Validation set")
plt.plot(np.sqrt(train_errors), "r--", linewidth=2, label="Training set")
plt.legend(loc="upper right", fontsize=14)
plt.xlabel("Epoch", fontsize=14)
plt.ylabel("RMSE", fontsize=14)
plt.show()
best_epoch