《蜥蜴书》_讲义及源码解读_04

线性回归

线性回归预测模型

$$\begin{gathered} \hat{y}={\theta }_0+{\theta }_1{x}_1+{\theta }_2{x}_2+\cdots +{\theta }_n{x}_n \\ \hat{y}:表示预测结果 \\ n：表示特征的个数 \\ {x}_i:表示第i个特征的值,x_0恒为1 \\ {\theta }_n:表示第n个参数（包括偏置项{\theta }_0和特征权重值{\theta }_1\cdots {\theta }_n） \end{gathered}$$

线性回归预测模型（向量形式）

$$\begin{gathered} \hat{y}=h_{\theta }(x)={\theta }^{T}x \\ {\theta }^T:表示向量的转置（行向量变为了列向量） \\ X:为每个样本中特征值的向量形式，包括x_0到x_n，而且x_0恒为1 \\ {\theta }^T\cdot X:表示{\theta }^T和X的点积 \\ {h}_{\theta }:表示参数为的假设函数 \end{gathered}$$

线性回归模型的MSE 损失函数

$$\begin{gathered} MSE(X,h_{\theta })=\frac{1}{m}\sum _{i=1}^m(\hat{y}-y{)}^2=MSE(\theta ) \\ y:是一个向量，样本实例中的标签向量，其包含了y_1到y_m的值 \end{gathered}$$

标准方程(The Normal Equation，正规方程)

为了找到最小化损失函数的 值，可以采用公式解，换句话说，就是可以通过解正规方程直接得到最后的结果。
$$\begin{gathered} \hat{\theta }=(X^T\cdot X{)}^{-1}\cdot X^T\cdot y \\ \hat{\theta }:指最小化损失\theta 的值 \\ y:是一个向量，样本实例中的标签向量，其包含了y_1到y_m的值 \\ 当(X^T\cdot X)是不可逆的（奇异的）， \\ 可以使用如下方法（伪逆，奇异值分解的知识）: \\ \hat{\theta }=X^{+}\cdot y \\ {X}^{+}:X的伪逆 \end{gathered}$$

示例

此例，可以看做求解目标函数，$y=3x_0+4x_1$ 其中 $\hat{\theta }=({\theta }_0,{\theta }_1),最优解应该是：(3,4)$

import numpy as np

X = 2 * np.random.rand(100, 1) #均匀分布  
y = 4 + 3 * X + np.random.randn(100, 1)  
# y =4+3x， x_0恒等于1，加入高斯噪声，

plt.plot(X, y, "b.")
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.axis([0, 2, 0, 15])
plt.show()

使用标准方程（正规方程）闭式求解 $\hat{\theta }$ 示例：

#将(100,1)的X,转换成(100,2) X_b =(X_0,X_1)  X_0全为1，X_1 为训练数据

X_b = np.c_[np.ones((100, 1)), X]  # add x_0 = 1 to each instance
theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)  
#如上正规方程求解 \theta  ,得：
#array([[4.21509616],
#       [2.77011339]])


X_new = np.array([[0], [2]])  #2行1列的，新的X1: 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) 
#用上面正规方程求解的到的参数\theta,来预测\hat y: y1,y2
#
plt.plot(X_new, y_predict, "r-")
plt.plot(X, y, "b.")
plt.axis([0, 2, 0, 15])
plt.show()

计算复杂度：

标准方程（正规方程）计算矩阵 $X^T\cdot X$ 的逆，它是一个$(n+1)(n+1)$ 的矩阵（$n$ 是特征的个数）。这样一个矩阵求逆的运算复杂度大约在$O(n^{2.4})$ 到 $O(n^3)$ 之间，取决于具体实现。即，若特征个数 n 翻倍，其计算时间大概会是原来的5-8倍。

使用LinearRegression模块的SVD奇异值分解的技术来计算，复杂度大约在$O(n^2)$ ，稍优

提示：当特征的个数较大的时候（例如：特征数量为 100000），以上方法的求解计算将会非常慢。

重点：内存不够呢？

使用Sklearn的线性回归模块LinearRegression闭式求解示例：

from sklearn.linear_model import LinearRegression

lin_reg = LinearRegression()
lin_reg.fit(X, y)
lin_reg.intercept_, lin_reg.coef_   #偏置和权重，即\theta_0 和 \theta_1

#(array([4.21509616]), array([[2.77011339]]))

#还是使用Numpy的线性代数库的最小二乘函数：np.linalg.lstsq
theta_best_svd, residuals, rank, s = np.linalg.lstsq(X_b, y, rcond=1e-6)
theta_best_svd
#或如下直接求伪逆
np.linalg.pinv(X_b).dot(y)

梯度下降

梯度下降是一种非常通用的优化算法，它能够很好地解决一系列最优解问题。梯度下降的整体思路是通过的迭代来逐渐调整参数使得损失函数达到最小值。

对损失函数求偏导

