线性回归在线拟合器 线性回归拟合度分析

所有代码块都是在Jupyter Notebook下进行调试运行，前后之间都相互关联。
文中所有代码块所涉及到的函数里面的详细参数均可通过scikit-learn官网API文档进行查阅，这里我只写下每行代码所实现的功能，参数的调整读者可以多进行试验调试。多动手！！！

主要内容：
线性回归方程实现
梯度下降效果
对比不同梯度下降策略
建模曲线分析
过拟合与欠拟合
正则化的作用
提前停止策略

一、线性回归

Ⅰ、参数直接求解法

前几章节中也都得出求解的方案，基于最小二乘法直接求解，但这是人工进行求解，为此引入了梯度下降法来训练机器学习。

回归方程，也就是最小二乘法求解，确实回归方程就是迭代更新即可。机器学习中核心的思想是迭代更新。

通过这些公式，我们需要把这个θ求解出来。

接下来，我们玩一个有意思的：将随机生成100个点，然后再进行线性拟合，拟合出一条较好的直线。

以下代码块都是基于Jupyter notebook 导包

import numpy as np
import os
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
plt.rcParams['axes.labelsize'] = 14
plt.rcParams['xtick.labelsize'] = 12
plt.rcParams['ytick.labelsize'] = 12
import warnings
warnings.filterwarnings('ignore')
np.random.seed(42)

# 随机生成100个点，随机构造一个方程y = 4 + 3*x + b
import numpy as np
X = 2*np.random.rand(100,1)#随机构造数据，选取100个点构建矩阵
y = 4+ 3*X +np.random.randn(100,1)#随机构造一个方程

# 显示一下这些点
plt.plot(X,y,'g.')#颜色设定和随机选取点
plt.xlabel('X_1')
plt.ylabel('y')
plt.axis([0,2,0,15])#x轴0-2和y轴0-15的取值范围
plt.show()

X_b = np.c_[np.ones((100,1)),X]#为了矩阵运算方便，需要再加一列全为1的向量，此时的X_b就是带有偏置项的矩阵
theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)#通过公式求解θ，.T转置，.dot相乘，inv求逆

theta_best#得到偏置项和权重项
"""
array([[4.21509616],
       [2.77011339]])
"""

X_new = np.array([[0],[2]])#定义测试数据，两个点就行，因为两点就能确定一条直线
X_new_b = np.c_[np.ones((2,1)),X_new]#训练的时候都加了一列1，测试的时候需要加一列1
y_predict = X_new_b.dot(theta_best)#当前得到的测试数据去乘以权重参数，得到最终的结果，相对于x×θ得到最后的y
y_predict#得出预测的结果值
"""
array([[4.21509616],
       [9.75532293]])
"""

X_new#当前的值
"""
array([[0],
       [2]])
"""

#将得到的结果进行绘制
plt.plot(X_new,y_predict,'b--')#由点x通过预测得到的y进行绘制
plt.plot(X,y,'g.')#上面的我们绘制的随机点
plt.axis([0,2,0,15])
plt.show()

最终即可得到通过训练一百个随机数，进行拟合得权重参数θ，从而绘制出线性拟合直线。

工具包实现求解权重和偏置参数
sklearn api文档

from sklearn.linear_model import LinearRegression
lin_reg = LinearRegression()
lin_reg.fit(X,y)
print (lin_reg.coef_)#权重参数
print (lin_reg.intercept_)#偏置参数
#很显然，跟上述的求解的theta_best一致
"""
[[2.77011339]]
[4.21509616]
"""

Ⅱ、预处理对结果的影响

选择一个初始值，沿着一个方向梯度按设置的步长走。

主要问题包括：

①步长的选取，过大容易发散，过小速度太慢

②局部最优解和全局最优解选取的问题

Ⅲ、标准化

标准化的作用：因为不同的特征取值范围不同，比如年龄和工资，取值范围肯定不同，故需要先对数据进行预处理，进行标准化归一化操作。拿到数据之后基本上都需要做一次标准化操作

