AI基础：线性回归及其最小二乘法和梯度下降法详细推导与代码示例

线性回归：Liner Regression
主要是回忆一下最小二乘和梯度下降

什么是线性回归

• 线性：两个变量之间的关系是一次函数关系的图象是直线，叫做线性。
• 非线性：两个变量之间的关系不是一次函数关系的图象不是直线，叫做非线性。
• 回归：人们在测量事物的时候因为客观条件所限，求得的都是测量值，而不是事物真实的值，为了能够得到真实值，无限次的进行测量，最后通过这些测量数据计算回归到真实值，这就是回归的由来。
• 回归又可理解为拟合，线性回归就是一批数据集拟合（回归）出一条直线

线性回归能做什么

• 对大量的观测数据进行处理，从而得到比较符合事物内部规律的数学表达式。
• 也就是说寻找到数据与数据之间的规律所在，从而就可以模拟出结果，也就是对结果进行预测。
• 解决的就是通过已知的数据得到未知的结果。例如：对房价的预测、判断信用评价、电影票房预估等。
• 且可以得到连续的预估值，连续就区别与离散，离散就是一个一个点

线性回归一般表达式

$$\hat{y}=wx+b$$

其中 w 是参数，b是偏置项

实际上我们训练模型的过程中就是为了计算出 w 和 b，只有w和b出来了才可以对新的x数据进行预测它的y是啥。

后面的篇幅都是在说如何根据历史数据x,y解出 w 和 b，就是计算出历史数据拟合的那条线的表达式。

如何计算（学习）参数w，b

我们找到的参数理应使得每一项通过w，b计算得到的y与实际y值的误差最小。有了这个关系，我们就可以进行计算符合这关系的w，b值

对于每个样本的误差：
$${\epsilon }_i=y_i-{\hat{y}}_i$$

通过w，b计算得到的y值和实际y值之间的差值就是损失值对吧，每一项x-y对应的损失值都是不一样的，这就有个关系表达式来表示每项x-y对应的损失值的关系，这个表达式就叫做损失函数（Loss Function）

针对任何模型求解问题，都是最终都是可以得到一组预测值y^ ，对比已有的真实值 y ，数据行数为 n ，可以将损失函数定义如下：
$$L=\frac{1}{n}\underset{i=1}{\stackrel{n}{\Sigma }}{\left({\hat{y}}_i-y_i\right)}^2$$
即预测值与真实值之间的平均的平方距离，统计中一般称其为MSE(mean square error)均方误差。把之前的函数式代入损失函数，并且将需要求解的参数w和b看做是函数L的自变量，可得
$$L\left(w,b\right)=\frac{1}{n}\underset{i=1}{\stackrel{n}{\Sigma }}{\left(w{x}_i+b-y_i\right)}^2$$
那我们的任务是求解最小化L时w和b的值，即核心目标优化式为
$$\left(w^{\ast },b^{\ast }\right)=arg\underset{\left(w,b\right)}{\min }\underset{i=1}{\stackrel{n}{\Sigma }}{\left(w{x}_i+b-y_i\right)}^2$$

求解损失函数最小化L时w和b值的方法：最小二乘法

最小二乘法 least square method

步骤:

1>直接对损失函数求参数的偏导
2>令导数为0，导数为0的时候二次函数取得最小值

求解 w 和 b 是使损失函数最小化的过程，在统计中，称为线性回归模型的最小二乘“参数估计”(parameter estimation)。我们可以将 L(w,b) 分别对 w 和 b 求导，并令导数为0进行求解，过程如下