$$\begin{gathered} MSE(X,h_{\theta })=\frac{1}{m}\sum _{i=1}^m(\hat{y}-y{)}^2=MSE(\theta ) \\ \frac{\partial }{\partial {\theta }_j}MSE(\theta )=\frac{2}{m}\sum _{i=1}^m({\theta }^T{x}^{(i)}-y^{(i)})x_j^{(i)} \\ {x}^{(i)},y^{(i)}:第i个样本实例的的特征数据和标签标签,共m个样本实例 \\ j：每一个实例样本的特征数量，0-n \\ 换一种写法(一次性计算梯度向量，即计算出关于所有\theta 的偏导)： \\ {\nabla }_{\theta }MSE(\theta )=\left(\begin{array}{c}\frac{\partial }{\partial {\theta }_0}MSE(\theta )\\ \frac{\partial }{\partial {\theta }_1}MSE(\theta )\\ ⋮\\ \frac{\partial }{\partial {\theta }_n}MSE(\theta )\end{array}\right)=\frac{2}{m}{X}^T(X\theta -y) \end{gathered}$$

梯度下降迭代公式

$${\theta }_{next}={\theta }_{current}-\eta {\nabla }_{\theta }MSE(\theta )$$ { }

批量梯度下降

示例代码：

eta = 0.1  # learning rate  \eta
n_iterations = 1000  # 手工指定迭代次数
m = 100              # 样本实例个数

theta = np.random.randn(2,1)  # 随机初始化\theta_0 和 \theta_1

for iteration in range(n_iterations):
    gradients = 2/m * X_b.T.dot(X_b.dot(theta) - y)  
    	#用MSE关于\theta 的求导公式求导，得到梯度
    theta = theta - eta * gradients  
    	#梯度下降迭代
    
print(theta)

#array([[4.21509616],
#       [2.77011339]])

讨论下梯度下降中的学习率\eta: $\eta$

theta_path_bgd = []  #记录\theta 变化情况

def plot_gradient_descent(theta, eta, theta_path=None):
    m = len(X_b)
    plt.plot(X, y, "b.")   #描原始的100个点
    n_iterations = 1000   #迭代1000次
    for iteration in range(n_iterations):
        if iteration < 10:                  #只绘制前10次迭代
            y_predict = X_new_b.dot(theta)    #用当前的\theta 向量，计算预测的y
            style = "b-" if iteration > 0 else "r--"  #等于0时即为初始时，描红色虚线，否则蓝线
            plt.plot(X_new, y_predict, style)     #描线
        gradients = 2/m * X_b.T.dot(X_b.dot(theta) - y)  #计算梯度
        theta = theta - eta * gradients   #更新\theta
        if theta_path is not None:   
            theta_path.append(theta)
    plt.xlabel("$x_1$", fontsize=18)
    plt.axis([0, 2, 0, 15])
    plt.title(r"$\eta = {}$".format(eta), fontsize=16)
    
np.random.seed(42)
theta = np.random.randn(2,1)  #随机初始化\theta

plt.figure(figsize=(10,4))
plt.subplot(131); plot_gradient_descent(theta, eta=0.02)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.subplot(132); plot_gradient_descent(theta, eta=0.1, theta_path=theta_path_bgd)
plt.subplot(133); plot_gradient_descent(theta, eta=0.5)

plt.show()

theta_path_bgd  
#得到 学习率 eta=0.1时，参数 \theta的变化情况。

学习率、迭代次数和收敛速率：

在左图中，学习率是最小的，算法几乎不能求出最后的结果，而且还会花费大量的时间。

在中间图中，学习率的表现看起来不错，仅仅几次迭代后，它就收敛到了最后的结果。

在右图中，学习率太大了，算法是发散的，跳过了所有的训练样本，同时每一步都离正确的结果越来越远。

为了找到一个好的学习率，你可以使用网格搜索。当然，你一般会限制迭代的次数，以便网格搜索在可接受时间内收敛。

关于迭代的次数：

• 如果它太小了，当算法停止的时候，你依然没有找到最优解。
• 如果它太大了，算法会非常的耗时同时后来的迭代参数也不会发生改变。

解决方法是：设置一个非常大的迭代次数，但是当梯度向量变得非常小的时候，结束迭代。早停技术early stopping

这里非常小的梯度向量指：小于一个值 \epsilon $ϵ$（称为容差）。

关于收敛速率：

当损失函数是凸函数，同时它的斜率不能突变（就像均方差损失函数那样），那么它的批量梯度下降算法固定学习率之后，它的收敛速率是$O(\frac{1}{ϵ})$ 。即，若将容差 缩小 10 倍，这个算法的迭代次数大约会变成原来的 10 倍。

随机梯度下降（SGD)

Stochastic Gradient Descent

批量梯度下降：要使用整个训练集来计算每一步的梯度，当训练集很大时，算法会非常慢。