(x-μ)/σ，其中μ为均值，x-μ将所有的数据向坐标原点靠拢，以原点为中心对称。σ为标准差，即数据的浮动范围，取值范围浮动大对应的标准差也大，除以σ使得数据取值范围差不多一致。
综合而言(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#更新参数

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

X_new_b.dot(theta)#得出预测的结果值
#很显然，跟上述的求解的theta_best一致
"""
array([[4.21509616],
       [9.75532293]])
"""

做个实验：不同学习率α对结果的影响

theta_path_bgd = []#保存theta值用于后续的对比试验
def plot_gradient_descent(theta,eta,theta_path = None):#theta_path指定是否要保存当前的实验
    m = len(X_b)#一共的样本个数
    plt.plot(X,y,'b.')#原始样本数据进行散点图显示
    n_iterations = 1000#迭代次数
    for iteration in range(n_iterations):#遍历每次迭代
        y_predict = X_new_b.dot(theta)#得到当前的预测值
        plt.plot(X_new,y_predict,'b-')#将预测的线性拟合结果进行绘制
        gradients = 2/m* X_b.T.dot(X_b.dot(theta)-y)
        theta = theta - eta*gradients
        if theta_path is not None:#不是空值就进行保存
            theta_path.append(theta)
    plt.xlabel('X_1')
    plt.axis([0,2,0,15])#取值范围
    plt.title('eta = {}'.format(eta))#把学习率当成图的名称

theta = np.random.randn(2,1)#随机初始化一个θ值，因为需要偏置参数和权重参数故这里是两个

plt.figure(figsize=(10,4))#指定当前图的大小
plt.subplot(131)#1列3个，首先画第1个子图
plot_gradient_descent(theta,eta = 0.02)
plt.subplot(132)#1列3个，画第2个子图
plot_gradient_descent(theta,eta = 0.1,theta_path=theta_path_bgd)
plt.subplot(133)#1列3个，画第3个子图
plot_gradient_descent(theta,eta = 0.5)
plt.show()#每根线都代表每次迭代得出的结果

每条线都代表这每次迭代的结果，很明显可以看出，学习率越小越慢，但效果好；学习率越大反而效果不好，学跑偏了都。

②，随机梯度下降

theta_path_sgd=[]#保存theta值用于后续的对比试验
m = len(X_b)#计算当前所有的样本
np.random.seed(42)
n_epochs = 50#指定迭代次数

t0 = 5#分子 随便定义
t1 = 50#分母 随便定义

def learning_schedule(t):#学习率进行动态衰减，传入迭代次数t
    return t0/(t1+t)#t0和t1不变，t是迭代次数，随着迭代次数的增加，整体会变小

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

for epoch in range(n_epochs):#每个epoch表示完整的迭代一次所有的样本
    for i in range(m):#一个epoch中完成了所有的样本迭代
        if epoch < 10 and i<10:
            y_predict = X_new_b.dot(theta)#拿到预测结果值
            plt.plot(X_new,y_predict,'r-')#将预测结果值进行展示
        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)#保存当前的theta值
        
plt.plot(X,y,'b.')#绘制原始数据，以蓝点的形式展示
plt.axis([0,2,0,15])#横纵坐标的取值范围
plt.show()#每次运行都会得到不同的结果

③，MiniBatch(小批量)梯度下降

theta_path_mgd=[]#保存theta值用于后续的对比试验
n_epochs = 50#指定迭代次数
minibatch = 16
theta = np.random.randn(2,1)#对theta进行随机初始化
#np.random.seed(0)#当然也可以指定随机种子，保证每次随机选取，最后的theta值的不变
t0, t1 = 200, 1000#指定衰减策略
def learning_schedule(t):
    return t0 / (t + t1)
np.random.seed(42)
t = 0
for epoch in range(n_epochs):
    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):#从第0个样本到第m个样本，每个的大小为minibatch
        t+=1#记录一下迭代次数
        xi = X_b_shuffled[i:i+minibatch]#从当前指定的位置开始，操作minibatch个
        yi = y_shuffled[i:i+minibatch]#标签也一样
        gradients = 2/minibatch* xi.T.dot(xi.dot(theta)-yi)#根据minibatch梯度计算公式计算梯度
        eta = learning_schedule(t)
        theta = theta-eta*gradients
        theta_path_mgd.append(theta)

theta
"""
array([[4.25490684],
       [2.80388785]])
"""

Ⅴ、三种梯度下降策略对比

#为了方便后续的操作，将保存过的梯度都转换为array格式
theta_path_bgd = np.array(theta_path_bgd)#批量梯度下降
theta_path_sgd = np.array(theta_path_sgd)#随机梯度下降
theta_path_mgd = np.array(theta_path_mgd)#MiniBatch梯度下降