$$\begin{gathered} \text{线性回归最小二乘求解} \\ L={\Sigma }_{i=1}^N{\left(w{x}_i+b-y_i\right)}^2 \\ \text{分别对}b\text{和}w\text{进行求导} \end{gathered}$$
$$\begin{gathered} \text{对}b\text{求导：}\frac{\partial L}{\partial b}={\Sigma }_{i=1}^{N}2\left(w{x}_i+b-y_i\right)\cdot 1 \\ =2{\Sigma }_{i=1}^N\left(w{x}_i+b-y_i\right) \\ =w{\Sigma }_{i=1}^N{x}_i+{\Sigma }_{i=1}^{N}b-{\Sigma }_{i=1}^N{y}_i \\ \text{令：}\bar{x}=\frac{1}{N}{\Sigma }_{i=1}^N{x}_i,\bar{y}=\frac{1}{N}{\Sigma }_{i=1}^N{y}_i \\ \text{则上式继续}=w\cdot N\cdot \bar{x}+N\cdot b-N\cdot \bar{y} \\ \text{令导数为}0\text{：}w\cdot N\cdot \bar{x}+N\cdot b-N\cdot \bar{y}=0 \\ \text{则}b=\bar{y}-w\bar{x} \end{gathered}$$

$$\begin{gathered} \text{对}w\text{求导：}\frac{\partial L}{\partial w}={\Sigma }_{i=1}^{N}2\left(w{x}_i+b-y_i\right)\cdot x_i \\ \text{把}b\text{代入}={\Sigma }_{i=1}^{N}2\left(w{x}_i+\bar{y}-w\bar{x}-y_i\right)\cdot x_i \\ ={\Sigma }_{i=1}^N{w{x}_i}^2+{\Sigma }_{i=1}^N\bar{y}\cdot x_i-{\Sigma }_{i=1}^{N}w\bar{x}\cdot x_i-{\Sigma }_{i=1}^N{y}_i\cdot x_i \\ \text{令导数为}0\text{：}w\left({\Sigma }_{i=1}^N{x_i}^2-{\Sigma }_{i=1}^N\bar{x}\cdot x_i\right)={\Sigma }_{i=1}^N{y}_i\cdot x_i-{\Sigma }_{i=1}^N\bar{y}\cdot x_i \\ \text{令}\overline{x^2}=\frac{1}{N}{\Sigma }_{i=1}^N{x_i}^2,\overline{xy}=\frac{1}{N}{\Sigma }_{i=1}^N{x}_i{y}_i \\ \text{则上式子继续：}w\cdot \left(N\cdot \overline{x^2}-N\cdot {\bar{x}}^2\right)=N\cdot \overline{xy}-N\cdot \bar{y}\cdot \bar{x} \\ w=\frac{\overline{xy}-\bar{y}\cdot \bar{x}}{\overline{x^2}-{\bar{x}}^2} \end{gathered}$$

其中几个均值的转换，在实践开发中是非常方便的，可以很容易的计算得到 (xy)的均值、x的均值、y的均值、x平方的均值、x均值的平方

这也是线性规划很迅速训练的原因，计算很快

代码实现

import numpy as np # 引用numpy库,主要用来做科学计算 
import matplotlib.pyplot as plt # 引用matplotlib库,主要用来画图 
# 创建数据集,把数据写入到numpy数组 
data = np.array([[152,51],[156,53],[160,54],[164,55], 
        [168,57],[172,60],[176,62],[180,65], 
        [184,69],[188,72]]) 
# 打印大小 
x, y = data[:,0], data[:,1] 
print (x.shape, y.shape) 
# 1. 手动实现一个线性回归算法
#  实现w和b参数, 这里w是斜率, b是偏移量 完全照抄上述公示推导的结论哈
w = (np.mean(x * y) - np.mean(y)*np.mean(x)) / (np.mean(x**2) - np.mean(x)**2)
b = np.mean(y) - w * (np.mean(x))
print ("通过手动实现的线性回归模型参数: %.5f %.5f"%(w,b)) 
# 模型函数 = w * x + b
print ("通过手动实现的线性回归参数的预测: %.5f"%(w * 190 + b)) 
# 对x每个值取一个预测值
y_pred = []
for i in x:
    y_pred.append(w * i + b)
print ("通过手动实现的线性回归x的所有预测值: ", y_pred) 


# 2. 使用sklearn来实现线性回归模型, 可以用来比较一下跟手动实现的结果 
from sklearn.linear_model import LinearRegression 
model = LinearRegression().fit(x.reshape(-1,1),y)
sk_w = model.coef_
sk_b = model.intercept_
print ("基于sklearn的线性回归模型参数:%.5f %.5f"%(sk_w, sk_b)) 
print("基于sklearn的线性回归预测：%.5f " %model.predict(np.array([190]).reshape(-1, 1))[0])

