关于线性回归的思考与部分代码

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

---

前言

本文仅作为记录自己在学习时获得的部分感悟用，如果能为各位朋友带来一点对于相关知识的思考，那便是我的荣幸。

（请注意：1、本文只是学习事情的自我思考，其中所写的内容或代码并不一定完全正确，如果您在学习过程中遇到问题或有其他更好的解法可以在评论或留言，或移步其他学习文档寻求帮助。2、在思考的过程中，由于个人水平有限，笔者也会有许多不理解的问题，因此如果文档中有错误或本文档无法解决您的疑问实属正常，欢迎指正，谢绝人身攻击。3、部分代码来源于网络，侵权即删。）

一、线性回归是什么？

在数学理论模型中，我们希望一个y值的确定能够通过一个函数和其中的一个或几个变量被清晰而准确的反应，例如对于函数y=2x+1，当x取一个确定的值时，y也为一个确定值，但在实际生活中，受到多方面因素的影响，并不是所有的y值都精准落在函数上的（如图1）
很明显图1上的散列点几乎可以拟合某一条形式为y=wx + b的函数，其中x为该函数中唯一的参数，而如何确定合适的w和b值，能让我没在新的x数据加入时以尽量小的误差去合理估计一个y值呢？这便是线性回归要解决的问题。

二、线性回归的理论基础

1.代价函数——Cost Function

在写下这篇博客的时候，我默认朋友们几乎都是和我一样刚刚接触机器学习或深度学习，那么以一个普通人的视角来看，使用最笨的方法，我们应该如何确定一组w和b值，使得构造出的函数拥有最小的误差呢？

如果我们是一台不知疲倦的机器，那么我们当然可以一点一点改变w和b的值，每次改变后依次计算每一个通过x值计算出的预测y值与实际y值的之间的误差值的绝对值（或平方值），并且将所有这样的差值相加起来除以参与预测数据的个数后得到一个最后的和，这便是在当前w和b值下的平均误差，在经历多次改变w和b值后，我们会得到多个这样的平均误差，挑选其中最小的，我们就可以得到相对来讲最理想的w和b的组合。

幸运的是，尽管我们无法完成这繁琐而枯燥的工作，但我们手边的计算机可以，对我们来说，需要的则是整理我们的思路，将上述内容整理成数学公式，转换成我们所熟悉的python代码。

以下是代价函数的数学计算式，其中第一个式子的Θ0和Θ1就是b和w，二者可以被视为一组数据，用一个1x2的矩阵表示：

2.梯度下降算法

1.使用代价函数找出最小值

好的，现在我们有了许多误差值（代价），那么我们应该如何获取这些误差值中最小的一个呢？很容易想到，我们只需要在计算误差值结束时加上一个一个if判断即可：

if 计算出的误差值<变量记录的误差值：
变量记录的误差值 = 计算出的误差值
变量记录的w = 参与本次计算的w
变量记录的b = 参与本次计算的b

（注意：该段为笔者自己的理解）类似的操作应当可以让我们很直白的找出计算过程中最低的误差并获得对应的w和b值，但是由此我们很容易发现一个问题，如果我们需要w和b的值尽可能地精准，那么我们每次对w和b施加的变化必须非常小（有时我们还需要考虑w在增大而b在减小才能达到最小误差的情况），那么计算出的误差个数势必成千上万甚至更多，一次一次会使得运行效率大幅度降低，因此我们需要一种能够帮助我们更效率地找到这组数据中的最小值。

2.梯度下降算法

1.算法引导

假想一种极端的情况：我们假设b=0，只考虑y=wx中w值的变化，如果已知有一点a处在常规坐标系的第一象限中（不落在坐标轴上），此时让w从0慢慢递增，那么y=wx与点a的最短距离将会慢慢缩小直至点a落在y=wx这条直线上，此时如果w继续增大，那么y=wx与点a的最短距离将会慢慢增大，此时我们记录下a与y=wx之间距离的变化并化成图像，那么理论上图像应当类似一条一元二次方程曲线，如果我们希望找到点a落在y=wx时w的值，很明显，我们只需要找到一阶导数使得y=0的w值。

回到我们的代价函数，当我们不断更改w和b的时候，如果我们对w和b的初始值足够合理，变化率足够小，那么我们得到的由多组误差值绘制成的曲线也应当基本拟合一条一元二次方程曲线并且这条曲线有最小值（直觉上笔者认为会有一些特殊情况脱离上述范畴，但前辈们的代码实际上完全规避了这个问题，因此只要能帮助朋友们理解梯度下降算法的存在意义和使用方法，我们姑且不考虑这些特殊情况）。