随机梯度下降: 在每一步的梯度计算上只随机选取训练集中的一个样本实例。并且仅基于该单个样本实例来计算梯度，算法变得非常快。由于每一次迭代，只需要在内存中有一个实例，这使随机梯度算法可以在大规模训练集上使用。

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

由于它的随机性，与批量梯度下降相比，其呈现出更多的不规律性：它到达最小值不是平缓的下降，损失函数会忽高忽低，只是在大体上呈下降趋势。随着时间的推移，它会非常的靠近最小值，但是它不会停止在一个值上，它会一直在这个值附近摆动。因此，当算法停止的时候，最后的参数还不错，但不是最优值。

随机性可好处在于可以逃离局部最优，但缺点是永远定位不出最小值。要解决这个问题。有一个好办法是：逐步降低学习率

代码示例：（因为是随机选取样本实例，所以，随机梯度下降前，要对训练集进行混洗，不能按标签排序）

theta_path_sgd = []
m = len(X_b)  # 添加了噪音的X(100,2) 训练数据的长度，即100
np.random.seed(42)

n_epochs = 50
t0, t1 = 5, 50  # learning schedule hyperparameters

def learning_schedule(t):
    return t0 / (t + t1)  #  控制学习率随t增加而下降，下降幅度不大于1/50?

theta = np.random.randn(2,1)  # 随机初始化\theta

for epoch in range(n_epochs):
    for i in range(m):
        if epoch == 0 and i < 20:                    #
            y_predict = X_new_b.dot(theta)           #
            style = "b-" if i > 0 else "r--"         #
            plt.plot(X_new, y_predict, style)        # 第一轮，用初始化的起始参数绘图，红色虚线
        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)       # epoch在（0,50）遍历，i在（1,100）遍历，100样本实例
        theta = theta - eta * gradients				 # 梯度下降	eta在变化，越来越小
        theta_path_sgd.append(theta)                 # 记录梯度变化

plt.plot(X, y, "b.")                                
plt.xlabel("$x_1$", fontsize=18)                    
plt.ylabel("$y$", rotation=0, fontsize=18)          
plt.axis([0, 2, 0, 15])                             
save_fig("sgd_plot")                               
plt.show()   
theta
#array([[4.21566571],
#       [2.85414642]])

使用SKlearn的SGDRegressor模块实现，示例代码

from sklearn.linear_model import SGDRegressor
#迭代1000次，eps为0.001，eta初始值为0.1 ，不使用任何正则化penalty
sgd_reg = SGDRegressor(max_iter=1000, tol=1e-3, penalty=None, eta0=0.1, random_state=42)
sgd_reg.fit(X, y.ravel())    #把列向量降维成一维数组，拉平。

sgd_reg.intercept_, sgd_reg.coef_

#(array([4.24365286]), array([2.8250878]))

小批量梯度下降

• 小型批量的随机实例集上计算梯度。

• 小批量梯度下降相比随机梯度下降的主要优点是：可以通过矩阵操作的硬件优化来提高性能，特别在使用GPU时。

代码示例：

theta_path_mgd = []

n_iterations = 50
minibatch_size = 20

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

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]   # 混洗之后的X，y
    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]    # 小批量 20个样本实例
        gradients = 2/minibatch_size * xi.T.dot(xi.dot(theta) - yi)  # 计算小批量梯度
        eta = learning_schedule(t)       # 学习率变换
        theta = theta - eta * gradients  # 梯度下降
        theta_path_mgd.append(theta)     # 记录梯度

三种梯度下降的总结

依前面示例代码记录的梯度下降过程的变量theta_path_bgd ,theta_path_sgd 和theta_path_mgd, 绘图总结

theta_path_bgd = np.array(theta_path_bgd)
theta_path_sgd = np.array(theta_path_sgd)
theta_path_mgd = np.array(theta_path_mgd)

plt.figure(figsize=(7,4))
plt.plot(theta_path_sgd[:, 0], theta_path_sgd[:, 1], "r-s", linewidth=1, label="Stochastic")
plt.plot(theta_path_mgd[:, 0], theta_path_mgd[:, 1], "g-+", linewidth=2, label="Mini-batch")
plt.plot(theta_path_bgd[:, 0], theta_path_bgd[:, 1], "b-o", linewidth=3, label="Batch")
plt.legend(loc="upper left", fontsize=16)
plt.xlabel(r"$\theta_0$", fontsize=20)
plt.ylabel(r"$\theta_1$   ", fontsize=20, rotation=0)
plt.axis([2.5, 4.5, 2.3, 3.9])

plt.show()

多项式回归Polynomial regression

大多数实际项目中的数据关系不是简单的线性关系，可以使用线性模型来拟合非线性数据。

方法就是对每个特征进行加权后添加为一个新特征，然后在这个扩展特征上训练一个线性模型， 这种方法称为多项式回归。

示例：

import numpy as np
import numpy.random as rnd