plt.figure(figsize=(12,6))#指定画图的最终大小
plt.plot(theta_path_sgd[:,0],theta_path_sgd[:,1],'r-s',linewidth=1,label='SGD')#随机梯度下降，对两个参数分别展示
plt.plot(theta_path_mgd[:,0],theta_path_mgd[:,1],'g-+',linewidth=2,label='MINIGD')#MiniBatch梯度下降
plt.plot(theta_path_bgd[:,0],theta_path_bgd[:,1],'b-o',linewidth=3,label='BGD')#批量梯度下降
plt.legend(loc='upper left')#label位置进行设置放在左上角
plt.axis([3.5,4.5,2.0,4.0])
plt.show()
"""
蓝色为批量梯度下降，直接奔着结果去的，不过很慢
绿色为MiniBatch梯度下降，随机选取其中的几个样本为参数，因为每个样本都不同，故有点跌宕起伏的感觉（常用，参数常用64、128、256，在速度能容忍的前提下参数越大越好）
红色为随机梯度下降，速度很快，跟个小傻子似的，好坏程度完全取决于随机样本选取的好坏
"""

二、多项式回归

import numpy as np
import os
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
plt.rcParams['axes.labelsize'] = 14
plt.rcParams['xtick.labelsize'] = 12
plt.rcParams['ytick.labelsize'] = 12
import warnings
warnings.filterwarnings('ignore')
np.random.seed(42)

m = 100#样本数据个数
X = 6*np.random.rand(m,1) - 3#横坐标的取值范围[-3,3]
y = 0.5*X**2+X+np.random.randn(m,1)#自己随便定义一个函数方程 y=0.5* x² + x + b，其中b为高斯随机抖动

plt.plot(X,y,'b.')#绘制随机选取的离散点
plt.xlabel('X_1')#横坐标名称
plt.ylabel('y')#纵坐标名称
plt.axis([-3,3,-5,10])#x、y的取值范围
plt.show()

from sklearn.preprocessing import PolynomialFeatures#导入工具包
poly_features = PolynomialFeatures(degree = 2,include_bias = False)#对函数进行实例化操作
X_poly = poly_features.fit_transform(X)#fit执行所有的计算操作，transform将所有的计算结果进行整合回传
X[0]#当前的x值
"""
array([2.82919615])
"""

X_poly[0]#[x,x²]
"""
array([2.82919615, 8.00435083])
"""

2.82919615 ** 2#很显然就是x²值
"""
8.004350855174822
"""

from sklearn.linear_model import LinearRegression#导包。开始训练
lin_reg = LinearRegression()
lin_reg.fit(X_poly,y)
print (lin_reg.coef_)#权重参数，即回归方程y=1.1087967x + 0.53435278x² - 0.03765461
print (lin_reg.intercept_)#与刚开始自定义的方程相比可知拟合效果还是很不错的
"""
[[1.10879671 0.53435287]]
[-0.03765461]
"""

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--',label='prediction')#再通过得出的预测值画曲线
plt.axis([-3,3,-5,10])#显示的时候限制一下取值范围
plt.legend()#加上标签
plt.show()

Ⅰ、根据不同degree值(不同多项式)进行拟合

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler#标准化操作
plt.figure(figsize=(12,6))
for style,width,degree in (('g-',1,100),('b--',1,2),('r-+',1,1)):#颜色、宽度、degree值
    poly_features = PolynomialFeatures(degree = degree,include_bias = False)#对函数进行实例化操作
    std = StandardScaler()#标准化 实例化
    lin_reg = LinearRegression()#线性回归 实例化
    polynomial_reg = Pipeline([('poly_features',poly_features),
             ('StandardScaler',std),
             ('lin_reg',lin_reg)])#车间三部曲
    polynomial_reg.fit(X,y)#传入数据
    y_new_2 = polynomial_reg.predict(X_new)#对测试点预测
    plt.plot(X_new,y_new_2,style,label = 'degree   '+str(degree),linewidth = width)
plt.plot(X,y,'b.')#画原始数据点
plt.axis([-3,3,-5,10])#显示的时候限制一下取值范围
plt.legend()
plt.show()

很显然，绿色函数degree为100，函数很复杂，已经过拟合了，为了尽可能的去满足所有的点，已经没必要了。

Ⅱ、不同样本数量对结果的影响