那么现在万事俱备，我们该如何引导编译器一步步找到由多组误差值绘制成的一元二次曲线的一阶导数为0的点呢？当然，我们不可能继续让编译器从左到右让编译器依次查询横坐标对应的误差值，这实际上完全没有使用梯度下降算法。如果朋友们有兴趣，大可以在此暂停自己思考。

那么为了方便朋友们更好地理解，我希望大家能看一下下面将要引出的案例。

2.梯度下降算法案例

首先我们需要导入两个个常用的python包：numpy和matplotlib.pyplot

import numpy as np
import matplotlib.pyplot as plt

现在假设我们已知一条一元二次方程曲线：y=x^2+2x+5，它的一阶导数我们很容易就能得到：y’ = 2x + 2，取最低点x=-1时y=4，编写显示它们的代码：

fig, ax = plt.subplots()
temp = np.linspace(-10, 10, 200)
vir = temp ** 2 + 2 * temp + 5
cir = temp * 2 + 2
zero = 0*temp
ax.plot(temp, vir, linestyle="-")
ax.plot(temp, cir, linestyle="--")
ax.plot(temp,zero,linestyle="-",color="black")
ax.scatter(-1, 4,color="k")
plt.show()

以上代码除2-5行外均为绘图代码，因此我们只简单介绍一下2-5行代码
第二行实际上只是告诉编译器创建一个空间，起点为-10到终点为10，在这个空间内创建200个点。这段代码的意义是让这200个点作为x值参与接下来的运算。
第三行、第四行和第五行则对应着图像上的曲线、斜线与y值为0的直线。

仔细观察此图像，我们假设初始点为x=3，那么点就需要从曲线右侧一点一点向最低点逼近，此时一阶导数值y’>0；再次假设初始点为x=-3，那么点就需要从曲线左侧一点一点向最低点逼近，此时一阶导数值y’<0

如果能做到每次循环能够做到：
___/ x-一小段距离
x= |
___\ x+一小段距离
那么x不就能慢慢向最低点逼近了吗？

“幸运”的是，最低点的导数为0，那么减去一阶导数为0左边的一阶导数值最终会做减法，而减去一阶导数为0右边的一阶导数值其实会做加法
那么我们对于x的迭代式可以更新为如下式：
___/ x-ay’(x)
x= |
___\ x-ay’(x)
因为一阶导数或许会很大，因此我们让y’(x)乘上步长a并且保证a是一个很小的小数，这样就可以让x一点一点地逼近最低点（如果不乘步长a，那么假设我们初始点设为x=3,此处y’(3)=8，则x=3-8=-5，相对来讲其实我们离最低点反而更远了）。继续编写代码如下：

inter = 1000
alpha = 0.01
x = 5
for i in range(inter):
    x = x - alpha*(2*x+2)

以上inter为迭代次数，在一定范围内，迭代次数越大，得到的结果越精准，但超出范围后会使x在最低点左右不断来回。
alpha为此次案例中的步长a，一般来说a的选择应该谨慎，如果选取过大就会出现之前提到的离最低点反而更远的情况；如果选取过小那么一方面在迭代结束前可能x无法很好地逼近最低点，而延长迭代次数又会使运行效率下降。
x的初始值在本案例中选择为5，尽管x的选择一般是随意的，但在实际应用中应当尽可能的估计初始值以方便编译器的运算。
for循环中即是我们之前提到的x = x - a*y’(x)，执行1000次即要求x向最低点逼近1000次。

在编写完上述代码后，编译器已经可以做到通过不断更新x值逼近最低点了，考虑到我们之前对最低点(-1,4)是直接在图上绘制出来，现在对代码进行改写，让最终呈现在图上的是经过计算后由程序得出的点：

fig, ax = plt.subplots()
temp = np.linspace(-10, 10, 200)
vir = temp ** 2 + 2 * temp + 5
cir = temp * 2 + 2
zero = 0*temp
ax.plot(temp, vir, linestyle="-")
ax.plot(temp, cir, linestyle="--")
ax.plot(temp,zero,linestyle="-",color="black")
ilter = 1000
alpha = 0.01
x = 3
for i in range(ilter):
    x = x - alpha * (2 * x + 2)

ax.scatter(x, x ** 2 + 2 * x + 5,color="k")
print(x)
print(x ** 2 + 2 * x + 5)
plt.show()

最终，我们能得到如下图片以及编译器的输出：

3.小结与问题