np.random.seed(42)
m = 100
X = 6 * np.random.rand(m, 1) - 3   # 放大6倍 在-3，数据分布在-3到+3之间
y = 0.5 * X**2 + X + 2 + np.random.randn(m, 1)  # y = 0.5x^2 + x + 2 + 高斯噪声  

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

使用SKlearn 中的PolynomialFeatures 将训练数据中每个特征的平方（二次多项式）添加为新特征。（此例只有一个特征：即 n=1）

from sklearn.preprocessing import PolynomialFeatures
poly_features = PolynomialFeatures(degree=2, include_bias=False) 
#degree 2阶，此例中只有1个特征x，对应产生1个新特征x^2,
#若多个特征(a,b),2阶degree,则产生 a^2 ,b^2 和ab 这个些新特征；3阶degree，则产生，a^2 ,a^3, b^2, b^3, ab,ab^2 ,a^2 b

X_poly = poly_features.fit_transform(X)  #得到预处理后的新训练集

X[0],X_poly[0]
#如下，产生新特征 0.5666....
#array([-0.75275929]), array([-0.75275929,  0.56664654])

lin_reg = LinearRegression()
lin_reg.fit(X_poly, y)            # 对预处理后的训练集进行线性回归，训练样本m=100, 特征个数n=2
lin_reg.intercept_, lin_reg.coef_  
#(array([1.78134581]), array([[0.93366893, 0.56456263]]))
#效果还行，偏置 \theta_0  1.78==2,  权重系数：\theta_1 == 1 ,\theta_2  ==0.5

X_new=np.linspace(-3, 3, 100).reshape(100, 1)# 随机产生（-3,3）100个测试用例
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()

学习曲线

先看示例：

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

for style, width, degree in (("g-", 1, 300), ("b--", 2, 2), ("r-+", 2, 1)):
    polybig_features = PolynomialFeatures(degree=degree, include_bias=False)  
    # degree=300,2,1  预处理数据 这里的绿色，细线，300阶有点扩张，尝试下5阶，6阶
    std_scaler = StandardScaler()
    lin_reg = LinearRegression()
    polynomial_regression = Pipeline([
            ("poly_features", polybig_features),
            ("std_scaler", std_scaler), #数据标准化
            ("lin_reg", lin_reg),
        ])
    # 定义流水线，对多项式特征和标准化后的数据进行线性拟合训练
    polynomial_regression.fit(X, y)   
    y_newbig = polynomial_regression.predict(X_new)  #用测试数据预测，得到3组不同的 y向量，绘图
    plt.plot(X_new, y_newbig, style, label=str(degree), linewidth=width)

plt.plot(X, y, "b.", linewidth=3)
plt.legend(loc="upper left")
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.axis([-3, 3, 0, 10])
plt.show()

结论：高阶多项式回归模型在训练集上严重过拟合了，而线性模型则欠拟合。

当然此例中，二次模型有着较好的泛化能力。因为此数据本来就是二次函数加高斯噪声模拟出来的。

在不知道你要预测的函数的阶是多少时，这么办？

有一种办法，就是学习曲线：这个曲线绘制的是模型在训练集和验证集上关于训练集大小（或训练迭代）的性能函数，也就是：画出模型在训练集上的表现，同时画出以训练集规模为自变量的训练集函数。为了得到图像，需要在训练集的不同规模子集上进行多次训练。

使用普通线性模型示例：

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)
    # 8成训练集，2成验证集 
    train_errors, val_errors = [], []
    for m in range(1, len(X_train)):
        model.fit(X_train[:m], y_train[:m])   # m 从1 - 训练集长度。 拟合
        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)   # 
    plt.xlabel("Training set size", fontsize=14) # 
    plt.ylabel("RMSE", fontsize=14)              # 
    
    
lin_reg = LinearRegression()  #普通线性回归
plot_learning_curves(lin_reg, X, y)
plt.axis([0, 80, 0, 3])                        
plt.show()       

由图结论：欠拟合，RMSE很高，并且随着训练数据增加，RMSE不会降低（对于验证集和训练集）

使用10阶多项式+线性模型示例

from sklearn.pipeline import Pipeline

polynomial_regression = Pipeline([
        ("poly_features", PolynomialFeatures(degree=10, include_bias=False)),
        ("lin_reg", LinearRegression()),
    ])
#10阶多项式特征后的线性回归
plot_learning_curves(polynomial_regression, X, y)  
plt.axis([0, 80, 0, 3])           # 
save_fig("learning_curves_plot")  # 
plt.show()                        # wn

由图结论：误差MSE降低很多。曲线之间的间隙较大（训练数据的性能 更优于验证数据），过拟合的现象，但随着训练数据量增大，间隙慢慢变小。

正则化模型（Regularized models）

降低模型的过拟合的好方法是正则化这个模型（即限制它）：模型有越少的自由度，就越难以拟合数据。

例如：

• 正则化一个多项式模型，简单方法就是减少多项式的阶数。
• 正则化一个线性模型，典型实现就是约束模型中参数的权重。

