一元线性回归（原理+python代码实现）,看完就能上手实践

问题背景

某公司每月都会举行<步行领奖金>的活动，具体发多少奖金由公司决定，公司有一套自己的计算规则，并且不会事先告知会发多少奖金。能够确定的是，步数越多，奖金越多。现在你想参加这个月的活动，想要预先算一下能领多少钱。你能拿到公司前几个月这个活动的数据，包括两列，一列是步数，一列是奖金金额。

解决方法

一元线性回归

什么是一元线性回归？

一元线性回归是理解两个变量（自变量和因变量，步数就是自变量，奖金金额就是因变量）之间关系的一种方法。

一元指的是只有一个自变量（或者叫特征、输入变量）。

线性指的是两个变量之间的关系是线性的。

一元线性回归的目标是找到一个直线，来拟合所有数据（<步行领奖金>活动的历史数据），使得通过这个直线，可以预测因变量（奖金金额）的值。

假如这个直线的方程式这样表示：

$y=wx+b$

这里，x是自变量，y是因变量，w是直线的斜率，b是直线的截距。

要想确定这个直线，就是要确定这里的w和b。

那么，我们来一步一步分析

什么样的直线是我们想要的直线？当然是可以"最好"的拟合所有数据（活动的历史数据）的直线。

什么是最好的拟合？怎么衡量最好的拟合所有数据呢？所有数据的真实值和预测值相差的最少。这里的真实值就是历史奖金金额，预测值就是通过直线预测的奖金金额。

那么问题来到了，用什么方法来衡量所有数据的真实值和预测值相差的大小呢？用下面这个公式：

$J(w,b)=\frac{1}{m}{\sum }_{i=1}^m{\left({\hat{y}}_i-y_i\right)}^2$

这里，m是数据集的数量，$\widehat{y_i}$是预测值，$y_i$是真实值。

公式的意思是：计算所有数据的真实值和预测值差的平方和，再取平均值。在机器学习中，这个公式被称作"代价函数"。

为什么有平方？因为衡量相差大小用的是距离，和正负没有关系。

为什么要求和？因为要拟合所有数据。

为什么要除以m取平均？为了避免数据量的多少影响代价函数的值。你想，如果没有最前边的 $\frac{1}{m}$，那么是不是数据越多，代价函数的值就越大。而我们最终的目标是不是要找到某个直线，使得代价函数的值最小。

到这里，接下来应该思考的就是如何找到使代价函数最小的那个直线？也就是要找到使代价函数最小的w和b。

最笨的办法是我们拿不同的w和b依次去试，看哪个组合求出来的J是最小的，就留下它。但这个办法显然不行，不知道要试多少次才算找到了最优解。所以，"梯度下降"方法就要登场了。

梯度下降

梯度下降算法就是根据函数的斜率（梯度），找到函数值下降最快的方向，然后朝着那个方向进行减小。重复这个过程，直到找到函数的局部最小值。

举个例子：想象你在一座大山上，想要找到山底。但是你眼睛被蒙上了，所以你不能直接看到山底。你的任务是找到最快的方法下山。这时候你会怎么做呢？只能通过脚下的坡度感受到山的陡峭程度。找到最陡峭的那个方向，然后走一小步，不断重复这个过程，直到最终走到山底。

在写梯度下降算法公式之前，我稍微改动一下代价函数J，但是改动并不会影响最终结果。

$J(w,b)=\frac{1}{2m}{\sum }_{i=1}^m{\left({\hat{y}}_i-y_i\right)}^2$

其实改动不大，就是乘了个$\frac{1}{2}$，为什么这么做呢？后面你看完公式就知道了，就是为了方便公式化简。虽然乘了0.5，但是不会影响最终结果。因为只要找到使J最小的w和b就好了。J是多少不重要。

梯度下降公式

$w=w-\alpha \frac{\partial }{\partial w}J(w,b)$
$b=b-\alpha \frac{\partial }{\partial b}J(w,b)$

这里，$\alpha$叫学习率，决定了下山步子的大小。以下是两个偏导数的推导过程。

$\begin{array}{rl}\frac{\partial }{\partial w}J(w,b)& =\frac{\partial }{\partial w}\frac{1}{2m}\sum _{i=1}^m{\left({\hat{y}}_i-y_i\right)}^2\\ \frac{\partial }{\partial w}J(w,b)& =\frac{\partial }{\partial w}\frac{1}{2m}\sum _{i=1}^m{\left({\hat{y}}_i-\left(w{x}_i+b\right)\right)}^2\\ \frac{\partial }{\partial w}J(w,b)& =2\times \frac{1}{2m}\sum _{i=1}^m\left({\hat{y}}_i-\left(w{x}_i+b\right)\right)\times -x_i\\ \frac{\partial }{\partial w}J(w,b)& =\frac{1}{m}\sum _{i=1}^m\left(\left(w{x}_i+b\right)-{\hat{y}}_i\right)x_i\\ \frac{\partial }{\partial w}J(w,b)& =\frac{1}{m}\sum _{i=1}^m\left(y_i-{\hat{y}}_i\right)x_i\end{array}$