测试结果

(10,) (10,)
通过手动实现的线性回归模型参数: 0.57576 -38.07879
通过手动实现的线性回归参数的预测: 71.31515
通过手动实现的线性回归x的所有预测值:  [49.43636363636362, 51.739393939393935, 54.04242424242423, 56.345454545454544, 58.64848484848484, 60.95151515151514, 63.25454545454545, 65.55757575757575, 67.86060606060606, 70.16363636363636]
基于sklearn的线性回归模型参数:0.57576 -38.07879
基于sklearn的线性回归预测：71.31515 

画图查看拟合

#绘图进行比较
plt.figure(figsize=(10,7))       #画布大小 
plt.scatter(x,y,label='target',color='blue')      #目标取值
plt.plot(x,y_pred ,label='preds',color='red')        #预测取值
plt.legend(loc='upper right')  #线条显示位置
plt.show()

求解损失函数最小化L时w和b值的方法：梯度下降法

梯度下降 gradient descent 非常常用的方式

梯度下降核心内容是对自变量进行不断的更新（针对w和b求偏导），使得目标函数不断逼近最小值的过程

梯度下降再细分又有随机梯度下降SGD、批量梯度下降BGD、MBGD小批量梯度下降，这里暂时不详细描述，后面再聊

目标：寻找某个凸函数的最小值。

凸函数： 某段区间只有一个低谷的函数

含义：朝着导数变小的方向逐步更新参数（这里的逐步，每一步，SGD就是每一步随机取样本，BGD就是每一步批量的使用样本，其他同理）
梯度：导数的方向

步骤：（这个步骤是梯度下降的逐步下降步骤，想像下盲人下山一步一步的走的过程。代码是梯度下降每次更新参数的过程比较清晰）

• 随机取一个值 x0；
• 算出对应的 f(x0)；
• 计算 f(x0)处损失函数 f(x)的导数；
• 从 f(x0)开始，沿着该处导数的方向增加𝛼，取值为x1，也就是说：|x0-x1|/𝛼的值为f(x) 在f(x0)的斜率；
• 重复循环2-4步，直到达到退出条件(比如设定循环N次后停止)

$$\begin{gathered} w\leftarrow w-\alpha \frac{\partial L}{\partial w} \\ b\leftarrow b-\alpha \frac{\partial L}{\partial b} \end{gathered}$$

关键超参数：步长的设置，步长，就是𝛼，每次我走多长，走太长可能跳过了低估，走太短计算就很慢

关于w和b的求导（梯度）先放这里，方便参看：注意系数2化简划掉了
$$\begin{gathered} L={\Sigma }_{i=1}^N{\left(w{x}_i+b-y_i\right)}^2 \\ \text{对}b\text{求导：}\frac{\partial L}{\partial b}={\Sigma }_{i=1}^N\left(w{x}_i+b-y_i\right) \\ \text{对}w\text{求导：}\frac{\partial L}{\partial w}={\Sigma }_{i=1}^N\left(w{x}_i+b-y_i\right)\cdot x_i \end{gathered}$$

这里的难点就在于：局部极小值 。解决：优化方法有一堆，凸优化啥的，这个是一个大课题，这里不展开了，总之知道这里就是逐步逼近最小值，从而得到对应的w和b

后面再开个梯度下降的专题聊聊

代码实现

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# 测试数据（仅供娱乐。。。。。）
data = np.array([[152,51],[156,53],[160,54],[164,55], 
        [168,57],[172,60],[176,62],[180,65], 
        [184,69],[188,72]]) 

x, y = data[:,0], data[:,1] 
# 画图查看测试数据分别
plt.figure(figsize=(5,3))
# plt.subplot(221)
plt.scatter(x,y,label='target')
plt.show()