三种不同约束权重的方法：Ridge 回归，Lasso 回归和 Elastic Net。

岭回归（Ridge）

岭回归（也称为 Tikhonov 正则化）是线性回归的正则化版：在损失函数上直接加上一个正则项

复习下普通线性模型的损失函数：
$$\begin{gathered} MSE(X,h_{\theta })=\frac{1}{m}\sum _{i=1}^m(\hat{y}-y{)}^2=MSE(\theta ) \\ y:是一个向量，样本实例中的标签向量，其包含了y_1到y_m的值 \end{gathered}$$
岭回归模型的损失函数：
$$\begin{gathered} J(\theta )=MSE(\theta )+\alpha \frac{1}{2}\sum _{j=1}^n{\theta }_j^2 \\ 超参数\alpha 为0时，还原成普通线性回归模型。 \\ 超参数\alpha 很大时，所有权重系数\theta 越接近0，变成经过数据均值的一条水平线。 \\ n:为特征数量，注意不包括{\theta }_0,也就是说，偏置项没有被正则化。 \end{gathered}$$
注意，在岭回归之前，要对数据进行规范化、标准化，也就是缩放数据

岭回归的闭式求解公式：对比线性回归的标准方程（正规方程）
$$\begin{gathered} \hat{\theta }=(X^{T}X+\alpha A{)}^{-1}{X}^{T}y \\ A:单位矩阵，主对角线为1，其余为0 \end{gathered}$$

示例代码：使用Sklearn的模块Ridge闭式求解

np.random.seed(42)
m = 20
X = 3 * np.random.rand(m, 1)                    # 20 的演示数据，范围（0-3）
y = 1 + 0.5 * X + np.random.randn(m, 1) / 1.5   # y= 1+0.5x  高斯噪声，训练数据
X_new = np.linspace(0, 3, 100).reshape(100, 1)  # 测试数据

from sklearn.linear_model import Ridge
ridge_reg = Ridge(alpha=1, solver="cholesky", random_state=42)  # cholesky 是岭回归闭式解公式的变体
ridge_reg.fit(X, y)  #训练
ridge_reg.intercept_, ridge_reg.coef_  #偏置，权重系数

#(array([1.00645006]), array([[0.3628467]]))

ridge_reg.predict([[1.5]]) # 预测示例

示例代码：对比线性模型和多项式模型，各种 $\alpha$ 的岭回归正则化。

from sklearn.linear_model import Ridge

def plot_model(model_class, polynomial, alphas, **model_kargs):
    for alpha, style in zip(alphas, ("b-", "g--", "r:")):
        model = model_class(alpha, **model_kargs) if alpha > 0 else LinearRegression()  
        # \alpha=0 就是普通线性模型
        # 如果是多项式模型，则10阶，缩放，再岭回归
        if polynomial:    
            model = Pipeline([
                    ("poly_features", PolynomialFeatures(degree=10, include_bias=False)),
                    ("std_scaler", StandardScaler()),
                    ("regul_reg", model),
                ])
         # 如果不是多项式模型，直接岭回归
        model.fit(X, y)  # 训练拟合     
        y_new_regul = model.predict(X_new)  # 用测试集预测，得到描图数据
        lw = 2 if alpha > 0 else 1   # \alpha=0 对应蓝色，细线，否则是绿色、红色粗线，线宽为2
        plt.plot(X_new, y_new_regul, style, linewidth=lw, label=r"$\alpha = {}$".format(alpha))
    plt.plot(X, y, "b.", linewidth=3)  #描训练数据
    plt.legend(loc="upper left", fontsize=15)
    plt.xlabel("$x_1$", fontsize=18)
    plt.axis([0, 3, 0, 4])

plt.figure(figsize=(8,4))
plt.subplot(121)
plot_model(Ridge, polynomial=False, alphas=(0, 10, 100), random_state=42)  # 不使用多项式模型
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.subplot(122)
plot_model(Ridge, polynomial=True, alphas=(0, 10**-5, 1), random_state=42) # 使用多项式模型
plt.show()

示例代码：使用Sklearn的模块SGD随机梯度下降，penalty参数

sgd_reg = SGDRegressor(penalty="l2", max_iter=1000, tol=1e-3, random_state=42)   
#l2 ,损失函数中增加曾正则项，权重向量l2范数的二分之一，即\alpha=1的岭回归
sgd_reg.fit(X, y.ravel())  #
sgd_reg.predict([[1.5]])

Lasso回归

也称 Least Absolute Shrinkage，或者Selection Operator Regression，是另一种正则化版的线性回归：就像岭回归那样，它也在损失函数上添加了一个正则化项，但是它使用权重向量的 $l_1$ 范数而不是权重向量 $l_2$ 范数平方的一半。
$$\begin{gathered} J(\theta )=MSE(\theta )+\alpha \sum _{j=1}^n|{\theta }_j| \\ 超参数\alpha :，为0时，还原成普通线性回归模型， \\ 超参数\alpha :很大时，所有权重系数\theta 越接近0，变成过数据均值的一条水平线。 \\ n:为特征数量，注意不包括{\theta }_0,也就是说，偏置项没有被正则化。 \end{gathered}$$