另一个偏导数：

$\begin{array}{rl}& \frac{\partial }{\partial b}J(w,b)=\frac{\partial }{\partial b}\frac{1}{2m}\sum _{i=1}^m{\left({\hat{y}}_i-y_i\right)}^2\\ & \frac{\partial }{\partial b}J(w,b)=\frac{\partial }{\partial b}\frac{1}{2m}\sum _{i=1}^m{\left({\hat{y}}_i-\left(w{x}_i+b\right)\right)}^2\\ & \frac{\partial }{\partial b}J(w,b)=2\times \frac{1}{2m}\sum _{i=1}^m\left({\hat{y}}_i-\left(w{x}_i+b\right)\right)\times -1\\ & \frac{\partial }{\partial b}J(w,b)=\frac{1}{m}\sum _{i=1}^m\left(\left(w{x}_i+b\right)-{\hat{y}}_i\right)\\ & \frac{\partial }{\partial b}J(w,b)=\frac{1}{m}\sum _{i=1}^m\left(y_i-{\hat{y}}_i\right)\end{array}$

到此，公式推导完毕，接下来就是具体实践了，还记得本篇文章最开始的那个问题吗，以下代码就使用一元线性回归来解决这个问题。代码中有丰富的注释，保证看完就能理解。

注：以下所有代码的运行环境是（anaconda1.12 python3.11）

x = [34, 82, 76, 14, 66, 30, 68, 13, 56, 64, 98, 80, 33, 36, 36, 92, 18,37, 91, 33, 44, 32,  2, 16, 28, 35, 10, 47, 93,  9,  9, 74, 84, 69, 33, 88, 33, 28, 31, 51, 32,  2,  6, 25, 12, 16, 87,  8, 63, 40]
y = [94, 248, 193, 60, 177, 99, 165, 33, 128, 245, 254, 199, 126, 96, 144, 302, 13, 114, 324, 61, 108, 115, 48, 58, 121, 93, 69, 109, 283, 61, 54, 252, 301, 179, 135, 289, 100, 128, 83, 166, 69, -19, 10, 126, 92, 62, 291, 72, 205, 159]
x = np.array(x) # 转换为ndarray数据类型
y = np.array(y) # 转换为ndarray数据类型

def J(x,y,w,b): # 代价函数
    m = len(x) # 数量
    dif = w*x+b - y # 真实值和预测值的差
    dif_2 = dif*dif # 真实值和预测值的差的平方
    return np.sum(dif_2)/(2*m) # 真实值和预测值的差的平方和，再取平均值

def gradientDescent(x,y,w,b,alpha,iterations): # 梯度下降
    m = len(x) # 数量
    cost = [] # 保存代价函数求的值
    for i in range(iterations):
        w = w - alpha*(np.sum(w*x+b-y)/m) # 梯度下降法更新w
        b = b - alpha*(np.sum((w*x+b-y)*x)/m) # 梯度下降法更新b
        cost.append(J(x,y,w,b)) # 代价函数求的值存储到cost列表中
    return w,b,cost
 
result = gradientDescent(x,y,2,1,0.0001,500)

在运行代码的过程中，结合下面这两个图来不断调整w,b,alpha,iterations各个参数的值。直到找到最优解。注意这个过程需要反复调试。

plt.plot(np.arange(count), result[2]) # 折线图
plt.show()

• x轴表示迭代步数。
• y轴表示代价函数J的值。

plt.plot(x, result[0]*x+result[1], color='r') # 折线图
plt.scatter(x, y, color='c') # 散点图
plt.show()

• x轴表示步数
• y轴表示奖金金额
• 散点图表示的是原始数据的分布情况
• 折线图表示的是拟合出的直线

最终得到的直线方程式为：y=2.40x+36.96（保留了两位小数）

第二种代码实现是借助sklearn库实现的，直接调用函数接口即可，算法细节封装到了函数里。对于初学者，我不建议直接调库，而是要按照上面的方法，一步一步写代码，这对理解算法非常有帮助。如果已经非常熟练了，可以使用sklearn库实现。

from sklearn.linear_model import LinearRegression # 导入线性回归接口
model = LinearRegression() # 模型搭建
x = np.array([x]).reshape(50,1) # 模型训练要求输入的参数必须是二维的，所以将x的维度改为（50，1）
model.fit(x,y) # 模型训练

# 绘制原始数据的散点图，以及，一元线性方程
plt.scatter(x,y) # 散点图
plt.plot(x,model.predict(x), color='red') # 折线图
plt.show()

• x轴表示步数
• y轴表示奖金金额
• 散点图表示的是原始数据的分布情况
• 折线图表示的是拟合出的直线