那么现在恭喜朋友们完成了本案例，我在某天发呆时突然觉得部分理解了梯度下降算法，并尝试了上述的例子。
现在或许大家能够理解算法中“梯度”的含义：x是经由计算一段一段下降的，不必遍历过程中每一个点。
当然，这种办法其实存在着一些问题：如果代价函数拟合的图像是一个倒着的麦当劳图标怎么办（类似于一个"W"）？此时图像的一阶导数存在两个零点，那么根据初始点选择的不同，函数会得到不同的结果，类似这样的问题被称为“局部最优”，顾名思义，在局部范围内我们找到了最优点，但当范围扩大，这个点或许并不是最优的。如果你询问一些ai：“如何解决线性回归的局部最优问题？”有些ai会回答你：线性回归一般不会出现局部最优问题，仔细想想我们确实很难让拟合一条直线的一些点再去拟合另一条直线并且保证拟合效果相当好，但ai仍然给了我们解决办法：正则表达式，并且正则表达式的代码极为简洁、处理效率也更高。然而很可惜，这部分属于线性代数内容，而笔者大学生涯中的历次线性代数成绩几乎可以拟合一条y=0.1x的直线，x甚至是运气而不是学习时长，因此笔者给大家介绍正则表达式的内容无异于张飞雕花，幸而任何一个ai都可以相对详细地介绍正则表达式，简单搜索也能找出许多水平高出笔者一大截的博主，有兴趣的朋友可以自行查阅。

3.案例讲解与代码实例

1.获取数据

通过以上学习，我们了解了代价函数与梯度下降算法，是时候看看详细的案例了，在那之前我们需要一些数据，朋友们可以直接生成一些拟合某条一元一次方程图像的点，之后选取其中的点让它们的y值随机加减一些值，类似于使用pytorch的如下代码：

def synthetic_data(w, b, num_examples):
    x = torch.normal(0, 1, (num_examples, len(w)))
    y = torch.matmul(x, w) + b
    y += torch.normal(0, 0.01, y.shape)
    return x, y.reshape((-1, 1))

但我们默认大家在初识线性回归时还没有接触或使用过pytorch，因此如果大家需要数据而又懒得自己编写，可以直接创建一个名为"ex1data1"的txt文件，将如下数据按照97行2列的形式复制上去并保存在大家需要编写线性回归代码的python项目目录下：
6.1101,17.592
5.5277,9.1302
8.5186,13.662
7.0032,11.854
5.8598,6.8233
8.3829,11.886
7.4764,4.3483
8.5781,12
6.4862,6.5987
5.0546,3.8166
5.7107,3.2522
14.164,15.505
5.734,3.1551
8.4084,7.2258
5.6407,0.71618
5.3794,3.5129
6.3654,5.3048
5.1301,0.56077
6.4296,3.6518
7.0708,5.3893
6.1891,3.1386
20.27,21.767
5.4901,4.263
6.3261,5.1875
5.5649,3.0825
18.945,22.638
12.828,13.501
10.957,7.0467
13.176,14.692
22.203,24.147
5.2524,-1.22
6.5894,5.9966
9.2482,12.134
5.8918,1.8495
8.2111,6.5426
7.9334,4.5623
8.0959,4.1164
5.6063,3.3928
12.836,10.117
6.3534,5.4974
5.4069,0.55657
6.8825,3.9115
11.708,5.3854
5.7737,2.4406
7.8247,6.7318
7.0931,1.0463
5.0702,5.1337
5.8014,1.844
11.7,8.0043
5.5416,1.0179
7.5402,6.7504
5.3077,1.8396
7.4239,4.2885
7.6031,4.9981
6.3328,1.4233
6.3589,-1.4211
6.2742,2.4756
5.6397,4.6042
9.3102,3.9624
9.4536,5.4141
8.8254,5.1694
5.1793,-0.74279
21.279,17.929
14.908,12.054
18.959,17.054
7.2182,4.8852
8.2951,5.7442
10.236,7.7754
5.4994,1.0173
20.341,20.992
10.136,6.6799
7.3345,4.0259
6.0062,1.2784
7.2259,3.3411
5.0269,-2.6807
6.5479,0.29678
7.5386,3.8845
5.0365,5.7014
10.274,6.7526
5.1077,2.0576
5.7292,0.47953
5.1884,0.20421
6.3557,0.67861
9.7687,7.5435
6.5159,5.3436
8.5172,4.2415
9.1802,6.7981
6.002,0.92695
5.5204,0.152
5.0594,2.8214
5.7077,1.8451
7.6366,4.2959
5.8707,7.2029
5.3054,1.9869
8.2934,0.14454
13.394,9.0551
5.4369,0.61705
请注意，这些数据应当有97行2列，请记住这两个数据，方便在后面进行数据预处理的时候理解。

2.导入库并录入数据

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