示例代码

from sklearn.linear_model import Lasso

plt.figure(figsize=(8,4))
plt.subplot(121)
plot_model(Lasso, polynomial=False, alphas=(0, 0.1, 1), random_state=42)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.subplot(122)
plot_model(Lasso, polynomial=True, alphas=(0, 10**-7, 1), random_state=42)

plt.show()

Lasso 回归的一个重要特征是它倾向于完全消除最不重要的特征的权重（即将它们设置为零）。

如右图中的绿色虚线所示（$\alpha =10^{-7}$时），曲线看起来像一条二次曲线，而且几乎是线性的。这是因为所有的高阶多项特征都被设置为零。即：Lasso 回归自动的进行特征选择，同时输出一个稀疏模型（即，具有很少的非零权重）。

弹性网络（ElasticNet）

弹性网络介于 Ridge 回归和 Lasso 回归之间。它的正则项是 Ridge 回归和 Lasso 回归正则项的简单混合，同时你可以控制它们的混合率 ，当$r=0$ 时，弹性网络就是 Ridge 回归，当 $r=1$时，其就是 Lasso 回归，公式如下：
$$J(\theta )=MSE(\theta )+\alpha r\sum _{j=1}^n|{\theta }_j|+\alpha \frac{1-r}{2}\sum _{j=1}^n{\theta }_j^2$$
示例：

from sklearn.linear_model import ElasticNet
elastic_net = ElasticNet(alpha=0.1, l1_ratio=0.5, random_state=42)  # l1_ratio  即公式中的混合比 r
elastic_net.fit(X, y)
elastic_net.predict([[1.5]])

如何选择正则化模型

如何选择岭回归，Lasso 回归，弹性网络呢？

• 一般来说有一点正则项表现更好，因此通常你应该避免使用简单的纯线性回归。岭回归是一个很好的首选项。

• 假如训练数据实际相关的特征仅有少数几个，那就更倾向于选择 Lasso 和弹性网络。因为，这两个正则化模型能够将无用特征的权重降为零。

• 一般来说，弹性网络的表现要比 Lasso 好，因为当特征数量比样本的数量大的时候，或者特征之间有很强的相关性时，Lasso 可能会表现的不规律。

提前停止（Early Stopping）

随着训练的进行，算法一直学习，它在训练集上的预测误差（RMSE）自然而然的下降。然而一段时间后，验证误差停止下降反而开始上升。这意味着模型在训练集上开始出现过拟合。

通过早期停止算法一旦验证错误达到最小值，便提早停止训练。这种简单有效的正则化方法被Geoffrey Hinton 称为“完美的免费午餐”。

示例代码

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)
#分割训练集和验证集，各一半。

from copy import deepcopy

poly_scaler = Pipeline([
        ("poly_features", PolynomialFeatures(degree=90, include_bias=False)), #  90阶多项式模型
        ("std_scaler", StandardScaler())
    ])

X_train_poly_scaled = poly_scaler.fit_transform(X_train)
X_val_poly_scaled = poly_scaler.transform(X_val)  # 标准化缩放数据

#epoch 1000次，先拿到随机梯度下降模型中最佳的迭代次数epoch,即best_epoch
sgd_reg = SGDRegressor(max_iter=1, tol=-np.infty, warm_start=True,  #为了绘图演示，设置了早停止之后继续
                       penalty=None, learning_rate="constant", eta0=0.0005, random_state=42)

minimum_val_error = float("inf")
best_epoch = None
best_model = None
for epoch in range(1000):
    sgd_reg.fit(X_train_poly_scaled, y_train)  			  # 训练集数据拟合，随循环在停止迭代的地方继续迭代
    y_val_predict = sgd_reg.predict(X_val_poly_scaled)    # 用验证集数据预测
    val_error = mean_squared_error(y_val, y_val_predict)  # 计算验证集的MSE
    if val_error < minimum_val_error:
        minimum_val_error = val_error
        best_epoch = epoch
        best_model = deepcopy(sgd_reg) 
        # 验证集测试MSE不降反升的话，提前深拷贝模型，为了保留好最佳
        
#epoch500次，重新跑一次训练，迭代次数500     
sgd_reg = SGDRegressor(max_iter=1, tol=-np.infty, warm_start=True,
                       penalty=None, learning_rate="constant", eta0=0.0005, 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, y_train_predict))  # 记录训练集MSE
    val_errors.append(mean_squared_error(y_val, y_val_predict))        # 记录验证集MSE