# ===========================梯度下降训练定义====================
# 记录每次梯度计算的mse值，方便画图观察
test_epochs = []
test_mse = []
# GD(每次都使用全量数据进行计算)
def gradient_descent(x, y, b, w, lr, epochs):
    global test_epochs
    global test_mse

    # 循环epochs次
    for i in range(epochs):
        b_grad = 0
        w_grad = 0
        # 计算梯度的总和再求平均（套用导数后的公式）
        for j in range(0, len(x)):
            b_grad += (((w * x[j]) + b) - y[j])
            w_grad += (((w * x[j]) + b) - y[j]) * x[j]
        # 更新b和w
        b = b - (lr * b_grad)
        w = w - (lr * w_grad)
        # 打印每次迭代后的损失
        cur_mse = compute_error(b, w, x, y)
        # print("epochs {0} , b = {1}, w = {2}, error = {3}".format(i, b, w, cur_mse))
        test_epochs.append(i)
        test_mse.append(cur_mse)

    return b, w

# 计算MSE
def compute_error(b, w, x, y):
    totalError = 0
    for i in range(0, len(x)):
        totalError += (y[i] - (w * x[i] + b)) ** 2
    return totalError

# 训练
# 初始化学习率learning rate，参数b和w并设置最大迭代次数epochs
lr = 0.000001
b = 0
w = 0
epochs = 50
b_result, w_result = gradient_descent(x, y, b, w, lr, epochs)
print("="*50)
print("梯度下降计算得到的参数：w:{0},b:{1}".format(w_result, b_result))
# 画图查看损失值变化
# plt.subplot(222)
plt.figure(figsize=(5,3))
plt.plot(test_epochs,test_mse)
plt.xlabel("epochs")
plt.ylabel("MSE")
plt.show()

# ==============================结果测试==========================
# 使用计算得到的w，b进行测试
# 对x每个值取一个预测值
y_pred = []
for i in x:
    y_pred.append(w_result * i + b_result)
print ("通过手动实现的GD线性回归x的所有预测值: ", y_pred) 

# 画图查看差异
plt.figure(figsize=(5,3))       #画布大小 
plt.scatter(x,y,label='target',color='blue')      #目标取值
plt.plot(x,y_pred ,label='preds',color='red')        #预测取值
plt.legend(loc='upper right')  #线条显示位置
plt.show()

综上，从测试过程来看，学习率和迭代次数非常影响最终效果和计算压力，这俩超参数重点关注

多项式的回归

• 原理：前面提到，根据样本的趋势来决定假设函数，当趋势较为复杂时，线性模型就不适用了。就是说一条直线无法拟合样本趋势了
• 多项式回归可以理解为多个线性回归组织起来来拟合复杂的非线性的样本趋势

一元多项式回归：

公式
$$h_{\theta }\left(x\right)={\theta }_0+{\theta }_1{x}_1+{\theta }_2{x_2}^2+{\theta }_3{x_3}^3+{\theta }_4{x_4}^4...$$
这里不深入探讨多项式回归损失函数那些了，上面详细的聊是清晰线性回归过程，多项式回归其实就是多维特征。

代码实现

Sklearn里的求解

import numpy as np 
from sklearn.linear_model import LinearRegression 
# 生成样本数据, 特征维度为2 
X = np.array([[1, 1], [1, 2], [2, 2], [2, 3]]) 
# y = 1 * x_0 + 2 * x_1 + 3 
y = np.dot(X, np.array([1, 2])) + 3 
# 先使用sklearn自带的库来解决 
model = LinearRegression().fit(X, y) 
# 打印参数以及偏移量（bias） 
print ("基于sklearn的线性回归模型参数为 coef: ", model.coef_, " intercept: %.5f" %model.intercept_) 

输出：
基于sklearn的线性回归模型参数为 coef: [1. 2.] intercept: 3.00000

过拟合、欠拟合、正则化

这里只略微提一下，这仨内容很多，暂不展开

• 过拟合，测试数据拟合的非常好，误差很小，真实去用时，效果很差
• 欠拟合就是拟合的不好呗，误差还是很大（如上GD做的测试拟合效果很差）
• 正则化就是解决过拟合问题的方式：原理就是给计算出来的参数增加惩罚项，让它不要发挥那么大作用。详细可以BD下或者关注后续文章吧