path = "ex1data1.txt"
data = pd.read_csv(path, header=None, names=['Population','Profit'])
#plt.figure(figsize=(12, 8))
#plt.scatter(data.Population, data.Profit)
#plt.show()

导入常用的三个包后，我们将文件的路径赋值到path中，在对data赋值时获取path所指向的文件，将两列数据分别命名为“Population”和“Profit”，header的作用稍微复杂一些，简单来讲header参数主要用于指定 CSV 文件中用作列名的行，如果header的值为None，数据框的列名会自动变为 0 和 1，但因为我们使用names给列进行了命名，相当于对0和1进行了覆盖。
被注释掉的三行代码用于可视化data所读取到的数据，如果读取正确，那么数据中的点应如下图散列：

3.数据预处理

这部分的代码有些抽象，如果您的手边有笔记本，那么建议您在读代码的同时将出现的变量形状写出来，或者您可以参考我在每行代码后的注释理解，它们代表了在执行了对应行的代码后变量的形状。

data.insert(0,"Ones",1)#97*3
cols = data.shape[1]#3
X = data.iloc[:,0:cols-1]#97*2
Y = data.iloc[:, cols-1:cols]#97*1
theta = np.matrix(np.array([1,1]))#1*2

X = np.matrix(X.values)
Y = np.matrix(Y.values)

我们来依次解释每一行代码：
第一行在data数据中插入了一列新的数据，插入列的位置在第一列，列名为“Ones”，列中所有的值均为1。于是data形状变成了97x3
第二行创建了一个名为cols的变量，将其赋值为data的列数，因为目前data有3列，因此cols值为3。

第三行创建X变量，对data进行切片处理，获取data的所有行与第一第二列，即名为"Ones"且所有值均为1的列和名为"Population"的列，将这两列数据全部赋值给X。
第四列创建Y变量，用同样的方法获取了data的最后一列，即名为“Profit”的列。
这里我们稍微介绍一些iloc函数：iloc函数是定义在pandas包中的函数，主要作用就是对数据进行切片，一般使用时会传入两组数据，以逗号为分界，前面代表需要获取的行，逗号后代表需要获取的列，使用时前后都需要有冒号的参与，冒号前代表从哪个行/列开始，冒号后则代表结束，如果冒号前后没有值，则代表获取所有行/列。需要注意冒号前后所代表的数字是左闭右开的，例如在对X赋值的过程中，第一个冒号前后无数字，代表需要获取data的所有行传递给X；第二个冒号前后为0和cols-1即为0和2，这代表获取data的第1列和第2列，即列号为0、1的两列赋值给X。这也解释了为什么在赋值Y的列时第二个冒号后面是3，但我们根本没有列号为3的数据。

第五行代码创建了一个名为theta的变量，使用np.array直接确定theta的形状为1x2，值均为1，并且将其“矩阵化”后赋值给theta，关于矩阵化的函数matrix我们目前不用了解太多，只需要知道矩阵化后才可以进行矩阵运算即可。而关于theta的作用，我们稍后做出解释。
类似的，之后的代码分别对X和Y都进行了矩阵化处理。

4.编写代价函数

根据我们之前的公式，代价函数可以编写如下：

def Cost(X,Y,theta):
    inner = np.power(X @ theta.T-Y,2)
    return np.sum(inner)/2*(len(X))

根据上述代码，大家应该可以理解theta的作用：theta其实就是b和w的值，我们在之前对他们都初始化值为1。
X@theta代表了预测值，97x2的矩阵乘以2x1的矩阵后会得到一个91x1的矩阵，这就是预测y值，减去实际的y值后平方后取平均值再除以2（这里除以2不影响最终结果的比较，同时可以在之后的求导式子更加简单）。

5.编写梯度下降函数

def gradientDescent(X, Y, alpha, iterations, theta):
    for i in range(iterations):
        h = X @ theta.T#97*1
        error = h - Y#97*1
        gradient = (X.T @ error) / X.shape[0]#2*1 = theta.T-->2*1
        new_theta = theta.T - alpha * gradient#2*1
        theta = new_theta.T#1*2

        #if i % 100 == 0:
            #cost = Cost(X, Y, theta)
            #print(f"Iteration {i}, Cost: {cost}")
    return theta#1*2