best_epoch = np.argmin(val_errors)
best_val_rmse = np.sqrt(val_errors[best_epoch])
#在途中标记注释最佳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  # 
plt.plot([0, n_epochs], [best_val_rmse, best_val_rmse], "k:", linewidth=2)  
# 绘制最小MSE水平线,下移了0.03
plt.plot(np.sqrt(val_errors), "b-", linewidth=3, label="Validation set")    
# 描绘验证机MSE变化曲线
plt.plot(np.sqrt(train_errors), "r--", linewidth=2, label="Training set")
# 描绘训练集MSE变化曲线
plt.legend(loc="upper right", fontsize=14)
plt.xlabel("Epoch", fontsize=14)
plt.ylabel("RMSE", fontsize=14)
plt.show()
best_epoch, best_model
#(239,
# SGDRegressor(eta0=0.0005, learning_rate='constant', max_iter=1, penalty=None,
#              random_state=42, tol=-inf, warm_start=True))

逻辑回归Logistic

Logistic 回归（也称为 Logit 回归）通常用于估计一个实例属于某个特定类别的概率（例如，这电子邮件是垃圾邮件的概率是多少？）。 如果估计的概率大于 50%，那么模型预测这个实例属于当前类（称为正类，标记为“1”），反之预测它不属于当前类（即它属于负类 ，标记为“0”）。 这样便成为了一个二元分类器。

线性回归预测模型公式（向量形式）：
$$\begin{gathered} \hat{y}=h_{\theta }(x)={\theta }^{T}x \\ \theta :表示模型的参数向量包括偏置项和特征权重值{\theta }_1到{\theta }_2 \\ {\theta }^T:表示向量的转置（行向量变为了列向量） \end{gathered}$$

就像线性回归模型一样，Logistic 回归模型计算输入特征的加权和（加上偏差项），但它不像线性回归模型那样直接输出结果，而是把结果输入logistic()函数进行二次加工后进行输出(数理逻辑值)。

逻辑回归预测模型公式

$$\begin{gathered} \hat{P}=h_{\theta }(X)=\sigma (X^T\theta ) \\ \sigma ()函数：大名鼎鼎的S型函数： \\ \sigma (t)=\frac{1}{1+e^{-t}} \end{gathered}$$

逻辑回归的损失函数公式

$$\begin{gathered} J(\theta )=-\frac{1}{m}\sum _{i=1}^m[y^{(i)}log({\hat{p}}^{(i)})+(1-y^{(i)})log(1-{\hat{p}}^{(i)})] \\ m:样本实例数量 \end{gathered}$$

逻辑回归损失函数的偏导数方程：

$$\begin{gathered} \frac{\partial }{\partial {\theta }_j}J(\theta )=\frac{1}{m}\sum _{i=1}^m(\sigma ({\theta }^T{x}^{(i)})-y^{(i)})x_j^{(i)} \\ i:m个样本实例数量中第i个 \\ j:n个特征数量，第j个特征 \\ \theta :参数向量({\theta }_0,thet{a}_1,\cdots ,thet{a}_n),{\theta }_0是偏置 \end{gathered}$$

决策边界

用著名的鸢尾花数据集来示例：

from sklearn import datasets
iris = datasets.load_iris()  # 鸢尾花数据集

print(iris.data.shape)
print(iris.data[:10])
print(iris.feature_names)
print(iris.target.shape)
print(iris.target[:10])
list(iris.keys())  #观察数据

X = iris["data"][:, 3:]  # 简化实验，只保留第三列特征，即petal width
y = (iris["target"] == 2).astype(np.int)  # 1 if Iris virginica, else 0，  如果是第2号类别，维吉尼亚鸢尾花，转成1，否则为0

from sklearn.linear_model import LogisticRegression
log_reg = LogisticRegression(solver="lbfgs", random_state=42)
log_reg.fit(X, y)

X_new = np.linspace(0, 3, 1000).reshape(-1, 1)
y_proba = log_reg.predict_proba(X_new)
decision_boundary = X_new[y_proba[:, 1] >= 0.5][0]

plt.figure(figsize=(8, 3))
plt.plot(X[y==0], y[y==0], "bs")
plt.plot(X[y==1], y[y==1], "g^")
plt.plot([decision_boundary, decision_boundary], [-1, 2], "k:", linewidth=2)
plt.plot(X_new, y_proba[:, 1], "g-", linewidth=2, label="Iris virginica")
plt.plot(X_new, y_proba[:, 0], "b--", linewidth=2, label="Not Iris virginica")
plt.text(decision_boundary+0.02, 0.15, "Decision  boundary", fontsize=14, color="k", ha="center")
plt.arrow(decision_boundary, 0.08, -0.3, 0, head_width=0.05, head_length=0.1, fc='b', ec='b')
plt.arrow(decision_boundary, 0.92, 0.3, 0, head_width=0.05, head_length=0.1, fc='g', ec='g')
plt.xlabel("Petal width (cm)", fontsize=14)
plt.ylabel("Probability", fontsize=14)
plt.legend(loc="center left", fontsize=14)
plt.axis([0, 3, -0.02, 1.02])
save_fig("logistic_regression_plot")
plt.show()