from sklearn.metrics import mean_squared_error#MSE均方误差
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=100)#测试集20%、训练集80%、随机种子为100，也就是每次数据切分都是相同的方式
    train_errors,val_errors = [],[]#这里没有使用到测试集，使用了是训练集和验证集，保存
    for m in range(1,len(X_train)):#从1个样本训练到使用所有的测试集数据样本训练
        model.fit(X_train[:m],y_train[:m])#训练
        y_train_predict = model.predict(X_train[:m])#得出训练集的结果，依次训练1个2个直到所有样本都训练
        y_val_predict = model.predict(X_val)#得出验证集的结果，验证集选择全部的，否则不公平
        train_errors.append(mean_squared_error(y_train[:m],y_train_predict[:m]))#使用m个样本进行训练，到时候预测也得是m个
        val_errors.append(mean_squared_error(y_val,y_val_predict))#使用所有样本进行验证，到时候预测个数也得一致
    plt.plot(np.sqrt(train_errors),'r-+',linewidth = 2,label = 'train_error')
    plt.plot(np.sqrt(val_errors),'b-',linewidth = 3,label = 'val_error')
    plt.xlabel('Trainsing set size')
    plt.ylabel('RMSE')
    plt.legend()#显示label

lin_reg = LinearRegression()#线性回归
plot_learning_curves(lin_reg,X,y)
plt.axis([0,80,0,3.3])#x、y轴取值范围
plt.show()

图中展示的是error，当训练集样本数量较少时，训练集error较小效果较好，但验证集validation的error较大，两者的差值较大，效果不好。

数据量越少，训练集的效果会越好，但是实际测试效果很一般。实际做模型的时候需要参考测试集和验证集的效果。

Ⅲ、多项式回归的过拟合风险

当degree 过大，也就是多项式太高，过拟合太严重。

polynomial_reg = Pipeline([('poly_features',PolynomialFeatures(degree = 25,include_bias = False)),
             ('lin_reg',LinearRegression())])
plot_learning_curves(polynomial_reg,X,y)
plt.axis([0,80,0,5])
plt.show()

训练集和验证集相差太明显了，可以等价于训练集为平时的刷题，验证集为期中期末考试。平常刷题错误率很低，最后的期中期末考试却错了一塌糊涂，这就是过拟合。

三、正则化

出现过拟合咋办？正则化专门解决过拟合问题的。

正则化实际上就是对权重参数进行惩罚，让权重参数尽可能平滑一些，有两种不同的方法来进行正则化惩罚：

θ1和θ2与x相乘结果都一样，要你你会选择哪个？小傻子都知道，会选择θ2，因为θ2考虑了三个元素，而θ1只考虑了第一个元素。

此时正则化公式就出现了：为了就是对每个θ参数进行选取最优

Ⅰ、岭回归

其中J(θ)为损失函数，当然越小越好。前面的均方误差是一样的，后面的为整数，越小越好。带入公式可得，θ2效果最佳。

from sklearn.linear_model import Ridge
np.random.seed(42)#构建随机种子
m = 20#指定样本数量
X = 3*np.random.rand(m,1)#x的取值范围[0,3]之间选取20个点就行了
y = 0.5 * X +np.random.randn(m,1)/1.5 +1#自己定义一个函数方程
X_new = np.linspace(0,3,100).reshape(100,1)#测试数据

def plot_model(model_calss,polynomial,alphas,**model_kargs):
    for alpha,style in zip(alphas,('b-','g--','r:')):#选取多个α值
        model = model_calss(alpha,**model_kargs)#对model实例化
        if polynomial:
            model = Pipeline([('poly_features',PolynomialFeatures(degree = 10,include_bias = False)),
             ('StandardScaler',StandardScaler()),
             ('lin_reg',model)])
        model.fit(X,y)
        y_new_regul = model.predict(X_new)#预测一下当前得到的结果
        lw = 2 if alpha > 0 else 1#指定线条宽度
        plt.plot(X_new,y_new_regul,style,linewidth = lw,label = 'alpha = {}'.format(alpha))
    plt.plot(X,y,'b.',linewidth =3)
    plt.legend()#显示标签

plt.figure(figsize=(14,6))#设置一个大图展示
plt.subplot(121)#两个子图的第一个
plot_model(Ridge,polynomial=False,alphas = (0,10,100))#不做多项式拟合操作  α的取值表示正则化的程度
plt.subplot(122)#两个子图的第二个
plot_model(Ridge,polynomial=True,alphas = (0,10**-5,1))#做多项式拟合操作   α的取值表示正则化的程度
plt.show()

Ⅱ、lasso

from sklearn.linear_model import Lasso

plt.figure(figsize=(14,6))
plt.subplot(121)
plot_model(Lasso,polynomial=False,alphas = (0,0.1,1))
plt.subplot(122)
plot_model(Lasso,polynomial=True,alphas = (0,10**-1,1))
plt.show()