代码中的alpha代表步长，iterations代表迭代次数，h代表预测y值，error代表着误差，gradient则是计算出的error（代价）这一点的导数，new_theta和theta则只是为了在正确接收theta的更新值的同时让其形状与参与运算前相同。
在函数中我们创建了一个for循环，循环次数即我们需要的迭代次数，在每一次循环中计算出误差值error，根据公式（如下图）计算得出代价函数的导数——梯度，根据梯度更新导数并进行下一次计算。
如果按照代码一步一步画出运行中各个变量的形状，会发现h是97x1的矩阵，error是97x1的矩阵，X.T是2x97的矩阵，而X.T@error则是一个2x1的矩阵，我们原本传入的theta实际上是1x2的，因此代码中我新加入了new_theta保存更新计算后的theta值，并在使用新的theta前转置回原来的形状。
如果仔细思考这部分的代码，你会发现实际上我们并不需要代价函数参与运算，事实上如果如上述代码所示的if语句注释掉，我们之前所编写的Cost函数在代码中并不会被调用。这是因为当我们在gradientDescent函数中计算出的error值已经可以代表代价函数中的某个代价了，而利用公式则可以直接求出这个点的导数并对theta进行更新，因此其实并不需要显式的求出整个代价函数，在实际应用中，代价函数一般用于可视化函数图像并帮助我们找出它的导数，如果你希望让cost函数起到一些作用，将上面代码中的注释符去掉即可，你可以看到每一百次循环后的总误差值，这可以方便你对代码进行调试以便让你找到更合适的iterations。

6.设置参数、调用函数以及绘图

那么万事俱备，我们接下来只需要设定好alpha和iterations的值、调用gradientDescent函数并绘制图像即可。
如果您在之前没有学习过matplotlib，请放宽心，这些代码对matplotlib库的要求并不高，您大可在现在将它们复制粘贴到您的编译器中，之后按照您的喜好询问ai每个函数的作用，或者留作以后进行学习，不过出于方便，这里仍然给出每个函数作用的简单注释。

alpha = 0.01  # 学习率
iterations = 3000  # 迭代次数

theta = gradientDescent(X, Y, alpha, iterations, theta)  #获得更新结束后的theta值

x = np.linspace(data.Population.min(), data.Population.max(), 100)  #创建一个线性空间，空间中最小值是data.Population的最小值，最大值是data.Population的最大值，中间插入100个点充当x
f = theta[0,0] + theta[0,1] * x  #根据x和计算得出的theta构造一条f函数，稍后直接将这条函数展示在图像上

plt.figure(figsize=(12, 8))  #设置图像框的长和高
plt.xlabel('Population')  #设置x坐标名称
plt.ylabel('Profit')  #设置y坐标名称
l1 = plt.plot(x, f, label='Prediction', color='red')  #绘制l1曲线，其中横坐标取x的值，纵坐标取f的值，标签名为：Prediction，颜色设置为红色
l2 = plt.scatter(data.Population, data.Profit, label='Traing Data', )  #绘制各个数据点，其中横坐标为data.Population，纵坐标为data.Profit，标签名为：Traing Data
plt.legend(loc='best')  #添加一个图例，并自动选择最佳的位置放置图例。
plt.title('Predicted Profit vs Population Size')  #设置图像框的标题
plt.show()  #展示图像框

如果您的代码还保留着本小节第二部分中的

#plt.figure(figsize=(12, 8))
#plt.scatter(data.Population, data.Profit)
#plt.show()

并且没有将其注释掉，那么我建议您将其删去或注释即可，否则会生成两个图像框影响您观看。
最终的结果图应当如下图所示：

---

总结

终于，我们完成了我们的第一个线性回归项目，这对我们的学习来说是一个好的开始，尽管过程磕磕绊绊，我们仍然坚持了下来，但我们仍然留下了一些问题：
1.您是否想过上述的所有代码中，我们每次处理代价函数时的横坐标都是什么呢？是x还是theta？这似乎很容易导致思维混乱，也许这其实是个无伤大雅或者实际上无比简单的问题，但依然困扰了我很久，直到现在都没有理清楚。
2.关于我们之前提到过的正则表达式，实际上它的代码是这样的：

def normalEqn(X, Y):
    theta = np.linalg.inv(X.T@X)@X.T@Y
    return theta

theta = normalEqn(X, Y)
computeCost(X, Y, theta.reshape(1, -1))

你可以看到，它的代码简单到可怕，但运行效率却极高：在我们的gradientDescent代码中，iterations被设置为3000，而调用正则表达式时只需要1000次迭代即可达到最佳，但其原理我并未进行深入了解，仅仅只是学会线性回归就已经消磨掉了我的大部分精力，如果您有兴趣，可在之后进行了解。

以上就是笔者个人对线性回归的理解，也许不尽正确，但如果能对您的学习起到一点帮助，那这些文字就并非毫无意义。另外如果文章中出现了错误和纰漏，请您拨冗留言或私信，感激不尽！