decision_boundary
#array([1.66066066])

SOFTMAX回归

Logistic 回归模型可以直接推广到支持多类别分类，不必组合和训练多个二分类器，其称为Softmax 回归或多类别Logistic 回归.

某类 $K$ 的Softmax 得分：
$$\begin{gathered} S_k(X)=X^T{\theta }^{(k)} \\ 每个类都有自己特定的参数向量{\theta }^{(k)},所有这些向量作为行存储在参数矩阵\Theta 中 \\ k:分类数量 \\ X:样本实例 \\ m:样本实例的数量 \\ n:样本特征的数量 \\ \Theta :k行n列的参数矩阵 \end{gathered}$$
一旦通过样本实例$X$ 计算了每个类的得分，就可以通过Softmax函数来估计样本实例属于第 $K$ 类的概率 ：

Softmax 函数公式：
$$\begin{gathered} {\hat{p}}_k=\sigma (S_k(X))=\frac{e^{S_k(X)}}{{\sum }_{i=1}^k{e}^{S_k(X)}} \\ 对某个实例样本，计算第k类的得分的指数，然后对它们进行归一化（除以所有类上得分的指数的总和） \end{gathered}$$
和 Logistic 回归分类器一样，Softmax 回归分类器将估计概率最高（它只是得分最高的类）的那类作为预测结果。

Softmax 回归分类预测：
$$\hat{y}=argmax\sigma (S_k(X))$$
Softmax 回归分类器一次只能预测一个类（即它是多类的，但不是多输出的），因此它只能用于判断互斥的类别，如不同类型的植物。 你不能用它来识别一张照片中的多个人。

Softmax 回归的损失函数：交叉熵
$$\begin{gathered} J(\Theta )=-\frac{1}{m}\sum _{i=1}^m\sum _{k=1}^K{y}_k^{(i)}log({\hat{p}}_k^{(i)}) \\ {y}_k^{(i)}:第i个样本实例属于类k的概率，一般是1或者0 \\ 当（k=2）时，等价于逻辑回归的损失函数 \end{gathered}$$
Softmax 回归的损失函数相对于${\theta }_k$的梯度
$$\begin{gathered} {\nabla }_{\theta }(k)J(\Theta )=-\frac{1}{m}\sum _{i=1}^m({\hat{p}}_k^{(i)}-y_k^{(i)})x^{(i)} \\ 可以计算每一类的梯度向量，然后使用梯度下降（或者其他的优化算法） \\ 找到使得损失函数达到最小值的参数矩阵\Theta \end{gathered}$$
示例：

X = iris["data"][:, (2, 3)]  # petal length, petal width
y = iris["target"]

softmax_reg = LogisticRegression(multi_class="multinomial",solver="lbfgs", C=10, random_state=42)
softmax_reg.fit(X, y)
#multi_class ：分类方式选择参数，可选参数为ovr和multinomial，默认为ovr。
#ovr即前面提到的one-vs-rest(OvR)，
#而multinomial即前面提到的many-vs-many(MvM)
#将multi_class切换成multinomial后，即转换为Softmax回归！！！

#solver 优化算法选择参数，只有五个可选参数，即newton-cg,lbfgs,liblinear,sag,saga。默认为liblinear。
#lbfgs：牛顿法的一种，即最小二乘法，利用损失函数二阶导数矩阵即海森矩阵来迭代优化损失函数。即，l2范数。

x0, x1 = np.meshgrid(
        np.linspace(0, 8, 500).reshape(-1, 1),
        np.linspace(0, 3.5, 200).reshape(-1, 1),
    )
X_new = np.c_[x0.ravel(), x1.ravel()]


y_proba = softmax_reg.predict_proba(X_new)
y_predict = softmax_reg.predict(X_new)

zz1 = y_proba[:, 1].reshape(x0.shape)
zz = y_predict.reshape(x0.shape)

plt.figure(figsize=(10, 4))
plt.plot(X[y==2, 0], X[y==2, 1], "g^", label="Iris virginica")
plt.plot(X[y==1, 0], X[y==1, 1], "bs", label="Iris versicolor")
plt.plot(X[y==0, 0], X[y==0, 1], "yo", label="Iris setosa")

from matplotlib.colors import ListedColormap
custom_cmap = ListedColormap(['#fafab0','#9898ff','#a0faa0'])

plt.contourf(x0, x1, zz, cmap=custom_cmap)
contour = plt.contour(x0, x1, zz1, cmap=plt.cm.brg)
plt.clabel(contour, inline=1, fontsize=12)
plt.xlabel("Petal length", fontsize=14)
plt.ylabel("Petal width", fontsize=14)
plt.legend(loc="center left", fontsize=14)
plt.axis([0, 7, 0, 3.5])

plt.show()